Decohesion finite element with enriched basis functions for delamination

Composites Science and Technology 69 (2009) 2616–2624

Contents lists available at ScienceDirect

Composites Science and Technology journal homepage: www.elsevier.com/locate/compscitech

Decohesion finite element with enriched basis functions for delamination I. Guiamatsia *, J.K. Ankersen, G.A.O. Davies, L. Iannucci Imperial College London, Dpt. Aeronautics, Prince Consort Rd., South Kensington SW7 2AZ, UK

a r t i c l e

i n f o

Article history: Received 12 May 2009 Received in revised form 20 July 2009 Accepted 2 August 2009 Available online 8 August 2009 Keywords: C. Delamination B. Modelling C. Damage mechanics Mesh size

a b s t r a c t In order to address the stringent mesh size requirement of the cohesive element, currently 0.5 mm or less, the mixed-mode interface finite element is enriched with the analytical solution of an idealized beam on elastic foundation. Both the interface and the solid continuum elements are enriched; the partition of unity (PU) method is utilized to obtain the enhanced interpolation, implemented with two user elements (UEL) in the commercial package ABAQUS/Standard. The new formulation yields predictions in excellent agreement with theoretical and experimental results for a typical mode I delamination benchmark, with elements 5 mm in size. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction Delamination is the most common mode of damage, and arguably one of the most fatal for laminated composite structures. The simulation of delamination with a finite element approach is eased by the knowledge of the potential damage path – the interface – and is usually performed with one of the following two approaches:  The virtual crack closure technique (VCCT), Rybicki and Kanninen [1], allows to calculate the energy released at the crack tip with a single calculation, based on the assumption of self-similar propagation. The modes I and II energies can be distinguished and the energy is given as the product of internal reactions at the crack tip and the difference between the displacements of selected nodes in the delaminated area before the tip. The released energy is compared to the critical value or fracture toughness Gc , and the crack is advanced. The crack advancement used to require the generation of a new finite element (FE) model and involved computationally costly remeshing operations, and was the greatest limitation of the technique. Krueger [2] gives an extensive review of the method and its finite element (FE) implementation. In recent years, however, VCCT has recaptured the interest of researchers and interesting solutions to avoid remeshing have surfaced: the VCCT-based interface element, Mabson et al. [3] and reverse contact in ABAQUS/Standard [10].

* Corresponding author. Tel.: +44 20 7594 5118; fax: +44 20 7584 8120. E-mail addresses: [email protected], [email protected] (I. Guiamatsia). 0266-3538/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compscitech.2009.08.002

 The cohesive zone model (CZM) was first introduced by Dugdale [8] and Barenblatt [9], with the idea of cohesive tractions acting within the plastic deformation zone ahead of a crack tip. The finite element (FE) implementation by Hillerborg et al. [4] was improved by Schellekens and De Borst [5] who were the first to propose a separate element for the interface. Further contributions from Wisnom et al. [23], Allix et al. [12,13], Mi et al. [20], Camanho [14] and others were mainly related to the consideration of mixed-mode propagation. Progress in non-linear solution procedures, Crisfield [18], have successfully addressed the numerical issues pertaining to the negative tangent stiffness caused by the interface material softening. The last battlefield, however, is that of the cohesive element size. Typically, the elements have to be less or equal to 0.5 mm in size in order to achieve convergence, as illustrated by the mesh sensitivity analysis shown in Fig. 1a, for a typical double cantilever beam (DCB) benchmark described in Ref. [7]. In contrast, VCCT, Fig. 1b, yields excellent results with much larger elements. The stringent mesh size requirement finds its rationale in the inherent concept of the cohesive element itself: because the cohesive law is only active when damage occurs, any interface element that is larger than the realistic cohesive zone (of the order of 1 mm) must be able to allow for active (damaged) and inactive (undamaged) regions within the same element. Fig. 2 plots the normal transverse stress for several points along the interface of a DCB specimen, showing steep gradients in the stress field ahead of the fracture process zone – within a region of approximately 4–5 mm – that cannot be accurately captured with a linear element.

2617

I. Guiamatsia et al. / Composites Science and Technology 69 (2009) 2616–2624

(a) Cohesive interface elements

(b) VCCT 120 linear theory 5mm vcct 2mm vcct 1mm vcct 0.5mm vcct

Reaction Force (N)

Reaction Force (N)

100 80 60 40 20 0 0

2

Applied displacement (mm)

4

6

8

10

Applied displacement (mm)

Fig. 1. Mesh sensitivity analysis for DCB Benchmark. (a) ABAQUS generic interface element. (b) ABAQUS VCCT based on reverse contact.

The aim of this work is to improve upon the current punitive cohesive element size requirement. This problem has previously been tackled by Turon et al. [16], who proposed the use of an artificially low interface strength. Although the strategy is able to bring the mesh size up to approximately 1 mm, the solution has notorious drawbacks such as inaccurate prediction of the limit point of the global load versus displacement curve. Additionally, problem-dependent tuning is required: Alfano and Crisfield [24] examined the difference in ‘ideal’ strengths between linear and quadratic elements, ABAQUS User’s Manual [10] reports results obtained for the two-crack DCB problem [11] and Harper and Hallett [17] recommend caution when using the strategy for mode II fracture. In this work, we propose to enrich the element’s basis functions with an analytical solution for the displacement near the crack tip, in order to improve the prediction and alleviate the mesh size requirement. It is worth reminding the reader that, compared to VCCT, the cohesive zone model is conceptually superior: it is a more generally applicable approach, capable of handling both initiation and propagation, with no initial flaw requirement. This paper is organized into three sections as follows: first, the formulation for the current linear interface element is presented, then the process of PU enrichment of a solid continuum element, as well as a cohesive interface element, is explained, and finally results obtained with a larger, enriched element are shown for a double cantilever beam test.

4.00E+07

2. Interface element formulation The complexity of the formulation of interface damage stems from the fact that the interface between plies is not a geometrically distinguishable entity, but rather a modelling artifact. In order to use continuum damage mechanics (CDM) concepts for the simulation of delamination, the interface is represented as a material that is initially elastic, then elastically damaged. 2.1. Element topology Because the constitutive equations for the interface are expressed as traction versus separation, the equivalent of strain for these elements is the separation between the nodes on the upper surface and those on the lower surface of the element: the element is hence directional. An interface element with linear topology is illustrated in Fig. 3. For the sake of simplicity, equations are written for a twodimensional (2D) model, but are easily extendable to the threedimensional (3D) case. In a 2D model, the interface element has a one dimensional formulation: the separation in the material direction i; i ¼ 1; 2 is written, in matrix form, as:

di ¼ ½ N10

N20

N 20

8 9 ui > > > > > 1i > = < u2 i ¼ 1; 2 N 10  > ui3 > > > > ; : i > u4

fracture process zone

3.00E+07

(MPa)

2.00E+07

crack tip

1.00E+07

~5mm 0.00E+00 0

0.002

0.004

0.006

0.008

0.01

0.012

0.014

0.016

0.018

(m)

-1.00E+07 -2.00E+07

-3.00E+07 Fig. 2. Variation of the transverse normal stress

r22 along a DCB specimen.

0.02

ð1Þ

2618

I. Guiamatsia et al. / Composites Science and Technology 69 (2009) 2616–2624

upper surface 3

4 +

+

+

+

-

-

-

-

1'

2' 2

1

lower surface

small / zero thickness

Fig. 3. Linear interface element. Left: embedded at an interface. Right: topology and node numbering.

In the above, N 10 ¼ N1 ¼ N4 and N 20 ¼ N 2 ¼ N 3 are the interpolation functions for a two-node linear element, and uij is the displacement in the ith direction at node j. As shown in Fig. 3, the element is connected on the upper and lower side to elements representing laminate’s plies.1 2.2. Constitutive equations Fracture toughness is the critical material property for fracture; it is embedded within the constitutive law of the interface element, as the area under a traction versus separation curve. The bilinear shape, illustrated in Fig. 4, is the most widely used, though other softening shapes exist. The plot expresses the function qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiof the stress in terms of the mixed-mode separation d ¼

d2I þ d2II . Before the initial separation

d0 ¼ r0 =K is reached, the material is purely elastic and the slope of the curve is the penalty stiffnessK that may, either be made as high as numerically possible or be assigned values based equating the interface to a thin resin rich layer. The difference between the two options is important to outline:  In the first case – infinite penalty stiffness – the initiation separation d0  df is zero, for all practical means. Hence, the element is a pure fracture element and is only driven by the definition of the fracture toughness. The concept of ‘fracture’ element is however slightly misleading as it implies that the interface strength, r0 , bears no influence on the outcome of the analysis, although the opposite has been observed, Turon et al. [16], Alfano and Crisfield [24].  In the second case – realistic penalty stiffness2 – it is assumed that the interface is a resin rich thin layer, which usually mimics the real structure. This approach is recommended here, along with using a realistic interface strength, for the sake of consistency. After initiation, the stress linearly reduces, and the stiffness is now calculated as Kð1  dÞ until d ¼ 1 when df is reached. In order to use a single damage variable d, the separation corresponding to two delamination modes – mode I and shear (modes II q and III) – ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi are combined as a single mixed-mode separation d ¼

d2I þ d2II .

This development is well detailed in Camanho et al. [14], Turon et al. [15], as well as Jiang et al. [22]. Some notes are worth making at this point: 1. Turon’s and Jiang’s models first differ by the constant proposed to represent the crack mixity ratio; we have, respectively.

1 These can be two-dimensional plane strain or plane stress continuum solids or in three-dimensional models, continuum or shell elements can be used. 2 Note that realistic values of interface stiffness are comparable, because of the very small realistic interface thickness, to values equivalent to 50 times the half beam stiffness, as proposed in Turon et al. [16].

dI dI þ dII dI ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 dI þ d2II

bTuron ¼ cosJiang

It seems that b performs better when the Benzeggagh–Kenane [21] failure criterion is used, while cos is more appropriate for the classical power law, especially with regard to the determination of the initial separation d0 . 2. Turon’s model introduced a new definition of the damage initiation, based on the B–K fracture criterion, replacing the traditional strength criterion. 3. Most formulations propose a degradation of the stiffness matrix

ri ¼ Dij dj ; Dii ¼ Kð1  dÞ

ð2Þ

and a hyperbolic damage evolution law:

d¼

df ðd  d0 Þ dðdf  d0 Þ

ð3Þ

The value of the damage parameter increases very rapidly from 0, which is required by the fact that the penalty stiffness is very high and needs to be brought down to a realistic value as soon as damage initiates. Jiang’s model, however, degrades the stresses and a linear evolution of the damage variable is proposed:

ri ¼ r0i ð1  dÞ i ¼ 1; 2

d¼

d  d0 df  d0

ð4Þ

Although simple, this formulation shows more difficulties with convergence in an implicit formulation, mainly caused by the non-symmetric tangent stiffness matrix. Turon’s model is utilized for the improvement proposed in this work. The commercial analysis package ABAQUS/Standard is used for this development and for validation purposes, an interface element based on Turon’s formulation and implemented via a user element (UEL) for ABAQUS/Standard is initially compared to ABAQUS generic interface element. Fig. 5 shows the plot of the predicted traction against separation of a double cantilever beam specimen, using elements of size 0.5 mm. At this point, a quick digression into the issue of integration is worthy of mention. The question of Gauss versus nodal integration has been dealt with in the past, Schellekens and De Borst [6], often to the conclusion that nodal quadrature yields better stability. However, it appears that nodal integration also impedes on convergence, as demonstrated in Fig. 6 which shows the prediction of the previous DCB benchmark using 1 mm mesh size for both integration schemes. Note that these results are obtained with the same UEL [15] and are also compared with the ABAQUS generic interface element. The prediction obtained with Gauss integration is in very good agreement with the theoretical solution, unlike that obtained with a nodal integration scheme (ABAQUS generic element and UEL).

2619

mode II

I. Guiamatsia et al. / Composites Science and Technology 69 (2009) 2616–2624

dfII

S0

e

stress S (MPa)

od

ed

df

d0II

d0

m

ix

m

d0

df

d0I

dfI

mode I

separation d (mm) Fig. 4. Interface constitutive law. Left: traction versus separation. Right: mixed-mode separation.

120

Reaction Force (N)

100

80

ABAQUS

UEL

Beam theory

60

40

20

0 0

0.002

0.004

0.006

0.008

0.01

0.012

0.014

0.016

0.018

0.02

Separation (m) Fig. 5. Reaction force against separation of benchmark DCB.

3. Enrichment Simple elements, hence interpolation functions are preferred in FE models not only because of their ability to represent any function accurately over sufficiently small domain intervals, but also because the enforcement of continuity, compatibility, and conformity requirements is generally straightforward. However, in cases where the shape of the solution is known, it may be advantageous to integrate that knowledge within the formulation in order to reduce the problem’s size, typically via enriching the shape functions. 3.1. Beam on elastic foundation As mentioned, the aim is to represent the steep variation of stress ahead of the crack tip by enriching the displacement interpolation functions. The idealization of the double cantilever beam (DCB) as a semi-infinite, isotropic, beam on elastic foundation, illustrated in Fig. 7, will provide the required analytical solution.3 The elastic foundation represents the midplane of the interface in a DCB problem. The equilibrium equation is: 3 This shape should also be appropriate for any thin plates debonding or delaminating.

EI

@ 4 v ðxÞ þ K v ðxÞ ¼ 0 @x4

ð5Þ

subject to boundary conditions:

v ð1Þ ¼ 0;

@ 3 v ð0Þ ¼0 @x3 3

In the above, E and I ¼ Bd are the Young’s modulus and second mo12 ment of area of the beam cross section; K ¼ Eti B is the elastic stiffness of the spring, with Ei the Young’s modulus of the resin, and t is half the interface thickness.4 Eq. (5) has a known general solution of the form:

v ¼ ebx ½A cosðbxÞ þ B sinðbxÞ þ ebx ½C cosðbxÞ þ D sinðbxÞ; which reduces in this specific case to

v ¼ v 0e

bx

½cosðbxÞ  sinðbxÞ;

rffiffiffiffiffiffiffi K 4 b¼ 4EI

ð6Þ

For the case of the DCB benchmark described in Ref. [7] – which will be used throughout this paper, the material properties of the 4

Typical values of the interface thickness are between 0.01 mm and 0.1 mm.

2620

I. Guiamatsia et al. / Composites Science and Technology 69 (2009) 2616–2624

120

100

Gauss_1mm

Nodal_1mm

ABQ_1mm

Beam Theory

Reaction (N)

80

60

40

20

0 0

0.005

0.01

0.015

0.02

Displacement (m) Fig. 6. Reaction force against separation of benchmark DCB – 1 mm mesh size.

v(x) B t

K d x

Fig. 7. DCB idealization as a semi-infinite beam on elastic foundation.

idealized beam on elastic foundation, E ¼ 125 GPa;Ei ¼ 5 GPa;B ¼ 30 mm;d ¼ 1:5 mm;t ¼ 0:05 mm yield b  920 m1 .5 3.2. PU enrichment In the partition of unity (PU) development, Babuska and Melenk [25] define a consistent approach for enhancing the basic isoparametric shape functions with the known solution, while guaranteeing the conformity of the approximation. Denoting by N i the standard shape functions, and assuming that each node i is enriched by J(i) shape functions, the proposed representation of the solution (displacements) u is:

uðxÞ ¼ Ni ðxÞðui þ uij /j ðxÞÞ;

i ¼ 1; M; j ¼ 1; JðiÞ

ð7Þ

where x and u are, respectively, the position and displacement vector; and ui ; uij are the standard and enriched displacement degrees of freedom of the element, M is the number of nodes of the element, J(i) is the number of functions enriching node i, and /j ðxÞ is the enrichment function. A typical PU enrichment strategy only involve those nodes affected by the enrichment function. This suggests, for a propagating crack, that (1) the crack front location should be recorded, and (2) the topology of some elements is modified during the course of the analysis. Such an option is not attractive with respect to implementation within most commercial packages, and it is preferable to impose the enrichment to all nodes affected by the potential crack path. This is realized with two user element subroutines (UEL): the first (continuum) element to model the core material and the second (cohesive) to model the interface. The development of these enriched formulations revealed a certain number of obstacles worthy of mention: 5 The range is between 900 m1 and 1600 m1 for a resin thickness between 0.1 mm and 0.01 mm, respectively.

1. An enriched linear, four-node element yields wrong predictions for thin beam bending, due to shear locking. The simple solution of reduced integration is, however, not applicable due to the presence of the higher-order enrichment. 2. Incompatible modes are strictly non-conforming, and may not be included along with the enrichment as proposed in Eq. (7), due to an inconsistency between the two enhancement strategies. 3. It does not seem possible to use, through the use of point-topoint constraints or surface tie constraints, an enriched solid element along with the standard linear cohesive element, and conversely, an enriched interface element with a generic continuum element. The discrepancy between the number of degrees of freedom of the nodes on the continuum element and the corresponding nodes on the cohesive elements is not permitted by the software used for this implementation, ABAQUS/Standard. Ongoing investigation is looking at overcome this limitation in order to explore the possibility of enriching only the interface element, and thus reduce the computational cost associated with the strategy proposed. 4. Some improvement in convergence is observed by using an enrichment function that vanishes at the opposite edge of the element; instead of using the analytical solution found above:

/ðÞ ¼ ebx ðcosðbxÞ  sinðbxÞÞ;

x¼

ðx2  x1 Þ ð þ 1Þ 2

ð8Þ

the enrichment function is modified to:

0

/ðÞ ¼ eb0 x ðcosðb0 x0 Þ  sinðb0 x0 ÞÞ;

x0 ¼  þ 1; b0 ¼

5p 8 ð9Þ

In Eq. (14), x1 and x2 are the respective coordinates of the first and second nodes of the enriched element. b0 is now strictly dependent on the element size, while b was a material property. Fig. 8 shows that the new enrichment function is only significantly different from the desired interpolation for elements of size smaller than 2 mm and larger than 6 mm – for the figure, b was set to 1000 m1 . In the current implementation, the parameter b remains an input to the model, and a cutoff element size is calculated. This allows to determine if the element size is too small, in which case no enrichment is utilized, i.e. a kinematic constraint is imposed on the enriched degrees of freedom. Because both the continuum and interface elements are to be enriched, an accurate representation of beam bending is as

2621

I. Guiamatsia et al. / Composites Science and Technology 69 (2009) 2616–2624 1

Normalized enrichment function

0.8

0.6

boef 6mm 5mm 3mm

0.4

2mm 0.5mm 0.2

0

-0.2

-0.4

0

1

2

3

4

5

6

Distance from crack tip (m)

x 10

-3

Fig. 8. Enrichment curves for element-size dependent b0 , versus material-based b (boef).

; g are the coordinates of the biunit square, the indices from 1 to 8 correspond to node numbers as per Fig. 9, N i are classic nodal shape functions for the eight-node serendipity element, and the enrichment functions LEð; gÞ and REð; gÞ are applied, respectively, to the left (nodes 1,4,8) and the right (nodes 2,3,6) edges of the element, are given as:

much an issue as the correct prediction of damage evolution. The solution adopted is based on quadratic shape functions in order to circumvent the shear locking problem. 3.3. Continuum element The enrichment is applied to nodes on both sides of the element (along the longitudinal direction x), a total of 22 degrees of freedom (dof), for the solid element, Fig. 9. Only the transverse displacement, v is enriched for mode I propagation. The detailed interpolation scheme can be written, in a two-dimensional reference frame (x, y), as:

LEð; gÞ ¼ MðgÞwðÞ; MðgÞ ¼ ð1  aÞg2 þ a; b0 xr

/ðÞ ¼ e

v ðx; yÞ ¼ N1 ð; gÞ½v 1 þ LEð; gÞv 1  þ N2 ð; gÞ½v 2 þ REð; gÞv 2  þ N3 ð; gÞ½v 3 þ REð; gÞv 3  þ N 4 ð; gÞ½v 4 þ LEð; gÞv 4  þ N5 ð; gÞv 5 þ N6 ð; gÞ½v 6 þ REð; gÞv 6  þ N7 ð; gÞv 7 þ N8 ð; gÞ½v 8 þ LEð; gÞv 8  ð11Þ

10 11 12

10 11 12

4 1 1 2 3

7 8 9

15 16

5 13 14

Enriched node

ðcosðb0 xr Þ  sinðb0 xr ÞÞ;

ð13Þ

xr ¼  þ 1

ð14Þ

7 8 9

18 19

4

20 21 22

3

6

ð12Þ xl ¼  þ 1

As shown in Fig. 10, the enrichment is the product of the analytical solution of a beam on an elastic foundation (active along the longitudinal axis), multiplied by some function with a variation in the transverse direction, MðgÞ (here, a ¼ 0:2). This was done to prevent rank deficiency that could arise from a uni-dimensional enrichment. The results presented here assume geometrical linearity. It is also worth noting that this element is directional in the formulation presented, with node 1 being assigned to a corner on the left side.

ð10Þ

zero thickness

a<1

wðÞ ¼ eb0 xl ðcosðb0 xl Þ  sinðb0 xl ÞÞ;

uðx; yÞ ¼ N 1 ð; gÞu1 þ N2 ð; gÞu2 þ N3 ð; gÞu3 þ N4 ð; gÞu4 þ N5 ð; gÞu5 þ N6 ð; gÞu6 þ N7 ð; gÞu7 þ N8 ð; gÞu8

REð; gÞ ¼ MðgÞ/ðÞ;

8

3 15 16 17

6

8

2 4 5 6

1 1 2 3

5 13 14

2 4 5 6

Fig. 9. Proposed elements with enrichment of the transverse displacement only. Left: cohesive element. Right: solid element.

2622

I. Guiamatsia et al. / Composites Science and Technology 69 (2009) 2616–2624

Fig. 10. Shape of the enrichment function for the 2D plane strain/stress element.

Table 1 Bending modelling capabilities of various ABAQUS elements. Number of elements through thickness

ABAQUS element

1

CPE4 – fully integrated, four-node linear CPE4I – incompatible mode, four-node linear CPE4R – reduced integration, four-node linear

4

P

CPE8 – fully integrated, eight-node quadratic UEL, eight-node quadratic, with enrichment

0.7967 0.1253 2.99E3 0.1263 0.1263

CPE4 – fully integrated, four-node linear CPE4I – incompatible mode, four-node linear CPE4R – reduced integration, four-node linear CPE8 – Fully integrated, eight-node quadratic UEL, eight-node quadratic, with enrichment

0.8705 0.1264 0.1203 0.1257 0.1258

As mentioned, the good quality of the solid element in modelling bending is crucial to the performance for the intended use for beam delamination. A beam of rectangular section, clamped at one end and loaded with an applied transverse displacement v at the other end, is considered in order to verify the bending behavior of the new enriched solid element. In this example, we have: Young’s modulus E = 150,000 MPa, beam thickness h = 0.5 mm, and width B = 0.2 mm. Ten elements are used in the longitudinal direction, and either 1 or 4 elements are used in the transverse direction. The analytical solution for the reaction force at the location of applied displacement is known:

P¼

3EIv L3

¼ 0:1172N

Table 1 compares the results obtained with the current element to other generic elements available in ABAQUS/Standard, and shows a good accuracy for the fully integrated enriched element, comparable to the incompatible mode linear element. 3.4. Cohesive element Analogous to the solid element, the enrichment is applied to the left and right side of a six-node quadratic cohesive element, resulting in a 16 dof interface element, as shown in Fig. 9.

If N i0 ðÞ; i ¼ 1; 3 are isoparametric shape functions for a onedimensional, three-node element, then the separation is formulated as follows:

DuðxÞ ¼ N10 ðÞðu4  u1 Þ þ N20 ðÞðu3  u2 Þ þ N30 ðÞðu6  u5 Þ ð15Þ Dv ðxÞ ¼ N 10 ðÞ½ðv 4  v 1 Þ þ wðÞðv 4  v 1 Þ þ N20 ðÞ½ðv 3  v 2 Þ þ /ðÞðv 3  v 2 Þ þ N 30 ðÞ½ðv 6  v 5 Þ ð16Þ The enrichment functions wðxÞ and /ðxÞ are given by Eqs. (13) and (14). An eight-point Gaussian quadrature is utilized in the cohesive element, while five Gauss points are sufficient for the continuum element. This high-order integration is necessary to benefit from the high-gradients enrichment function, especially with both damaged and undamaged zones set to appear within the same, larger element. Similar to the solid element, this cohesive element is two-way directional:6 Node 1 is assigned to the lower left corner of the element. 4. Results The double cantilever beam (DCB) test described in reference [7] is used for the validation of the current implementation. A mesh sensitivity study is performed, with element sizes of 1 mm, 2 mm and 5 mm. The results are reported for the quadratic formulation without enrichment (Fig. 11) and the proposed quadratic formulation with enrichment (Fig. 12). It is interesting to note that, although the results for the 5 mm mesh remain above the theoretical solution, the model accurately predicts the delamination of the beam over its entire length, for the displacement applied. Further improvement of the solution for the 5 mm mesh will be sought through techniques such as viscous regularization, not linked to the enrichment procedure itself. 4.1. Computational cost Cohesive element PU enrichment has been utilized before, Crisfield and Alfano [19] used cubic and quintic polynomials to 6 As opposed to the standard, non-enriched element that is only oriented with respect to the transverse direction.

2623

I. Guiamatsia et al. / Composites Science and Technology 69 (2009) 2616–2624

160 140

1mm

2mm

5mm

Reaction Force (N)

120 100 80 60 40 20 0 0

0.005

0.01

0.015

0.02

Separation (m) Fig. 11. DCB benchmark: quadratic interpolation, without enrichment.

120

Reaction Force (N)

100

80

60

40

1mm

20

2mm

5mm

0 0

0.002

0.004

0.006

0.008

0.01

0.012

0.014

0.016

0.018

0.02

Separation (m) Fig. 12. DCB benchmark: quadratic interpolation, with enrichment.

improve the non-linear solution convergence. The main advantage of the topology proposed here is that there is no need to keep record of neither the crack front or the direction of propagation, due to the symmetry of the enrichment functions of each element. However, there is a significant increase of the number of degrees of freedom per element, and the computational performance of the enriched formulation needs to be assessed. If we consider a 2 mm segment of the double cantilever beam, the total number of degrees of freedom 44 dof for a single element7 through the length, using the topology proposed here. In contrast, the segment only involves a total of 40 dof if linear, 0.5 mm cohesive elements are used for the interface and continuum incompatible mode elements for the beam halves. However, for the DCB benchmark considered throughout this paper, the processing cost measured as the wallclock time is 828 s for the 0.5 mm, linear model, compared to 524 s for the 2 mm proposed element. This seems to suggest that the presence of the enriching functions also improves the convergence of the non-linear solution.

7 One enriched cohesive element sandwiched between two enriched continuum elements.

5. Concluding remarks This paper addressed the mesh size issue in the simulation of delamination with the CZM approach by proposing two enriched elements: a continuum solid element for the plies, and a cohesive element for the interface. The enrichment function is the analytical solution of a semi-infinite beam on elastic foundation, which mimics the high gradients of stress field ahead of the crack tip. The enhanced elements have a quadratic formulation, as a baseline; they are formulated as partition of unity (PU) finite elements and implement as user elements (UEL), compatible with the finite element software ABAQUS/Standard. Several issues were noted during this development, including (1) the inability to use a linear interpolation element as a baseline for enrichment, as well as (2) the need to enforce a two-dimensional variation of the enrichment functions in the solid element formulation. Excellent results were obtained for element sizes of up to 5 mm, which represents a tenfold improvement with respect to the current situation. On the other hand, given that the average number of degrees of freedom per node is 2.75 for the current formulation, as opposed to 2 for the standard linear interpolation, meshes of at least 2 mm in

2624

I. Guiamatsia et al. / Composites Science and Technology 69 (2009) 2616–2624

size are needed in order to make up for the additional number of degrees of freedom incurred with the enriched formulation. However, it should be mentioned that the computational cost is not solely measured by the number of degrees of freedom: for the cases presented here, the average solution time for the enriched formulation with elements 2 mm in size was significantly lower than the standard linear formulation with elements of size 0.5 mm. This indicates that the enrichment also equipped the element with better convergence properties. Further development of this work will extend this formulation for three-dimensional models, and look at implementation with an explicit solver. This should enable greater use of cohesive interface elements within large scale models, as this have been prevented by the small mesh size requirement. Acknowledgements This work is part of a project looking at the development of efficient numerical approaches for the analysis of damage in composite structures with geometrical discontinuities. The joint support of EPSRC and DSTL through grant EP/E023967/1, is gratefully acknowledged. References [1] Rybicki EF, Kanninen EF. A finite element calculation of stress intensity factors by a modified crack closure integral. Eng Fract Mech 1977;9:931–8. [2] Krueger R. Virtual crack closure technique: history approach and applications. Appl Mech Rev 2004;57(2):109–44. [3] Mabson GE, Deobald LR, Dopker B, Hoyt DM, Baylor JS. Fracture interface elements for static and fatigue analysis. In: 16th International conference on composite materials, 2007, Kyoto, Japan. [4] Hillerborg A, Modeer M, Petersson P-E. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 1976;6:773–82. [5] Schellekens JCJ, De Borst R. A non-linear finite-element approach for the analysis of mode I free-edge delamination in composites. Int J Solids Struct 1993;30(9):1239–53. [6] Schellekens JCJ, De Borst R. On the numerical integration of interface elements. Int J Numer Meth Eng 1993;36:43–66.

[7] Davies GAO. Benchmarks for composite delamination. In: Benchmarks for composites. NAFEMS, Imperial College London; 2002. [8] Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids 1960;8:100–4. [9] Barenblatt GI. The mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech 1962;7:55–129. [10] DS Simulia, ABAQUS/STANDARD Analysis user’s manual, v.6.8; 2008. [11] Robinson P, Besant T, Hitchings D. Delamination growth prediction using a finite element approach. In: 2nd ESIS TC4 conference on polymers and composites, 1999, Les Diablerets, Switzerland. [12] Allix O, Ladeveze P. Interlaminar interface modelling for the prediction of delamination. Compos Struct 1992;22:235–42. [13] Allix O, Ladeveze P, Corigliano A. Damage analysis of interlaminar fracture specimens. Compos Struct 1995;31:61–74. [14] Camanho PP, Davila CG, Moura M. Numerical simulation of mixed-mode progressive delamination in composite materials. J Compos Mater 2003;37:1415–38. [15] Turon A, Camanho PP, Costa J, Davila CG. A damage model for the simulation of delamination in advanced composites under variable loading. Mech Mater 2006;38:1072–89. [16] Turon A, Davila 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–82. [17] Harper PW, Hallett SR. Cohesive zone length in numerical simulations of composite delamination. Eng Fract Mech 2008;75(16):4774–92. [18] Crisfield MA. A fast incremental/iterative solution procedure that handles ‘snap-through’. Comput Struct 1981;13:55–62. [19] Crisfield MA, Alfano G. Adaptive hierarchical enrichment for delamination fracture using a decohesive zone model. Int J Numer Meth Eng 2002;54:1369–90. [20] Mi Y, Crisfield MA, Davies GAO, Hellweg HB. Progressive delamination using interface elements. J Compos Mater 1997;32:1246–72. [21] Benzeggagh MK, Kenane M. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/eposy with mixed-mode bending apparatus. Compos Sci Technol 1996;56:439–49. [22] Jiang WG, Hallett SR, Green BG, Wisnom MR. A concise interface constitutive law dor analysis of delamination and splitting in composite materials and its application to scaled notched tensile specimens. Int J Numer Meth Eng 2006;69:1982–95. [23] Wisnom MR, Petrossian ZJ, Jones MI. Interlaminar failure of unidirectional glass/epoxy due to combined through thickness shear and tension. Composites Part A: Appl Sci Manuf 1996;27A:921–9. [24] Alfano G, Crisfield MA. Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues. Int J Numer Meth Eng 2001;50:1701–36. [25] Melenk JM, Babuska I. The partition of unity finite element method: basic theory and applications. Austin: The University of Texas at Austin; 1996. 32p.