- Email: [email protected]
rigid substrate system
Mechanics Research Communications 38 (2011) 283–287
Contents lists available at ScienceDirect
Mechanics Research Communications journal homepage: www.elsevier.com/locate/mechrescom
Elastic thermal stresses in a circular overlay/rigid substrate system Z.F. Yuan, H.M. Yin ∗ Department of Civil Engineering and Engineering Mechanics, Columbia University, 610 Seeley W. Mudd, 500 West 120th Street, New York, NY 10027, USA
a r t i c l e
i n f o
Article history: Received 1 November 2010 Received in revised form 17 February 2011 Keywords: Thermal stresses Layered materials Thermoelasticity Axisymmetry Thin film
a b s t r a c t When an overlay/substrate system is subjected to a temperature change, thermal stress is induced due to the difference of thermal expansion coefficients between the overlay and substrate. This paper studies elastic thermal stress distribution of a circular overlay bonded to a rigid substrate due to temperature variation. Under the plane assumption, the radial displacement in the overlay is obtained in a series form. Using a weak form boundary condition, one can obtain an approximate solution in closed form for the thermoelastic field. This approximate solution gives a description of displacement along both radial and thickness directions. Elastic thermal stress distributions are demonstrated and compared with the finite element results. An asymptotic analysis of the solution is given when the thickness of a thin film is much less than its radius. The closed-form solution is useful for stress-strain analyses and design of circular layered systems such as protective coatings and microelectromechanical components. Application of this solution in modeling spiral cracking of circular thin film/substrate is underway. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction Overlay/substrate systems have been widely used in engineering structures (Timm et al., 2003; Leung and Tung, 2006; Yin et al., 2007a), electronic packaging (Suhir, 1994; Yin et al., 2008), surface coatings (Shevchuk and Silberschmidt, 2005; Yin et al., 2007b), and material testing and modeling components (Agrawal and Raj, 1989; Hutchinson and Suo, 1992; Porter and Hills, 2001). When these systems are subjected to an ambient temperature change or a service temperature load, thermal stress will be induced due to the different coefficients of thermal expansion (CTEs) between the overlay and substrate, which may cause thermal cracking within the layers or delamination along the interface. Advances in design and maintenance of these material systems require the quantitative understanding of the thermal stress distribution as part of the structural analysis process. The crack within the thin film has different patterns, some of which exhibit a very interesting phenomenon. Yin et al. (2007a,b) studied a parallel thermal cracking of asphalt overlay bonded to rigid substrate, and the cracking is nearly uniformly distributed. The solution was extended to general materials (Yin et al., 2008) and used for investigation of fracture initiation and saturation (Yin, 2010a,b). Spiral cracking was observed in the drying process of thin films by Leung et al. (2001) and Neda et al. (2002). Sendova and Willis (2003) presented four different crack patterns in sili-
∗ Corresponding author. Tel.: +1 212 851 1648; fax: +1 212 854 6267. E-mail address: [email protected] (H.M. Yin). 0093-6413/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2011.04.005
cate sol–gel films such as straight-line, curved closed-loop, crescent (sawtooth) periodic, and spiral cracks. Crack patterns in thin films may change with the material stiffness and fracture toughness of the materials, interface condition, and loading conditions. Stress and strain analyses will be valuable to understand the fracture mechanisms of an overly/substrate system and predict its cracking behavior under a certain working condition. Some theoretical research has been established for cracking propagation. Suo (1989) developed the solutions to singularities such as point force, point moment, edge dislocation and transformation strain spot in bonded elastic blocks of dissimilar materials. Thouless (1995) derived an approximate analytical solution for the minimum spacing that can exist between a series of parallel cracks propagating in a thin film and investigated mechanism of development and relaxation of stresses in films. Although fracture mechanics based modeling can provide local mechanical fields with a high fidelity, it is not straightforward to formulate thermal stress field in the whole system. The first theoretical formula for evaluation of thermal stresses in an overlay/substrate system was proposed by Stoney (1909). It has been widely used to predict stresses of thin films with the system curvatures and later extended to multilayer structures by many investigators, for example, see Freund and Suresh (2004) and Feng et al. (2008). Obviously, this school of semi-empirical approaches cannot be applied if the substrate is so stiff or so thick that no bowing occurs (Suhir, 1994). Suhir proposed the interfacial compliance concept (Suhir, 1989, 1991) and could address these special cases (Suhir, 1994). However, most existing models assumed a uniform or linear stress distribution in the thickness direction and reduced the two dimensional (2D) problem into a one dimensional (1D) problem, so that the
284
Z.F. Yuan, H.M. Yin / Mechanics Research Communications 38 (2011) 283–287
z
R
Substituting Eq. (6) into Eqs. (2)–(4), and then into Eq. (1), we obtain the governing equation as
E1 ,ν 1 , α1 , h
E0 ,ν 0 , α 0 , H
r
r 2 ur,rr + rur,r − ur 1 − 1 + ur,zz = 0 2 r2
(7)
Using the method of separation of variables, we can write the displacement in the following form as
T Fig. 1. A circular thin film boned to a rigid substrate and subjected to a temperature change T.
ur = R(r)Z(z)
stress variation in the thickness direction could not be directly studied (Suhir, 1994, 2000; Xia and Hutchinson, 2000; Timm et al., 2003). However, because the thermal stress transferred across the interface in mainly through the shear stress, which vanishes at the surface, a simple assumption of uniform or linear stress may not reflect the actual stress distribution. Yin et al. (2007a, 2008) have developed a 2D analytical model to explicitly predict the stress distribution for rectangular overlay/substrate systems. The purpose of this paper is to extend that solution to circular overlay/substrate systems. A compliant overlay bonded to a rigid substrate and subjected to uniform temperature change is considered as seen in Fig. 1, which is very common in polymer coating systems.
r 2 Rrr
(8)
Then, Eq. (7) can be rewritten as 1 − 1 Zzz + rRr − R =− = c2 2 Z r2R
where c is a constant to be determined from boundary conditions. Then we obtain Zzz +
2c 2 Z=0 1 − 1
(c 2 r 2 + 1)R − rRr − r 2 Rrr = 0
rr,r + zr,z
(1)
where is the coordinate in the circumferential direction. In the z direction, because no constraint is applied, it is assumed that zz = 0. Then the isotropic material constitutive law provides rr = =
E1 1 − 12 E1
[εrr + 1 ε − (1 + 1 )˛T ]
[1 εrr + ε − (1 + 1 )˛T ] 2
1 − 1
(2)
(3)
2 c 1 − 1
(12)
where A and B are to be determined by the boundary conditions. Because the substrate is so stiff compared with the overlay that it can be approximated as a perfectly rigid material. Therefore, the deformation of the overlay is completely constrained by the substrate. Notice that although the substrate has a thermal strain, because the thermal stresses are induced by the mismatch of CTEs between the overlay and the substrate, the effect of the thermal strain in the substrate has been included in Eqs. (2) and (3) and will not have effect on the stress analysis. Therefore, the mechanical displacement is considered to be zero. In addition, the shear stress on the surface of the overlay is zero. We can write ur (r, 0) = 0, and ur,z (r, h) = 0
(13)
Eq. (13) implies B = 0 and di =
1 2
+i
h
, with i = 0, 1, 2, ...
(14)
so that a series form of solution for Z(z) is obtained as Zi (z) = Ai sin(di z)
(15)
For each di with i = 0,1,2,.., we can obtain
ci =
1 − 1 di 2
(16)
The general solution for Eq. (11) can be written as
and zr =
(11)
The general solution for Eq. (10) is written as
2. Basic formulation
rr − =0 + r
(10)
and
Z = A sin(dz) + B cos(dz), with d =
Consider a circular overlay with thickness h, Young’s modulus E1 , Poisson’s ratio v1 , the coefficient of thermal expansion (CTE) ˛1 and radius R fully bonded to a rigid substrate with thickness H, Young’s modulus E0 E1 , Poisson’s ratio v0 , and CTE ˛0 . The overlay/substrate system is subjected to an ambient temperature change T, as illustrated in Fig. 1. Because the overlay and the substrate have different CTEs, whose difference is defined as ˛ ¯ = ˛1 − ˛0 , a thermal stress is induced in the overlay. Assume the material is thermoelastically isotropic. The thermal stresses can be solved by an axi-symmetric formulation. The coordinate r is set up along the interface between two layers, and the coordinate z is along the axisymmetric axis in the thickness direction. The equilibrium equation in r direction is written
(9)
E1 zr 2(1 + 1 )
(4)
Ri (r) = Ci I1 (ci r) + Di K1 (ci r)
(17)
where I1 (ci r) and K1 (ci r) are the modified Bessel functions. Using the axisymmetric condition
Because the substrate provides a strong constraint for the overlay, the thermal stress is mainly induced in the plane perpendicular to the z direction. Therefore, a plane assumption is introduced that all the points in a plane normal to the z direction keep in the same plane after the deformation, i.e.:
We obtain Di = 0. Therefore, the solution for Eq. (7) can be written in a series form as
uz (r, z) = uz (z)
ur (r, z) =
(5)
Using Eq. (5), we can obtain the strain-displacement relation for this axisymmetric problem, εrr = ur,r ; ε =
ur ; zr = ur,z r
(6)
ur (0, z) = 0
(18)
∞
i=0
Ai sin(di z)I1 (ci r)
(19)
where Ai can be determined by the stress free boundary condition at r = R, i.e. rr (R, z) = 0
(20)
Z.F. Yuan, H.M. Yin / Mechanics Research Communications 38 (2011) 283–287
Substituting Eq. (6) into Eq. (2) and using Eq. (20), we obtain ur,r (R, z) +
1 ur (R, z) ¯ = (1 + 1 )˛T R
= ci I0 (ci r) −
= (1 + 1 )˛T ¯
(22)
Inserting Eq. (19) into Eq. (21) and using I1 (ci r)/r, we can write ∞
Ai sin(di z) ci I0 (ci R) −
i=0
(1 − 1 )I1 (ci R) R
Extending the integral domain to (0,2h) and using the orthogonality of trigonometric functions, we can write
Ai
(1 − 1 )I1 (ci R) ci I0 (ci R) − R
1 = h
2h
(1 + 1 )˛T ¯ sin(di z)dz (23) z=0
Simplifying Eq. (23), we can analytically solve the coefficient Ai as
Ai =
¯ 2(1 + 1 )˛T
(1/2 + i) ci I0 (ci R) − (1 − 1 )I1 (ci R)/R
(24)
Overall, the elastic solution for displacement as Eq. (19) is obtained with parameter Ai ,di and ci given in Eqs. (24), (14) and (16), respectively. Then we can solve the stress field from Eqs. (2)–(4) with (6). A series form solution may reach a higher precision to the boundary condition at the circumferential surface of the thin film, but not very straightforward to use. In addition, due to the approximation of Eq. (5), the series form solution cannot guarantee a higher accuracy of the solution, which will be revisited in Section 3. In the sense of approximation, if only the first item in Eq. (14) is adopted, we can write ur (r, z) = A sin(dz)I1 (cr); d = ; and c = 2h
1 − 1 2 2h
(25)
Using this approximation, we cannot make the boundary condition in Eq. (20) exactly satisfied at each point. However, we can set the resultant force to be zero, so that we obtain: A=
(1 + 1 )˛RT ¯ 2[cRI0 (cR) − (1 − 1 )I1 (cR)]
(26)
Then we obtain a closed form solution with Eqs. (25) and (26). In summary, the elastic fields are given as ur (r, z) = rr = =
¯ (1 + 1 )˛RT sin(dz)I1 (cr) 2[cRI0 (cR) − (1 − 1 )I1 (cR)]
¯ E1 ˛T 1 − 1
R cRI (cr) − (1 − )I (cr) 0 1 1
2r cRI0 (cR) − (1 − 1 )I1 (cR)
sin(dz) − 1
R crI (cr) + (1 − )I (cr) E1 ˛T ¯ 1 0 1 1
1 − 1
Table 1 Material properties defined in FEM simulation (Suhir, 1994).
(21) I1 (ci r)
2r cRI0 (cR) − (1 − 1 )I1 (cR)
(27)
sin(dz) − 1
285
Young’s modulus kgf/cm2 Poisson’s ratio CTE (1/◦ C)
Thin film
Substrate
28 × 103 0.4 34.8 × 10−6
1.23 × 106 0.24 0
free, and the outer boundary and bottom of the substrate are fixed. We applied a temperature change of T = − 280 ◦ C. Following we will compare the theoretical predictions and FEM results of displacement and stresses. The geometric parameters of the thin film/substrate system are as follows: radius of R = 5, and thickness h = 1, as shown in Fig. 1. Figs. 2–4 illustrate the distributions of displacement ur and normal stresses rr and along the surface of the thin film. In Fig. 2, ur decreases as r increases, which means the thin film will shrink inwards. Theoretical predictions of displacement ur through Eq. (27) excellently agree with the FEM results. However, the series form solution of Eq. (19) considerably deviates from FEM results. In Fig. 3, theoretical predictions of rr and agree with the FEM results very well at the central range. However, close to the edge, FEM results of rr will approach zero, which means it is a stress-free boundary; whereas the theoretical solution monotonically decrease due to the weak boundary condition used. On the other hand, the series form solution in Eq. (22) can satisfy the stressfree boundary condition, but the accuracy in the central range is lost. The similar tendency of stress distribution of can also be found except that it is not zero at the free edge. In Figs. 2–3, the closed form solution provides very good agreement with FEM results except the points close to the edge. Although the series form solution can satisfy the stress-free boundary condition at the edge, the overall accuracy is even worse than the closed form solution. The reason to cause this problem is that the plane assumption in Eq. (5) imposes an artificial constraint on the displacement field in the vertical direction. Enforcement of the stress boundary condition at the edge amplifies the negative effect of this assumption. Further investigation of this problem in the context of mathematical physics is underway. Fortunately, the simple closed form solution provides a very good prediction of the overall elastic field. The shearing stress along the interface is of interest in design and analysis of thin film/substrate system. Fig. 4 illustrates the shearing stress rz along the interface. The closed-form solution provides a good agreement with the FEM results except the range close to the
(28)
(29)
and zr =
¯ 2 R I1 (cr) E1 ˛T cos(dz) 2 4h cRI0 (cR) − (1 − 1 )I1 (cR)
(30)
3. Comparison with finite element simulation Because the formulation developed in the last section is still an approximate solution, verification of its accuracy is important to its application. In this section, an FEM model will be developed for verification of the analytical solution. The model contains the thin film and substrate. We use ABAQUS to run the model. The material properties of thin film and substrate are given in Table 1. Notice that we will focus on the difference of CTE between thin film and substrate, which is a resource of mismatch stress. Therefore, the CTE of substrate is purposely chosen at zero. In this model, the boundaries of the thin film except the interface of substrate are
Fig. 2. Comparison of ur at top surface of thin film with radius R = 5.
286
Z.F. Yuan, H.M. Yin / Mechanics Research Communications 38 (2011) 283–287
Fig. 3. Comparison of rr (left) and (right) at top surface of thin film with radius R = 5.
To compare the closed form solution with the FEM results in z direction, Fig. 5 illustrate the distributions of ur and rz along z direction of thin film at r = R/2. The closed form solution provides an excellent agreement with the FEM results. Notice that at the free edge, we expect a considerable difference between them and thus do not show it in the z direction. In Figs. 2–5, we find a very good fit of the closed form solution and FEM results, but a considerable difference of the series form solution from the FEM results in the central range of the overlay. Although there is a slight difference at the outer part, Eq. (28) gives a precise prediction of turning point of compressive zone and tensile zone. In detail, from Fig. 3, we find at the inner part, rr is positive which denotes a tensile stress, and at the outer part about r > 4, rr is negative which denotes a compressive stress. Both FEM and Eq. (28) predict such kind of turning point, which is very important to spiral cracking development. Overall, the FEM results agree with the closed form solution very well. 4. Discussion Fig. 4. Comparison of rz at interface of thin film with radius R = 5.
edge, which is supposed to exhibit a singular effect of the shearing stress distribution. Due to the weak form boundary condition used, the proposed solution cannot capture the singularity. However, in an average sense, it provides a very good approximation.
The proposed model provides a closed form approximate solution instead of an exact solution. The assumptions and approximations adopted may impose some limitations on the application of the formulation. In summary, the following two main assumptions are used in the derivation: first, because the stress in the z direction is free, all the points in the same plane perpendicular
Fig. 5. Comparison of ur (left) and rz (right) along z direction at r = R/2 with radius R = 5.
Z.F. Yuan, H.M. Yin / Mechanics Research Communications 38 (2011) 283–287
to z axis are assumed to keep in the same plane during deformation. Therefore, the governing equation of ur is obtained. This assumption imposes a limitation of this model that the thickness and stiffness of the substrate should be large enough to make the assumption of rigid substrate applicable and to keep the plane flatten. In addition, due to this assumption, the equilibrium in the z direction cannot be rigorously satisfied at each point. Secondly, to obtain the closed-form solution, the boundary condition along the free edge surface is approximated as a weak form boundary condition that the resultant force is zero. However, it makes the stress field around the fracture surface inaccurate. The closed form solution cannot capture the singularity of stress distribution. If the local stress around the fracture surface is emphasized, the series form solution should be used. However, it is not guaranteed that the series form solution will provide higher accuracy than the closed-form solution does because of the first assumption. Generally, the normal stress in the z direction ( zz ), which is also called peeling stress, will not be zero in the neighborhood of the free edge. However, the first assumption implies zz = 0, which further affects the accuracy of this model for the stress distribution close to the free edge. Because the closed-form solution has assured the balance of the resultant force across the thickness in the radial direction and satisfied the governing equations, and is easy to apply. In addition, the substrate is assumed to be rigid in this paper. For many applications, elastic substrates may be more realistic. In these cases, an elastic overlay fully bonded to an elastic substrate needs to be considered. 5. Conclusions A series form 2D solution is derived for a circular overlay bonded to a rigid substrate and subjected to a temperature change. Using a weak form boundary condition, a closed form solution is obtained. Compared with Suhir’s 1D solution, this solution can predict the stress variation in the thickness direction. Take an average of the stress distribution, the proposed solution can be reasonably reduced to a form which is very similar to Suhir’s solution. The explicit thermal stress solution is useful to structure design and thermoelastic analysis for circular overlay/rigid substrate systems. The assumptions and limitations of this model are discussed and analyzed. Future research efforts are emphasized on releasing some of these assumptions and combining this model with numerical methods to conduct fracture simulation.
287
Acknowledgment This work is sponsored by the National Science Foundation CMMI 0954717 and the Department of Homeland Security CU091155, whose support is gratefully acknowledged. References Agrawal, D.C., Raj, R., 1989. Measurement of the ultimate shear-strength of a metal ceramic interface. Acta Metall. 37 (4), 1265–1270. Feng, X., Huang, Y., et al., 2008. Multi-layer thin films/substrate system subjected to non-uniform misfit strains. Int. J. Solid Struct. 45 (13), 3688–3698. Freund, L.B., Suresh, S., 2004. Thin Film Materials; Stress Defect Formation and Surface Evalution. Cambridge University Press, Cambridge, UK. Hutchinson, J.W., Suo, Z., 1992. Mixed-mode cracking in layered materials. Advances Appl. Mech., vol. 29. Academic Press Inc., San Diego, pp. 63–191. Leung, C.K.Y., Tung, W.K., 2006. Three-parameter model for debonding of FRP plate from concrete substrate. J. Eng. Mech. 132 (5), 509–518. Leung, K.T., Jozsa, L., et al., 2001. Pattern formation – spiral cracks without twisting. Nature 410 (6825), 166–1166. Neda, Z., Leung, K.T., et al., 2002. Spiral cracks in drying precipitates. Phys. Rev. Lett. 88 (9), 4. Porter, M.I., Hills, D.A., 2001. Adhesive contact between a rigid punch and a halfplane via a thin, soft interlayer. J. Eng. Mech. 127 (2), 176–179. Sendova, M., Willis, K., 2003. Spiral and curved periodic crack patterns in sol–gel films. Appl. Phys. A-Mater. Sci. Proc. 76 (6), 957–959. Shevchuk, V.A., Silberschmidt, V.V., 2005. Semi-analytical analysis of thermally induced damage in thin ceramic coatings. Int. J. Solid Struct. 42 (16–17), 4738–4757. Stoney, G.G., 1909. The tension of metallic films deposited by electrolysis. Proc. Roy. Soc. Lond., Ser. A 82 (553), 172–175. Suhir, E., 1989. Interfacial stresses in bimetal thermostats. J. Appl. Mech. 56 (3), 595–600. Suhir, E., 1991. Structural Analysis in Microelectronic and Fiber-optic Systems. van Nostrand Reinhold, New York. Suhir, E., 1994. Approximate evaluation of the elastic thermal stresses in a thin film fabricated on a very thick circular substrate. J. Electron. Package 116, 171–176. Suhir, E., 2000. Predicted thermally induced stresses in, and the bow of, a circular substrate/thin-film structure. J. Appl. Phys. 88 (5), 2363–2370. Suo, Z.G., 1989. Singularities interacting with interfaces and cracks. Int. J. Solid Struct. 25 (10), 1133–1142. Thouless, M.D., 1995. Modeling the development and relaxation of stresses in films. Annu. Rev. Mater. Sci. 25, 69–96. Timm, D.H., Guzina, B.B., et al., 2003. Prediction of thermal crack spacing. Int. J. Solid Struct. 40 (1), 125–142. Xia, Z.C., Hutchinson, J.W., 2000. Crack patterns in thin films. J. Mech. Phys. Solids 48 (6–7), 1107–1131. Yin, H.M., 2010a. Fracture saturation and critical thickness in layered materials. Int. J. Solid Struct. 47, 1007–1015. Yin, H.M., 2010b. Opening-mode cracking in asphalt pavements: crack initiation and saturation. Road Mater. Pavement Des. 11, 435–457. Yin, H., Paulino, G., et al., 2008. An explicit elastic solution for a brittle film with periodic cracks. Int. J. Fract. 153 (1), 39–52. Yin, H.M., Buttlar, W.G., et al., 2007a. Periodic thermal cracking in an asphalt overlay bonded to a rigid pavement. J. Transport. Eng. 133, 39–46. Yin, H.M., Paulino, G.H., et al., 2007b. Micromechanics-based thermoelastic model for functionally graded particulate materials with particle interactions. J. Mech. Phys. Solids 55 (1), 132–160.
Contents lists available at ScienceDirect
Mechanics Research Communications journal homepage: www.elsevier.com/locate/mechrescom
Elastic thermal stresses in a circular overlay/rigid substrate system Z.F. Yuan, H.M. Yin ∗ Department of Civil Engineering and Engineering Mechanics, Columbia University, 610 Seeley W. Mudd, 500 West 120th Street, New York, NY 10027, USA
a r t i c l e
i n f o
Article history: Received 1 November 2010 Received in revised form 17 February 2011 Keywords: Thermal stresses Layered materials Thermoelasticity Axisymmetry Thin film
a b s t r a c t When an overlay/substrate system is subjected to a temperature change, thermal stress is induced due to the difference of thermal expansion coefficients between the overlay and substrate. This paper studies elastic thermal stress distribution of a circular overlay bonded to a rigid substrate due to temperature variation. Under the plane assumption, the radial displacement in the overlay is obtained in a series form. Using a weak form boundary condition, one can obtain an approximate solution in closed form for the thermoelastic field. This approximate solution gives a description of displacement along both radial and thickness directions. Elastic thermal stress distributions are demonstrated and compared with the finite element results. An asymptotic analysis of the solution is given when the thickness of a thin film is much less than its radius. The closed-form solution is useful for stress-strain analyses and design of circular layered systems such as protective coatings and microelectromechanical components. Application of this solution in modeling spiral cracking of circular thin film/substrate is underway. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction Overlay/substrate systems have been widely used in engineering structures (Timm et al., 2003; Leung and Tung, 2006; Yin et al., 2007a), electronic packaging (Suhir, 1994; Yin et al., 2008), surface coatings (Shevchuk and Silberschmidt, 2005; Yin et al., 2007b), and material testing and modeling components (Agrawal and Raj, 1989; Hutchinson and Suo, 1992; Porter and Hills, 2001). When these systems are subjected to an ambient temperature change or a service temperature load, thermal stress will be induced due to the different coefficients of thermal expansion (CTEs) between the overlay and substrate, which may cause thermal cracking within the layers or delamination along the interface. Advances in design and maintenance of these material systems require the quantitative understanding of the thermal stress distribution as part of the structural analysis process. The crack within the thin film has different patterns, some of which exhibit a very interesting phenomenon. Yin et al. (2007a,b) studied a parallel thermal cracking of asphalt overlay bonded to rigid substrate, and the cracking is nearly uniformly distributed. The solution was extended to general materials (Yin et al., 2008) and used for investigation of fracture initiation and saturation (Yin, 2010a,b). Spiral cracking was observed in the drying process of thin films by Leung et al. (2001) and Neda et al. (2002). Sendova and Willis (2003) presented four different crack patterns in sili-
∗ Corresponding author. Tel.: +1 212 851 1648; fax: +1 212 854 6267. E-mail address: [email protected] (H.M. Yin). 0093-6413/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2011.04.005
cate sol–gel films such as straight-line, curved closed-loop, crescent (sawtooth) periodic, and spiral cracks. Crack patterns in thin films may change with the material stiffness and fracture toughness of the materials, interface condition, and loading conditions. Stress and strain analyses will be valuable to understand the fracture mechanisms of an overly/substrate system and predict its cracking behavior under a certain working condition. Some theoretical research has been established for cracking propagation. Suo (1989) developed the solutions to singularities such as point force, point moment, edge dislocation and transformation strain spot in bonded elastic blocks of dissimilar materials. Thouless (1995) derived an approximate analytical solution for the minimum spacing that can exist between a series of parallel cracks propagating in a thin film and investigated mechanism of development and relaxation of stresses in films. Although fracture mechanics based modeling can provide local mechanical fields with a high fidelity, it is not straightforward to formulate thermal stress field in the whole system. The first theoretical formula for evaluation of thermal stresses in an overlay/substrate system was proposed by Stoney (1909). It has been widely used to predict stresses of thin films with the system curvatures and later extended to multilayer structures by many investigators, for example, see Freund and Suresh (2004) and Feng et al. (2008). Obviously, this school of semi-empirical approaches cannot be applied if the substrate is so stiff or so thick that no bowing occurs (Suhir, 1994). Suhir proposed the interfacial compliance concept (Suhir, 1989, 1991) and could address these special cases (Suhir, 1994). However, most existing models assumed a uniform or linear stress distribution in the thickness direction and reduced the two dimensional (2D) problem into a one dimensional (1D) problem, so that the
284
Z.F. Yuan, H.M. Yin / Mechanics Research Communications 38 (2011) 283–287
z
R
Substituting Eq. (6) into Eqs. (2)–(4), and then into Eq. (1), we obtain the governing equation as
E1 ,ν 1 , α1 , h
E0 ,ν 0 , α 0 , H
r
r 2 ur,rr + rur,r − ur 1 − 1 + ur,zz = 0 2 r2
(7)
Using the method of separation of variables, we can write the displacement in the following form as
T Fig. 1. A circular thin film boned to a rigid substrate and subjected to a temperature change T.
ur = R(r)Z(z)
stress variation in the thickness direction could not be directly studied (Suhir, 1994, 2000; Xia and Hutchinson, 2000; Timm et al., 2003). However, because the thermal stress transferred across the interface in mainly through the shear stress, which vanishes at the surface, a simple assumption of uniform or linear stress may not reflect the actual stress distribution. Yin et al. (2007a, 2008) have developed a 2D analytical model to explicitly predict the stress distribution for rectangular overlay/substrate systems. The purpose of this paper is to extend that solution to circular overlay/substrate systems. A compliant overlay bonded to a rigid substrate and subjected to uniform temperature change is considered as seen in Fig. 1, which is very common in polymer coating systems.
r 2 Rrr
(8)
Then, Eq. (7) can be rewritten as 1 − 1 Zzz + rRr − R =− = c2 2 Z r2R
where c is a constant to be determined from boundary conditions. Then we obtain Zzz +
2c 2 Z=0 1 − 1
(c 2 r 2 + 1)R − rRr − r 2 Rrr = 0
rr,r + zr,z
(1)
where is the coordinate in the circumferential direction. In the z direction, because no constraint is applied, it is assumed that zz = 0. Then the isotropic material constitutive law provides rr = =
E1 1 − 12 E1
[εrr + 1 ε − (1 + 1 )˛T ]
[1 εrr + ε − (1 + 1 )˛T ] 2
1 − 1
(2)
(3)
2 c 1 − 1
(12)
where A and B are to be determined by the boundary conditions. Because the substrate is so stiff compared with the overlay that it can be approximated as a perfectly rigid material. Therefore, the deformation of the overlay is completely constrained by the substrate. Notice that although the substrate has a thermal strain, because the thermal stresses are induced by the mismatch of CTEs between the overlay and the substrate, the effect of the thermal strain in the substrate has been included in Eqs. (2) and (3) and will not have effect on the stress analysis. Therefore, the mechanical displacement is considered to be zero. In addition, the shear stress on the surface of the overlay is zero. We can write ur (r, 0) = 0, and ur,z (r, h) = 0
(13)
Eq. (13) implies B = 0 and di =
1 2
+i
h
, with i = 0, 1, 2, ...
(14)
so that a series form of solution for Z(z) is obtained as Zi (z) = Ai sin(di z)
(15)
For each di with i = 0,1,2,.., we can obtain
ci =
1 − 1 di 2
(16)
The general solution for Eq. (11) can be written as
and zr =
(11)
The general solution for Eq. (10) is written as
2. Basic formulation
rr − =0 + r
(10)
and
Z = A sin(dz) + B cos(dz), with d =
Consider a circular overlay with thickness h, Young’s modulus E1 , Poisson’s ratio v1 , the coefficient of thermal expansion (CTE) ˛1 and radius R fully bonded to a rigid substrate with thickness H, Young’s modulus E0 E1 , Poisson’s ratio v0 , and CTE ˛0 . The overlay/substrate system is subjected to an ambient temperature change T, as illustrated in Fig. 1. Because the overlay and the substrate have different CTEs, whose difference is defined as ˛ ¯ = ˛1 − ˛0 , a thermal stress is induced in the overlay. Assume the material is thermoelastically isotropic. The thermal stresses can be solved by an axi-symmetric formulation. The coordinate r is set up along the interface between two layers, and the coordinate z is along the axisymmetric axis in the thickness direction. The equilibrium equation in r direction is written
(9)
E1 zr 2(1 + 1 )
(4)
Ri (r) = Ci I1 (ci r) + Di K1 (ci r)
(17)
where I1 (ci r) and K1 (ci r) are the modified Bessel functions. Using the axisymmetric condition
Because the substrate provides a strong constraint for the overlay, the thermal stress is mainly induced in the plane perpendicular to the z direction. Therefore, a plane assumption is introduced that all the points in a plane normal to the z direction keep in the same plane after the deformation, i.e.:
We obtain Di = 0. Therefore, the solution for Eq. (7) can be written in a series form as
uz (r, z) = uz (z)
ur (r, z) =
(5)
Using Eq. (5), we can obtain the strain-displacement relation for this axisymmetric problem, εrr = ur,r ; ε =
ur ; zr = ur,z r
(6)
ur (0, z) = 0
(18)
∞
i=0
Ai sin(di z)I1 (ci r)
(19)
where Ai can be determined by the stress free boundary condition at r = R, i.e. rr (R, z) = 0
(20)
Z.F. Yuan, H.M. Yin / Mechanics Research Communications 38 (2011) 283–287
Substituting Eq. (6) into Eq. (2) and using Eq. (20), we obtain ur,r (R, z) +
1 ur (R, z) ¯ = (1 + 1 )˛T R
= ci I0 (ci r) −
= (1 + 1 )˛T ¯
(22)
Inserting Eq. (19) into Eq. (21) and using I1 (ci r)/r, we can write ∞
Ai sin(di z) ci I0 (ci R) −
i=0
(1 − 1 )I1 (ci R) R
Extending the integral domain to (0,2h) and using the orthogonality of trigonometric functions, we can write
Ai
(1 − 1 )I1 (ci R) ci I0 (ci R) − R
1 = h
2h
(1 + 1 )˛T ¯ sin(di z)dz (23) z=0
Simplifying Eq. (23), we can analytically solve the coefficient Ai as
Ai =
¯ 2(1 + 1 )˛T
(1/2 + i) ci I0 (ci R) − (1 − 1 )I1 (ci R)/R
(24)
Overall, the elastic solution for displacement as Eq. (19) is obtained with parameter Ai ,di and ci given in Eqs. (24), (14) and (16), respectively. Then we can solve the stress field from Eqs. (2)–(4) with (6). A series form solution may reach a higher precision to the boundary condition at the circumferential surface of the thin film, but not very straightforward to use. In addition, due to the approximation of Eq. (5), the series form solution cannot guarantee a higher accuracy of the solution, which will be revisited in Section 3. In the sense of approximation, if only the first item in Eq. (14) is adopted, we can write ur (r, z) = A sin(dz)I1 (cr); d = ; and c = 2h
1 − 1 2 2h
(25)
Using this approximation, we cannot make the boundary condition in Eq. (20) exactly satisfied at each point. However, we can set the resultant force to be zero, so that we obtain: A=
(1 + 1 )˛RT ¯ 2[cRI0 (cR) − (1 − 1 )I1 (cR)]
(26)
Then we obtain a closed form solution with Eqs. (25) and (26). In summary, the elastic fields are given as ur (r, z) = rr = =
¯ (1 + 1 )˛RT sin(dz)I1 (cr) 2[cRI0 (cR) − (1 − 1 )I1 (cR)]
¯ E1 ˛T 1 − 1
R cRI (cr) − (1 − )I (cr) 0 1 1
2r cRI0 (cR) − (1 − 1 )I1 (cR)
sin(dz) − 1
R crI (cr) + (1 − )I (cr) E1 ˛T ¯ 1 0 1 1
1 − 1
Table 1 Material properties defined in FEM simulation (Suhir, 1994).
(21) I1 (ci r)
2r cRI0 (cR) − (1 − 1 )I1 (cR)
(27)
sin(dz) − 1
285
Young’s modulus kgf/cm2 Poisson’s ratio CTE (1/◦ C)
Thin film
Substrate
28 × 103 0.4 34.8 × 10−6
1.23 × 106 0.24 0
free, and the outer boundary and bottom of the substrate are fixed. We applied a temperature change of T = − 280 ◦ C. Following we will compare the theoretical predictions and FEM results of displacement and stresses. The geometric parameters of the thin film/substrate system are as follows: radius of R = 5, and thickness h = 1, as shown in Fig. 1. Figs. 2–4 illustrate the distributions of displacement ur and normal stresses rr and along the surface of the thin film. In Fig. 2, ur decreases as r increases, which means the thin film will shrink inwards. Theoretical predictions of displacement ur through Eq. (27) excellently agree with the FEM results. However, the series form solution of Eq. (19) considerably deviates from FEM results. In Fig. 3, theoretical predictions of rr and agree with the FEM results very well at the central range. However, close to the edge, FEM results of rr will approach zero, which means it is a stress-free boundary; whereas the theoretical solution monotonically decrease due to the weak boundary condition used. On the other hand, the series form solution in Eq. (22) can satisfy the stressfree boundary condition, but the accuracy in the central range is lost. The similar tendency of stress distribution of can also be found except that it is not zero at the free edge. In Figs. 2–3, the closed form solution provides very good agreement with FEM results except the points close to the edge. Although the series form solution can satisfy the stress-free boundary condition at the edge, the overall accuracy is even worse than the closed form solution. The reason to cause this problem is that the plane assumption in Eq. (5) imposes an artificial constraint on the displacement field in the vertical direction. Enforcement of the stress boundary condition at the edge amplifies the negative effect of this assumption. Further investigation of this problem in the context of mathematical physics is underway. Fortunately, the simple closed form solution provides a very good prediction of the overall elastic field. The shearing stress along the interface is of interest in design and analysis of thin film/substrate system. Fig. 4 illustrates the shearing stress rz along the interface. The closed-form solution provides a good agreement with the FEM results except the range close to the
(28)
(29)
and zr =
¯ 2 R I1 (cr) E1 ˛T cos(dz) 2 4h cRI0 (cR) − (1 − 1 )I1 (cR)
(30)
3. Comparison with finite element simulation Because the formulation developed in the last section is still an approximate solution, verification of its accuracy is important to its application. In this section, an FEM model will be developed for verification of the analytical solution. The model contains the thin film and substrate. We use ABAQUS to run the model. The material properties of thin film and substrate are given in Table 1. Notice that we will focus on the difference of CTE between thin film and substrate, which is a resource of mismatch stress. Therefore, the CTE of substrate is purposely chosen at zero. In this model, the boundaries of the thin film except the interface of substrate are
Fig. 2. Comparison of ur at top surface of thin film with radius R = 5.
286
Z.F. Yuan, H.M. Yin / Mechanics Research Communications 38 (2011) 283–287
Fig. 3. Comparison of rr (left) and (right) at top surface of thin film with radius R = 5.
To compare the closed form solution with the FEM results in z direction, Fig. 5 illustrate the distributions of ur and rz along z direction of thin film at r = R/2. The closed form solution provides an excellent agreement with the FEM results. Notice that at the free edge, we expect a considerable difference between them and thus do not show it in the z direction. In Figs. 2–5, we find a very good fit of the closed form solution and FEM results, but a considerable difference of the series form solution from the FEM results in the central range of the overlay. Although there is a slight difference at the outer part, Eq. (28) gives a precise prediction of turning point of compressive zone and tensile zone. In detail, from Fig. 3, we find at the inner part, rr is positive which denotes a tensile stress, and at the outer part about r > 4, rr is negative which denotes a compressive stress. Both FEM and Eq. (28) predict such kind of turning point, which is very important to spiral cracking development. Overall, the FEM results agree with the closed form solution very well. 4. Discussion Fig. 4. Comparison of rz at interface of thin film with radius R = 5.
edge, which is supposed to exhibit a singular effect of the shearing stress distribution. Due to the weak form boundary condition used, the proposed solution cannot capture the singularity. However, in an average sense, it provides a very good approximation.
The proposed model provides a closed form approximate solution instead of an exact solution. The assumptions and approximations adopted may impose some limitations on the application of the formulation. In summary, the following two main assumptions are used in the derivation: first, because the stress in the z direction is free, all the points in the same plane perpendicular
Fig. 5. Comparison of ur (left) and rz (right) along z direction at r = R/2 with radius R = 5.
Z.F. Yuan, H.M. Yin / Mechanics Research Communications 38 (2011) 283–287
to z axis are assumed to keep in the same plane during deformation. Therefore, the governing equation of ur is obtained. This assumption imposes a limitation of this model that the thickness and stiffness of the substrate should be large enough to make the assumption of rigid substrate applicable and to keep the plane flatten. In addition, due to this assumption, the equilibrium in the z direction cannot be rigorously satisfied at each point. Secondly, to obtain the closed-form solution, the boundary condition along the free edge surface is approximated as a weak form boundary condition that the resultant force is zero. However, it makes the stress field around the fracture surface inaccurate. The closed form solution cannot capture the singularity of stress distribution. If the local stress around the fracture surface is emphasized, the series form solution should be used. However, it is not guaranteed that the series form solution will provide higher accuracy than the closed-form solution does because of the first assumption. Generally, the normal stress in the z direction ( zz ), which is also called peeling stress, will not be zero in the neighborhood of the free edge. However, the first assumption implies zz = 0, which further affects the accuracy of this model for the stress distribution close to the free edge. Because the closed-form solution has assured the balance of the resultant force across the thickness in the radial direction and satisfied the governing equations, and is easy to apply. In addition, the substrate is assumed to be rigid in this paper. For many applications, elastic substrates may be more realistic. In these cases, an elastic overlay fully bonded to an elastic substrate needs to be considered. 5. Conclusions A series form 2D solution is derived for a circular overlay bonded to a rigid substrate and subjected to a temperature change. Using a weak form boundary condition, a closed form solution is obtained. Compared with Suhir’s 1D solution, this solution can predict the stress variation in the thickness direction. Take an average of the stress distribution, the proposed solution can be reasonably reduced to a form which is very similar to Suhir’s solution. The explicit thermal stress solution is useful to structure design and thermoelastic analysis for circular overlay/rigid substrate systems. The assumptions and limitations of this model are discussed and analyzed. Future research efforts are emphasized on releasing some of these assumptions and combining this model with numerical methods to conduct fracture simulation.
287
Acknowledgment This work is sponsored by the National Science Foundation CMMI 0954717 and the Department of Homeland Security CU091155, whose support is gratefully acknowledged. References Agrawal, D.C., Raj, R., 1989. Measurement of the ultimate shear-strength of a metal ceramic interface. Acta Metall. 37 (4), 1265–1270. Feng, X., Huang, Y., et al., 2008. Multi-layer thin films/substrate system subjected to non-uniform misfit strains. Int. J. Solid Struct. 45 (13), 3688–3698. Freund, L.B., Suresh, S., 2004. Thin Film Materials; Stress Defect Formation and Surface Evalution. Cambridge University Press, Cambridge, UK. Hutchinson, J.W., Suo, Z., 1992. Mixed-mode cracking in layered materials. Advances Appl. Mech., vol. 29. Academic Press Inc., San Diego, pp. 63–191. Leung, C.K.Y., Tung, W.K., 2006. Three-parameter model for debonding of FRP plate from concrete substrate. J. Eng. Mech. 132 (5), 509–518. Leung, K.T., Jozsa, L., et al., 2001. Pattern formation – spiral cracks without twisting. Nature 410 (6825), 166–1166. Neda, Z., Leung, K.T., et al., 2002. Spiral cracks in drying precipitates. Phys. Rev. Lett. 88 (9), 4. Porter, M.I., Hills, D.A., 2001. Adhesive contact between a rigid punch and a halfplane via a thin, soft interlayer. J. Eng. Mech. 127 (2), 176–179. Sendova, M., Willis, K., 2003. Spiral and curved periodic crack patterns in sol–gel films. Appl. Phys. A-Mater. Sci. Proc. 76 (6), 957–959. Shevchuk, V.A., Silberschmidt, V.V., 2005. Semi-analytical analysis of thermally induced damage in thin ceramic coatings. Int. J. Solid Struct. 42 (16–17), 4738–4757. Stoney, G.G., 1909. The tension of metallic films deposited by electrolysis. Proc. Roy. Soc. Lond., Ser. A 82 (553), 172–175. Suhir, E., 1989. Interfacial stresses in bimetal thermostats. J. Appl. Mech. 56 (3), 595–600. Suhir, E., 1991. Structural Analysis in Microelectronic and Fiber-optic Systems. van Nostrand Reinhold, New York. Suhir, E., 1994. Approximate evaluation of the elastic thermal stresses in a thin film fabricated on a very thick circular substrate. J. Electron. Package 116, 171–176. Suhir, E., 2000. Predicted thermally induced stresses in, and the bow of, a circular substrate/thin-film structure. J. Appl. Phys. 88 (5), 2363–2370. Suo, Z.G., 1989. Singularities interacting with interfaces and cracks. Int. J. Solid Struct. 25 (10), 1133–1142. Thouless, M.D., 1995. Modeling the development and relaxation of stresses in films. Annu. Rev. Mater. Sci. 25, 69–96. Timm, D.H., Guzina, B.B., et al., 2003. Prediction of thermal crack spacing. Int. J. Solid Struct. 40 (1), 125–142. Xia, Z.C., Hutchinson, J.W., 2000. Crack patterns in thin films. J. Mech. Phys. Solids 48 (6–7), 1107–1131. Yin, H.M., 2010a. Fracture saturation and critical thickness in layered materials. Int. J. Solid Struct. 47, 1007–1015. Yin, H.M., 2010b. Opening-mode cracking in asphalt pavements: crack initiation and saturation. Road Mater. Pavement Des. 11, 435–457. Yin, H., Paulino, G., et al., 2008. An explicit elastic solution for a brittle film with periodic cracks. Int. J. Fract. 153 (1), 39–52. Yin, H.M., Buttlar, W.G., et al., 2007a. Periodic thermal cracking in an asphalt overlay bonded to a rigid pavement. J. Transport. Eng. 133, 39–46. Yin, H.M., Paulino, G.H., et al., 2007b. Micromechanics-based thermoelastic model for functionally graded particulate materials with particle interactions. J. Mech. Phys. Solids 55 (1), 132–160.