New analytical model for thermomechanical responses of multi-layered structures with imperfect interfaces

A new analytical model is developed for thermomechanical responses of multi-layered structures with an arbitrary number of layers and subjected to general thermal and mechanical loading. The formulation is based on an extended Bernoulli–Euler beam theory and a slip-interface model. The former includes Poisson’s effect and covers both the plane stress and plane strain deformations, and the latter allows slipping between two adjacent layers but no jump in the normal displacement or traction. An analytical solution for a multi-layered structure under general thermomechanical loading is derived by using a new approach that first determines one interfacial shear stress and the curvature of the deformed structure. To illustrate the newly developed model, three example problems for two-, three- and five-layer structures respectively are analytically solved by directly applying the new model. In all three cases, the solutions are obtained in closed-form expressions by considering both temperature changes and mechanical loads including body forces, distributed normal and shear stresses on the top and bottom surfaces, and normal forces, transverse shear forces and bending moments at the two ends, unlike existing ones. It is shown that the current solution for two-layer structures recovers an existing solution without considering Poisson’s effect and mechanical loading and the classical solution of Timoshenko for perfectly bonded bi-metal thermostats as two special cases. The closed-form solution for five-layer structures with imperfect interfaces is derived here for the first time. In addition, numerical results are provided for five- and seven-layer transistor stacks to quantitatively demonstrate the new model. It is found that the current results for the five-layer transistor stack agree well with those obtained by others, thereby further validating the new model.


Introduction
Multi-layered structures have been widely used in electronics, aerospace, construction, manufacturing, and power generation industries, which are often subjected to thermomechanical loads.Various models have been developed to describe responses of such structures, including laminate theories (e.g., [29,32]), layer-wise analyses (e.g., [1,4,24]), and shear-lag models (e.g., [8,14,26,40]).These models are typically based on the assumption of perfectly-bonded or imperfect interfaces between adjacent layers.For the former, both the traction and displacement are taken to be continuous across each interface (e.g., [48]), while for the latter the continuity conditions are imposed only for the traction and the displacement field can be discontinuous across the same interface (e.g., [9]).Imperfect interface-based models are physically more representative, but they suffer from the following major drawbacks: (1) boundary value problems (BVPs) formulated using these models depend on the local equilibrium for each layer, global equilibrium of the entire multilayer structure, and continuity (compatibility) conditions at interfaces; (2) such a BVP is governed by a reduced system of coupled differential equations, whose orders increase with the number of layers.As a result, analytical solutions based on these imperfect interface models have been limited to structures with only two or three layers (e.g., [5, 6, 13, 15-18, 20, 22, 23, 25, 34, 35, 38, 42, 43, 46, 47]).
In the current study, a new analytical model is developed for thermomechanical responses of multi-layered structures with an arbitrary number of layers by using an extended Bernoulli-Euler beam theory and a slipinterface model.The former accounts for Poisson's effect and covers both the plane stress and plane strain deformation modes, and the latter allows slipping between two adjacent layers (i.e., a tangential displacement jump) but no jump in the normal displacement or traction.
The rest of this paper proceeds as follows.In Sect.2, the new analytical model is formulated, and a general solution is derived for a multi-layered structure, which is enabled by using a new approach proposed here that first determines one interfacial shear stress along with the curvature of the deformed structure.To illustrate the newly developed model, three example problems for two-, three-and five-layer structures respectively are analytically solved in Sect. 3 by directly applying the new model.In all three cases, the solutions are obtained in closed-form expressions.The current solution for two-layer structures recovers those of Liu and Chen [22] for bi-layer beams with an imperfect interface and Timoshenko [39] for perfectly bonded bi-metal thermostats as two special cases.The closed-form solution for five-layer structures with imperfect interfaces is obtained here for the first time.In this section, numerical results are also provided for five-and seven-layer transistor stacks and compared with the data from Pao and Eisele [30] for a five-layer transistor stack to verify and demonstrate the new model.A summary is presented in Sect. 4.

New model for a multi-layered structure with interfacial slipping
Consider a multi-layered structure of a uniform rectangular cross-section, which consists of k homogeneous, isotropic, linear elastic thin layers of uniform thickness, as shown in Fig. 1.The k bonded layers are made from dissimilar materials, and the k − 1 interfaces are imperfect (weak), allowing adjacent layers to slip over each other while staying in contact all the time.Each layer satisfies an extended Bernoulli-Euler beam theory that incorporates Poisson's effect and covers both the plane stress and plane strain deformation modes.This type of structure can well represent electronic packaging assemblies (e.g., [30,38,41]) and composite beam-columns with partial interaction (e.g., [18,21,27]).
The multi-layered structure undergoes axial and flexural deformations due to the mechanical loads applied at its boundaries and a temperature change T T -T 0 , where T 0 is a reference temperature.In Fig. 1, a global coordinate system {x, y, z}, with its origin placed at the center of the mid-plane of the first layer from the bottom, is adopted to describe the entire structure, while a local coordinate system {x, y, z m } with its origin located at the center of the mid-plane of the layer is used for each layer, where m denotes the mth layer.The choice of this global coordinate system is not unique, since the origin for the global coordinate system can be placed at any point on the symmetry axis z of the cross section (e.g., [22,23,44]) without affecting the determination of the curvature of the deformed structure from simultaneously enforcing the equilibrium Fig. 1 a Multi-layered structure with a uniform rectangular cross-section subjected to thermomechanical loads; b deformation of a differential element (dx) of the mth layer equations for forces in the x-and z-directions and moments in the y-direction.Note that the z-or z m -axis is also the symmetry axis in the width direction.In addition, z m zd m , where d m is the distance between the centroidal axis of the mth layer and the centroidal axis of the first layer, as shown in Fig. 1.

Kinematic relations
In the current study, each layer is regarded as a Bernoulli-Euler beam satisfying the following kinematic relations (see Fig. 1b): where u (m)  x , u (m)   y and u (m)   z are, respectively, the x-, y-and z-components of the displacement vector u of a point (x, 0, z m ) on the vertical symmetry plane of the mth layer, and u (m)  0 and w (m) are, respectively, the axial and vertical displacement components of the corresponding point (x, 0, 0) on the centroidal axis of the mth layer, the superscript/subscript m indicates the mth layer, and () , x ≡ d()/dx.
From Eq. ( 1), the strain components in each layer can be readily obtained as where are the axial strain on the centroidal axis and the curvature of the mth layer, respectively.

Slip-interface model
Each of the k − 1 interfaces is taken to be imperfect, which allows slipping between the two adjacent layers (i.e., a tangential displacement jump in the x-direction) but no jump in the normal (z-direction) displacement or traction.This type of slip-interface model has been widely used in designing composite beam-columns in civil engineering (e.g., [18,21,22,27]).It can be viewed as a special case of the spring-layer model for imperfect interfaces known in the mechanics community (e.g., [2,7,9,19,40]).
According to the slip-interface model (e.g., [9,18,23]), where τ xz is the shear stress acting on the mth interface (between layers m and m + 1) (see Fig. 2a), t m and t m+1 are, respectively, the thickness of the mth and (m + 1)th layers, and K m is the stiffness constant of the mth interface (with a unit of Pa/m).Clearly, Eq. (3) indicates that K m → ∞ when the two layers are perfectly bonded (and hence the interface is perfect), and K m 0 when the two layers are completely separated.Using Eq. (1) in Eq. (3) gives Note that the slip-interface model with u (m)   x 0 (while u (m) z 0) adopted here differs from those implemented for perfectly-bonded interfaces, where u (m)   x u (m)   z 0 (e.g., [4,25]), or spring-layer imperfect interfaces, where the interfacial traction vector is related to the jump in the displacement vector u (m) through the elastic stiffness of the interface (treated as a spring layer) (e.g., [7,19,40]).Based on this slip-interface model, w (m) (x) w(x), ϕ (m) (x) ϕ(x) w(x) ,x x , ( 5 ) where w and ϕ are, respectively, the deflection and curvature of the structure at the x cross-section, which are independent of z m or z.This z-independence results from the kinematic relations in the Bernoulli-Euler beam theory (see Eq. ( 1)) and the continuity condition of no jump in the vertical displacement across each interface, which is embodied in the slip-interface model.

Constitutive equations
The material of each layer is assumed to be isotropic.Then, applying Hooke's law yields the axial normal stress component in each layer as, upon using Eq.(2a), where T is the temperature change that is the same for all layers and uniform along the x-direction, and c (m) 11 and d 11 are, respectively, the elastic stiffness constants and thermal expansion coefficients given by (e.g., [3,35]) for plane stress deformations, and for plane strain deformations.In Eqs.(7a) and (7b), E (m) , ν (m) and α (m) are, respectively, Young's modulus, Poisson's ratio and the coefficient of thermal expansion of the material of the mth layer.Note that Eq. (7b) can be either directly derived from the plane strain conditions or simply obtained from Eq. (7a) using the well-known transformation relations for E (m) , ν (m) and α (m) between the plane stress and plane strain cases (e.g., [33]).When Poisson's effect is suppressed by setting ν (m) 0, Eqs.(7a) and (7b) both reduce to Equation (7c) has been used in existing studies based on the Bernoulli-Euler beam theory (e.g., [8,22,23]), which did not consider the Poisson effect.From Eqs. (2a) and ( 6), the resultant normal force N in the mth (with m ∈ {1, 2, …, k}) layer can be obtained as where with b being the width of the beam structure, which is uniform along the beam length in the x-direction.

Equilibrium analysis
Consider the equilibrium of a differential element of the multi-layered structure, whose length is dx. Figure 2a shows the free body diagram (FBD) of the mth layer in the differential element, with the normal and shear stresses on the two interfaces indicated.The FBD for the differential element dx of the structure without separating the layers is displayed in Fig. 2b.From the FBD in Fig. 2a, the equilibrium equations of the mth layer can be readily obtained as, upon enforcing the force balances in the x-and z m -directions and the moment balance in the y-direction, x and M (m) x are, respectively, the resultant normal force, transverse shear force and bending moment on the x cross-section of the mth layer, q (m) x and q (m) z are, respectively, the xand z m -components of the body force acting in the mth layer (with a unit of N/m 2 ), σ for any x ∈ (− L/2, L/2).Hence, Eqs.(10a) and (11) will be used as the two equilibrium equations for the mth layer in terms of the resultants N (m) x and M (m) x .The boundary conditions (BCs) for the mth layer (associated with Eqs.(10a-c)) are given by where the overhead bar represents the prescribed value and use has been made of Eq. ( 5).These boundary conditions are required for the determination of N (m) x and M (m) x from the equilibrium equations listed in Eqs.(10a) and (11), which is demonstrated in Sect.3.
Similarly, from the FBD in Fig. 2b, the equilibrium equations for the entire multi-layered structure can be obtained as xz are, respectively, the normal and shear stresses applied on the top and bottom surfaces of the structure.In reaching Eq. (13c), the moment is taken about the y-centroidal axis of the first layer on the (x + dx) cross-section.
Note that Eqs.(13b) and (13c) can be combined to give From Fig. 3, it follows that which are obtained from the equivalence of the normal force, shear force and bending moment on the x cross-section.Note that the values of d m in Eqs.(13c), ( 15) and (16c) (see Fig. 1a) are given by In addition, the moment balance leading to Eq. ( 16c) is about the y-centroidal axis of the first layer on the x cross-section.From Eqs. (10a,b), (11), (14a,b) and (16a-c), it can be readily shown that satisfying the local equilibrium equations in Eqs.(10a,b) and (11) ensures that the global equilibrium equations listed in Eqs.(13a,b) and (15) are automatically satisfied for the current multi-layered structure under the general thermal and mechanical loads prescribed.
The boundary conditions for the entire structure (associated with Eqs.(13a-c)) read xz V or w (1) w at x ±L/2, Mx M or w (1)   ,x w ,x at x ±L/2, (18a-c) where
Substituting Eq. ( 27) into the (k -1)th (the last) equation in Eq. ( 26) leads to the following 2(k − 1)th order ODE involving τ (1) xz, x and ϕ ,xx : (28) where the first five constant coefficients in each of the A i and B i groups are listed in Eqs.(A1) and (A2) in "Appendix", and the inhomogeneous term H 1 is given by with Note that A i and B i in Eq. ( 28) depend on the number of layers and the geometrical and material properties of the layers, and H 1 depends on the mechanical loads q (m) x , τ (0)  xz and τ (k)  xz additionally, as shown in Eqs. ( 25), (29a), (29b), (A1) and (A2).
Next, using Eq. ( 27) in Eq. ( 23) results in the following 2(k − 2)th order ODE involving τ (1)  xz, x and ϕ ,xx : xz,x dx 2 + C 0 τ (1)   xz,x where the first four constant coefficients in each of the C i and D i groups are listed in Eqs.(A3) and (A4) in "Appendix", and the inhomogeneous term H 2 is given by Similarly, C i and D i in Eq. ( 30) depend on the number of layers and the geometrical and material properties of the layers, and H 2 depends on the mechanical loads q zz additionally, as shown in Eqs. ( 25), (29b), ( 31), (A3) and (A4).
Solving Eqs. ( 28) and (30) simultaneously will lead to the determination of τ (1)  xz, x and ϕ ,xx .For the case with H 1 andH 2 being constants, the general solution of this equation system has the form: where λ n is the nth root of the following 2(k -1)th-degree polynomial characteristic equation: n is a constant depending on λ n , and C n are constants to be determined from the following 2(k − 1) boundary conditions: , which are obtained directly from Eqs. ( 19) and (12a,c), with N (m) and M (m) being the axial force and bending moment acting at the two ends of the mth (m ∈ {1, 2, …, k}) layer.Using Eqs.(32a,b) and ( 27) in Eq. ( 34) will lead to a set of 2(k -1) algebraic equations, which can be solved to obtain the constants C n .From Eqs. ( 27) and (32a), the interfacial stresses τ (m)  xz (m ∈ {1, 2, …, k -1}) can be determined by directly integrating τ (m)  xz, x , which gives where x x and bending moments M (m) x x can be acquired from the equilibrium equations in Eqs.(10a) and (11) as where From Eqs. ( 5), (8b), ( 9) and (21b), it follows that With τ xz, x obtained in Eqs. ( 27) and (32a) and ϕ , x x given in Eq. (32b), the interfacial normal stresses σ zz (m ∈ {1, 2, …, k -1}) can be determined sequentially from Eq. ( 37), starting from m 1 through m k -1, with σ (0) zz prescribed on the bottom surface of the structure (see Figs. 2b and 3).
The curvature ϕ can be readily obtained from Eq. (32b) as, after integrating it twice, where J 1 and J 2 are two integration constants to be determined by applying the boundary conditions of Mx M at x ±L/2 in Eq. (18c), with Mx obtained from Eqs. ( 5), (8b), ( 9) and (16c) as

Examples
To illustrate the new analytical model developed in Sect. 2 for multi-layered structures with an arbitrary number of layers and subjected to general thermal and mechanical loading, three example problems for two-, threeand five-layer structures respectively are solved here by directly applying the new model.In all three cases, the solutions are obtained in closed-form expressions.For the two-layer case, the current solution recovers the one provided by Liu and Chen [22] without considering Poisson's effect and any mechanical load as a special case and reduces to the classical solution of Timoshenko [39] for bi-metal thermostats with perfect bonding and subjected to a uniform temperature change.For the five-layer case, the closed-form solution is derived here for the first time, which can be directly used to study electronic packaging assemblies and composite joints with imperfect interfaces (e.g., [8,30,44]).In this case, the number of layers k 2, and there is only one interface.From Eq. ( 26), it follows that τ (1)   xz,x x x β (1) 0 τ (1)  xz,x + (1) ϕ ,x x + η (1)  x + β (1) where xz and τ (2 xz are the shear stresses applied on the bottom and top surfaces of the structure, as shown in Fig. 4, and x K 1 c (1) which are directly obtained from Eq. ( 25).With k 2, Eq. ( 23) becomes (1) Substituting Eq. ( 42) into Eq.(40) gives xz,x 12 (1)   c (1) xz,x + η (1)  x + β (1) This is a second-order inhomogeneous ODE for τ xz, x with constant coefficients.When the body force components (i.e., q (m) x and q zz ) and the first derivative of the distributed shear stress (i.e., τ xz, x ) prescribed on the top and bottom surfaces are constants, the inhomogeneous term on the right-hand side of Eq. ( 43) becomes a constant upon using Eq. ( 41).As a result, the general solution of Eq. ( 43) can be readily obtained as τ (1)   xz,x a 1 e γ 1 x + a 2 e γ 2 x − 12 (1) qz − σ where a 1 and a 2 are constants, and γ 1 and γ 2 are the roots of the characteristic equation (a quadratic one) of the homogeneous part of Eq. ( 43) given by in which use has been made of Eq. ( 41).
For the special case with ν (m) 0 (m ∈ {1, 2}), q 0 and M (2) 0, Eq. (52) simplifies to, with the help of Eq. (7c), The curvature of the deformed two-layer structure obtained in Eq. (53a) is identical to that derived in Liu and Chen [22] by using a different approach for the special case without considering Poisson's effect, interfacial normal stress, body forces, distributed normal and shear stresses prescribed on the top and bottom surfaces, and normal forces, transverse shear forces and bending moments applied at the two ends.This provides a validation of the current new model.In addition, using Eq.(53a) in Eq. ( 6) gives, with the help of Eqs.(2a), ( 5), (7c), (8a) and (12a), the axial normal stress in layer 2 in this special case as A comparison of Eq. (53b) with Eq. (1) in Noyan et al. [28], which is based on a shear-lag model, shows that the functional relationship between σ (2) x x and x given by the current solution for the film (top layer) with two free-free ends is the same as that presented in Noyan et al. [28], thereby further validating the new model.
Equation (53a) explicitly shows that the bending deformation (as measured by the curvature ϕ) for the special case with only a temperature change (but without any mechanical load) is purely due to the mismatch of the coefficients of thermal expansion of the two layers that are imperfectly bonded with 0 < K 1 < ∞.If the two layers are completely separated, then K 1 0 and hence ϕ 0 according to Eqs. ( 45) and (53a).That is, there will be no bending deformation induced by the temperature change if the two layers are not bonded at all.However, this is no longer the case when a mechanical load is applied in addition to the temperature change.It is seen from Eq. ( 52) that with any of qz , σ xz, x , M (1) and M (2) being non-zero, ϕ 0 even if K 1 0. When the two layers are perfectly bonded, K 1 → ∞ and Eq. ( 52) reduces to 11 11 11 t 3  1 +c (2) b c (1) which is the curvature of the perfectly bonded two-layer structure under the general thermal and mechanical loading.If ν (m) 0 (m ∈ {1, 2}), q 0 and M (2) 0, Eq. (54) further simplifies to, upon using Eq.(7c), which is the same as that initially derived by Timoshenko [39] for perfectly bonded bi-metal thermostats subjected to a uniform temperature change (see also Zhang and Xing [45] for an alternative derivation of Eq. ( 55)).This indicates that the classical constant curvature formula for perfectly bonded bi-layer structures first obtained by Timoshenko [39] is recovered by the current two-layer solution as a special case.
In addition, for perfectly bonded two-layer structures subjected to both a temperature change and applied mechanical loads, Eq. (54) shows that there will always be bending deformations (with ϕ 0) even if there is no temperature change (i.e., T 0) or no mismatch in the CTE for the two layers (i.e., α (1) − α (2) 0).This differs from what is exhibited by two-layer structures subjected to a temperature change only, for which the curvature is given by Eq. (55), which is a special case of the current formula derived in Eq. (54) for perfectly bonded two-layer structures under general thermal and mechanical loading.

Three-layer structure
Consider a three-layer structure subjected to a uniform temperature change T and mechanical loads on the top and bottom surfaces and at the two ends x ±L/2, as shown in Fig. 5.

Five-layer structure
Consider a five-layer structure subjected to a uniform temperature change T and mechanical loads on the top and bottom surfaces and at the two ends x ±L/2, as shown in Fig. 6.
When d < 0, Eq. ( 92) has one real root and two conjugated complex roots that can be obtained using the Cardano formula, which give the three roots of Eq. (90) as, upon using Eq.(93a), where χ 1 in Eq. (96a) is the real root, which will be used in Eq. ( 89) to obtain ω 2 j (j ∈ {1, 2, 3, 4}), and thereby the eight roots ω i of Eq. ( 85).
When d 0, Eq. ( 92) has a triple root of 0 if 0 or a single root of ξ 1 3 and a double root of Then, it follows from these results and Eq.(93a) that the roots of Eq. ( 90) in this case are given by Using any of the three values of χ listed in Eq. (97a) or (97b) in Eq. (89) will yield ω 2 j (j ∈ {1, 2, 3, 4}) and thus the eight roots ω i of Eq. (85).
The rest of the quantities, including N (m) x , M (m) x and u (m) 0 (m ∈ {1, 2, 3, 4, 5}), for the current five-layer structure can be determined by following procedures similar to those used in Sect.3.2 for the three-layer structure.

Numerical Results
To quantitatively illustrate the new analytical model developed in Sect. 2 and the closed-form solutions obtained above in this section, two transistor stacks subjected to thermomechanical loading are considered as numerical examples in this sub-section.
A typical industrial transistor stack has a multi-layer structure containing a semiconductor (e.g., silicon (Si)) layer, a thermal barrier/spreader (e.g., beryllium oxide (BeO)), and a metal substrate (e.g., aluminum (Al) or copper (Cu)), which are connected through eutectic solder joints (e.g., lead/tin alloy (Pb/Sn)) (e.g., [30]).Two transistor stacks are analyzed here, one of which is five-layer and the other is seven-layer, as shown in Figs.7a and 7b.The difference is that only one Pb/Sn joint is used in the former, while three Pb/Sn layers are employed in the latter.A thermomechanical analysis has been conducted for the five-layer transistor stack displayed in Fig. 7a by Pao and Eisele [30], who extended the bi-metal thermostat model (e.g., [31,36,37,39]) to determine interfacial shear and peel (normal) stresses in five-layer joints and transistors.Pao and Eisele [30] obtained a system of four coupled second-order ODEs for the interfacial stresses, which was converted to an eigenvalue problem and subsequently solved numerically.Their analysis is limited to five-layer structures and confined to cases with temperature changes (through uniform heating or cooling) but without any mechanical load, which is different from the current model that can be applied to structures with an arbitrary number of layers and subjected to both mechanical loads and temperature changes.
In order to directly compare with the results reported in Pao and Eisele [30], the geometrical and material parameters used by them are adopted for the five-layer transistor stack considered in the current numerical example.The same material properties and geometrical dimensions are also employed for the seven-layer transistor stack which contains two more Pb/Sn layers.The parameter values for both the five-and seven-layer transistor stacks are listed in Table 1.
As a first approximation, the stiffness constant K m of the mth interface between layers m and m + 1 is taken to be the average of the shear modulus values of the two adjacent layers over the average thickness of the two layers (i.e., t (m) i (t m + t m+1 )/2) given by Note that this expression remains the same for both the plane stress and plane strain cases.Equation ( 106) is used to obtain the interfacial stiffness constants K m for the five-layer (with m ∈ {1, 2, 3, 4}) and seven-layer (with m ∈ {1, 2, 3, 4, 5, 6}) stacks considered here from the values of Young's modulus E (m) , Poisson's ratio ν (m) and layer thickness t m listed in Table 1.
In addition, the material constants c 11 for each layer in the two example problems are determined from the values of E (m) , ν (m) and α (m) given in Table 1 using Eq.(7a) for the plane stress case.

Model verification
The newly developed analytical model is verified against the results of Pao and Eisele [30] for a five-layer transistor stack consisting of a Si chip, a BeO thermal spreader, an eutectic solder (lead/tin alloy (Pb/Sn)), a Cu interlayer, and an Al substrate, as shown in Fig. 7a.This five-layer stack is subjected to a temperature change of T −65 • C (uniform cooling) without any applied mechanical load (i.e., σ 0).In addition, the body forces are all neglected (i.e., q ) and the axial force, transverse shear force and bending moment at the two ends of each layer are taken to be zero-valued (i.e., N xz (m ∈ {1, 2, 3, 4}) varying with x (the distance from the midspan) obtained using Eqs.(99a) and (101a-c) are plotted in Fig. 7c, where the results presented in Pao and Eisele [30] are also displayed for comparison.It is seen from Fig. 7c that the interfacial shear stresses τ (m) xz (m ∈ {1, 2, 3, 4}) predicted by the current model agree well with those provided in Pao and Eisele [30], which verifies and supports the newly developed model.

Transistor stacks under thermal and mechanical loads
To investigate the responses of the five-layer transistor stack under both thermal and mechanical loads, the stack is subjected to a uniform temperature rise of T +65 • C and a mechanical load σ ( 5) zz −20MPa (a uniform compression) applied on the top surface.The other mechanical loads are taken to be zero-valued, i.e., Figure 8 shows the mechanical responses of the five-layer transistor stack under the prescribed thermomechanical loading.The numerical results displayed in Fig. 8 are obtained using Eqs.(99a) and (101a-c) for τ zz both change sharply near the two ends of the stack, indicating a strong edge (boundary) effect.In addition, the interfacial shear stress at the BeO-Si interface, i.e., τ (4) xz , is found to be the lowest due to the small CTE mismatch between the BeO and Si layers, while the highest interfacial shear stress (excluding the regions with the edge effect) is developed at the BeO-Pb/Sn interface, i.e., τ (3) xz , because of the large CTE mismatch between the BeO and Pb/Sn layers.Furthermore, the interfacial normal stress at the BeO-Si interface, i.e., σ (4) zz , is seen to be the smallest, while the interfacial normal stress at the Cu-Al interface, i.e., σ (1) zz , is the largest.From Figs. 8d and 8e, it is observed that the BeO layer (thermal spreader) exhibits the smallest axial strain ε (4) 0 and axial displacement u (4) 0 among all layers.This is due to the low CTE and high elastic modulus of BeO, as shown in Table 1.In contrast, the highest axial strain ε  Figures 8c-8f show that all layers are bent due to the applied mechanical load σ (5)  zz and the large CTE mismatch between the bottom three layers (Al, Cu, and Pb/Sn) and the top two layers (BeO and Si).In addition, Fig. 8e shows that each of the Al, Cu, Pb/Sn and BeO layers exhibits an axial stretch, while the Si layer displays an axial compression.The transition from the stretch to the compression indicates the existence of a neutral surface, which is located in the solder (Pb/Sn) layer where both the bending moment and axial normal force are seen to be almost zero from Figs. 8c and 8f.
To further understand how the curvature ϕ and deflection w of the transistor stack vary with applied mechanical loads, ϕ and w of the five-layer stack subjected to the temperature rise T +65 • C and applied mechanical load σ (5)  zz of different magnitudes are plotted in Figs.9a and 9b, respectively.The numerical results shown in Figs.9a and 9b are obtained using Eqs.( 5) and (105), with σ 0 (m ∈ {1, 2 , 3, 4, 5}) at x ±L/2.It is clear from Figs. 9a and 9b that both the curvature ϕ and deflection w (magnitude) increase with the increase of the applied load σ (5)  zz (in magnitude).zz of different magnitudes; c curvature ϕ and d deflection w of the seven-layer transistor stack induced by the temperature rise T +65 • C and applied mechanical load σ (7)  zz of different magnitudes Finally, to demonstrate the application of the newly developed model for stacks with more than five layers, the seven-layer transistor stack shown in Fig. 7b is considered.In this transistor stack, the Pb/Sn solder joint is used between every two adjacent layers, as indicated in Fig. 7b.The transistor stack is subjected to a temperature rise of T +65 • C and a mechanical load σ (7)  zz varying from 0 to −20 MPa on the top surface.The curvature and deflection of this seven-layer transistor stack are plotted in Figs.9c, and 9d, respectively.The numerical values displayed in these figures are obtained using Eqs.( 5) and (38) in the general analytical model, with σ It is observed from Figs. 9c and 9d that both ϕ and w (magnitude) of the seven-layer stack increase with the increase of the applied mechanical load σ (7)  zz (in magnitude), which is similar to what is exhibited by the five-layer transistor stack (see Figs. 9a and 9b).However, compared to the five-layer stack, the overall stiffness of the seven-layer stack is enhanced due to the inclusion of two additional solder (Pb/Sn) layers, and, as a result, the curvature ϕ and deflection w (magnitude) of the latter are lower than those of the former under the same thermal and mechanical loads, as shown in Fig. 9.

Summary
A new analytical model for multi-layered structures under general thermomechanical loading is developed using an extended Bernoulli-Euler beam theory and a slip-interface model.It accounts for Poisson's effect and allows for slipping between two adjacent layers.An analytical solution for a multi-layered structure with an arbitrary number of layers and subjected to a temperature change and mechanical loads on the top and bottom surfaces and at the two ends is obtained here for the first time by following a new approach that first determines one interfacial shear stress along with the curvature of the deformed structure.
Three example problems for two-, three-and five-layer structures respectively are analytically solved by directly applying the new model.Closed-form solutions are derived for all three problems, each of which includes Poisson's effect and considers a uniform temperature change and general mechanical loads, unlike existing solutions.The current solution for two-layer structures recovers the solution obtained by Liu and Chen [22] without considering Poisson's effect and mechanical loading and the classical solution of Timoshenko [39] for perfectly bonded bi-metal thermostats as two special cases.This verifies the new model.In addition, the closed-form solution for five-layer structures with imperfect interfaces is derived here for the first time, which can find direct applications in designing electronic packaging assemblies and composite joints.Moreover, all expressions of the general and specific solutions are derived in terms of physical parameters and variables, which can be easily adjusted to fit different problems.
To quantitatively illustrate the new model and current analytical solutions, numerical results are presented for physical problems of five-and seven-layer transistor stacks subjected to prescribed thermomechanical loading.The interfacial shear stresses predicted by the current model for the five-layer transistor stack are seen to agree very well with those obtained by Pao and Eisele [30] using an extended bi-metal thermostat model, which further validates the newly developed analytical model.The general and specific solutions obtained in the current study can be applied to analyze multi-layered stack structures with imperfect interfaces that are commonly used in electronic assemblies, composite beam-columns, and all-solid-state batteries.
Finally, it should be mentioned that the current formulation is limited to multi-layered structures with a uniform temperature field.It can be extended to cases with non-uniform temperature distributions and/or creep deformations, which is currently being pursued.

Fig. 2
Fig. 2 FBDs of the differential element dx of the structure: a mth layer; b entire structure xz are, respectively, the normal and shear stress components on the mth interface, and σ (m−1) zz and τ (m−1) xz are, respectively, the normal and shear stress components on the (m -1)th interface.From Eqs. (8a) and (10a), it is clearly seen that the axial force N (m) x and normal stress σ (m) x x are directly linked to the interfacial bonding constraints through the interfacial shear stresses τ (m) xz and τ (m−1) xz .Note that Eqs.(10b) and (10c) can be combined to yield

Fig. 3
Fig. 3 Equivalence of the resultant normal force, shear force and bending moment between local and global systems

Fig. 9 a
Fig. 9 a Curvature ϕ and b deflection w of the five-layer transistor stack induced by the temperature rise T +65 • C and applied mechanical load σ (5) 7}) at x ±L/2.

Table 1
Material and geometrical parameters for the five-layer and seven-layer transistors