A modified model for simulation of mode I delamination growth in laminated composite materials

Accepted Manuscript A modified model for simulation of mode I delamination growth in laminated composite materials M.M. Shokrieh, Z. Daneshjoo, M. Fakoor PII: DOI: Reference:

S0167-8442(15)30132-4 http://dx.doi.org/10.1016/j.tafmec.2015.12.012 TAFMEC 1655

To appear in:

Theoretical and Applied Fracture Mechanics

Please cite this article as: M.M. Shokrieh, Z. Daneshjoo, M. Fakoor, A modified model for simulation of mode I delamination growth in laminated composite materials, Theoretical and Applied Fracture Mechanics (2015), doi: http://dx.doi.org/10.1016/j.tafmec.2015.12.012

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.

A modified model for simulation of mode I delamination growth in laminated composite materials M. M. Shokrieha1, Z. Daneshjooa and M. Fakoorb a

Composites Research Laboratory, Center of Excellence in Experimental Solid Mechanics and Dynamics, School of Mechanical Engineering, Iran University of Science and Technology, Tehran, Iran, [email protected], [email protected] b Faculty of New Sciences and Technologies, University of Tehran, Tehran, Iran, [email protected]

Abstract In this paper, a modified model has been proposed for simulation of mode I delamination growth in laminated composite materials. The proposed model called “multi-linear cohesive zone model” utilizes spring elements to consider the residual strength of the damaged material in the fracture process zone. As the crack propagates, stiffness characteristics of spring elements change gradually according to the traction-separation curve. In this model, unlike the previous ones, all possible softening behaviors of the traction-separation curve, even nonlinear behavior, could be taken into account. An analytical equation based on the energy approach has also been proposed for extraction of mode I delamination R-curve. The multi-linear cohesive zone model along with the R-curve extracted by the model considers the fracture process zone effects. A comparison of the R-curve obtained by the new model with the available experimental data shows the capabilities of the proposed method. Keywords Delamination, fracture process zone, R-curve, composite materials, traction-separation curve

1. Introduction Linear elastic fracture mechanics (LEFM) has been developed to study a variety of fracture problems. It has been proved that if the fracture process zone (FPZ) ahead of the crack tip is small, then LEFM approach will be effective in predicting crack growth [1–3]. However, in many cases of crack growth in composite materials and structures made of other quasi-brittle materials, because of toughening mechanisms in the FPZ (such as microcracking, fiber bridging, coalescence of voids and other mechanisms), the FPZ may be large. So, the LEFM is not suitable fracture analysis of these materials [4, 5]. Therefore, nonlinear fracture models that can accurately evaluate the effects of toughening mechanisms have been considered to investigate the crack growth in such materials [6-8]. Cohesive zone model (CZM) is a nonlinear fracture 1

Corresponding Author Composites Research Laboratory, Center of Excellence in Experimental Solid Mechanics and Dynamics, School of Mechanical Engineering, Iran University of Science and Technology, Tehran, Iran, 16846-13114.

mechanics method for analysis [9, 10] in which, the crack should not only overcome the crack resistance of the material, but also should neutralize the stresses generated by the cohesive zone [11]. Hillerborg [12] showed that if properties and behavior of material are well understood, the CZM is a general model that is applicable to different types of materials. He showed that this model has been mainly applied to ordinary concrete, fiber reinforced concrete, rock, carbon/epoxy composites, glass/epoxy composites, wood and metals. Many variations of CZMs (such as polynomial, triangular, exponential, trapezoidal, linear/bilinear, multi-section, etc.) have been proposed and applied to predict fracture behaviors of materials [13–16]. It has been shown that different softening behavior of traction-separation curves leads to different results [17, 18]. By considering that the traction-separation curves of various materials are different, Hillerborg [12] stated that for the proper use of the CZM for a variety of different materials, it is required to have a sufficient knowledge about the material properties, the FPZ and a procedure which covers all FPZ effects [19, 20]. It is well known that the "crack growth resistance curve" or R-curve is a suitable concept for quantifying the FPZ effects [21, 22]. Recent research of Sorensen et al. [23] and Tamuz et al. [24] demonstrated that a trilinear CZM can better predict the R-curve behavior. It has already been shown that as a result of several types of toughening mechanisms with different strain energy release rates in the FPZ of composite materials, the bilinear CZM does not cover all FPZ toughening events. So, it is required to propose a multi-linear CZM [25-27]. A review of published works available in the literatures reveals that the effects of considering FPZ on the delamination growth in laminated composite materials are not studied well. For a real simulation of the delamination growth, it is necessary to have a method capable of implementing any arbitrary shape of softening behavior in the damaged area. However, such a versatile method cannot be found in previous researches. Considering the importance of modeling of different forms of traction-separation curves, in this paper, a modified simulation model is proposed in which all possible softening behaviors of a traction-separation curve are taken into account. The proposed model called the multi-linear cohesive zone model (ML-CZM), using variable stiffness springs, considers FPZ effects to simulate the initiation and propagation of mode I delamination in laminated composites. For this purpose, ML-CZM along with an extracted analytical equation for R-curve has been employed. The proposed model was verified by the available experimental data.

2. Simulation of initiation and propagation of mode I delamination 2.1. Description of the multi-linear cohesive zone model

In quasi-brittle composite materials, particularly in fibrous composites, due to the toughening mechanisms such as fiber bridging, the crack growth phenomenon is associated with the creation of a FPZ at the crack tip. It is important to know that FPZ is capable to carry the load. Moreover, this zone at the crack tip of composites and quasi-brittle materials absorbs the fracture energy,

delays the fracture of the specimen and prevents the catastrophic failure. Thus, effects of FPZ must be considered in the modeling and calculation of the crack growth [28]. For this purpose, in the present research, a model called the multi-linear cohesive zone (ML-CZ) model is proposed. The ML-CZ model utilizes spring elements to consider the residual strength of the damaged material in the FPZ. As the crack propagates, stiffness characteristics of spring elements change gradually according to the traction-separation curve. Therefore, the mentioned toughening effects will be simulated by variable stiffness spring elements in the model. Since the delamination phenomenon takes place between the plies, the crack propagation path is usually known in advance. An example of the configuration of spring elements along the crack propagation path for a double cantilever beam (DCB) specimen and the fiber bridging zone of a unidirectional DCB specimen [27] are shown in Figs. 1(a) and (b). A general traction-separation curve is considered as shown in Fig. 2. The algorithm of the MLCZ model for simulation of initiation and propagation of mode I delamination in laminated composites is shown in Fig. 3. Also, a computer code has been provided in accordance with this algorithm. As shown in Fig. 3, the computation is performed in a step by step approach, i.e. a unit load is applied to the specimen and the system of equations is solved under this loading condition. Forces and displacements in the model are obtained from the numerical solution. Since the analysis is linear elastic, by using the definition of a proportionality coefficient, the amount of load that provides a force equivalent to maximum tensile strength in a spring with a maximum stress can be calculated. In this case, the nodes on both sides of the spring have reached a maximum stress according to the traction-separation curve. The maximum stress is defined for the FPZ, in accordance with the traction-separation curve ( Sut in Fig. 2). In this case, at the first step, a proportionality coefficient ( 11 ) is defined. Then, it is assumed that a FPZ is created in the material. It means that the spring with the maximum tensile strength still has the loading bearing capacity. The amount of this capacity depends on the amount of the crack opening length. Therefore, the stiffness of the described spring should be changed based on the traction-separation curve. In the next steps, unlike the first step, two events are expected to be occurred; 1) Another spring reaches the maximum tensile strength, or 2) A break point (the border between the two lines with different slopes) in the tractionseparation curve is achieved by another spring in the FPZ. Both of these events should be considered. After calculating all proportionality coefficients at each step, the smallest absolute value of the proportionality coefficient is chosen and the related event will be selected. After that, the stiffness of the spring with lowest proportionality coefficient must be changed based on the traction-separation curve. When the loading capacity in a spring goes to zero, the spring is removed from the model and the spring stiffness is considered to be zero. It is important to know that all springs in the FPZ in every step of the calculations are examined, especially in the case of various break points in the traction-separation curve or the specimen with complex geometry. This process continues until the last spring in the FPZ reaches the final break point in the traction-separation curve.

The main focus of this paper is on mode I delamination in laminated composites. The initiation and propagation of mode I delamination are predicted by the maximum strain energy release rate criterion. According to this criterion, the crack growth will happen when the energy release rate (GR ) is greater than or equal to a critical value (Gc ) . Depending on how the GR changes with the crack length, crack growth can be stable or unstable. Therefore, in the analytical modeling of delamination, the most important point is the estimation and calculation of GR . In the next section, using the results of the ML-CZ model, a comprehensive analytical relation for extraction of mode I delamination R-curve is developed. The ML-CZ model along with this analytical formulation for the R-curve can consider FPZ effects using the variable stiffness of spring elements and simulate the initiation and propagation of mode I delamination of laminated composites.

2.2. R-curve behavior of laminated composites A convenient tool for quantifying the FPZ effects is the crack growth resistance curve or Rcurve. The R-curve is a plot of the crack growth resistance GR as a function of crack length. An example of such a curve is shown in Fig. 4. Davila et al. [25] obtained a mathematical expression for the R-curve behavior based on an approximate method proposed by Foote et al. [25]. With the assumption that bridging stresses are not sensitive to crack-opening profile and depend only on the distance from the crack tip. Foote et al. [29] developed an approximate solution for the Rcurve related to the DCB specimen by a power softening law equation and simplified that for the bridging stress with a linear softening as follows: a  x GR (a)  Gi   c c (1  )dx a  lFPZ (1) lFPZ lFPZ 0 where  c is the critical displacement for the softening law and  c is the ultimate strength of the bridging tractions and a  lFPZ . For a linear softening:

c 

2Gb c

(2)

By substitution of Eq. (2) into Eq. (1) and performing an integration, Davila et al. [25] obtained an expression for the R-curve denoted as GRNL that is quadratic in terms of a when a  lFPZ : a a  Gi  Gb NL (2  NL ) lFPZ lFPZ G (a)    Gi  Gb  NL R

NL a  lFPZ

a  l

(3)

NL FPZ

where Gi is the initial fracture toughness, Gb is the energy release rate related to bridging effect, NL is the length of the FPZ for a nonlinear form of R-curve. A a is crack extension and lFPZ

simple expression for the R-curve of a large scale bridged crack with analysis of different forms of bridging stress was presented by Suo et al. [30]: a  L a  lFPZ Gi  Gb L L lFPZ GR (a)   (4) L  G G a  lFPZ i b  L where lFPZ was defined as follows:

L lFPZ L

Ez Gb

(5)

 c2

where Ez is elastic modulus and  c is the tensile strength. Regarding the prediction of the maximum traction, it should be equal to the transverse tensile strength of a transversely isotropic unidirectional composite ( YT ). For a uniform bridging stress, the FPZ coefficient,  L , calculated by Bao and Suo [31] is a constant (  L  0.393 ). However, for a linear softening bridging stress, the FPZ coefficient  L is a function of the crack extension a and according to Bao and Suo [31], ranges from 0.393 for short cracks to 0.732 for steady-state propagation. Bao nonlinear model for linear softening can be linearized as follows. By equating Eqs. (3) and (4), the following relation between  and  L was obtained by Davila et al. [25]:

 L (a) 

 2

a lFPZ

(6)

where  is constant and  L is a function of a . Eq. (4) is linearized using a constant value for

 L regardless of the softening law. An adequate constant value of  L was obtained by Davila et al. [25] by evaluating Eq. (6) at a  lFPZ 2 : 2 2 L , l FPZ  l FPZ 3 3 Davila et al. [25] defined the following expression for an R-curve that results from the superposition of two linear softening laws: a q a a (1  q) a GRNL (a)  qGC (2  )  (1  q)GC (2  ) lC p lC lC (1  p) lC

L  

G1 if

p a  lc q

G2

if

(7)

(8)

1 p a  l 1 q c

where lc   EGc  c2 is the length of the FPZ for a linear softening. They showed that Eq. (8) satisfies the following conditions: 1. If q  p , i.e. having  c1   c 2 in Fig. 5, then, Eq. (8) will be independent of q and p ; and it reverts to Eq. (9). 2. When a  lc p q , the first term on the right hand side of Eq. (8) will be equal to G1 .

3. When a  lc (1  p) (1  q) , the second term on the right hand side of Eq. (8) will be equal

to G2 . According to Fig. 5, Davila et al. [25] showed that if the 1 and 2 linear softening laws are ordered such that p q  (1  p) (1  q) , then the length of the FPZ for bilinear softening will be defined as the crack extension necessary to reach the steady-state value of GR  Gc is: NL lFPZ 

1 p lc 1 q

if

p q  (1  p) (1  q )

(9)

if

p q  (1  p) (1  q)

(10)

otherwise: NL lFPZ 

p lc q

NL where lFPZ is the length of FPZ in a nonlinear form of R-curve that results from the superposition

of the two linear softening laws (Fig. 5). Davila et al. [25] extracted a simpler R-curve expression using superposition of two linear softening laws, in Eq. (4) and assuming  L  2 3  :

G LR (a)  MIN [G1; qGc

3 a 3 a ]  MIN [G 2 ;(1  q)Gc ] 2 lc 2 lc

(11)

Davila et al. [25] showed that if the 1 and 2 linear softening laws are ordered such that p q  (1  p) (1  q) , then:

 3 a a 2 p Gc if   2 l l 3q c c   3 a 2 p a 2 1  p G LR (a)  G1  (1  q) Gc if   2 l 3 q l 3 1 q c c   a 2 1  p Gc if   (12) lc 3 1  q  Also, according to Eq. (7) (relation between the length of linear and nonlinear FPZs), the length of the FPZ in linearized R-curve with softening law is [25]: L lFPZ 

2 1 p lc 3 1 q

if

p q  (1  p) (1  q)

(13)

2.3. Analytical equation for extraction of mode I delamination R-curve In this section, the variation of the fracture toughness versus the crack length for initiation and propagation of mode I delamination of a laminated composites is calculated using the elastic potential energy of springs that simulated the FPZ in the ML-CZ model. Due to existence of the FPZ in the crack tip of composite materials, the variation of the fracture toughness via crack extension is usually described by the concept of resistance curves (R-curves), i.e., the crack growth resistance G R as a function of crack extension a as:

GR (a)  Gi  GFPZ

(14)

where a is the crack extension, G i is the strain energy release rate at the initiation of delamination, G FPZ is the FPZ delamination toughness and G R is the strain energy release rate (equal to the propagation of delamination toughness). So, according to the maximum strain energy release rate criterion, a crack will grow when: GR (a)  Gi  GFPZ  Gc (15) where G c is the critical strain energy release rate. The R-curve can be divided into three regions: a) The initiation of delamination region G i . b) The FPZ and raised toughening mechanisms in this zone G FPZ . c) The propagation delamination region G ss . In the following, the proportion of strain energy release rate in each of these three regions will be calculated in order to derive an analytical model for the R-curve.

2.3.1. Fracture toughness of initiation of delamination ( G i ) According to the ML-CZ model, the initiation of delamination occurs when the first spring in the crack tip reaches the final break point in the traction-separation curve. If the traction-separation curve is considered in a general form as shown in Fig. 2, fracture toughness in the initiation of delamination can be obtained from the following equation: n 1 Gi   ki wi 2 (16) i 1 2 where n is the number of break points in the softening zone of the traction-separation curve, k i is the stiffness of spring elements (the slope of the traction-separation curve in the softening zone) and w i is the difference between the opening (separation) of break points in the tractionseparation curve.

2.3.2. Fracture resistance of the process zone ( G FPZ ) In experimental procedures to measure the R-curve, it is observed that after the crack growths, the fracture toughness increases with a steep slope to a stable value. The reason of this increasing trend is the existence of the toughening mechanisms in the FPZ such as the fiber bridging at the back of the crack front that cause a resistance to the crack growth. In order to calculate the strain energy release rate of the toughening mechanisms in the FPZ, it is needed to determine the number of involved spring elements versus the crack tip opening. Also, it must be clarified that which break points in the traction-separation curve have been passed by a specific spring and each spring is located on the slope of the softening zone of the curve. Then, G FPZ can be calculated by the summation of the elastic potential energy of these springs at any moment.

2.3.3. Fracture toughness of delamination in the propagation state ( Gss ) The R-curve flattening at higher crack lengths indicates that the effect of toughening mechanisms is limited and the increase of the fracture toughness with the crack length is not permanent and after a fully development of the bridging zone, the fracture toughness reaches a constant value. In the ML-CZ model, the crack grows stably after the last spring in the FPZ reaches the final break point in the traction-separation curve. So, the steady-state value of GR is defined as follows: n 1 Gc  Gss  m ki wi 2 (17) i 1 2 where m is the number of springs in the FPZ. According to the above explanations, R-curve behavior of unidirectional DCB specimens can be expressed as a mathematical equation. Evaluation of linear and polynomial functions illustrates that the simplest function that can fit the

2 experimental data in the FPZ is a quadratic function. Thus, GR (a)  A(a)  B(a)  C is defined for the FPZ. To obtain the unknown coefficients A, B and C, the following three conditions must be applied: a) The initiation point of the FPZ ( 0, Gi ).

b) The end point of the FPZ ( aFPZ , Gss ). c) The R-curve has a downward curvature in this zone and must be tangent to ( aFPZ , Gss ) at the end point. By applying these three conditions, the following equation for the R-curve behavior of DCB specimens is achieved:  n 1 2 [ i 1 2 ki wi ]  Gi   n 1 a n 1 a 2 2(m  1)[ i 1 ki wi 2 ]( )  (1  m)[ i 1 ki wi 2 ]( )   2 aFPZ 2 aFPZ GR (a)    GFPZ  n 1  2 m[ i 1 2 ki wi ]    Gi GFPZ Gss

0  a  aFPZ

(18)

a  aFPZ

where n is the number of break points in the softening region of the traction-separation curve, k i is the stiffness of spring elements (the slope of the traction-separation curve in the softening zone), wi is the difference between the opening (separation) of break points in the tractionseparation curve, m is the number of springs in the FPZ, a is the crack extension (a  a0 ) and

aFPZ is the crack extension in the FPZ (aFPZ  a0 ) . The main advantage of the proposed analytical versatile equation for R-curve (Eq. 18) is that it only requires the material traction-separation curve, length of the FPZ and the number of springs

in the FPZ. The analytical equation can be used in numerical modeling of delamination based on fracture mechanics (for example, numerical simulation was carried out in the previous section).

3. Implementation of multi-linear cohesive zone (ML-CZ) model A user developed Matlab code in a combination with the commercial software has been used to implement the ML-CZ model for the simulation of mode I delamination growth in laminated composites. A finite element model was prepared by ANSYS14.0 to simulate a DCB specimen having 150  25  4.2 mm dimensions and a pre-crack length of 35 mm (Fig. 6). SOLID46 and COMBIN14 elements were used to model a double cantilever beam and the springs, respectively. Spring elements are configured along the crack propagation path. To make the finite element model of the DCB specimen, a hexahedral mapped mesh was used. Hence, the finite element model of the specimen consists of one element along the thickness, ten elements in the direction of the width and thirty elements along the length of the specimen. According to the results of the mesh sensitivity analysis, the convergence can be obtained by using at least twenty-three spring elements through the crack growth path. The traction-separation curve of the DCB specimen was shown in Fig. 7. This trilinear curve has been obtained from the experimental R-curve of the DCB specimen. This method is fully described by Heidari-Rarani et al. [26]. According to the experimental R-curve (Fig. 8), the length of the FPZ has been estimated about 8.6 mm. So, during the implementation of ML-CZ model, at least three spring elements are considered in the FPZ and the stop point for the presented algorithm will be set after the final fracture (i.e., fracture of last spring element in the FPZ). It is important to note that after the separation of the last spring element, the crack growth is unstable. Finally, the extracted load-displacement curve has been considered as the output of the ML-CZ model.

4. Results and discussion The ML-CZ model has been validated using available experimental results provided by HeidariRarani et al. [26]. They used DCB specimens made of 24-ply unidirectional E-glass/epoxy laminate. The mechanical properties of the material are summarized in Table 1. A Teflon film with a thickness of 17 µm has been used to create the initial crack. The DCB specimen has a total length of L  150 mm, a width of b  25 mm and initial crack lengths of a0  35 mm. The experimental load-displacement and corresponding R-curve of the DCB specimen obtained by Heidari-Rarani et al. [26] is illustrated in Fig. 8. More details about the experiments can be found in the method proposed by Shokrieh et al. [27]. Table 1 Mechanical properties of unidirectional E-glass/ML506.

Ex (GPa)

E y  Ez (GPa)

 xy   zy   zx

Gxy  Gxz (GPa)

G yz (GPa)

33.5

10.23

0.27

4.26

10.65

The results of simulation of ML-CZ model are presented in form of the load-displacement curve. As it can be seen, the result of the model is in a good agreement with the available experimental data. Also, softening behavior of the material after the delamination initiation has been simulated correctly (Fig. 9). In order to highlight the ability of ML-CZ model in considering the multi linear softening law and its accuracy, a comparison between ML-CZ model with a trilinear softening law and available experimental results provided by Morais et al. [18] has been done. They used DCB specimens were obtained from unidirectional  024  laminates. Laminated plies were cut from high strength carbon fibre (T300)/toughened epoxy. More details about the experiments can be found in [18]. In this case, three different models as shown in Fig.10 were used for softening    laws. These traction-separation curves were obtained from [18] with  u  45 MPa, which was the ply transverse strength, and k  106 N/mm3. The first model used a linear softening with GIc  1000 J/m2. For the second case, a bilinear softening with GIc,i   u b,1 2  200 J/m2, GIc,b   b,1 u 2  800 J/m2 and  b1  1 MPa was assumed. The last model was a trilinear

softening law with  b 2  0.1 MPa. In Fig. 11, load–displacement curves obtained by simulations with linear, bilinear and trilinear softening laws using the ML-CZ model are compared with the results of experiments. The comparison indicates that there are many differences between the p   curve obtained considering linear softening law and the experimental result. In this case, the specimen was damaged by subjecting it to a load much smaller than the actual one and clearly it could not consider the load carrying capacity of FPZ. The result of the bilinear

softening law is in a better agreement with the experimental data. However, Fig. 11 shows that it is possible to obtain a well fit with the trilinear softening law. It is confirmed that a multi linear softening law can predict the R-curve behavior and provides more accurate fitting by considering the effects of the FPZ. In the following, to evaluate the results of the proposed analytical equation, R-curve behavior of a double cantilever beam (DCB) is calculated using Eq. (18) and finally is compared with the experimental data obtained by Heidari-Rarani et al. [26] and Shokrieh et al. [27]. A comparison between the magnitudes of delamination toughness calculated with the analytical equation (Eq. (18)) and the experimental values is given in Table 2. Table 2 A comparison between analytical delamination toughness values with experimental data of [0◦]24 unidirectional DCB specimens.

a(mm)

GR ( J m2 )

Error %

0 2aFPZ  3.44 5 aFPZ  4.30 2 3aFPZ  5.16 5 4aFPZ  6.88 5

aFPZ  8.60

Experimental [26]

Analytical

163.344

186.205

13.99

393.728

424.547

7.83

446.150

465.510

4.34

472.898

499.029

5.52

527.039

543.718

3.16

557.000

558.615

0.29

According to the experimental R-curve (Fig. 8.(b)), the value of the initiation delamination toughness is equal to 98.2 J/m2 while this value found by the model is equal to G=163.344 J/m2. However, in this figure from G=98.2 J/m2 to G=163.344 J/m2, no delamination growth is observed in experimental results (Fig. 8.(b)). The reason for this difference is different definitions for the initiation delamination toughness used by the proposed analytical R-curve and the experiments. In the proposed model, the amount of energy required to increase the crack length is considered as the initiation delamination toughness, while, in experimental tests, the Rcurve is plotted from the moment that the energy is released. Therefore, in order to have the same definitions for both experimental and analytical methods, for a  0 in Table 2, G=163.344 J/m2 considered as the initiation delamination toughness and was compared with the G i obtained by the proposed analytical model. A comparison between the results of the proposed analytical R-curve and experimental data was shown in Fig. 12. As can be seen, this difference has no effect on the estimation of the R-curve behavior during the delamination initiation and propagation.

4.1. Influence of specimen geometry and relative crack size on the R-curve behavior of unidirectional DCB specimens In this section the effect of the thickness and initial crack length on the R-curve behavior of DCB specimens was examined. Thus, the R-curve behavior of two DCB specimens made of 18-ply and 24-ply unidirectional E-glass/epoxy laminate with four different initial crack lengths a0  35, 45,55,65 mm were examined using the proposed analytical equation. The analytical Rcurves of [0o]18 and [0o]24 unidirectional DCB specimens with different initial crack lengths are shown in Figs. 13 and 14, respectively. In all these curves, three regions including the initiation delamination point, the FPZ and the propagation delamination region were observed. Comparison of these curves shows that the variation in the initiation fracture toughness in the interval of 35  a0  65 mm is negligible. Also, the effect of the thickness on the initiation

fracture toughness is negligible. The initiation fracture toughness values in the range of (8  a0 2h  20) remain constant and equal to their mean value Gi  170.03 J/m2. Moreover, these curves show that almost all unidirectional DCB specimens have an identical slope in the FPZ. The values of FPZ length obtained from the above analytical R-curves for the unidirectional DCB specimen are summarized in Table 3. According to this table, it can be concluded that the changes of the FPZ length values in the range of (8  a0 2h  20) is negligible. The average length of the FPZ is equal to 9.6 mm. Furthermore, an accurate and complete observation of the FPZ (especially the fiber bridging zone in DCB specimens) during the test shows that the bridged fibers between the two arms of the DCB specimen remain intact up to the complete development of the FPZ. However, when the fracture toughness reaches a steady-state, fibers at the beginning of the FPZ start to break. The loss of these fibers does not mean that the length of the bridging zone is reduced. But, on the opposite side of the FPZ (just behind the delamination front) fibers bridging happens and in spite of the crack growth, the energy balance is restored. Table 3 Average length of the FPZ for [0o]18 and [0o]24 unidirectional DCB specimens. Length of FPZ (mm)

a0 (mm)

0  18

0  24

35 45 55 65

10.40 9.90 9.85 9.30

8.60 9.50 9.73 9.25

Finally, by a comparison of these curves, it is observed that change of the initial crack length does not affect the steady state fracture toughness values. However, the steady state fracture toughness unlike the initiation fracture toughness shows sensitivity to the thickness more than the initial crack length. Sue et al. [30] showed that the steady state fracture toughness did not dependent on the thickness of the specimen and it can be assumed as a material property for the DCB specimen. The average value of the steady state fracture toughness is equal to 510.06 J/m2. All analytical and experimental R-curves for unidirectional DCB specimens are illustrated in Fig. 15(a) and (b). It was observed that the FPZ is in a large scale in the unidirectional DCB specimen made of E-glass/epoxy, so we must consider the R curve effects. In the analytical solution, GR started from an average value of 170.03 J/m2 at a  0 mm and when the crack length reaches 9.6 mm the value of GR is 510.06 J/m2. Moreover, in experiments, GR started from an average value of 103.2  8.3 J/m2 at a  0 and when the crack reaches 9.6 mm, the value of GR is

534.7  39.6 J/m2 [27]. A comparison of these results confirms that the analytical R-curves are compatible with the experimental data for different thicknesses and initial cracks lengths. So, the

proposed analytical equation in this paper can extract the R-curve behavior of DCB specimens very well.

5. Conclusions In the present research, a modified model called the multi-linear cohesive zone (ML-CZ) model has been proposed to consider the FPZ effects on mode I delamination behavior of laminated composites. Using the spring elements with variable stiffness according to the traction-separation curve, this model simulates the residual strength of the material in the FPZ. This model, unlike the existing models and commercial codes, considers all possible softening behaviors of the traction-separation curve (n-linear CZM). The exponential behavior of the softening region of the traction-separation curve can be modeled by increasing the number of breaks. The results of simulation presented in forms of load-displacement curves are in good agreement with the available experimental data and also correctly show the softening behavior of the material after the delamination initiation. Using the results of the ML-CZ model, an analytical equation has been proposed to extract the mode I delamination R-curve. The ML-CZ model along with the R-curve equation can estimate the initiation and propagation of mode I delamination. This analytical equation calculates the initiation delamination toughness ( G i ), the propagation delamination toughness ( G ss ) and the FPZ delamination toughness ( G FPZ ) using the elastic potential energy of spring elements that simulate the FPZ along the delamination growth. The main advantage of the proposed analytical equation for R-curve is that it only needs the traction-separation curve of the material, the FPZ length and the number of spring elements in the FPZ. The R-curves obtained by the proposed model in this study are in good agreement with experimental ones. References [1] R.W. Hertzberg, Deformation and fracture mechanics of engineering materials, Fourth

Edition. John Wiley & Sons, INC, (1983). [2] M. Fan, D.K. Yi, Z.M. Xiao, “Fracture analysis for a sub-interface Zener–Stroh crack in a bimaterial plate under small-scale yielding condition”, Theoretical and Applied Fracture Mechanics, Vol. 76, pp. 60–66, (2015). [3] D. Xie, S.B. Biggers Jr, “Strain energy release rate calculation for a moving delamination front of arbitrary shape based on virtual crack closure technique, Part I: Formulation and validation”, Engineering Fracture Mechanics, Vol. 73, No. 6, pp. 771–785, (2006). [4] A. Urala, V. R. Krishnanb, K. D. Papoulia “A cohesive zone model for fatigue crack growth allowing for crack retardation,” International Journal of Solids and Structures, Vol. 46, pp. 2453–2462, (2009).

[5] M.F.S.F. Moura and J.A.G. Chousal, “Cohesive and continuum damage model applied to

fracture characterization of bonded joints”, International Journal of Mechanical Sciences, vol. 48, No. 5, pp. 493-503, (2006). [6] D. S. Dugdale, “Yielding of steel sheets containing slits”, Journal of the mechanics and physics of solids, Vol. 8, pp. 100-108 (1960). [7] G. Wang, S.F. Li, “A penny-shaped cohesive crack model for material damage”, Theoretical and Applied Fracture Mechanics, Vol. 42, No.3, pp. 303-316, (2004). [8] P. Qiao, Y. Chena, “Cohesive fracture simulation and failure modes of FRP–concrete bonded interfaces”, Theoretical and Applied Fracture Mechanics, Vol. 49, No.2, pp. 213-225, (2008). [9] A. Hillerborg, “Analysis of fracture by means of the fictitious crack model, particularly for fiber reinforcement concrete”, International Journal of Cement Composites, Vol. 2, pp. 177184, (1980). [10] B. Ren, S. Li, “Multiscale modeling and prediction of bonded joint failure by using an adhesive process zone model”, Theoretical and Applied Fracture Mechanics, Vol. 72, pp. 7688, (2014). [11] F. Tin-Loi and H. Li, “Numerical simulations of quasibrittle fracture processes using the discrete cohesive crack model”, International Journal of Mechanical Sciences, Vol. 42, pp. 367-379, (2000). [12] A. Hillerborg, “Application of the fictitious crack model to different types of materials”, International Journal of Fracture, Vol. 51, pp. 95-102, (1991). [13] Tvergaard V, Hutchinson JW. The relation between crack growth resistance and fracture parameters in elastic–plastic solids. J Mech Phys Solids 1992; 40:1377–97. [14] Xu XP, Needleman A. Numerical simulation of fast crack growth in brittle solids. J Mech Phys Solids 1994; 42:1397–434. [15] Camacho GT, Ortiz M. Computational modeling of impact damage in brittle materials. Int J Solids Struct 1996; 33:2899–938. [16] Li S, Thouless MD, Waas AM, Schroeder JA, Zavattieri PD. Mixed-mode cohesive-zone models for fracture of an adhesively bonded polymer–matrix composite. Engng Fract Mech 2006; 73:64–78. [17] G. Alfano, “On the influence of the shape of the interface law on the application of cohesivezone models”, Comp. Sci. Technol., Vol. 66, pp. 723–730, (2006). [18] A.B. Morais, A .B. Pereira, “Application of the effective crack method to mode I and mode II interlaminar fracture of carbon/epoxy unidirectional laminates”, Comp. Part A, Vol. 38, pp. 785–794, (2007). [19] G.S. Amrutharaja, K.Y. Lama, B. Cotterell, “Fracture process zone concept and delamination of composite laminates”, Theoretical and Applied Fracture Mechanics, Vol. 24, No.1, pp. 5764, (1995). [20] A. R. Gowhari Anaraki, M. Fakoor, “Mixed mode fracture criterion for wood based on a reinforcement microcrack damage model”, Materials Science and Engineering, Vol. 527, pp. 7184–7191, (2010).

[21] G. Gatalanoti, A. Arteiro, M. Hayati, P.P. Camanho, “Determination of the mode I crack

resistance curve of polymer composites using the size-effect law”, Engineering Fracture Mechanics, Vol. 118, pp. 49-65, (2014). [22] Y.W. Mai, “Cohesive zone and crack-resistance R-curve of cementitious materials and their fiber-reinforced composites”, Engineering Fracture Mechanics, Vol. 69, pp. 219-234, (2002). [23] L. Sorensen, J. Botsis, Th. Gmür, L. Humbert, “Bridging tractions in mode I delamination: Measurements and simulations”, Composites Science and Technology, Vol. 68, pp. 2350– 2358, (2008). [24] V. Tamuzs, S. Tarasovs, U. Vilks, “Progressive delamination and fiber bridging modeling in double cantilever beam composite specimens”, Eng Fract Mech, Vol. 68, pp. 513–25, (2001). [25] C. G. Dávila, CA. Rose, PP. Camanho, “A procedure for superposing linear cohesive laws to represent multiple damage mechanisms in the fracture of composites”, International Journal of Fracture, Vol. 158, No. 2, pp. 211-223, (2009). [26] M. Heidari-Rarani, M. M. Shokrieh, P. Camanho, “Finite element modeling of mode I delamination growth in laminated DCB specimens with R-curve effects”, Composites: Part B, Vol. 45, pp. 897–903, (2013). [27] M. M. Shokrieh, M. Heidari-Rarani, M. R. Ayatollahi “Delamination R-curve as a material property of unidirectional glass/epoxy composites”, Materials and Design, Vol. 34, pp. 211– 8, (2012). [28] I. Smith, E. Landis and M. Gong, Fracture and Fatigue in Wood, John Wiley & Sons Ltd, England, (2003). [29] R. M. L. Foote, Y-W. Mai, B. Cotterell, “Crack growth resistance curves in strain-softening materials”, J Mech Phys Solids, Vol.34, pp. 593–607, (1986). [30] Z. Suo, G. Bao, B. Fan, “Delamination R-curve phenomena due to damage”, J. Mech. Phys. Solids, Vol. 40, No. 1, pp. 1-16, (1992). [31] G. Bao, Z. Suo, “Remarks on crack-bridging concepts”, Appl Mech Rev, Vol. 45, pp. 355– 366, (1992).

Fig. 1. (a) Configuration of spring elements along the crack propagation path in a DCB specimen (b) Fiber bridging zone in a unidirectional DCB specimen [27].

Fig. 2. Traction-separation curve (   w ) with multi-linear breaks in the softening zone.

Start Increase applying load : Fapplying×∆11 determination of σij for all spring elements Update stiffness of first spring: k11 =(∆σ1/∆w1) ×A1

Input Traction- Separation Curve: determination of wi σi kp wut Sut n (1≤ i ≤n+1) (1)

No

The last spring in the FPZ attains the final break point of the traction- separation curve?

j=j+1

(9) Yes

(5)

Input FE Modeling (applying Funit): determination of initial stiffness of spring elements (ki ), number of spring elements (K), Ai (1≤ i ≤K), σij i: counter of spring elements j: counter of step (2)

Calculation of ∆σi ∆wi (1≤ i ≤n): ∆ σi =σ(i+1)-σi (3) ∆wi =w(i+1)-wi Calculation of proportionality coefficient of first spring in the first step: ∆11 =Sut/σ11 , (i=1, j=1) (4)

For j≥2 Calculation of proportionality coefficients of 1st to jth spring elements: ∆jj=(1/σjj)×[Sut- ∑n=1:(j-1) ∆nn σjn] ∆ij=(1/σij)×[ σ(Last Break)- ∑n=1:(j-1) ∆nn σin] i
Choose of Minimum ∆ in the each step: T=min (∆jj, ∆ij) (7)

Increase applying load : Fapplying×T determination of σij for all spring elements Update stiffness of spring having minimum ∆ (8)

Calculation of Load (Fj) & Opening Displacement (δj) in the each step: Fj =∑n=1:j ∆n Funit δj =∑n=1:j ∆n δn (10) Output F-δ Curve

(11)

Calculation of Gi GFPZ and Gss: Gi=∑i=1:n ½ ki ∆wi2 Gss=m ∑i=1:n ½ ki ∆wi2

(12)

m is the number of springs in FPZ Calculation of GR(∆a) (13) Output R- Curve

Fig. 3. The algorithm of the ML-CZ model.

Fig. 4. The R-curve for a composite laminate under large-scale fiber bridging.

End

(14)

Fig. 5. Superposing of two bilinear CZMs [25].

Fig. 6. (a) Geometry of the DCB specimen, (b) spring elements and (c) the DCB finite element model.

Fig. 7. Traction-separation curve of a DCB specimen.

Fig. 8. Load–displacement and R-curve of a DCB specimen with a0  35 mm [26].

Fig. 9. Load-displacement curve obtained by ML-CZ model in comparison with the experimental data.

Fig. 10. Traction-separation curve of a DCB specimen used the ML-CZ model: linear (LS), bilinear (BS) and trilinear (TS) softening [18].

Fig. 11. Load-displacement curve obtained by the ML-CZ model in a comparison with the experimental data.

Fig. 12. R-curves obtained by the proposed analytical Eq. (18) and experiment.

Fig. 13. The R-curve of [0o]18 unidirectional DCB specimens with different initial crack lengths.

Fig. 14. The R-curve of [0o]24 unidirectional DCB specimens with different initial crack lengths.

Fig. 15. A comparison of (a) analytical and (b) experimental R-curve behavior of unidirectional DCB specimens with different thickness and initial crack lengths.

   

A modified model has been proposed for simulation of mode I delamination growth in laminated composite materials An analytical equation based on the energy approach has been proposed for extraction of mode I delamination R-curve The new proposed model along with the extracted R-curve considers the fracture process zone effects. A comparison of the extracted R-curve with available experimental data shows the capabilities of the proposed method.