- Email: [email protected]
Energy release rates for delamination of composite beams
theoretical and
applied fracture mechanics ELSEVIER
Theoretical and AppliedFracture Mechanics 25 (1996) 225-232
Energy release rates for delamination of composite beams F. Fraternali * Department of Civil Engineering, University of Salerno, 84084 Fisciano (SA), Italy
Abstract
This paper deals with the formulation of a mechanical/numerical model for analyzing the delamination effects of layered composite beams. The laminate is modelled through a multiple-beam model and interfacial constitutive laws are obtained by introducing interlaminar bilateral and unilateral springs. Delamination growth is described by employing the classical energy release rate criterion. A path-following procedure with delamination growth control is presented for the numerical analysis of the given model. Numerical results on delamination buckling and growth in compressed beams are given and comparisons with simplified theories are established.
1. I n t r o d u c t i o n
The loss of adhesion between the laminae of a composite laminate is well recognized as a phenomenon which may drastically reduce the mechanical performance of the material. Often, it originates with the material manufacturing process, due to flaws or imperfections in the interlaminar surfaces. The delamination growth is essentially produced by the action of normal and tangential interlaminar stresses. Such stresses can appear either in the linear-elastic field, due to particular loading conditions or a marked curvature of the structure, and in the field of large or moderately large displacements, as it happens, for example, in the presence of local or global buckling phenomena. The technical literature on the subject of delamination initiation and growth is wide and collects studies mainly developed during the last twenty years. The first, pioneeristic contributions analyzed the problem of buckling and growth of near surface or
* E-mail: [email protected].
symmetric delaminations in one-dimensional or two-dimensional structures [1-4]. They modelled the structure behavior through simplified theories (such as the so-called thin film approximation) and described delamination growth by means of the energy release rate criterion. More refined theories have been proposed in recent years. Most of them describe the delaminated structure as an assembly of several one-dimensional or two-dimensional elements. Mainly, the incipient propagation of initial defects is investigated. The delamination growth is described by employing either fracture mechanics [5-8] or interface models [9-14]. The present work deals with the formulation of a one-dimensional mechanical model for the analysis of delamination problems in composite laminated beams. The proposed model generalizes previous ones oriented to stability analysis and interlaminar stress calculation [ 15-19]. It includes geometric nonlinearities due to moderately large rotations and shear strains. The delamination growth is modelled by combining the energy release rate criterion with in-
0167-8442/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved. PH S0167-8442(96)00024-9
F. Fraternali / Theoretical and Applied Fracture Mechanics 25 (1996) 225-232
226
terrace models. It is stated in terms of the generalized stresses and displacements of the proposed model. The work is completed by finite element applications on buckling and growth of delaminations in compressed composite beams.
2. The mechanical model Consider the equilibrium problem of a plane composite beam, which contains a single delamination as shown in Fig. 1. It is useful to model such a structure as an assembly of N beams: B 0), B (2). . . . . B (N), which correspond to packages of one or more plies. Herein, such a model is referred to as multiple-beam model and the elements B (*) as sublaminates. The delamination separates B (') and B ("+ ~) along the entire beam width b. The mechanical behavior of each individual sublaminate can be described by the one-dimensional theory proposed in [15], that here is just briefly sketched. It accounts for moderately large rotations of the cross-sections, moderately large shear strains, small extensional strains and small warping strains. According to such assumptions and disregarding warping effects, the displacement field of the generic sublaminate (say the kth one) is the following:
u (*) = ~,(*) + (~b(*) +2~
I~,",
A[ *))
"',
(~)
XJk)
Here: XI *), X~k) and are local coordinates on B (k) (Fig. 1); a~ k), A~k) and A(3k) are the corresponding base vectors; v (k) = v(k)(X~ k)) is the centerline displacement field; ~ ( , ) = q~(k)(Xjk) ) is a field of skew tensors. Hereafter, the delaminated portion of the interface OB(') NOB ~'+ ~) will be denoted by J2 as shown in Fig. 1. Constitutive laws between interlaminar stresses and interfacial displacement jumps can be obtained by introducing continuous distributions of transverse and axial elastic springs over the undelaminated interfaces and continuous distributions of unilateral transverse springs over J2. See Fig. 2. Expressing the stiffnesses K of such spring distributions (interface stiffnesses) in the form K = 1/'q, and letting "r/ (penalty parameter) approach zero, it is possible to reproduce the following constraints between the sublaminate displacement fields: [ u ~*)] = 0 over the undelaminated interfaces,
(2)
[u~ ")] _> 0 over a . Here, by definition, [f(,)] = f ( , + I)_f(k).
2.1. Delamination growth The one-dimensional model here discussed allows for an accurate prediction of the interlaminar stress field [ 17-19]. Nevertheless, in studying delamination growth, it seems to be more reasonable to fit to such a model an energetic criterion than a local stress criterion.
B(~) 0
b
x~ x1~ •
•
I
h
I
h (k)
1 ×2
L
B (k) -
---i-
x(k) 2 Fig. I. Multiple beam model
3 x(k)
(k)
Xl
I
•
~-~
b
I I
F. Fraternali / Theoretical and Applied Fracture Mechanics 25 (1996) 225-232 (1/~
/
=1
/
/I
I
,
~ = l
/t
0
if
[ U2(") 3
< 0 ,
if
[ U2(")3
~ 0 .
/B
227
(n)
,
[
J Fig. 2. I n t e r f a c i a l s p r i n g s .
Upon assuming that the material of the beam is linearly elastic and orthotropic (see [15]), in the present work the classical Griffit fracture criterion is adopted. The association of such a criterion with the interface model defined above is particularly convenient to follow the shape transformation of the delaminated structure. Let the axial coordinate of the delamination tip T,~ ( a = 1, 2) be denoted by X~~) (Fig. 3). Applicating energy balance to delamination growth, there results the following variational equation (see e.g. [20,21]):
617+ Gcb6X~ ~) = 0 (4) where: 6X(3~) is a virtual displacement of the delam-
A powerful tool for the computation of G,~ is the integral [22]: OU i
OXi OX3
where: 7~ is a path that encircles the delamination tip; W is the strain-energy-density (per unit volume); N is the outward normal to %; S is the second Piola-Kirchhoff stress tensor. By using the displacement field in Eq. (1) and the paths shown in Fig. 3, Eq. (6) can be reduced to the form: -
(fracture energy). ment, that is the quantity (5)
(7)
q¢ being the function defined by the following equation: N
1 6H b aX~ ~)
G~
ds
(6)
ination tip; 6H is the corresponding variation of the total potential energy; Gc is a material property One is hence interested in computing the energy release rate associated with a virtual tip displace-
Ouh Ouh )
N-1
~b(X3) = Y'. W(k)+½ ) " bK[u<*)] 2 k=l
k=l
N [T(k )dv~) (I')
0")
I
-- k~-'~ =l ~
I
X3
dX s
I
f)
T,
+ N (k) X3 II
+ W)
__ +p(k) dX 3
u-%
Y3 - - ~ (2"1
X3
I
Fig. 3. P a t h s o f G - i n t e g r a l s .
Here: W (k) is the strain-energy-density (per unit length) of B{k); ~o{k) = q~{k)A~k)is the vector associ-
228
F. F r a t e r n a l i / T h e o r e t i c a l a n d A p p l i e d F r a c t u r e M e c h a n i c s 2 5 ( 1 9 9 6 ) 2 2 5 - 2 3 2
ated with tO(k); T (~), N ~k), M (k) and p(k) are the
generalized stresses of B ~k): a _ b / 2 J _ h ( k ) / 2 S23
X(": f b/2 fh'k'/2 S33dXl"dX " - b/2
(10)
- h( k ) / 2
M(k) = fb/2 fh(k)/2 S33Xlk'dXlk dX "-
ill)
b/2 a- h(k)/2
p(,i fh/ r,,,/: S.( XI") dXl*' =
•" - b / 2 "
k) (12)
- h(k)/2
3. Finite element approximation
(R(U, *, X (l,, X3(2))}
In this section, a finite element approximation of the non-linear problem stated in the previous section is briefly discussed. It extends the numerical procedure proposed in [17] to include delamination effects. Consider a discretization of the beam by four-node isoparametric elements and an approximation with Lagrangian cubic polynomials of the geometry and the displacement field of each element [23]. Such approximations and a standard assembly procedure lead one to write the equilibrium equations of the present model in the following discrete form [17]:
R(U,A)=
(
_' )
Ku(U ) + ~ K P ( U )
where (U o, Ao) is a previously computed equilibrium point and A s a prescribed arc-length (for the solution technique of Eq. (14) as in [24]). Once a critical equilibrium point for delamination growth is detected (G~ _> G c a n d / o r G 2 _> Gc), the procedure here proposed leaves the arc-length method to employ a continuation procedure with delamination growth control. Suppose, for example, to find both G t > G c and G 2 > G c. The present procedure suggests that Eq. (14) be substituted with a new extended system, which admits U, A and the actual coordinates X~ ~), X~2) of the delamination tips as unknown (the latter are allowed to coincide with Gauss points of the discrete model). Such a system is the following:
t
c,(u)-Cc cdu)-c
= 0
(15)
where [u(j,")] is the largest component of the displacement jump [u (")] and /~ a fixed value. The solution technique of Eq. (15) consists of a bordering algorithm similar to those usually employed in stability analysis of discrete systems [25,26].
U-AQ=0
(13) Note that U is the (global) nodal displacement vector, A the load multiplier, Q the nodal force vector, K v the (secant) stiffness matrix of the unconstrained problem (disconnected sublaminates); (1/-1/) Kp the contribution to the stiffness due to interlaminar springs (penalty stiffness matrix), and R the residual vector. Up to delamination growth, the well known arclength method [24] is useful for a step-by-step/iterative solution of Eq. (13). It leads one to solve the
4. Numerical results To show the features of the model here presented, numerical simulations of delamination buckling and growth in compressed composite beams were developed. The two-beam model (TBM) of Fig. 4 was analyzed by setting L = 0.17 m, h = 3 mm, b = 10 mm.
Ih
b(')
B(1) P
extended system:
R*(V,A)=
-as =0 (14)
L. . . . . F--
~__ _-q L
4
Fig. 4. T w o - b e a m model used in numerical calculations.
F. Fraternali / Theoretical and Applied Fracture Mechanics 25 (1996) 225-232
Concerning the aspect ratios h(l)/h and f/L, three models were examined. One (a) considers a very thin delamination: h°)/h = 0.01, ,elL = 0.25, the second (b) a moderately thin delamination: h(l)/h = 0.10, f / L = 0.25 and the third a thick symmetric delamination: h(l)/h = I l L = 0.50. The material properties of unidirectional (0 °) graphite/epoxy [27] were used: Young's moduli: E l = E 2 = 11.5 GPa, E 3 = 141.4 GPa; shear moduli: G I 2 = 3.4 GPa, G31 = G32 = 6.0 GPa; Poisson's ratios: vl2 = 0 . 4 3 , v31 1'32 = 0 . 2 8 ; Critical (mode I) energy release rate: Gic = 130 N / m .
The following will use the notations PD and P~) for the value of the applied load which produces a local buckling of the delaminated zone (delamination buckling load) and the portion of PD acting on the delamination (P~Dl) ----PD hO)/h), respectively. The thin film approximation admits analytic solutions for P~D~) and the energy release rate G in post-buckling equilibrium conditions. They can be expressed as follows [2]: ~2E3bh(|)3
3f2 G=
=
Throughout the numerical applications, a uniform mesh of 40 four-node elements was employed to represent the beam axis. The central-point deflection of B°)(v°)), the applied load (P), the energy release rate (G = G l = G 2) and the penalty parameter (7/, Section 2) were normalized as follows: vO)
/' _
p
Pc'
GL 4
G = E3h5,
-~= ~K v
(16)
Here: / is the delamination length (Fig. 4); Pc is the Eulerian critical load of the undelaminated beam:
~r2E3bh3 Pc =
3L2
(17)
K~ is the largest (translational) element of K v (Section 3). First, delamination growth was ignored and a comparison between the present theory and the socalled thin film approximation (TFA) was established. The latter is a simplified theory, conceived for the analysis of near-surface delaminations. By disregarding the interaction between the delaminated zone and the subgrade, it models the first one as a clamped beam of length / (Fig. 4) loaded by a part p(l) of the global load (due to material homogeneity, the relation p(l) = ph(l)/h holds in the present case).
229
18f4
(18)
'
~
-1
~ +3 I p(Dl)
(19)
Numerical solutions for PD were obtained through the present theory by considering decreasing values of the parameter ~. To this end, the numerical procedure given in [26] was employed (see also [16,17]). The present results corresponding to the above cases a and b (thin and moderately thin delaminations) are compared in Table 1 with the TFA analytic solutions. Fig. 5 shows the pattern of the curves P - ~(~) for each of the three examined cases. They were obtained by employing the TBM and the path-following procedure described in Section 3. The curves corresponding to the models of a and b show a similar shape. Here, delamination buckling occurs at a load level markedly smaller than Pc (see also Table 1). The post-critical path is very stiff and presents a second bifurcation point, where two paths branch off leading to global buckling. They follow the opening (a', b') and the closure (g', b") of the delamination, respectively. In particular, the latter asymptotically approaches the post-critical equilibrium path of the undelaminated beam.
Table 1 Comparison between the present theory (TBM) and the thin film approximation (TFA) in terms of the normalized delamination buckling load P D / P c ( f / L = 0.25)
hO)/ h TBM ~ = 1 0 -I 0.01 0.10
TFA "~=10 -3
~=10 -5
"~=10 -6
0.0015842 0.0015997 0.0016000 0.0016000 0.0016000 0.1405 0.15391 0.15831 0.15848 0.16000
230
F. Fraternali / Theoretical and Applied Fracture Mechanics 25 (1996) 225-232
~ ~___-V
1.00
o"
1.0
'-n
b' I ................. ,,.-r'" T o', b'
o",~b'" ]
0.80
/
s'
0.8
0.6C
G_
........
)/h=0.01
06
0.40
a
0.20
Ic i i I - O: hO)/h=O,01, t/L=0.25 , I '4 . . . . b: hll)/h=O.lO* ~/L=O 25 - - - c: h°)/h=0 50, ~/L=0.50
........ "...... i
O_ 0.4 /
Present Theory
0.2
.....
t/L=0.25
0.00
0.0 -010
-0106
' -0.I02 v°)
0.42
I
0.06
0.10
I O. 1
0.0
0.i2
/
Fig. 5. Load-displacementcurves of beams with thin (a), moderately thin (b) and thick (c) delaminations. The behavior of the model with a symmetric delamination (model c) is different. As a matter of fact, in such a case, the value of PD is close to that of Pc (it was found: PD = 0"89Pc)" Moreover, the post-critical equilibrium path is flat and global buckling does not occur. Often, symmetric delaminations are analyzed by modelling the delaminated parts of the structure as two elements clamped on the boundary [3,11]. Due to the coincidence of the aspect ratios h(D/l and h/L, such a theory leads to the result PD = PC in the case under consideration. The appreciable difference in the value of P-D between the TBM and the above simplified model shows that the displacements on the boundary of the delaminated elements cannot really be known a priori. Results similar to those of Fig. 5 are given in [6] for delaminated plates. The dependence of the energy release rate upon the load in the post-buckling range is shown in Fig. 6, which presents the results of a comparative analysis between the present theory and the thin film approximation. As before, such a comparison considers the models of a and b (the equilibrium paths denoted by a' and b' in Fig. 5 were followed after delamination buckling). The quoted results underline that the above theories give markedly different values of G, unless the delamination is very thin and the load is considerably smaller than the global buckling load Pc. Since the actual value of G c seems to be reached during global buckling (the used material has G1c =
Thin Film Approx;motTon
0.13
GL 4
/
0.14
0.5
(E.3h s)
Fig. 6. Comparison between computed energy release rate and thin film approximation.
3E3hS/L4),
one is led to conclude that the thin film approximation is not adequate to study delamination growth, even when the delamination is very thin. The delamination growth was studied with reference to the models of b and c. By considering a wide range of Glc = GcL4/(E3 hS) values, and letting the delamination length / reach a value equal to 0.9L, several load-displacement curves were obtained. They are shown in Figs. 7 and 8, where fo denotes the initial value of [ . In the model of b (Fig. 7), it can be observed that the nature of delamination growth depends on the value of the parameter Gc- It is indeed first unstable and then stable for G c a~ 1 and first stable and then unstable for G c > 0.2. It is interesting to note that the ultimate load Pu of the delaminated beam is about equal to 0.7 Pc in the latter case. Gc=3.0 080
05
c~ 0.60 ~" 0..
""~'~
0.40
0.20
V
_ ~
0.1 Gc = Gc
C/(%h~)
h0)/h = 0.10 eo/L 0.00
' --0.12
= 0.25 '
I --0.09
~ '
~
I --0.06
'
v O)
'
I ' -0.03
&~][ ~
Ii 0.00
'
' 0.03
/ eo
Fig. 7. Influenceof delamination growth on load-displacementfor
model b.
F. Fraternali / Theoretical and Applied Fracture Mechanics 25 (1996) 225-232 1.0
References Gc
= Oc L*//(E.Ih ~)
/////11
0.8
c~
0.6
rl
0.4 3.0 0.2
231
1.0
0.5 0.3 0.2 0.1 0.01
T
0.0 -0.035
'
I
-0.025
i
I
-0.015
I
-0.005
v (~) / eo
Fig. 8. Influence of delamination growth on load-displacement for model c.
The delamination growth in the model of c is, on the other hand, always unstable (Fig. 8). Here, the decrease in the compressive load-carrying-capacity of the structure is very strong and the value of Pv is close to 0._3Pc, independently of the value of the parameter G c.
5. Conclusions A one-dimensional mechanical model for the analysis of geometrically non-linear delamination problems in composite beams has been presented. A finite element approximation of the proposed model has been formulated, too. Numerical results on delamination buckling and growth in layered composite beams have been given and comparisons with the analytic solutions of simplified theories have been established. The presented results allow one to conclude on the importance of a global modelling of the delaminated structure behavior, even though the delamination is very thin or symmetric. Use of the present theory for the study of the evolution of multiple delaminations in straight and curved beams awaits attention.
Acknowledgements The author gratefully acknowledges the support of the Italian Ministry of University and Scientific Research (MURST 40%) in the present research.
[1] L.M. Kachanov, Failure of Composite Materials due to Delamination, Polym. Mech. 5, 918-922 (1976). [2] H. Chai, C.D. Babcock and W.G. Knauss, One Dimensional Modelling of Failure in Laminated Plates by Delamination Buckling, Int. J. Solids Struct. 17, 1069-1083 (1981). [3] W.J. Bottega and A. Maewal, Delamination Buckling and Growth in Laminates, J. Appl. Mech. 50, 184-189 (1983). [4] A.G. Evans and J.W. Hutchinson, On the Mechanics of Delamination and Spalling in Compressed Films, Int. J. Solids Struct. 20, 455-466 (1984). [5] B. Storakers and B. Anderson, Non-Linear Plate Theory Applied to Delamination in Composites, J. Mech. Phys. Solids 36, 689-718 (1988). [6] B. Cochelin and M. Potier-Ferry, A Numerical Model for Buckling and Growth of Delamination in Composite Laminates, Comput. Methods Appl. Mech. Eng. 89, 361-380 (1991). [7] E.J. Barbero and J.N. Reddy, An Application of the Generalized Plate Theory to Delamination Buckling, Proc. of the Am. Soc. of Compos. (ASC), VPI and SU, Blacksburg (VA). [8] B.V. Sankar, A Finite Element for Modelling Delaminations in Composite Beams, Comput. Struct. 38, 239-246 (1991). [9] L. Ascione and D. Bruno, On the Delamination Problem of Two-Layer Plates, in: Unilateral Problems in Structural Analysis, ed. G. Del Piero and E. Maceri (Springer-Vedag, 1985). [10] A. Grimaldi and J.N. Reddy, On Delamination in Plates: A Unilateral Contact Approach, in: Unilateral Problems in Structural Analysis, ed. G. Del Piero and F. Maceri (Springer-Vedag, 1985). [11] D. Bruno and A. Grimaldi, Delamination Failure of Layered Composite Plates Loaded in Compression, Int. J. Solids Struct. 26, (1990). [12] P.D. Panagiotopoulos and G. Stavroulakis, The Delamination Effect in Laminated Von Karman Plates under Unilateral Boundary Conditions - A Variational-Hemivariational Inequality Approach, J. Elasticity 23, 69-96 (1990). [13] O. Allix and P. Ladev~ze, Interlaminar Interface Modelling for the Prediction of Delamination, Int. J. Compos. Struct. 22, 235-242 (1992). [14] A. Corigliano, Formulation, Identification and Use of Interface Models in the Numerical Analysis of Composite Delamination, Int. J. Solids Struct. 30, 2779-2811 ('1993). [15] L. Ascione and F. Frateruali, A Moderate Rotation Theory of Laminated Composite C~rved Beams, Int. J. Eng. AnaL Des. 1, 161-176 (1994). [16] L. Ascione and F. Fraternali, A Finite Element Analysis of the Stability of Bimodular Composite Curved Beams, Int. J. Eng. Anal. Des. 1,315-334 (1994). [17] F. Frateruali and G. Bilotti, Non-Linear-Elastic Stress Analysis in Curved Composite Beams, to be published in Comput. Struct. (1996). [18] L. Ascione and F. Frateruali, A Penalty Model for the Analysis of Laminated Curved Beams, Comput. Struct. 45, 985-999 (1992).
232
F. Fraternali / Theoretical and Applied Fracture Mechanics 25 (1996) 225-232
[19] F. Fratemali and J.N. Reddy, A Penalty Model lor the Analysis of Laminated Composite Shells, Int. J. Solids Struct. 30, 3337-3355 (I 992). [20] A. Carpinteri, Scienza delle Costruzioni, Vol. I1 (Pitagora, Bologna, 1992) ch. 20. [21] A. Di Tommaso, Fondamenti di Seienza delle Costruzioni, Part II (P~tron, Bologna, 1993) ch. 6. [22] G.C. Sih, in: Dynamic Aspects of Crack Propagation; Inelastic Behavior of Solids, ed. M.J Kannmen, W.F. Adler, A.R. Rosenfield and R.J. Jaffee (McGraw-Hill, 1970) pp. 607-639. [23] J.N. Reddy, An Introduction to the Finite Element Method, 2nd Ed. (McGraw-Hill, New York, 1992). [24] M.A. Crisfield, A Fast Incremental/Iterative Solution Proce-
dure that Handles Snap Through, Comput. Struct. 13, 55-62 (1981). [25] G. Moore and A. Spence, The Calculation of Turning Points of Nonlinear Equations, SIAM J. Numer. Anal. 17, 567-576 (1980). [26] J.C. Simo and P. Wriggers, A General Procedure for the Direct Computation of Turning and Bifurcation Points, Int. J. Numer. Methods Eng. 30, 155-176 (1990). [27] K.Y. Rhee, An Application of the Single Specimen Technique to a Double Cantilever Beam Specimen to Determine Gt, from a Single Test Record, Compos. Struct. 26, 109-113 (1993).
applied fracture mechanics ELSEVIER
Theoretical and AppliedFracture Mechanics 25 (1996) 225-232
Energy release rates for delamination of composite beams F. Fraternali * Department of Civil Engineering, University of Salerno, 84084 Fisciano (SA), Italy
Abstract
This paper deals with the formulation of a mechanical/numerical model for analyzing the delamination effects of layered composite beams. The laminate is modelled through a multiple-beam model and interfacial constitutive laws are obtained by introducing interlaminar bilateral and unilateral springs. Delamination growth is described by employing the classical energy release rate criterion. A path-following procedure with delamination growth control is presented for the numerical analysis of the given model. Numerical results on delamination buckling and growth in compressed beams are given and comparisons with simplified theories are established.
1. I n t r o d u c t i o n
The loss of adhesion between the laminae of a composite laminate is well recognized as a phenomenon which may drastically reduce the mechanical performance of the material. Often, it originates with the material manufacturing process, due to flaws or imperfections in the interlaminar surfaces. The delamination growth is essentially produced by the action of normal and tangential interlaminar stresses. Such stresses can appear either in the linear-elastic field, due to particular loading conditions or a marked curvature of the structure, and in the field of large or moderately large displacements, as it happens, for example, in the presence of local or global buckling phenomena. The technical literature on the subject of delamination initiation and growth is wide and collects studies mainly developed during the last twenty years. The first, pioneeristic contributions analyzed the problem of buckling and growth of near surface or
* E-mail: [email protected].
symmetric delaminations in one-dimensional or two-dimensional structures [1-4]. They modelled the structure behavior through simplified theories (such as the so-called thin film approximation) and described delamination growth by means of the energy release rate criterion. More refined theories have been proposed in recent years. Most of them describe the delaminated structure as an assembly of several one-dimensional or two-dimensional elements. Mainly, the incipient propagation of initial defects is investigated. The delamination growth is described by employing either fracture mechanics [5-8] or interface models [9-14]. The present work deals with the formulation of a one-dimensional mechanical model for the analysis of delamination problems in composite laminated beams. The proposed model generalizes previous ones oriented to stability analysis and interlaminar stress calculation [ 15-19]. It includes geometric nonlinearities due to moderately large rotations and shear strains. The delamination growth is modelled by combining the energy release rate criterion with in-
0167-8442/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved. PH S0167-8442(96)00024-9
F. Fraternali / Theoretical and Applied Fracture Mechanics 25 (1996) 225-232
226
terrace models. It is stated in terms of the generalized stresses and displacements of the proposed model. The work is completed by finite element applications on buckling and growth of delaminations in compressed composite beams.
2. The mechanical model Consider the equilibrium problem of a plane composite beam, which contains a single delamination as shown in Fig. 1. It is useful to model such a structure as an assembly of N beams: B 0), B (2). . . . . B (N), which correspond to packages of one or more plies. Herein, such a model is referred to as multiple-beam model and the elements B (*) as sublaminates. The delamination separates B (') and B ("+ ~) along the entire beam width b. The mechanical behavior of each individual sublaminate can be described by the one-dimensional theory proposed in [15], that here is just briefly sketched. It accounts for moderately large rotations of the cross-sections, moderately large shear strains, small extensional strains and small warping strains. According to such assumptions and disregarding warping effects, the displacement field of the generic sublaminate (say the kth one) is the following:
u (*) = ~,(*) + (~b(*) +2~
I~,",
A[ *))
"',
(~)
XJk)
Here: XI *), X~k) and are local coordinates on B (k) (Fig. 1); a~ k), A~k) and A(3k) are the corresponding base vectors; v (k) = v(k)(X~ k)) is the centerline displacement field; ~ ( , ) = q~(k)(Xjk) ) is a field of skew tensors. Hereafter, the delaminated portion of the interface OB(') NOB ~'+ ~) will be denoted by J2 as shown in Fig. 1. Constitutive laws between interlaminar stresses and interfacial displacement jumps can be obtained by introducing continuous distributions of transverse and axial elastic springs over the undelaminated interfaces and continuous distributions of unilateral transverse springs over J2. See Fig. 2. Expressing the stiffnesses K of such spring distributions (interface stiffnesses) in the form K = 1/'q, and letting "r/ (penalty parameter) approach zero, it is possible to reproduce the following constraints between the sublaminate displacement fields: [ u ~*)] = 0 over the undelaminated interfaces,
(2)
[u~ ")] _> 0 over a . Here, by definition, [f(,)] = f ( , + I)_f(k).
2.1. Delamination growth The one-dimensional model here discussed allows for an accurate prediction of the interlaminar stress field [ 17-19]. Nevertheless, in studying delamination growth, it seems to be more reasonable to fit to such a model an energetic criterion than a local stress criterion.
B(~) 0
b
x~ x1~ •
•
I
h
I
h (k)
1 ×2
L
B (k) -
---i-
x(k) 2 Fig. I. Multiple beam model
3 x(k)
(k)
Xl
I
•
~-~
b
I I
F. Fraternali / Theoretical and Applied Fracture Mechanics 25 (1996) 225-232 (1/~
/
=1
/
/I
I
,
~ = l
/t
0
if
[ U2(") 3
< 0 ,
if
[ U2(")3
~ 0 .
/B
227
(n)
,
[
J Fig. 2. I n t e r f a c i a l s p r i n g s .
Upon assuming that the material of the beam is linearly elastic and orthotropic (see [15]), in the present work the classical Griffit fracture criterion is adopted. The association of such a criterion with the interface model defined above is particularly convenient to follow the shape transformation of the delaminated structure. Let the axial coordinate of the delamination tip T,~ ( a = 1, 2) be denoted by X~~) (Fig. 3). Applicating energy balance to delamination growth, there results the following variational equation (see e.g. [20,21]):
617+ Gcb6X~ ~) = 0 (4) where: 6X(3~) is a virtual displacement of the delam-
A powerful tool for the computation of G,~ is the integral [22]: OU i
OXi OX3
where: 7~ is a path that encircles the delamination tip; W is the strain-energy-density (per unit volume); N is the outward normal to %; S is the second Piola-Kirchhoff stress tensor. By using the displacement field in Eq. (1) and the paths shown in Fig. 3, Eq. (6) can be reduced to the form: -
(fracture energy). ment, that is the quantity (5)
(7)
q¢ being the function defined by the following equation: N
1 6H b aX~ ~)
G~
ds
(6)
ination tip; 6H is the corresponding variation of the total potential energy; Gc is a material property One is hence interested in computing the energy release rate associated with a virtual tip displace-
Ouh Ouh )
N-1
~b(X3) = Y'. W(k)+½ ) " bK[u<*)] 2 k=l
k=l
N [T(k )dv~) (I')
0")
I
-- k~-'~ =l ~
I
X3
dX s
I
f)
T,
+ N (k) X3 II
+ W)
__ +p(k) dX 3
u-%
Y3 - - ~ (2"1
X3
I
Fig. 3. P a t h s o f G - i n t e g r a l s .
Here: W (k) is the strain-energy-density (per unit length) of B{k); ~o{k) = q~{k)A~k)is the vector associ-
228
F. F r a t e r n a l i / T h e o r e t i c a l a n d A p p l i e d F r a c t u r e M e c h a n i c s 2 5 ( 1 9 9 6 ) 2 2 5 - 2 3 2
ated with tO(k); T (~), N ~k), M (k) and p(k) are the
generalized stresses of B ~k): a _ b / 2 J _ h ( k ) / 2 S23
X(": f b/2 fh'k'/2 S33dXl"dX " - b/2
(10)
- h( k ) / 2
M(k) = fb/2 fh(k)/2 S33Xlk'dXlk dX "-
ill)
b/2 a- h(k)/2
p(,i fh/ r,,,/: S.( XI") dXl*' =
•" - b / 2 "
k) (12)
- h(k)/2
3. Finite element approximation
(R(U, *, X (l,, X3(2))}
In this section, a finite element approximation of the non-linear problem stated in the previous section is briefly discussed. It extends the numerical procedure proposed in [17] to include delamination effects. Consider a discretization of the beam by four-node isoparametric elements and an approximation with Lagrangian cubic polynomials of the geometry and the displacement field of each element [23]. Such approximations and a standard assembly procedure lead one to write the equilibrium equations of the present model in the following discrete form [17]:
R(U,A)=
(
_' )
Ku(U ) + ~ K P ( U )
where (U o, Ao) is a previously computed equilibrium point and A s a prescribed arc-length (for the solution technique of Eq. (14) as in [24]). Once a critical equilibrium point for delamination growth is detected (G~ _> G c a n d / o r G 2 _> Gc), the procedure here proposed leaves the arc-length method to employ a continuation procedure with delamination growth control. Suppose, for example, to find both G t > G c and G 2 > G c. The present procedure suggests that Eq. (14) be substituted with a new extended system, which admits U, A and the actual coordinates X~ ~), X~2) of the delamination tips as unknown (the latter are allowed to coincide with Gauss points of the discrete model). Such a system is the following:
t
c,(u)-Cc cdu)-c
= 0
(15)
where [u(j,")] is the largest component of the displacement jump [u (")] and /~ a fixed value. The solution technique of Eq. (15) consists of a bordering algorithm similar to those usually employed in stability analysis of discrete systems [25,26].
U-AQ=0
(13) Note that U is the (global) nodal displacement vector, A the load multiplier, Q the nodal force vector, K v the (secant) stiffness matrix of the unconstrained problem (disconnected sublaminates); (1/-1/) Kp the contribution to the stiffness due to interlaminar springs (penalty stiffness matrix), and R the residual vector. Up to delamination growth, the well known arclength method [24] is useful for a step-by-step/iterative solution of Eq. (13). It leads one to solve the
4. Numerical results To show the features of the model here presented, numerical simulations of delamination buckling and growth in compressed composite beams were developed. The two-beam model (TBM) of Fig. 4 was analyzed by setting L = 0.17 m, h = 3 mm, b = 10 mm.
Ih
b(')
B(1) P
extended system:
R*(V,A)=
-as =0 (14)
L. . . . . F--
~__ _-q L
4
Fig. 4. T w o - b e a m model used in numerical calculations.
F. Fraternali / Theoretical and Applied Fracture Mechanics 25 (1996) 225-232
Concerning the aspect ratios h(l)/h and f/L, three models were examined. One (a) considers a very thin delamination: h°)/h = 0.01, ,elL = 0.25, the second (b) a moderately thin delamination: h(l)/h = 0.10, f / L = 0.25 and the third a thick symmetric delamination: h(l)/h = I l L = 0.50. The material properties of unidirectional (0 °) graphite/epoxy [27] were used: Young's moduli: E l = E 2 = 11.5 GPa, E 3 = 141.4 GPa; shear moduli: G I 2 = 3.4 GPa, G31 = G32 = 6.0 GPa; Poisson's ratios: vl2 = 0 . 4 3 , v31 1'32 = 0 . 2 8 ; Critical (mode I) energy release rate: Gic = 130 N / m .
The following will use the notations PD and P~) for the value of the applied load which produces a local buckling of the delaminated zone (delamination buckling load) and the portion of PD acting on the delamination (P~Dl) ----PD hO)/h), respectively. The thin film approximation admits analytic solutions for P~D~) and the energy release rate G in post-buckling equilibrium conditions. They can be expressed as follows [2]: ~2E3bh(|)3
3f2 G=
=
Throughout the numerical applications, a uniform mesh of 40 four-node elements was employed to represent the beam axis. The central-point deflection of B°)(v°)), the applied load (P), the energy release rate (G = G l = G 2) and the penalty parameter (7/, Section 2) were normalized as follows: vO)
/' _
p
Pc'
GL 4
G = E3h5,
-~= ~K v
(16)
Here: / is the delamination length (Fig. 4); Pc is the Eulerian critical load of the undelaminated beam:
~r2E3bh3 Pc =
3L2
(17)
K~ is the largest (translational) element of K v (Section 3). First, delamination growth was ignored and a comparison between the present theory and the socalled thin film approximation (TFA) was established. The latter is a simplified theory, conceived for the analysis of near-surface delaminations. By disregarding the interaction between the delaminated zone and the subgrade, it models the first one as a clamped beam of length / (Fig. 4) loaded by a part p(l) of the global load (due to material homogeneity, the relation p(l) = ph(l)/h holds in the present case).
229
18f4
(18)
'
~
-1
~ +3 I p(Dl)
(19)
Numerical solutions for PD were obtained through the present theory by considering decreasing values of the parameter ~. To this end, the numerical procedure given in [26] was employed (see also [16,17]). The present results corresponding to the above cases a and b (thin and moderately thin delaminations) are compared in Table 1 with the TFA analytic solutions. Fig. 5 shows the pattern of the curves P - ~(~) for each of the three examined cases. They were obtained by employing the TBM and the path-following procedure described in Section 3. The curves corresponding to the models of a and b show a similar shape. Here, delamination buckling occurs at a load level markedly smaller than Pc (see also Table 1). The post-critical path is very stiff and presents a second bifurcation point, where two paths branch off leading to global buckling. They follow the opening (a', b') and the closure (g', b") of the delamination, respectively. In particular, the latter asymptotically approaches the post-critical equilibrium path of the undelaminated beam.
Table 1 Comparison between the present theory (TBM) and the thin film approximation (TFA) in terms of the normalized delamination buckling load P D / P c ( f / L = 0.25)
hO)/ h TBM ~ = 1 0 -I 0.01 0.10
TFA "~=10 -3
~=10 -5
"~=10 -6
0.0015842 0.0015997 0.0016000 0.0016000 0.0016000 0.1405 0.15391 0.15831 0.15848 0.16000
230
F. Fraternali / Theoretical and Applied Fracture Mechanics 25 (1996) 225-232
~ ~___-V
1.00
o"
1.0
'-n
b' I ................. ,,.-r'" T o', b'
o",~b'" ]
0.80
/
s'
0.8
0.6C
G_
........
)/h=0.01
06
0.40
a
0.20
Ic i i I - O: hO)/h=O,01, t/L=0.25 , I '4 . . . . b: hll)/h=O.lO* ~/L=O 25 - - - c: h°)/h=0 50, ~/L=0.50
........ "...... i
O_ 0.4 /
Present Theory
0.2
.....
t/L=0.25
0.00
0.0 -010
-0106
' -0.I02 v°)
0.42
I
0.06
0.10
I O. 1
0.0
0.i2
/
Fig. 5. Load-displacementcurves of beams with thin (a), moderately thin (b) and thick (c) delaminations. The behavior of the model with a symmetric delamination (model c) is different. As a matter of fact, in such a case, the value of PD is close to that of Pc (it was found: PD = 0"89Pc)" Moreover, the post-critical equilibrium path is flat and global buckling does not occur. Often, symmetric delaminations are analyzed by modelling the delaminated parts of the structure as two elements clamped on the boundary [3,11]. Due to the coincidence of the aspect ratios h(D/l and h/L, such a theory leads to the result PD = PC in the case under consideration. The appreciable difference in the value of P-D between the TBM and the above simplified model shows that the displacements on the boundary of the delaminated elements cannot really be known a priori. Results similar to those of Fig. 5 are given in [6] for delaminated plates. The dependence of the energy release rate upon the load in the post-buckling range is shown in Fig. 6, which presents the results of a comparative analysis between the present theory and the thin film approximation. As before, such a comparison considers the models of a and b (the equilibrium paths denoted by a' and b' in Fig. 5 were followed after delamination buckling). The quoted results underline that the above theories give markedly different values of G, unless the delamination is very thin and the load is considerably smaller than the global buckling load Pc. Since the actual value of G c seems to be reached during global buckling (the used material has G1c =
Thin Film Approx;motTon
0.13
GL 4
/
0.14
0.5
(E.3h s)
Fig. 6. Comparison between computed energy release rate and thin film approximation.
3E3hS/L4),
one is led to conclude that the thin film approximation is not adequate to study delamination growth, even when the delamination is very thin. The delamination growth was studied with reference to the models of b and c. By considering a wide range of Glc = GcL4/(E3 hS) values, and letting the delamination length / reach a value equal to 0.9L, several load-displacement curves were obtained. They are shown in Figs. 7 and 8, where fo denotes the initial value of [ . In the model of b (Fig. 7), it can be observed that the nature of delamination growth depends on the value of the parameter Gc- It is indeed first unstable and then stable for G c a~ 1 and first stable and then unstable for G c > 0.2. It is interesting to note that the ultimate load Pu of the delaminated beam is about equal to 0.7 Pc in the latter case. Gc=3.0 080
05
c~ 0.60 ~" 0..
""~'~
0.40
0.20
V
_ ~
0.1 Gc = Gc
C/(%h~)
h0)/h = 0.10 eo/L 0.00
' --0.12
= 0.25 '
I --0.09
~ '
~
I --0.06
'
v O)
'
I ' -0.03
&~][ ~
Ii 0.00
'
' 0.03
/ eo
Fig. 7. Influenceof delamination growth on load-displacementfor
model b.
F. Fraternali / Theoretical and Applied Fracture Mechanics 25 (1996) 225-232 1.0
References Gc
= Oc L*//(E.Ih ~)
/////11
0.8
c~
0.6
rl
0.4 3.0 0.2
231
1.0
0.5 0.3 0.2 0.1 0.01
T
0.0 -0.035
'
I
-0.025
i
I
-0.015
I
-0.005
v (~) / eo
Fig. 8. Influence of delamination growth on load-displacement for model c.
The delamination growth in the model of c is, on the other hand, always unstable (Fig. 8). Here, the decrease in the compressive load-carrying-capacity of the structure is very strong and the value of Pv is close to 0._3Pc, independently of the value of the parameter G c.
5. Conclusions A one-dimensional mechanical model for the analysis of geometrically non-linear delamination problems in composite beams has been presented. A finite element approximation of the proposed model has been formulated, too. Numerical results on delamination buckling and growth in layered composite beams have been given and comparisons with the analytic solutions of simplified theories have been established. The presented results allow one to conclude on the importance of a global modelling of the delaminated structure behavior, even though the delamination is very thin or symmetric. Use of the present theory for the study of the evolution of multiple delaminations in straight and curved beams awaits attention.
Acknowledgements The author gratefully acknowledges the support of the Italian Ministry of University and Scientific Research (MURST 40%) in the present research.
[1] L.M. Kachanov, Failure of Composite Materials due to Delamination, Polym. Mech. 5, 918-922 (1976). [2] H. Chai, C.D. Babcock and W.G. Knauss, One Dimensional Modelling of Failure in Laminated Plates by Delamination Buckling, Int. J. Solids Struct. 17, 1069-1083 (1981). [3] W.J. Bottega and A. Maewal, Delamination Buckling and Growth in Laminates, J. Appl. Mech. 50, 184-189 (1983). [4] A.G. Evans and J.W. Hutchinson, On the Mechanics of Delamination and Spalling in Compressed Films, Int. J. Solids Struct. 20, 455-466 (1984). [5] B. Storakers and B. Anderson, Non-Linear Plate Theory Applied to Delamination in Composites, J. Mech. Phys. Solids 36, 689-718 (1988). [6] B. Cochelin and M. Potier-Ferry, A Numerical Model for Buckling and Growth of Delamination in Composite Laminates, Comput. Methods Appl. Mech. Eng. 89, 361-380 (1991). [7] E.J. Barbero and J.N. Reddy, An Application of the Generalized Plate Theory to Delamination Buckling, Proc. of the Am. Soc. of Compos. (ASC), VPI and SU, Blacksburg (VA). [8] B.V. Sankar, A Finite Element for Modelling Delaminations in Composite Beams, Comput. Struct. 38, 239-246 (1991). [9] L. Ascione and D. Bruno, On the Delamination Problem of Two-Layer Plates, in: Unilateral Problems in Structural Analysis, ed. G. Del Piero and E. Maceri (Springer-Vedag, 1985). [10] A. Grimaldi and J.N. Reddy, On Delamination in Plates: A Unilateral Contact Approach, in: Unilateral Problems in Structural Analysis, ed. G. Del Piero and F. Maceri (Springer-Vedag, 1985). [11] D. Bruno and A. Grimaldi, Delamination Failure of Layered Composite Plates Loaded in Compression, Int. J. Solids Struct. 26, (1990). [12] P.D. Panagiotopoulos and G. Stavroulakis, The Delamination Effect in Laminated Von Karman Plates under Unilateral Boundary Conditions - A Variational-Hemivariational Inequality Approach, J. Elasticity 23, 69-96 (1990). [13] O. Allix and P. Ladev~ze, Interlaminar Interface Modelling for the Prediction of Delamination, Int. J. Compos. Struct. 22, 235-242 (1992). [14] A. Corigliano, Formulation, Identification and Use of Interface Models in the Numerical Analysis of Composite Delamination, Int. J. Solids Struct. 30, 2779-2811 ('1993). [15] L. Ascione and F. Frateruali, A Moderate Rotation Theory of Laminated Composite C~rved Beams, Int. J. Eng. AnaL Des. 1, 161-176 (1994). [16] L. Ascione and F. Fraternali, A Finite Element Analysis of the Stability of Bimodular Composite Curved Beams, Int. J. Eng. Anal. Des. 1,315-334 (1994). [17] F. Frateruali and G. Bilotti, Non-Linear-Elastic Stress Analysis in Curved Composite Beams, to be published in Comput. Struct. (1996). [18] L. Ascione and F. Frateruali, A Penalty Model for the Analysis of Laminated Curved Beams, Comput. Struct. 45, 985-999 (1992).
232
F. Fraternali / Theoretical and Applied Fracture Mechanics 25 (1996) 225-232
[19] F. Fratemali and J.N. Reddy, A Penalty Model lor the Analysis of Laminated Composite Shells, Int. J. Solids Struct. 30, 3337-3355 (I 992). [20] A. Carpinteri, Scienza delle Costruzioni, Vol. I1 (Pitagora, Bologna, 1992) ch. 20. [21] A. Di Tommaso, Fondamenti di Seienza delle Costruzioni, Part II (P~tron, Bologna, 1993) ch. 6. [22] G.C. Sih, in: Dynamic Aspects of Crack Propagation; Inelastic Behavior of Solids, ed. M.J Kannmen, W.F. Adler, A.R. Rosenfield and R.J. Jaffee (McGraw-Hill, 1970) pp. 607-639. [23] J.N. Reddy, An Introduction to the Finite Element Method, 2nd Ed. (McGraw-Hill, New York, 1992). [24] M.A. Crisfield, A Fast Incremental/Iterative Solution Proce-
dure that Handles Snap Through, Comput. Struct. 13, 55-62 (1981). [25] G. Moore and A. Spence, The Calculation of Turning Points of Nonlinear Equations, SIAM J. Numer. Anal. 17, 567-576 (1980). [26] J.C. Simo and P. Wriggers, A General Procedure for the Direct Computation of Turning and Bifurcation Points, Int. J. Numer. Methods Eng. 30, 155-176 (1990). [27] K.Y. Rhee, An Application of the Single Specimen Technique to a Double Cantilever Beam Specimen to Determine Gt, from a Single Test Record, Compos. Struct. 26, 109-113 (1993).