- Email: [email protected]
A nonlinear cohesive model for mixed-mode delamination of composite laminates
Accepted Manuscript A nonlinear cohesive model for mixed-mode delamination of composite lami‐ nates P.F. Liu, M.M. Islam PII: DOI: Reference:
S0263-8223(13)00269-9 http://dx.doi.org/10.1016/j.compstruct.2013.05.049 COST 5194
To appear in:
Composite Structures
Please cite this article as: Liu, P.F., Islam, M.M., A nonlinear cohesive model for mixed-mode delamination of composite laminates, Composite Structures (2013), doi: http://dx.doi.org/10.1016/j.compstruct.2013.05.049
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
A nonlinear cohesive model for mixed-mode delamination of composite laminates
P.F.Liu a,b,*, M.M.Islam a a. Department of Civil and Environmental Engineering, Princeton University, Princeton NJ 08544, USA b. Institute of Chemical Machinery and Process Equipment, Zhejiang University, Hangzhou 310027, China * Corresponding author. Tel.: +01(609)258 1619 E-mail address: [email protected] and [email protected] (P.F.Liu) _______________________________________________________________________________ Abstract An anisotropic nonlinear cohesive model based on continuum damage mechanics is proposed to predict the mixed-mode delamination initiation and growth of composite laminates numerically. A damaged exponential traction-displacement jump relationship is proposed for the cohesive model and an exponential damage evolution law on the evolvement of displacement jump is proposed to describe the irreversible delamination crack propagation. The implementation of cohesive model uses zero-thickness cohesive element and middle-plane isoparametric interpolation technique. The mesh size effect is solved by regulating the delamination fracture toughness by the cohesive length, and the numerical convergence problem for the snap instability is addressed by viscous regularization. Four representative delamination cases of composite laminates are used to validate the proposed cohesive model: (a) double cantilever beam (Mode-I), (b) end notch flexure (Mode-II), (c) mixed-mode bending fracture, and (d) fixed ratio mixed mode. The effects of the cohesive strengths and mesh sizes on the global load-displacement response for four cases are studied. In addition, the numerical results using proposed cohesive model are also compared with the results obtained by the virtual crack closure technique and other existing
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
cohesive models. The proposed cohesive model is demonstrated to predict the mixed-mode delamination response of composite laminates accurately and effectively. Keywords: Composite laminates, Delamination, Cohesive modelling, Finite element analysis
1 Introduction Carbon fiber reinforced polymer laminated composites have been increasingly used in areas of the aerospace, construction and new energy use due to the high stiffness-to-weight ratio. These laminated composites are generally manufactured by deploying the fiber layups in different orientations for different layers. These stacked angle-ply layers are expected to resist external loads from different directions. However, the laminated properties of composites introduce complicated failure mechanisms. In most cases, the interlaminar delamination of composites caused by weak shear strength of the polymer matrix leads to the degradation of the physical and mechanical performance of composites. Besides, the interlaminar delamination is usually strongly coupled with the intralaminar transverse matrix cracking [1]. The delamination problem of composites is essentially represented by the crack initiation and propagation. The energy dissipation near the crack tip due to crack propagation is the main reason for the catastrophic material failure, and the mechanical response of composites is largely affected by the interlaminar fracture properties. To analyze the interlaminar fracture properties, some advanced fracture theories have been developed. For example, the virtual crack closure technique (VCCT) [2-4] and the cohesive theory, which are increasingly applied to the research on the delamination problems of composites. While the VCCT is very efficient in calculating the energy release rate for crack propagation within the linear elastic fracture mechanics, it cannot predict the
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
crack initiation. In contrast, Dugdale and Barenblatt [5,6] originally proposed the cohesive process zone conception near the crack tip that described discrete fracture as a material separation across the surface. In the cohesive model, the nonlinear fracture properties are controlled by an assumed traction-displacement jump law. One of the major advantages of the cohesive model is that it can predict both the crack initiation and propagation effectively. Besides, it can avoid the complicated stress singularity at the crack tip. In addition, the width of the failure zone in this discrete representation is approximated as zero, and the length scale associated with failure is infinitely smaller than the length scale associated with the structure [7]. Some cohesive models have been developed to predict the delamination crack initiation and propagation of composite laminates. Particularly, the numerical application of cohesive models in the delamination problem of composites was reviewed [8-12]. On one hand, the triangle cohesive models were proposed for the mixed-mode delamination of composites by many researchers [13-22]. Specially, a discrete triangle cohesive zone model for the delamination problem of composites was proposed [23]. Although the triangle cohesive models have achieved success in predicting the delamination growth of composites, they cannot represent true nonlinear shape of cohesive zone such as the fiber bridging toughening mechanisms in fiber-reinforced composites. On the other hand, the potential-based cohesive models for the mixed-mode fracture were proposed [24-26]. Yet, the numerical convergence and mesh size effect for the finite element implementation for these cohesive models have not yet been solved efficiently. Currently, the exponential-type cohesive models [24-26] and polynomial cohesive models [27] have been used to represent the nonlinear cohesive properties. The advantage of the exponential-type cohesive models is that the exponential function is continuous for the
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
traction-displacement jump cohesive relationship and their derivatives. It could be the reason for the exponential cohesive model to be the most popular for the mathematical derivation and finite element implementation. In this paper, an original nonlinear cohesive model for the mixed-mode delamination of composite laminates is proposed. This cohesive model starts from the Helmholtz free energy and introduces the internal damage to describe the irreversible damage evolution properties during the delamination of composites. The numerical convergence problem using finite element analysis is solved by introducing viscous effect and the mesh size effect is solved by following the thought of the crack band model [28]. The four typical delamination cases: (a) double cantilever beam (DCB), (b) end notch flexure (ENF), (c) mixed-mode bending fracture (MMB) and (d) fixed ratio mixed mode (FRMM) are used to validate the cohesive model by the global load-displacement response. For the MMB and FRMM cases, the damage evolution locuses between individual failure mode are also discussed. The paper is organized as follows. First, we introduce the boundary value problem with an interface discontinuity for linear and nonlinear deformations, respectively; Second, we establish the nonlinear cohesive models for the single-mode and mixed-mode delamination of composite laminates, and emphasize the mesh size effect and convergence problem for implementing the proposed cohesive model using finite element analysis; Third, we introduce the VCCT for the delamination growth of composite laminates by comparison; Finally, we validate the proposed cohesive model and VCCT for the DCB, ENF, MMB and FRMM cases respectively, and also explore the effects of the cohesive strengths on the global load-displacement response and damage evolution locuses.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
2 Boundary value problem with an interface discontinuity Let us consider a bounded domain
3
occupied by an initially linear elastic body, as
shown in Figure 1. The boundary of the domain with the outward normal n is subdivided into either Neumann or Dirichlet boundary conditions defined on t and u , respectively,
t and u
such that u
internal discontinuity d such that
t =Ø. Further, is divided into and by an
in the reference configuration. The
displacement field u across d is discontinuous but continuous on and . Neglecting the inertia and body force, the strong form is given by the balance of linear momentum and the boundary conditions
0 , u ub ,
on u
n Tb ,
on
t
on
d
in
N t u , where ub : u
Tb : t
3
3
(1)
is the field of prescribed displacement on the Dirichlet boundary and
is the field of prescribed traction on the Neumann boundary. The Cauchy stress
tensor is denoted by
. The Cauchy cohesive traction t which acts on the internal boundary
d depends on the displacement jump u at the interface d . The reference normal vector
N points from to . The virtual work equation for the body with volume V containing the cohesive interface d is given by the following weak form
X u : dV t u dA u Tb dA d
t
(2)
3 Cohesive modeling for mixed-mode delamination of composite laminates The exponential cohesive traction-displacement jump relationship governing the interface
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
failure is used, as shown in Figure 2. The expression for each failure mode is given by
ti etic / ui
c
ui exp ui / ui
c
,
(i 1, 2, 3)
(3)
where i=1, 2 and 3 represent the normal and two tangential directions, respectively. The fracture toughness for each fracture mode is equal to the area under the tractiondisplacement jump curve
0
0
t1d u1 et1c u1
c
G1c
t2 d u2 t3d u3 et2c u2
(4) c
0
et3c u3
c
G2c G3c
(5)
where G1c , G2c and G3c are the fracture toughness for fracture mode-I, II and III, respectively. 3.1 Single-mode delamination model based on continuum damage mechanics Based on Eq.(3), an anisotropic damaged cohesive traction-displacement jump relationship is proposed
ti etic / ui
c
u 1 ui exp i c , (i 1, 2, 3) ui 1 di
(6)
where d is the damage vector. The Helmholtz free energy is written as (t , u , d , T ) ( , T ), ( ) e
d
d
1 2
a
(7)
where and are the free energy representing the elastic deformation and damage evolvement e
d
of interface, respectively. T is the temperature and is a vector of internal variable. a is a material parameter. According to the thermodynamic theory, the relationships between the thermodynamic conjugate forces
t
Y , B and internal variables d , are given by
, Y , B a u d
(8)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
The Clausius-Duhem dissipation inequality imposed by the continuum thermodynamics on the evolution laws is [29] 1
t u Y d B T q 0
(9)
T
where is the power of dissipation due to damage. q is the heat flux. If the dissipation reaches the maximum, the damage evolution law is given by
F d F d 0 d d , 0 d Y Y B B
(10)
where F d (Y , d , B) is the damage potential function and is the damage consistency parameter. d
They are assumed to obey the Kuhn-Tucker loading/uploading conditions [29]
d 0, F d 0, d F d 0
(11)
Here, a damage potential function F d (Y , d , B) is proposed for each failure mode if the temperature effect is neglected
F d (Y , d , B) Y Yc B Bc
(12)
From Eqs.(11) and (12), the damage criterion is not satisfied and d = 0 holds at F 0 d
while the further damage appears and d > 0 holds at F 0 . Thus, the consistency condition d
F d 0 for the damage evolution by using Eq.(12) derives the consistency parameter d .
F d . F d B Y d Y/ Y B a d
(13)
where is the unit vector. Here, an exponential function for the conjugate force Y is proposed when the displacement jump based failure criterion is assumed
u c 1 exp u f u Y / a 1 exp(1) where
c
u and
u
f
c
(14)
are the initial and critical failure displacement jump, respectively. By
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
combining Eqs.(11)-(14), the damage variable d is found to be equal to the right-hand formula of Eq.(14). The progression of damage is possible if the actual displacement jump u is equal to the evolving internal variable and increases for different delamination modes. The damage
u is smaller than or
remains constant if
conditions for the load function g ( u
d > 0, if d = 0, if
et c / u et
c
are given by
u > 0,
Damage evolution
u or
u 0,
No damage evolution
f
0
) u
u and
f
The critical displacement jump u u
u is not increasing. Thus, the damage evolution
c
u et c
c
u
c
u
f
in Eq.(14) is numerically solved by Eq.(16)
d u exp u
u exp u / u
(15)
c
f
/ u
c
G
where is a constant approximating unity. For example,
(16) c
u
f
11.76 u is obtained at c
=0.9999 by combining Eqs.(4) or (5) and (16). 3.2 Mixed-mode delamination model based on continuum damage mechanics The traction-displacement jump relationship for mixed-mode delamination is assumed by
t et c / u
c 2 1
c
1 u u exp c u 1 d
t t t
where t
u
c
u1
2
c 2 2
u2
c 2 3
2
and
u3
2
(17)
u
c
u
c 2
1
u2
u
c 2
c 2
3
.
and denotes the MacAuley bracket.
For mixed-mode delamination, the anisotropic damage variable d in Eq.(6) is replaced by the scalar damage variable d in Eq.(18) u c 1 exp u f u d 1 exp(1)
c
The damage evolution conditions for mixed-mode delamination are given by
(18)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
d > 0, if u and u > 0, d = 0, if u or u 0,
Damage evolution
(19)
No damage evolution
The mixed-mode interlaminar fracture toughness G c is assumed to obey the B-K law [30]
G G G G G shear GT c
c 1
c 2
c 1
(20)
where Gshear G2 G3 and GT G1 G2 G3 . is a constant depending on the mode mixity. f
u
The mixed-mode critical failure displacement jump
can be numerically solved by
combining Eqs.(20) and (21)
u
0
f
et c / u
et u c
u exp u / u
c
c
et
c
u
c
u
d u exp u
c
f
f
/ u
c
G
(21) c
tan
The tangent stiffness Dij for mixed-mode interface fracture is derived from Eqs.(6) and (18) by considering the interpenetration problem between the delaminated surface
Dijtan
ti uj
gij 1 f g1 j gij 1 f g1 j
/ uj 1 uj / uj
uj uj
uj 1 uj
c
1 u j / u j
t cj e
c
/ uj
c c uj 1 d uj 1 1 d d , c 2 u j u j 1 d
tc / uj j u j
0
c
exp 1
uj
c
0 u u
,
1
u
c
u 0 or
u u u u
(22)
where
d d u , uj u u j
d u
1 f
u
c
i 1, 2, 3
u c exp u f u 1 exp(1)
ui u u i 1 g1i ui u ui
c
c
(23)
(24)
(25)
f
f
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
and f is a penalty factor and g ij is the Kronecker dalta.
1 is for damage loading and u
0 is for damage uploading. u
3.3 Numerical considerations for cohesive models using finite element analysis The finite element implementation of cohesive model sees Segurado and LLorca [31]. The finite element analysis of proposed cohesive model is performed using ABAQUS-UEL (User element subroutine). The zero-thickness interface element is used and the middle-plane d is used for numerical interpolation. In the following, the mesh size effect and numerical convergence using finite element analysis are addressed respectively. 3.3.1 Solution for the mesh size effect Yang and Cox [12] pointed out a necessary condition for mesh-independent numerical results is that the length of a cohesive zone should be larger than the size of the solid element on either side of crack. For the triangular cohesive model, Turon et al. [32] explored the size effect by introducing the cohesive characteristic length l cz and adjusting the cohesive strength. Yet, the problem for the nonlinear cohesive model is more complicated due to variable interface stiffness. Here, the thought of crack band model [28] is introduced into the exponential cohesive model to solve the size effect. The objective response of the finite element analysis is ensured by regularizing the fracture toughness G1c in Eq.(4) and G2c in Eq.(5)
g c G c / l cz
(26)
c
where g is the energy dissipation per unit volume. The cohesive length l cz was evaluated that is defined as one over which stress increases up to the maximum interfacial strength some distance ahead of the crack tip [5,6,33-38]. For different fracture modes, the generalized formula is given by
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
l cz
EG c
tc
2
(27)
where E is Young’s modulus and is a constant depending on different cohesive models and failure cases. For practical delamination cases, the cohesive lengths are evaluated by adjusting the values for t c and . Here, two points should be noted: ① the cohesive length for mode-II is larger than that for mode-I since the delaminated fracture toughness for mode-II is generally larger than that for mode-I; ② for the mixed-mode delamination, the cohesive length for mode-I and mode-II is taken respectively. 3.3.2 Solution for the convergence problem The convergence problem arises from the softening effect of cohesive model. Essentially, the negative eigenvalues in the stiffness matrix lead to the ill-conditioned finite element equations and unstability of numerical solution. In order to solve the convergence problem, Camanho et al. [15], and Alfano and Crisfield [17] used the local-control arc-length and line searches, Hamitouche et al. [39] introduced the damping effect into cohesive zone model. Specially, Gao and Bower [40] introduced the viscous effect into the Xu and Needleman’s cohesive model [24]. Recently, Liu et al. [41] gave a review on the convergence problem and found the viscous damping technique shows higher calculation efficiency than the arc-length algorithm since the latter imposes additional constraint equations between the displacement increment and the arc-length radius. Here, the viscous damping method is used for stabilizing the numerical solution.
4 VCCT for predicting the delamination growth of composite laminates Figure 3 shows crack propagation from the crack tip node i to j with the increment crack length Δa. Node i is separated into two nodes i1 and i2 after crack propagation. For node i, the
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
relative displacements in three directions (x, y, z) are Δuix, Δuiy and uiz after propagation and the node forces before propagation are Fix, Fiy and Fiz. The energy release rate due to crack propagation is given by [2,3]
G GI GII GIII 1 a 0 2 S lim
a
0
Fix (a)uix (a)da Fiy (a)uiy (a)da Fiz (a) uiz (a)da 0 0 a
a
(28)
where GI , GII and GIII are energy release rate for the mode-I, mode-II and mode-III, and S is the new area generated due to the crack propagation length Δa. Assume a two-step method is adopted by using the node force before crack propagation and node relative displacement after crack propagation, the G due to crack propagation on the crack closure surface is calculated as [2,3]
G
1
Fi x u i x Fi y u i y F 2S
u
iz
iz
(29)
Often, the two-step method can be approximately substituted by the one-step method if the mesh sizes on the crack closure surface are sufficiently small. The relative displacement after crack propagation can be substituted by the relative displacements between nearest node pairs before crack propagation [42]. For example, the relative displacements in three directions for the node pair i1 and i2 can be approximately substituted by those for the node pair e1 and e2. In this analysis, the one-step method is used for the calculation of energy release rate. The master-slave node technique is used to simulate crack propagation in the finite element analysis. The technique divides the node pair at the same position into two nodes by releasing the coupling freedom degree if the critical energy release rate G is larger than the equivalent one Gc according to the B-K law in Eq.(20). The finite element implementation of the VCCT has been integrated into ABAQUS finite element software as a module. In this research, the delamination analysis of composite laminates using the VCCT adopts this module.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
5 Numerical examples and discussions The DCB, ENF, MMB and FRMM cases are used to testify the cohesive model, in which the DCB case belongs to the mode-I delamination type, the ENF case belongs to the mode-II delamination type, and the MMB and FRMM cases belong to the mixed-mode delamination type. In order to study the mesh size effect on the load-displacement response using the cohesive model and VCCT, four mesh sizes with different element side length along the direction of crack propagation h=0.015mm, 0.03mm, 0.1mm and 0.5mm are applied. The number of bulk elements is 2E5, 1E5, 3E4 and 6E3, respectively. The Newton-Cotes integration is used [43] and two integration points are set at each direction for each cohesive element, and the displacement, traction and damage variable as the status variables that depend on the solution are monitored. The contact pairs between the delaminated surface are initially set to prevent the normal interpenetration. The Newton-Raphson algorithm is used and the time increment is properly controlled. in Eq.(27) is taken 0.5 and 0.1 respectively for the mode-I and mode-II delamination. 5.1 Numerical results for the DCB delamination case Figure 4 shows the DCB case with boundary and loads. Two loads are symmetrically and oppositely applied to the upper lamina and low lamina of specimen. The material and size parameters are listed in Table 1 [23]. Figure 5 shows three mesh models for the DCB case and the number of cohesive interface elements is 4668, 700 and 140 for h=0.015mm, 0.1mm and 0.5mm mesh sizes, respectively. Figure 6 shows the crack length a-crack opening displacement 2w curves for the DCB case for three mesh models at t1c 60MPa . The ratio a/2w changes from about 12.5 to 4.5. The delamination crack propagates quickly at the beginning, but then enters the unstable
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
delamination stage until the end of the specimen. Figure 7 shows the global load-displacement curves by considering the effects of the cohesive strengths and mesh sizes. The load first increases and then decreases with the displacement. It can be seen the curves at low cohesive strengths using cohesive model deviate from the true results, but the curve at t1c 60MPa is in good agreement with those by the VCCT and discrete cohesive model at 5.7MPa cohesive strength by [23]. Yet, the peak load using the VCCT is largest. The error may arise from the evaluation of energy release rate as the failure criterion for the VCCT. Spurious numerical oscillation appears using both the cohesive model and VCCT at the softening stage for coarse mesh size, but coarse and fine mesh models lead to consistent results, showing the proposed cohesive model is insensitive to the mesh size. The main reason for numerical oscillation arises from the sudden release of strain energy in the intralaminar materials which leads to instantaneous failure of the element, and thus spurious solution jump appears in the unstable snap-back delamination process. Since the practical cohesive model must reflect the stress distribution associated with damage mechanisms occurring ahead of the physical crack tip [44], the numerical oscillation shows the coarse mesh size cannot accurately capture the cohesive zone stress distribution. However, the curve becomes smooth at the decreasing stage when the mesh size h is equal to the cohesive length l1cz . Xie and Waas [23] also explored the effect of the bias mesh on the curve and found it can more easily lead to oscillation using coarse and bias mesh sizes. According to the present numerical results, only the element side length h along the direction of crack propagation shows large influence on the numerical results. Besides, the initial crack propagation is shown to appear much earlier than the peak load. This conclusion is also consistent with that obtained by Chen et al. [14]. For the mode-I delamination, no convergence problem occurs without the viscous technique.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
5.2 Numerical results for the ENF delamination case Figure 8 shows the ENF case with boundary and load. Three-point bending with an interlaminar initial crack is generally used to test the mode-II fracture toughness. The contact pair is introduced into the delaminated surface. The material and size parameters are listed in Table 2 [23]. The number of cohesive interface elements is 2334, 700 and 140 for h=0.03mm, 0.1mm and 0.5mm mesh sizes, respectively. Figure 9 shows the global load-displacement curves. By using low cohesive strengths and coarse mesh sizes, the load-displacement curves by the cohesive model differ from the results using the VCCT and the discrete cohesive model at 57MPa cohesive strength by Xie and Waas [23]. The load-displacement curve experiences three stages: ① the load first increases with displacement; ② a large ‘snap’ phenomenon appears indicating the unstable crack propagation; ③ and then the load increases rapidly again. Goyal et al. [25] classified this case into the unstable delamination growth at a0 0.35L for the ENF specimen. It can also be seen the curves are basically consistent using three mesh models. Different from the DCB case, the convergence becomes more difficult using the Newton-Raphson algorithm for the ENF case which requires the viscous regularization and linear search algorithm as the cohesive strength increases to 500MPa. The convergence problem due to the increase of the cohesive strength and interface stiffness is also found by Chen et al. [14] and Turon et al. [19]. In addition, Feih [45] found the restriction on the element size can be improved by choosing a larger number of integration points instead of reducing the element sizes. Yet, this still requires sufficiently fine mesh sizes since the combination of a large number of integration points and coarse mesh sizes may also lead to numerical oscillation. The oscillatory performance of interface elements in terms of the selection of integration points is specially studied by Schellekens and de Borst [43]. They
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
showed Gaussian integration can easily lead to the oscillation of traction when large traction gradients are present over an element. Thus, the Newton-Cotes integration is used in this research. 5.3 Numerical results for the MMB delamination case Figure 10 shows the MMB case with boundary and loads. The experimental equipment of the MMB is initially proposed by Reeder and Crews [46] to predict the mixed-mode I/II delamination fracture toughness. By loading with a lever, a single applied load produces bending loads on the specimen. Different mode I/II ratios can be obtained by adjusting the length e of the lever. The material and size parameters are listed in Table 3 [13,23]. The number of cohesive interface elements is 700 and 140 for h=0.1mm and 0.5mm mesh sizes, respectively. Figure 11 shows the load-displacement curves for the MMB case. Similar to the DCB and ENF cases, the curve shows small numerical oscillation at the softening stage for coarse mesh sizes using the cohesive model and VCCT. As the cohesive strength t increases, the initial stiffness et c u of cohesive c
c
interface increases rapidly and differs for normal and shear failure. The initial normal interface stiffness reaches about 5E4N/mm at t1c 500MPa . As the cohesive strength increases to 500MPa, the curve is close to that obtained by the VCCT. Same to the ENF case, the curve for the MMB case also exhibits three distinct stages, in which the sudden ‘snap’ effect shows unstable delamination properties. The delamination appears just at the beginning of loading, much earlier than the moment at which the peak load reaches. The numerical results by the proposed cohesive model are also compared with the results using the interface delamination model and cohesive model [7,13,23]. In general, these curves show good agreement except small errors at the last increasing stage. Compared with the B-K law used by the present work, the power law is given by [13,22, 23,25,27]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
G1 G2 c c 1 G1 G2
(30)
where 1 and 2 denote the linear and quadratic failure criteria, respectively. The power law in Eq.(30) fails to accurately capture the dependence of the mixed-mode fracture toughness on the mode ratio [15]. In contrast, the difficulty using the B-K law is to determine the coefficient in Eq.(20) which yet depends on the mode-mixity and the phase angle between the tangent and normal displacement jumps [30]. Turon et al. [18] explored the effect of different mode ratios on the load-displacement curves and got the value 2.284 in their cohesive models by applying a least-square fit to the experimental data for AS4/PEEK carbon fiber composites. In this research, the determination for is to calculate Eq.(18). We found
u
f
/ u
c
varies little as changes. For example,
u
u f
f
/ u
c
in
c
/ u changes
from 2.1 to 2.3 as varies from 1.5 to 2.5 at t1c t2c 500MPa . Thus, 2.0 is used for the present work. Figure 12(a) shows the damage variable d-mixed mode displacement jump
u curves for three cohesive strengths and Figure 12(b) shows the damage evolution locuses at t1c t2c 500MPa . As the cohesive strength increases, the ability for the interface to resist the
delamination becomes stronger. The approximate elliptic damage evolution locus according to the present cohesive model is similar to that obtained by Turon et al. [18]. The only difference is that the present cohesive model requires higher cohesive strengths than that by Turon et al. [18] in order to get accurate numerical results since a nonlinear exponential cohesive law instead of a simpler triangle relationship is used. Similar to the proposed cohesive model, the conclusion for the requirement for the cohesive strength to increase to high values such as 500MPa in order to lead to accurate results is also obtained by the potential-based exponential cohesive model [26]. Besides, the numerical convergence is also a large challenge for both the triangle and exponential
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
cohesive models for softening problems. For the exponential cohesive models, high cohesive strength adds the difficulty of convergence. Yet, the convergence problem for the triangle cohesive models lies in the sharp and discontinuous transition at the peak traction. Although the artificial viscous method is an ideal one for improving the convergence, the viscous coefficients should be properly chosen, similar to the analysis by Gao and Bower [40]. 5.4 Numerical results for the FRMM delamination case Figure 13 shows the FRMM case with boundary and load. The material and size parameters are listed in Table 4 [14]. Different from the MMB case, the mode-I and II delamination fracture toughness for the FRMM case is different. The number of cohesive interface elements is 1834, 550 and 110 for h=0.03mm, 0.1mm and 0.5mm mesh sizes, respectively. Figure 14 shows the load-displacement curves for the FRMM case. Figure 15 shows the damage evolution locuses. Same to the DCB, ENF and MMB cases, the delamination starts much earlier than the peak load. Similar to the ENF and MMB cases, three stages appear in the process of delamination, accompanied by the sudden softening effect. It can be seen all curves for the DCB, ENF, MMB and FRMM cases are close to the true results as the cohesive strength increases. This conclusion is c
c
also obtained by Xie and Waas [23]. At t1 t2 500MPa , the curve using proposed cohesive model is basically consistent with the results obtained by the experiments, the VCCT and Chen et al. [14]. Also, they studied the effect of parameter in Eq.(30) and found the results using the quadratic failure criterion 2 are more close to the experimental results. In addition, Chandra et al. [8] gave an overview on the selection of cohesive strength (or displacement jump) for current cohesive models and applied cases. The traction-displacement cohesive properties between bi-material interface can be determined by inverse optimization computation by measuring the
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
crack-tip opening displacements directly [47-48]. Beisdes, such as the infra-red crack opening interferometry measurements [49] and laser treatment technique [50] are also used to determine the cohesive properties and to improve the bond toughness of interface. In terms of the sensitivity analysis, the effects of cohesive parameters including the cohesive fracture energy, cohesive shape and cohesive strengths on the mixed-mode delamination mechanisms using cohesive models are studied. The cohesive models as a phenomenological method is expected to expound the true bonding properties for various bi-material interface by adjusting the cohesive zone parameters.
6 Concluding remarks In recent years, the cohesive theory has been increasingly used in the research on the delamination problem of composite laminates. Various cohesive models were developed to predict the delamination initiation and growth. However, until now no universal and numerically robust cohesive model exists for the efficient and accurate simulation for the interface failure. In this research, an anisotropic exponential cohesive model based on continuum damage mechanics is proposed to predict the mixed-mode delamination response of composites. By performing numerical simulation on the DCB, ENF, MMB and FRMM delamination cases, the global load-displacement curves using proposed cohesive model are shown to be in good agreement with the results obtained by using the VCCT and existing cohesive models after evaluating the cohesive characteristic length accurately. It is demonstrated the nonlinear cohesive model can solve the mesh size effect and convergence problems effectively. Using the proposed cohesive model, even the coarse mesh sizes can lead to consistent global load-displacement response. The proposed nonlinear cohesive model is expected to be applied to more delamination problems.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Acknowledgement (1) The author Dr. P.F. Liu at Zhejiang University as a Visiting Scholar at Princeton University would like to thank Prof. Jean H. Prévost at Princeton University for helpful academic guidance and discussions. (2)This research is supported by Natural Science Funding of Zhejiang Province of China (Number: LY13E050002).
References [1] Ladeveze P, Lubineau G, Marsal D. Towards a bridge between the micro- and the mesomechanics of delamination for laminated composites. Compos Sci Technol 2006; 66:698-712. [2] Rybicki EF, Kanninen MF. A finite element calculation of stress intensity factors by a modified crack closure integral. Eng Fract Mech 1977; 9: 931-938. [3] Krueger R. Virtual crack closure technique: history, approach, and applications. Appl Mech Rev 2004; 57: 109-143. [4] Liu PF, Hou SJ, Chu JK, Hu XY, Zhou CL, YL Liu, et al. Finite element analysis of postbuckling and delamination of composite laminates using virtual crack closure technique. Compos Struct 2011; 93:1549-1560. [5] Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids 1960; 8: 100-104. [6] Barenblatt I. Mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech 1962; 7: 55-129. [7] Mosler J, Scheider I. A thermodynamically and variationally consistent class of damage-type cohesive models. J Mech Phys Solids 2011; 59: 1647-1666. [8] Chandra N, Li H, Shet C, Ghonem H. Some issues in the application of cohesive zone models for metal-ceramic interfaces. Inter J Solids Struct 2002; 39: 2827-2855. [9] Elices M, Guinea GV, Gómez J, Planas J. The cohesive zone model advantages, limitations and challenge. Eng Fract Mech 2002; 69: 137-163.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[10] Cornec A, Scheider I, Schwalbe KH. On the practical application of the cohesive model. Eng Fract Mech 2003; 70 :1963-1987. [11] de Borst R. Numerical aspects of cohesive-zone models. Eng Fract Mech 2003; 70: 1743-1757. [12] Yang QD, Cox BN. Cohesive models for damage evolution in laminated composites. Int J Fract 2005; 133: 107-137. [13] Mi Y, Crisfield MA, Davies GA. Progressive delamination using interface elements. J Compos Mater 1998; 32:1246-1272. [14] Chen J, Crisfield M, Kinloch AJ, Busso EP, Matthews FL, Qiu Y. Predicting progressive delamination of composite material specimens via interface elements. Mech Compos Mater Struct 1999; 6:301-317. [15] Camanho PP, Davila CG, De Moura MF. Numerical simulation of mixed-mode progressive delamination in the composite materials. J Compos Mater 2003; 37:1415-1438. [16] Alfano G, Crisfield MA. Finite element interface models for the delamination analysis of laminated composites: mechanics and computational issues. Int J Numer Meth Eng 2001; 50: 1701-1736. [17] Alfano G, Crisfield MA. Solution strategies for the delamination analysis based on a combination of local-control arc-length and line searches. Int J Numer Meth Eng 2003; 58: 999-1048. [18] Turon A, Camanho PP, Costaa J, D´avila CG. A damage model for the simulation of delamination in advanced composites under variable-mode loading. Mech Mater 2006; 38: 1072-1089. [19] Turon A, Camanho PP, Costa J, Renart J. Accurate simulation of delamination growth under mixed-mode loading using cohesive elements: Definition of interlaminar strengths and elastic stiffness. Compos Struct 2010; 92:1857-1864. [20] Swaminathan S, Pagano NJ, Ghosh S. Analysis of interfacial debonding in three-dimensional composite microstructures. J Eng Mater Technol 2006; 128: 96-106. [21] Bui VQ, Iannucci L, Robinson P, Pinho ST. A coupled mixed-mode delamination model for laminated composites. J Compos Mater 2011; 45: 1717-1729. [22] Jiang WG, Hallett SR, Green BG, Wisnom MR. A concise interface constitutive law for
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
analysis of delamination and splitting in composite materials and its application to scaled notched tensile specimens. Int J Numer Meth Eng 2007; 69:1982-1995. [23] Xie D, Waas AM. Discrete cohesive zone model for mixed-mode fracture using finite element analysis. Eng Fract Mech 2006;73:1783-1796. [24] Xu XP, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids 1994;42:1397-1434. [25] Goyal VK, Johnson ER, Davila CG. Irreversible constitutive law for modeling the delamination process using interfacial surface discontinuities. Compos Struct 2004; 65: 289-305. [26] Park K, Paulino GH, Roesler JR. A unified potential-based cohesive model of mixed-mode fracture. J Mech Phys Solids 2009; 57: 891-908. [27] El-Sayed S, Sridharan S. Predicting and tracking interlaminar crack growth in composites using a cohesive layer model. Compos Part B 2001; 32:545-553. [28] Bažant ZP, Oh B. Crack band theory for fracture of concrete. Mater Struct 1983;16:155-177. [29] Lemaitre J, Chaboche, JL, 1990. Mechanics of Solid Materials. Cambridge University Press, Cambridge. [30] Benzeggagh ML, Kenane M. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Compos Sci Technol 1996; 56: 439-449. [31] Segurado J, LLorca J. A new three-dimensional interface finite element to simulate fracture in composites. Inter J Solids Struct 2004; 41: 2977-2993. [32] Turon A, Da´vila CG, Camanho PP, Costa J. An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Eng Fract Mech 2007; 74: 1665-1682. [33] Irwin GR. Plastic zone near a crack and fracture toughness. Proceedings of the seventh Sagamore Ordnance materials conference, vol. IV. New York: Syracuse University; 1960. p. 63-78. [34] Hillerborg A, Mode´er M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement Concr Res 1976; 6: 773-782.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[35] Rice JR. The mechanics of earthquake rupture. In: Dziewonski AM, Boschhi E, editors. Physics of the Earth’s interior. Proceedings of the international school of physics ‘‘Enrico Fermi’’, Course 78, 1979. Amsterdam: Italian Physical Society/North-Holland; 1980. p. 555-649. [36] Bao G, Suo Z. Remarks on crack-bridging concepts. Appl Mech Rev 1992; 24, 355-366. [37] Cox BN, Marshall. Concepts for bridged cracks in fracture and fatigue. Acta Metallurgica et Materialia 1994; 42: 341-363. [38] Harper PW, Hallett SR. Cohesive zone length in numerical simulations of composite delamination. Eng Fract Mech 2008;75:4774-4792. [39] Hamitouche L, Tarfaoui M, Vautrin A. An interface debonding law subject to viscous regularization for avoiding instability: Application to the delamination problems. Eng Fract Mech 2008; 75: 3084-3100. [40] Gao YF, Bower AF. A simple technique for avoiding convergence problems in finite element simulations of crack nucleation and growth on cohesive interfaces. Model Simul Mater Sci Eng 2004; 12(3): 453-463. [41] Liu PF, Zheng JY. Recent developments on damage modeling and finite element analysis for composite laminates: a review. Mater Des 2010; 31: 3825-3834. [42] Xie D, Biggers SB. Progressive crack growth analysis using interface element based on the virtual crack closure technique. Finite Elem Anal Design 2006; 42: 977-984. [43] Schellekens JCJ, de Borst R. On the Numerical Integration of Interface Elements. Inter J Numer Meth Eng 1993; 36: 43-66. [44] Schellekens JCJ, de Borst R. A nonlinear finite-element approach for the analysis of mode-I free edge delamination in composites. Int J Solids Struct 1993;30:1239-1253. [45] Feih S. Development of a user element in ABAQUS for modelling of cohesive laws in composite structures. Risø National Laboratory Roskilde, Denmark, Risø-R-1501, 2005. [46] Reeder JR, Crews Jr JR. Mixed-Mode Bending Method for Delamination Testing. AIAA J 1990; 28: 1270-1276. [47] Jacobsen TK, Sorensen BF. Mode I intra-laminar crack growth in composites-modelling of R-curves from measured bridging laws. Compos Part A 2001; 32: 1-11. [48] Shen B, Stanciulescu I, Paulino GH. Inverse computation of cohesive fracture properties from
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
displacement fields. Inver Probl Sci Eng 2012; 18: 1103-1128. [49] Gowrishankar S, M HX, Liechti K M, Huang R. A comparison of direct and iterative methods for determining traction-separation relations. Int J Fract 2012;177: 109-128. [50] Alfano M, Lubineau G, Furgiuele F, Paulino GH. On the enhancement of bond toughness for Al-epoxy T-peel joints with laser surface treated surfaces. Int J Fract 2011;171:139-150.
Figure captions: Figure 1 Boundary value problem with an interface discontinuity Figure 2 Exponential traction-displacement cohesive law Figure 3 Schematic representation of the VCCT: (a) before crack propagation and (b) after crack propagation Figure 4 The DCB delamination model with boundary and loads Figure 5 Three mesh models for the DCB case Figure 6 The crack length-crack opening displacement curves for the DCB case for three mesh models at t1c 60MPa Figure 7 The load-displacement curves for the DCB case: (a) effects of the cohesive strengths using coarse mesh size h 33l1cz 0.5mm and (b) effects of the mesh sizes at
t1c 60MPa . The curves using proposed cohesive model are also compared with the results obtained using the VCCT and existing cohesive models Figure 8 The ENF delamination model with boundary and load Figure 9 The load-displacement curves for the ENF case: (a) effects of the cohesive strengths using coarse mesh size h 33l2cz 0.5mm and (b) effects of the mesh sizes at
t2c 500MPa . The curves using proposed cohesive model are also compared with the results obtained using the VCCT and existing cohesive models Figure 10 The MMB delamination model with boundary and loads Figure 11 The load-displacement curves for the MMB case: (a) effects of the cohesive strengths using coarse mesh size h 5l1cz 0.5mm and (b) effects of the mesh sizes at
t1c t2c 500MPa . The curves using proposed cohesive model are also compared with the results obtained using the VCCT and existing cohesive models Figure 12 Damage evolution of cohesive elements for the MMB case: (a) damage variable
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
d-mixed mode displacement jump u curve for three cohesive strengths and (b) damage evolution locus at t1c t2c 500MPa . Figure 13 The FRMM delamination model with boundary and load Figure 14
The load-displacement curves for the FRMM case: (a) effects of the cohesive strengths using coarse mesh h 12.5l1cz 0.5mm and (b) effects of the mesh sizes at
t1c t2c 500MPa . The curves using proposed cohesive model are also compared with the results using the VCCT, existing cohesive model and experiment Figure 15 Damage evolution of cohesive elements for the FRMM case: (a) damage variable d-mixed mode displacement jump u curve for three cohesive strengths and (b) damage evolution locus at t1c t2c 500MPa
Table 1 Material and geometry parameters for the DCB case [23] E2 (GPa) G12 (GPa)
E1 (GPa)
135.9
9
5.2
c
v12
G1 ( N / mm) L(mm) a0 (mm) 2h(mm) B(mm)
0.24
0.28
100
30
3
10
Table 2 Material and geometry parameters for the ENF case [23] E1 (GPa)
70
E2 (GPa)
70
c
v12
G2 ( N / mm) L(mm) a0 (mm) 2h(mm) B(mm)
0.24
1.45
100
30
3
10
Table 3 Material and geometry parameters for the MMB case[13,23] c
c
E1 (GPa) E2 (GPa) G12 (GPa) v12 G1 G2 ( N / mm) L(mm) e(mm) a0 (mm) 2h(mm) B(mm)
135.9
9
5.2
0.24
4
100
50
30
3
10
Table 4 Material and geometry parameters for the FRMM case [14] c
c
E1 (GPa) E2 (GPa) G12 (GPa) v12 G1 ( N / mm) G2 ( N / mm) L(mm) a0 (mm) 2h(mm) B(mm)
176
8
6
0.27
0.257
0.856
105
45
3.1
24
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
S0263-8223(13)00269-9 http://dx.doi.org/10.1016/j.compstruct.2013.05.049 COST 5194
To appear in:
Composite Structures
Please cite this article as: Liu, P.F., Islam, M.M., A nonlinear cohesive model for mixed-mode delamination of composite laminates, Composite Structures (2013), doi: http://dx.doi.org/10.1016/j.compstruct.2013.05.049
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
A nonlinear cohesive model for mixed-mode delamination of composite laminates
P.F.Liu a,b,*, M.M.Islam a a. Department of Civil and Environmental Engineering, Princeton University, Princeton NJ 08544, USA b. Institute of Chemical Machinery and Process Equipment, Zhejiang University, Hangzhou 310027, China * Corresponding author. Tel.: +01(609)258 1619 E-mail address: [email protected] and [email protected] (P.F.Liu) _______________________________________________________________________________ Abstract An anisotropic nonlinear cohesive model based on continuum damage mechanics is proposed to predict the mixed-mode delamination initiation and growth of composite laminates numerically. A damaged exponential traction-displacement jump relationship is proposed for the cohesive model and an exponential damage evolution law on the evolvement of displacement jump is proposed to describe the irreversible delamination crack propagation. The implementation of cohesive model uses zero-thickness cohesive element and middle-plane isoparametric interpolation technique. The mesh size effect is solved by regulating the delamination fracture toughness by the cohesive length, and the numerical convergence problem for the snap instability is addressed by viscous regularization. Four representative delamination cases of composite laminates are used to validate the proposed cohesive model: (a) double cantilever beam (Mode-I), (b) end notch flexure (Mode-II), (c) mixed-mode bending fracture, and (d) fixed ratio mixed mode. The effects of the cohesive strengths and mesh sizes on the global load-displacement response for four cases are studied. In addition, the numerical results using proposed cohesive model are also compared with the results obtained by the virtual crack closure technique and other existing
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
cohesive models. The proposed cohesive model is demonstrated to predict the mixed-mode delamination response of composite laminates accurately and effectively. Keywords: Composite laminates, Delamination, Cohesive modelling, Finite element analysis
1 Introduction Carbon fiber reinforced polymer laminated composites have been increasingly used in areas of the aerospace, construction and new energy use due to the high stiffness-to-weight ratio. These laminated composites are generally manufactured by deploying the fiber layups in different orientations for different layers. These stacked angle-ply layers are expected to resist external loads from different directions. However, the laminated properties of composites introduce complicated failure mechanisms. In most cases, the interlaminar delamination of composites caused by weak shear strength of the polymer matrix leads to the degradation of the physical and mechanical performance of composites. Besides, the interlaminar delamination is usually strongly coupled with the intralaminar transverse matrix cracking [1]. The delamination problem of composites is essentially represented by the crack initiation and propagation. The energy dissipation near the crack tip due to crack propagation is the main reason for the catastrophic material failure, and the mechanical response of composites is largely affected by the interlaminar fracture properties. To analyze the interlaminar fracture properties, some advanced fracture theories have been developed. For example, the virtual crack closure technique (VCCT) [2-4] and the cohesive theory, which are increasingly applied to the research on the delamination problems of composites. While the VCCT is very efficient in calculating the energy release rate for crack propagation within the linear elastic fracture mechanics, it cannot predict the
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
crack initiation. In contrast, Dugdale and Barenblatt [5,6] originally proposed the cohesive process zone conception near the crack tip that described discrete fracture as a material separation across the surface. In the cohesive model, the nonlinear fracture properties are controlled by an assumed traction-displacement jump law. One of the major advantages of the cohesive model is that it can predict both the crack initiation and propagation effectively. Besides, it can avoid the complicated stress singularity at the crack tip. In addition, the width of the failure zone in this discrete representation is approximated as zero, and the length scale associated with failure is infinitely smaller than the length scale associated with the structure [7]. Some cohesive models have been developed to predict the delamination crack initiation and propagation of composite laminates. Particularly, the numerical application of cohesive models in the delamination problem of composites was reviewed [8-12]. On one hand, the triangle cohesive models were proposed for the mixed-mode delamination of composites by many researchers [13-22]. Specially, a discrete triangle cohesive zone model for the delamination problem of composites was proposed [23]. Although the triangle cohesive models have achieved success in predicting the delamination growth of composites, they cannot represent true nonlinear shape of cohesive zone such as the fiber bridging toughening mechanisms in fiber-reinforced composites. On the other hand, the potential-based cohesive models for the mixed-mode fracture were proposed [24-26]. Yet, the numerical convergence and mesh size effect for the finite element implementation for these cohesive models have not yet been solved efficiently. Currently, the exponential-type cohesive models [24-26] and polynomial cohesive models [27] have been used to represent the nonlinear cohesive properties. The advantage of the exponential-type cohesive models is that the exponential function is continuous for the
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
traction-displacement jump cohesive relationship and their derivatives. It could be the reason for the exponential cohesive model to be the most popular for the mathematical derivation and finite element implementation. In this paper, an original nonlinear cohesive model for the mixed-mode delamination of composite laminates is proposed. This cohesive model starts from the Helmholtz free energy and introduces the internal damage to describe the irreversible damage evolution properties during the delamination of composites. The numerical convergence problem using finite element analysis is solved by introducing viscous effect and the mesh size effect is solved by following the thought of the crack band model [28]. The four typical delamination cases: (a) double cantilever beam (DCB), (b) end notch flexure (ENF), (c) mixed-mode bending fracture (MMB) and (d) fixed ratio mixed mode (FRMM) are used to validate the cohesive model by the global load-displacement response. For the MMB and FRMM cases, the damage evolution locuses between individual failure mode are also discussed. The paper is organized as follows. First, we introduce the boundary value problem with an interface discontinuity for linear and nonlinear deformations, respectively; Second, we establish the nonlinear cohesive models for the single-mode and mixed-mode delamination of composite laminates, and emphasize the mesh size effect and convergence problem for implementing the proposed cohesive model using finite element analysis; Third, we introduce the VCCT for the delamination growth of composite laminates by comparison; Finally, we validate the proposed cohesive model and VCCT for the DCB, ENF, MMB and FRMM cases respectively, and also explore the effects of the cohesive strengths on the global load-displacement response and damage evolution locuses.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
2 Boundary value problem with an interface discontinuity Let us consider a bounded domain
3
occupied by an initially linear elastic body, as
shown in Figure 1. The boundary of the domain with the outward normal n is subdivided into either Neumann or Dirichlet boundary conditions defined on t and u , respectively,
t and u
such that u
internal discontinuity d such that
t =Ø. Further, is divided into and by an
in the reference configuration. The
displacement field u across d is discontinuous but continuous on and . Neglecting the inertia and body force, the strong form is given by the balance of linear momentum and the boundary conditions
0 , u ub ,
on u
n Tb ,
on
t
on
d
in
N t u , where ub : u
Tb : t
3
3
(1)
is the field of prescribed displacement on the Dirichlet boundary and
is the field of prescribed traction on the Neumann boundary. The Cauchy stress
tensor is denoted by
. The Cauchy cohesive traction t which acts on the internal boundary
d depends on the displacement jump u at the interface d . The reference normal vector
N points from to . The virtual work equation for the body with volume V containing the cohesive interface d is given by the following weak form
X u : dV t u dA u Tb dA d
t
(2)
3 Cohesive modeling for mixed-mode delamination of composite laminates The exponential cohesive traction-displacement jump relationship governing the interface
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
failure is used, as shown in Figure 2. The expression for each failure mode is given by
ti etic / ui
c
ui exp ui / ui
c
,
(i 1, 2, 3)
(3)
where i=1, 2 and 3 represent the normal and two tangential directions, respectively. The fracture toughness for each fracture mode is equal to the area under the tractiondisplacement jump curve
0
0
t1d u1 et1c u1
c
G1c
t2 d u2 t3d u3 et2c u2
(4) c
0
et3c u3
c
G2c G3c
(5)
where G1c , G2c and G3c are the fracture toughness for fracture mode-I, II and III, respectively. 3.1 Single-mode delamination model based on continuum damage mechanics Based on Eq.(3), an anisotropic damaged cohesive traction-displacement jump relationship is proposed
ti etic / ui
c
u 1 ui exp i c , (i 1, 2, 3) ui 1 di
(6)
where d is the damage vector. The Helmholtz free energy is written as (t , u , d , T ) ( , T ), ( ) e
d
d
1 2
a
(7)
where and are the free energy representing the elastic deformation and damage evolvement e
d
of interface, respectively. T is the temperature and is a vector of internal variable. a is a material parameter. According to the thermodynamic theory, the relationships between the thermodynamic conjugate forces
t
Y , B and internal variables d , are given by
, Y , B a u d
(8)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
The Clausius-Duhem dissipation inequality imposed by the continuum thermodynamics on the evolution laws is [29] 1
t u Y d B T q 0
(9)
T
where is the power of dissipation due to damage. q is the heat flux. If the dissipation reaches the maximum, the damage evolution law is given by
F d F d 0 d d , 0 d Y Y B B
(10)
where F d (Y , d , B) is the damage potential function and is the damage consistency parameter. d
They are assumed to obey the Kuhn-Tucker loading/uploading conditions [29]
d 0, F d 0, d F d 0
(11)
Here, a damage potential function F d (Y , d , B) is proposed for each failure mode if the temperature effect is neglected
F d (Y , d , B) Y Yc B Bc
(12)
From Eqs.(11) and (12), the damage criterion is not satisfied and d = 0 holds at F 0 d
while the further damage appears and d > 0 holds at F 0 . Thus, the consistency condition d
F d 0 for the damage evolution by using Eq.(12) derives the consistency parameter d .
F d . F d B Y d Y/ Y B a d
(13)
where is the unit vector. Here, an exponential function for the conjugate force Y is proposed when the displacement jump based failure criterion is assumed
u c 1 exp u f u Y / a 1 exp(1) where
c
u and
u
f
c
(14)
are the initial and critical failure displacement jump, respectively. By
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
combining Eqs.(11)-(14), the damage variable d is found to be equal to the right-hand formula of Eq.(14). The progression of damage is possible if the actual displacement jump u is equal to the evolving internal variable and increases for different delamination modes. The damage
u is smaller than or
remains constant if
conditions for the load function g ( u
d > 0, if d = 0, if
et c / u et
c
are given by
u > 0,
Damage evolution
u or
u 0,
No damage evolution
f
0
) u
u and
f
The critical displacement jump u u
u is not increasing. Thus, the damage evolution
c
u et c
c
u
c
u
f
in Eq.(14) is numerically solved by Eq.(16)
d u exp u
u exp u / u
(15)
c
f
/ u
c
G
where is a constant approximating unity. For example,
(16) c
u
f
11.76 u is obtained at c
=0.9999 by combining Eqs.(4) or (5) and (16). 3.2 Mixed-mode delamination model based on continuum damage mechanics The traction-displacement jump relationship for mixed-mode delamination is assumed by
t et c / u
c 2 1
c
1 u u exp c u 1 d
t t t
where t
u
c
u1
2
c 2 2
u2
c 2 3
2
and
u3
2
(17)
u
c
u
c 2
1
u2
u
c 2
c 2
3
.
and denotes the MacAuley bracket.
For mixed-mode delamination, the anisotropic damage variable d in Eq.(6) is replaced by the scalar damage variable d in Eq.(18) u c 1 exp u f u d 1 exp(1)
c
The damage evolution conditions for mixed-mode delamination are given by
(18)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
d > 0, if u and u > 0, d = 0, if u or u 0,
Damage evolution
(19)
No damage evolution
The mixed-mode interlaminar fracture toughness G c is assumed to obey the B-K law [30]
G G G G G shear GT c
c 1
c 2
c 1
(20)
where Gshear G2 G3 and GT G1 G2 G3 . is a constant depending on the mode mixity. f
u
The mixed-mode critical failure displacement jump
can be numerically solved by
combining Eqs.(20) and (21)
u
0
f
et c / u
et u c
u exp u / u
c
c
et
c
u
c
u
d u exp u
c
f
f
/ u
c
G
(21) c
tan
The tangent stiffness Dij for mixed-mode interface fracture is derived from Eqs.(6) and (18) by considering the interpenetration problem between the delaminated surface
Dijtan
ti uj
gij 1 f g1 j gij 1 f g1 j
/ uj 1 uj / uj
uj uj
uj 1 uj
c
1 u j / u j
t cj e
c
/ uj
c c uj 1 d uj 1 1 d d , c 2 u j u j 1 d
tc / uj j u j
0
c
exp 1
uj
c
0 u u
,
1
u
c
u 0 or
u u u u
(22)
where
d d u , uj u u j
d u
1 f
u
c
i 1, 2, 3
u c exp u f u 1 exp(1)
ui u u i 1 g1i ui u ui
c
c
(23)
(24)
(25)
f
f
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
and f is a penalty factor and g ij is the Kronecker dalta.
1 is for damage loading and u
0 is for damage uploading. u
3.3 Numerical considerations for cohesive models using finite element analysis The finite element implementation of cohesive model sees Segurado and LLorca [31]. The finite element analysis of proposed cohesive model is performed using ABAQUS-UEL (User element subroutine). The zero-thickness interface element is used and the middle-plane d is used for numerical interpolation. In the following, the mesh size effect and numerical convergence using finite element analysis are addressed respectively. 3.3.1 Solution for the mesh size effect Yang and Cox [12] pointed out a necessary condition for mesh-independent numerical results is that the length of a cohesive zone should be larger than the size of the solid element on either side of crack. For the triangular cohesive model, Turon et al. [32] explored the size effect by introducing the cohesive characteristic length l cz and adjusting the cohesive strength. Yet, the problem for the nonlinear cohesive model is more complicated due to variable interface stiffness. Here, the thought of crack band model [28] is introduced into the exponential cohesive model to solve the size effect. The objective response of the finite element analysis is ensured by regularizing the fracture toughness G1c in Eq.(4) and G2c in Eq.(5)
g c G c / l cz
(26)
c
where g is the energy dissipation per unit volume. The cohesive length l cz was evaluated that is defined as one over which stress increases up to the maximum interfacial strength some distance ahead of the crack tip [5,6,33-38]. For different fracture modes, the generalized formula is given by
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
l cz
EG c
tc
2
(27)
where E is Young’s modulus and is a constant depending on different cohesive models and failure cases. For practical delamination cases, the cohesive lengths are evaluated by adjusting the values for t c and . Here, two points should be noted: ① the cohesive length for mode-II is larger than that for mode-I since the delaminated fracture toughness for mode-II is generally larger than that for mode-I; ② for the mixed-mode delamination, the cohesive length for mode-I and mode-II is taken respectively. 3.3.2 Solution for the convergence problem The convergence problem arises from the softening effect of cohesive model. Essentially, the negative eigenvalues in the stiffness matrix lead to the ill-conditioned finite element equations and unstability of numerical solution. In order to solve the convergence problem, Camanho et al. [15], and Alfano and Crisfield [17] used the local-control arc-length and line searches, Hamitouche et al. [39] introduced the damping effect into cohesive zone model. Specially, Gao and Bower [40] introduced the viscous effect into the Xu and Needleman’s cohesive model [24]. Recently, Liu et al. [41] gave a review on the convergence problem and found the viscous damping technique shows higher calculation efficiency than the arc-length algorithm since the latter imposes additional constraint equations between the displacement increment and the arc-length radius. Here, the viscous damping method is used for stabilizing the numerical solution.
4 VCCT for predicting the delamination growth of composite laminates Figure 3 shows crack propagation from the crack tip node i to j with the increment crack length Δa. Node i is separated into two nodes i1 and i2 after crack propagation. For node i, the
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
relative displacements in three directions (x, y, z) are Δuix, Δuiy and uiz after propagation and the node forces before propagation are Fix, Fiy and Fiz. The energy release rate due to crack propagation is given by [2,3]
G GI GII GIII 1 a 0 2 S lim
a
0
Fix (a)uix (a)da Fiy (a)uiy (a)da Fiz (a) uiz (a)da 0 0 a
a
(28)
where GI , GII and GIII are energy release rate for the mode-I, mode-II and mode-III, and S is the new area generated due to the crack propagation length Δa. Assume a two-step method is adopted by using the node force before crack propagation and node relative displacement after crack propagation, the G due to crack propagation on the crack closure surface is calculated as [2,3]
G
1
Fi x u i x Fi y u i y F 2S
u
iz
iz
(29)
Often, the two-step method can be approximately substituted by the one-step method if the mesh sizes on the crack closure surface are sufficiently small. The relative displacement after crack propagation can be substituted by the relative displacements between nearest node pairs before crack propagation [42]. For example, the relative displacements in three directions for the node pair i1 and i2 can be approximately substituted by those for the node pair e1 and e2. In this analysis, the one-step method is used for the calculation of energy release rate. The master-slave node technique is used to simulate crack propagation in the finite element analysis. The technique divides the node pair at the same position into two nodes by releasing the coupling freedom degree if the critical energy release rate G is larger than the equivalent one Gc according to the B-K law in Eq.(20). The finite element implementation of the VCCT has been integrated into ABAQUS finite element software as a module. In this research, the delamination analysis of composite laminates using the VCCT adopts this module.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
5 Numerical examples and discussions The DCB, ENF, MMB and FRMM cases are used to testify the cohesive model, in which the DCB case belongs to the mode-I delamination type, the ENF case belongs to the mode-II delamination type, and the MMB and FRMM cases belong to the mixed-mode delamination type. In order to study the mesh size effect on the load-displacement response using the cohesive model and VCCT, four mesh sizes with different element side length along the direction of crack propagation h=0.015mm, 0.03mm, 0.1mm and 0.5mm are applied. The number of bulk elements is 2E5, 1E5, 3E4 and 6E3, respectively. The Newton-Cotes integration is used [43] and two integration points are set at each direction for each cohesive element, and the displacement, traction and damage variable as the status variables that depend on the solution are monitored. The contact pairs between the delaminated surface are initially set to prevent the normal interpenetration. The Newton-Raphson algorithm is used and the time increment is properly controlled. in Eq.(27) is taken 0.5 and 0.1 respectively for the mode-I and mode-II delamination. 5.1 Numerical results for the DCB delamination case Figure 4 shows the DCB case with boundary and loads. Two loads are symmetrically and oppositely applied to the upper lamina and low lamina of specimen. The material and size parameters are listed in Table 1 [23]. Figure 5 shows three mesh models for the DCB case and the number of cohesive interface elements is 4668, 700 and 140 for h=0.015mm, 0.1mm and 0.5mm mesh sizes, respectively. Figure 6 shows the crack length a-crack opening displacement 2w curves for the DCB case for three mesh models at t1c 60MPa . The ratio a/2w changes from about 12.5 to 4.5. The delamination crack propagates quickly at the beginning, but then enters the unstable
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
delamination stage until the end of the specimen. Figure 7 shows the global load-displacement curves by considering the effects of the cohesive strengths and mesh sizes. The load first increases and then decreases with the displacement. It can be seen the curves at low cohesive strengths using cohesive model deviate from the true results, but the curve at t1c 60MPa is in good agreement with those by the VCCT and discrete cohesive model at 5.7MPa cohesive strength by [23]. Yet, the peak load using the VCCT is largest. The error may arise from the evaluation of energy release rate as the failure criterion for the VCCT. Spurious numerical oscillation appears using both the cohesive model and VCCT at the softening stage for coarse mesh size, but coarse and fine mesh models lead to consistent results, showing the proposed cohesive model is insensitive to the mesh size. The main reason for numerical oscillation arises from the sudden release of strain energy in the intralaminar materials which leads to instantaneous failure of the element, and thus spurious solution jump appears in the unstable snap-back delamination process. Since the practical cohesive model must reflect the stress distribution associated with damage mechanisms occurring ahead of the physical crack tip [44], the numerical oscillation shows the coarse mesh size cannot accurately capture the cohesive zone stress distribution. However, the curve becomes smooth at the decreasing stage when the mesh size h is equal to the cohesive length l1cz . Xie and Waas [23] also explored the effect of the bias mesh on the curve and found it can more easily lead to oscillation using coarse and bias mesh sizes. According to the present numerical results, only the element side length h along the direction of crack propagation shows large influence on the numerical results. Besides, the initial crack propagation is shown to appear much earlier than the peak load. This conclusion is also consistent with that obtained by Chen et al. [14]. For the mode-I delamination, no convergence problem occurs without the viscous technique.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
5.2 Numerical results for the ENF delamination case Figure 8 shows the ENF case with boundary and load. Three-point bending with an interlaminar initial crack is generally used to test the mode-II fracture toughness. The contact pair is introduced into the delaminated surface. The material and size parameters are listed in Table 2 [23]. The number of cohesive interface elements is 2334, 700 and 140 for h=0.03mm, 0.1mm and 0.5mm mesh sizes, respectively. Figure 9 shows the global load-displacement curves. By using low cohesive strengths and coarse mesh sizes, the load-displacement curves by the cohesive model differ from the results using the VCCT and the discrete cohesive model at 57MPa cohesive strength by Xie and Waas [23]. The load-displacement curve experiences three stages: ① the load first increases with displacement; ② a large ‘snap’ phenomenon appears indicating the unstable crack propagation; ③ and then the load increases rapidly again. Goyal et al. [25] classified this case into the unstable delamination growth at a0 0.35L for the ENF specimen. It can also be seen the curves are basically consistent using three mesh models. Different from the DCB case, the convergence becomes more difficult using the Newton-Raphson algorithm for the ENF case which requires the viscous regularization and linear search algorithm as the cohesive strength increases to 500MPa. The convergence problem due to the increase of the cohesive strength and interface stiffness is also found by Chen et al. [14] and Turon et al. [19]. In addition, Feih [45] found the restriction on the element size can be improved by choosing a larger number of integration points instead of reducing the element sizes. Yet, this still requires sufficiently fine mesh sizes since the combination of a large number of integration points and coarse mesh sizes may also lead to numerical oscillation. The oscillatory performance of interface elements in terms of the selection of integration points is specially studied by Schellekens and de Borst [43]. They
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
showed Gaussian integration can easily lead to the oscillation of traction when large traction gradients are present over an element. Thus, the Newton-Cotes integration is used in this research. 5.3 Numerical results for the MMB delamination case Figure 10 shows the MMB case with boundary and loads. The experimental equipment of the MMB is initially proposed by Reeder and Crews [46] to predict the mixed-mode I/II delamination fracture toughness. By loading with a lever, a single applied load produces bending loads on the specimen. Different mode I/II ratios can be obtained by adjusting the length e of the lever. The material and size parameters are listed in Table 3 [13,23]. The number of cohesive interface elements is 700 and 140 for h=0.1mm and 0.5mm mesh sizes, respectively. Figure 11 shows the load-displacement curves for the MMB case. Similar to the DCB and ENF cases, the curve shows small numerical oscillation at the softening stage for coarse mesh sizes using the cohesive model and VCCT. As the cohesive strength t increases, the initial stiffness et c u of cohesive c
c
interface increases rapidly and differs for normal and shear failure. The initial normal interface stiffness reaches about 5E4N/mm at t1c 500MPa . As the cohesive strength increases to 500MPa, the curve is close to that obtained by the VCCT. Same to the ENF case, the curve for the MMB case also exhibits three distinct stages, in which the sudden ‘snap’ effect shows unstable delamination properties. The delamination appears just at the beginning of loading, much earlier than the moment at which the peak load reaches. The numerical results by the proposed cohesive model are also compared with the results using the interface delamination model and cohesive model [7,13,23]. In general, these curves show good agreement except small errors at the last increasing stage. Compared with the B-K law used by the present work, the power law is given by [13,22, 23,25,27]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
G1 G2 c c 1 G1 G2
(30)
where 1 and 2 denote the linear and quadratic failure criteria, respectively. The power law in Eq.(30) fails to accurately capture the dependence of the mixed-mode fracture toughness on the mode ratio [15]. In contrast, the difficulty using the B-K law is to determine the coefficient in Eq.(20) which yet depends on the mode-mixity and the phase angle between the tangent and normal displacement jumps [30]. Turon et al. [18] explored the effect of different mode ratios on the load-displacement curves and got the value 2.284 in their cohesive models by applying a least-square fit to the experimental data for AS4/PEEK carbon fiber composites. In this research, the determination for is to calculate Eq.(18). We found
u
f
/ u
c
varies little as changes. For example,
u
u f
f
/ u
c
in
c
/ u changes
from 2.1 to 2.3 as varies from 1.5 to 2.5 at t1c t2c 500MPa . Thus, 2.0 is used for the present work. Figure 12(a) shows the damage variable d-mixed mode displacement jump
u curves for three cohesive strengths and Figure 12(b) shows the damage evolution locuses at t1c t2c 500MPa . As the cohesive strength increases, the ability for the interface to resist the
delamination becomes stronger. The approximate elliptic damage evolution locus according to the present cohesive model is similar to that obtained by Turon et al. [18]. The only difference is that the present cohesive model requires higher cohesive strengths than that by Turon et al. [18] in order to get accurate numerical results since a nonlinear exponential cohesive law instead of a simpler triangle relationship is used. Similar to the proposed cohesive model, the conclusion for the requirement for the cohesive strength to increase to high values such as 500MPa in order to lead to accurate results is also obtained by the potential-based exponential cohesive model [26]. Besides, the numerical convergence is also a large challenge for both the triangle and exponential
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
cohesive models for softening problems. For the exponential cohesive models, high cohesive strength adds the difficulty of convergence. Yet, the convergence problem for the triangle cohesive models lies in the sharp and discontinuous transition at the peak traction. Although the artificial viscous method is an ideal one for improving the convergence, the viscous coefficients should be properly chosen, similar to the analysis by Gao and Bower [40]. 5.4 Numerical results for the FRMM delamination case Figure 13 shows the FRMM case with boundary and load. The material and size parameters are listed in Table 4 [14]. Different from the MMB case, the mode-I and II delamination fracture toughness for the FRMM case is different. The number of cohesive interface elements is 1834, 550 and 110 for h=0.03mm, 0.1mm and 0.5mm mesh sizes, respectively. Figure 14 shows the load-displacement curves for the FRMM case. Figure 15 shows the damage evolution locuses. Same to the DCB, ENF and MMB cases, the delamination starts much earlier than the peak load. Similar to the ENF and MMB cases, three stages appear in the process of delamination, accompanied by the sudden softening effect. It can be seen all curves for the DCB, ENF, MMB and FRMM cases are close to the true results as the cohesive strength increases. This conclusion is c
c
also obtained by Xie and Waas [23]. At t1 t2 500MPa , the curve using proposed cohesive model is basically consistent with the results obtained by the experiments, the VCCT and Chen et al. [14]. Also, they studied the effect of parameter in Eq.(30) and found the results using the quadratic failure criterion 2 are more close to the experimental results. In addition, Chandra et al. [8] gave an overview on the selection of cohesive strength (or displacement jump) for current cohesive models and applied cases. The traction-displacement cohesive properties between bi-material interface can be determined by inverse optimization computation by measuring the
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
crack-tip opening displacements directly [47-48]. Beisdes, such as the infra-red crack opening interferometry measurements [49] and laser treatment technique [50] are also used to determine the cohesive properties and to improve the bond toughness of interface. In terms of the sensitivity analysis, the effects of cohesive parameters including the cohesive fracture energy, cohesive shape and cohesive strengths on the mixed-mode delamination mechanisms using cohesive models are studied. The cohesive models as a phenomenological method is expected to expound the true bonding properties for various bi-material interface by adjusting the cohesive zone parameters.
6 Concluding remarks In recent years, the cohesive theory has been increasingly used in the research on the delamination problem of composite laminates. Various cohesive models were developed to predict the delamination initiation and growth. However, until now no universal and numerically robust cohesive model exists for the efficient and accurate simulation for the interface failure. In this research, an anisotropic exponential cohesive model based on continuum damage mechanics is proposed to predict the mixed-mode delamination response of composites. By performing numerical simulation on the DCB, ENF, MMB and FRMM delamination cases, the global load-displacement curves using proposed cohesive model are shown to be in good agreement with the results obtained by using the VCCT and existing cohesive models after evaluating the cohesive characteristic length accurately. It is demonstrated the nonlinear cohesive model can solve the mesh size effect and convergence problems effectively. Using the proposed cohesive model, even the coarse mesh sizes can lead to consistent global load-displacement response. The proposed nonlinear cohesive model is expected to be applied to more delamination problems.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Acknowledgement (1) The author Dr. P.F. Liu at Zhejiang University as a Visiting Scholar at Princeton University would like to thank Prof. Jean H. Prévost at Princeton University for helpful academic guidance and discussions. (2)This research is supported by Natural Science Funding of Zhejiang Province of China (Number: LY13E050002).
References [1] Ladeveze P, Lubineau G, Marsal D. Towards a bridge between the micro- and the mesomechanics of delamination for laminated composites. Compos Sci Technol 2006; 66:698-712. [2] Rybicki EF, Kanninen MF. A finite element calculation of stress intensity factors by a modified crack closure integral. Eng Fract Mech 1977; 9: 931-938. [3] Krueger R. Virtual crack closure technique: history, approach, and applications. Appl Mech Rev 2004; 57: 109-143. [4] Liu PF, Hou SJ, Chu JK, Hu XY, Zhou CL, YL Liu, et al. Finite element analysis of postbuckling and delamination of composite laminates using virtual crack closure technique. Compos Struct 2011; 93:1549-1560. [5] Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids 1960; 8: 100-104. [6] Barenblatt I. Mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech 1962; 7: 55-129. [7] Mosler J, Scheider I. A thermodynamically and variationally consistent class of damage-type cohesive models. J Mech Phys Solids 2011; 59: 1647-1666. [8] Chandra N, Li H, Shet C, Ghonem H. Some issues in the application of cohesive zone models for metal-ceramic interfaces. Inter J Solids Struct 2002; 39: 2827-2855. [9] Elices M, Guinea GV, Gómez J, Planas J. The cohesive zone model advantages, limitations and challenge. Eng Fract Mech 2002; 69: 137-163.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[10] Cornec A, Scheider I, Schwalbe KH. On the practical application of the cohesive model. Eng Fract Mech 2003; 70 :1963-1987. [11] de Borst R. Numerical aspects of cohesive-zone models. Eng Fract Mech 2003; 70: 1743-1757. [12] Yang QD, Cox BN. Cohesive models for damage evolution in laminated composites. Int J Fract 2005; 133: 107-137. [13] Mi Y, Crisfield MA, Davies GA. Progressive delamination using interface elements. J Compos Mater 1998; 32:1246-1272. [14] Chen J, Crisfield M, Kinloch AJ, Busso EP, Matthews FL, Qiu Y. Predicting progressive delamination of composite material specimens via interface elements. Mech Compos Mater Struct 1999; 6:301-317. [15] Camanho PP, Davila CG, De Moura MF. Numerical simulation of mixed-mode progressive delamination in the composite materials. J Compos Mater 2003; 37:1415-1438. [16] Alfano G, Crisfield MA. Finite element interface models for the delamination analysis of laminated composites: mechanics and computational issues. Int J Numer Meth Eng 2001; 50: 1701-1736. [17] Alfano G, Crisfield MA. Solution strategies for the delamination analysis based on a combination of local-control arc-length and line searches. Int J Numer Meth Eng 2003; 58: 999-1048. [18] Turon A, Camanho PP, Costaa J, D´avila CG. A damage model for the simulation of delamination in advanced composites under variable-mode loading. Mech Mater 2006; 38: 1072-1089. [19] Turon A, Camanho PP, Costa J, Renart J. Accurate simulation of delamination growth under mixed-mode loading using cohesive elements: Definition of interlaminar strengths and elastic stiffness. Compos Struct 2010; 92:1857-1864. [20] Swaminathan S, Pagano NJ, Ghosh S. Analysis of interfacial debonding in three-dimensional composite microstructures. J Eng Mater Technol 2006; 128: 96-106. [21] Bui VQ, Iannucci L, Robinson P, Pinho ST. A coupled mixed-mode delamination model for laminated composites. J Compos Mater 2011; 45: 1717-1729. [22] Jiang WG, Hallett SR, Green BG, Wisnom MR. A concise interface constitutive law for
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
analysis of delamination and splitting in composite materials and its application to scaled notched tensile specimens. Int J Numer Meth Eng 2007; 69:1982-1995. [23] Xie D, Waas AM. Discrete cohesive zone model for mixed-mode fracture using finite element analysis. Eng Fract Mech 2006;73:1783-1796. [24] Xu XP, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids 1994;42:1397-1434. [25] Goyal VK, Johnson ER, Davila CG. Irreversible constitutive law for modeling the delamination process using interfacial surface discontinuities. Compos Struct 2004; 65: 289-305. [26] Park K, Paulino GH, Roesler JR. A unified potential-based cohesive model of mixed-mode fracture. J Mech Phys Solids 2009; 57: 891-908. [27] El-Sayed S, Sridharan S. Predicting and tracking interlaminar crack growth in composites using a cohesive layer model. Compos Part B 2001; 32:545-553. [28] Bažant ZP, Oh B. Crack band theory for fracture of concrete. Mater Struct 1983;16:155-177. [29] Lemaitre J, Chaboche, JL, 1990. Mechanics of Solid Materials. Cambridge University Press, Cambridge. [30] Benzeggagh ML, Kenane M. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Compos Sci Technol 1996; 56: 439-449. [31] Segurado J, LLorca J. A new three-dimensional interface finite element to simulate fracture in composites. Inter J Solids Struct 2004; 41: 2977-2993. [32] Turon A, Da´vila CG, Camanho PP, Costa J. An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Eng Fract Mech 2007; 74: 1665-1682. [33] Irwin GR. Plastic zone near a crack and fracture toughness. Proceedings of the seventh Sagamore Ordnance materials conference, vol. IV. New York: Syracuse University; 1960. p. 63-78. [34] Hillerborg A, Mode´er M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement Concr Res 1976; 6: 773-782.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[35] Rice JR. The mechanics of earthquake rupture. In: Dziewonski AM, Boschhi E, editors. Physics of the Earth’s interior. Proceedings of the international school of physics ‘‘Enrico Fermi’’, Course 78, 1979. Amsterdam: Italian Physical Society/North-Holland; 1980. p. 555-649. [36] Bao G, Suo Z. Remarks on crack-bridging concepts. Appl Mech Rev 1992; 24, 355-366. [37] Cox BN, Marshall. Concepts for bridged cracks in fracture and fatigue. Acta Metallurgica et Materialia 1994; 42: 341-363. [38] Harper PW, Hallett SR. Cohesive zone length in numerical simulations of composite delamination. Eng Fract Mech 2008;75:4774-4792. [39] Hamitouche L, Tarfaoui M, Vautrin A. An interface debonding law subject to viscous regularization for avoiding instability: Application to the delamination problems. Eng Fract Mech 2008; 75: 3084-3100. [40] Gao YF, Bower AF. A simple technique for avoiding convergence problems in finite element simulations of crack nucleation and growth on cohesive interfaces. Model Simul Mater Sci Eng 2004; 12(3): 453-463. [41] Liu PF, Zheng JY. Recent developments on damage modeling and finite element analysis for composite laminates: a review. Mater Des 2010; 31: 3825-3834. [42] Xie D, Biggers SB. Progressive crack growth analysis using interface element based on the virtual crack closure technique. Finite Elem Anal Design 2006; 42: 977-984. [43] Schellekens JCJ, de Borst R. On the Numerical Integration of Interface Elements. Inter J Numer Meth Eng 1993; 36: 43-66. [44] Schellekens JCJ, de Borst R. A nonlinear finite-element approach for the analysis of mode-I free edge delamination in composites. Int J Solids Struct 1993;30:1239-1253. [45] Feih S. Development of a user element in ABAQUS for modelling of cohesive laws in composite structures. Risø National Laboratory Roskilde, Denmark, Risø-R-1501, 2005. [46] Reeder JR, Crews Jr JR. Mixed-Mode Bending Method for Delamination Testing. AIAA J 1990; 28: 1270-1276. [47] Jacobsen TK, Sorensen BF. Mode I intra-laminar crack growth in composites-modelling of R-curves from measured bridging laws. Compos Part A 2001; 32: 1-11. [48] Shen B, Stanciulescu I, Paulino GH. Inverse computation of cohesive fracture properties from
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
displacement fields. Inver Probl Sci Eng 2012; 18: 1103-1128. [49] Gowrishankar S, M HX, Liechti K M, Huang R. A comparison of direct and iterative methods for determining traction-separation relations. Int J Fract 2012;177: 109-128. [50] Alfano M, Lubineau G, Furgiuele F, Paulino GH. On the enhancement of bond toughness for Al-epoxy T-peel joints with laser surface treated surfaces. Int J Fract 2011;171:139-150.
Figure captions: Figure 1 Boundary value problem with an interface discontinuity Figure 2 Exponential traction-displacement cohesive law Figure 3 Schematic representation of the VCCT: (a) before crack propagation and (b) after crack propagation Figure 4 The DCB delamination model with boundary and loads Figure 5 Three mesh models for the DCB case Figure 6 The crack length-crack opening displacement curves for the DCB case for three mesh models at t1c 60MPa Figure 7 The load-displacement curves for the DCB case: (a) effects of the cohesive strengths using coarse mesh size h 33l1cz 0.5mm and (b) effects of the mesh sizes at
t1c 60MPa . The curves using proposed cohesive model are also compared with the results obtained using the VCCT and existing cohesive models Figure 8 The ENF delamination model with boundary and load Figure 9 The load-displacement curves for the ENF case: (a) effects of the cohesive strengths using coarse mesh size h 33l2cz 0.5mm and (b) effects of the mesh sizes at
t2c 500MPa . The curves using proposed cohesive model are also compared with the results obtained using the VCCT and existing cohesive models Figure 10 The MMB delamination model with boundary and loads Figure 11 The load-displacement curves for the MMB case: (a) effects of the cohesive strengths using coarse mesh size h 5l1cz 0.5mm and (b) effects of the mesh sizes at
t1c t2c 500MPa . The curves using proposed cohesive model are also compared with the results obtained using the VCCT and existing cohesive models Figure 12 Damage evolution of cohesive elements for the MMB case: (a) damage variable
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
d-mixed mode displacement jump u curve for three cohesive strengths and (b) damage evolution locus at t1c t2c 500MPa . Figure 13 The FRMM delamination model with boundary and load Figure 14
The load-displacement curves for the FRMM case: (a) effects of the cohesive strengths using coarse mesh h 12.5l1cz 0.5mm and (b) effects of the mesh sizes at
t1c t2c 500MPa . The curves using proposed cohesive model are also compared with the results using the VCCT, existing cohesive model and experiment Figure 15 Damage evolution of cohesive elements for the FRMM case: (a) damage variable d-mixed mode displacement jump u curve for three cohesive strengths and (b) damage evolution locus at t1c t2c 500MPa
Table 1 Material and geometry parameters for the DCB case [23] E2 (GPa) G12 (GPa)
E1 (GPa)
135.9
9
5.2
c
v12
G1 ( N / mm) L(mm) a0 (mm) 2h(mm) B(mm)
0.24
0.28
100
30
3
10
Table 2 Material and geometry parameters for the ENF case [23] E1 (GPa)
70
E2 (GPa)
70
c
v12
G2 ( N / mm) L(mm) a0 (mm) 2h(mm) B(mm)
0.24
1.45
100
30
3
10
Table 3 Material and geometry parameters for the MMB case[13,23] c
c
E1 (GPa) E2 (GPa) G12 (GPa) v12 G1 G2 ( N / mm) L(mm) e(mm) a0 (mm) 2h(mm) B(mm)
135.9
9
5.2
0.24
4
100
50
30
3
10
Table 4 Material and geometry parameters for the FRMM case [14] c
c
E1 (GPa) E2 (GPa) G12 (GPa) v12 G1 ( N / mm) G2 ( N / mm) L(mm) a0 (mm) 2h(mm) B(mm)
176
8
6
0.27
0.257
0.856
105
45
3.1
24
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure
Figure