Energy Release Rates for Interface Cracks in Multilayered Structures

This paper examines the evolution of the interfacial deflection energy release rates in multilayered structures under four-point bending. The J-integral and the extended finite element method (XFEM) are adopted to investigate the evolution of the interfacial deflection energy release rates of composite structures. Numerical results not only verify the accuracy of analytical solutions for the steady-state interfacial deflection energy release rate, but also provide the evolutionary history of the interfacial deflection energy release rate under different crack lengths. In addition, non-dimensional parametric analyses are performed to discuss the effects of normalized ratios of the crack length, the elastic modulus, and the thickness on the interfacial deflection energy release rate. The results demonstrate that the elastic modulus ratio and thickness ratio have a distinct influence on the interfacial deflection energy release rate for multilayered beams. Furthermore, an unstable interfacial crack tends to occur for elastic multilayer beams with higher elastic modulus on the upper sub-beam under bending moments. The unstable interfacial fracture shows a decreasing interfacial deflection energy release rate with an increasing interfacial crack length.


Introduction
Multilayered structures have been widely adopted in microelectronic, semiconductors, optical, protective coating, and composite structures. Different types of bonding, such as inter diffusion bonds, chemical interaction bonds, thermal-compression bonds, and adhesive bonds, are developed in manufacturing these advanced devices [Abdelhadi, Ladani and Razmi (2011)]. The functionality and reliability of these multilayered devices largely depend on the ability of the bonds at the interfaces to sustain mechanical and thermal stresses during the fabrication process and the service life. As a result, it is inevitable to develop different measured techniques when assessing the interfacial adhesion of multilayered structures [Dannenberg (1961); William (1969); Vossen (1978); Hinkley (1983); Valli (1986) ;Steinmann, Tardy and Hintermann (1987); Westerlind, Larsson and Rigdahl (1987); Stone, LaFontaine, Alexopoulos et al. (1988); Charalambides, Lund, Evans et al. (1989); Bull (1997) ;Zhou, Zhang and Cai (2002); Modi and Sitaraman (2004) ;Zhao, Sun, Zhu et al. (2010); Gadelrab, Chiesa, Hecker et al. (2012)]. One popular technique to measure the interface adhesion strength of a bi-layered system with dissimilar materials is the four-point bending adhesion test, which was proposed by Charalambides et al. [Charalambides, Lund, Evans et al. (1989)]. In the adhesion test, a notch is cut from the top of the specimen to provide an initiation point of the interfacial crack, as shown in Fig. 1. The interfacial crack begins growth on the bonding interface between the upper layer and the lower layer from the notch tip and then propagates along the concerned interface under load. In addition to the reliable and convenient measurements over a wide range of multilayered specimens [Klingbell and Beuth (1997); An, Qin and Li (2011)], the most attractive character of the adhesion test is that a steady-state interfacial deflection energy release rate can be measured under the four-point bending. The terminology of "steadystate" implies that the interfacial deflection energy release rate remains the same value under the applied loadings when the crack length exceeds a lower limit and is within two inner supports. This special property is useful in the experimental procedure because engineers or scientists do not need to accurately monitor the crack length during a test. Charalambides et al. [Charalambides, Lund, Evans et al. (1989)] developed an analytical formulation for the steady-state interfacial deflection energy release rate in a bi-layered system. Klingbell et al. [Klingbell and Beuth (1997)] extended this adhesion test under the four-point bending to trilayered systems and developed the corresponding analytical solutions. Hsueh et al. [Hsueh, Tuan and Wei (2009)] further derived the analytical solutions of the steady-state interfacial deflection energy release rate for multilayered systems in the adhesion test. However, to our knowledge, few research studies investigate the critical length (lower limit) of interface cracks corresponding to the steady-state deflection energy release rate of multilayered structures in the adhesion test. In addition, the evolution of interfacial deflection energy release rates in multilayered structures is seldom found in existing literature. This study follows the linear elastic fracture mechanics (LEFM) and the extended finite element method (XFEM), using the J-integral method, to investigate the evolution of the interfacial deflection energy release rate on the interfaces of composite beams. Numerical results not only provide the evolutionary history of the interfacial deflection energy release rate with different crack lengths, but also help to determine the critical length (lower limit) of the interface cracks to reach a steady-state interfacial deflection energy release rate for different combinations of material and sectional properties. To simplify the analyses, we do not include the residual stress-induced delamination in multilayered systems. In addition, small-scale yielding is assumed, i.e., the plastic deformation near the crack tip is not considered. In other words, the present study focuses solely on interfacial cracking resulting from external loads under four-point bending in elastic multilayered systems.

Methods
We adopted different methods, including the closed-form solution, the J-integral method based on LEFM, and XFEM, to estimate the interfacial deflection cracking energy release rates for multilayered systems in the four-point bending test. The following sections describe the details of these methods.

Analytical solution
A multilayered beam which consists of n layers is taken into consideration. Each layer has individual thickness ti and elastic modulus Ei. In this study, the subscript i represents the layer number and the first layer (i=1) is at the bottom of the multilayered beam. The coordinate system is defined such that the origin (z=0) is located at the bottom of the multilayered beam and the top of the beam is at z=hn, as shown in Fig. 2. The transformation between the height hi and the thickness ti of the i -th layer is given by: Without losing generality, the interfacial crack can propagate on the interface between the m -th layer and the (m+1) -th layer. The corresponding steady-state interfacial deflection energy release rate is obtained by Hsueh et al. [Hsueh, Tuan and Wei (2009) where P is the external applied load, L is the distance between the roller supports and the applied loads (as shown in Fig. 2), and b is the width of the composite beam. Notably, Eq.
(2) was derived for plane stress conditions. The elastic modulus Ei in Eqs. (2), (4) and (5) should be replaced by Ei/(1-vi 2 ) for plane strain conditions, where νi is the Poisson's ratio of the i -th material.

J-integral method
LEFM considers three distinct fracture types, modes I, II, and III, which encompass all possible ways a crack tip may deform. In LEFM, the stress field around the crack tip is described through the stress intensity factor, K. The stress intensity factor, K, which depends on the geometries of the specimen, the size and location of the crack, and the applied external loads is defined from the elastic stresses close to a sharp crack under remote loading. For an interfacial crack between two dissimilar materials, with a crack length of 2a, the singular stress on the interface at distance r in front of the crack tip can be written as [Sun and Jih (1987)]: where KI and KII represent the stress intensity factors for fracture modes I and II, respectively. The parameter ε is related to the material properties of the two dissimilar materials in adjacent layers and can be defined by 3 4 for the case of plane strain ( 1, 2) (3 ) / (1 ) for the case of plane stress ( 1, 2) where µ and ν are the shear modulus and Poisson's ratio of materials, respectively.
In addition to the local stress intensity factors, the energy-based fracture toughness is an alternative global parameter to quantify the characteristics of cracks. The energy release rate, G, is the negative value of the variation in the potential energy per unit length of the crack extension under the same load, and is given by: where U is the potential energy and a is the length of the crack. Theoretically speaking, the energy release rate of the considered interfacial crack is obtained from the change of the compliance under different crack lengths with infinitesimal crack length increment. However, this computation becomes difficult when complicated stress fields of interfacial cracks are considered in the four-point bending adhesion test. Therefore, the J-integral method is adopted to estimate the global fracture toughness, i.e., the energy release rate.
The J-integral is path-independent when the line integral contour is around the crack tip. A common definition of the J-integral is given by Rice [Rice (1968)]: where Γdenotes the possible contour route around the crack tip, w is the strain energy density. In addition, t and u represent the traction and the corresponding displacement vector, respectively, and ds is an infinitesimal piece of the route length of contour Γ. For linear-elastic materials, the magnitude of J-integral equals the value of the energy release rate associated with crack advance. In this study, concentric circular contours surrounding the crack tip, as shown in Fig. 3, are adopted as the integral paths. The energy release rate is then obtained by averaging these J-values from different circular contours.

XFEM
Compared to other finite element-based crack simulation methods, XFEM enables an accurate prediction of the crack initiation and growth along arbitrary paths without a prespecified direction [Belytschko and Black (1999)]. It is well-accepted that XFEM has several advantages, such as (1) prediction of the crack propagation along an arbitrary direction, (2) independent of mesh for crack initiation and propagation, (3) easier definition of the initial crack, (4) simpler mesh refinement to improve the convergence rates, and (5) wide applicability in static and implicit dynamic analyses [Mohammadi (2008)].
In addition to the polynomial function in the standard finite element methods, XFEM adopts extra interpolation functions, which include the Heaviside step function and the crack tip asymptotic function. The former represents the discontinuous displacement fields across the crack faces and the later considers the singularity of stress fields around the crack tip. As a result, XFEM integrates the discontinuous geometries and the singular stress (and strain) fields due to cracks for the displacement interpolation. The displacement interpolation in XFEM can be written as [Fries and Baydoun (2012)]: where u h (x) represents the displacement at a point x in the domain, and NI(x) and uI are the conventional polynomial shape function and the continuous nodal displacement used in the standard finite element method, respectively. In addition, H(x) is the Heaviside step function, aI is the nodal enrichment degree of freedom to describe the jump discontinuity across the crack surfaces, Fα(x) is the asymptotic enrichment function for the crack tip, and bI α is the additional nodal degree of freedom associated with the asymptotic enrichment function. The asymptotic enrichment function for the crack tip is [Lecampion (2009) where r and θ denote the position of the nodes of the element containing the crack tip in terms of the local polar coordinate system. The origin of the local polar coordinate system is set at the crack tip, while θ=0 is set to be tangent to the crack at the tip.
The most promising development in the XFEM is the description of cracks. To facilitate the treatments of cracks, XFEM permits a crack to be located in the element interior and the corresponding mesh is not required to conform to the crack geometry. These advantages are achieved by implementing the level set method (LSM) which enables the representation of surfaces in interface tracking problems [Hansbo and Hansbo (2004); Song, Areias and Belytschko (2006)]. For a real number function, a level set of is the sub-set in the design domain at which the function is equal to a pre-specified value. Generally, two orthogonal signed functions, Φ and Ψ, are required for a complete description of the cracks. The nodal values of these two functions, Φ and Ψ, represent the distance from the crack face to the node and the perpendicular distance from the crack front to the the node, respectively. The interaction of these two level sets, Φ=0 and Ψ=0, gives the crack tip, as shown in Fig. 4. At each iteration, crack growth is modeled by updating these two functions, Φ and Ψ, at nodes which are belonged to elements cut by the crack.

Numerical results
To verify the analytical solution for the interfacial deflection energy release rate Gs versus the crack length and to determine the critical crack length for which the interfacial deflection energy release rate of an interfacial crack reaches a steady-state condition in multilayered systems, finite element models, based on LEFM and XFEM, are constructed using the software package Abaqus. This study considers a composite specimen with length of 64 mm, height of 4 mm, and width of 1 mm. Two concentrated loads of 2,000 N each are applied at 2 mm away from the two ends of the beam. The distance between the outer loading line and inner support is 10 mm, and the distance between the two inner supports is 40 mm. To simplify the analyses, only alternate layered systems where the layer management is ABAB in sequence are investigated. In other words, the first layer on the bottom of the composite beam and other odd-number layers are made of material A, while the even-number layers are made of material B. The elastic modulus of material A is 64 GPa, while the Poisson's ratio of both materials is 0.21. In addition, non-dimensional parameters related to the elastic modulus ratio, the layer thicknesses ratio, the crack lengths ratio, and the normalized energy release rates are introduced and defined as follows:

Bilayered system
For bilayered systems, with n=2 and m=1 in Eqs.
(2)-(4), an interfacial crack propagates along the interface between the top and bottom layers. Figs. 5 and 6 demonstrate the evolution of the energy release rates of the interfacial cracks, obtained by LEFM and XFEM, with various values of elastic modulus ratio, RE, for thickness ratios Rt=0.25 and Rt=4.0 respectively. The dashed lines in Figs. 5 and 6 represent the analytical solution of the steady-state deflection energy release rates for interfacial cracks. One can observe that the energy release rates from LEFM and XFEM show good agreements with each other and converge to the analytical solutions for the steady-state interfacial deflection energy release rates. Some discrepancy in the energy release rates from LEFM and XFEM result mainly from mesh qualities. In addition, the interfacial deflection energy release rate achieves a steady-state rate at different crack length ratios according to the normalized elastic modulus RE and the normalized thickness Rt. In the case of Rt=0.25, the steady-state energy release rates are achieved when the crack length ratio, Ra, is greater than 0.25. However, in the case of Rt=4.0, the energy release rate achieves a steady-state value when the crack length ratio, Ra, is greater than 0.5. It is also found that the deflection energy release rates of the interfacial cracks in the case of Rt=4.0 are much larger than those for Rt=0.25. The results imply that the increase of thickness of the top-layer sub-beam could prevent the propagation of the interfacial crack, corresponding to a large value of the energy release rate. In such cases, the elastic modulus of the top layer sub-beam only slightly affects the steady-state deflection energy release rate of the interfacial crack. However, in the case of Rt=0.25, where the thickness of the top layer sub-beam is smaller than that of the bottom layer sub-beam, the elastic modulus of the top layer sub-beam significantly affects the steady-state interfacial deflection energy release rate.
Notably, the interfacial fractures shows decreasing deflection energy release rates with increasing the interfacial crack length for RE=5.0 and Rt=0.25. The evolution of the deflection energy release rate with the decreasing trend implies that the interfacial crack will propagate automatically without increasing external loads until reaching a steady state. This phenomenon, for which the evolution of the deflection energy release rate has a decreasing trend, is contrary to that of a stable crack and becomes more obvious for the stiffer and thinner top layer sub-beams in bi-layered systems.

Tri-layered systems
For tri-layered systems, with n=3 and m=1 in Eqs.
(2)-(4), the interfacial crack propagates along the interface between the middle and bottom layers. Figs. 7(a) and 7(b) demonstrate the evolution of the deflection energy release rates of the interfacial cracks for different thickness ratios Rt and crack length ratio Ra under RE=0.2 and RE=5.0, respectively. From Fig. 7(a), one can observe that the deflection energy release rates increase with increasing crack length ratios and reach a steady-state value when Ra=0.15 for RE=0.2. In addition, the increase of the thickness ratio Rt also increases the deflection energy release rate, which indicates that a thicker upper sub-beam provides greater adhesion strength to the interfacial crack. However, the deflection energy release rates show a decreasing trend with an increasing crack length ratio in the case of RE=5.0. The results imply that such interfacial cracks are unstable when the elastic modulus of the upper sub-beam is larger than that of the lower sub-beam. Since the layer management is ABA in the considered tri-layered systems, the top layer provides some resistance to the interfacial cracks between the middle and bottom layers. Thus, the interfacial crack becomes unstable when the thickness ratio Rt increases, as shown in Fig. 7(b). For tri-layered systems, with n=3 and m=2 in Eqs.
(2)-(4), the interfacial crack propagates along the interface between the middle and top layers. In such a case, the middle layer plays the role of the substrate. Figs. 8(a) and 8(b) demonstrate the evolution of the deflection energy release rates of the interfacial cracks for different thickness ratios Rt and crack length ratios Ra under RE=0.2 and RE=5.0, respectively. In the case of RE=0.2, the top layer (material A) has a larger elastic modulus than that of the middle layer. Therefore, one can observe that the deflection energy release rates decrease as the crack length ratio increases and reaches a steady state after Ra=0.15 in Fig. 8(a). In addition, the increase in the thickness ratio Rt slightly decreases the deflection energy release rate and the gradient of the interfacial deflection energy release rate with respect to the crack length ratio Ra appears to be less steep at higher Rt values, which indicates that the instability has reduced slightly (but still unstable) for the thinner upper sub-beam. Nevertheless, the evolution of the deflection energy release rates show a rising trend with an increasing crack length ratio in the case of RE=5.0 in Fig. 8(b), which result from the upper sub-beam having a smaller elastic modulus than that of the lower sub-beam for the interfacial crack. These numerical results imply that such interfacial cracks are stable. One can observe that the energy release rates reach their steady-state values after Ra=0.10.

Conclusions
This study investigated the interfacial defection energy release rates in multilayer structures under the four-point bending adhesion test. Finite element models, were constructed using the software package Abaqus to determine the interfacial deflection energy release rate. The numerical results have good agreement with the analytical solution for the steady-state deflection energy release rate. Furthermore, the critical crack length for which the interfacial deflection energy release rate achieves a steady state is obtained from our simulations. The numerical results, from LEFM and XFEM, show good agreement with the analytical solutions. In addition, from the evolutionary histories of the interfacial deflection energy release rates under different crack lengths, the interfacial crack lengths should be larger than half of the system height in order to reach their steady-state condition. In addition, non-dimensional parametric analyses for interfacial deflection energy release rate were also performed to discuss the effects of the ratios of the crack length, the elastic modulus, and the thickness on the interfacial deflection energy release rate in this study.
Numerical simulations indicate that increasing the elastic modulus of the upper layer in multilayered structures may enhance the interfacial deflection energy release rate. However, unstable interfacial cracks occur when the upper layer has a larger elastic modulus than that of the lower layer in multilayer systems. Although only interfacial cracks of bi-and tri-layered structures are simulated in this study, it is worth emphasizing that these results can be extended to multilayered structures.