Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 45, 4, pp. 739-752, Warsaw 2007 INFLUENCE OF ANISOTROPY ON THE ENERGY RELEASE RATE GI FOR HIGHLY ORTHOTROPIC MATERIALS Paweł Grzegorz Kossakowski Kielce University of Technology, Faculty of Civil and Environmental Engineering, Kielce, Poland e-mail: kossak@tu.kielce.pl The paper presents results of a numerical analysis concerning the ener- gy release rate, GI, for highly orthotropicmaterials such as composites, laminates or wood. The values of GI were calculated using the Adina v. 8.1 Finite Element Method (FEM) program. Different material mo- delswere considered to establish the influence of anisotropyon GI. Two- dimensional (2D) and three-dimensional (3D) isotropic and anisotropic models were employed to study the performance of a Double Cantile- ver Beam (DCB) with various crack length-to-thickness ratios. It was reported that the smaller the ratio, the bigger the difference between the energy release rates GI calculated for the isotropic and anisotropic (transversal isotropic and orthotropic) material models. Thus, it is im- portant that a fracture or fracture toughness analysis should be based on the transversal isotropic and orthotropic models and it should take into account anisotropy. Key words: fracture toughness, pinewood, highly orthotropic materials, mode I loading, energy release rate GI Notations a – crack length B – specimen width Ei – elastic modulus in the i-direction (i=1,2,3) GI – mode I energy release rate GIc – critical value of mode I energy release rate Gi,j – shear modulus in the ij-orientation (i,j=1,2,3; i 6= j) H – specimen half-thickness h – cantilever thickness i,j – indices 740 P.G. Kossakowski P – load Pc – critical load U – potential energy νi,j – Poisson’s ratio in the ij-orientation (i,j =1,2,3; i 6= j) EL,ET ,ER – elastic modulus of wood in the L-, T- and R-direction, respectively GLT ,GLR,GTR – shear modulus of wood in the LT-,LR- and TR-orientation, respectively νLT ,νTL,νLR,νRL,νTR,νRT – Poisson’s ratio ofwood in the LT-,TL-,LR-, RL-, TR- and RT-orientation, respectively 1. Introduction A variety of inhomogeneous materials are applied in structural engineering. Themost commonare composites, laminates andwood.As structural elements made of such materials are prone to cracking, they need to be designed for particular service conditions. Their strength and stiffness decrease with crack extension, which may eventually lead to failure. One of the most frequent loading modes that occur in orthotropic struc- tural elements is the opening, denoted by mode I. In order to avoid fracture, it is essential to assess the crack driving force and fracture toughness. This requires determination of the energy release rate GI, which defines the load causing delamination. The energy release rate can be determined analytically as a function of the load P. Its critical value, GIc, is a parameter that cha- racterizes fracture toughness of amaterial formode I loading. It is established on the basis of GI for the critical load Pc determined experimentally at the crack initiation. As analytical solutions are suitable only for simple loading cases, specimens and crack configurations, their application is limited. More general solutions are obtained by means of numerical methods. Various nu- merical techniques have been used to calculate GI. They were described, for example, by Hellen (1975), Parks (1977), Nikishkov and Vaynshtok (1980), Haber and Koh (1985), and Seweryn (1998). It should be noted that, unlike analytical solutions, numerical techniques are reported to be very useful in fracture analyses. They make it possible to determine the energy release ra- te of a structural element, GI, for any geometry, any material (for instance, isotropic and anisotropic), and any crack location. This paper is concerned with the influence of anisotropy on the energy release rate, GI. It is possible to determine how anisotropy affects fracture to- Influence of anisotropy on the energy... 741 ughness of highly orthotropicmaterials suchas composites, laminates orwood. It enables us to usepropermaterialmodels in fracture toughness examinations and calculations of GI for any structural element containing a crack. In the analysis, wood was selected to act as a model of an inhomogeneous material. The analysis results are found to be universal and, thus, applica- ble to other anisotropic materials with less complex structure, for instance, composites and laminates. The principles of fracture mechanics have been determined for isotropic materials mainly. Being a complex and highly anisotropic material, wood is definedon thebasis of an orthotropicmodel.As can be seen inFig.1, there are threemaindirections inwood: longitudinal (L), tangential (T) andradial (R). Fig. 1. Orthotropic material model of wood Consequently, there are six main systems of crack propagation in wood denoted by TL, TR,RL,RT ,LR and LT (Fig.2). The first letter stands for the direction normal to the crack plane, while the other denotes the direction of crack propagation. In order to characterize the fracturemechanics of wood formode I loading, it is necessary to determine six critical values of the energy release rate, GIc. Fig. 2. Systems of crack propagation in wood 742 P.G. Kossakowski A fracture toughness analysis of an inhomogeneous material requires ta- king into account anisotropy. A majority of authors involved in the research into themode I loading treatwoodas an isotropicmaterial, for example,Atack et al. (1961), Porter (1964), DeBaise et al. (1966), Johnson (1973), Schniewind and Lyon (1973), Schniewind andCenteno (1973), Schniewind (1977), Petter- son and Bodig (1983), Yeh and Schniewind (1992), Ando et al. (1992), and Ando andOhta (1995). In order to determine fracture toughness of wood they often applied the same solutions as those for metal testing. Some authors, for instance,Wu (1967), Stanzl-Tschegg et al. (1995, 1996), Reiterer et al. (2000), and Kossakowski (2004), studied wood anisotropy basing on the orthotropic model. To estimate fracture toughness of wood, it is necessary to specify under what conditions the isotropic model is used. This problem has been analyzed by several authors, yet the focuswasonanother fractureparameter of anisotro- picmaterials, the so-called stress intensity factor.According toPatton-Mallory andCramer (1987), the isotropic stress intensity can be used to determine the orthotropic stress intensity under certain circumstances. For infinite bodies, where the crack surface is not subjected to any local loads, the isotropic and orthotropic stress intensity factors were the same (Cook and Rau, 1974; To- min, 1972;Williams andBirch, 1976). For a finite element with real boundary conditions, a change in thematerial or geometry causes that the isotropic and orthotropic stress intensity factors are different (Tomin, 1971). For a finite orthotropic plate with a center crack, subjected to tensile isotropic stress, the intensity factors lied within 10 percent of the orthotropic values (Bowie and Freese, 1972;Ghandi, 1972). However, for specimens subjected to uniform ten- sion and pure bending, the isotropic and orthotropic stress intensity factors differed less than 7 percent (Walsh, 1972). Mandel et al. (1974) found the difference between isotropic and orthotropic stress intensity to be as much as 25 percent for an edge crack in a cantilever beam specimen. It seems that the isotropic stress intensity factors canbe employed to estimate orthotropic stress intensity, yet they depend on the calculation accuracy and problem type. 2. Numerical procedures TheAdinaversion8.1FiniteElementMethod(FEM)programwasused for the numerical analysis. The calculations were made for Double Cantilever Beam (DCB) specimens subjected to mode I loading (Fig.3). Influence of anisotropy on the energy... 743 Fig. 3. DCB specimen Two- (2D) and three-dimensional (3D) elements were considered. The geo- metry and boundary conditions were the same for the 2D and 3D numerical models. As the DCB specimen is symmetrical, only half was modeled for nu- merical calculations. In order to establish the effect of the element geometry on the energy release rate for different material models, it was assumed that the cantilever thickness, h, and the element thickness, 2H, were different. Four types of 2D and 3D DCB elements were examined. The dimensions of the numerical models with the cantilever thickness ranging from 4.5mm up to 9.0mm are presented in Table 1. Table 1.Dimensions and loadings of DCBmodels Element No. B [mm] H [mm] h [mm] a [mm] a/h a/H P [N] 2D-1 18 4.5 3.6 80.0 22.22 17.78 16.46 2D-2 18 6.0 5.1 80.0 15.69 13.33 26.75 2D-3 18 7.5 6.6 80.0 12.12 10.67 38.00 2D-4 18 9.0 8.1 80.0 9.88 8.89 50.00 3D-1 18 4.5 3.6 80.0 22.22 17.78 16.46 3D-2 18 6.0 5.1 80.0 15.69 13.33 26.75 3D-3 18 7.5 6.6 80.0 12.12 10.67 38.00 3D-4 18 9.0 8.1 80.0 9.88 8.89 50.00 The fracture analysis module was used to perform a linear-elastic fracture mechanics analysis. A stationary crack with a singularity in the crack tip was modeled for all the elements.Middle nodes in the vicinity of the crack tipwere moved to the so-called Quarter-Point Location. All the two-dimensional (2D) specimens were modeled using eight-node elements in a plane stress system. Accordingly, twenty-node elements were used to build three-dimensional (3D) models. The numerical models of the DCB specimens for 2D and 3D analyses are shown in Figs. 4 and 5. 744 P.G. Kossakowski Fig. 4. Two-dimensional (2D) DCB element (dimensions in mm) Fig. 5. Three-dimensional (3D) DCB element (dimensions in mm) Thevirtual crack extensionmethodwasused tocalculate theenergy release rates directly on the basis of changes in the potential energy of the element. In theAdinaManual (Adina 1999), one reads that the virtual crack extension method evaluates the J-integral for a given body using the difference in the total potential energybetween twoconfigurationswith slightly different cracks. The energy release rates, GI, were calculated on the basis of the value of the J-parameter adopted directly from the Adina program. Studying the results of his previous research (Kossakowski, 2004), the au- thor assumed that the critical value of the energy release rate for pinewo- od, GIc, was 228 J/m 2. In the analysis, the loads P were considered to be the same for all the examined elements. The load values had to be calibrated so that the energy release rates, GI, were in the elastic range. They had to be smaller than the critical energy release rates, i.e. GI < GIc = 228J/m 2. The loads P were calibrated to obtain a constant value of GI ≈ 209J/m 2 for the transversal-isotropicmaterialmodel, which is used in thePolish standards (PN-EN 338) as the basic material model of wood. Table 1 shows values of the loads P for particular models. Influence of anisotropy on the energy... 745 3. Material models In order to analyze the influence of different material models on the values of GI for pinewood, the isotropic and anisotropic, i.e. transversal isotropic and orthotropic models, were considered. First, the elastic constants were calculated for pinewood using the ortho- tropic material model. Three main directions were taken into account: longi- tudinal (L), tangential (T), and radial (R). The values of the elastic moduli, EL, ET , ER, and shear moduli, GLT , GLR, GTR, were determined in accor- dance with the Polish standards (PN-EN 380, PN-EN 384 and PN-EN 408), while Poisson’s ratios, νLT , νLR, νRT , were established on the basis of the sta- tistical methods proposed byBodig andGoodman (1973). The values of νTL, νRL, νTR were calculated from the following constants EL, ET , ER, GLT , GLR, GTR, νLT , νLR, νRT , assuming the symmetry of the compliance matrix of wood. Table 2 presents the values of all the 12 elastic constants obtained for pinewood (Kossakowski, 2004). Table 2.Elastic constants of pinewood Elastic Value constant EL 6919MPa ET 271 MPa ER 450 MPa GLT 262 MPa GLR 354 MPa GTR 33.8 MPa νLT 0.387897 νLR 0.374817 νRT 0.462337 νTL 0.015187 νRL 0.024377 νTR 0.278327 Four material models were analyzed for all the numerical 2D and 3D ele- ments: • isotropic (ISO), described by 2 elastic constants: E1 and ν12 • transversal-isotropic (TRANS), describedby 5 elastic constants: E1,E2, G12, ν12 and ν23 746 P.G. Kossakowski • orthotropic for the TLorientation (TL),describedby9elastic constants: E1,E2,E3,G12,G13,G23, ν12, ν13 and ν23 • orthotropic for the RLorientation (RL), describedby9elastic constants: E1,E2,E3,G12,G13,G23, ν12, ν13 and ν23. For the isotropic and transversal-isotropic models, the elastic modulus E1 was assumed to be equal to EL. The other constants, i.e. elastic and shear moduli andPoisson’s ratioswere themeanvalues for the isotropic plane TRon the basis of elastic constants presented inTable 1 prepared for the orthotropic model.For example, the elasticmodulus E2 for the transversal-isotropicmodel was equal to the mean value of moduli ET and ER, i.e. E2 =(ET +ER)/2. Inwooden structural elements, cracks usuallypropagate in the longitudinal direction (L), so two crack propagation systems, TL and RL, were conside- red for the orthotropic material models. In this case, elastic constants were determined basing on the values presented in Table 1. 4. Results The analysis shows that the energy release rates, GI, are dependent on the material model, whether two- and three-dimensional DCB ones. The values of GI for all the considered material models are presented in Fig.6. Fig. 6. Changes in GI for different material models; (a) 2Dmodels; b) 3Dmodels As can be seen, there are several relationships between the energy release rates GI, material models and geometry of elements. Influence of anisotropy on the energy... 747 Firstly, the values of GI are the same for the orthotropic (TL and RL) as well as transversal-isotropic (TRANS) 2D and 3D elements. In the isotropic (ISO) 2D models, the values of GI are slightly higher than for 3D elements, but the differences are small, i.e. 1-1.5%. Another relationship is the influence of the element geometry, i.e. the ratio of the crack length a (cantilever length) to the cantilever thickness h on the values of GI. For the isotropic model, GI decreases linearly along with the element thickness. In this case, the difference between the values of GI is 9.6% for models 2D-1 and 2D-4, and 10.1% for models 3D-1 and 3D-4. The ratio a/h is 22.22 and 9.88, respectively. A similar interdependence is observed in the case of an orthotropic element in the orientation RL, but the difference is smaller, 1% only. For an orthotropic element in the orientation TL, a reverse relationship was reported; the values of GI increase with the element thick- ness, and the difference is small – 1.4%. Because the load P was calibrated to the constant values of GI for the transversal-isotropic material model, GI is independent of the element geometry. For the isotropic material model, the energy release rates, GI, are affected by the specimen geometry, thus for the anisotropic material models there is slight influence of geometry on the ener- gy release rates GI. It can be assumed that the energy release rates GI are independent of geometry whenever anisotropy is taken into consideration. Fig. 7. Changes in the ratio of GI determined for the anisotropic material models (TRANS, TL and RL) to GI determined for the isotropic material model (ISO): (a) 2Dmodels; (b) 3Dmodels It seems important to study the values of the energy release rates, GI, determined for the anisotropic models, as they are higher than those reported for the isotropic models (Figs. 6 and 7). The differences between the ratios 748 P.G. Kossakowski of GI determined for the anisotropic material models to GI determined for the isotropic material models change linearly from 7.1-23.1% for 2D elements and 8.5-24.0% for 3Dnumerical elements (Fig.7). Thus, the differences betwe- en the energy release rates, GI, determined for the anisotropic and isotropic material models increase if the a/h ratio decreases. The information is very important as it can be used in practice. In the case of elements containing a relatively short crack in comparison to its thickness, it is essential to ta- ke anisotropy into account so as to properly determine fracture toughness. Assuming that the difference between GI determined for the isotropic and anisotropic models is 5%, it is possible to estimate the corresponding a/h ra- tio, thus a/h=25.When a/h> 25, the difference between GI determined for the isotropic and anisotropic models should be less than 5%. In such a case, it is possible to use the isotropic material model instead. For a/h < 25, the differences between GI determined for the isotropic and anisotropic models are significant and higher than 5%. In order to calculate fracture toughness properly and accurately, material anisotropy should be taken into considera- tion. These relationships can be used for wood and other highly anisotropic materials, especially when high accuracy of results is required. This is impor- tant in the case of laminates, and some composites, when the scatter of elastic constants is small. As fracture toughness is dependent on thematerial model, it is necessary that, in calculations, the material anisotropy should be taken into consideration. As can be seen, the values of GI for the transversal-isotropic and orthotro- picmodels are similar (Fig. 6). It canbeassumed that in technical calculations it is possible to use the transversal-isotropic model described by 5 elastic con- stants: E1,E2,G12, ν12 and ν23. As mentioned above, it is possible to determine the critical values of the energy release rate, GIc, by employing a numerical method and basing on the critical load, Pc. Then, the material fracture toughness can be established also numerically. The relationships described above can be applied in a frac- ture toughness analysis of pinewood, other species of wood and other highly orthotropic materials such as composites and laminates. 5. Conclusions The following conclusions have been drawn basing on the results of the nume- rical analysis of the influence of anisotropy on the energy release rate, GI, for pinewood. Influence of anisotropy on the energy... 749 • The numerical analysis showed that the energy release rate GI depends on thematerial model. • Similar values of the energy release rates GI are obtained for both aniso- tropic models, i.e. the orthotropic TL and RL and transversal-isotropic TRANS ones. They are higher than those reported for the isotropicmo- dels. • It can be assumed that the values of the energy release rate, GI, are independent of the geometry of the anisotropic material models. Howe- ver, for the isotropicmodel, GI is dependent on the specimen geometry, i.e. the ratio of the crack length, a, to the cantilever thickness, h. The energy release rate, GI, for the isotropic model decreases when there is a fall in the a/h ratio. • Differencesbetween theenergy release rates GI observed for theanisotro- pic and isotropicmaterialmodels increasewhen the a/h ratio decreases. If a/h< 25, thematerial anisotropy should be taken into consideration, otherwise fracture toughness cannot be calculated precisely. • Thevalues of GI for the transversal-isotropic and orthotropicmodels are similar, thus for practical reasons, the calculations canbemadeusing the transversal-isotropicmodel.Accordingly, the number of elastic constants required for describing the material is reduced to 5 only: E1, E2, G12, ν12 and ν23. References 1. ADINA,1999,TheoryandModellingGuide.Volume I:ADINA,ADINASystem OnlineManuals, Report ARD 99-7 2. Ando K., Ohta M., 1995, Relationships between the morphology of micro- fractures ofwood and the acoustic emission characteristics,Mokuzai Gakkaishi, 41, 640-646 3. Ando K., Sato K., Fushitani M., 1992, Fracture toughness and acoustic emission characteristics of wood II. Effects of grain angle,Mokuzai Gakkaishi, 38, 342-349 4. Atack D., May W.D., Morris E.L., Sproule R.N., 1961, The energy of tensile and cleavage fracture of black spruce,Tappi, 44, 555-567 5. Bodig J., Goodman J.R., 1973, Prediction of elastic parameters for wood, Wood Science, 5, 249-264 750 P.G. Kossakowski 6. Bowie O.L., Freese C.E., 1972, Central crack in plane orthotropic rectan- gular sheet, Journal of Fracture Mechanics, 8, 49-58 7. Cook T.S., Rau C.A. Jr., 1974, A critical rewiev of anisotropic fracture mechanics, In:Prospects of FractureMechanics, Sih, vanElst, andBroek,Edit., 509-523 8. DeBaise G.R., Porter A.W., PentoneyR.E., 1966,Morphology andme- chanics of wood fracture,Materials Research and Standards, 6, 493-499 9. GhandiK.R., 1972,Analysis of an inclined crack centrally placed in an ortho- tropic rectangular plate, Journal of Strain Analysis, 7, 157-163 10. HaberR.B., KohH.M., 1985,Explicit Expressions for EnergyReleaseRates Using Virtual Crack Extensions, International Journal for Numerical Methods in Engineering, 21, 301-315 11. Hellen T.K., 1975, On the method of virtual crack extensions, International Journal of Numerical Methods in Engineering, 9, 187-207 12. JohnsonJ.A., 1973,Crack initiation inwoodplates,Wood Science,6, 151-158 13. Kossakowski P.G., 2004,An Analysis of Mixed Mode Fracture Toughness of Pinewood Beam Elements, PhD. Thesis, Faculty of Civil and Environmental Engineering, Kielce University of Technology, Kielce 14. Mandel J.F., McGarry F.J., Wann S.S., Im J., 1974, Stress intensity factors for anisotropic fracture test specimens of several geometries, Journal of Composite Materials, 8, 106-116 15. Nikishkov G.P., Vaynshtok V.A., 1980, Metod virtualnogo rosta treshiny dlja otredelenija koeffitsientov intensivnosti naprijazenij KI i KII, Problemy Procznosti, 6, 26-30 16. ParksD.M., 1977, The virtual crack extensionmethod for nonlinearmaterial behavior, Computer Methods in Applied Mechanics and Engineering, 12, 353- 364 17. Patton-Mallory M., Cramer S.M., 1987, Fracture mechanics: a tool for predicting wood component strength,Forest Products Journal, 37, 39-47 18. Petterson R., Bodig J., 1983, Prediction of fracture toughness of conifers, Wood and Fiber Science, 15, 302-316 19. PN-EN 338:1999 Structural timber – Strength classes 20. PN-EN 380:1998 Timber structures – Test methods – General principles for static load testing 21. PN-EN 384:1999 Structural timber – Determination of characteristic values of mechanical properties and density 22. PN-EN 408:1998 Timber structures – Structural timber and glued laminated timber – Determination of some physical andmechanical properties Influence of anisotropy on the energy... 751 23. Porter A.W., 1964, On the mechanics of fracture in wood, Forest Products Journal, 14, 325-331 24. ReitererA., Stanzl-Tschegg S.E.,TscheggE.K., 2000,Mode I fracture and acoustic emission of softwood and hardwood,Wood Science and Technolo- gy, 34, 417-430 25. Seweryn A., 1998, Numerical methods for calculation of stress intensity factors, 17-th Symposium of Fatigue of Materials and Structures, Bydgoszcz- Pieczyska, 301-308 26. Schniewind A.P., Lyon D.E., 1973, A fracture mechanics approach to the tensile strength perpendicular to grain of dimension lumber,Wood Science and Technology, 7, 45-59 27. Schniewind A.P., Centeno J.C., 1973, Fracture toughness and duration of load factor I. Six principal systemsof crackpropagationand the duration factor for cracks propagating parallel to grain,Wood and Fiber, 5, 152-159 28. Schniewind A.P., 1977, Fracture toughness and duration of load factor II. Duration factor for crackspropagatingperpendicular-to-grain,Wood andFiber, 9, 216-226 29. Stanzl-Tschegg S.E., Tan D.M., Tschegg E.K., 1995, New splittingme- thod for wood fracture characterization, Wood Science and Technology, 29, 31-50 30. Stanzl-TscheggS.E.,TanD.M.,TscheggE.K., 1996,Fracture resistance to the crack propagation in wood, International Journal of Fracture, 75, 347- 356 31. Tomin M., 1971, Influence of wood orthotropy on basic equations of linear fracture mechanics,Drevarsky Vyskum, 16, 219-230 32. Tomin M., 1972, Influence of anisotropy on fracture toughness of wood,Wood Science, 5, 118-121 33. Walsh P.F., 1972, Linear fracture mechanics in orthotropic materials, Engi- neering Fracture Mechanics, 5, 533-541 34. Williams J.G., Birch M.W., 1976,Mixedmode fracture in anisotropic me- dia,Cracks and Fracture, ASTM 601, 125-137 35. WuE.M., 1967,Application of fracturemechanics to anisotropic plates, Jour- nal of Applied Mechanics, 34, 967-974 36. Yeh B., Schniewind A.P., 1992, Elasto-plastic fracture mechanics of wood using the J-integral method,Wood and Fiber Science, 24, 364-376 752 P.G. Kossakowski Wpływ anizotropii na współczynnik uwalniania energii GI dla materiałów wysokoortotropowych Streszczenie W artykule przedstawionowyniki analizy numerycznej dotyczącej współczynnika uwalniania energii GI dla materiałów wysokoortotropowych takich jak kompozyty, laminaty czy drewno.Współczynnik GI obliczano przy użyciu programuAdina v. 8.1 opartego na metodzie elementów skończonych (MES). Określając wpływ anizotro- pii na GI, przyjmowano różne modele materiałowe. Podczas analizy użyto dwuwy- miarowych (2D) i trójwymiarowych (3D) modeli numerycznych próbek podwójnie wspornikowych (ang. Double Cantilever Beam – DCB) o różnych proporcjach dłu- gości pęknięcia do grubości elementu. Zauważono, że im te proporcje są mniejsze, tym większe są różnice pomiędzy współczynnikami GI obliczanymi przy założeniu izotropowego i anizotropowych modeli materiałowych (transwersalnie izotropowych i ortotropowych). Dlatego też analiza pękania czy odporności na pękanie powinna być oparta namodelach transwersalnie izotropowych lub ortotropowych oraz powin- na być uwzględniana anizotropiamateriału. Manuscript received March 9, 2007; accepted for print April 25, 2007