A modified mode I cohesive zone model for the delamination growth in DCB laminates with the effect of fiber bridging

Journal Pre-proof

A modified mode I cohesive zone model for the delamination growth in DCB laminates with the effect of fiber bridging Yu Gong , Yixin Hou , Libin Zhao , Wangchang Li , Jianyu Zhang , Ning Hu PII: DOI: Reference:

S0020-7403(19)32998-4 https://doi.org/10.1016/j.ijmecsci.2020.105514 MS 105514

To appear in:

International Journal of Mechanical Sciences

Received date: Revised date: Accepted date:

12 August 2019 15 January 2020 5 February 2020

Please cite this article as: Yu Gong , Yixin Hou , Libin Zhao , Wangchang Li , Jianyu Zhang , Ning Hu , A modified mode I cohesive zone model for the delamination growth in DCB laminates with the effect of fiber bridging, International Journal of Mechanical Sciences (2020), doi: https://doi.org/10.1016/j.ijmecsci.2020.105514

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2020 Published by Elsevier Ltd.

Highlights 

A three-linear CZM is proposed based on the microscopic failure mechanism of the delamination in composite laminates.



Bridging stress is experimentally obtained by J-integral method and implemented into the new CZM using a user-subroutine UMAT.



A sensitivity analysis shows that the interfacial strength has little effect on the simulated results.

1

A modified mode I cohesive zone model for the delamination growth in DCB laminates with the effect of fiber bridging Yu Gonga,b,c,, Yixin Houd, Libin Zhaoe, *, Wangchang Lif, Jianyu Zhanga, *, Ning Hua a

b

College of Aerospace Engineering, Chongqing University, Chongqing 400044, China

Key Laboratory of Fundamental Science for National Defence of Aeronautical Digital Manufacturing Process, Shenyang Aerospace University, Shenyang 110136, China

c

Chongqing Key Laboratory of Heterogeneous Material Mechanics, Chongqing University, Chongqing 400044, China d

AECC Sichuan Gas Turbine Establishment, Mianyang 621000, China e

f

School of Astronautics, Beihang University, Beijing 100191, China

College of Materials Science and Engineering, Zhejiang University of Technology, Hangzhou, 310014, China

Abstract: Fiber bridging has a significant influence on the delamination propagation behavior in multidirectional composite laminates. Traditional pure mode I bilinear cohesive zone models (CZM) do not consider the effect of fiber bridging and result in an inaccurate simulation on the delamination behavior. This study proposed a physical-based three-linear CZM superposed by two bilinear CZMs, which represent two different phenomena including the quasi-brittle matrix fracture characterized by a higher peak stress and a shorter critical opening displacement, and the fiber bridging characterized by a lower peak stress and a longer critical opening displacement, respectively. The three-linear CZM was implemented in the commercial FE software using a user-subroutine UMAT. Double Cantilever Beam (DCB) tests on the multidirectional composite laminates with 0°/5° and 45°/-45° interfaces conducted in our previous studies are shown to have large-scale fiber bridging in mode I delamination and are used to provide experimental data for calibrating the new CZM. Good agreements between the predicted and tested results can be 

Corresponding author. E-mail address: [email protected] (Y. Gong); [email protected] (L. Zhao); [email protected] (J. Zhang). 2

achieved by adopting this new CZM, demonstrating its applicability on predicting the mode I delamination behavior in composite laminates with the effect of fiber bridging. Keywords: Carbon fiber; Laminate; Delamination; Finite element analysis (FEA)

1 Introduction One weakness of the widely used laminated fiber-reinforced composites is the risk of delamination, due to the lack of reinforcement in the thickness direction. Delamination is regarded as one of the most dangerous failure mechanisms in composite structure [1-3], as it develops inside of the material rather than obvious on the surface and finally results in a significant decrease of mechanical properties [4,5]. The characterization and prediction for delamination assume a prominent position in the design of composite structures. Designs of new laminated composite structures, currently based on expensive and time-consuming experiments, can greatly benefit from advanced numerical techniques. Some fracture mechanics methods have been implemented in FEM codes, such as the Virtual Crack closure Technique (VCCT) [6] to predict the delamination growth. However, there are several drawbacks which limit its application, including the inability to predict delamination onset and the need to pre-define the delamination path. There are also some recently developed numerical modelling technologies for simulating delamination in composite structures, e.g. cohesive zone model (CZM), embedded finite element method (E-FEM) [7], variational multiscale cohesive method (VMCM) [8], augment finite element method (A-FEM) [9], floating node method (FNM) [10], extended finite element method (XFEM) [11-13], extended cohesive damage model (ECDM) [14] and continuum decohesive finite element (CDFE) [15]. A comprehensive assessment on the VCCT, CZM and XFEM for the delamination growth modelling in composites and a detailed discussion on their advantages and disadvantages have been conducted by Heidari-Rarani and Sayedain [16]. The CZM can circumvent the aforementioned limitations and thus is an appealing alternative for modelling the delamination behavior. Different cohesive laws, which present the relation between 3

tractions acting on separated crack faces in a cohesive zone near the crack tip and the interfacial separation) have been proposed in the last four decades, such as bilinear, exponential, trapezoidal and cubic forms [17,18]. Turon et al. [19] proposed a thermodynamically consistent model for the simulation of delamination based on damage mechanics, where a novel constitutive equation is developed to account for crack closure effects to avoid interfacial penetration. Lindgaard el al. [20] proposed a mixed-mode CZM defined by the modified B-K criterion, where an adaptive numerical integration scheme was used for improved accuracy and convergence of elements. Alfano [21] conducted a benchmark study and comparison of prevalent CZM and concluded that the shape of the CZM law has trivial effect on the delamination analysis. However, this is only true for case of small-scale fracture process zone, where the size of fracture process zone is small compared with the crack length and other physical dimensions of specimens [18]. Composite multidirectional laminates are mostly used in practical engineering applications and delamination propagates in the interface between plies with different fiber orientations [22-24]. In composite unidirectional and multidirectional laminates, the large-scale fiber bridging usually occurs with the dimensions comparable to the specimen's size [25,26]. The intact bridging fibers act as closing tractions applied on the faces and result in an increased fracture toughness with the delamination growth, i.e. the resistance curve (R-curve) [27-31], which is dependent on the specimen's geometry (such as the thickness of specimen) and should not be considered as a material property [32-35]. In such case, the shape of the CZM law becomes important for delamination growth analysis and traditional bilinear CZMs cannot capture the delamination behavior well due to the ignorance of fiber bridging effect [17,18,25,36,37]. To include the effect of large-scale fiber bridging, Gong et al. [36,37] proposed to integrate the R-curve into the traditional B-K and power law criterions. Numerical simulations give good results compared to experimental ones for the standard MMB delamination tests. However, the proposed criterions are only effective for the one-dimensional delamination growth while not capable for two-dimensional delamination prediction. Developing new CZMs to present failure mechanisms on 4

both small and large scales is another approach. This can be achieved by adding complexity to the shape of the cohesive law [18]. In terms of this, Sorensen et al. [38] and Tamuzs et al. [39] applied a three-part cohesive law. The first part of the law is a bilinear traction-separation relation accounting for crack initiation and the second part characterizing the bridging traction distribution with an exponential decaying function. Jensen et al. [18] developed a multi-linear CZM with an arbitrary number of line segments enabling simulation of delamination in composite materials with R-curve effects. Sørensen et al. [40] determined the cohesive law by measuring the J-integral and end-opening of the cohesive zone of DCB specimens loaded with pure bending moments, which agreed well with mirco-mechanical modelling for the fiber bridging problem. Airoldi and Dávila [28] applied two procedures for determining the shape parameters of cohesive laws required to predict delamination growth behavior. The first method extracted the cohesive parameters from an accurate experimental R-curve by a new semi-analytical equation. The second method identified material parameter by minimizing the error between the averaged experimental and numerical load-displacement curves. Compared with the first method, the second one is advantageous when inaccuracies are brought by the fiber bridging in measuring the experimental fracture toughness. Canal et al. [41] proposed a new CZM to assess the thickness scaling effect associated with the fiber bridging. The new CZM was developed through a computation strategy utilizing multi-scale modeling approach which allows mapping micromechanical simulations into continuum finite element representations for macro-scale analyses. In addition, an embedded cell model of the DCB specimen was applied which explicitly accounts for the physical behavior including the delamination and the bridging bundles on the fracture plane. Farmand-Ashtiani et al. [42] measured the strain by means of embedded fiber bragg grating sensors and applied an iterative approach to identify the bridging law in the DCB specimens of different thickness. The identified bridging law is employed in a CZM to predict delamination in unidirectional DCB specimens. Daneshjoo et al. [43] proposed a novel mixed mode I/II micromechanical bridging model based on the calculation of the delamination bridging zone energy absorption to account for the fiber bridging effects. Gutkin et 5

al. [44] proposed a semi-analytical methodology for identifying the parameters of tri-linear cohesive laws from experimental R-curves for DCB and Compact Tension specimens. The cohesive parameters are obtained using assumed crack profiles and numerical analyses are used to calculate specimen compliances at different crack lengths. The methodology does not consider the effect of the bridging tractions on the assumed displacement profiles and the error in the toughness obtained by conventional data reduction schemes. Hansen et al. [45] proposed a CZM for simulating multi-scale failure mechanisms in the glass-epoxy laminates using a combination of bilinear and higher-order polynomial functions. Others attempted to use tri-linear cohesive laws, wherein one line segment is added to the traditional bilinear cohesive law to present a simple linear distribution of the bridging stress. For example, Dávila et al. [46] proposed a procedure for superposing two bilinear cohesive laws to approximate the experimentally determined R-curve. In the same way, Heidari-Rarani et al. [17] proposed a three-linear CZM with physical background, which is the superposition of two bilinear CZMs, presenting the quasi-brittle matrix fracture characterized by a higher peak stress and a shorter critical opening displacement, and the fiber bridging characterized by a lower peak stress and a longer critical opening displacement, respectively. This method is efficient due to its dependency on the R-curve while no requirement to measure the crack tip opening experimentally. The goal of this study is to propose a three-linear CZM with simple parameters required, which is based on the microscopic mechanism of the delamination failure. The organization of the paper is as followings: In Section 2, the traditional three-linear CZM in the literatures is briefly introduced and discussed. Afterwards, a new physical-based three-linear CZM is proposed. Tests used for validating the new CZM are introduced in Section 3. Section 4 establishes a numerical delamination model with the implementation of the new CZM. Comparisons between the experimental and numerical results obtained in composite laminates are made to verify the applicability of this new CZM. Finally, the influence of the interfacial strength on the numerical results is analyzed.

6

2 New three-linear cohesive zone model It has been shown that material softening and the R-curve are directly related to each other and that bilinear laws cannot accurately represent toughening mechanisms causing an R-curve response [17]. The method of superposition of two bilinear bridging laws is useful to obtain the three-linear CZM accounting for the toughening mechanisms, which avoids the need to develop cohesive elements with complex bridging laws. In this section, the traditional three-linear CZM in the literature is briefly introduced and discussed. A new physical-based three-linear CZM is subsequently proposed and its constitutive law is derived. 2.1 Traditional three-linear CZM in the literatures Fig. 1 shows the traditional three-linear CZM proposed in Refs. [17,46]. It is assumed that the two bilinear CZMs peak at the same displacement jump δ0. The sum of two arbitrary bilinear CZMs describing bridging phenomena results in a three-linear CZM. The blue line denotes a bilinear CZM with a short critical opening displacement δf1, characterizing the quasi-brittle matrix fracture or micro-cracking formation. The red line denotes a bilinear CZM with a longer critical opening displacement δf2, characterizing the fiber bridging. Each bilinear CZM includes three critical parameters: the initial interfacial stiffness K, the interfacial strength σ and the fracture toughness G. As illustrated in Ref. [46], the bilinear CZM can be described by the proportions of interfacial strength  01 =n 0 ,  02 = 1-n   0 , and the fracture toughness Gini=mGprop and Gbri=(1-m)Gprop with 0≤n, m≤1, so that Gprop=Gini+Gbri and σ0=  01   02 . The m and n are defined as the fracture toughness ratio and the strength ratio, respectively. Their value can be calculated by the following Eqs. (1) and (2) [47]. m

7

Gini Gprop

(1)

 2 (1  m) Ez  Gprop m 1 m ,  1   2 lcz 0 n 1 n  3 n  2  m Ez  Gprop , m  1  m   02 n 1 n  3 lcz

(2)

where γ is an empirical parameter which depends on the model and the suggested value is 0.884, lcz is the length of the fracture process zone. Ez is the through-the-thickness Young’s modulus. For transversely isotropic materials, Ez is assumed to be equal to the transverse Young’s modulus (E2). Once the m and n are determined, the bridging strength σb, which is the maximum bridging stress in the fiber bridging zone, can be calculated by Eq. (3) and the CZM constitutive law can thus be obtained. m   1 n  b = 0 1  n  1    n 1 m  

(3)

The CZM obtained by the superposition is efficient for the simulation of delamination propagation with the large-scale fiber bridging. However, in the traditional three-linear CZM, it is assumed that the two bilinear CZMs peak at the same displacement jump δ0. This assumption indeed does not conform to the realistic microscopic mechanism of delamination failure, and this will be discussed in the following Section 2.2. 2.2 Constitutive law of a new cohesive zone model Fig. 2(a) shows the new three-linear CZM. The blue and red lines denote a bilinear CZM characterizing the quasi brittle matrix fracture and the fiber bridging, respectively. At the initial stage (0≤δ<δ0), both the interface and fibers bear tensile loading simultaneously. However, the fibers are not peeled-off from the matrix and cannot bridge the delamination faces. In this case, the interface traction is dominant and much higher than the bridging traction. As the crack opening displacement increases, δ0≤δ<δb, the interface gradually damages with a decreasing load capability, as shown in Fig. 2(b). Because a new delamination still not occurs, the bridging fibers are wrapped up by matrix and the bridging traction is relatively small. When the crack opening displacement arrives the critical damage onset displacement of the interface (δ=δb), a new delamination surface is 8

created and the matrix starts to spalling [48]. The peeled-off bridging fibers hinder delamination propagation by restraining the fracture surfaces at this point, as shown in Fig. 2(c). When δb<δ<δf, the bridging fibers cause extra growth resistance individually. The new delamination surface is opened under tensile loading, damaged fibers carry loads until the final failure displacement δf is arrived, at which point the bridging fibers are broken. Different from the assumption in the traditional three-linear CZM that the two bilinear CZMs peak at the same displacement jump, the bridging traction peaks at the complete failure displacement of interface in this new CZM. This means the damage onset displacement of bridging fibers δb is same with the complete failure displacement of interface, which is consistent with the micro failure mechanism during the delamination growth in composite laminates with the effect of fiber bridging. The three-linear, softening constitutive behavior for single-mode loading shown in Fig. 2 can be defined as:

 = 1  d  Dij0  0,   0  1  K AB  1   0  ,      0 b   K0     d   1  K BC 1   f  ,     b f    K0      1,   f where d is a global damage variable. Dij0 is the initial stiffness tensor and defined as:

Dij0 = ij K0

(4)

(5)

(6)

The scalar parameter K0 is a penalty stiffness, which is equal to the sum of initial stiffness of CZM characterizing quasi brittle matrix fracture, K1, and that of CZM characterizing the fiber bridging, K2. KAB and KBC are the slope of AB and BC lines, respectively. They can be obtained by Eqs. (7)(8).

9

K AB 

 b   0  K 2 0  b  0

K BC 

b b  f

(7)

(8)

The parameters in the three-linear CZM include the initial fracture toughness Gini, the steady-state fracture toughness Gprop, the interfacial strength σ0 and the bridging strength σb, and the initial stiffness K1 and K2. The initial and steady-state fracture toughness can be determined by standard delamination tests. The bridging strength σb can also be experimentally determined by J-integral method using DCB tests. A referenced value of the initial stiffness K1 can be analytically prescribed from a general linear-elastic CZM as: K1=2Ez/h [49], h is the half thickness of a DCB specimen. The finally adopted value of K1 should be as high as possible (preferably infinite which is not possible for computational reasons) in order to ensure compatibility of the undamaged interface. As mentioned earlier, in the new CZM, it is assumed that the bridging traction peaks at the complete failure displacement of interface. Under this condition, there are only one independent parameter required to be numerically determined, i.e. the interfacial strengths σ0. According to Ref. [22], the interfacial strength is difficult to be measured by experiments. An important conclusion from this reference is that the value of the interfacial strength should be 55~65% of the interlaminar strength of the matrix, which is a useful guidance for choosing the interfacial strength in this study. 3 DCB test and results DCB tests are conducted to provide test data for validating the new CZM. Multidirectional laminates made of T700/QY9511 bismaleimide prepregs were manufactured. The material properties [50] are E11 = 130 GPa, E22 = E33 = 10.4 GPa and ν12 = ν13 = 0.3, and the tensile strength of the matrix is 100MPa. Specimens with two kinds of lay-up 016//(+5/-5/06)S and (+45/-45/06)S//(-45/+45/06)S were designed based on the classical laminated plate theory, where the symbol ‘//’ denotes the position of the artificial pre-crack introduced during the fabrication process. The 0°/5° interface was selected to alleviate the degree of fiber bridging, which exhibits very 10

similar mechanical properties with the unidirectional laminates. The 45°/-45° interface is very commonly applied in practical engineering and hence studied here. The values of Dc for both kinds of specimens are smaller than the upper limit value 0.25 [51] and values of Bt are also very low. Thus the bending-twisting coupling is limited and obvious non-uniformity of SERR along the width direction can be avoided. A 35mm-long and 25μm-thick Teflon film was inserted in the middle plane of laminates during the lay-up process to achieve an artificial pre-crack. After the curing process, the cured plates were cut by a diamond saw into specimens with specific geometry dimensions of 180 mm-long, 25 mm-wide and 4.16 mm-thick. Prior to DCB tests, the cutted specimens were inspected by the ultrasonically C-scanned technique to ensure defect-free specimens being used. The edges of each specimen were coated with a thin layer of white correction fluid in order to enhance visibility of the crack tip during the delamination growth process. It is worthwhile to point out that the use of white correction fluid may has the undesired effect that the delamination may propagate behind the thin white layer influencing the precise measurement of the crack tip location. All DCB tests were conducted according to the ASTM standard D5528-13 [52] and on a MTS 880 servo-hydraulic test machine. The displacement control mode was adopted at a low loading rate of 0.1mm/min. During the test, the displacement and load data at the loading point were automatically recorded by sensors with a pre-defined interval. An instrumented travelling microscope (JCXE-DK) was used to monitor and locate the position of the crack tip, which allowed a non-contact measurement of delamination length at the edge of the tested specimen with a precision of 0.01 mm. Three specimens were tested in DCB tests for each interface. The fracture toughness are calculated by the modified beam theory as the following equation:

GC 

3PC C F 2b  a    N '

(9)

where a, b, PC and C are the delamination length, the specimen width, the applied load and 11

displacement, respectively. F is a correction factor for considering a large displacement, N´ is a correction factor for considering the load-block effect, which is one due to the quick-mounted hinge used in this study. And  is a correction for crack tip displacement and rotation. The experimental results for this study had already been reported in our previous studies [22]. The relations between fracture toughness GC versus crack growth length ∆a for both interfaces are shown in Fig. 3(a), where significant R-curve behavior resulted by bridging fibers can be observed. The GC gradually increased with the crack growing at the initial stage until arriving a constant value after a certain distance. To quantitatively characterize the relations between the fracture toughness GC and ∆a, an empirical model [22,27,53] was adopted here as follows: 2   a   a  G +2  Gprop  Gini     a  lcz  -  Gprop  Gini    ； GC (a )   ini  lcz   lcz   a  lcz  Gprop ；

(10)

where ∆a is crack growth length. The lcz is defined as the length of the transient region where the R-curve develops from an initial value to reach a steady state value. The value of lcz can be experimentally determined by fitting the fracture toughness data via Eq. (10) with the best fitting degree. It mathematically means that the slope of the quadratic equation should be zero at Δa=lcz, where the steady-state condition is satisfied. The detailed values of Gini, Gprop and lcz for both interfaces are listed in Table 1. All specimens exhibit the same value of Gini (100 J/m2), indicating its independence on the ply orientation. However, the steady-state fracture toughness and the length of fiber bridging zone of both interfaces are different between each other. The quadratic model fits all GC data well with correlation coefficient r2 higher than 0.9, which indicates the high reliability of the GC values. The load versus displacement curves of DCB specimens for two kinds of specimens are presented in Fig. 3(b), where obvious difference can be observed for both interfaces. The structural stiffness, the initial damage load and the ultimate load are higher in the specimen with 0°/5° interface than in the specimen with 45°/-45° interface. The curve of 0°/5° interface is smoother after arriving at the ultimate load, while the 45°/-45° interface exhibited obvious stick-slip 12

behaviors. As presented in Fig. 4 showing the delamination paths [54], some delamination migrations occur in the specimens with the 45°/-45° interface and result in the stick-slip behaviors, meaning that delamination propagation partly happens at another interface than the original. However, only delamination strictly along the interface is modeled in the numerical simulation because it is usually deemed that the fiber bridging is the dominant factor which contributes to the most energy absorption [54,55]. More accurate simulation on the matrix cracking and the interactions between delamination and matrix cracking will be pursued in the next phase. 4 Implementation of CZM in FE model In this section, the numerical procedure for the simulation of the delamination propagation using the new three-linear CZM are described in the commercial finite element package ABAQUS ®, for implicit finite element analysis. The FE model is shown in Fig. 5. The three dimensional finite element mesh was prepared to simulate the DCB specimen with the 8-node solid brick element (C3D8). Every ply was individually modeled using one element through-the-thickness in order to consider the effect of fiber angle near the delamination plane. The 8-node cohesive elements (COH3D8) with 10-μm-thickness were pre-arranged along the middle plane of the laminates from the tip of the initial pre-crack to the end of the specimen. The DCB specimen was divided into three sections along the length with different mesh refinements. A refined mesh with 0.2mm length was used in the region where the crack was expected to propagate to ensure the independence of predicted results on the mesh size. The other regions along the length had solid elements with 4 mm length. There were eight elements through-the width of each arm. To implement the proposed new CZM in the cohesive elements, user subroutines were developed in the commercial software ABAQUS®. The upper and lower arms and interfaces were connected by tie constraints, which can ensure better mesh adaptability. For the cohesive elements, initial interfacial stiffness K1 of 1×1015 N/m3 and viscosity coefficient of 1×10-5 were adopted as suggested in Ref. [22]. The FE analysis was conducted under displacement control. Therefore, the loading points of the specimen were given a z-wise displacement, while the other end of the specimen was 13

fixed with rotation allowed. The implementation of new CZM includes three basic modules: 1) inputs of material parameters, material parameters will be transmitted to the user subroutine; 2) calculate the values of δ0, δb and δf, then obtain the damage variable d based on Eqs. (5)(7)(8) and stiffness matrix, finally update the stress matrix; 3) calculate the Jacobian matrix of the elements, update the Jacobian matrix and state variables. 5 Results and discussion The value of the interfacial strength is an important parameter in the new CZM. The suggested values of σ0 in the literatures are usually in the range of 35~75MPa [22]. In addition, the study from Ye and Chen [56] also illustrated that the σ0 is about 55–65% of the tensile strength of the matrix. In this study, the tensile strength of the QY9511 is 100MPa. It may suggest that a suitable value of the interfacial strength should be around 60MPa, which is also consistent with the conclusion in Ref. [56]. Therefore, the value of σ0 is adopted as 60Mpa and a sensitivity analysis of its effect on numerical results will be conducted in the following Section 5.2. In order to determine the σb, the J-integral method is adopted here. The crack opening displacement (COD) at the initial pre-crack should be recorded during the tests. Fig. 6 shows the relationship between the fracture toughness Gc and the COD at the initial pre-crack tip δ*. A non-linear fitting relation between them can be observed. The Gc dramatically increases with the δ* firstly until a slow increase trend exhibiting in the last stage. An exponential fitting function [57] as Eq. (11) is used to fit the G-δ* data. This is because that, as mentioned in Refs. [38,58,59], the bridging stress decreases monotonically from a maximum at the crack tip to zero at the end of bridging zone, which is also the case for the studied specimens here. Thus the exponential function is assumed to fit the G-δ* data in this study, which is also recommended by Shokrieh et al. [57,60]. The value of r2 is higher than 0.96 for both interfaces. G(  ) =Ga (1  e *



* a

)  Gb (1  e



* b

)  Gtip

(11)

where Ga, Gb, δa and δb are fitting parameters. The corresponding bridging stress distribution 14

equation is: *

*

dG( * ) Ga   a Gb  b  b ( )   ( )e  ( )e d * a b *

(12)

The bridging stress distribution can be calculated via Eq. (12). As shown in Fig. 7, it can be found that the bridging strength increases with the interfacial fiber angle. The bridging strength of specimens with 45°/-45° interface is 2.1MPa, which is higher than that of specimens with 0°/5° interface (1.0MPa). 5.1 Numerical results from the new model The numerical results using the traditional bilinear CZM are shown in Fig. 8 for both interfaces. It can be observed that numerical load versus displacement curves are totally different with the experimental ones, except for good agreements of the initial structural stiffness. The numerical ones exhibit rapidly linear increasing followed by a gradually non-linear reduction. In addition, it can be found that the ultimate load obtained by using the steady-state fracture toughness in the traditional bilinear CZM is higher than that by using the initial fracture toughness. Above results give evidence that the traditional bilinear CZM is not capable for simulating the propagated behaviors in the multidirectional DCB specimens. Comparisons between the numerical results obtained by the traditional three-linear CZM and the experimental results are also made. Using the traditional three-linear CZM, σ0=50Mpa, σb=0.4Mpa and σ0=55Mpa, σb=0.5Mpa are applied for the delamination simulation of the 0°/5° and 45°/-45° interfaces, respectively. The obtained numerical results are shown in Fig. 9. It can be seen that the traditional trilinear CZM is also capable of capturing the mode I delamination growth behavior. For the 0°/5° interface, some differences exist between the numerical and experimental results after the delamination onset. The value of σ0 is lower than the transverse tensile strength. The conclusions are consistent with those from Ref. [17]. Fig. 10 shows the predicted load versus displacement curves by using the new proposed CZM. It can be observed that the initial slopes of the numerical curves are in good agreements and keep 15

constant until the damage initiation denoted by the beginning of the degraded stiffness. The simulated stiffness is gradually and smoothly decreased with the delamination growth. Good consistent propagation phase with the experimental one can be easily observed in all figures. In addition, the predicted ultimate loads of the two kinds of DCB specimens and their comparisons with the experimental data are listed in Table 2, marked by bold text. The numerical errors are no more than 10% showing the high accuracy of the ultimate load predictions. Overall, the predictions by the new three-linear CZM show a good accordance with the experimental outcomes on the initial stiffness and ultimate loads as well as the delamination propagation behavior. Therefore, it follows that the new three-linear CZM proposed in Section 2.2 is efficient for accurately simulating the delamination of DCB specimens with the effect of fiber bridging. 5.2 The effect of interfacial strength Considering the realistic value of the interfacial strength are difficult to be obtained by experiments, a sensitivity analysis of the interfacial strength on the numerical results is conducted in this study. Four kinds of value (50MPa, 60MPa, 70MPa and 80MPa) are chosen for investigation. Fig. 11 presents the numerical load versus displacement curves originated from the new three-linear CZM for the 0°/5° and 45°/-45° interfaces. Comparisons between predictions and the corresponding test results are also shown in these figures. The detailed ultimate loads from the predictions and tests are listed in Table 2. For both interfaces, it can be seen that similar predicted results can be obtained when the interfacial strength is in the range of 50~80Mpa, which are all in good agreement with the test results. Take a closer look at the values of ultimate load, it can be found that the interfacial strength has little influence on the numerical results for 50Mpa≤σ0≤80Mpa. This is a helpful conclusion for the application of the new CZM as it indicates more applicability of this CZM. It is worth mentioning that in Ref. [61] an interesting sensitivity analysis is made on the all input parameters of the used models, which is not the main topic here and should be the subject for further researches. Currently the ability of the presented cohesive law to represent the experimental data is only 16

compared with the standard bilinear model and the traditional three-linear model rather than other state-of-the-art cohesive law shapes scoped for modelling fiber bridging effects, e.g. the exponential, multi-linear. This is because the goal of this study is not to achieve to a cohesive law with a better performance. The traditional trilinear CZM can give a closed relation for bridging stress versus n and m. However, in the proposed new CZM, the bridging stress is required to be determined by tests, which is a fundamental drawback. It may give very large bridging stress at the crack tip due to the sensitivity of bridging law to the very small crack opening displacement [62]. In the future, it is meaningful to determine the bridging stress by analytical methods with only basic experimental load and displacement data required. 6 Conclusions A physical-based three-linear cohesive zone model has been proposed to describe the effect of fiber bridging of the mode I delamination in multidirectional composite laminates. This model reflects the microcosmic mechanism of fiber bridging, and is implemented through the use-subroutine UMAT in the commercial finite element software ABAQUS®. To provide experimental data for the validation of the applicability of the proposed three-linear CZM, DCB tests on the multidirectional composite laminates with two kinds of interface (0°/5° and 45°/-45°) are conducted. The bridging strength is determined by the J-integral method. The simulations of the DCB tests are in good agreement with the experiment results for both interfaces, with the relative errors between the numerical and experimental results less than 10%, which illustrates the applicability of the new CZM on predicting the delamination behavior in multidirectional laminates with large-scale fiber bridging. Furthermore, a sensitivity analysis of the interfacial strength on the modelling results is carried out. It is shown that the interfacial strength has little effect on the simulated results, which again illustrates the applicability of the proposed new CZM. Though this new model is applied to the specific case of fiber reinforced composite laminates, it would be interesting to study its applicability for other material systems.

17

Acknowledgements The research work is supported by the National Natural Science Foundation of China (Project nos. 11902054, 11872131, 11772028, 11372020, 11572058 and U1864208), the Chongqing Natural Science

Foundation

(Projects

no.

cstc2018jcyjAX0235,

cstc2019jsyj-yzysbAX0015

and

cstc2019jscx-zdztzxX0028), the Fundamental Research Funds for the Central Universities (Project no. 2019CDXYHK0001) and the Key Laboratory of Fundamental Science for National Defence of Aeronautical Digital Manufacturing Process of Shenyang Aerospace University (Project no. SHSYS2018001).

Declaration of interests

☑ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References [1] Y Liu, C Zhang. A critical plane-based model for mixed-mode delamination growth rate prediction under fatigue cyclic loadings. Composites Part B: Engineering 2018; 139: 185-94. [2] OR Shah, M Tarfaoui. Determination of mode I & mode II strain energy release rates in composite foam core sandwiches. An experimental study of the composite foam core interfacial fracture resistance. Composites Part B: Engineering 2017; 111: 134-42. [3] Y Liv, G Guillamet, J Costa, EV González, L Marín, JA Mayugo. Experimental study into compression after impact strength of laminates with conventional and nonconventional ply orientations. Composites Part B: Engineering 2017; 126: 133-42. [4] VPW Shim, LM Yang. Characterization of the residual mechanical properties of woven fabric reinforced composites after low-velocity impact. Int J Mech Sci 2005; 47: 647-65. [5] MR Parlapalli, D Shu. Buckling of composite beams with two non-overlapping delaminations: Lower and upper bounds. Int J Mech Sci 2007; 49: 793-805. [6] R Krueger. Virtual crack closure technique: History, approach, and applications. Appl Mech Rev 2004; 57: 109-43. [7] C Annavarapu, M Hautefeuille, JE Dolbow. Stable imposition of stiff constraints in explicit dynamics for embedded finite element methods. Int J Numer Meth Eng 2012; 92: 206-28. [8] J Mergheim. A variational multiscale method to model crack propagation at finite strains. Int J Numer Meth Eng 2009; 80: 269-89. [9] D Ling, Q Yang, B Cox. An augmented finite element method for modeling arbitrary discontinuities in composite materials. Int J Fracture 2009; 156: 53-73. 18

[10] J Zhi, B Chen, T Tay. Geometrically nonlinear analysis of matrix cracking and delamination in composites with floating node method. Comput Mech 2018. [11] NA Abdullah, JL Curiel-Sosa, ZA Taylor, B Tafazzolimoghaddam, JL Martinez Vicente, C Zhang. Transversal crack and delamination of laminates using XFEM. Compos Struct 2017; 173: 78-85. [12] XF Hu, BY Chen, M Tirvaudey, V Tan, TE Tay. Integrated XFEM-CE analysis of delamination migration in multi-directional composite laminates. Composites Part A: Applied Science and Manufacturing 2016; 90: 161-73. [13] N Sukumar, DL Chopp, N Moës, T Belytschko. Modeling holes and inclusions by level sets in the extended finite-element method. Comput Method Appl M 2001; 190: 6183-200. [14] X Li, J Chen. An extended cohesive damage model for simulating multicrack propagation in fibre composites. Compos Struct 2016; 143: 1-8. [15] P Prabhakar, AM Waas. A novel continuum-decohesive finite element for modeling in-plane fracture in fiber reinforced composites. Compos Sci Technol 2013; 83: 1-10. [16] M Heidari-Rarani, M Sayedain. Finite element modeling strategies for 2D and 3D delamination propagation in composite DCB specimens using VCCT, CZM and XFEM approaches. Theor Appl Fract Mec 2019; 103: 102246. [17] M Heidari-Rarani, MM Shokrieh, PP Camanho. Finite element modeling of mode I delamination growth in laminated DCB specimens with R-curve effects. Composites Part B: Engineering 2013; 45: 897-903. [18] SM Jensen, MJ Martos, BLV Bak, E Lindgaard. Formulation of a mixed-mode multilinear cohesive zone law in an interface finite element for modelling delamination with R-curve effects. Compos Struct 2019; 216: 477-86. [19] A Turon, PP Camanho, J Costa, CG Dávila. A damage model for the simulation of delamination in advanced composites under variable-mode loading. Mech Mater 2006; 38: 1072-89. [20] E Lindgaard, BLV Bak, JA Glud, J Sjølund, ET Christensen. A user programmed cohesive zone finite element for ANSYS Mechanical. Eng Fract Mech 2017; 180: 229-39. [21] G Alfano. On the influence of the shape of the interface law on the application of cohesive-zone models. Compos Sci Technol 2006; 66: 723-30. [22] L Zhao, Y Gong, J Zhang, Y Chen, B Fei. Simulation of delamination growth in multidirectional laminates under mode I and mixed mode I/II loadings using cohesive elements. Compos Struct 2014; 116: 509-22. [23] L Zhao, Y Wang, J Zhang, Y Gong, N Hu, N Li. XFEM-based model for simulating zigzag delamination growth in laminated composites under mode I loading. Compos Struct 2017; 160: 1155-62. [24] Y Gong, B Zhang, SR Hallett. Delamination migration in multidirectional composite laminates under mode I quasi-static and fatigue loading. Compos Struct 2018; 189: 160-76. [25] LP Canal, G Pappas, J Botsis. Large scale fiber bridging in mode I intralaminar fracture. An embedded cell approach. Compos Sci Technol 2016; 126: 52-9. [26] Y Gong, B Zhang, L Zhao, J Zhang, N Hu, C Zhang. R-curve behaviour of the mixed-mode I/II delamination in carbon/epoxy laminates with unidirectional and multidirectional interfaces. Compos Struct 2019; 223: 110949. [27] L Zhao, Y Wang, J Zhang, Y Gong, Z Lu, N Hu, J Xu. An interface-dependent model of plateau fracture toughness in multidirectional CFRP laminates under mode I loading. Composites Part B: Engineering 2017; 131: 196-208. [28] A Airoldi, CG Dávila. Identification of material parameters for modelling delamination in the presence of fibre bridging. Compos Struct 2012; 94: 3240-9. [29] A Ghadirdokht, M Heidari-Rarani. Delamination R-curve behavior of curved composite laminates. Composites Part B: Engineering 2019; 175: 107139. [30] MM Shokrieh, M Salamat-talab, M Heidari-Rarani. Effect of interface fiber angle on the R-curve behavior of E-glass/epoxy DCB specimens. Theor Appl Fract Mec 2016. 19

[31] Y Gong, L Zhao, J Zhang, N Hu, C Zhang. An insight into three approaches for determining fatigue delamination resistance in DCB tests on composite laminates. Compos Part B-Eng 2019; 176: 107206. [32] BF Sørensen, TK Jacobsen. Large-scale bridging in composites: R-curves and bridging laws. Composites Part A: Applied Science and Manufacturing 1998; 29: 1443-51. [33] G Frossard, J Cugnoni, T Gmür, J Botsis. Mode I interlaminar fracture of carbon epoxy laminates: Effects of ply thickness. Composites Part A: Applied Science and Manufacturing 2016; 91: 1-8. [34] Z Suo, G Bao, B Fan. Delamination R-curve phenomena due to damage. J Mech Phys Solids 1992; 40: 1-16. [35] SM Spearing, AG Evans. The role of fiber bridging in the delamination resistance of fiber-reinforced composites. Acta metallurgica et materialia 1992; 40: 2191-9. [36] Y Gong, L Zhao, J Zhang, Y Wang, N Hu. Delamination propagation criterion including the effect of fiber bridging for mixed-mode I/II delamination in CFRP multidirectional laminates. Compos Sci Technol 2017; 151: 302-9. [37] Y Gong, L Zhao, J Zhang, N Hu. An improved power law criterion for the delamination propagation with the effect of large-scale fiber bridging in composite multidirectional laminates. Compos Struct 2018; 184: 961-8. [38] L Sørensen, J Botsis, T Gmür, L Humbert. Bridging tractions in mode I delamination: measurements and simulations. Compos Sci Technol 2008; 68: 2350-8. [39] V Tamuzs, S Tarasovs, U Vilks. Progressive delamination and fiber bridging modeling in double cantilever beam composite specimens. Eng Fract Mech 2001; 68: 513-25. [40] BF Sørensen, TK Jacobsen. Determination of cohesive laws by the J integral approach. Eng Fract Mech 2003; 70: 1841-58. [41] LP Canal, M Alfano, J Botsis. A multi-scale based cohesive zone model for the analysis of thickness scaling effect in fiber bridging. Compos Sci Technol 2017; 139: 90-8. [42] E Farmand-Ashtiani, J Cugnoni, J Botsis. Specimen thickness dependence of large scale fiber bridging in mode I interlaminar fracture of carbon epoxy composite. Int J Solids Struct 2015; 55: 58-65. [43] Z Daneshjoo, MM Shokrieh, M Fakoor. A micromechanical model for prediction of mixed mode I/II delamination of laminated composites considering fiber bridging effects. Theor Appl Fract Mec 2018; 94: 46-56. [44] R Gutkin, ML Laffan, ST Pinho, P Robinson, PT Curtis. Modelling the R-curve effect and its specimen-dependence. Int J Solids Struct 2011; 48: 1767-77. [45] S Giannis, K Hansen, RH Martin.Accounting for the R-curve effects on the mode I fatigue delamination growth characterisation of unidirectional composites. Proceedings of the American Society for composites2010. p. 20-2. [46] CG Dávila, CA Rose, PP Camanho. A procedure for superposing linear cohesive laws to represent multiple damage mechanisms in the fracture of composites. Int J Fracture 2009; 158: 211-23. [47] M Heidari-Rarani, AR Ghasemi. Appropriate shape of cohesive zone model for delamination propagation in ENF specimens with R-curve effects. Theor Appl Fract Mec 2017; 90: 174-81. [48] Z Daneshjoo, MM Shokrieh, M Fakoor. A micromechanical model for prediction of mixed mode I/II delamination of laminated composites considering fiber bridging effects. Theor Appl Fract Mec 2018; 94: 46-56. [49] MM Shokrieh, M Heidari-Rarani, MR Ayatollahi. Calculation of GI for a multidirectional composite double cantilever beam on two-parametric elastic foundation. Aerosp Sci Technol 2011; 15: 534-43. [50] J Zhang, L Peng, L Zhao, B Fei. Fatigue delamination growth rates and thresholds of composite laminates under mixed mode loading. Int J Fatigue 2012; 40: 7-15. [51] BD Davidson, R Kruger, M König. Three dimensional analysis and resulting design 20

recommendations for unidirectional and multidirectional end-notched flexure tests. J Compos Mater 1995; 29: 2108-33. [52] ASTM. Standard Test Method for Mode I Interlaminar Fracture Toughness of Unidirectional Fiber-Reinforced Polymer Matrix Composites. ASTM International, West Conshohocken PA, 2007, DOI: 10.1520/D5528_D5528M. [53] MM Shokrieh, M Heidari-Rarani, MR Ayatollahi. Delamination R-curve as a material property of unidirectional glass/epoxy composites. Mater Design 2012; 34: 211-8. [54] L Peng, J Zhang, L Zhao, R Bao, H Yang, B Fei. Mode I delamination growth of multidirectional composite laminates under fatigue loading. J Compos Mater 2011; 45: 1077-90. [55] L Yao, RC Alderliesten, M Zhao, R Benedictus. Discussion on the use of the strain energy release rate for fatigue delamination characterization. Composites Part A: Applied Science and Manufacturing 2014; 66: 65-72. [56] Q Ye, P Chen. Prediction of the cohesive strength for numerically simulating composite delamination via CZM-based FEM. Composites Part B: Engineering 2011; 42: 1076-83. [57] MM Shokrieh, M Salamat-talab, M Heidari-Rarani. Dependency of bridging traction of DCB composite specimen on interface fiber angle. Theor Appl Fract Mec 2017; 90: 22-32. [58] S Stutz, J Cugnoni, J Botsis. Studies of mode I delamination in monotonic and fatigue loading using FBG wavelength multiplexing and numerical analysis. Compos Sci Technol 2011; 71: 443-9. [59] S Stutz, J Cugnoni, J Botsis. Crack-fiber sensor interaction and characterization of the bridging tractions in mode I delamination. Eng Fract Mech 2011; 78: 890-900. [60] MM Shokrieh, M Salamat-talab, M Heidari-Rarani. Effect of initial crack length on the measured bridging law of unidirectional E-glass/epoxy double cantilever beam specimens. Mater Design 2014; 55: 605-11. [61] N Chandra, H Li, C Shet, H Ghonem. Some issues in the application of cohesive zone models for metal-ceramic interfaces. Int J Solids Struct 2002; 39: 2827-55. [62] DM Montenegro, G Pappas, J Botsis, M Zogg, K Wegener. A comparative study of mode I delamination behavior of unidirectional glass fiber-reinforced polymers with epoxy and polyurethane matrices using two methods. Eng Fract Mech 2019; 206: 485-500. Figure Captions Fig. 1. Sketch of the traditional three-linear CZM in literatures [17]. Fig. 2. (a) Sketch of the new three-linear CZM, (b) interfacial damage and (c) peeled-off bridging fibers. Fig. 3. Experimental results: (a) R-curves and (b) Typical load versus displacement curves. Fig. 4. Delamination path pictures of tested specimens with (a) 0°/5° interface and (b) 45°/-45° interface. Fig. 5. The FE model of DCB tests. Fig. 6. The relation between fracture toughness and crack opening displacement at the pre-crack tip. Fig. 7. The bridging stress distribution in tested specimens. Fig. 8. Numerical results from the traditional bilinear CZM model for the (a) 0°/5° interface and (b) 21

45°/-45° interface. Fig. 9. Numerical results from the traditional three-linear CZM model for the (a) 0°/5° interface and (b) 45°/-45° interface. Fig. 10. Numerical results from the new CZM model for the (a) 0°/5° interface and (b) 45°/-45° interface. Fig. 11. Sensitivity analysis of the interfacial strength on the simulated results for the (a) 0°/5° interface and (b) 45°/-45° interface.

22

Fig. 1. Sketch of the traditional three-linear CZM in literatures [17].

Fig. 2. (a) Sketch of the new three-linear CZM, (b) interfacial damage and (c) peeled-off bridging fibers.

23

Fig. 3. Experimental results: (a) R-curves and (b) Typical load versus displacement curves.

Fig. 4. Delamination path pictures of tested specimens with (a) 0°/5° interface and (b) 45°/-45° interface.

24

Fig. 5. The FE model of DCB tests.

Fig. 6. The relation between fracture toughness and crack opening displacement at the pre-crack tip.

25

Fig. 7. The bridging stress distribution in tested specimens.

Fig. 8. Numerical results from the traditional bilinear CZM model for the (a) 0°/5° interface and (b) 45°/-45° interface.

26

Fig. 9. Numerical results from the traditional three-linear CZM model for the (a) 0°/5° interface and (b) 45°/-45° interface.

Fig. 10. Numerical results from the new CZM model for the (a) 0°/5° interface and (b) 45°/-45° interface.

27

Fig. 11. Sensitivity analysis of the interfacial strength on the simulated results for the (a) 0°/5° interface and (b) 45°/-45° interface.

28

Table Captions Table 1 Detailed values of the fracture toughness, the length of fiber bridging zone and other input parameters in FE models. Table 2 Ultimate loads for both interfaces under different values of the interfacial strength. Table 1 Detailed values of the fracture toughness, the length of fiber bridging zone and other input parameters in FE models. Ply sequence

Interface

016//(+5/-5/06)S 0°/5° (+45/-45/06)S//(-45/+45/06)S 45°/-45°

Gini (J/m2) 100

Gprop (J/m2) 825 890

lcz (mm) 50 26

K1 (N/m3)

μ

1×1015 1×10-5

σ0 σb (MPa) (MPa) 1.0 60 2.1

Table 2 Ultimate loads for both interfaces under different values of the interfacial strength. σb (MPa) Experimental ultimate loads (N) 0°/5° Numerical ultimate loads (N) Relative errors (%) Experimental ultimate loads (N) 45°/-45° Numerical ultimate loads (N) Relative errors (%) Note: S.D. is the standard deviation.

50 60* 104.2(S.D.=2.6) 100.4 105.9 -3.6 1.6 91.3 (S.D.=0.4) 83.6 82.7 -8.4 -9.4

Interface

70

80

102.9 -1.2

104.9 0.7

84.7 -7.2

85.8 -6.0

Yu Gong: Methodology, Formal Analysis, Writing – Original Draft, Writing – Review & Editing. Yixin Hou: Methodology, Software, Formal Analysis, Writing – Original Draft. Libin Zhao: Conceptualization, Methodology, Formal Analysis, Writing – Original Draft, Writing – Review & Editing, Supervision, Project Administration. Wangchang Li: Investigation, Resources. Jianyu Zhang: Validation, Data Curation, Resources, Writing – Review & Editing, Project Administration. 29

Ning Hu: Resources, Writing – Review & Editing.

30

Graphical abstract

31