- Email: [email protected]
Three-dimensional fracture analysis in transversely isotropic solids
Engineering Analysis with Boundary Elements 20 (1997) 287-298 © 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain Plh S 0 9 5 5 - 7 9 9 7 ( 9 7 ) 0 0 0 7 0 - 2
0955-7997/97/$17.00
ELSEVIER
Three-dimensional fracture analysis in transversely isotropic solids A. Sfiez, M. P. Ariza & J. Dominguez* Escuela Superior de Ingenieros, Universidad de Sevilla, Av. Reina Mercedes s/n, 41012-Sevilla, Spain
(Received 21 August 1997; accepted 5 September 1997)
In this paper a boundary element formulation for three-dimensional crack problems in transversely isotropic bodies is presented. Quarter-point and singular quarter-point elements are implemented in a quadratic isoparametric element context. The point load fundamental solution for transversely isotropic media is implemented. Numerical solutions to several three-dimensional crack problems are obtained. The accuracy and robustness of the present approach for the analysis of fracture mechanics problems in transversely isotropic bodies are shown by comparison of some of the results obtained with existing analytical solutions. The approach is shown to be a simple and useful tool for the evaluation of stress intensity factors in transversely isotropic media. © 1998 Elsevier Science Ltd. All rights reserved Keywords: Cracks, fracture mechanics, stress intensity factor, boundary elements,
transversely isotropic bodies. 1 INTRODUCTION
dimensional anisotropic bodies 3 years previously. More recently, Rajiyah and Atluri 3 generalized the analytical solution for the SIF of a flat elliptical crack, to arbitrary crack face loading. The previous solution obtained by Kassir and Sih I was limited to the case of constant and linear variations of tractions on the crack face. Rajiyah and Atluri 3 employed their generalized solution with the finite element alternating technique to analyze embedded and surface elliptical cracks in transversely isotropic bodies of finite dimensions where the crack plane is perpendicular to the material axis. The number of boundary element papers dedicated to the analysis of fracture mechanics problems in non-isotropic materials is rather small. Among these papers, the early work of Snyder and C r u s e 4 and the more recent one of Chan and Cruse 5 should be mentioned. Cruse and co-workers developed a boundary element approach for the analysis of two-dimensional crack problems in fully anisotropic bodies. They employed a special fundamental solution which incorporates the presence of a traction-free crack in an infinite anisotropic body. The use of such fundamental solution makes it unnecessary to represent the crack surface in the boundary element discretization. The stress intensity factors are computed using a path independent integral once the boundary displacements and tractions have been obtained. In 1992, Tan and Gao 6 presented the quadratic
Transversely isotropic elastic materials are those with an axis of symmetry such that all directions perpendicular to this axis are equivalent. In other words, any plane perpendicular to the axis is a plane of isotropy. Crystals of the hexagonal system are transversely isotropic solids. There is also a significant number of fiber-reinforced composites which show this kind of behavior from a macroscopic point of view. They are composite solids made of unidirectionally oriented fibers, the fiber diameter and the fiber spacing being much smaller than the dimensions of the body. The basic fracture mechanics concepts for transversely isotropic solids were set 30 years ago by Kassir and Sih. I They showed that the stress singularity of the order ~r near the periphery of three-dimensional cracks, well known for isotropic materials, remains in the case of transversely isotropic solids. The magnitude of the local stresses may also be described in this case in terms of stress intensity factors (SIFs). Kassir and Sih obtained expressions of stresses and displacements near the crack front for cracks of arbitrary shape in a plane perpendicular to the material axis of symmetry. Sih et al. 2 had obtained expressions for stress and displacement fields near the tips of cracks in two*To whom correspondence should be addressed. 287
288
A. & i e : et al.
isoparametric formulation for fully anisotropic twodimensional problems. They used quarter-point singular traction crack-tip elements to compute stress intensity factors directly from the traction nodal values at the crack tip as done by Martfnez and Domfnguez7 for isotropic materials. The results obtained in that paper show great accuracy and little mesh dependency, All the previously mentioned boundary element papers deal with two-dimensional problems. To the authors' knowledge, the only three-dimensional boundary element formulation for fracture mechanics in non-isotropic materials was presented by lshikawa. ~ This author studied crack problems in three-dimensional transversely isotropic and fully anisotropic bodies. The approach is based on the subdomain technique and the use of quarter-point elements in combination with transition elements placed next to them. No singularity is introduced to represent the tractions at the crack front. The SIF values are computed from a twodimensional approximation of the crack opening displacements. The expression for the displacement near the crack tip in terms of the SIF in an anisotropic two-dimensional plate is used in planes normal to the crack front. Only results for mode-! SIF are given by Ishikawa. s In the present paper, the nine node singular quarter-point quadratic element for three-dimensional crack problems, as presented by the authors, 9 is extended to transversely isotropic materials. The crack displacement variation is represented by quarter-point elements which are obtained by locating the mid-side nodes of a quadratic element at one quarter of the element width. The singular tractions over the elements next to the crack front are represented by singular functions obtained by dividing the quadratic shape functions by the square root of the distance to the front. This fact constitutes one of the main differences with respect to Ishikawa's approach, s which contains a tractions representation unable to reproduce the actual singular traction distribution. The integration approach over the boundary elements used in the present paper for transversely isotropic media is the same as given in Ref. 9 for isotropic materials since the singularities of the fundamental solution and of the stress field in the proximity of the crack front are of the same type in both cases. Values of the SIF are computed as traction system unknowns at the crack front nodes. They are also computed by making the obtained values of the crack opening displacement at the quarter-point nodes inside the crack, equal to their analytical expression at the same points, which is a known function of the SIFs. The actual threedimensional expressions given by Kassir and Sih ~ are used. As a general rule, the SIF values obtained from crack opening displacements are less accurate and more sensitive to the discretization than those obtained from the traction crack-front nodal values. Nevertheless, in the studies presented in this paper, both approaches yield very accurate results. In any case, having a second approach for stress intensity factor computations provides a very useful crosscheck on the present numerical approach. Numerical examples with different materials and
geometries are presented in this paper. The computed values of the SIF are compared with other exact or approximate results when they exist. The accuracy of the results obtained tot three-dimensional crack problems in transversely isotropic materials is very good.
2 TRANSVERSELY ISOTROPIC MATERIALS
It is useful to review the basic relations governing the behavior of elastic transversely isotropic materials. The strain tensor is defined, as usual, in terms of displacement derivatives 1 8ii = ~(ui,/+ Ui.i)
(1)
The stress-strain relation for a general elastic anisotropic material can be written as G i = Cii~/ek:
(2)
where Ci/k ~ are the elastic constants. There are 81 possible constants in eqn (2) but only 21 are independent (see for instance Ref. m). When a material is transversely isotropic and the axis of symmetry coincides with a Cartesian axis, for instance x~, the stress-strain relation can be written in terms of five independent constants in the following matrix form
= ]
Cll
Cj2
Ct:~
0
0
0
Cj2
Cj]
Cj3
0
0
0
Cj3
Cr3
C~
0
0
0
0
0
0
C44
0
0
0
0
0
0
C44
0
0
0
0
0
0
Cij - Cj2 2
[
811 82 8~
(3) 28 28 28 where Cu are elastic constants in a condensed notation in terms of only two subindices. lsotropic materials can be considered as particular cases of transversely isotropic materials with
Cll =C33=X+2/.t; Ci2=Ci3=)k; C44~#
(4)
The stress-strain relation can be written in the inverse form as eli = aij,C#~kl
(5)
where a6k / are elastic compliances such that the matrix of
289
3-D fracture analysis in transversely isotropic solids elastic compliances is the inverse of the matrix of elastic constants. In the case of transversely isotropic materials, the matrix of elastic compliances may be written in terms of engineering material constants as follows II
'Ell
the origin of coordinates of an infinite transversely isotropic three-dimensional domain was obtained by Pan and Chou 12 in 1976. They represented the displacements in terms of three potential functions. Assuming that the x3-axis is the material axis of symmetry, the displacements at an arbitrary point x are as follows:
E22 3.2 U n i t p o i n t l o a d in the x3 d i r e c t i o n
E33 2E31
2e32
* ~{ 113o~~-
i 1 Eli
- vl2 Ei t
- v31 E33
0
0
0
- vl~
1
- v31
0
0
0
Ell
Ell
E33
-vt3 Etl
-vL3 Etl
1 E33
0
0
0
0
GI3
0
0
0
0
• 17-11
0
0
0
Gl3
1
0
0
0
0
0
0
-G12
O31
• O"12.
1
v~3 Ell
v3] and Gi2 E33
in the I direction, Gu is and vu is the Poisson's Only five of the seven are independent. They
Eli
(8)
Cllni + C44v2Ai'] 1
-1- v~(A i -k- Bi) ~ C44 R2 -q- Ci i x2
(9)
where ut*~stands for the displacement along the k direction due to a unit point load in the l direction, R 2 = x~ + x 2
R~ : R 2 + Zai Zi = ViX3
(6) where E//is the Young's modulus the shear modulus in the/J-plane, ratio relating directions / and J. engineering constants in eqn (6) satisfy the relations
ui(Aiq-Bi)XaZ--~--~) ~ : 1 , 2
022
0"32 0
2{(
"~3 = -- Z
0"33 1
viAixa, RiRi
Ri : Ri + zi va:{
( P I 3 - C 1 3 ) (1PCI 33 3/+C44 C I23 + 24C 4 4 ) }
- ( - l ) c ~ { (P13+CI3)(PI3" 4C3"3 C44-CI3- 2C44)} 1/2 c~= 1,2 (7)
2(1 +vl2) Cll - C12
v3 ~---
2C44
3 BOUNDARY E L E M E N T F O R M U L A T I O N The displacement integral representation and the boundary element discretization for transversely isotropic elastic bodies is the same as for the isotropic case. The only change in the boundary element formulation is that a new fundamental solution, corresponding to a point load in an infinite transversely isotropic domain, must be used. This change obviously transforms the expressions to be integrated over the boundary elements but does not alter the integration approach since the singularities appearing in the transversely isotropic fundamental solution are of the same kind as those in the isotropic fundamental solution.
PI3 : V/Cll C33 and A~ and Bi are given by:
1. PI3 -- C13 - 2C44 ~ 0
CI3+C44 vIA I : - v2A 2 : 47rC33C44(v2 - Vl) Bi = - A i
2. PI3-C]3-2C44=0 A] = A 2 = 0
3.1 F u n d a m e n t a l s o l u t i o n
B I = B 2 -The solution to the problem of a unit point load applied at
CI3+C44 167rCll C44
290
A. S S e : et al.
3.3 U n i t p o i n t load in the x~ direction
i:
the strain components and substi.tution of the stress-strain relation (3). Explicit expressions for the stress tensor components can be found in Pan and Chou's paper. J2 Tractions are immediately obtained by projection over any surface defined by its normal. The singularity contained in the fundamental solution tractions is of the type ( l/r 2) when r tends to zero: i.e., the same as lbr the isotropic case.
Rt /I
+ 2viBi
RiR2.i~2 } + D
110)
{ u~o :- =
-rl.r~ , -rlx~ ] - ui(A" ' - Bi') R~ ' - - 2uiBi RiR;" J
+ D - -X I X2
(II)
R3R*32
1t73 :
I ":i --U i2 li(A i , -B, , ) XR~
CIIA i -- C44u~B i 2.rl
4 T H R E E D I M E N S I O N A L C R A C K S IN TRANSVERSELY ISOTROPIC SOLIDS
According to the work of Kassir and Sih j the stress intensity factors K j (J = I, II, Ill) for cracks in a plane perpendicular to the material axis of symmetry can be defined in terms of the stress distribution in the usual way K I = lim ~ ( o - ) 0 :
}
r~[)
Cli/u[-C4a
CI3"~-C44 D--
1
4"xv:~(744
and A i' and B i' are given by:
I. PI3 Ci3 2C44 ~ 0 -
-
-
-
I ) ~ G 3 -- C44
r~0
KH, = lim ~/(2"n-r)(o,:)0=ll
(12) where
~: KH = lim " ~ ' ) ( ~ , , : ) 0 : , , ; ""
where the stress components refer to a natural local coordinate system with origin at the arbitrary point P on the crack front where the SIF is computed (Fig. 1 ). The n-axis is directed along the normal to the crack front; the t-axis along the tangent to the crack front; and the z-axis perpendicular to the crack plane. The polar coordinates r and 0 are defined in the nz-plane. Displacement and stress components in the proximity of the crack front for a crack in an infinite body loaded symmetrically with respect to the crack plane can be written LIS:
",,-
KI , ~ r / ~ ) ( ]
n/2
isin20)l/2lJ/~ :
A 2 =B~':
":
,h
I + m2
2. PI3 - CI3 - 2C44 = 0 A n' = A 3 _ r ~__
m
B i' =B.'
1
167rCll
[cos0+(cos-0+n;
Jsm-0)l/-] j/2
} +0(r) (14)
uf = 0(r)
1 • 167rC44vi
3.4 U n i t p o i n t load in the x z direction
The displacements are the same as given by eqns (10)-(12) when indices I and 2 are interchanged. The displacements due to a unit point load applied at the origin, given by eqns (8)-(12) are singular as (I/r) when distance r between the load and the observation point tends to zero. This singularity is of the same type as that appearing in the fundamental solution for isotropic media. The stress at any point can be easily derived by differentiation of the above displacement expressions to obtain
" *
ack Boundary
n "~
Fig. 1. Local coordinates at crack boundary point P.
3-D fracture analysis in transversely isotropic solids
291
KI
Kjsin 0 ( ~
G,z = 2(nli2 _ n~12)v/(r~r) { [(COS20 "1- n l -
ID.74 nl'(+ml)
.[cosO+(cos20+nllsin20)ll2]
I sin20)1/2 __ COS 0] 1/2
-I/2 (C0S20 +
n~ 1sin20) - t/2]ln
_ [(cos20 +//2
1sin20)l/2
1172
n~/2(l + m2)
-- COS 0] I/2(COS20 q- n 2 I sin20) - 1/2 } +
0(r o)
a,z = 0(r °)
lcos 0 + (cos20 + n2- t sin20)1/2] - I/2} _4_O(r)
~,,, = O(r°)
and
where n~ and n2 are the two solutions of the quadratic equation of the material properties
KI Unn --
C11C44 n2 + [C13(C13 +
2C44) -
C I I C 3 3 ] n + C33C44 = 0
Cl3ml ~ I----n C I J [cos 0 + (COS20 + n i- I sin 20) 1/2l I/2 (1 + ml)nll/2
(16) (Ci3 + C44)nj, j = l , 2
(cos20 + n ~ I sin20) - 1/2
C33 _ C44nj
mj--
and
CII n2 [cos I sin20) l/2] i/2 (1 -4- m2)n I/2 ~ 0 + (cos20 + n2-
CI 3m2 -
fit = C44(n lt/2 -- n~/2)
For a loading state which is skew-symmetric with respect to the crack plane, the displacement and stress components in the proximity of the crack front are: 1
( COS20 q- nZ Isine0) - I/2 } q_ O(r0) K!
u,,= -
ut t --
Ci3ml -- Ci2nl lc°s0+¢c°s20+n
(17)
-
Kll(n,n.,/,,~sln0(,r/--r~)l/2 ~ - - - ~ v "{[(1 + m , ) n , ] - ' C44(nl 2 _ nil )
[COS 0 q'- (C0S20 q-
l sin20)l/2]t/2
n 1 I sin20)1/2] - I/2 _ [(1 q- m2)n21 - I
[COS 0 + (COS20 "4- n 2
(cos20 + n 1-J sin20)- J/2 _ Cj 3rn2 - Ci2n2 (1 + m2)n I/2
Lsin20) I/2] -
1/2 } -4- O(r)
Kmn~nsin 0 v / ( r - ~ ) /l t
C44 (m2 - m I )
[cos 0 + (cos20 + n f I sin20)l/2] I/2
1 + m2[co s 0 + (cos20 + nl- I sin20) 1/2] - I/2 (cos20 + n2 t sin20) - t/2 } + 0(/i)
o-__--
KI C44
{nll/2[cos 0 + (cos20 + n I- I sin20)l/2] 1/2 (COS20 + n l- t sin20)- t/2 _ n~/2 [COS 0 + (COS 2 0 q- n 2 I sin 20)1/2 ] 1/2 • 2 -I/2 (COS2 0 + n 2- I sin 0) } + 0(r °)
//1
(15)
1 + ml [cos 0 + (cos20 + n f t sin20) In] - In ~ + 0(r) ) n2 (18)
Kll (Vq77~ oQC44(m I -- m2) ml(1n j/~ + m l ) r .Icos O + (cos 20 + n i sin20)ln] i/2
292
A. S a e z et al.
and Kmn~/2(Ci) - Ci
O,tI
KH
~l(cos20+ n f
I sine0)1:2 _ cos 0l i/2
2C44 k/('ar)
2oe2C44(m I -//12) ~(rlr) {(1 + m2)(Cl3m I - - CIIIll
)Hi
I[(cos20
q- I1[
(cos20 + n~ t sin'0) - .2 + O(r o)
I sin:O)I/e
where - - COS
0]1/2(COS20-~11! Isin20)
I/2
I1; -- (1 + m l ) ( C j 3 m 2 -- CI 172)n2 I - cos 0l I/-~
[(cos~.O+nflsin:O),/2
0~ (cos20 + n-, I sin20)
(19)
20Q C4a (m I - 1112) Ci2nl)n
I I
× [(cos20 + hi- I sin:O) I/2 - cos OIIn (cosRO+nl isin20)
I/2
- (1 + ml)(G3m2 - Ci2n2)n2 [(cos20 + n2 1sine0)1/2 _ cos 01 i/2
X (COS20+112 Isin20) I/2} + 0 ( r 0) K H(nj n 2)1/2
--2(n~/2_nln)xfOrr) {[cos20 + n ( I sin20)l/2 - cos 01": lcos20+nl --
sin'0)
i/2
[ ( c o s 2 0 + n 2 I s i•n - t,~,1/~ /) I
"~
•
n
(cos'0 + n2 sln~0) at:=
Kill
1/2}
- - COS
011/2
+ O(r u) - Isin:O)t:-II:-
2 X/(,a-r) ] • )
Ill]
I]
71-1}12) I1 n~/2
-- Ill')
,1~/~
5 DISCRETIZATION
The boundary element representation of three-dimensional crack problems in transversely isotropic bodies is done following the same steps as in fully isotropic media: (I) the domain under consideration is divided into subdomains by cutting it by a section through the crack; (2) each subdomain is analyzed by the BEM; and (3) the subdomains are coupled using equilibrium and compatibility relations. When the crack is in a plane of symmetry or skew-symmetry. the symmetry conditions can be directly applied over one subregion. The elements used are standard nine node or six node quadratic elements except for those with one side at the crack front. Among these elements, the ones inside the crack surface (non-singular traction) are quarter-point elements and those extending inside the material are singular quarter-point elements. These two types of elements have been used by the authors for similar problems in isotropic materials. 9 The quarter-point element (Fig. 2) is a quadrilateral and has its nine nodes on a plane. The displacement representation over the element when quadratic shape functions are used and the mid-side nodes are located at one quarter of the width, can be written as:
,,, =,,,!
,
~__tcosO+(cos20-kn~
(cos20 + n~- sin-0)
(l ~-/ll])(l
i/2 } + 0(r o)
Kn
{(1 q-m2)(Ci3m I -
2(744 (Clt - C ~ : )
i/2 + 0(r o)
t20/
where ? is the distance to the crack front, L the element 7
Kll (nt u2) J/2
6
o , , : = 2(n~/--- n l / 2 ) . ~ r r )
{n ~ I/2[cos 0 + (cos:'0 + n ~- sin 20) i/:, 1I/2 (cose0+n~ -Isin20)
J/:--n~ 1/2
[cos O + (cos20 + n2 Isin20)J/2 i 1/2
xl
3 4
5
I
(cos20 + n £ l sin20)
I/2 } +
O(r o) Fig. 2. Ninc node quarter-point quadratic elemcnt.
3-D fracture analysis in transversely isotropic solids width, and a'i" are polynomial functions of the nodal values of ui and the element coordinate ~l along the crack front. The displacement variation in this quarterpoint element is able to reproduce the displacement behavior in the proximity of the front of a crack in a transversely isotropic material, which according to eqns (14) and (18), have a variation of the type x/~ (notice that the distance to the crack front, r, used in the general equations, has been changed to ? in the BE formulation because r is used for the distance to the collocation point). The singular quarter-point elements have one side along the crack front and are part of the boundary extending inside the material. They are able to represent the (l/x/if) behavior of the stress near the crack front given by eqns (15) and (19) and have the same geometry and displacement representation as the quarter-point elements. In the tractions representation, the standard shape functions are divided by I/X/~/L. Thus, variation with distance ? to the crack front becomes:
Pi =
a) ( 1 / v / N L ) + ~ + a ~ ( X / ~ )
(21)
where a~" are polynomial functions of the nodal values/3i of the tractions representation, and the element coordinate ~ ~. For instance, along the line between nodes 1 and 7 (Fig. 2) eqn (21) gives the variation of pi with the distance ? to the crack front, with a~ =/3], a 2 = -/37 + 4/3~ - 3/31 and fi{=2/3~-4p~. Therefore, comparing eqn (21) with eqn (15) or eqn (19), it is concluded that the nodal values/3] represent, except for a constant, the SIF values at node 1. The same can be said with respect to nodes 2 and 3. The integration over the boundary elements of the fundamental solution tractions or displacements times the shape functions, is done numerically. The same approach as for isotropic materials is followed. This procedure was described by the authors in Ref. 9. The only transformation required is the change of the fundamental solution expressions. The singularities contained in the transversely isotropic materials fundamental solution are of the same type as those in the isotropic fundamental solution and the singularity contained in the near crack tractions representation is also the same. A more involved discussion on the numerical integration of the fundamental solution for cases of general anisotropy can be found in Schclar.~3
6 STRESS I N T E N S I T Y F A C T O R S C O M P U T A T I O N Once the displacement and traction nodal values have been computed by solution of the boundary element system of equations, the SIF can be obtained directly from the traction nodal values at the crack front. Assuming that z is the material axis of symmetry perpendicular to the crack plane and the t-axis is tangent to the crack front line at the point where a node k is located (Fig. 1), the SIF components
293
at this point are K, = / ~ 2 X / ~
(22)
where/3~ is the T (z, n, t) component of the traction nodal value at node k, and L the width of the singular quarterpoint element to which k belongs. The expressions in eqn (22) are obtained by simple comparison of the (l/v/?) term of the representation of p=, p . and p, near the crack front (eqn (21)) and the ( l I v e ) term of the limit for ~--~0 and 0 = 0 of o= (eqn (15)), a,: (eqn (19)) and at: (eqn (19)), respectively. In order to have a second procedure for SIF evaluation, an expression in terms of the crack opening displacement at quarter-point nodes is also obtained. Given a node on the crack front where the SIF is desired, the computed displacement at the quarter-point node next to it, is made equal to the expression of the displacement at the same point in terms of the SIF. To do so, expressions (14) and (18) must be particularized for quarter-point nodes. In case of symmetric loading, u. can be obtained from eqn (14) when 0 = 7r and ? = L / 4 as
u=(L/4, 70 = K I
1 +rnl
l ~m
(23)
and
2[~
131u:(LI4, 7r)
K I ----- , - [ml/(1 + m l ) -
m21(1 + m 2 ) ]
(24)
where nl, n>/31, ml and m2 are given in eqns (16) and (17). In case of skew-symmetric loading, the values of Ktl and Kill can easily be obtained from the expressions of un and u= for 0=Tr and ?=L/4 in eqn (18). Thus,
13,u~(L/4__~)
Km = ~
[C44(m2- mJ)/V/~3] ut(L/4' Tr) [(1 +mz~)/--V/~--~(5~m,)/--~n2]
(26)
where n3 = 2 C 4 J ( C I I - C12). In the case of cracks under a general loading state, the above SIF expressions remain valid with ui (L/4, lr) representing one half of the 'i' (z, n, t) component of the crack opening displacement at the corresponding quarter-point nodes.
7 NUMERICAL EXAMPLES
To validate the technique presented in this paper, several three-dimensional crack problems in transversely isotropic
A. S6ez et al.
294
media have been analyzed. In all cases the boundary is discretized using standard quadratic quadrilateral or triangular elements except for the two rows in contact with the crack front. The elements of the row next to the front and inside the crack (0 = 7r) are quarter-point elements. Those in the row next to the front and outside the crack (0 = 0) are singular-quarter point elements.
7.1 Example 1: Prismatic plate with central crack under uniform traction This problem has been solved in two dimensions tbr isotropic and orthotopic materials. The geometry of the problem is shown in Fig. 3(a). Owing to its symmetry, only one-quarter of the problem has been discretized. Symmetry boundary conditions are applied to the discretized subregion. Plane strain conditions are simulated by prescribing zero normal displacements over the two larger faces of the plate. A plate width t = w/4 is assumed for the present threedimensional model. Two different BE meshes are used to solve the problem, one consists of 100 elements (Fig. 3(b)) and the other of 148 elements. There are eight elements in half of the crack upper or lower surface in the first case and 24 in the second. There are two singular quarter-point elements and two quarter-point elements in the first model and four of each type in the second, i n order to compare with a known two-dimensional solution for isotropic media ]4 the
O"
"";;;"
2h, li"
(a)
,"' ~ ' ~ ' ~ / ""--...L
I
,
./<
(b) Fig. 3. Prismatic plate with central crack: (at geometry; (b) 100 element mesh.
material properties are first assumed as those of an isotropic material (CjI = C33 = X + 2~, C]2 = C]3 = k, C44 = #) with a Poisson's ratio v = 0.3 (material A). Orthotropic materials with the same properties in two of the three orthogonal directions are transversely isotropic. Two of these materials (B and C) for which the problem was solved by Bowie and Freese ~5 are also assumed. Material B has the following properties: Cr] = 12.043 GPa, Cj2 = 4.351 GPa, Cl.~ = 4.918 GPa, C33 = 22.951 GPa and C44 = 4.762 GPa; and material C: Ci] = 11.178 GPa, CF2 = 3.485 GPa, CI3 = 4.399 GPa, C33 = 102.639 GPa and C~4 = 5.88 GPa. The normalized mode-I SIF [FI = Kl/a ~(Tra)] computed using the tractions formula (22) is denoted by Fit and the one computed using the quarter-point displacement formula (24) is denoted by F~u. In both cases the values correspond to nodes in the mid-plane of the plate. Results are shown in Table 1 in comparison with those obtained by Isida 14 for material A and Bowie and Freese 15 for materials B and C, using two-dimensional models. The SIF values obtained using the proposed BE technique are accurate in all cases and show very little dependency on the mesh size or the tbrmula used to evaluate the SIF. Results obtained at other planes parallel to the mid-plane are almost identical to those presented.
7.2 Example 2: Penny shape crack in infinite domain This is a classical problem Ibr which there are analytical solutions both for isotropic and transversely isotropic media. The infinite medium is represented in the present study as a cylinder with a crack in the mid plane. The diameter of the cylinder is assumed to be ten times that of the crack and its length equal to the diameter. The BE mesh used to represent one-eighth of the problem is shown in Fig. 4. It contains 155 elements. The crack surface is represented by the four innermost rows of elements. Six elements are quarter-point and six are singular quarterpoint. The same crack surface discretization will be used in Example 4. Fig. 5 presents a closer view of this part of the discretization. Two different materials are assumed, one is isotropic with a Poisson's ratio v --- 0.3 and the second is a transversely isotropic graphite-epoxy with the following properties: Cll = 13.92MPa, C]2 = 6-92MPa, Cl~ = 6-44 MPa, C3~ = 160.7 MPa and C44 :- 7.07 MPa. Two loading states are assumed: one is a uniform traction which produces a mode-I SIF, and the second is a shear traction in two opposite directions ( _+ x2) over the crack laces which yields mode-II and mode-III SIFs. The SIF analytic solution for both materials when the crack is under mode-I load is Fl=K[/a~(Tra)=2/Tr= 0.637. The values computed using the BE mesh of Fig. 4 with the tractions formula (22) and the displacement formula (24) are: Fit -- 0.638 and Flu = 0.630, respectively, for the isotropic material and Fh = 0.643 and Flu = 0.634 for the graphite-epoxy. All are within 1% of the analytic solution. The analytic solution for mode-II and mode-III SIFs of a penny shape crack in a boundless isotropic domain under
3-D fracture analysis in transversely isotropic solids
295
Table 1. Normalized mode-I SIF computed using the tractions formula (Fit) and the quarter-point displacement formula (Flu). Results obtained by lsida 14 for material A and Bowie and Freese Is for materials B and C (Fir,r) are also shown
Material
A B C
100 elements
148 elements
Flref
Flu
Fit
Flu
Fit
1-319 1.44 1.86
1.336 1.45 1-84
1.327 1.44 1.85
1.336 1.45 1.84
shear loading are:16 F~I = Kulr~;(Tra) = [4sin 0]/[(2 - v)Tr] and F m = Km/r((Tra) = [4(1 - v)cos 0]/[(2 - v)Tr]. In the case of transversely isotropic materials the solution has the same form as for the isotropic case 16 with the only change of substituting 1 - N* instead of the Poisson's ratio v, with N* = [n 3 l/2(m 2 -- ml)]/[(l + ml)(1 + m2)(nl- I/2 _ n2- 1/2)]. For the graphite-epoxy: 1 - N* = 0.431. The values of F~ and Fill computed from the BE analysis using the tractions formula (22) are represented in Fig. 6 versus the position indicated by 0 along the crack front. The 0 angle is measured from the direction perpendicular to the shear traction over the crack surface. The computed results reproduce very accurately the analytical solution for both materials, the isotropic and the transversely isotropic graphite-epoxy.
7.3 Example 3. Elliptical crack in infinite domain
The analytical solution for the Kl SIF of elliptical cracks in infinite isotropic media was obtained by Irwin.17 Expressions for Kll and Km were derived by Kassir and Sih. Is The analytical solution for the SIFs of elliptical cracks in transversely isotropic boundless domains can be found in Ref. 16. In the present section the problem of an elliptical crack in a plate, the dimensions of which are very large compared with the size of the crack, is solved using the BEM approach presented in the previous sections. The same two materials of the previous example are assumed. The
Fig. 4. Model for one-eighth of a cylinder with central penny shape crack (155 elements),
1.334 1.46 1-85
results obtained are compared with the known analytical solutions for the infinite domain. The ellipse aspect ratio is (Fig. 7) b/a = 0.4 and the plate dimensions b/t = 0.2, a/w = 0.1, h/w = 1. The BE discretization of one-eighth of the problem is shown in Fig. 7. It consist of 297 elements. Symmetry boundary conditions are applied to the BE model when the plate is under uniform traction (mode-I) and skew-symmetry conditions when the crack surfaces are subject to a uniform shear traction r along the a axis of the ellipse (modes II and III). The stress intensity factor values are computed from the traction nodal values at the crack front (eqn (22)). Fig. 8 shows results for Kl computed for the two materials and different positions along the crack front. The position is denoted by the angle ~b measured with respect to the largest radius of the ellipse. The computed results are in very good agreement with the analytical solution which in the case of mode-I SIF is the same for both materials. Results for Kn and Knl obtained when the crack is excited by opposite uniform shear tractions over the two faces along the a axis are shown in Fig. 9. The computed values are compared with the analytical solution for the crack in a boundless domain. The agreement between the analytical solutionsl6-18 and the present results is very good both for the isotropic material and for the graphite-epoxy. 7.4 Example 4. Semi-circular surface crack
Consider a semi-circular crack on the surface of a half-space under the effects of a uniform traction in the direction perpendicular to the crack surface. This problem has been
Fig. 5. Boundary discretization for semi-circular crack on the surface of a finite thickness plate (w/a = 10).
296
A. S a e c et al.
1.0 0.8 ~ , 0.6
-5b
a
~ = 0.4
o Graphite-Epoxy
0.0 0.0
.................. 0.2 0.4
0.6
0.8
1.0
~/,
(a) 0.8
--A~=iyfic o Oral~ite-E~ry
0.6
M O2
0.0 0.0 (b)
0.2
0.4
0.6
0.8
1.0
Fig. 7. Boundary discretization for one-eighth of an elliptical crack in a very large plate (297 elements).
20/w
Fig. 6. Normalized SIFs for penny shape crack under uniform shear: (a) mode-II; (b) mode-Ill.
solved for isotropic materials by lsida et al. ~9 The geometry of the problem is represented in the present paper as a semicircular crack in the mid-plane of a semi-cylinder which dimensions are much larger than those of the crack, The same BE mesh (Fig. 4) as for the penny shape crack in infinite domain (Example 2) is used. The only difference is that the mesh now represents one-quarter of the problem and free surface boundary conditions are prescribed over one of the plane faces. Results obtained for different positions along the crack front for an isotropic material with v = 0.3 and for the g r a p h i t e - e p o x y (with the axis of symmetry perpendicular to the crack surface) are given in Fig. 10. The position is denoted by the angle 0 measured from the half-space free surface. It can be seen from the figure that the agreement between the BE and the results of Isida et al. J9 is very good. It is worth noticing that the 1/~ r behavior of the stress is not valid at the free surface (0 --- 0) of the half-space and therefore the SIF concept is meaningless at that point. A second semi-circular surface crack problem is studied. The crack is on the surface of a boundless plate of finite thickness subject to uniform traction. Raju and Newman 2° obtained an approximate numerical solution for the mode-1 stress intensity factor of this problem. A BE model representing a finite dimension square plate is used. To assess the
influence of plate dimensions on Kb several models are analyzed. In all cases the ratio between the crack radius a and the plate thickness t is a/t = 0.4. Plate length and width values from h = w = 10a to h - - w = 2.5a, with the two materials used in the previous examples, are considered. The BE mesh used for w / a = I 0 is shown in Fig. 5. The discretization of the crack zone is the same in all cases. Fig. 11 (a) shows the computed mode-I SIF versus the position along the crack front for the case w/a = 5. The results are compared with the infinite plate analytical solution. BE results obtained with larger dimensions of the plate are very close to those in Fig. 11 (a) and also very close to Raju and Newman's solution 2°, as shown in the figure.
1.0 0.9 0.8 0.7 0.6
o Graphite-Epoxy
0.5 0.4 0.0
'
,
0.2
,
J
,
i
0.4
i
i
0.6
,
~
,
i
0.8
,
1.0
~/~
Fig. 8. Normalized mode-I SIF for elliptical crack in a very large plate under unitbrm traction.
3-D fracture analysis in transversely isotropic solids 0.8
297
1.0
--Analytic x t~ropic 0.6 -
0.9 ~
1
o Ga'aphite-Epoxy
~o.4 -
"~ 0.7
O. R
0.6
0.0
• . . . . . . . . . . . . . 0.2 0.4 0.6
(a)
'
0.8
''
0.4
1.0
~1,
I
,o Graphic-Epoxy
0.0
. . . . . . . 0.2
0.4
(a)
0.6
0.8
1
1.0
2~/w
°
0.9 0.6
o
- - R a j u & N e w m a n (1979 x ~u'ol~c
05 •
R
°
o
o
o
o
X
X
X
X
M
0.8
0.7 . /
0.2
oo 0.0
ot,o ,
0.6 t
o G'raphite..Epox3¢ •
................ 0.2 0.4
(b)
0.8
1.0
2¢/x
When the dimensions of the plate are reduced to values smaller than w/a = 5, the effects of the dimensions of the plate and of the type of material on the SIF become very important. As an example, results for w/a = 2.5 are shown in Fig. 11 (b). The differences between the computed values for a crack in an isotropic small plate, a crack in a transversely isotropic small plate and the general solution of a crack in a large plate, are very important. 1.0 0.9 -'- 0.8 X
0.7
- - h i d a et a]. 0 9 8 4 )
x Isotmpic
o Graphite-Epoxy
0.4 0.0
. . .
i.
0.2
i
,
,
,
0.4
0.4
(b)
Fig. 9. Normalized SIFs for elliptical crack under uniform shear: (a) mode-lI; (b)mode-Ill.
~ " 0.6
i,
0.6
,,
i,
0.8
x
ho~.'©
0.4
0.6
t
,,,,
0.6
. . . . .
0979) 1 /
........ 0.0
0.2
/
o.s
" "
1.0
2#/,
Fig. 11. Normalized mode-I SIF for semi-circular crack under uniform traction on the surface of a plate: (a) w/a = 5; (b) w/a = 2-5. 8 CONCLUSIONS A boundary element formulation using quadratic isoparametric elements, quarter-point elements and singular quarter-point elements for three-dimensional crack problems in transversely isotropic solids has been successfully implemented in the present paper. Stress intensity factors have been evaluated as system unknowns and also as functions of the computed nodal displacements. Solutions to several three-dimensional crack problems in transversely isotropic bodies have been presented. Numerical solutions computed in the present study for rectangular, penny shape and elliptical cracks were shown to be in very good agreement with those in the literature. The proposed approach is very robust and easy to implement. It follows the same scheme as the one developed for isotropic materials by the same authors 9 in a previous paper. The results computed are very accurate even with a relatively small number of elements to represent the crack and its proximity.
,,
1.0
~/* Fig. 10. Normalized mode-I SIF for semi-circular crack under uniform traction on the surface of a half-space.
ACKNOWLEDGEMENTS Financial
support for this work by the Comisi6n
298
A. S~ie= et al.
lnterministerial de Ciencia y Tecnologfa of Spain (PB931191 and PB93-1184) is gratefully acknowledged.
REFERENCES 1. Kassir, M. K. and Sih, G. C. Three-dimensional stresses around elliptical cracks in transversely isotropic solids. Eng. Fract. Mech., 1968, 1, 327-345. 2. Sih, G. C., Paris, P. C. and Irwin, G. R. On cracks in rectilinearly anisotropic bodies. Int. J. Fract. Mech., 1965, I, 189-203. 3. Rajiyah, H. and Atluri, S. N. Analysis of embedded and surface elliptical flaws in transversely isotropic bodies by Finite Element alternating method. J. Appl, Mech., 1991, 58, 435443. 4. Snyder, M. D. and Cruse, T. A. Boundary integral equation analysis of cracked anisotropic plates. Int. J. Fract. , 1975, 11, 315-328. 5. Chan, K. S. and Cruse, T. A. Stress intensity lactors for anisotropic compact-tension specimens with inclined cracks. Eng. Fracl. Mech., 1986, 23, 863-874. 6. Tan, C. L. and Gao, Y. L. Boundary element analysis of plane anisotropic bodies with stress concentrations and cracks. Compos. Struct. , 1992, 20, 17-28. 7. Martinez, J. and Doml'nguez, J. On the use of quarter-point boundary elements for stress intensity factors computation. Int. J. Namer. Meth. Eng., 1984, 20, 1941 - 1950. 8. Ishikawa, H. Application of the boundary element melhod In anisotropic crack problems. In Advances in Boun&uv Element Methods for Fracture Mechanics, ed. M. H. Aliabadi & C. A. Brebbia. Computational Mechanics Publications, Southampton and Elsevier Applied Science, London, 1990, pp. 269-292.
9. Ariza, M. P., Sfiez, A. & Domfnguez, J. A singular element for three-dimensional fracture mechanics analysis. Eng. Anal. Bounda O, Elements, 1997, 20, 275-285. 10. Love, A. E. H. A Treatise on the Mathematical Theoo, qf Elasticity. Dover, New York, 1944. I I. Lekhnitskii, S. G. Theory of Elasticity ql an Anisotropic Elastic Bo&,. Holden-Day, San Francisco, 1963. 12. Pan, Y. C. and Chou, T. W. Point force solution for an inlinite transversely isotropic solid. J. Appl. Mech., ASME, 1976, 43, 608-612. 13. Schclar, N. A. Anisotropic Analysis using Boun&to, Flements. Topics in Engineering, 20, Computational Mechanics Publications, Southampton and Boston, 1994. 14. lsida, M. Arbitrary loading problems of doubly symmetric regions containing a central crack. Eng. Fract. Mech.. 1975, 7, 505-514. 15. Bowie, O. L. and Freese, C. E. Central crack in plane ortholopic rectangular sheet. Int. J. Fract. Mech., 1972, 8, 49-58. 16. Kassir, M. K. & Sih, G. C. Three-dimensional crack problems. In Mechanics q/Fraclure, ed. G. C. Sih. Noordhoff, The Netherlands, 1975. 17. Irwin, G. R. Crack-extension Ii)rce for a part-through crack in a plate. ,/. Appl. Mech.. ASME, 1962, 29, 651-654. 18. Kassir, M. K. and Sih, G. C. Three-dimensional stress distribution around an elliptical crack under arbitrary loadings. J. Appl. Mech., ASME, 1966, 33, 601-61 I. 19. lsida, M., Noguchi, H. and Yoshida, T. Tension and bending of finite thickne~,s plate with a semi-elliptical surface crack, Int. J. Fract. , 1984, 26, 157- 188. 20. Raju, I. S. and Newman, J. C. SIFs for a wide range of sentielliptical surface cracks in finite-thickness plates. Eng. Fract. Mech., 1979, 11,817-829
0955-7997/97/$17.00
ELSEVIER
Three-dimensional fracture analysis in transversely isotropic solids A. Sfiez, M. P. Ariza & J. Dominguez* Escuela Superior de Ingenieros, Universidad de Sevilla, Av. Reina Mercedes s/n, 41012-Sevilla, Spain
(Received 21 August 1997; accepted 5 September 1997)
In this paper a boundary element formulation for three-dimensional crack problems in transversely isotropic bodies is presented. Quarter-point and singular quarter-point elements are implemented in a quadratic isoparametric element context. The point load fundamental solution for transversely isotropic media is implemented. Numerical solutions to several three-dimensional crack problems are obtained. The accuracy and robustness of the present approach for the analysis of fracture mechanics problems in transversely isotropic bodies are shown by comparison of some of the results obtained with existing analytical solutions. The approach is shown to be a simple and useful tool for the evaluation of stress intensity factors in transversely isotropic media. © 1998 Elsevier Science Ltd. All rights reserved Keywords: Cracks, fracture mechanics, stress intensity factor, boundary elements,
transversely isotropic bodies. 1 INTRODUCTION
dimensional anisotropic bodies 3 years previously. More recently, Rajiyah and Atluri 3 generalized the analytical solution for the SIF of a flat elliptical crack, to arbitrary crack face loading. The previous solution obtained by Kassir and Sih I was limited to the case of constant and linear variations of tractions on the crack face. Rajiyah and Atluri 3 employed their generalized solution with the finite element alternating technique to analyze embedded and surface elliptical cracks in transversely isotropic bodies of finite dimensions where the crack plane is perpendicular to the material axis. The number of boundary element papers dedicated to the analysis of fracture mechanics problems in non-isotropic materials is rather small. Among these papers, the early work of Snyder and C r u s e 4 and the more recent one of Chan and Cruse 5 should be mentioned. Cruse and co-workers developed a boundary element approach for the analysis of two-dimensional crack problems in fully anisotropic bodies. They employed a special fundamental solution which incorporates the presence of a traction-free crack in an infinite anisotropic body. The use of such fundamental solution makes it unnecessary to represent the crack surface in the boundary element discretization. The stress intensity factors are computed using a path independent integral once the boundary displacements and tractions have been obtained. In 1992, Tan and Gao 6 presented the quadratic
Transversely isotropic elastic materials are those with an axis of symmetry such that all directions perpendicular to this axis are equivalent. In other words, any plane perpendicular to the axis is a plane of isotropy. Crystals of the hexagonal system are transversely isotropic solids. There is also a significant number of fiber-reinforced composites which show this kind of behavior from a macroscopic point of view. They are composite solids made of unidirectionally oriented fibers, the fiber diameter and the fiber spacing being much smaller than the dimensions of the body. The basic fracture mechanics concepts for transversely isotropic solids were set 30 years ago by Kassir and Sih. I They showed that the stress singularity of the order ~r near the periphery of three-dimensional cracks, well known for isotropic materials, remains in the case of transversely isotropic solids. The magnitude of the local stresses may also be described in this case in terms of stress intensity factors (SIFs). Kassir and Sih obtained expressions of stresses and displacements near the crack front for cracks of arbitrary shape in a plane perpendicular to the material axis of symmetry. Sih et al. 2 had obtained expressions for stress and displacement fields near the tips of cracks in two*To whom correspondence should be addressed. 287
288
A. & i e : et al.
isoparametric formulation for fully anisotropic twodimensional problems. They used quarter-point singular traction crack-tip elements to compute stress intensity factors directly from the traction nodal values at the crack tip as done by Martfnez and Domfnguez7 for isotropic materials. The results obtained in that paper show great accuracy and little mesh dependency, All the previously mentioned boundary element papers deal with two-dimensional problems. To the authors' knowledge, the only three-dimensional boundary element formulation for fracture mechanics in non-isotropic materials was presented by lshikawa. ~ This author studied crack problems in three-dimensional transversely isotropic and fully anisotropic bodies. The approach is based on the subdomain technique and the use of quarter-point elements in combination with transition elements placed next to them. No singularity is introduced to represent the tractions at the crack front. The SIF values are computed from a twodimensional approximation of the crack opening displacements. The expression for the displacement near the crack tip in terms of the SIF in an anisotropic two-dimensional plate is used in planes normal to the crack front. Only results for mode-! SIF are given by Ishikawa. s In the present paper, the nine node singular quarter-point quadratic element for three-dimensional crack problems, as presented by the authors, 9 is extended to transversely isotropic materials. The crack displacement variation is represented by quarter-point elements which are obtained by locating the mid-side nodes of a quadratic element at one quarter of the element width. The singular tractions over the elements next to the crack front are represented by singular functions obtained by dividing the quadratic shape functions by the square root of the distance to the front. This fact constitutes one of the main differences with respect to Ishikawa's approach, s which contains a tractions representation unable to reproduce the actual singular traction distribution. The integration approach over the boundary elements used in the present paper for transversely isotropic media is the same as given in Ref. 9 for isotropic materials since the singularities of the fundamental solution and of the stress field in the proximity of the crack front are of the same type in both cases. Values of the SIF are computed as traction system unknowns at the crack front nodes. They are also computed by making the obtained values of the crack opening displacement at the quarter-point nodes inside the crack, equal to their analytical expression at the same points, which is a known function of the SIFs. The actual threedimensional expressions given by Kassir and Sih ~ are used. As a general rule, the SIF values obtained from crack opening displacements are less accurate and more sensitive to the discretization than those obtained from the traction crack-front nodal values. Nevertheless, in the studies presented in this paper, both approaches yield very accurate results. In any case, having a second approach for stress intensity factor computations provides a very useful crosscheck on the present numerical approach. Numerical examples with different materials and
geometries are presented in this paper. The computed values of the SIF are compared with other exact or approximate results when they exist. The accuracy of the results obtained tot three-dimensional crack problems in transversely isotropic materials is very good.
2 TRANSVERSELY ISOTROPIC MATERIALS
It is useful to review the basic relations governing the behavior of elastic transversely isotropic materials. The strain tensor is defined, as usual, in terms of displacement derivatives 1 8ii = ~(ui,/+ Ui.i)
(1)
The stress-strain relation for a general elastic anisotropic material can be written as G i = Cii~/ek:
(2)
where Ci/k ~ are the elastic constants. There are 81 possible constants in eqn (2) but only 21 are independent (see for instance Ref. m). When a material is transversely isotropic and the axis of symmetry coincides with a Cartesian axis, for instance x~, the stress-strain relation can be written in terms of five independent constants in the following matrix form
= ]
Cll
Cj2
Ct:~
0
0
0
Cj2
Cj]
Cj3
0
0
0
Cj3
Cr3
C~
0
0
0
0
0
0
C44
0
0
0
0
0
0
C44
0
0
0
0
0
0
Cij - Cj2 2
[
811 82 8~
(3) 28 28 28 where Cu are elastic constants in a condensed notation in terms of only two subindices. lsotropic materials can be considered as particular cases of transversely isotropic materials with
Cll =C33=X+2/.t; Ci2=Ci3=)k; C44~#
(4)
The stress-strain relation can be written in the inverse form as eli = aij,C#~kl
(5)
where a6k / are elastic compliances such that the matrix of
289
3-D fracture analysis in transversely isotropic solids elastic compliances is the inverse of the matrix of elastic constants. In the case of transversely isotropic materials, the matrix of elastic compliances may be written in terms of engineering material constants as follows II
'Ell
the origin of coordinates of an infinite transversely isotropic three-dimensional domain was obtained by Pan and Chou 12 in 1976. They represented the displacements in terms of three potential functions. Assuming that the x3-axis is the material axis of symmetry, the displacements at an arbitrary point x are as follows:
E22 3.2 U n i t p o i n t l o a d in the x3 d i r e c t i o n
E33 2E31
2e32
* ~{ 113o~~-
i 1 Eli
- vl2 Ei t
- v31 E33
0
0
0
- vl~
1
- v31
0
0
0
Ell
Ell
E33
-vt3 Etl
-vL3 Etl
1 E33
0
0
0
0
GI3
0
0
0
0
• 17-11
0
0
0
Gl3
1
0
0
0
0
0
0
-G12
O31
• O"12.
1
v~3 Ell
v3] and Gi2 E33
in the I direction, Gu is and vu is the Poisson's Only five of the seven are independent. They
Eli
(8)
Cllni + C44v2Ai'] 1
-1- v~(A i -k- Bi) ~ C44 R2 -q- Ci i x2
(9)
where ut*~stands for the displacement along the k direction due to a unit point load in the l direction, R 2 = x~ + x 2
R~ : R 2 + Zai Zi = ViX3
(6) where E//is the Young's modulus the shear modulus in the/J-plane, ratio relating directions / and J. engineering constants in eqn (6) satisfy the relations
ui(Aiq-Bi)XaZ--~--~) ~ : 1 , 2
022
0"32 0
2{(
"~3 = -- Z
0"33 1
viAixa, RiRi
Ri : Ri + zi va:{
( P I 3 - C 1 3 ) (1PCI 33 3/+C44 C I23 + 24C 4 4 ) }
- ( - l ) c ~ { (P13+CI3)(PI3" 4C3"3 C44-CI3- 2C44)} 1/2 c~= 1,2 (7)
2(1 +vl2) Cll - C12
v3 ~---
2C44
3 BOUNDARY E L E M E N T F O R M U L A T I O N The displacement integral representation and the boundary element discretization for transversely isotropic elastic bodies is the same as for the isotropic case. The only change in the boundary element formulation is that a new fundamental solution, corresponding to a point load in an infinite transversely isotropic domain, must be used. This change obviously transforms the expressions to be integrated over the boundary elements but does not alter the integration approach since the singularities appearing in the transversely isotropic fundamental solution are of the same kind as those in the isotropic fundamental solution.
PI3 : V/Cll C33 and A~ and Bi are given by:
1. PI3 -- C13 - 2C44 ~ 0
CI3+C44 vIA I : - v2A 2 : 47rC33C44(v2 - Vl) Bi = - A i
2. PI3-C]3-2C44=0 A] = A 2 = 0
3.1 F u n d a m e n t a l s o l u t i o n
B I = B 2 -The solution to the problem of a unit point load applied at
CI3+C44 167rCll C44
290
A. S S e : et al.
3.3 U n i t p o i n t load in the x~ direction
i:
the strain components and substi.tution of the stress-strain relation (3). Explicit expressions for the stress tensor components can be found in Pan and Chou's paper. J2 Tractions are immediately obtained by projection over any surface defined by its normal. The singularity contained in the fundamental solution tractions is of the type ( l/r 2) when r tends to zero: i.e., the same as lbr the isotropic case.
Rt /I
+ 2viBi
RiR2.i~2 } + D
110)
{ u~o :- =
-rl.r~ , -rlx~ ] - ui(A" ' - Bi') R~ ' - - 2uiBi RiR;" J
+ D - -X I X2
(II)
R3R*32
1t73 :
I ":i --U i2 li(A i , -B, , ) XR~
CIIA i -- C44u~B i 2.rl
4 T H R E E D I M E N S I O N A L C R A C K S IN TRANSVERSELY ISOTROPIC SOLIDS
According to the work of Kassir and Sih j the stress intensity factors K j (J = I, II, Ill) for cracks in a plane perpendicular to the material axis of symmetry can be defined in terms of the stress distribution in the usual way K I = lim ~ ( o - ) 0 :
}
r~[)
Cli/u[-C4a
CI3"~-C44 D--
1
4"xv:~(744
and A i' and B i' are given by:
I. PI3 Ci3 2C44 ~ 0 -
-
-
-
I ) ~ G 3 -- C44
r~0
KH, = lim ~/(2"n-r)(o,:)0=ll
(12) where
~: KH = lim " ~ ' ) ( ~ , , : ) 0 : , , ; ""
where the stress components refer to a natural local coordinate system with origin at the arbitrary point P on the crack front where the SIF is computed (Fig. 1 ). The n-axis is directed along the normal to the crack front; the t-axis along the tangent to the crack front; and the z-axis perpendicular to the crack plane. The polar coordinates r and 0 are defined in the nz-plane. Displacement and stress components in the proximity of the crack front for a crack in an infinite body loaded symmetrically with respect to the crack plane can be written LIS:
",,-
KI , ~ r / ~ ) ( ]
n/2
isin20)l/2lJ/~ :
A 2 =B~':
":
,h
I + m2
2. PI3 - CI3 - 2C44 = 0 A n' = A 3 _ r ~__
m
B i' =B.'
1
167rCll
[cos0+(cos-0+n;
Jsm-0)l/-] j/2
} +0(r) (14)
uf = 0(r)
1 • 167rC44vi
3.4 U n i t p o i n t load in the x z direction
The displacements are the same as given by eqns (10)-(12) when indices I and 2 are interchanged. The displacements due to a unit point load applied at the origin, given by eqns (8)-(12) are singular as (I/r) when distance r between the load and the observation point tends to zero. This singularity is of the same type as that appearing in the fundamental solution for isotropic media. The stress at any point can be easily derived by differentiation of the above displacement expressions to obtain
" *
ack Boundary
n "~
Fig. 1. Local coordinates at crack boundary point P.
3-D fracture analysis in transversely isotropic solids
291
KI
Kjsin 0 ( ~
G,z = 2(nli2 _ n~12)v/(r~r) { [(COS20 "1- n l -
ID.74 nl'(+ml)
.[cosO+(cos20+nllsin20)ll2]
I sin20)1/2 __ COS 0] 1/2
-I/2 (C0S20 +
n~ 1sin20) - t/2]ln
_ [(cos20 +//2
1sin20)l/2
1172
n~/2(l + m2)
-- COS 0] I/2(COS20 q- n 2 I sin20) - 1/2 } +
0(r o)
a,z = 0(r °)
lcos 0 + (cos20 + n2- t sin20)1/2] - I/2} _4_O(r)
~,,, = O(r°)
and
where n~ and n2 are the two solutions of the quadratic equation of the material properties
KI Unn --
C11C44 n2 + [C13(C13 +
2C44) -
C I I C 3 3 ] n + C33C44 = 0
Cl3ml ~ I----n C I J [cos 0 + (COS20 + n i- I sin 20) 1/2l I/2 (1 + ml)nll/2
(16) (Ci3 + C44)nj, j = l , 2
(cos20 + n ~ I sin20) - 1/2
C33 _ C44nj
mj--
and
CII n2 [cos I sin20) l/2] i/2 (1 -4- m2)n I/2 ~ 0 + (cos20 + n2-
CI 3m2 -
fit = C44(n lt/2 -- n~/2)
For a loading state which is skew-symmetric with respect to the crack plane, the displacement and stress components in the proximity of the crack front are: 1
( COS20 q- nZ Isine0) - I/2 } q_ O(r0) K!
u,,= -
ut t --
Ci3ml -- Ci2nl lc°s0+¢c°s20+n
(17)
-
Kll(n,n.,/,,~sln0(,r/--r~)l/2 ~ - - - ~ v "{[(1 + m , ) n , ] - ' C44(nl 2 _ nil )
[COS 0 q'- (C0S20 q-
l sin20)l/2]t/2
n 1 I sin20)1/2] - I/2 _ [(1 q- m2)n21 - I
[COS 0 + (COS20 "4- n 2
(cos20 + n 1-J sin20)- J/2 _ Cj 3rn2 - Ci2n2 (1 + m2)n I/2
Lsin20) I/2] -
1/2 } -4- O(r)
Kmn~nsin 0 v / ( r - ~ ) /l t
C44 (m2 - m I )
[cos 0 + (cos20 + n f I sin20)l/2] I/2
1 + m2[co s 0 + (cos20 + nl- I sin20) 1/2] - I/2 (cos20 + n2 t sin20) - t/2 } + 0(/i)
o-__--
KI C44
{nll/2[cos 0 + (cos20 + n I- I sin20)l/2] 1/2 (COS20 + n l- t sin20)- t/2 _ n~/2 [COS 0 + (COS 2 0 q- n 2 I sin 20)1/2 ] 1/2 • 2 -I/2 (COS2 0 + n 2- I sin 0) } + 0(r °)
//1
(15)
1 + ml [cos 0 + (cos20 + n f t sin20) In] - In ~ + 0(r) ) n2 (18)
Kll (Vq77~ oQC44(m I -- m2) ml(1n j/~ + m l ) r .Icos O + (cos 20 + n i sin20)ln] i/2
292
A. S a e z et al.
and Kmn~/2(Ci) - Ci
O,tI
KH
~l(cos20+ n f
I sine0)1:2 _ cos 0l i/2
2C44 k/('ar)
2oe2C44(m I -//12) ~(rlr) {(1 + m2)(Cl3m I - - CIIIll
)Hi
I[(cos20
q- I1[
(cos20 + n~ t sin'0) - .2 + O(r o)
I sin:O)I/e
where - - COS
0]1/2(COS20-~11! Isin20)
I/2
I1; -- (1 + m l ) ( C j 3 m 2 -- CI 172)n2 I - cos 0l I/-~
[(cos~.O+nflsin:O),/2
0~ (cos20 + n-, I sin20)
(19)
20Q C4a (m I - 1112) Ci2nl)n
I I
× [(cos20 + hi- I sin:O) I/2 - cos OIIn (cosRO+nl isin20)
I/2
- (1 + ml)(G3m2 - Ci2n2)n2 [(cos20 + n2 1sine0)1/2 _ cos 01 i/2
X (COS20+112 Isin20) I/2} + 0 ( r 0) K H(nj n 2)1/2
--2(n~/2_nln)xfOrr) {[cos20 + n ( I sin20)l/2 - cos 01": lcos20+nl --
sin'0)
i/2
[ ( c o s 2 0 + n 2 I s i•n - t,~,1/~ /) I
"~
•
n
(cos'0 + n2 sln~0) at:=
Kill
1/2}
- - COS
011/2
+ O(r u) - Isin:O)t:-II:-
2 X/(,a-r) ] • )
Ill]
I]
71-1}12) I1 n~/2
-- Ill')
,1~/~
5 DISCRETIZATION
The boundary element representation of three-dimensional crack problems in transversely isotropic bodies is done following the same steps as in fully isotropic media: (I) the domain under consideration is divided into subdomains by cutting it by a section through the crack; (2) each subdomain is analyzed by the BEM; and (3) the subdomains are coupled using equilibrium and compatibility relations. When the crack is in a plane of symmetry or skew-symmetry. the symmetry conditions can be directly applied over one subregion. The elements used are standard nine node or six node quadratic elements except for those with one side at the crack front. Among these elements, the ones inside the crack surface (non-singular traction) are quarter-point elements and those extending inside the material are singular quarter-point elements. These two types of elements have been used by the authors for similar problems in isotropic materials. 9 The quarter-point element (Fig. 2) is a quadrilateral and has its nine nodes on a plane. The displacement representation over the element when quadratic shape functions are used and the mid-side nodes are located at one quarter of the width, can be written as:
,,, =,,,!
,
~__tcosO+(cos20-kn~
(cos20 + n~- sin-0)
(l ~-/ll])(l
i/2 } + 0(r o)
Kn
{(1 q-m2)(Ci3m I -
2(744 (Clt - C ~ : )
i/2 + 0(r o)
t20/
where ? is the distance to the crack front, L the element 7
Kll (nt u2) J/2
6
o , , : = 2(n~/--- n l / 2 ) . ~ r r )
{n ~ I/2[cos 0 + (cos:'0 + n ~- sin 20) i/:, 1I/2 (cose0+n~ -Isin20)
J/:--n~ 1/2
[cos O + (cos20 + n2 Isin20)J/2 i 1/2
xl
3 4
5
I
(cos20 + n £ l sin20)
I/2 } +
O(r o) Fig. 2. Ninc node quarter-point quadratic elemcnt.
3-D fracture analysis in transversely isotropic solids width, and a'i" are polynomial functions of the nodal values of ui and the element coordinate ~l along the crack front. The displacement variation in this quarterpoint element is able to reproduce the displacement behavior in the proximity of the front of a crack in a transversely isotropic material, which according to eqns (14) and (18), have a variation of the type x/~ (notice that the distance to the crack front, r, used in the general equations, has been changed to ? in the BE formulation because r is used for the distance to the collocation point). The singular quarter-point elements have one side along the crack front and are part of the boundary extending inside the material. They are able to represent the (l/x/if) behavior of the stress near the crack front given by eqns (15) and (19) and have the same geometry and displacement representation as the quarter-point elements. In the tractions representation, the standard shape functions are divided by I/X/~/L. Thus, variation with distance ? to the crack front becomes:
Pi =
a) ( 1 / v / N L ) + ~ + a ~ ( X / ~ )
(21)
where a~" are polynomial functions of the nodal values/3i of the tractions representation, and the element coordinate ~ ~. For instance, along the line between nodes 1 and 7 (Fig. 2) eqn (21) gives the variation of pi with the distance ? to the crack front, with a~ =/3], a 2 = -/37 + 4/3~ - 3/31 and fi{=2/3~-4p~. Therefore, comparing eqn (21) with eqn (15) or eqn (19), it is concluded that the nodal values/3] represent, except for a constant, the SIF values at node 1. The same can be said with respect to nodes 2 and 3. The integration over the boundary elements of the fundamental solution tractions or displacements times the shape functions, is done numerically. The same approach as for isotropic materials is followed. This procedure was described by the authors in Ref. 9. The only transformation required is the change of the fundamental solution expressions. The singularities contained in the transversely isotropic materials fundamental solution are of the same type as those in the isotropic fundamental solution and the singularity contained in the near crack tractions representation is also the same. A more involved discussion on the numerical integration of the fundamental solution for cases of general anisotropy can be found in Schclar.~3
6 STRESS I N T E N S I T Y F A C T O R S C O M P U T A T I O N Once the displacement and traction nodal values have been computed by solution of the boundary element system of equations, the SIF can be obtained directly from the traction nodal values at the crack front. Assuming that z is the material axis of symmetry perpendicular to the crack plane and the t-axis is tangent to the crack front line at the point where a node k is located (Fig. 1), the SIF components
293
at this point are K, = / ~ 2 X / ~
(22)
where/3~ is the T (z, n, t) component of the traction nodal value at node k, and L the width of the singular quarterpoint element to which k belongs. The expressions in eqn (22) are obtained by simple comparison of the (l/v/?) term of the representation of p=, p . and p, near the crack front (eqn (21)) and the ( l I v e ) term of the limit for ~--~0 and 0 = 0 of o= (eqn (15)), a,: (eqn (19)) and at: (eqn (19)), respectively. In order to have a second procedure for SIF evaluation, an expression in terms of the crack opening displacement at quarter-point nodes is also obtained. Given a node on the crack front where the SIF is desired, the computed displacement at the quarter-point node next to it, is made equal to the expression of the displacement at the same point in terms of the SIF. To do so, expressions (14) and (18) must be particularized for quarter-point nodes. In case of symmetric loading, u. can be obtained from eqn (14) when 0 = 7r and ? = L / 4 as
u=(L/4, 70 = K I
1 +rnl
l ~m
(23)
and
2[~
131u:(LI4, 7r)
K I ----- , - [ml/(1 + m l ) -
m21(1 + m 2 ) ]
(24)
where nl, n>/31, ml and m2 are given in eqns (16) and (17). In case of skew-symmetric loading, the values of Ktl and Kill can easily be obtained from the expressions of un and u= for 0=Tr and ?=L/4 in eqn (18). Thus,
13,u~(L/4__~)
Km = ~
[C44(m2- mJ)/V/~3] ut(L/4' Tr) [(1 +mz~)/--V/~--~(5~m,)/--~n2]
(26)
where n3 = 2 C 4 J ( C I I - C12). In the case of cracks under a general loading state, the above SIF expressions remain valid with ui (L/4, lr) representing one half of the 'i' (z, n, t) component of the crack opening displacement at the corresponding quarter-point nodes.
7 NUMERICAL EXAMPLES
To validate the technique presented in this paper, several three-dimensional crack problems in transversely isotropic
A. S6ez et al.
294
media have been analyzed. In all cases the boundary is discretized using standard quadratic quadrilateral or triangular elements except for the two rows in contact with the crack front. The elements of the row next to the front and inside the crack (0 = 7r) are quarter-point elements. Those in the row next to the front and outside the crack (0 = 0) are singular-quarter point elements.
7.1 Example 1: Prismatic plate with central crack under uniform traction This problem has been solved in two dimensions tbr isotropic and orthotopic materials. The geometry of the problem is shown in Fig. 3(a). Owing to its symmetry, only one-quarter of the problem has been discretized. Symmetry boundary conditions are applied to the discretized subregion. Plane strain conditions are simulated by prescribing zero normal displacements over the two larger faces of the plate. A plate width t = w/4 is assumed for the present threedimensional model. Two different BE meshes are used to solve the problem, one consists of 100 elements (Fig. 3(b)) and the other of 148 elements. There are eight elements in half of the crack upper or lower surface in the first case and 24 in the second. There are two singular quarter-point elements and two quarter-point elements in the first model and four of each type in the second, i n order to compare with a known two-dimensional solution for isotropic media ]4 the
O"
"";;;"
2h, li"
(a)
,"' ~ ' ~ ' ~ / ""--...L
I
,
./<
(b) Fig. 3. Prismatic plate with central crack: (at geometry; (b) 100 element mesh.
material properties are first assumed as those of an isotropic material (CjI = C33 = X + 2~, C]2 = C]3 = k, C44 = #) with a Poisson's ratio v = 0.3 (material A). Orthotropic materials with the same properties in two of the three orthogonal directions are transversely isotropic. Two of these materials (B and C) for which the problem was solved by Bowie and Freese ~5 are also assumed. Material B has the following properties: Cr] = 12.043 GPa, Cj2 = 4.351 GPa, Cl.~ = 4.918 GPa, C33 = 22.951 GPa and C44 = 4.762 GPa; and material C: Ci] = 11.178 GPa, CF2 = 3.485 GPa, CI3 = 4.399 GPa, C33 = 102.639 GPa and C~4 = 5.88 GPa. The normalized mode-I SIF [FI = Kl/a ~(Tra)] computed using the tractions formula (22) is denoted by Fit and the one computed using the quarter-point displacement formula (24) is denoted by F~u. In both cases the values correspond to nodes in the mid-plane of the plate. Results are shown in Table 1 in comparison with those obtained by Isida 14 for material A and Bowie and Freese 15 for materials B and C, using two-dimensional models. The SIF values obtained using the proposed BE technique are accurate in all cases and show very little dependency on the mesh size or the tbrmula used to evaluate the SIF. Results obtained at other planes parallel to the mid-plane are almost identical to those presented.
7.2 Example 2: Penny shape crack in infinite domain This is a classical problem Ibr which there are analytical solutions both for isotropic and transversely isotropic media. The infinite medium is represented in the present study as a cylinder with a crack in the mid plane. The diameter of the cylinder is assumed to be ten times that of the crack and its length equal to the diameter. The BE mesh used to represent one-eighth of the problem is shown in Fig. 4. It contains 155 elements. The crack surface is represented by the four innermost rows of elements. Six elements are quarter-point and six are singular quarterpoint. The same crack surface discretization will be used in Example 4. Fig. 5 presents a closer view of this part of the discretization. Two different materials are assumed, one is isotropic with a Poisson's ratio v --- 0.3 and the second is a transversely isotropic graphite-epoxy with the following properties: Cll = 13.92MPa, C]2 = 6-92MPa, Cl~ = 6-44 MPa, C3~ = 160.7 MPa and C44 :- 7.07 MPa. Two loading states are assumed: one is a uniform traction which produces a mode-I SIF, and the second is a shear traction in two opposite directions ( _+ x2) over the crack laces which yields mode-II and mode-III SIFs. The SIF analytic solution for both materials when the crack is under mode-I load is Fl=K[/a~(Tra)=2/Tr= 0.637. The values computed using the BE mesh of Fig. 4 with the tractions formula (22) and the displacement formula (24) are: Fit -- 0.638 and Flu = 0.630, respectively, for the isotropic material and Fh = 0.643 and Flu = 0.634 for the graphite-epoxy. All are within 1% of the analytic solution. The analytic solution for mode-II and mode-III SIFs of a penny shape crack in a boundless isotropic domain under
3-D fracture analysis in transversely isotropic solids
295
Table 1. Normalized mode-I SIF computed using the tractions formula (Fit) and the quarter-point displacement formula (Flu). Results obtained by lsida 14 for material A and Bowie and Freese Is for materials B and C (Fir,r) are also shown
Material
A B C
100 elements
148 elements
Flref
Flu
Fit
Flu
Fit
1-319 1.44 1.86
1.336 1.45 1-84
1.327 1.44 1.85
1.336 1.45 1.84
shear loading are:16 F~I = Kulr~;(Tra) = [4sin 0]/[(2 - v)Tr] and F m = Km/r((Tra) = [4(1 - v)cos 0]/[(2 - v)Tr]. In the case of transversely isotropic materials the solution has the same form as for the isotropic case 16 with the only change of substituting 1 - N* instead of the Poisson's ratio v, with N* = [n 3 l/2(m 2 -- ml)]/[(l + ml)(1 + m2)(nl- I/2 _ n2- 1/2)]. For the graphite-epoxy: 1 - N* = 0.431. The values of F~ and Fill computed from the BE analysis using the tractions formula (22) are represented in Fig. 6 versus the position indicated by 0 along the crack front. The 0 angle is measured from the direction perpendicular to the shear traction over the crack surface. The computed results reproduce very accurately the analytical solution for both materials, the isotropic and the transversely isotropic graphite-epoxy.
7.3 Example 3. Elliptical crack in infinite domain
The analytical solution for the Kl SIF of elliptical cracks in infinite isotropic media was obtained by Irwin.17 Expressions for Kll and Km were derived by Kassir and Sih. Is The analytical solution for the SIFs of elliptical cracks in transversely isotropic boundless domains can be found in Ref. 16. In the present section the problem of an elliptical crack in a plate, the dimensions of which are very large compared with the size of the crack, is solved using the BEM approach presented in the previous sections. The same two materials of the previous example are assumed. The
Fig. 4. Model for one-eighth of a cylinder with central penny shape crack (155 elements),
1.334 1.46 1-85
results obtained are compared with the known analytical solutions for the infinite domain. The ellipse aspect ratio is (Fig. 7) b/a = 0.4 and the plate dimensions b/t = 0.2, a/w = 0.1, h/w = 1. The BE discretization of one-eighth of the problem is shown in Fig. 7. It consist of 297 elements. Symmetry boundary conditions are applied to the BE model when the plate is under uniform traction (mode-I) and skew-symmetry conditions when the crack surfaces are subject to a uniform shear traction r along the a axis of the ellipse (modes II and III). The stress intensity factor values are computed from the traction nodal values at the crack front (eqn (22)). Fig. 8 shows results for Kl computed for the two materials and different positions along the crack front. The position is denoted by the angle ~b measured with respect to the largest radius of the ellipse. The computed results are in very good agreement with the analytical solution which in the case of mode-I SIF is the same for both materials. Results for Kn and Knl obtained when the crack is excited by opposite uniform shear tractions over the two faces along the a axis are shown in Fig. 9. The computed values are compared with the analytical solution for the crack in a boundless domain. The agreement between the analytical solutionsl6-18 and the present results is very good both for the isotropic material and for the graphite-epoxy. 7.4 Example 4. Semi-circular surface crack
Consider a semi-circular crack on the surface of a half-space under the effects of a uniform traction in the direction perpendicular to the crack surface. This problem has been
Fig. 5. Boundary discretization for semi-circular crack on the surface of a finite thickness plate (w/a = 10).
296
A. S a e c et al.
1.0 0.8 ~ , 0.6
-5b
a
~ = 0.4
o Graphite-Epoxy
0.0 0.0
.................. 0.2 0.4
0.6
0.8
1.0
~/,
(a) 0.8
--A~=iyfic o Oral~ite-E~ry
0.6
M O2
0.0 0.0 (b)
0.2
0.4
0.6
0.8
1.0
Fig. 7. Boundary discretization for one-eighth of an elliptical crack in a very large plate (297 elements).
20/w
Fig. 6. Normalized SIFs for penny shape crack under uniform shear: (a) mode-II; (b) mode-Ill.
solved for isotropic materials by lsida et al. ~9 The geometry of the problem is represented in the present paper as a semicircular crack in the mid-plane of a semi-cylinder which dimensions are much larger than those of the crack, The same BE mesh (Fig. 4) as for the penny shape crack in infinite domain (Example 2) is used. The only difference is that the mesh now represents one-quarter of the problem and free surface boundary conditions are prescribed over one of the plane faces. Results obtained for different positions along the crack front for an isotropic material with v = 0.3 and for the g r a p h i t e - e p o x y (with the axis of symmetry perpendicular to the crack surface) are given in Fig. 10. The position is denoted by the angle 0 measured from the half-space free surface. It can be seen from the figure that the agreement between the BE and the results of Isida et al. J9 is very good. It is worth noticing that the 1/~ r behavior of the stress is not valid at the free surface (0 --- 0) of the half-space and therefore the SIF concept is meaningless at that point. A second semi-circular surface crack problem is studied. The crack is on the surface of a boundless plate of finite thickness subject to uniform traction. Raju and Newman 2° obtained an approximate numerical solution for the mode-1 stress intensity factor of this problem. A BE model representing a finite dimension square plate is used. To assess the
influence of plate dimensions on Kb several models are analyzed. In all cases the ratio between the crack radius a and the plate thickness t is a/t = 0.4. Plate length and width values from h = w = 10a to h - - w = 2.5a, with the two materials used in the previous examples, are considered. The BE mesh used for w / a = I 0 is shown in Fig. 5. The discretization of the crack zone is the same in all cases. Fig. 11 (a) shows the computed mode-I SIF versus the position along the crack front for the case w/a = 5. The results are compared with the infinite plate analytical solution. BE results obtained with larger dimensions of the plate are very close to those in Fig. 11 (a) and also very close to Raju and Newman's solution 2°, as shown in the figure.
1.0 0.9 0.8 0.7 0.6
o Graphite-Epoxy
0.5 0.4 0.0
'
,
0.2
,
J
,
i
0.4
i
i
0.6
,
~
,
i
0.8
,
1.0
~/~
Fig. 8. Normalized mode-I SIF for elliptical crack in a very large plate under unitbrm traction.
3-D fracture analysis in transversely isotropic solids 0.8
297
1.0
--Analytic x t~ropic 0.6 -
0.9 ~
1
o Ga'aphite-Epoxy
~o.4 -
"~ 0.7
O. R
0.6
0.0
• . . . . . . . . . . . . . 0.2 0.4 0.6
(a)
'
0.8
''
0.4
1.0
~1,
I
,o Graphic-Epoxy
0.0
. . . . . . . 0.2
0.4
(a)
0.6
0.8
1
1.0
2~/w
°
0.9 0.6
o
- - R a j u & N e w m a n (1979 x ~u'ol~c
05 •
R
°
o
o
o
o
X
X
X
X
M
0.8
0.7 . /
0.2
oo 0.0
ot,o ,
0.6 t
o G'raphite..Epox3¢ •
................ 0.2 0.4
(b)
0.8
1.0
2¢/x
When the dimensions of the plate are reduced to values smaller than w/a = 5, the effects of the dimensions of the plate and of the type of material on the SIF become very important. As an example, results for w/a = 2.5 are shown in Fig. 11 (b). The differences between the computed values for a crack in an isotropic small plate, a crack in a transversely isotropic small plate and the general solution of a crack in a large plate, are very important. 1.0 0.9 -'- 0.8 X
0.7
- - h i d a et a]. 0 9 8 4 )
x Isotmpic
o Graphite-Epoxy
0.4 0.0
. . .
i.
0.2
i
,
,
,
0.4
0.4
(b)
Fig. 9. Normalized SIFs for elliptical crack under uniform shear: (a) mode-lI; (b)mode-Ill.
~ " 0.6
i,
0.6
,,
i,
0.8
x
ho~.'©
0.4
0.6
t
,,,,
0.6
. . . . .
0979) 1 /
........ 0.0
0.2
/
o.s
" "
1.0
2#/,
Fig. 11. Normalized mode-I SIF for semi-circular crack under uniform traction on the surface of a plate: (a) w/a = 5; (b) w/a = 2-5. 8 CONCLUSIONS A boundary element formulation using quadratic isoparametric elements, quarter-point elements and singular quarter-point elements for three-dimensional crack problems in transversely isotropic solids has been successfully implemented in the present paper. Stress intensity factors have been evaluated as system unknowns and also as functions of the computed nodal displacements. Solutions to several three-dimensional crack problems in transversely isotropic bodies have been presented. Numerical solutions computed in the present study for rectangular, penny shape and elliptical cracks were shown to be in very good agreement with those in the literature. The proposed approach is very robust and easy to implement. It follows the same scheme as the one developed for isotropic materials by the same authors 9 in a previous paper. The results computed are very accurate even with a relatively small number of elements to represent the crack and its proximity.
,,
1.0
~/* Fig. 10. Normalized mode-I SIF for semi-circular crack under uniform traction on the surface of a half-space.
ACKNOWLEDGEMENTS Financial
support for this work by the Comisi6n
298
A. S~ie= et al.
lnterministerial de Ciencia y Tecnologfa of Spain (PB931191 and PB93-1184) is gratefully acknowledged.
REFERENCES 1. Kassir, M. K. and Sih, G. C. Three-dimensional stresses around elliptical cracks in transversely isotropic solids. Eng. Fract. Mech., 1968, 1, 327-345. 2. Sih, G. C., Paris, P. C. and Irwin, G. R. On cracks in rectilinearly anisotropic bodies. Int. J. Fract. Mech., 1965, I, 189-203. 3. Rajiyah, H. and Atluri, S. N. Analysis of embedded and surface elliptical flaws in transversely isotropic bodies by Finite Element alternating method. J. Appl, Mech., 1991, 58, 435443. 4. Snyder, M. D. and Cruse, T. A. Boundary integral equation analysis of cracked anisotropic plates. Int. J. Fract. , 1975, 11, 315-328. 5. Chan, K. S. and Cruse, T. A. Stress intensity lactors for anisotropic compact-tension specimens with inclined cracks. Eng. Fracl. Mech., 1986, 23, 863-874. 6. Tan, C. L. and Gao, Y. L. Boundary element analysis of plane anisotropic bodies with stress concentrations and cracks. Compos. Struct. , 1992, 20, 17-28. 7. Martinez, J. and Doml'nguez, J. On the use of quarter-point boundary elements for stress intensity factors computation. Int. J. Namer. Meth. Eng., 1984, 20, 1941 - 1950. 8. Ishikawa, H. Application of the boundary element melhod In anisotropic crack problems. In Advances in Boun&uv Element Methods for Fracture Mechanics, ed. M. H. Aliabadi & C. A. Brebbia. Computational Mechanics Publications, Southampton and Elsevier Applied Science, London, 1990, pp. 269-292.
9. Ariza, M. P., Sfiez, A. & Domfnguez, J. A singular element for three-dimensional fracture mechanics analysis. Eng. Anal. Bounda O, Elements, 1997, 20, 275-285. 10. Love, A. E. H. A Treatise on the Mathematical Theoo, qf Elasticity. Dover, New York, 1944. I I. Lekhnitskii, S. G. Theory of Elasticity ql an Anisotropic Elastic Bo&,. Holden-Day, San Francisco, 1963. 12. Pan, Y. C. and Chou, T. W. Point force solution for an inlinite transversely isotropic solid. J. Appl. Mech., ASME, 1976, 43, 608-612. 13. Schclar, N. A. Anisotropic Analysis using Boun&to, Flements. Topics in Engineering, 20, Computational Mechanics Publications, Southampton and Boston, 1994. 14. lsida, M. Arbitrary loading problems of doubly symmetric regions containing a central crack. Eng. Fract. Mech.. 1975, 7, 505-514. 15. Bowie, O. L. and Freese, C. E. Central crack in plane ortholopic rectangular sheet. Int. J. Fract. Mech., 1972, 8, 49-58. 16. Kassir, M. K. & Sih, G. C. Three-dimensional crack problems. In Mechanics q/Fraclure, ed. G. C. Sih. Noordhoff, The Netherlands, 1975. 17. Irwin, G. R. Crack-extension Ii)rce for a part-through crack in a plate. ,/. Appl. Mech.. ASME, 1962, 29, 651-654. 18. Kassir, M. K. and Sih, G. C. Three-dimensional stress distribution around an elliptical crack under arbitrary loadings. J. Appl. Mech., ASME, 1966, 33, 601-61 I. 19. lsida, M., Noguchi, H. and Yoshida, T. Tension and bending of finite thickne~,s plate with a semi-elliptical surface crack, Int. J. Fract. , 1984, 26, 157- 188. 20. Raju, I. S. and Newman, J. C. SIFs for a wide range of sentielliptical surface cracks in finite-thickness plates. Eng. Fract. Mech., 1979, 11,817-829