We were among the first to use geometric models, based on Non-Uniform Rational B-Splines or NURBS, directly for analysis eliminating the meshing step needed in finite element analysis [1,2]. The use of geometric models for analysis is commonly referred to, as Isogeometric Analysis (IGA, [3]), and is a popular strategy for directly integrating CAD and CAE without needing an intervening mesh.

A challenging problem in IGA is stitching together, parametrically described spline patches. We developed a parametric stitching, or p-stitching, procedure for coupling patches with assured, arbitrary smoothness at the interface between the subdomains, including problems with sharp changes in gradient, as at dissimilar material interfaces [4]. The coupling procedure relies on Enriched Isogeometric Analysis (EIGA) previously developed by our group in 2012 [5].

Figure 1: Isoparametric descriptions of the NURBS geometrical surface and behavior for predictions of optimal shape of a droplet constrained between a cylindrical surface and a circular pad [1].

Using IGA to study crack initiation and propagation in semiconductor packaging

Automated crack initiation and propagation along arbitrary paths across interfaces in multilayer materials was first demonstrated by our group in 2012 [6]. In the figure below, simulations of arbitrary crack initiation and propagation in micron-scale semiconductor chip dielectric stacks using enriched isogeometric analysis is illustrated [7]. We later used EIGA for singular as well as discontinuous enrichments [8] based on asymptotic analysis of singularities in multi-material wedge corners and crack tips [9].

EIGA, mentioned earlier, utilizes measures of distance (algebraic level sets) constructed from the parametric enrichments; it explicitly preserves the CAD geometry of the interface allowing ease in calculating the curvature and normal that drive the interface motion under physical forces. The algebraic level sets are constructed from parametric boundaries using the theory of resultants of algebraic equations [10,11]. The algebraic level sets enable one to avoid the vexing CAD problem of needing to numerically calculate the intersection between complex three-dimensional parametric surfaces and stitching the surfaces together to form the common B-rep solid model (Figure 3).

Furthermore, a rigorous technique was developed by our group to recover implicitizations of parametric surfaces with trivially singular Dixon resultants and to carryout behavioral analysis of complex intersecting closed regions bounded by parametric surfaces without surface-surface intersection calculations [12].

Interfacial phenomena are of critical importance in many fields, however the conditions governing phase evolution at the interface are often unclear. We derived the thermodynamic configurational force associated with a moving interface with multiple diffusing species and arbitrary surface stress [13], inspired by the work of Truesdell and Toupin [14] and Gurtin [15]. The mass, momentum, energy balance as well as the second law condition were derived on the evolving phase interface subject to mechanical loads, heat and multiple diffusing species. The derived second law condition naturally extends the Eshelby energy–momentum tensor associated with moving cracks to include species diffusion terms (Eq. (43) in [13]), and naturally yields the interface equivalents of Fick’s (Eq. (41)) and Fourier’s laws (Eq. (42)). The second law restriction was then used to derive the condition for the growth and nucleation of new phases in a body undergoing finite deformation subject to inhomogeneous as well as anisotropic interface stress, and multiple diffusing species.

Figure 4: Explicit interface EIGA simulations of void and Cu6Sn5 intermetallic compound (IMC) growth in solder microbumps [17] (a) experimental observations of void growth after 1000 hours at 150 oC (b) initial mass concentrations of Sn in the model (c) mass concentrations of Sn after 200 hours at 150 oC (d) electric potential in the domain and void evolution after 200 hours at 150 oC under a current density of 𝟑 × 𝟏𝟎𝟓 𝑨/𝒄𝒎𝟐. The model accurately captures diffusion dominated IMC growth as well as the greater IMC growth at anode with addition of current.

We also posed a configurational optimization problem and derived the sensitivity of an arbitrary objective to arbitrary motions of one or more finite-sized heterogeneities inserted into a homogeneous domain [18]. In the derivation, they naturally obtained the definition of a generalized Eshelby energy-momentum tensor for arbitrary objectives and expressed the sensitivities as surface integrals with jump terms across the heterogeneity boundaries that vanished on homogeneous domains yielding generalized conservation laws for arbitrary objectives. They then showed that the specific path- independent forms of the sensitivity of the objective to arbitrary translation, rotation or scaling of the inserted heterogeneities naturally yielded the J-, L- and M-integrals of fracture mechanics when the objective is strain energy. The theory was computationally implemented using EIGA to optimally identify best/worst-case locations for line cracks that are inserted into the domain as well as to optimally mitigate the risk of fracture due to a crack at its worst-case location by sequentially inserting and optimizing the configurations of circular/elliptical stiff/soft inclusions (Figure 5). The simulations required singular as well as discontinuous enrichments mentioned earlier for modeling the cracks [8].

Solder System Constitutive Behavior and Fatigue

Solder joints are ubiquitous in microelectronic assemblies. They operate at high homologous temperatures, experience significant rate-dependent plastic deformation, and develop fatigue cracks that are large relative to the size of the joint. A significant effort in the Subbarayan group is to characterize the constitutive and fatigue behavior of a wide variety of solder alloys under a wide range of strain rates (Figure 6, [19,20,21]).

Figure 6: Micro- and nano-precision closed-loop mechanical testers with capacitively sensed feedback were built and used to characterize viscoplastic constitutive behavior and cyclic damage accumulation in solder joints varying in size from 100s of micrometers to 10s of micrometers.

Using point wise material damage described by a Weibull distribution of equivalent plastic strain, crack fronts and fatigue lives were predicted accurately to within 20% relative to experimental observations [22]. However, a Weibull damage description lacks physical foundation. So, we developed a maximum entropy theory for fatigue damage accumulation at crack tips in ductile solids based on the concepts of statistical mechanics and information theory (Figure 7, [23,24]).

Figure 7 A comparison of simulated and experimentally observed crack fronts in a Sn3.8Ag0.5Cu solder joint. The max entropy damage model used a single experimentally determined parameter to simulate thermal cycling fatigue crack initiation and propagation [24]. The predicted life was accurate to within 10% of the experimentally measured life.

Characterization of Phase Evolution

As solder joints shrink in size in response to general microelectronics drive towards transistor scaling, the intermetallic phases occupy a significant fraction of the solder joint. They are also prone to new failure mechanisms as illustrated in Figure 4. The existing characterization procedures rely on cross-sectioning and polishing to observe the phases, which interrupts the test. So, an ongoing effort is to design nanoprobe Blech-like test structures using semiconductor fabrication techniques that allow in situ visualization of phase evolution without cross-sectioning (Figure 8, [25]).

Figure 8 Semiconductor fabrication techniques were used to create test structures for thermal aging and electromigration characterization of microscale solder joints [25]. Tests can be carried out inside an SEM using the four-probe test setup.

Related modeling work includes phase field multi-physics simulations of Cu6Sn5 intermetallic growth in micron-scale solder joints (Figure 9, [26]) and phase field electromigration simulations [2].

Figure 9 Phase field simulations of the dynamics and stability of Cu6Sn5 interfacial growth [26]. The model accounted for bulk diffusion of copper through the IMC, reaction at the interface followed by surface diffusion of Cu6Sn5 to yield the scalloped interfacial morphology. (a) scalloped morphology of Cu6Sn5 in as reflowed solder joints (b) problem domain and boundary conditions (c) initial morphology and (d) predicted scalloped morphology resulting from the competition between bulk diffusion of Cu, rate of reaction at the interface and the surface diffusion of IMC.

Thin Film Buckling and Wrinkling

Thin films are common to microelectronic devices, and during processing, they are prone to wrinkling and debonding. Characterizing the adhesion energy of the films is a significant challenge at the present time. Subbarayan and his students were the first to evaluate thermally induced buckling leading to debond in thin films (Figure 10) to estimate the interfacial fracture toughness [28]. They characterized the bond between aluminum and Su8 films and estimated the yield strength of the aluminum film using these experiments [29].

Particulate composites are common in many applications including in thermal interface materials used to remove heat from the chips in microelectronic packages. The interface materials are engineered to maximally pack particles in a silicone or epoxy matrix to enhance thermal conductivity while allowing a compliant system that mechanically decouples the chip from the heat spreader. The near-percolation thermal conduction in these systems often requires explicit simulations of microstructural ensembles. Simulations that don’t explicitly account for particle-particle interactions typically don’t agree well with experiments when the volume loading of the particles exceeds 30%. Subbarayan and his students developed full-filed NURBS-based models of realistic three-dimensional microstructures to identify the near percolation behavior in these systems [30]. The models were validated against experimental measurements of the effective thermal conductivity as well as classical analytical models (Figure 11).

Figure 11 Comparison between explicit, NURBS-based isogeometric simulations (of effective thermal conductivity) of particle filled thermal interface materials against experimental measurements as well as classical analytical models. Each data point in the plot corresponds to multiple simulations of random particle arrangements of the same volume fraction [30]. Numerical simulations matched the experimental measurements very accurately, demonstrating that particle geometric arrangements are sufficient to capture the observed behavior, without needing an (often assumed) interfacial resistance parameter.

They then developed random network models that capture the near percolation behavior to within 10% accuracy of the full-field models, while requiring three-orders of magnitude lower computational effort [31]. The models also relied on sophisticated particle packing algorithms including for ellipsoidal particles ([32], Figure 12).

Figure 12 (Left) The developed algorithms can pack 200,000 ellipsoidal particles in a few hours on a desktop computer [32].

Figure 13 (Right) NEMD simulation system of silica nanowire used to estimate size- dependent thermal conductivity [33].

Prof. Subbarayan was among the first to propose material density distribution to optimally design topologies of structures [34] popularly referred as SIMP (solid isotropic microstructure with penalization for intermediate densities). He proposed a Pareto-optimal problem for load-bearing bones (Figure 14) that traded metabolic cost measured by mass against risk of fracture as measured by stored energy. Bones have a known relationship between elastic modulus and apparent density of the form 𝐸 = 𝑐𝜌3 resulting in a topological optimization problem that was efficiently solved using variational sensitivity analysis.

Figure 14 Simulated optimal material distribution in a femur resulting from a trade-off between the objectives of mass and stored energy [34]. This study is one of the earliest examples of topology optimization as well as of using variational sensitivity analysis.

Other examples of Subbarayan’s optimization research include Pareto-optimal formulations and machine learning models for mapping color from monitors to printing devices [35, 36] (Figure 15) and applications to topological design of micro-channels for squeeze flow optimization of thermal interface materials in semiconductor chips [37]. Subbarayan was among the first to systematically evaluate the utility of machine learning models as global approximations for optimal design. He demonstrated that machine learning models are useful only if the cost of training is unimportant relative to cost of use [38,39].