- Email: [email protected]
Decomposed crack closure integrals for estimation of SIF variations
Engineering Fracture Mechanics 63 (1999) 165±178
Decomposed crack closure integrals for estimation of SIF variations Ripudaman Singh a,*, S.K. Patel b, B. Dattaguru c a
Karta Technologies, Inc., 1892 Grandstand, San Antonio, TX 78238, USA b Gas Turbine Research Establishment, Bangalore, 560017, India c Indian Institute of Science, Bangalore, 560012, India
Received 1 April 1998; received in revised form 20 December 1998; accepted 4 January 1999
Abstract The Modi®ed Crack Closure Integral is a popular, convenient and very eective post-processing tool for extracting the stress intensity factors (SIF) from ®nite element analysis (FEA) of cracked structures in individual or mixed-mode situations. In 3D, the current practice is to estimate the strain energy release rate over the brick element and mark the corresponding SIF to the center of the brick edge representing the crack front. This suppresses the information regarding SIF variations within the element. This paper presents a novel concept of decomposing the nodal forces, nodal displacements and the strain energy release rate into their constituent parts for extracting the SIF variations consistent with the element shape functions. Decomposed Crack Closure Integral (DCCI) has been formulated for linearly varying SIF in an 8-noded brick and parabolically varying SIF in a 20-noded brick element. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Finite element methods; Brick elements; 3D Crack closure integrals; Stress intensity factor variation
1. Introduction Damaged structures call for numerical analysis to estimate residual strength. Fatigue damage in thick structural parts normally manifests in the form of 3D cracks which may be through, part through or embedded. Accurate estimation of stress intensity factors along the crack front is important for reliable strength and fatigue life predictions and in turn for damage tolerance * Corresponding author. Tel.: +1-201-681-9102; fax: +1-201-681-9198. E-mail address: [email protected] (R. Singh) 0013-7944/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 3 - 7 9 4 4 ( 9 9 ) 0 0 0 3 1 - 4
166
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
analysis/design of high performance structures. Finite element analysis (FEA) of cracked structures with various fracture mechanics based post processor tools is being extensively used for such applications. Crack Closure Integral (CCI) as given by Irwin [1] has been a very useful tool for Strain Energy Release Rate (SERR) estimation and then Stress Intensity Factor (SIF) extraction at crack tip/front in engineering structures. The SERR can be estimated using CCI concept by considering an incremental crack extension and evaluating the work done to close the crack to original con®guration. Rybicki and Kanninen [2] formulated crack closure integral (MCCI) wherein, SERR can be estimated from the product of crack tip nodal force and crack opening displacement at ®rst node in crack tip wake, from single ®nite element analysis using four-noded quadratic elements. This simple formulation is element dependent as it was mathematically derived from the basic shape functions of the four-noded quadratic element. Over the last decade, various mathematical expressions have been developed for SERR estimation, involving nodal displacements and forces across the crack tip/front, each consistent with a speci®c element type in FEM [3±9]. Most of these formulations are summations of product of nodal forces and displacements, representing work done to close the crack over one element length/face. Closing the crack length/face displacements in various directions provides SERRs/SIFs in dierent modes. For brick elements in 3D FEA, the SERRs in three directions are estimated over the element face and the corresponding SIFs in three modes are attached to the center of the element edge coincident with the crack front [4,5]. This is exact if SIFs are constant along the crack front. However, for a 3D crack, the SIFs will generally vary along the crack front. Thus, the direct MCCI formulations will provide root mean square SIF values over the elemental width, at its center section. To capture the SIF variations along the crack front, one would essentially need a ®ner mesh. Narayana et al. [8,9] proposed a general procedure for derivation of MCCI expressions for 3D problems with cracks and mathematically con®rmed the MCCI formulae in practice. Based on their approach, they introduced a concept of subarea integration for brick elements to capture SERRs at several segments of the element. Once again, an RMS value of SIF over this segment is attached to the center point of the segment. The overall accuracy and variations within the element would now depend upon the number of segments per element. None of the CCI formulations to date can de®ne the variations in SIF. In 3D FEA, a varying SIF along the elemental edge representing the crack front, as an obvious consequence of varying stress/displacement ®elds, will show up as varying nodal displacements and forces. It should be possible to . Decompose the nodal forces and displacements to represent the constant and perturbing components of stress and displacement ®elds along a direction parallel to the crack front . Decompose the elemental strain energy release rate into constant and perturbing components . Formulate Decomposed Crack Closure Integral (DCCI) to estimate the constant and the perturbing components of the SIF, over elemental width The scope of this paper covers DCCI formulation for linear SIF variation in an eight-noded brick and quadratic SIF variation in a 20-noded brick element.
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
167
2. DCCI for an eight-noded brick element Consider the scenario as shown in Fig. 1a for a conventional eight-noded parallelepiped brick element. The forces (F1 and F2 from elements 1 and 2) and displacements (v1 and v2 from element 1), shown along `y'-direction, are the only ones that contribute to the Mode-I strain energy release rate [8] from the upper surface. We can see that the crack face is opening as well as rotating (if v1 $ v2). The total strain energy associated with crack closure will have two components: . Closure strain energy (Gclo) and . Rotational strain energy (Grot). The stress and displacement ®elds with the conventional shape functions for this element [8] are bilinear in x and z, on the faces in parametric space. They have the form s,v a0 a1 x a2 z a3 xz The nodal forces (F1 and F2) and displacements (v1 and v2) can be decomposed into two parts: . Mean: Fm and vm as obtained from stress and displacement state independent of z i.e. (a0+a1x ). (Refer to Fig. 1b), and . Alternating: Fa and va as obtained from stress and displacement state varying linearly with z i.e. z(a2+a3x ). (Refer to Fig. 1c.) For the eight-noded element, these force and displacement components are obtained as 1 1 vm v1 v2 va v1 ÿ v2 2 2
1 Fm F1 F2 2
1 Fa F1 ÿ F2 2
A similar set of forces (F'1 and F'2, not shown) and displacements (v1' and v'2, not shown) exist for the lower surface. Two components of the strain energy release rate associated with the two deformation constituents from both the surfaces are obtained as: Gclo
1 1 2Fm vm 2F 0 m v 0 m Grot 2Fa va 2F 0 a v 0 a 2Da 2Da
Force equilibrium demands that F'i=ÿFi for i=1, 2, m, a. Let us represent the gap between the nodes to be closed as Vi=viÿvi' for i=1, 2, m, a. The energy components become Gclo
1 F1 F2 V1 V2 4Da
Grot
1 F1 ÿ F2 V1 ÿ V2 4Da
Total strain energy released (Gtot=Gclo+Grot), turns out to be Gtot
1 F1 V1 F2 V2 2Da
168
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
Fig. 1. DCCI formulation from an 8-noded brick element.
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
169
This matches with the form of expression used by Shivakumar et al. [5] and Narayana et al. [8]. For an eight-noded brick, it is a reasonable assumption that the SIF varies linearly along the elemental edge coinciding with the crack front. This can also be expressed in terms of a mean (Km) and alternating (Ka) terms. For any location z from the mid-plane, we can write K z Km zKa Such that at the ends, z= 2 1, and K=Km 2 Ka. In fact, this is consistent with the element shape functions. From the well known SERR±SIF relationship, G
1 2 1 K K 2m 2zKm Ka z2 K 2a M M
where the material-state index `M' is a function of material properties (E, v ), the state of stress (plane stress or plane strain) and the mode of crack opening(-I, -II, -III). The total strain energy released, over the element can be obtained by integration, 1 1 1 1 2 21 2 K m K a Gclo Grot K dz Gtot M ÿ1 2 M 3 Equating the mean and alternating components, we get p p Km MGclo Ka 3MGrot Combining this with the expression of SERR components, DCCI formulation for an eightnoded brick element becomes: r 1 Km M F1 F2 V1 V2 4Da r 1 Ka 3M F1 ÿ F2 V1 ÿ V2 4Da The SIF values obtained by taking a square root of energy will have to be assigned a positive or a negative sign. The mean SIF Km is always positive, and the alternating SIF Ka will have the sign of (F1ÿF2). It may be noted that Ka cannot exist in isolation, as it would imply overlapping of elements when the crack is loaded. For a feasible solution the gaps have to be positive, i.e. Vir0 for i=1, 2. The above expressions correspond to Fig. 1, that is Mode-I. The subscript 'I' was conveniently dropped to make the expressions more readable. Mode -II, -III expressions would be identical, except that the directions of the nodal forces and displacements would be along `x' and `z'-axis, respectively, and the material state index `M' is speci®c to the fracture mode. The DCCI provides SIF variation along the elemental edge in all three modes. At any node, the two SIF values from the two adjacent elemental edges may not exactly match with each other, since the stresses are not necessarily continuous in FEA for this element. In such a
170
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
situation one might need to take an average SIF value based on contributions from the two neighboring elements on lines similar to those adopted for nodal stresses. The above equations for DCCI based SIF variations are also applicable to the six-noded brick element obtained from collapsing the eight-noded brick. The forces `F1 and F2', have to be obtained from contribution of all the elements on one side of the crack faces being virtually closed. 2.1. Numerical example Single sided patch repairs to cracked aluminum sheets are common to aerospace applications for the ease of carrying out in-situ external repairs to airframes. Such repairs demand
Fig. 2. Analysis of an unsymmetric patch repair.
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
171
geometrically non-linear analysis and produce a variation of SIF within the sheet thickness even for straight crack fronts. The variations are known to be close to linear for relatively thin sheets [10]. DCCI can be a useful tool to obtain SIF variations with as low as one element across the thickness of the sheet. Consider a strip of thin elastic aluminum sheet 200 40 1.5 mm, with an edge crack of length 15 mm. The repair is a 26 26 1.4 mm glass epoxy patch, bonded by 0.2 mm thick ®lm-epoxy adhesive, on one side of the sheet, as shown in Fig. 2a. The sheet is subjected to a remote uni-axial tensile stress of sy=58.33 MPa. Due to symmetry about x-axis, only one half needs to be analyzed. The adhesive is treated as cracked along with the sheet. The material properties are: aluminum E=70 GPa, n=0.32; ®lm adhesive E=2.2 GPa, n=0.32; glass-epoxy Es=38.6 and 8.27 GPa, ns=0.168 and 0.035, Gs=4.14 and 3.14 GPa. Non-linear ®nite element analysis with 8-noded 24-dof brick element is carried out using commercial code NISA on an IBM workstation. The element size at the crack tip is 0.5 0.5 mm. Multi-brick model with three, two and three bricks across the thickness of the cracked sheet, adhesive and patch, respectively, are considered to capture the thickness-wise variation in SIF using conventional MCCI and the proposed DCCI formulations. Single brick model with one brick each for the sheet, adhesive and patch is modeled to demonstrate the power of DCCI over MCCI. The two mesh con®gurations are schematically shown in Fig. 2b. Fig. 2c shows ÿ point SIF values at mid element locations obtained using MCCI (marked by `+' and `O'), and linearly varying SIF values using DCCI (marked by the dotted and solid lines) for both mesh resolutions used. The MCCI values are RMS average over the elemental width. Very clearly, the MCCI fails to capture the SIF variations across the width of the element. DCCI brings out these SIF variations. It also demonstrates that single element DCCI gives SIF variation as good as three element MCCI at almost 1/5 the computational cost. Numerically, the DCCI based nodal SIF values plotted in Fig. 2c are a mean of the SIF values obtained from two neighboring elements, due to the fore mentioned reason of stress discontinuity. The actual values dier by 3.8% and 1.8% from the mean value, at the two corresponding nodes. The detailed values are listed in Table 1, to make this section complete and exact. Table 1 SIF values corresponding to problem of Fig. 2 MCCI z Three bricks (eight-noded)
One brick (eight-noded)
DCCI KRMS
Km
Ka
z
K
0.00 0.50 0.50 1.00 1.00 1.50 0.00 1.50
4.34 8.34 7.73 11.19 10.80 13.24 4.73 13.43
0.25
6.45
6.34
2.00
0.75
9.51
9.46
1.73
1.25
12.04
12.02
1.22
0.75
9.42
9.08
4.35
172
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
Fig. 3. DCCI formulation from a 20-noded brick element.
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
173
3. DCCI for a 20-noded brick element Consider the scenario as shown in Fig. 3a for a conventional 20-noded parallelepiped brick element. The forces (F1, F2, F3, F4 and F5 from elements 1 and 2) and displacements (v1, v2, v3, v4 and v5 from element 1), shown along `y'-direction, are the only ones that contribute to the Mode-I strain energy release [9] from the upper surface. We can see that the crack face is opening, rotating as well as ¯exing. The total strain energy associated with crack closure will have three components: . Closure strain energy (Gclo), . Rotational strain energy (Grot), and . Flexural strain energy (G¯ex). The stress and displacement ®elds with the conventional shape functions for this element [9] are parabolic in x and z, on the faces in parametric space. They have the form s,v a0 a1 x a2 z a3 xz a4 x2 a5 z2 a6 x2 z a7 xz2 The nodal forces (F1, F2, F3, F4 and F5) and displacements (v1, v2, v3, v4 and v5) can be decomposed into three parts: . Mean: Fmo, Fm1, Fm2 and vm1, vm2 as obtained from stress and displacement state independent of z i.e. (a0+a1x+a4x 2). (Refer to Fig. 3b), . Alternating: Fa1, Fa2 and va1, va2 as obtained from stress and displacement state varying linearly with z i.e. z(a2+a3x+a6x 2). (Refer to Fig. 3c), and . Quadratic: Fq and vq as obtained from stress and displacement state varying parabolically with z i.e. z 2(a5+a7x ). (Refer to Fig. 3 d). The force and displacement components for the 20-noded element are obtained as 1 Fm0 2F1 2F2 2F3 F4 F5 3
1 vm1 4v1 v2 v3 6
1 1 Fm1 F1 F2 F3 ÿ F4 ÿ F5 vm2 2v1 ÿ v2 ÿ v3 3v4 3v5 6 6 1 1 Fm2 F4 F5 va1 v2 ÿ v3 2 2 1 1 Fa1 F2 ÿ F3 va2 v4 ÿ v5 2 2 1 1 Fa2 F4 ÿ F5 vq ÿ2v1 v2 v3 2 6
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
174
A similar set of forces (F'1, F'2, F'3, F'4 and F'5, not shown) and displacements (v1' , v2' , v3' , v4' and v5' , not shown) exist for the lower surface. Three components of the strain energy release rate associated with the three deformation constituents from both the surfaces are obtained as: Gclo
1 Fmo vm1 2Fm1 vm1 2Fm2 vm2 F 0 m0 v 0 m1 2F 0 m1 v 0 m1 2F 0 m2 v 0 m2 2Da
Grot
1 2Fa1 va1 2Fa2 va2 2F 0 a1 v 0 a1 2F 0 a2 v 0 a2 2Da
Gflex
1 6Fq vq 6F 0 q v 0 q 2Da
Force equilibrium demands that F'i=ÿFi for i=1, 2, 3, 4, 5, m0, m1, m2, a1, a2, q. Let us represent the gap between the nodes to be closed as Vi=viÿv 'i for i=1, 2, 3, 4, 5, m1, m2, a1, a2, q. The energy components become Gclo
1 F1 F2 F3 4V1 V2 V3 F4 F5 2V1 ÿ V2 ÿ V3 3V4 3V5 6Da
Grot
1 F2 ÿ F3 V2 ÿ V3 F4 ÿ F5 V4 ÿ V5 2Da
Gflex
1 ÿF1 2F2 2F3 F4 F5 ÿ2V1 V2 V3 6Da
Total strain energy released (Gtot=Gclo+Grot+G¯ex), turns out to be Gtot
1 F1 V1 F2 V2 F3 V3 F4 V4 F5 V5 2Da
This matches with the form of expression used by Shivakumar et al. [5] and Narayana et al. [8]. For a 20-noded brick, it is a reasonable assumption that the SIF varies quadratically along the elemental edge coinciding with the crack front. This can also be expressed in terms of a mean (Km), alternating (Ka) and quadratic (Kq) terms. For any location z from the mid-plane, we can write K z Km zKa 3z2 ÿ 1Kq Such that at the ends, z=21, and K=Km2Ka+2 Kq, and at the center, z=0, K=KmÿKq. In fact this is consistent with element shape functions. The total strain energy released over the element can once again be obtained by integration, 1 1 1 1 2 4 2 21 2 K m K a K q Gclo Grot Gflex K dz Gtot M ÿ1 2 M 3 5
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
175
Equating the mean, alternating and ¯exural components, we get r p p 5 MGflex Km MGclo Ka 3MGrot Kq 4 Combining this with the expression for SERR components, DCCI formulation for a 20-noded brick comes out to be: r 1 Km M F1 F2 F3 4V1 V2 V3 F4 F5 2V1 ÿ V2 ÿ V3 3V4 3V5 6Da
Fig. 4. Analysis of a symmetric patch repair.
176
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
r 1 F2 ÿ F3 V2 ÿ V3 F4 ÿ F5 V4 ÿ V5 Ka 3M 2Da
r 5 1 Kq M ÿF1 2F2 2F3 F4 F5 ÿ2V1 V2 V3 4 6Da The SIF values obtained by taking a square root of energy will have to be assigned a positive or a negative sign. The mean SIF Km is always positive, the alternating SIF Ka will have the sign of Fa1, and the ¯exural SIF Kq will have a sign of Fq. It may be noted that Ka and Kq cannot exist in isolation, as they would imply overlapping of elements when the crack is loaded. For a feasible solution the gaps have to be positive, i.e. Vir0 for i=1, 2, 3, 4, 5. The above expressions correspond to Fig. 3, that is Mode-I. The subscript `I' was conveniently dropped to make the expressions more readable. Mode-II, -III expressions would be identical, except that the directions of the nodal forces and displacements would be along `x' and `z'-axis respectively, and the material state index `M' is speci®c to the fracture mode. The DCCI provides SIF variation along the elemental edge in all three modes. At any node, the two SIF values from the two adjacent elemental edges may not exactly match with each other, since the stresses are not necessarily continuous in FEA for this element. In such a situation one might need to take an average SIF value based on contributions from the two neighboring elements on lines similar to those adopted for stresses. The above equations for DCCI based SIF variations are also applicable to the 15-noded brick element obtained from collapsing the 20-noded brick. The forces `F1, F2, F3, F4 and F5', have to be obtained from contribution of all the elements on one side of the crack faces being virtually closed. 3.1. Numerical example Symmetric patch repair to a thin cracked sheet produces a symmetric variation of SIF within the sheet thickness [11]. DCCI can be a useful tool to get a parabolic approximation to SIF variations with as low as single 20-noded brick element across the thickness of the sheet. Consider a change in repair con®guration of Fig. 2a from a single 1.4 mm patch to two symmetric patches 0.7 mm thick each, as shown in Fig. 4a. The element size at the crack tip is retained as 0.5 0.5 mm. Multi-brick model with 3, 2 and 3 8-noded and 20-noded bricks across the thickness of the cracked sheet, adhesive and patch respectively, are considered to capture the thickness-wise variations of SIF within the sheet. The con®guration is modeled with a single 20-noded brick each for the sheet, adhesive and patch, to demonstrate the power of DCCI. The two mesh con®gurations are schematically shown in Fig. 4b. Fig. 4c shows ÿ point SIF values at mid element locations obtained using MCCI (marked by `+', `O' and `D'), and parabolically varying SIF values using DCCI (marked by the dotted, dashed and solid lines), for both mesh resolutions used. The MCCI values are RMS average over the elemental width. Very clearly, the MCCI fails to capture the SIF variations across the width of the element. DCCI brings out these SIF variations. It also demonstrates that single element DCCI gives SIF variation as good as three element MCCI approximation at almost 1/ 3 the computational cost. Numerically, the DCCI based nodal SIF values plotted in Fig. 4c are
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
177
a mean of the SIF values obtained from two neighboring elements, due to the fore mentioned reason of stress discontinuity. The values dier by 0.6% from the mean value, for both 8- and 20-noded brick elements. The detailed values are listed in Table 2 to make this section complete and exact. In fact the single 20-noded brick DCCI based solution appears to be the most superior of all. It may be noted here that a 20-noded brick, which is believed to provide poor stress recovery, appears to behave very well, when it comes to energy based SIF estimation. 4. Closing remarks The application of MCCI with 3D FEA using brick elements need not be restricted to evaluation of total elemental strain energy release rate for an RMS average SIF estimate at the elemental edge center. The nodal forces and displacements along each of the three orthogonal directions can be conveniently decomposed to extract the variations in three modes of SIF along the elemental edge (crack front). The decomposed crack closure integral (DCCI) formulated in this paper can bring out linear variation of SIF from a 6/8-noded brick and a parabolic SIF variation from a 15/20-noded brick element with as much ease. This can help reduce the ®nite element mesh size, particularly for problems where SIF variations along the crack front are of signi®cance and large mesh sizes are computationally expensive. Acknowledgements Discussions with Professor T. S. Ramamurthy and Dr K. N. Raju at the Indian Institute of Science were helpful in streamlining some of the thoughts. Special thanks are due to Mr Turaga Umamaheswar for his help in computing and Mr Bhavik Jhaveri with type setting.
Table 2 SIF values corresponding to problem of Fig. 4 MCCI
Three bricks (8-noded)
Three bricks (20-noded)
One brick (20-noded)
DCCI
z
KRMS
Km
Ka
0.25
3.14
3.14
0.14
±
0.75
3.26
3.26
0.00
±
0.25
3.12
3.12
0.16
ÿ0.035
0.75
3.26
3.26
0.00
ÿ0.016
0.75
3.15
3.15
0.00
ÿ0.12
Kq
z
K
0.00 0.50 0.50 0.75 0.00 0.50 0.50 0.75 0.00 0.75
3.00 3.28 3.26 3.26 2.89 3.21 3.23 3.27 2.91 3.27
178
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
References [1] Broek D. Elementary engineering fracture mechanics. Nordo, 1974. [2] Rybicki EF, Kanninen MF. A ®nite element calculation of stress intensity factors by a modi®ed crack closure integral. Engng Frac Mech 1977;9:931±8. [3] Buchholz F. Improved formulae for the FE calculation of strain energy release rate by the modi®ed crack closure integral method. In: Proceedings of the Fourth World Congress and Exhibition in Finite Element Methods, Interlaken, 1984. p. 650±9. [4] Raju IS. Simple formulas for strain energy release rates with higher order and singular ®nite elements, NASA CR 178186 1986. [5] Shivakumar KN, Tan PW, Newman JC. A virtual crack closure technique for calculating stress intensity factors for cracked three dimensional bodies. Int J Frac 1988;56:R43±R50. [6] Ramamurthy TS, Krishnamurthy T, Narayana KB, Vijaykumar K, Dattaguru B. Modi®ed crack closure integral method with quarter point elements. Mech Res Commun 1986;13:179±86. [7] Viswanath S, Lakshminarayana HV, Ravindranath BR. A modi®ed crack closure integral method for calculating stress intensity factors for cracked plates subjected to bending. Int J Frac 1989;41:R45±R50. [8] Narayana KB, Dattaguru B, Ramamurthy TS, Vijaykumar K. A general procedure for modi®ed crack closure integral in 3D problems with cracks. Engng Frac Mech 1994;48:167±76. [9] Narayana KB, George S, Dattaguru B, Ramamurthy TS, Vijaykumar K. Modi®ed crack closure integral for 3D problems using 20-noded brick elements. Fatigue Frac Engng Mater Struct 1994;17:145±57. [10] Sun CT, Klug J, Arendt C. Analysis of cracked plates repaired with bonded composite patches. AIAA Journal 1996;34:369±74. [11] Kumar AM, Ripudaman Singh. 3D ®nite element modeling of a composite patch repair. In: Proceedings of the International Conference on Fracture, ICF-9, (April), 1997. p. 2159±66.
Decomposed crack closure integrals for estimation of SIF variations Ripudaman Singh a,*, S.K. Patel b, B. Dattaguru c a
Karta Technologies, Inc., 1892 Grandstand, San Antonio, TX 78238, USA b Gas Turbine Research Establishment, Bangalore, 560017, India c Indian Institute of Science, Bangalore, 560012, India
Received 1 April 1998; received in revised form 20 December 1998; accepted 4 January 1999
Abstract The Modi®ed Crack Closure Integral is a popular, convenient and very eective post-processing tool for extracting the stress intensity factors (SIF) from ®nite element analysis (FEA) of cracked structures in individual or mixed-mode situations. In 3D, the current practice is to estimate the strain energy release rate over the brick element and mark the corresponding SIF to the center of the brick edge representing the crack front. This suppresses the information regarding SIF variations within the element. This paper presents a novel concept of decomposing the nodal forces, nodal displacements and the strain energy release rate into their constituent parts for extracting the SIF variations consistent with the element shape functions. Decomposed Crack Closure Integral (DCCI) has been formulated for linearly varying SIF in an 8-noded brick and parabolically varying SIF in a 20-noded brick element. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Finite element methods; Brick elements; 3D Crack closure integrals; Stress intensity factor variation
1. Introduction Damaged structures call for numerical analysis to estimate residual strength. Fatigue damage in thick structural parts normally manifests in the form of 3D cracks which may be through, part through or embedded. Accurate estimation of stress intensity factors along the crack front is important for reliable strength and fatigue life predictions and in turn for damage tolerance * Corresponding author. Tel.: +1-201-681-9102; fax: +1-201-681-9198. E-mail address: [email protected] (R. Singh) 0013-7944/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 3 - 7 9 4 4 ( 9 9 ) 0 0 0 3 1 - 4
166
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
analysis/design of high performance structures. Finite element analysis (FEA) of cracked structures with various fracture mechanics based post processor tools is being extensively used for such applications. Crack Closure Integral (CCI) as given by Irwin [1] has been a very useful tool for Strain Energy Release Rate (SERR) estimation and then Stress Intensity Factor (SIF) extraction at crack tip/front in engineering structures. The SERR can be estimated using CCI concept by considering an incremental crack extension and evaluating the work done to close the crack to original con®guration. Rybicki and Kanninen [2] formulated crack closure integral (MCCI) wherein, SERR can be estimated from the product of crack tip nodal force and crack opening displacement at ®rst node in crack tip wake, from single ®nite element analysis using four-noded quadratic elements. This simple formulation is element dependent as it was mathematically derived from the basic shape functions of the four-noded quadratic element. Over the last decade, various mathematical expressions have been developed for SERR estimation, involving nodal displacements and forces across the crack tip/front, each consistent with a speci®c element type in FEM [3±9]. Most of these formulations are summations of product of nodal forces and displacements, representing work done to close the crack over one element length/face. Closing the crack length/face displacements in various directions provides SERRs/SIFs in dierent modes. For brick elements in 3D FEA, the SERRs in three directions are estimated over the element face and the corresponding SIFs in three modes are attached to the center of the element edge coincident with the crack front [4,5]. This is exact if SIFs are constant along the crack front. However, for a 3D crack, the SIFs will generally vary along the crack front. Thus, the direct MCCI formulations will provide root mean square SIF values over the elemental width, at its center section. To capture the SIF variations along the crack front, one would essentially need a ®ner mesh. Narayana et al. [8,9] proposed a general procedure for derivation of MCCI expressions for 3D problems with cracks and mathematically con®rmed the MCCI formulae in practice. Based on their approach, they introduced a concept of subarea integration for brick elements to capture SERRs at several segments of the element. Once again, an RMS value of SIF over this segment is attached to the center point of the segment. The overall accuracy and variations within the element would now depend upon the number of segments per element. None of the CCI formulations to date can de®ne the variations in SIF. In 3D FEA, a varying SIF along the elemental edge representing the crack front, as an obvious consequence of varying stress/displacement ®elds, will show up as varying nodal displacements and forces. It should be possible to . Decompose the nodal forces and displacements to represent the constant and perturbing components of stress and displacement ®elds along a direction parallel to the crack front . Decompose the elemental strain energy release rate into constant and perturbing components . Formulate Decomposed Crack Closure Integral (DCCI) to estimate the constant and the perturbing components of the SIF, over elemental width The scope of this paper covers DCCI formulation for linear SIF variation in an eight-noded brick and quadratic SIF variation in a 20-noded brick element.
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
167
2. DCCI for an eight-noded brick element Consider the scenario as shown in Fig. 1a for a conventional eight-noded parallelepiped brick element. The forces (F1 and F2 from elements 1 and 2) and displacements (v1 and v2 from element 1), shown along `y'-direction, are the only ones that contribute to the Mode-I strain energy release rate [8] from the upper surface. We can see that the crack face is opening as well as rotating (if v1 $ v2). The total strain energy associated with crack closure will have two components: . Closure strain energy (Gclo) and . Rotational strain energy (Grot). The stress and displacement ®elds with the conventional shape functions for this element [8] are bilinear in x and z, on the faces in parametric space. They have the form s,v a0 a1 x a2 z a3 xz The nodal forces (F1 and F2) and displacements (v1 and v2) can be decomposed into two parts: . Mean: Fm and vm as obtained from stress and displacement state independent of z i.e. (a0+a1x ). (Refer to Fig. 1b), and . Alternating: Fa and va as obtained from stress and displacement state varying linearly with z i.e. z(a2+a3x ). (Refer to Fig. 1c.) For the eight-noded element, these force and displacement components are obtained as 1 1 vm v1 v2 va v1 ÿ v2 2 2
1 Fm F1 F2 2
1 Fa F1 ÿ F2 2
A similar set of forces (F'1 and F'2, not shown) and displacements (v1' and v'2, not shown) exist for the lower surface. Two components of the strain energy release rate associated with the two deformation constituents from both the surfaces are obtained as: Gclo
1 1 2Fm vm 2F 0 m v 0 m Grot 2Fa va 2F 0 a v 0 a 2Da 2Da
Force equilibrium demands that F'i=ÿFi for i=1, 2, m, a. Let us represent the gap between the nodes to be closed as Vi=viÿvi' for i=1, 2, m, a. The energy components become Gclo
1 F1 F2 V1 V2 4Da
Grot
1 F1 ÿ F2 V1 ÿ V2 4Da
Total strain energy released (Gtot=Gclo+Grot), turns out to be Gtot
1 F1 V1 F2 V2 2Da
168
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
Fig. 1. DCCI formulation from an 8-noded brick element.
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
169
This matches with the form of expression used by Shivakumar et al. [5] and Narayana et al. [8]. For an eight-noded brick, it is a reasonable assumption that the SIF varies linearly along the elemental edge coinciding with the crack front. This can also be expressed in terms of a mean (Km) and alternating (Ka) terms. For any location z from the mid-plane, we can write K z Km zKa Such that at the ends, z= 2 1, and K=Km 2 Ka. In fact, this is consistent with the element shape functions. From the well known SERR±SIF relationship, G
1 2 1 K K 2m 2zKm Ka z2 K 2a M M
where the material-state index `M' is a function of material properties (E, v ), the state of stress (plane stress or plane strain) and the mode of crack opening(-I, -II, -III). The total strain energy released, over the element can be obtained by integration, 1 1 1 1 2 21 2 K m K a Gclo Grot K dz Gtot M ÿ1 2 M 3 Equating the mean and alternating components, we get p p Km MGclo Ka 3MGrot Combining this with the expression of SERR components, DCCI formulation for an eightnoded brick element becomes: r 1 Km M F1 F2 V1 V2 4Da r 1 Ka 3M F1 ÿ F2 V1 ÿ V2 4Da The SIF values obtained by taking a square root of energy will have to be assigned a positive or a negative sign. The mean SIF Km is always positive, and the alternating SIF Ka will have the sign of (F1ÿF2). It may be noted that Ka cannot exist in isolation, as it would imply overlapping of elements when the crack is loaded. For a feasible solution the gaps have to be positive, i.e. Vir0 for i=1, 2. The above expressions correspond to Fig. 1, that is Mode-I. The subscript 'I' was conveniently dropped to make the expressions more readable. Mode -II, -III expressions would be identical, except that the directions of the nodal forces and displacements would be along `x' and `z'-axis, respectively, and the material state index `M' is speci®c to the fracture mode. The DCCI provides SIF variation along the elemental edge in all three modes. At any node, the two SIF values from the two adjacent elemental edges may not exactly match with each other, since the stresses are not necessarily continuous in FEA for this element. In such a
170
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
situation one might need to take an average SIF value based on contributions from the two neighboring elements on lines similar to those adopted for nodal stresses. The above equations for DCCI based SIF variations are also applicable to the six-noded brick element obtained from collapsing the eight-noded brick. The forces `F1 and F2', have to be obtained from contribution of all the elements on one side of the crack faces being virtually closed. 2.1. Numerical example Single sided patch repairs to cracked aluminum sheets are common to aerospace applications for the ease of carrying out in-situ external repairs to airframes. Such repairs demand
Fig. 2. Analysis of an unsymmetric patch repair.
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
171
geometrically non-linear analysis and produce a variation of SIF within the sheet thickness even for straight crack fronts. The variations are known to be close to linear for relatively thin sheets [10]. DCCI can be a useful tool to obtain SIF variations with as low as one element across the thickness of the sheet. Consider a strip of thin elastic aluminum sheet 200 40 1.5 mm, with an edge crack of length 15 mm. The repair is a 26 26 1.4 mm glass epoxy patch, bonded by 0.2 mm thick ®lm-epoxy adhesive, on one side of the sheet, as shown in Fig. 2a. The sheet is subjected to a remote uni-axial tensile stress of sy=58.33 MPa. Due to symmetry about x-axis, only one half needs to be analyzed. The adhesive is treated as cracked along with the sheet. The material properties are: aluminum E=70 GPa, n=0.32; ®lm adhesive E=2.2 GPa, n=0.32; glass-epoxy Es=38.6 and 8.27 GPa, ns=0.168 and 0.035, Gs=4.14 and 3.14 GPa. Non-linear ®nite element analysis with 8-noded 24-dof brick element is carried out using commercial code NISA on an IBM workstation. The element size at the crack tip is 0.5 0.5 mm. Multi-brick model with three, two and three bricks across the thickness of the cracked sheet, adhesive and patch, respectively, are considered to capture the thickness-wise variation in SIF using conventional MCCI and the proposed DCCI formulations. Single brick model with one brick each for the sheet, adhesive and patch is modeled to demonstrate the power of DCCI over MCCI. The two mesh con®gurations are schematically shown in Fig. 2b. Fig. 2c shows ÿ point SIF values at mid element locations obtained using MCCI (marked by `+' and `O'), and linearly varying SIF values using DCCI (marked by the dotted and solid lines) for both mesh resolutions used. The MCCI values are RMS average over the elemental width. Very clearly, the MCCI fails to capture the SIF variations across the width of the element. DCCI brings out these SIF variations. It also demonstrates that single element DCCI gives SIF variation as good as three element MCCI at almost 1/5 the computational cost. Numerically, the DCCI based nodal SIF values plotted in Fig. 2c are a mean of the SIF values obtained from two neighboring elements, due to the fore mentioned reason of stress discontinuity. The actual values dier by 3.8% and 1.8% from the mean value, at the two corresponding nodes. The detailed values are listed in Table 1, to make this section complete and exact. Table 1 SIF values corresponding to problem of Fig. 2 MCCI z Three bricks (eight-noded)
One brick (eight-noded)
DCCI KRMS
Km
Ka
z
K
0.00 0.50 0.50 1.00 1.00 1.50 0.00 1.50
4.34 8.34 7.73 11.19 10.80 13.24 4.73 13.43
0.25
6.45
6.34
2.00
0.75
9.51
9.46
1.73
1.25
12.04
12.02
1.22
0.75
9.42
9.08
4.35
172
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
Fig. 3. DCCI formulation from a 20-noded brick element.
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
173
3. DCCI for a 20-noded brick element Consider the scenario as shown in Fig. 3a for a conventional 20-noded parallelepiped brick element. The forces (F1, F2, F3, F4 and F5 from elements 1 and 2) and displacements (v1, v2, v3, v4 and v5 from element 1), shown along `y'-direction, are the only ones that contribute to the Mode-I strain energy release [9] from the upper surface. We can see that the crack face is opening, rotating as well as ¯exing. The total strain energy associated with crack closure will have three components: . Closure strain energy (Gclo), . Rotational strain energy (Grot), and . Flexural strain energy (G¯ex). The stress and displacement ®elds with the conventional shape functions for this element [9] are parabolic in x and z, on the faces in parametric space. They have the form s,v a0 a1 x a2 z a3 xz a4 x2 a5 z2 a6 x2 z a7 xz2 The nodal forces (F1, F2, F3, F4 and F5) and displacements (v1, v2, v3, v4 and v5) can be decomposed into three parts: . Mean: Fmo, Fm1, Fm2 and vm1, vm2 as obtained from stress and displacement state independent of z i.e. (a0+a1x+a4x 2). (Refer to Fig. 3b), . Alternating: Fa1, Fa2 and va1, va2 as obtained from stress and displacement state varying linearly with z i.e. z(a2+a3x+a6x 2). (Refer to Fig. 3c), and . Quadratic: Fq and vq as obtained from stress and displacement state varying parabolically with z i.e. z 2(a5+a7x ). (Refer to Fig. 3 d). The force and displacement components for the 20-noded element are obtained as 1 Fm0 2F1 2F2 2F3 F4 F5 3
1 vm1 4v1 v2 v3 6
1 1 Fm1 F1 F2 F3 ÿ F4 ÿ F5 vm2 2v1 ÿ v2 ÿ v3 3v4 3v5 6 6 1 1 Fm2 F4 F5 va1 v2 ÿ v3 2 2 1 1 Fa1 F2 ÿ F3 va2 v4 ÿ v5 2 2 1 1 Fa2 F4 ÿ F5 vq ÿ2v1 v2 v3 2 6
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
174
A similar set of forces (F'1, F'2, F'3, F'4 and F'5, not shown) and displacements (v1' , v2' , v3' , v4' and v5' , not shown) exist for the lower surface. Three components of the strain energy release rate associated with the three deformation constituents from both the surfaces are obtained as: Gclo
1 Fmo vm1 2Fm1 vm1 2Fm2 vm2 F 0 m0 v 0 m1 2F 0 m1 v 0 m1 2F 0 m2 v 0 m2 2Da
Grot
1 2Fa1 va1 2Fa2 va2 2F 0 a1 v 0 a1 2F 0 a2 v 0 a2 2Da
Gflex
1 6Fq vq 6F 0 q v 0 q 2Da
Force equilibrium demands that F'i=ÿFi for i=1, 2, 3, 4, 5, m0, m1, m2, a1, a2, q. Let us represent the gap between the nodes to be closed as Vi=viÿv 'i for i=1, 2, 3, 4, 5, m1, m2, a1, a2, q. The energy components become Gclo
1 F1 F2 F3 4V1 V2 V3 F4 F5 2V1 ÿ V2 ÿ V3 3V4 3V5 6Da
Grot
1 F2 ÿ F3 V2 ÿ V3 F4 ÿ F5 V4 ÿ V5 2Da
Gflex
1 ÿF1 2F2 2F3 F4 F5 ÿ2V1 V2 V3 6Da
Total strain energy released (Gtot=Gclo+Grot+G¯ex), turns out to be Gtot
1 F1 V1 F2 V2 F3 V3 F4 V4 F5 V5 2Da
This matches with the form of expression used by Shivakumar et al. [5] and Narayana et al. [8]. For a 20-noded brick, it is a reasonable assumption that the SIF varies quadratically along the elemental edge coinciding with the crack front. This can also be expressed in terms of a mean (Km), alternating (Ka) and quadratic (Kq) terms. For any location z from the mid-plane, we can write K z Km zKa 3z2 ÿ 1Kq Such that at the ends, z=21, and K=Km2Ka+2 Kq, and at the center, z=0, K=KmÿKq. In fact this is consistent with element shape functions. The total strain energy released over the element can once again be obtained by integration, 1 1 1 1 2 4 2 21 2 K m K a K q Gclo Grot Gflex K dz Gtot M ÿ1 2 M 3 5
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
175
Equating the mean, alternating and ¯exural components, we get r p p 5 MGflex Km MGclo Ka 3MGrot Kq 4 Combining this with the expression for SERR components, DCCI formulation for a 20-noded brick comes out to be: r 1 Km M F1 F2 F3 4V1 V2 V3 F4 F5 2V1 ÿ V2 ÿ V3 3V4 3V5 6Da
Fig. 4. Analysis of a symmetric patch repair.
176
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
r 1 F2 ÿ F3 V2 ÿ V3 F4 ÿ F5 V4 ÿ V5 Ka 3M 2Da
r 5 1 Kq M ÿF1 2F2 2F3 F4 F5 ÿ2V1 V2 V3 4 6Da The SIF values obtained by taking a square root of energy will have to be assigned a positive or a negative sign. The mean SIF Km is always positive, the alternating SIF Ka will have the sign of Fa1, and the ¯exural SIF Kq will have a sign of Fq. It may be noted that Ka and Kq cannot exist in isolation, as they would imply overlapping of elements when the crack is loaded. For a feasible solution the gaps have to be positive, i.e. Vir0 for i=1, 2, 3, 4, 5. The above expressions correspond to Fig. 3, that is Mode-I. The subscript `I' was conveniently dropped to make the expressions more readable. Mode-II, -III expressions would be identical, except that the directions of the nodal forces and displacements would be along `x' and `z'-axis respectively, and the material state index `M' is speci®c to the fracture mode. The DCCI provides SIF variation along the elemental edge in all three modes. At any node, the two SIF values from the two adjacent elemental edges may not exactly match with each other, since the stresses are not necessarily continuous in FEA for this element. In such a situation one might need to take an average SIF value based on contributions from the two neighboring elements on lines similar to those adopted for stresses. The above equations for DCCI based SIF variations are also applicable to the 15-noded brick element obtained from collapsing the 20-noded brick. The forces `F1, F2, F3, F4 and F5', have to be obtained from contribution of all the elements on one side of the crack faces being virtually closed. 3.1. Numerical example Symmetric patch repair to a thin cracked sheet produces a symmetric variation of SIF within the sheet thickness [11]. DCCI can be a useful tool to get a parabolic approximation to SIF variations with as low as single 20-noded brick element across the thickness of the sheet. Consider a change in repair con®guration of Fig. 2a from a single 1.4 mm patch to two symmetric patches 0.7 mm thick each, as shown in Fig. 4a. The element size at the crack tip is retained as 0.5 0.5 mm. Multi-brick model with 3, 2 and 3 8-noded and 20-noded bricks across the thickness of the cracked sheet, adhesive and patch respectively, are considered to capture the thickness-wise variations of SIF within the sheet. The con®guration is modeled with a single 20-noded brick each for the sheet, adhesive and patch, to demonstrate the power of DCCI. The two mesh con®gurations are schematically shown in Fig. 4b. Fig. 4c shows ÿ point SIF values at mid element locations obtained using MCCI (marked by `+', `O' and `D'), and parabolically varying SIF values using DCCI (marked by the dotted, dashed and solid lines), for both mesh resolutions used. The MCCI values are RMS average over the elemental width. Very clearly, the MCCI fails to capture the SIF variations across the width of the element. DCCI brings out these SIF variations. It also demonstrates that single element DCCI gives SIF variation as good as three element MCCI approximation at almost 1/ 3 the computational cost. Numerically, the DCCI based nodal SIF values plotted in Fig. 4c are
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
177
a mean of the SIF values obtained from two neighboring elements, due to the fore mentioned reason of stress discontinuity. The values dier by 0.6% from the mean value, for both 8- and 20-noded brick elements. The detailed values are listed in Table 2 to make this section complete and exact. In fact the single 20-noded brick DCCI based solution appears to be the most superior of all. It may be noted here that a 20-noded brick, which is believed to provide poor stress recovery, appears to behave very well, when it comes to energy based SIF estimation. 4. Closing remarks The application of MCCI with 3D FEA using brick elements need not be restricted to evaluation of total elemental strain energy release rate for an RMS average SIF estimate at the elemental edge center. The nodal forces and displacements along each of the three orthogonal directions can be conveniently decomposed to extract the variations in three modes of SIF along the elemental edge (crack front). The decomposed crack closure integral (DCCI) formulated in this paper can bring out linear variation of SIF from a 6/8-noded brick and a parabolic SIF variation from a 15/20-noded brick element with as much ease. This can help reduce the ®nite element mesh size, particularly for problems where SIF variations along the crack front are of signi®cance and large mesh sizes are computationally expensive. Acknowledgements Discussions with Professor T. S. Ramamurthy and Dr K. N. Raju at the Indian Institute of Science were helpful in streamlining some of the thoughts. Special thanks are due to Mr Turaga Umamaheswar for his help in computing and Mr Bhavik Jhaveri with type setting.
Table 2 SIF values corresponding to problem of Fig. 4 MCCI
Three bricks (8-noded)
Three bricks (20-noded)
One brick (20-noded)
DCCI
z
KRMS
Km
Ka
0.25
3.14
3.14
0.14
±
0.75
3.26
3.26
0.00
±
0.25
3.12
3.12
0.16
ÿ0.035
0.75
3.26
3.26
0.00
ÿ0.016
0.75
3.15
3.15
0.00
ÿ0.12
Kq
z
K
0.00 0.50 0.50 0.75 0.00 0.50 0.50 0.75 0.00 0.75
3.00 3.28 3.26 3.26 2.89 3.21 3.23 3.27 2.91 3.27
178
R. Singh et al. / Engineering Fracture Mechanics 63 (1999) 165±178
References [1] Broek D. Elementary engineering fracture mechanics. Nordo, 1974. [2] Rybicki EF, Kanninen MF. A ®nite element calculation of stress intensity factors by a modi®ed crack closure integral. Engng Frac Mech 1977;9:931±8. [3] Buchholz F. Improved formulae for the FE calculation of strain energy release rate by the modi®ed crack closure integral method. In: Proceedings of the Fourth World Congress and Exhibition in Finite Element Methods, Interlaken, 1984. p. 650±9. [4] Raju IS. Simple formulas for strain energy release rates with higher order and singular ®nite elements, NASA CR 178186 1986. [5] Shivakumar KN, Tan PW, Newman JC. A virtual crack closure technique for calculating stress intensity factors for cracked three dimensional bodies. Int J Frac 1988;56:R43±R50. [6] Ramamurthy TS, Krishnamurthy T, Narayana KB, Vijaykumar K, Dattaguru B. Modi®ed crack closure integral method with quarter point elements. Mech Res Commun 1986;13:179±86. [7] Viswanath S, Lakshminarayana HV, Ravindranath BR. A modi®ed crack closure integral method for calculating stress intensity factors for cracked plates subjected to bending. Int J Frac 1989;41:R45±R50. [8] Narayana KB, Dattaguru B, Ramamurthy TS, Vijaykumar K. A general procedure for modi®ed crack closure integral in 3D problems with cracks. Engng Frac Mech 1994;48:167±76. [9] Narayana KB, George S, Dattaguru B, Ramamurthy TS, Vijaykumar K. Modi®ed crack closure integral for 3D problems using 20-noded brick elements. Fatigue Frac Engng Mater Struct 1994;17:145±57. [10] Sun CT, Klug J, Arendt C. Analysis of cracked plates repaired with bonded composite patches. AIAA Journal 1996;34:369±74. [11] Kumar AM, Ripudaman Singh. 3D ®nite element modeling of a composite patch repair. In: Proceedings of the International Conference on Fracture, ICF-9, (April), 1997. p. 2159±66.