Micromechanics-based progressive failure analysis of fibre-reinforced composites with non-iterative element-failure method

Computers and Structures 89 (2011) 1103–1116

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Micromechanics-based progressive failure analysis of fibre-reinforced composites with non-iterative element-failure method X.S. Sun ⇑, V.B.C. Tan, T.E. Tay Department of Mechanical Engineering, National University of Singapore, 9 Engineering Drive 1, Singapore 117576, Singapore

a r t i c l e

i n f o

Article history: Received 8 June 2010 Accepted 2 December 2010 Available online 6 January 2011 Keywords: Fibre-reinforced composite Progressive failure analysis Micromechanics of failure (MMF) criterion Non-iterative element-failure method (NEFM)

a b s t r a c t This paper proposes a micromechanics-based progressive failure analysis strategy for fibre-reinforced composites. A non-iterative element-failure method (NEFM) is developed and integrated with a micromechanics-based failure criterion. The micromechanics of failure (MMF) is further investigated and modified to improve predictions. The NEFM is focused on implementation of damaged elements with direct solution and implicit nodal forces. Progressive failure of fibre-reinforced laminates is analysed with the micromechanics-based NEFM and the results are compared with those from the traditional finite element method or experiments. The comparison demonstrates computational efficiency and validity of the micromechanics-based NEFM in progressive failure analysis of composites. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Fibre-reinforced composite materials are being extensively used in modern engineering practices due to their high strength to weight ratio compared to traditional metal materials. Usually the failure of composite materials is much more complicated than that of isotropic and homogeneous materials. Many researches on failure models and numerical methods for progressive failure analysis of composites have been conducted and published in literatures during the past decades [1–5]. As the failure of composites presents significant anisotropy and heterogeneity, various failure theories or criteria were developed for composite materials. Nineteen failure theories of the period were extensively assessed in a ‘‘World-Wide Failure Exercise (WWFE)’’ [6–9], where each participant theory was described in detail and then the corresponding prediction was compared with experimental results provided by the initiators of the WWFE. An overall comparison of all the participant failure theories against experimental results and an assessment of their predictive capabilities were presented with their validities and shortfalls. Among various failure theories, micromechanics-based failure criteria show more predictive capabilities because failure is determined at the constituent level in this kind of failure criteria. Typical examples of the micromechanics-based failure theories are the strain invariant failure theory (SIFT) by Gosse et al. [10,11] and the multicontinuum theory (MCT) by Garnich and Hansen [12]. The ⇑ Corresponding author. Tel.: +65 65162237; fax: +65 67791459. E-mail address: [email protected] (X.S. Sun). 0045-7949/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2010.12.003

strain-based SIFT criterion uses two strain invariants, i.e. the critical volumetric strain (the first strain invariant) and the critical equivalent strain (the von Mises strain which is a function of the first and second strain invariants), to identify failures in fibre or matrix. However, the critical values used in the SIFT criterion need to be determined by special tests on laminates with designed lay-ups. The MCT criterion is based on a continuum definition of stress at a point, where stress and strain of constituents (fibre and matrix) are averaged, and fibre and matrix failure modes are determined by separated conditions, respectively. Gotsis et al. [13] and Huang [14] also proposed other micromechanicsbase failure criteria, respectively. Recently, Ha et al. [15] developed a new micromechanics of failure (MMF) criterion for fibrereinforced composites, where fibre and matrix failures were determined by a maximum stress criterion and a modified von Mises failure criterion, respectively. Generally, micromechanicsbased failure criteria require some interaction information between the micro- and macro-levels, which should be obtained prior to failure analysis. In most cases of failure analysis for composites, the finite element method (FEM) is popularly used in the numerical analysis of failure progression for composites and numerous simulation methods integrating with failure models were developed [4,5]. Besides the traditional FEM, the finite strip method [16–18] and the element-failure method (EFM) [11,19,20] were also developed for progressive failure analysis of composites. In recent years, some other numerical methods, typically including the boundary element method (BEM) [21–24], the extended finite element method (X-FEM) [25,26], the enriched element-failure method (REFM) [27],

1104

X.S. Sun et al. / Computers and Structures 89 (2011) 1103–1116

the element-free or meshless method [28,29], among others, were applied to some special failure problems such as crack propagation or delamination of fibre-reinforced composites. These methods, however, are not suitable for general progressive failure analysis of composite structures with complex geometries or various layup configurations. Progressive failure analysis of composites is usually complicated and time-consuming. Conventional FEM-based numerical models for the progressive failure analysis of composites often adopts a scheme of material property degradation method (MPDM), which generally requires frequent assembly and decomposition of global stiffness matrix in each load increment in the progressive analysis. However, few elements may fail in some load increments (especially during the initial phase of failure), and the frequent assembly and decomposition of global stiffness matrix accordingly would reduce the computational efficiency. Tay et al. [11,19,20] proposed an element-failure method (EFM) for progressive failure analysis of composites based on the framework of FEM, where failed elements were explicitly applied external nodal forces determined by equilibrium of adjacent elements. In the EFM, solutions are obtained by equilibrium iteration of applied external nodal forces and the global stiffness matrix is kept constant throughout the analysis. Although the EFM does not require reformulation or decomposition of global stiffness matrix in each load increment, large numbers of equilibrium iterations usually have to be carried out to obtain the solutions, which would also reduce the computational efficiency. This paper develops a non-iterative element-failure method (NEFM) for failure analysis of composites based on the framework of FEM, and a micromechanics-based analysis strategy integrating NEFM with MMF is proposed for progressive failure analysis of fibre-reinforced composites. The MMF criterion is further investigated to improve the problems of premature matrix failure and under-prediction of shear strength, and the micromechanics analysis based on a representative volume element (RVE) or unit cell is elaborated. The basic equations in the NEFM are established from traditional FEM and the implementation of the micromechanics-based NEFM is presented. Examples are used to demonstrate applications of the proposed micromechanics-based NEFM to progressive failure analysis of fibre-reinforced composites, and the results from the proposed method are also compared with those from the traditional FEM (with MPDM) or experiments. 2. Micromechanics-based failure criterion for composites The micromechanics of failure (MMF) criterion was developed very recently for failure prediction of fibre-reinforced composites, which determines the failure initiation at the micromechanics level in the constituents of fibre and matrix [15]. In the MMF criterion, the failure at the macro level (ply level) is indicated by the failures of fibre and matrix at the micro level (constituent level). 2.1. Failure conditions and stress amplification factors Usually, two failure conditions based on micro stresses, i.e. the maximum stress failure condition and the modified von Mises failure condition, are used to determine the fibre failure and the matrix failure, respectively. These failure conditions constitute the MMF criterion and can be written as [15]: (1) Fibre failure condition:

r2f 1 T f Cf

 þ

1 1  T f Cf



rf 1 ¼ 1:

ð1Þ

(2) Matrix (inter-fibre) failure condition:



r2VM T mCm

þ

 1 1 I1 ¼ 1;  T m Cm

ð2Þ

where Tf, Cf, Tm and Cm are the tensile and compressive strengths of fibre and matrix, respectively, which are backcalculated from the ply strengths; rf1 is the micro longitudinal stress in the fibre; rVM and I1 are the micro von Mises stress and first stress invariant in the matrix; and

8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 > > < rVM ¼ I1  3I2 ; > I1 ¼ r1 þ r2 þ r3 ; > : I2 ¼ r1 r2 þ r2 r3 þ r3 r1  ðs212 þ s223 þ s213 Þ:

ð3Þ

Generally, the micro stress (constituent-level stress) r in both fibre and matrix can be calculated with the macro stress (ply-level  through some influence functions [30,31]. With a represtress) r sentative micromechanical block (or unit cell model), r can be computed through the mechanical and thermo-mechanical stress amplification factors, Mr and Ar [32]:

r ¼ M r r þ Ar DT;

ð4aÞ

where DT is the temperature increment, and

(

r ¼ ½r1 r2 r3 s12 s23 s13 T ; r ¼ ½r 1 r 2 r 3 s12 s23 s13 T ; 2

ð5Þ

M 12

M 13

0

M 15

6M 6 21 6 6 M 31 Mr ¼ 6 6 0 6 6 4 M 51

M 22

M 23

0

M 25

M 32

M 33

0

M 35

0

0

M 44

0

M 52

M 53

0

M 55

0 7 7 7 0 7 7 ; M46 7 7 7 0 5

0

0

M 64

0

M66 r

0 Ar ¼ ½ A1

A2

A3

0 A5

0 Tr :

0

3

M 11

ð6aÞ

ð7Þ

The stress amplification factors Mr and Ar in Eqs. (6a) and (7), which account for the macro and micro interactions, can be extracted by finite element analyses with representative micromechanical blocks (or unit cell models). In the application of the MMF criterion, some reference locations within the fibre and matrix materials of the unit cell models are selected, where the stress amplification factors are extracted and then used to calculate the micro stress. The typical unit cell models, the reference points in the models and the extraction of stress amplification factors are given in Appendix A. There are however two common problems in the abovementioned MMF criterion: (1) premature matrix failure under longitudinal tension (smaller matrix failure envelope), which typically occurs for glass/epoxy composites and (2) lower predicted shear strength. This is usually because the matrix under tension and the ply under shear often present significantly nonlinearity, while linear elastic stress–strain relationships are used in the micromechanics analysis. Here, a simplified approximation method is proposed to modify the MMF criterion and circumvent the abovementioned problems. In the modified MMF criterion, the additional interface failure condition for glass/epoxy composites used in the original MMF [15] is treated as a kind of inter-fibre (matrix) failure, and the failures of both carbon/epoxy and glass/epoxy composites are unified and determined by same conditions as defined in Eqs. (1) and (2). The basic idea of the approximation method is based on modifying the matrix mechanical stress amplification factors in order to compensate the nonlinear influence of the material

1105

X.S. Sun et al. / Computers and Structures 89 (2011) 1103–1116

properties. The matrix stress amplification factor Mr can be modified with some correction coefficients as the following:

2

c1 M11

6 M 21 6 6 6 M 31 0 Mr ¼ 6 6 0 6 6 4 M 51 0

Fig. 1. Flowchart for backcalculating Tm and Cm.

0

3

M12

M13

0

M 15

M22

M23

0

M 25

M32

M33

0

M 35

0 M52

0 M53

c4 M 44 0

0 M 55

7 7 7 0 7 7 ; c6 M46 7 7 7 0 5

0

0

c4 M 64

0

c6 M66 r

0

ð6bÞ

where c1, ck (k = 4, 6) are correction coefficients determined by the matrix tensile failure strain and the ply in-plane shear strength (see Appendix B), respectively. The failure in the fibre-reinforced composites is generally caused by mechanical loads rather than temperature. The micro thermal stresses (in fibre or matrix) caused by temperature are usually not very high and far less than the constituent strengths, which means the thermal-mechanical stress amplification factor Ar has only minor contribution to the total micro stresses and the major micro stresses are due to the mechanical stress amplification factor Mr. Accordingly, only Mr is modified in order to avoid the unnecessary complexity and low efficiency in the computation.

1  r  2 ; (b) Fig. 2. Failure envelopes of the modified MMF criterion in biaxial macro stress spaces (material: UD laminate of carbon/epoxy IM7/5250-4): (a) stress space r 1  s 12 and (c) stress space r 2  s 12 . stress space r

1106

X.S. Sun et al. / Computers and Structures 89 (2011) 1103–1116

This modification will not change the backcalculated values of Tm and Cm because only Mrj2 (j = 1, 2, . . . , 6) is used in the backcalculation (see Section 2.2). With the modified matrix stress amplification factors, the micro stress in the matrix can be rewritten as:

r ¼ M 0r r þ Ar DT:

ð4bÞ

2.2. Backcalculation of the constituent strengths In the MMF criterion, the constituent strengths including Tf, Cf, Tm and Cm can be backcalculated from the ply strengths with the unit cell modes and reference points. The constituent strengths of fibre, Tf and Cf, are backcalculated from the ply longitudinal tensile strength XT and compressive strength XC, respectively, which is straightforward and easy to obtain from Eqs. (1) and (4a) as:

8 < T f ¼ maxðM ðiÞ X T þ AðiÞ DTÞ; r11;f r1;f i ¼ 1; 2; 3 . . . : C f ¼ maxðj  MðiÞ X C þ AðiÞ DTjÞ; r11;f

ð8Þ

r1;f

where Mr11,f and Ar1,f are the fibre stress amplification factors, and i = 1, 2, 3. . . denotes reference points in the fibre of the unit cell model (see Appendix A). The constituent strengths of matrix, Tm and Cm, are backcaculated from the ply transverse tensile strength YT and compressive strength YC. This backcalculation, however, is not straightforward and achieved by iterative computation. Define the following two parameters:



bm ¼ C m =T m ;

ð9Þ

p ¼ req =T m ;

where req is an equivalent stress defined with Eq. (2) as the following:

req ¼

ðbm  1ÞI1 þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðbm  1Þ2 I21 þ 4bm r2VM 2bm

:

ð10Þ

Fig. 3. Flowchart for implementing the micromechanics-based NEFM.

With an initial value bm,0, e.g. bm,0 = 1, a uniaxial macro stress with the value of YT and YC, respectively, is applied to the unit cell models and the equivalent stress req at each reference point can be computed with Eqs. (3), (4a), and (10). Then two values of Tm corresponding to the two loading cases can be found and determined equal or not. If they are equal, the iteration terminates and the values of Tm and Cm are obtained; otherwise, the value of bm is changed for the next iteration. The flowchart of the iterative computation is illustrated in Fig. 1. It can be seen that the backcalculated constituent strengths in the MMF criterion usually reflect the influence of model geometry, fibre volume fraction, stress amplification factors and temperature. 2.3. Failure envelopes in macro stress spaces As mentioned above, failure prediction by the MMF criterion is based on micro stress at the constituent level. The interaction between the macro stress at the ply level and the micro stress at the constituent level is represented by the stress amplification factors. Therefore, failure envelopes in the macro stress spaces can be plotted with stress interaction between micro- and macro-levels. Fig. 2 shows the failure envelopes of the modified MMF criterion in typical biaxial macro stress spaces, where the original MMF criterion and the popular Tsai-Wu criterion are also used to compare with the modified MMF criterion. For carbon/epoxy composites, there is usually no much difference between the failure envelopes of the modified and original MMF criteria in the macro stress space r 1  r 2 , but the difference in shear-related stress spaces is much significant as shown in Fig. 2(b) and (c). The shear strength is under-predicted remarkably by the original MMF criterion, whereas

Fig. 4. Flowchart for micromechanics analysis based on unit cell models.

the under-prediction is well improved in the modified MMF criterion.

X.S. Sun et al. / Computers and Structures 89 (2011) 1103–1116

ðK 0  DKÞaI ¼ P I ;

3. Non-iterative element-failure method (NEFM) The non-iterative element-failure method (NEFM) is focused on implementation of damaged elements with direct solution and implicit external nodal forces applied to the damaged elements. Compared with the traditional FEM, NEFM does not require frequent reformulation and decomposition of global stiffness matrix in each load increment.

ð11Þ

where K0 is the initial global stiffness matrix without damage; DK is the stiffness matrix for damaged elements, which is updated with failure progression; aI and PI are the displacements and known loads in the Ith load increment. By left multiplying K 1 0 in Eq. (11), the basic equation in NEFM can be derived from Eq. (11) as:

~ I ¼ F~ I ; Ka 3.1. General equations in the NEFM

1107

ð12aÞ

where

For the progressive failure analysis under the framework of FEM, the discretized system equation in one load increment can be formulated as the following:

(

~ ¼ I  K 1 DK; K 0 ~ F I ¼ K 1 P I :

ð13Þ

0

In Eq. (13), only one-time LU decomposition (e.g. Cholesky factorization, K0 = LLT) in the initial analysis is required and then stored for use in subsequent computations. The right-hand term

Table 1 Material properties of the carbon/epoxy (IM7/5250-4) UD laminate. Material properties

Ply (IM7/ 5250-4)

Fibre (IM7)

Matrix (5250-4)

Fibre volume fraction Vf Longitudinal modulus E1 (GPa) Transverse moduli E2 = E3 (GPa) Shear moduli G12 = G13 (GPa) Shear modulus G23 (GPa) Poisson’s ratios v12 = v13 Poisson’s ratio v23 Longitudinal coefficient of thermal expansion a1 (/°C) Transverse coefficients of thermal expansion a2 = a3 (/°C)

0.62 172.4 10.3 5.52 3.45 0.32 0.40 0.11  106

1.0 276.0 27.6 138.0 7.8 0.30 0.80 0.036  106

0.0 3.45 3.45 1.28 1.28 0.35 0.35 46.8  106

7.610-6

5.04  106

46.8  106

Table 2 Uniaxial strengths of the carbon/epoxy (IM7/5250-4) UD laminate. Ply strengths

Values

Longitudinal tensile strength XT (MPa) Longitudinal compressive strength XC (MPa) Transverse tensile strength YT (MPa) Transverse compressive strength YC (MPa) Shear strengths S12 = S13 (MPa) Shear strength S23 (MPa)

2826.5 1620.0 65.5 248.0 122.0 122.0

Fig. 6. Comparison of computational efficiency.

Fig. 5. Example of crack propagation problem: (a) model and (b) mesh.

1108

X.S. Sun et al. / Computers and Structures 89 (2011) 1103–1116

Fig. 7. Comparison of computational results: (a) normal FEM; (b) MPDM and (c) NEFM.

Fig. 8. Example of open-hole-tension problem: (a) model; (b) Mesh 1 (400 elements per ply); (c) Mesh 2 (784 elements per ply) and (d) Mesh 3 (1296 elements per ply).

F~ I in Eq. (12a) can be easily obtained by forward and backward substitutions with the LU decomposition. Especially, under a proportional loading case which is common for many problems, i.e.

PI = cIP0, where P0 is the initial load and cI is a scalar parameter denoting loading factor in the Ith load increment, F~ I can be directly calculated by F~ I ¼ cI F~ 0 with F~ 0 ¼ K 1 0 P 0 , and only one-time

X.S. Sun et al. / Computers and Structures 89 (2011) 1103–1116

1109

forward and backward substitutions in the initial analysis are required. Additionally, K 1 0 is also not explicitly calculated and only columns related to damaged elements are required, which are obtained by forward and backward substitutions of the LU decomposition. For example, if Cholesky factorization is used in the LU decomposition of global stiffness matrix, i.e. K0 = LLT, the non-zero jth column related to damaged elements in K 1 0 can be calculated by the following equation: T 1 ½K 1 0 j ¼ L L I j ;

ð14Þ

where Ij is a unit array with zero rows except for the jth row, i.e.:

 ½Ii j ¼

1; if i ¼ j;

ð15Þ

0; if i–j:

As K0 is constant through the analysis, only columns in K 1 0 related to new damaged elements are calculated in each load increment and then stored for use in subsequent load increments. The multiplication K 1 0 DK in Eq. (13) accordingly takes the form as:

2

0 6  6  6 1 K 0 DK ¼ 6 6 0 6 4

   k1j



0

3

7    7 7    kij    0 7 ; 7 7    5    knj    0 ðnnÞ

0

ð16Þ

where only columns related to damaged elements are non-zero. ~ and F~ I are obtained, Eq. (12a) can be rearranged with After K the following formulation:

"

~d K ~u K

0

#

1

ad au

"

 ¼ I

F~ d F~ u

# ;

ð12bÞ

I

~ d and K ~ u are condensed matrices related to damaged and where K undamaged elements, respectively; the subscripts ‘‘d’’ and ‘‘u’’ denote damaged elements and undamaged elements, respectively, and

(

aI ¼ ½ ad F~ I ¼ ½ F~ d

au TI ; T F~ u I :

ð17Þ

Then the large-size problem described by Eqs. (12a) or (12b) can be reduced to a small-size problem described by the following equation:

~ d ad ¼ F~ d ; K

ð18Þ

where the solution is only focused on damaged elements and the problem size is dependent on the number of damaged elements. After the displacement related to damaged elements, ad, is obtained, the displacement related to undamaged elements, au, can be obtained straightforward with substituting ad into the following equation according to Eq. (12b):

~ u ad : au ¼ F~ u  K

ð19Þ

It can be seen that no iteration is required in solving Eq. (18), which is distinguished from the iterative EFM. The NEFM is also different with the general MPDM in treating damaged elements: the material properties of damaged elements are not degraded in the NEFM, but external nodal forces are implicitly applied to damaged elements and the damage is modeled through element failure caused by the nodal forces. The external nodal forces are implicit because they are not involved into solving system equations as in the iterative EFM. The implicit external nodal forces in the Ith load increment are determined by the displacement solution as the following:

F N;I ¼ K d ad;I ;

Fig. 9. Comparsion of results predicted by the MMF–NEFM with different meshes: (a) stress–strain curve and (b) ultimate strength.

ð20Þ

where Kd is the condensed stiffness matrix related to damaged elements. 3.2. Implementation of the micromechanics-based NEFM In the NEFM, the global stiffness matrix for damaged elements, DK, as indicated in Eq. (13), is assembled with the element stiffness matrix for damaged elements, DKe, which is determined with the following equations by the MMF criterion: ( e K ; if fibre fails; ð21Þ DK e ¼ K em ; if matrix fails; where

(

P

 BT DB; P  m B; K em ¼ BT D Ke ¼

ð22Þ

 and D  m are full and reduced global ply stiffnesses, respecand D tively, which are obtained with transformation of full and reduced local ply stiffnesses, D and Dm, respectively, and

1110

X.S. Sun et al. / Computers and Structures 89 (2011) 1103–1116

Fig. 10. Failure patterns predicted by the MMF–NEFM with different meshes (Red elements: fibre failure; Blue elements: matrix failure; White elements: no damage): (a) Mesh 1; (b) Mesh 2 and (c) Mesh 3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

 ðDm Þij ¼

0;

if i; j ¼ 1;

Dij ; if i; j–1:

ð23Þ

In order to keep the number of accumulated damaged elements not too large and speed up the computation of NEFM, the initial stiffness matrix without damage in Eq. (13), K0, can be replaced by a reassembled stiffness matrix after every J0 increments, K J0 , and DK is also replaced by that after every J0 increments, DK J0 , which is assembled only for new damaged elements in every J0 increments. Then the NEFM equation can be modified as:

 ~ J aI ¼ F~ J ; K 0 0

ð24Þ

I

and the reassembled stiffness matrix, K J0 , is determined by the following equation:

K J0 ¼

mX t md

Ke þ

md X

K ed ;

ð25Þ

where mt and md are total number of elements and number of previous accumulated damaged elements, respectively, and

( K ed

¼

K 0e ;

An in-house code is developed to implement progressive failure analysis by the micromechanics-based NEFM. The flowchart for implementing the micromechanics-based NEFM is illustrated in Fig. 3. A micromechanics analysis based on unit cell models is also carried out with the traditional FEM prior to the macro analysis in order to obtain the stress amplification factors and the backcalculated constituent strengths used in the MMF criterion, as shown in Fig. 4. Details of the micromechanics analysis can be referred to Appendix A. The data from the micromechanics analysis can be stored to use for other failure analyses with the same materials.

if fibre fails;

K 0em ; if matrix fails;

where dM  0 is a material degradation factor.

In this section, carbon/epoxy unidirectional (UD) laminate and multidirectional laminate [45/0/–45/90]s (quasi-isotropic laminate) are used to demonstrate the application of the proposed micromechanics-based NEFM in progressive failure analysis. The composite material properties and uniaxial ply strengths used in the demonstrations are listed in Tables 1 and 2, respectively.

ð26Þ

and the element stiffness matrix K 0e and K 0em are determined by the local ply stiffnesses D0ij and ðD0m Þij with appropriate transformations, which are defined as:

D0ij ¼ dM Dij  0;  if i; j ¼ 1; Dij ; ðD0m Þij ¼ dM Dij ; if i; j–1;

4. Examples of progressive failure analysis with the NEFM

ð27Þ ð28Þ

4.1. Inter-fibre crack propagation of UD laminates under transverse tension An example of inter-fibre crack propagation of a UD laminate under transverse tension is used to demonstrate the computational efficiency of the NEFM, which is also compared with that of MPDM. The model and mesh are shown in Fig. 5, where prescribed uniform displacement is applied along the transverse direction and the crack is modeled with failed elements (gray elements). Total

X.S. Sun et al. / Computers and Structures 89 (2011) 1103–1116

3690 3D linear brick elements with 20500 degree of freedoms (DOFs) are used in the analysis. The inter-fibre crack propagates along the fibre direction, which is dominated by matrix (inter-fibre) failure and causes a splitting failure of the UD laminate. The splitting failure progresses very quickly in an actual loading test. Nevertheless, the implementation of the micromechanics-based NEFM can be conceptually demonstrated with this example. Fig. 6 shows the comparison of computational efficiency of the NEFM and MPDM, with normalized computation time by that of MPDM, where NEFM_1 and NEFM_2 denote results from reassembling and decomposing global stiffness matrix every 5 and 10 load increments, respectively, and 3 elements fail in each load increment. It can be seen that the efficiency of NEFM is improved compared to that of MPDM. The transverse stress contour of the UD laminate after 10 load increments are illustrated in Fig. 7, and the result from the current NEFM is also compared with those from normal FEM and MPDM. The comparison shows similar results for the three methods but with different schemes on treating the failure. With failure progression, failed elements are deleted in the normal FEM model and material properties (stiffness) of failed elements are degraded to near zero in the MPDM model. In the NEFM model, however, failed elements are neither deleted nor degraded, but are implicitly applied with some nodal forces, which are indicated by large stress distribution in the failed elements as shown in Fig. 7(c).

1111

in the 90o ply. The three mesh sizes demonstrate similar failure patterns in each ply. The results in Figs. 9 and 10 show that the MMF–NEFM is not very sensitive to element size in the open-hole tension problem if the mesh is reasonable. Therefore, the Mesh 1 is used to demonstrate the comparison of the MMF–NEFM with other methods in the followings. The stress–strain curves and ultimate strengths predicted by different methods with Mesh 1 are shown in Fig. 11(a) and (b), respectively. As the material properties of damaged elements are not fully reduced to zero in the MPDM, both the residual stiffness after fibre failure and ultimate strength predicted by the MPDM are a little higher than those by the NEFM. Fig. 12 shows failure patterns in each ply predicted by different methods after fibre failure. All the four methods produce similar and reasonable failure patterns in each ply, which indicates that the proposed MMF–NEFM is applicable to progressive failure analysis of composites. The damage predicted by the MMF–NEFM is also compared with the X-ray image of damage from experiment, as shown in Fig. 13. It can be seen that the MMF–NEFM shows reasonable prediction of damage, which also demonstrates the validity of the MMF–NEFM in progressive failure analysis of fibre-reinforced composites.

4.2. Progressive failure of open-hole laminates under tension Progressive failure of a quasi-isotropic laminate [45/0/–45/90]s with open hole under tension is analysed with the proposed micromechanics-based NEFM. The dimension size of the laminate is L  W  H = 76.2  76.2  1 mm (each ply thickness is 0.125 mm), and the radius of the centered hole is r = 9.525 mm. Uniform displacement is applied to one side of the laminate. Fig. 8 shows the model and mesh of the open-hole laminate under tension. Three mesh sizes including meshes with 400, 784 and 1296 3D linear brick elements per ply, respectively, as shown in Fig. 8(b)–(d), are used to investigate the potential mesh sensitivity in the micromechanics-based NEFM (MMF–NEFM). Four analysis methods combining the MPDM and NEFM with the Tsai-Wu and MMF criteria, respectively, i.e. Tsai–MPDM, Tsai–NEFM, MMF– MPDM and MMF–NEFM, are employed to compare the numerical results. In the MPDM, a non-zero degradation factor d = 1.0  103 is used to reduce the material properties of damaged fibre or matrix. In the MMF criterion, a curing temperature DT = 155.6 °C is also considered to capture residual stresses in the constituents (fibre and matrix). The averaged stress and strain on the boundary with the prescribed displacement loads are traced to record the characteristic stress–strain curve. The stress–strain curves and ultimate strengths predicted by the MMF–NEFM with different meshes are shown in Fig. 9(a) and (b), respectively. It can be seen that the stress–strain curves before fibre failure predicted by the MMF– NEFM with different meshes change little with reduced element size. For a coarser mesh with larger element size, there is a fluctuation in the stress–strain curve after fibre failure. This is usually because the failure influence of larger elements on the stress redistribution of the structure is greater than that of smaller elements. The ultimate strengths predicted by the MMF–NEFM with three mesh sizes are similar and show good agreement with the experiment. Fig. 10 shows the failure patterns in each ply predicted by the MMF–NEFM with different meshes after fibre failure. Fibre failure will occur first in the 0o ply, and there are significant matrix failures in the 45o, 45o and 90o plies. Usually, there are few matrix failures in the 0o ply, whereas extensive matrix failures

Fig. 11. Comparsion of results predicted by different methods with Mesh 1: (a) stress–strain curve and (b) ultimate strength.

1112

X.S. Sun et al. / Computers and Structures 89 (2011) 1103–1116

Fig. 12. Failure patterns predicted by different methods with Mesh 1 (Red elements: fibre failure; Blue elements: matrix failure; White elements: no damage): (a) TsaiMPDM; (b) Tsai-NEFM; (c) MMF–MPDM and (d) MMF–NEFM. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

a

b

Fig. 13. Comparison of damage: (a) X-ray image of damage [33] and (b) predicted damage by the MMF–NEFM (with Mesh 1 and superposition of all plies).

X.S. Sun et al. / Computers and Structures 89 (2011) 1103–1116

5. Conclusions A micromechanics-based non-iterative element-failure method (NEFM), which is integrated with the micromechanics of failure (MMF) criterion, is developed for progressive failure analysis of fibre-reinforced composites. The MMF criterion is further investigated and modified to circumvent the shortcomings of premature matrix failure and under-prediction of shear strength in the original MMF criterion. The improvement of the modified MMF criterion is demonstrated by failure envelopes of fibre-reinforced composites in typical biaxial macro stress spaces. The micromechanics analysis is based on two typical representative volume elements (RVEs) or unit cell models: square and hexagonal models. The stress amplification factors, which bridge the macro and micro analyses, are extracted from the two unit cell models, and the constituent strengths of fibre and matrix used in the MMF criterion are backcalculated from the ply strengths through the micromechanics analysis. The NEFM is carried out under the framework of FEM and focused on implementation of damaged elements. The system equation is solved directly without iteration, and external nodal forces are implicitly applied to the damaged elements. The damage is accordingly modeled through element failure caused by the implicit external nodal forces. Compared with the material property degradation method (MPDM), the NEFM does not require frequent reformulation and decomposition of global stiffness matrix in each load increment in the progressive failure analysis. Inter-fibre crack propagation of UD laminates under transverse tension and progressive failure of open-hole laminates under tension are used to demonstrate the application of the NEFM, and the results are compared with those from the MPDM and experiments. The NEFM shows better computational efficiency compared to the MPDM when the number of failed elements is small in large finite element models. The ultimate strength and failure patterns of open-hole tension problem predicted by the micromechanics-based NEFM

1113

(MMF–NEFM) are in good agreement with those from experiments. Different mesh sizes and methods are also used to compare the results from the MMF–NEFM. The comparison shows the validity and reliability of the MMF–NEFM in progressive failure analysis of fibre-reinforced composites. Although only solid elements are used in the demonstrations, it is straightforward to extend the application of the NEFM with other element types such as plate element or shell element which are widely used in analysis of fibrereinforced laminates. Acknowledgements Financial support from US AFRL AOARD Grant No. 08-4026 and NUS Incentive Grant No. R-265-000-240-731 is gratefully acknowledged. Appendix A. Unit cell models and amplification factors for micromechanics analysis In the micromechanics analysis, an idealized periodical representative volume element (RVE) or a unit cell with appropriate periodical boundary conditions is used to model the ply structure with real fibre arrangements. The micro strain or stress, i.e. strain or stress at the constituent level (fibre or matrix), is calculated with amplification factors through some reference points in the unit cell model. Usually two typical unit cell models, i.e. square model and hexagonal model as shown in Fig. A1, are used in the analysis. The geometry sizes of the models are usually given unit value, p i.e. ffiffiffi L = W = 1 for both models, H = 1 for the square model and H ¼ 3 for the hexagonal model. Some reference points at the critical locations in the unit cell model are selected to calculate the micro stresses (in both fibre and matrix). These reference points usually represent the ‘‘critical positions’’ within the fibre and matrix materials, such as the fibre-matrix interface positions, inter-fibre positions and interstitial positions. Fig. A2 illustrates the reference points in the abovementioned two unit cell models. The mechanical and thermo-mechanical amplification factors, M and A, are extracted from the unit cell models (at the reference points as shown in Fig. A2), which are applied periodical boundary conditions and load cases, with traditional linear finite element analyses. Generally, M and A may be strain or stress type, i.e. strain or stress amplification factors, which can be computed with the following formulation:

 1 ; ½M A ¼ XX

ðA1Þ

where

( Fig. A1. Unit cell models for micromechanics analysis: (a) square mode and (b) hexagonal model.

 ¼ ½x I X X ¼ ½xI

 II x  III x  IV x V x  VI x  VII 77 ; x xII xIII xIV xV xVI xVII 67 ;

Fig. A2. Reference points in unit cell models: (a) matrix of square model; (b) matrix of hexagonal model and (c) fibre of both square and hexagonal models.

ðA2Þ

1114

X.S. Sun et al. / Computers and Structures 89 (2011) 1103–1116

 i ði ¼ I; II; . . . ; VIIÞ denote seven macro load case arrays, which and x may be assigned with strain type (displacement) or stress type (force); xi (i = I, II, . . . , VII) denote the micro strain or stress arrays at the reference points from the finite element analysis results un i , respectively; for strain-type loads, x  i and xi are assigned as der x followings:

(

12 c 23 c 13 DTT ;  i ¼ ½e1 e2 e3 c x i ¼ I; II; . . . ; VII xi ¼ ½e1 e2 e3 c12 c23 c13 T ; ðA3aÞ

where the strain components are only mechanical strains, and for stress-type loads:

(

1 r 2 r 3 s 12 s 23 s 13 DTT ;  i ¼ ½r x i ¼ I; II; . . . ; VII xi ¼ ½r1 r2 r3 s12 s23 s13 T : ðA3bÞ

In the micromechanics analysis based on the unit models in Fig. A1, it is assumed that each plane surface of the unit cell model is always kept a plane state during the deformation, and the corresponding boundary conditions are actually a special explicit form of periodic boundary conditions suitable for the parallelepiped unit cell models [32,34]. It can be seen that the macro strain or stress applied to the unit models is uniform with these boundary conditions. Usually, any seven valid cases of uniform macro load,  i ði ¼ I; II; . . . ; VII), among which any two cases are linearly x independent from the mathematical view point, can be used to compute the corresponding micro strains or stresses xi (i = I, II, . . . , VII) for extracting the strain or stress amplification

factors with Eq. (A1). Some special load cases, however, can simplify the computation procedures. Consider three cases of uniaxial tension load along 1-, 2- and 3-directions, respectively, three cases of simple shear load on planes perpendicular to 1-, 2- and 3-directions, respectively, and one case of pure thermal load, with unit value for all load components. Under these cases of macro load,  will become a unit matrix, i.e. X  ¼ I, which forms a kind of X macro ‘‘unit pivot load cases’’, and the strain or stress amplification factors will be easy to obtain from the micro strains or stresses with Eq. (A1). This kind of macro ‘‘unit pivot load cases’’, including the strain type and stress type, can be constructed by appropriate boundary conditions (prescribed displacements and forces if necessary) on the unit cell models (Fig. A1), as listed in Tables A1 and A2, respectively. If the stress-type macro loads (Table A2) are applied to the unit cell models, the stress amplification factors Mr and Ar can be directly obtained from Eq. (A1); otherwise, if the strain-type macro laods (Table A1) are applied, the strain amplification factors Me and Ae can be directly obtained from Eq. (A1) and the stress amplification factors Mr and Ar are calculated from the following formulations:

(

 M r ¼ DM e S;

ðA4Þ

Ar ¼ DAe ;

where D is the constituent (fibre or matrix) stiffness and S is the ply compliance. Alternatively, the thermo-mechanical stress amplification factor Ar can also be computed by other simple periodical boundary conditions without multiple-point constraints after the mechanical

Table A1 Boundary conditions for the macro strain-type unit pivot load cases.  i) Pivot load cases (x

Surface 1 (x1 = 0)

Surface 2 (x1 = L)

Surface 3 (x2 = W/2)

Surface 4 (x2 = W/2)

Surface 5 (x3 = H/2)

Surface 6 (x3 = H/2)

1-tension ( e1 ¼ 1) 2-tension ( e2 ¼ 1) 3-tension ( e3 ¼ 1) 12 ¼ 1) 12-shear (c

u1 = 0 u1 = 0 u1 = 0 u2 = 0 u3 = 0 u1 = 0

u1 = L u1 = 0 u1 = 0 u2 = 0 u3 = 0 u1 = 0

u2 = 0 u2 = 0 u2 = 0 u1 = 0

u2 = 0 u2 = W u2 = 0 u1 = W

u3 = 0 u3 = 0 u3 = 0 u3 = 0

u3 = 0 u3 = 0 u3 = H u3 = 0 u2 = H

u2 = 0 u3 = 0 u1 = MPCa F1 = 0

u1 = 0 u3 = 0 u2 = 0

u2 = 0

u2 = 0 u3 = 0 u1 = 0

u1 = 0 u3 = 0 u2 = 0

u1 = 0

u1 = H

u2 = 0

u2 = MPC F2 = 0

u3 = 0

u3 = MPC F3 = 0

23 ¼ 1) 23-shear (c 13 ¼ 1) 13-shear (c Temperature (DT = 1) a

Multiple-point constraints which make all the points at the surface have the same non-zero unknown displacement along the prescribed direction (similarly hereinafter).

Table A2 Boundary conditions for the macro stress-type unit pivot load cases.  i) Pivot load cases (x

Surface 1 (x1 = 0)

Surface 2 (x1 = L)

Surface 3 (x2 = W/2)

Surface 4 (x2 = W/2)

Surface 5 (x3 = H/2)

Surface 6 (x3 = H/2)

 1 ¼ 1) 1-tension (r

u1 = 0

u2 = 0

u1 = 0

 3 ¼ 1) 3-tension (r

u1 = 0

12 ¼ 1) 12-shear (s

u2 = 0

23 ¼ 1) 23-shear (s

u1 = 0

13 ¼ 1) 13-shear (s

u3 = 0

Temperature (DT = 1)

u1 = 0

u2 = MPC F2 = 0 u2 = MPC F2 = L  H u2 = MPC F2 = 0 u1 = MPC F1 = L  H u3 = MPC F3 = L  H u2 = MPC F2 = 0 u2 = MPC F2 = 0

u3 = 0

 2 ¼ 1) 2-tension (r

u1 = MPC F1 = HW u1 = MPC F1 = 0 u1 = MPC F1 = 0 u2 = MPC F2 = HW u1 = MPC F1 = 0 u3 = MPC F3 = HW u1 = MPC F1 = 0

u3 = MPC F3 = 0 u3 = MPC F3 = 0 u3 = MPC F3 = L  W u3 = MPC F3 = 0 u2 = MPC F2 = L  W u1 = MPC F1 = L  W u3 = MPC F3 = 0

u2 = 0 u2 = 0 u1 = 0 u3 = 0 u2 = 0 u2 = 0

u3 = 0 u3 = 0 u3 = 0 u2 = 0 u1 = 0 u3 = 0

1115

X.S. Sun et al. / Computers and Structures 89 (2011) 1103–1116

stress amplification factor Mr is obtained. Consider the boundary conditions that normal displacements on all surfaces are prescribed zero values, i.e. u1 = 0 on x1 = 0 and x1 = L, u2 = 0 on x2 = ±W/2, u3 = 0 on x3 = ±H/2. Under these boundary conditions and a unit temperature increment DT = 1, the micro stress array in the constituent (fibre or matrix) of the unit cell model, rFE, can be obtained from the finite element analysis, and then Ar can be computed with:

a ;  ¼ rFE þ M r D Ar ¼ rFE  M r r

ðA5Þ T

 is the ply stiffness; a 2 a  3 0 0 0  is the ply 1 a  ¼ ½a where D  eMech ¼ DT D a  ðDT ¼ 1Þ; eMech and r ¼D  thermal coefficient array; r are the macro mechanical strain and stress arrays, respectively, . with eMech ¼ DT a Appendix B. Correction coefficients in the stress amplification factors Under the tensile loading along the longitudinal direction, the macro and micro longitudinal strains are same (no longitudinal strain amplification). Therefore the macro longitudinal strain ea corresponding to the matrix tensile failure strain can be written as:

ea ¼ hðemT  em0 Þ;

ðA6Þ

where em0 = Ae1DT is the thermal strain in the matrix; Ae1 is the matrix thermal-mechanical strain amplification factor along the longitudinal direction; h = Tm/YmT is a reduced factor defined by Tm and the pure matrix tensile strength YmT; emT is the matrix tensile failure strain which can be chosen as emT = 0.04 for carbon/epoxy and emT = 0.05 for glass/epoxy. According to the elastic stress–strain relationship, the macro  a , can be written uniaxial stress along the longitudinal direction, r as:

r a ¼ E1 ea ¼ T m E1 =E0m ;

ðA7Þ E0m

where E1 is the ply longitudinal tension modulus and is the matrix secant modulus at the point corresponding to the tensile failure strain, and

E0m ¼ Y mT =ðemT  em0 Þ:

ðA8Þ

1 ¼ r  a , the micro stresses in Under the macro uniaxial stress r the matrix are calculated with Eq. (4b) and substituted into Eq. (2). Then the correction coefficient c1 in the modified stress amplification factors can be obtained by solving the quadratic equation c21 þ 2a0 c1 þ b0 ¼ 0 as:

c1 ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a20  b0  a0 ;

ðA9Þ

where

(

 a þ Ra Þ=ð2Mr11 r  a Þ; a0 ¼ ðR1 r  2a þ ðC m  T m ÞR1 r  a þ Rb =ðM r11 r  a Þ2 ; b 0 ¼ ½ R2 r

(

(

ðA10Þ

P

Ra ¼ ð3Ar1  3j¼1 Arj ÞDT þ ðC m  T m Þ;  a þ s 0 Þ DT  T m C m ; Rb ¼ ðs1 r R1 ¼ R2 ¼

ðA11Þ

P3

j¼2 M rj1 ;

1 2

P3

2 j¼2 ð3M rj1



P3

l¼2 M rj1 M rl1 Þ

þ 3M 2r51 ;

ðA12Þ

8 P  P > < s1 ¼ 3j¼2 3M rj1 Arj  3l¼1 M rj1 Arl þ 6M r51 Ar5 ; h P  i P P > : s0 ¼ 12 3j¼1 3A2rj  3l¼1 Arj Arl þ 3A2r5 DT þ ðC m  T m Þ 3j¼1 Arj ; ðA13Þ

and Mrjl and Arj (j = 1, 2, . . . , 6; l = 1, 2, . . . , 6) are the matrix stress amplification factors. As the ply failure under the shear stress is dominated by matrix, the correction coefficient ck (k = 4, 6) can be determined by the ply in-plane shear strengths S12 and S13 (usually S13 = S12), respectively. 12 ¼ S12 (or s 13 ¼ S13 ) is applied, When the macro uniaxial stress s the matrix will fail and its failure envelope should pass through the point corresponding to the macro uniaxial stress state in the macro stress space. Similar to the determination of c1, the micro 12 ¼ S12 (or s 13 ¼ S13 ) are stresses in the matrix calculated from s substituted into Eq. (2) and then the coefficient ck (k = 4, 6) is obtained:

ck ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T m C m  s 0 DT qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; S12 3ðM2r4k þ M2r6k Þ

k ¼ 4; 6;

ðA14Þ

where s0 is referred to the second formula in Eq. (A13). The above correction coefficients, c1, c4 and c6, are calculated corresponding to the stress amplification factors at each reference point, respectively.

References [1] Icardi U, Locatto S, Longo A. Assessment of recent theories for predicting failure of composite laminates. Appl Mech Rev 2007;60(2):76–86. [2] Garnich MR, Akula VMK. Review of degradation models for progressive failure analysis of fiber reinforced polymer. Appl Mech Rev 2009;62(1). No. 010801: 1-33. [3] Orifici AC, Herszberg I, Thomson RS. Review of methodologies for composite material modelling incorporating failure. Compos Struct 2008;86(13):194–210. [4] Tay TE, Liu G, Tan VBC, Sun XS, Pham DC. Progressive failure analysis of composites. J Compos Mater 2008;42(18):1921–66. [5] Zhang YX, Yang CH. Recent developments in finite element analysis for laminated composite plates. Compos Struct 2009;88(1):147–57. [6] Soden PD, Hinton MJ, Kaddour AS. A comparison of the predictive capabilities of current failure theories for composite laminates. Compos Sci Technol 1998;58(7):1225–54. [7] Hinton MJ, Kaddour AS, Soden PD. A comparison of the predictive capabilities of current failure theories for composite laminates, judged against experimental evidence. Compos Sci Technol 2002;62(12–13):1725–97. [8] Kaddour AS, Hinton MJ, Soden PD. A comparison of the predictive capabilities of current failure theories for composite laminates: additional contributions. Compos Sci Technol 2004;64(3–4):449–76. [9] Hinton MJ, Kaddour AS, Soden PD. A further assessment of the predictive capabilities of current failure theories for composite laminates: comparison with experimental evidence. Compos Sci Technol 2004;64(3-4):549–88. [10] Gosse JH, Christensen S. Strain invariant failure criteria for polymers in composite materials. In: Proceedings of the 42nd AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics, and materials conference and exhibit, AIAA2001-1184. Seattle, USA; April 16–19, 2001. [11] Tay TE, Tan SHN, Tan VBC, Gosse JH. Damage progression by the elementfailure method (EFM) and strain invariant failure theory (SIFT). Compos Sci Technol 2005;65(6):935–44. [12] Garnich MR, Hansen AC. A multicontinuum theory for thermalelastic finite element analysis of composite materials. J Compos Mater 1997;31(1):71–86. [13] Gotsis PK, Chamis CC, Minnetyan L. Prediction of composite laminate fracture: micromechanics and progressive fracture. Compos Sci Technol 1998;58(7):1137–50. [14] Huang ZM. A bridging model prediction of the tensile strength of composite laminates subjected to biaxial loads. Compos Sci Technol 2004;64(3– 4):395–448. [15] Ha SK, Jin KK, Huang Y. Micro-Mechanics of Failure (MMF) for continuous fiber reinforced composites. J Compos Mater 2008;42(18):1873–95. [16] Cheung MS, Akhras G, Li W. Progressive failure analysis of composite plates by the finite strip method. Comput Methods Appl Mech Eng 1995;124(1– 2):49–61. [17] Akhras G, Li WC. Progressive failure analysis of thick composite plates using the spline finite strip method. Compos Struct 2007;79(1):34–43. [18] Zahari R, EI-Zafrany A. Progressive failure analysis of composite laminated stiffened plates using the finite strip method. Compos Struct 2009;87(1):63–70. [19] Tay TE, Tan VBC, Tan SHN. Element-failure: An alternative to material property degradation method for progressive damage in composite. J Compos Mater 2005;39(18):1659–75. [20] Tay TE, Liu G, Tan VBC. Damage progression in open-hole tension laminates by the SIFT–EFM approach. J Compos Mater 2006;40(11):971–92.

1116

X.S. Sun et al. / Computers and Structures 89 (2011) 1103–1116

[21] Kwon YW, Eren H. Micromechanical study of interface stresses and failure in fibrous composites using boundary element method. Polym Polym Compos 2000;8(6):369–86. [22] Han R, Ingber MS, Schreyer HL. Progression of failure in fiber-reinforced materials. CMC-Comput Mater Continua 2006;4(3):163–76. [23] Blázquez A, Manticˇ V, París F, McCartney NL. Stress state characterization of delamination cracks in [0/90] symmetric laminates by BEM. Int J Solids Struct 2008;45(6):1632–62. [24] Blázquez A, Manticˇ V, París F, McCartney NL. BEM analysis of damage progress in 0/90 laminates. Eng Anal Bound Elem 2009;33(6):762–9. [25] Hettich T, Hund A, Ramm E. Modeling of failure in composites by X-FEM and level sets within a multiscale framework. Comput Methods Appl Mech Eng 2008;197(5):414–24. [26] Huynh DBP, Belytschko T. The extended finite element method for fracture in composite materials. Int J Num Methods Eng 2009;77(2):214–39. [27] Sun XS, Tan VBC, Liu G, Tay TE. An enriched element-failure method (REFM) for delamination analysis of composite structures. Int J Num Methods Eng 2009;79(6):639–66.

[28] Guiamatsia I, Falzon BG, Davies GAO, Iannucci L. Element-Free Galerkin modelling of composite damage. Compos Sci Technol 2009;69(15–16):2640–8. [29] Djeukou A, von Estorff O. Static and damage analyses of shear deformable laminated composite plates using the radial point interpolation method. In: Ferreira AJM et al., editors. Progress on meshless methods. Computational methods in applied sciences. Dordrecht: Springer; 2009. p. 199–216. [30] Dvorak GJ, Benveniste Y. On transformation strains and uniform fields in multiphase elastic media. Proc Roy Soc London Ser A 1992;437(1900): 291–310. [31] Dvorak GJ. Transformation field analysis of inelastic composite materials. Proc Roy Soc London Ser A 1992;437(1900):311–27. [32] Jin KK, Huang Y, Lee YH, Ha SK. Distribution of micro stresses and interfacial tractions in unidirectional composites. J Compos Mater 2008;42(18):1825–49. [33] Kim R, Sihn S. Effects of Stacking Sequence. In: 9th composites durability workshop (CDW-9). Shanghai, China; 1–2 November, 2004. [34] Xia Z, Zhang Y, Ellyin F. A unified periodical boundary conditions for representative volume elements of composites and applications. Int J Solids Struct 2003;40(8):1907–21.