《复合材料 Composites》课程教学资源（学习资料）第五章 陶瓷基复合材料_Prediction of the fracture toughness of a ceramic multilayer composite

Availableonlineatwww.sciencedirect.co ScienceDirect Acta materialia ELSEVIER Acta Materialia 55(2007)409-421 www.actamat-journals.com Prediction of the fracture toughness of a ceramic multilayer composite Modeling and experiments C R Chen a. Pascual b F D. Fischer .o. Kolednik d, *R. Danzer b Materials Center Leoben Forschung GmbH, Franz-Josef-Strasse 13, A-8700 Leoben, austria truktur und Funktionskeramik. Montanumiversitat Leoben. Peter. Tunner. Strasse 5.A-8700 Leoben. austria Institute of Mechanics, Montanumirersitat Leoben, Franz-Josef-Strasse 18, A-8700 d erich Schmid Institute of Materials Science, Austrian Academy of Sciences, Jahnstrasse 12, 4-8700 Leoben, Austria Received 24 March 2006: received in revised form 28 June 2006: accepted 2 July 2006 Available online 7 November 2006 Abstract equires experiments to measure the intrinsic fracture toughness of the phases and to determine the required material data. The numerical modeling includes a conventional finite element stress analysis and the calculation of the crack driving force based on the concept of configurational(material)forces. The procedure yields the fracture toughness of the composite as a function of the crack length. a bend- bar consisting of layers made of alumina and an alumina-zirconia composite is investigated The bar has a crack perpendicular to the interfaces. The spatial variations of both the thermal residual stresses and the elastic modulus induce shielding and anti-shielding effects on the crack, which are quantified. The numerically predicted fracture toughness is compared with the experimental values. o 2006 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved Keywords: Ceramics: Multilayer; Fracture toughness; Crack tip shielding: Thermal residual stresses 1. Introduction oxide fuel cells [8] but also in structural applications such as dental restorations or hip replacements [5] To improve the performance against brittle failure in A prerequisite for the application of multilayer ceramic materials, strategies have been developed in recent is the understanding of their resistance against crack initi- years to design tough and strong ceramics consisting of ation and propagation. Different authors have predicted multilayers. These strategies include the tailoring of weak the toughening effect of the residual stress state by means interfaces for crack deflection or the design of multilayers of the weight function method [1, 3, 9, 10]. They used the with compressive residual stresses in the outer layers [1]. classical weight function concept to calculate the stress The beneficial consequences of compressive stresses at the intensity factor, considering an inhomogeneous distribu rface are well known: increases in strength, fracture tion of the residual stresses in a homogeneous body (mostly oughness and reliability [2, 3]. Additionally, the wear resis- with the elastic modulus of the first layer). According to tance is enhanced [4]. Therefore, such ceramics receive seri- Fett et al. [11, 12] an approximate weight function method ous consideration for structural application [5-7]. Ceramic can also be applied to heterogeneous, graded or laminated multilayers have already found functional applications in materials with a variable Youngs modulus. substrates for low-load-bearing integrated circuits or solid The immanent inhomogeneity of the material, however causes implications which are not taken into account by the weight function method: spatially varying material proper- Tel: +43 114: fax: +43 3842804116. ties induce an additional crack driving force term. The nik(aunileoben ac at(O. Kolednik ) propagation of a crack in a direction orthogonal to the 1359-6454/$30.00 O 2006 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:l0.1016 actant200607.046

Prediction of the fracture toughness of a ceramic multilayer composite – Modeling and experiments C.R. Chen a , J. Pascual b , F.D. Fischer c , O. Kolednik d,*, R. Danzer b a Materials Center Leoben Forschung GmbH, Franz-Josef-Strasse 13, A-8700 Leoben, Austria b Institut fu¨r Struktur- und Funktionskeramik, Montanuniversita¨t Leoben, Peter-Tunner-Strasse 5, A-8700 Leoben, Austria c Institute of Mechanics, Montanuniversita¨t Leoben, Franz-Josef-Strasse 18, A-8700 Leoben, Austria d Erich Schmid Institute of Materials Science, Austrian Academy of Sciences, Jahnstrasse 12, A-8700 Leoben, Austria Received 24 March 2006; received in revised form 28 June 2006; accepted 2 July 2006 Available online 7 November 2006 Abstract A procedure to predict the fracture toughness of a ceramic multilayer composite made of different phases is presented. The procedure requires experiments to measure the intrinsic fracture toughness of the phases and to determine the required material data. The numerical modeling includes a conventional finite element stress analysis and the calculation of the crack driving force based on the concept of configurational (material) forces. The procedure yields the fracture toughness of the composite as a function of the crack length. A bend￾ing bar consisting of layers made of alumina and an alumina–zirconia composite is investigated. The bar has a crack perpendicular to the interfaces. The spatial variations of both the thermal residual stresses and the elastic modulus induce shielding and anti-shielding effects on the crack, which are quantified. The numerically predicted fracture toughness is compared with the experimental values. 2006 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Ceramics; Multilayer; Fracture toughness; Crack tip shielding; Thermal residual stresses 1. Introduction To improve the performance against brittle failure in ceramic materials, strategies have been developed in recent years to design tough and strong ceramics consisting of multilayers. These strategies include the tailoring of weak interfaces for crack deflection or the design of multilayers with compressive residual stresses in the outer layers [1]. The beneficial consequences of compressive stresses at the surface are well known: increases in strength, fracture toughness and reliability [2,3]. Additionally, the wear resis￾tance is enhanced [4]. Therefore, such ceramics receive seri￾ous consideration for structural application [5–7]. Ceramic multilayers have already found functional applications in substrates for low-load-bearing integrated circuits or solid oxide fuel cells [8], but also in structural applications such as dental restorations or hip replacements [5]. A prerequisite for the application of multilayer ceramics is the understanding of their resistance against crack initi￾ation and propagation. Different authors have predicted the toughening effect of the residual stress state by means of the weight function method [1,3,9,10]. They used the classical weight function concept to calculate the stress intensity factor, considering an inhomogeneous distribu￾tion of the residual stresses in a homogeneous body (mostly with the elastic modulus of the first layer). According to Fett et al. [11,12], an approximate weight function method can also be applied to heterogeneous, graded or laminated materials with a variable Young’s modulus. The immanent inhomogeneity of the material, however, causes implications which are not taken into account by the weight function method: spatially varying material proper￾ties induce an additional crack driving force term. The propagation of a crack in a direction orthogonal to the 1359-6454/$30.00 2006 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actamat.2006.07.046 * Corresponding author. Tel.: +43 3842 804 114; fax: +43 3842 804 116. E-mail address: kolednik@unileoben.ac.at (O. Kolednik). www.actamat-journals.com Acta Materialia 55 (2007) 409–421

C R Chen et al. Acta Materialia 55(2007)409-421 laminate planes can be promoted (anti-shielding) or diamond paste. The procedure yields a notch with a tip retarded(shielding) by the different elastic properties(elas- radius of about 15 um(measured at the lateral surfaces tic mismatch) of the laminae [1]. The spatially varying before fracture). This should assure a reliable fracture residual stress state can have a similar effect [13]. toughness measurement [15](see below ). Four-point-bend- The goals of this paper are, therefore, to present an ing tests were conducted in a commercial Zwick Z010 alternative procedure -the method of configurational (or machine under room conditions (34% relative humidity material) forces- to predict the fracture toughness of a and a temperature of 20C) ceramic multilayer composite which takes these spatial The fracture toughness Kic was determined from the inhomogeneities into account. These numerical predictions expression [14] will then be compared with experimental fracture tough less data. Furthermore, we will outline a generally applica- K F S,-s, 3 (1) ble concept of how to exploit the method of configurational Bh h 2(1-7) forces which would also be applicable for the design of lay- with ered composite materials with improved fracture resistance 2. Materials and processing Tapes of alumina(A)and an alumina-zirconia compos- Y=1.9887-13267 3.49-0687+1.35y2)(1-y) e(az were produced by tape casting at ISTEC Faenza, Italy. The AZ composite consists of 60 vol. Al_O3 and 40 vol %3Y-TZP, which is ZrO, with 3 mol %Y,O3. The In the above relations, F is the fracture load and S, sheets were alternatively stacked, forming a multilayer that (=20 mm) and S2(=10 mm) are the support span lengths was sintered at 1550C. As result, a seven-layer laminate The initial crack length, a, was measured after the tests, ta- was obtained with the structure A/AZ/A/AZ/A/AZ/A, in ken as an average of three measurements on the fracture which the thickness of the individual A- and AZ-layers surface are 190 and 220 um, respectively. More details about pro- Three specimens of the multilayer composite with vari- cessing can be found elsewhere, especially data about the ous initial crack lengths were tested(see Fig. Ib). In the composition of the slurry necessary to produce the tapes [6]. first specimen the crack tip was located in the middle of A complete characterization of the material layers was the first layer, in the second specimen it was shortly behind performed. The results are presented in Table 1. The elastic and in the third shortly beyond the interface with the sec- onstants, Youngs modulus E and Poisson's ratio v, were ond layer (interface 1). The initial geometries and the measured at room temperature by an impulse excitation results of the experiments are listed in Table 2. As well as chnique. The coefficient of thermal expansion( CTE)a the Kic values, the loads at fracture Ffr are also given. was measured between room temperature and the reference Additionally, approximate JIc values are listed which are temperature 1160C by means of a dilatometer (see calculated from the relation J=(1-22)K/E, using vol- below). The average CtE values are listed in Table 1. Also ume-averaged values of Poissons ratio, v=0. 25, and presented are the total thickness of the four A-layers (A, Youngs modulus, E=375 GPa nd the three Az-layers tAz in the multilayer are given The intrinsic fracture toughness Ko of the A and aZ which are needed in Appendix I material was determined by testing multilayered homoge- neous specimens consisting of only A-an 3. Fracture mechanics experiments respectively. Indentation toughness values following [16] were determined, and these are presented in Table 1. In The fracture toughness of the composite was measured addition, VAMAS-ESIS experiments were also performed on single edge V-notched beams, following the VAMAs- For the A-layers the VAMAS-ESIS procedure gave exactly ESIS procedure [14]. Bar-shaped specimens with length the same value as the indentation procedure: K0=3.8+0.3 L=28 mm, width h= 1.42 mm and thickness B 2.5 mm MPa vm. For the AZ-layers, however, the VAMAS-ESIS were cut from the original plates using a diamond saw(see values show a big scatter and higher mean value: KAZ Fig. la). The notches were machined in a home-made auto- 5.4+ 1.0 MPa vm compared with K0=4.3+0.1 matic device which uses a razor blade sprinkled with MPa vm from the indentation tests. There are two possi- Table l propertIes f(mm) E(GPa) x(10-6K-) Ko(MPa vm) AlO3(A 392±5 0.24±0.04 8.64±0.03 3.8±0.2 35±2 AlO -ZrO(AZ 0.66 305±4 0.26±0.03 9.24±0.02 4.3±0.1 7士3

laminate planes can be promoted (anti-shielding) or retarded (shielding) by the different elastic properties (elas￾tic mismatch) of the laminae [1]. The spatially varying residual stress state can have a similar effect [13]. The goals of this paper are, therefore, to present an alternative procedure – the method of configurational (or material) forces – to predict the fracture toughness of a ceramic multilayer composite which takes these spatial inhomogeneities into account. These numerical predictions will then be compared with experimental fracture tough￾ness data. Furthermore, we will outline a generally applica￾ble concept of how to exploit the method of configurational forces which would also be applicable for the design of lay￾ered composite materials with improved fracture resistance. 2. Materials and processing Tapes of alumina (A) and an alumina–zirconia compos￾ite (AZ) were produced by tape casting at ISTEC Faenza, Italy. The AZ composite consists of 60 vol.% Al2O3 and 40 vol.% 3Y-TZP, which is ZrO2 with 3 mol.% Y2O3. The sheets were alternatively stacked, forming a multilayer that was sintered at 1550 C. As result, a seven-layer laminate was obtained with the structure A/AZ/A/AZ/A/AZ/A, in which the thickness of the individual A- and AZ-layers are 190 and 220 lm, respectively. More details about pro￾cessing can be found elsewhere, especially data about the composition of the slurry necessary to produce the tapes [6]. A complete characterization of the material layers was performed. The results are presented in Table 1. The elastic constants, Young’s modulus E and Poisson’s ratio m, were measured at room temperature by an impulse excitation technique. The coefficient of thermal expansion (CTE) a was measured between room temperature and the reference temperature 1160 C by means of a dilatometer (see below). The average CTE values are listed in Table 1. Also presented are the total thickness of the four A-layers tA, and the three AZ-layers tAZ in the multilayer are given which are needed in Appendix 1. 3. Fracture mechanics experiments The fracture toughness of the composite was measured on single edge V-notched beams, following the VAMAS– ESIS procedure [14]. Bar-shaped specimens with length L = 28 mm, width h = 1.42 mm and thickness B 2.5 mm were cut from the original plates using a diamond saw (see Fig. 1a). The notches were machined in a home-made auto￾matic device which uses a razor blade sprinkled with diamond paste. The procedure yields a notch with a tip radius of about 15 lm (measured at the lateral surfaces before fracture). This should assure a reliable fracture toughness measurement [15] (see below). Four-point-bend￾ing tests were conducted in a commercial Zwick Z010 machine under room conditions (34% relative humidity and a temperature of 20 C). The fracture toughness KIC was determined from the expression [14] KIC ¼ F B ffiffiffi h p S1 S2 h 3 ffiffi c p 2ð1 cÞ 1:5 Y ð1Þ with c ¼ a h and Y ¼ 1:9887 1:326c ð3:49 0:68c þ 1:35c2Þcð1 cÞ ð1 þ cÞ 2 ð2Þ In the above relations, F is the fracture load and S1 (=20 mm) and S2 (=10 mm) are the support span lengths. The initial crack length, a, was measured after the tests, ta￾ken as an average of three measurements on the fracture surface. Three specimens of the multilayer composite with vari￾ous initial crack lengths were tested (see Fig. 1b). In the first specimen the crack tip was located in the middle of the first layer, in the second specimen it was shortly behind and in the third shortly beyond the interface with the sec￾ond layer (interface 1). The initial geometries and the results of the experiments are listed in Table 2. As well as the KIC values, the loads at fracture Ffr are also given. Additionally, approximate JIC values are listed which are calculated from the relation J ¼ ð1 m2ÞK2 =E, using vol￾ume-averaged values of Poisson’s ratio, m ¼ 0:25, and Young’s modulus, E ¼ 375 GPa. The intrinsic fracture toughness K0 of the A and AZ material was determined by testing multilayered homoge￾neous specimens consisting of only A- and AZ-layers, respectively. Indentation toughness values following [16] were determined, and these are presented in Table 1. In addition, VAMAS–ESIS experiments were also performed. For the A-layers the VAMAS–ESIS procedure gave exactly the same value as the indentation procedure: KA 0 ¼ 3:8  0:3 MPa pm. For the AZ-layers, however, the VAMAS–ESIS values show a big scatter and higher mean value: KAZ 0 ¼ 5:4  1:0 MPa pm compared with KAZ 0 ¼ 4:3  0:1 MPa pm from the indentation tests. There are two possi￾Table 1 Material properties Material t (mm) E (GPa) m (–) a (106 K1 ) K0 (MPa pm) J0 (J/m2 ) Al2O3 (A) 0.76 392 ± 5 0.24 ± 0.04 8.64 ± 0.03 3.8 ± 0.2 35 ± 2 Al2O3–ZrO2 (AZ) 0.66 305 ± 4 0.26 ± 0.03 9.24 ± 0.02 4.3 ± 0.1 57 ± 3 410 C.R. Chen et al. / Acta Materialia 55 (2007) 409–421

C R Chen et al. Acta Materialia 55(2007)409-421 919922 h=1.42 mm aX Cross-section F/2 symmetry plane Longitudinal section z=0 Fig. 1. Geometry of the four-point-bending test arrangement of the laminate specimen consisting of layers of A- and AZ-material Results of the fracture toughness tests on the multilayer composite Tip location Fr(N IC(M 117 n AZ 8.4 ble explanations for this difference. The first is that some lie in planes parallel to the y-z plane. The arrangement cutting-induced phase transformation occurs when produc- and thickness of the laminae can be taken from the long ing the notch. In the machine used, the force used while tudinal section at z=0(Fig. lb). In this longitudinal sec cutting the notch cannot be controlled. The second possible tion, an integration path I is marked surrounding a reason for the discrepancy is that the notch radius of about rectangle Q2far with the area hLy. The length Ly will be var 5 um is certainly small enough for the alumina, with its ied to obtain different paths T. The six interfaces, num- grain size of approximately 5 um, but it is quite large for bered from 1 to 6, intersect the integration path I at the AZ-microstructure, which has a grain size of approxi- y=0 and y= Ly. The normal unit vectors to the interfaces mately 0.7 um. This could cause an overestimation of the as well as to the integration path I are designated n. The fracture toughness. Therefore, to be on the safe side, we crack with variable length a is located in the plane y=0; decided to use the lower value of the indentation measure- the crack front is assumed to be parallel to the line x=0 ment in the paper. Fig. Ic presents a cross-section at y=0. The average Ko values and the corresponding Jo values are also listed in Table 1. It is seen that the fracture tough- 4.2. The residual stress state ness values of the composite are significantly larger than the corresponding intrinsic values of the A and Az-materials The specimens are fabricated at high temperatures. Due It should be noted here that the Jo values characterize the to the difference in the Cte, a cooling from the sintering fracture initiation toughness of the materials; in the ceramic temperature to the room temperature 20C introduces a ommunity the term"fracture energy is often used residual stress state(see e.g. [I]. At high temperatures, relaxation processes prevent the development of residual 4. Numerical modeling stresses. The main reason is that the az-material exhibits extensive plasticity between 1200C and the sintering tem- 4. 1. Description of the model perature [17]. Therefore, an upper reference temperature is the effective temperature The global setting is depicted in Fig. la. The laminate is calculated. This reference temperature was determined in beam is supported at a distance S, of 20 mm and loaded [18]as Tref=1160C, leading to an effective temperature by a pair of loads at a distance S2 of 10 mm. The laminae difference AT=-l140C

ble explanations for this difference. The first is that some cutting-induced phase transformation occurs when produc￾ing the notch. In the machine used, the force used while cutting the notch cannot be controlled. The second possible reason for the discrepancy is that the notch radius of about 15 lm is certainly small enough for the alumina, with its grain size of approximately 5 lm, but it is quite large for the AZ-microstructure, which has a grain size of approxi￾mately 0.7 lm. This could cause an overestimation of the fracture toughness. Therefore, to be on the safe side, we decided to use the lower value of the indentation measure￾ment in the paper. The average K0 values and the corresponding J0 values are also listed in Table 1. It is seen that the fracture tough￾ness values of the composite are significantly larger than the corresponding intrinsic values of the A- and AZ-materials. It should be noted here that the J0 values characterize the fracture initiation toughness of the materials; in the ceramic community the term ‘‘fracture energy’’ is often used. 4. Numerical modeling 4.1. Description of the model The global setting is depicted in Fig. 1a. The laminate beam is supported at a distance S1 of 20 mm and loaded by a pair of loads at a distance S2 of 10 mm. The laminae lie in planes parallel to the y z plane. The arrangement and thickness of the laminae can be taken from the longi￾tudinal section at z =0 (Fig. 1b). In this longitudinal sec￾tion, an integration path C is marked surrounding a rectangle Xfar with the area hLy. The length Ly will be var￾ied to obtain different paths C. The six interfaces, num￾bered from 1 to 6, intersect the integration path C at y = 0 and y = Ly. The normal unit vectors to the interfaces as well as to the integration path C are designated n. The crack with variable length a is located in the plane y = 0; the crack front is assumed to be parallel to the line x = 0. Fig. 1c presents a cross-section at y = 0. 4.2. The residual stress state The specimens are fabricated at high temperatures. Due to the difference in the CTE, a cooling from the sintering temperature to the room temperature 20 C introduces a residual stress state (see e.g. [1]). At high temperatures, relaxation processes prevent the development of residual stresses. The main reason is that the AZ-material exhibits extensive plasticity between 1200 C and the sintering tem￾perature [17]. Therefore, an upper reference temperature is introduced from which the effective temperature difference is calculated. This reference temperature was determined in [18] as Tref = 1160 C, leading to an effective temperature difference DT = 1140 C. Fig. 1. Geometry of the four-point-bending test arrangement of the laminate specimen consisting of layers of A- and AZ-material. Table 2 Results of the fracture toughness tests on the multilayer composite Specimen B (mm) h (mm) a (mm) Tip location Ffr (N) KIC (MPa pm) J appr IC ðJ=m2 Þ 1 2.72 1.42 0.10 In A 117 6.1 98 2 2.64 1.42 0.21 In AZ 110 8.4 188 3 2.64 1.42 0.18 In A 103 7.3 142 C.R. Chen et al. / Acta Materialia 55 (2007) 409–421 411

C R Chen et al. Acta Materialia 55(2007)409-421 The calculation procedure to find the residual stress and aaz= aAZ -A=0.6x 10-K- as the substitute state in laminate specimens has been heavily investigated CTE in the AZ-material. The second 2D computation used (see e.g. [19, 20]. We use a three-dimensional (3D) model the corresponding plane stress model of the specimen, consisting of a slice of length 8 in the y- direction and covering the area0≤x≤M2,0≤z≤B/2.4.3. Beam bending Symmetry conditions are applied along 2=0 and x= h/ 2, and the y-displacements are fixed on the bottom of the As outlined above, four-point-bend tests are performed layer. Only in the end regions far from the crack plane will on rather slender beams. The classical beam-bending the the residual stress state depend on the y-coordinate. Since ory could be used to evaluate the stress state the specimen is very long in y-direction, the residual stres- uncracked composite beam; for details, see e.g ses can be assumed to be independent of the y-coordinate. Since we need to determine the stress state in the specimen Therefore, an unknown but spatially constant displace- with a crack of length a, finite element calculations are per ment u] is assumed along the upper boundary y=8. The formed. The beam is replaced by a two-dimensional plane finite element program system ABAQUS(htp!∥ strain model covering the area0≤x≤ h and y≥0.Note www.hks.com)isengagedforthecomputationusingthattheplanestrainmodelcanbetreatedasaplanestress eight-node 3D solid elements model by replacing Youngs modulus E by E= E/(1-v) It is well known that the behavior of a crack is deter- Only x-displacements ux are allowed at y=0. The speci mined not only by the stresses perpendicular to the crack men is fixed at the point P, in the x-direction; the load F plane y=0. The in-plane and out-of-plane constraints also is applied at point P2(see Fig. 1). The mesh consists of play a role, i.e. the stresses in x-and z-directions. The mate- eight-node plane strain elements rial can freely move in x-direction, so that no correspond To model a realistic stress state of the fracture mechan- g residual stresses will appear In the z-direction residual ics specimens, the following procedure is applied.First, the stresses o, Ar will appear which infiuence the strain energy uncracked and unloaded specimen is subjected to a thermal loading by a temperature difference AT=-1140C to cal The materials of the A- and AZ-layers are modeled as culate the thermal residual stresses. Then in the unloaded linear-elastic. The corresponding material properties are specimen a crack of length a is introduced by a node release taken from Table 1. The longitudinal residual stresses o, Ar technique. Subsequently, the specimen is loaded by pre in the symmetry plane z=0 and in the side-surface plane scribing the load Fat the load application point. The final z= B/2 are plotted between 0<x<h/2 in Fig. 2. Fig. 3 stress and strain distribution within the specimen is used shows the variation of o, Ar from the midsection to the for the evaluation of the crack driving force, which is side-surface for sections in the middle of the various described beloy laminae In addition to the 3D computation, two simple 2D com- 4. 4. Calculation of the crack driving force putations were also performed. The first one used a plane strain model covering the area0≤x≤h,0≤y≤ s with The study of stress intensity factors, as well as of their assumes no displacement in the 2-direction, u,=0. To mechanics research in composite materials (see e. g Simha, Kolednik et al. [13, 26-28]. References to the open formulations can be taken from these extensive papers Specifically, Section 3 of [28] provides the corresponding equations which are reshaped below in the specific form a composites with constant material properties within each lamina 。-atz=1.25of3 D model The concept of configurational forces considers a mate- rial inhomogeneity as an additional defect in the material (besides the crack) which induces an additional contribu- tion to the crack driving force. This contribution has been plane stress b called the material inhomogeneity term Cinh. The thermo- dynamic force at the crack tip, denominated as the local near-tip crack driving force Jtip, is the sum of the nominally x[ mm] applied far-field crack driving force Far and the material inh [28]: Fig. 2. Thermal residual stresses o,. AT along the x-direction at 2=0 and 25 mm for a plane stress, plane strain and three-dimensional (3D)model Jtip=Far Ci

The calculation procedure to find the residual stress state in laminate specimens has been heavily investigated (see e.g. [19,20]). We use a three-dimensional (3D) model of the specimen, consisting of a slice of length d in the y￾direction and covering the area 0 6 x 6 h/2, 0 6 z 6 B/2. Symmetry conditions are applied along z = 0 and x = h/ 2, and the y-displacements are fixed on the bottom of the layer. Only in the end regions far from the crack plane will the residual stress state depend on the y-coordinate. Since the specimen is very long in y-direction, the residual stres￾ses can be assumed to be independent of the y-coordinate. Therefore, an unknown but spatially constant displace￾ment uy is assumed along the upper boundary y = d. The finite element program system ABAQUS (http:// www.hks.com) is engaged for the computation, using eight-node 3D solid elements. It is well known that the behavior of a crack is deter￾mined not only by the stresses perpendicular to the crack plane y = 0. The in-plane and out-of-plane constraints also play a role, i.e. the stresses in x- and z-directions. The mate￾rial can freely move in x-direction, so that no correspond￾ing residual stresses will appear. In the z-direction residual stresses rz,DT will appear which influence the strain energy density (see Appendix 1). The materials of the A- and AZ-layers are modeled as linear-elastic. The corresponding material properties are taken from Table 1. The longitudinal residual stresses ry,DT in the symmetry plane z = 0 and in the side-surface plane z = B/2 are plotted between 0 6 x 6 h/2 in Fig. 2. Fig. 3 shows the variation of ry,DT from the midsection to the side-surface for sections in the middle of the various laminae. In addition to the 3D computation, two simple 2D com￾putations were also performed. The first one used a plane strain model covering the area 0 6 x 6 h, 0 6 y 6 d with unit thickness in the z-direction. The plane strain model assumes no displacement in the z-direction, uz ” 0. To avoid any stresses due to the global shrinkage of the spec￾imen, we set a A ¼ 0 as the substitute CTE in the A-material and a AZ ¼ aAZ aA ¼ 0:6  106 K1 as the substitute CTE in the AZ-material. The second 2D computation used the corresponding plane stress model. 4.3. Beam bending As outlined above, four-point-bend tests are performed on rather slender beams. The classical beam-bending the￾ory could be used to evaluate the stress state in the uncracked composite beam; for details, see e.g. [21,22]. Since we need to determine the stress state in the specimen with a crack of length a, finite element calculations are per￾formed. The beam is replaced by a two-dimensional plane strain model covering the area 0 6 x 6 h and y P 0. Note that the plane strain model can be treated as a plane stress model by replacing Young’s modulus E by E* = E/(1 m 2 ). Only x-displacements ux are allowed at y = 0. The speci￾men is fixed at the point P1 in the x-direction; the load F is applied at point P2 (see Fig. 1). The mesh consists of eight-node plane strain elements. To model a realistic stress state of the fracture mechan￾ics specimens, the following procedure is applied. First, the uncracked and unloaded specimen is subjected to a thermal loading by a temperature difference DT = 1140 C to cal￾culate the thermal residual stresses. Then in the unloaded specimen a crack of length a is introduced by a node release technique. Subsequently, the specimen is loaded by pre￾scribing the load F at the load application point. The final stress and strain distribution within the specimen is used for the evaluation of the crack driving force, which is described below. 4.4. Calculation of the crack driving force The study of stress intensity factors, as well as of their relevance for crack growth, has been a topic of fracture mechanics research in composite materials (see e.g. [1,23–25]). In the current investigation the concept of con- figurational forces is used. Here we refer to the works by Simha, Kolednik et al. [13,26–28]. References to the open literature with respect to this concept and other related formulations can be taken from these extensive papers. Specifically, Section 3 of [28] provides the corresponding equations which are reshaped below in the specific form for composites with constant material properties within each lamina. The concept of configurational forces considers a mate￾rial inhomogeneity as an additional defect in the material (besides the crack) which induces an additional contribu￾tion to the crack driving force. This contribution has been called the material inhomogeneity term Cinh. The thermo￾dynamic force at the crack tip, denominated as the local, near-tip crack driving force Jtip, is the sum of the nominally applied far-field crack driving force Jfar and the material inhomogeneity term Cinh [28]: Jtip ¼ Jfar þ Cinh: ð3Þ Fig. 2. Thermal residual stresses ry,DT along the x-direction at z = 0 and 1.25 mm for a plane stress, plane strain and three-dimensional (3D) model. 412 C.R. Chen et al. / Acta Materialia 55 (2007) 409–421

C R Chen et al. Acta Materialia 55(2007)409-421 plane strain ri plane stress plane strain z[mm Fig 3. Thermal residual stresses , AT along the z-direction for various values of y for a plane stress, plane strain and three-dimensional (3D) model. Jfar is the classical J-integral of fracture mechanics. For a done by specifying the set of nodes on the interface as crack growing in the x-direction, Jfar is virtual crack tip nodes. Even a contour directly adjacent to the interface yields very accurate results. For the eval- (4) uation of the J-integral, the virtual crack growth direc tion must be specified; this is the(,0)-direction, as for The components t;(t;=fr,fy)of the traction vector t along he evaluation of Je the contour I follow from the stress tensor g with the rel-. Generally, both Far and Cinh depend on the crack length evant components ox, Oy, txy as I=g n. Note that the a. They also depend on Ly, but produce Jtip values stress components are the sum of their contributions due according to Eq .(3)which are independent of L,; for to bending and the residual stress state, e.g. 0,=0x, b+ details see Appendix 2 At. The components u;(u;=ux, ly) of the displacement vector u are differentiated with respect to the crack growth After the finite element stress analysis, the material i direction, i. e. the x-direction. The quantity o is the specific mogeneity term Cinh is calculated from Eq. (5)by a elastic strain energy and nx is the x-component of the unit processing procedure. The integration along the interface normal vector n to the integration path I is performed using the trapezoid formula. For this, the The material inhomogeneity term can be evaluated by node values of the stress and strain components and the strain energy density are taken, which are extrapolated val ues from the Gauss integration points. The far-field J-inte- Cnh=>Cinh i, Cinh =-2/(l,-(g)le)dy (5) gral Jar is calculated using the virtual crack extension method of ABAQUS. Then the near-tip crack driving force The jump [b] and the average (b)of a quantity b at an Jtip is calculated from Eq (3). The numerical results will be interface are defined presented in the following section [b]=b-b,(b)=(b+b) It should be noted that Sun and Wu[29]have calculated the effective crack driving force by replacing the region far where br and b denote the limiting values of the quantity b by subregions, each including only one layer, and applying on the right and left side of the interface, respectively. The the J-integral procedure for each individual layer. The index i refers to the individual interface; the integer I de- strength of the configurational forces concept lies in its gen notes the total number of interfaces in the specimen, in eral applicability. The material inhomogeneity can be our case I =6 either a sharp interface with a discrete jump of the material The following comments may be useful properties or a region where the material properties change continuously. The Cinh-evaluation procedure can be The multiplier 2 in Eqs. (4)and(5)points to the fact that applied to any arbitrary spatial distribution of these mate- only the upper half of a symmetric configuration with rial inhomogeneities in both elastic and elastic-plastic respect to I is considered. materials. In general, the evaluation of Cinh can be per The material inhomogeneity terms Cinh. i can be also formed very accurately. This enables us to evaluate Jtip as found via the J-integral calculation routine provided the sum of Far and Cinh more accurately than it would be by the finite element code by evaluating the J-integral possible from the calculation of Jtip using the conventional around the ith interface Jint. i[26]. In ABAQUS, this is J-evaluation procedures, especially for cases when the

Jfar is the classical J-integral of fracture mechanics. For a crack growing in the x-direction, Jfar is Jfar ¼ 2 Z C /nx ti oui ox ds: ð4Þ The components ti (ti = tx,ty) of the traction vector t along the contour C follow from the stress tensor r with the rel￾evant components rx, ry, sxy as t = r Æ n. Note that the stress components are the sum of their contributions due to bending and the residual stress state, e.g. ry = ry,b + ry,DT. The components ui (ui = ux,uy) of the displacement vector u are differentiated with respect to the crack growth direction, i.e. the x-direction. The quantity / is the specific elastic strain energy and nx is the x-component of the unit normal vector n to the integration path C. The material inhomogeneity term can be evaluated by [28] Cinh ¼ XI i¼1 Cinh;i; Cinh;i ¼ 2 Z Ly 0 ð½½/i hrii½½eiÞdy ð5Þ The jump [[b]] and the average Æbæ of a quantity b at an interface are defined as ½½b ¼ br bl; hbi¼ðbl þ brÞ=2; ð6Þ where br and bl denote the limiting values of the quantity b on the right and left side of the interface, respectively. The index i refers to the individual interface; the integer I de￾notes the total number of interfaces in the specimen, in our case I = 6. The following comments may be useful:  The multiplier 2 in Eqs. (4) and (5) points to the fact that only the upper half of a symmetric configuration with respect to C is considered.  The material inhomogeneity terms Cinh,i can be also found via the J-integral calculation routine provided by the finite element code by evaluating the J-integral around the ith interface Jint,i [26]. In ABAQUS, this is done by specifying the set of nodes on the interface as virtual crack tip nodes. Even a contour directly adjacent to the interface yields very accurate results. For the eval￾uation of the J-integral, the virtual crack growth direc￾tion must be specified; this is the (1, 0)-direction, as for the evaluation of Jfar.  Generally, both Jfar and Cinh depend on the crack length a. They also depend on Ly, but produce Jtip values according to Eq. (3) which are independent of Ly; for details see Appendix 2. After the finite element stress analysis, the material inho￾mogeneity term Cinh is calculated from Eq. (5) by a post￾processing procedure. The integration along the interface is performed using the trapezoid formula. For this, the node values of the stress and strain components and the strain energy density are taken, which are extrapolated val￾ues from the Gauss integration points. The far-field J-inte￾gral Jfar is calculated using the virtual crack extension method of ABAQUS. Then the near-tip crack driving force Jtip is calculated from Eq. (3). The numerical results will be presented in the following section. It should be noted that Sun and Wu [29] have calculated the effective crack driving force by replacing the region Xfar by subregions, each including only one layer, and applying the J-integral procedure for each individual layer. The strength of the configurational forces concept lies in its gen￾eral applicability. The material inhomogeneity can be either a sharp interface with a discrete jump of the material properties or a region where the material properties change continuously. The Cinh-evaluation procedure can be applied to any arbitrary spatial distribution of these mate￾rial inhomogeneities in both elastic and elastic-plastic materials. In general, the evaluation of Cinh can be per￾formed very accurately. This enables us to evaluate Jtip as the sum of Jfar and Cinh more accurately than it would be possible from the calculation of Jtip using the conventional J-evaluation procedures, especially for cases when the Fig. 3. Thermal residual stresses ry,DT along the z-direction for various values of y for a plane stress, plane strain and three-dimensional (3D) model. C.R. Chen et al. / Acta Materialia 55 (2007) 409–421 413