Simulation of Mode I Interlaminar Damage of a GFRP Composite Using Cohesive Laws in the Framework of the Equivalent LEFM R-Curve and an Optimised Algorithm.

This paper is focused on mode I delimitation of a unidirectional glass fibre reinforced polymer (GFRP) composite. The aim is to propose an accurate and simple characterisation of three cohesive zone models (CZM)-bilinear, trilinear, and potential-from the measurement of the load-displacement curve during a double cantilever beam experimental test. For that, a framework based on the equivalent linear elastic fracture mechanics (LEFM) R-curve is here proposed, which has never before been developed for a bilinear and a potential CZM. Besides, in order to validate this strategy, an optimisation algorithm for solving an inverse problem is also implemented. It is shown that the parameters' identification using the equivalent LEFM R-curve enables the same accuracy but reduces 72% the numerical efforts respect to a "blind fitting" (i.e., the optimisation algorithm). Therefore, even if optimisation techniques become popular at present due to their easy numerical implementation, strategies founded on physical models are still better solutions especially when evaluating the objective function is expensive as in mechanical problems.

1. Introduction

The use of structural composites in high performance applications, such as those present in the automotive or aerospace industry, has been steadily growing in during the last 50 years [1]. Despite the progress made on areas related to design, manufacturing and analysis of composites, one of the challenges that still remains is the accurate pren class="Disease">diction of the progressive failure of the composite due to delamination [2]. These failure mechanisms are intrinsically related to the hierarchical nature of the laminates. Micro (ply), meso (laminate), and macro (structural) length scales should be considered to correctly predict the mechanical response of the composite. In this context, multi-scale modelling approaches that inherently incorporate the different scales present in composites, are a natural and frequently used alternative to study their mechanical response [3,4]. Furthermore, the increase in the use of numerical techniques, such as multi-scale modelling strategies, observed in the recent years, has led to a boost on the use of virtual testing techniques in the design and optimisation process of new laminate materials [5,6]. Additional benefits of virtual testing techniques are the cost and time reductions in the design process compared to experimental testing [6].

Failure mechanisms in composite laminates can be categorised as intraply (e.g., fibre n class="Disease">fracture or matrix cracking) and interply (i.e., delamination) and, of course, complex interactions between them also can take place [7,8]. Moreover, it has been shown that the type of failure that is observed depends on the specimen size [9]. Delamination, defined as the crack propagation between adjacent plies, is one of the most relevant damage scenarios because it drastically reduces the mechanical strength and may lead to structural collapse [10,11,12,13]. The delamination fracture toughness is well established through the strain energy release rate , according respectively to the three fracture modes (, and ) and their combinations. The experimental procedures used to determine the delamination fracture toughness for each fracture mode are established in several standards [14,15,16,17]. Here, details on the specimens dimensions, loading protocol and data reduction are clearly outlined. The delamination occurs when the energy dissipated during fracture per unit of newly created surface is greater or equal than a critical strain energy release rate (i.e., the resistance to crack growth), which can be viewed as a material property. Taking into account a general body with constant thickness B and an initial crack length , under a loading P (N)ormal to the crack plane, the linear elastic fracture mechanics (LEFM) theory enables evaluating the fracture energy as follows [18]: where is the compliance or displacement to applied load P ratio.

The double cantilever beam (pan class="Chemical">DCB) test shown in Figure 1a and used in this work is the most widely used procedure for the measurement of mode I pan class="Disease">delamination fracture toughness. For the load-displacement curve (see Figure 1b) and considering a rigid foundation at the pan class="Chemical">crack ending, the compliance can be given by the beam theory as , where and E is the Young’s module. A corrected compliance taking into account an elastic foundation is presented in [19]. Finally, the relation linking P and during the propagation phase can be obtained from Equation (1). Further details can be found in [20].

Figure 1

Double cantilever beam test.

The R-curve shown in Figure 1c illustrates the variation of the material n class="Chemical">crack resistance with respect to the crack propagation length . Ideal brittle materials present a flat R-curve, as shown by the blue line in Figure 1c. In this case, and when the crack propagates, the energy release rate remains constant and equal to (for the sake of simplicity subindex is omitted in the following). This behaviour is also observed when the crack propagation process of a material is studied by means of the LEFM method. Quasibrittle materials have a rising R-curve with an initiation phase where the resistance to the crack growth increases. Then the critical strain energy release rate reaches a steady-state plateau for a critical crack extension denoted as , i.e., crack propagates in a self-similar steady way [21]. Rising R-curves require nonlinear fracture theories to describe the existence of a fracture process zone (FPZ) of length ahead of the crack tip (see Figure 2). The FPZ is where inelastic crack propagation mechanisms, such as fibre bridging and microcracks, take place. In fact, materials presenting large scale bridging have R-curves strongly dependent on the specimen’s geometry and therefore their constitutive damage model cannot be regarded as a material property [22,23]. More recently, the transition from 1D standard tests to 2D delamination scenarios has shown higher value of the fracture toughness for plates—due to stretching mechanisms affecting their stiffness—compared to the DCB specimens [24]. For all these quasibrittle behaviours, the R-curves need complex experimental setups for measuring the crack length, e.g., traveling microscope, crack gauge or video cameras. Another less expensive option is to use the equivalent LEFM [25], where the increase of the compliance can be related to the propagation of an equivalent LEFM crack. For any point of the experimental load-displacement curve in Figure 1b, a secant compliance is associated with each load; the corresponding equivalent crack extension is determined by solving the equation for and the crack growth resistance is then determined from Equation (1) [26].

Figure 2

Fracture process zone and equivalent LEFM: is defined as the distance between the tip of the stress-free crack and the point along the potential crack path where damage begins (schema inspired from [26]).

From a numerical point of view, there are a few techniques that can be used within a finite element (FE) method framework to address the delamination process in composites. Among the most frequently used strategies, there are adaptive remeshing procedures such as the virtual n class="Chemical">crack closure technique (VCCT) [27,28,29]; the use of enrichment functions near the crack tip such as the extended finite element method (X-FEM) [30,31,32,33]; or zero thickness interfaces with a continuum damage mechanics model such as the cohesive zone models (CZM) [34,35,36]. Advantages, limitations, and challenges of these three families of methods are discussed in [37]. If the crack path is known a priori, as in delamination of composite laminates, the CZM is the simplest and most accurate method. It has been widely used by researchers in recent decades for predicting both crack nucleation and propagation in composites [13,25,26,36,38,39,40,41,42,43,44,45,46,47,48,49,50,51]. The CZM method is defined by a constitutive law or softening function that relates the cohesive interface transfer traction f to the displacement jump w. The softening function is written in terms of un damage variable d, and characterised by a positive high initial stiffness and a maximum critical traction level . Upon reaching , the softening function is described by a negative tangent stiffness until a critical displacement jump is achieved. At this point, the system presents no more load-bearing capacity, i.e., . As depicted in Figure 3, different forms of softening laws have been introduced such as linear, bilinear, trilinear, trapezoidal or exponential [52]. The area under the entire traction-displacement jump curve is the fracture energy , i.e., the total energy required to completely separate the interface per unit area. These models are chosen according to a compromise between simple identification of its parameters and an accurate prediction of the crack propagation. In fact, due to the incapability to directly measure the curve, especially for materials with non-negligible FPZ, the model characterisation usually combines experimental data, theory (LEFM or J-Integral) or FE simulations. For example, it is possible to embed a fibre Bragg grating (FBG) sensor close to the crack tip and to measure the distributed strains [53,54,55], then an inverse method relates strains to the traction-separation curve. Digital image correlation (DIC) has also allowed inverse procedures combining full field cinematic data and FE simulations [56,57]. The approaches based on J-integral need experimental data such as the crack length, applied load or crack tip opening displacement (CTOD) [58,59,60,61]. In [62] a J-integral procedure is compared to an inverse optimisation scheme that minimises the difference between the experimental and simulated strains along the specimen. Although these techniques require very high resolution equipment to capture the CTOD or FPZ, which can cost a lot of money and be difficult to implement in specimens tested in a controlled environmental. Another option is to use inverse methods for minimising the residual between experimental and numerical curves [44,63,64,65,66]. However, despite their simplicity, they are very time-consuming because it is necessary several virtual tests to evaluate different values of each parameter of the softening function . Moreover, FE simulations can be very expensive if models are more accurate, such as 6–8 CPU hours for only one 3D DCB [67]. In order to be more efficient, an inverse method combined with a model based on a Dugdale’s condition [68] or closed-form analytical solutions [69] has been developed for identifying multilinear, piecewise constant or bilinear CZM. Nevertheless, for more sophisticated shapes of CZM, optimisation algorithms seem to be the only possible way to identify the parameters. In this work, another alternative for reducing numerical and experimental efforts is exploited, which is based on the equivalent LEFM R-curve and has been initially proposed for a trilinear CZM [26]. To the best of the authors’ knowledge, this methodology has never before been developed for identifying a bilinear and a potential CZM or compared to “blind fitting” (i.e., an optimisation algorithm).

Figure 3

Traction-separation laws: (a) trilinear; (b) bilinear; and (c) potential CZM.

2. Experimental Test

n class="Chemical">DCB specimens were manufactured through a vacuum infusion process, considering a fibre/resin weight ratio of and embedding a thin film at the mid-plane of the specimen for the pre-crack. Geometry (see Figure 1a) is specified in Table 1; an Epoxi 713 resin matrix with an E1174 hardener from the Chilean commical">pany Fibratec (Santiago, Chile) was employed, while the reinforcement is a unidirectional glass fibre from the German commical">pany P-D Interglas Technologies GmbH (now acquired by Porcher Industries, Eclose Badinieres, France), currently named as UD 220 g/m (i.e., with 207 g/m in the warp direction and 13 g/m in the fill direction). Then, the composite, previously characterised in [70], has the following elastic properties: GPa, GPa, (-), (-) and GPa, where direction is in the fibre direction and aligned through the specimen’s length.

Table 1

DCB specimen geometry according to [17].

Thickness h (mm)	Width B (mm)	Length D (mm)	Pre-Crack a_0 (mm)
2.72 ± 0.07	20.4 ± 0.08	124.68 ± 0.52	47

The pan class="Chemical">DCB test was carried out taking into account the ISO 15024 standard [17], under quasi-static conditions using a testing machine ZwickRoell (Ulm, Germany) provided with a 5 [kN] load cell and the testXpert testing software (see Figure 4). The load-displacement and resistance curves obtained from five experiments are drawn in Figure 5, whereas Table 2 summaries the critical energy release rate and the maximal applied load together with their standard deviations.

Figure 4

Experimental DCB tests: layout according to [17].

Figure 5

Experimental DCB tests: (a) load-displacement curve; (b) resistance curve.

Table 2

Experimental DCB tests: critical energy release rate and maximal applied load.

Specimen	Gcexp (N/m)	Pmaxexp (N)
1	981.65	31.08
2	926.60	32.16
3	1062.80	29.72
4	964.90	27.98
5	974.8	30.70
average value	982.15	30.33
standard deviation	49.85	1.58

3. Cohesive Zone Models

Traction-separation laws can be written in terms of a damage variable d, ranging from 0 to 1, for a healthy to a completely damaged interface point, respectively. It is important to notice that in the tridimensional case is a vector-valued function, while the initial stiffness is a second order tensor. However, because this paper is only related to the mode I delamination, variables are all scalars. Therefore, the initial stiffness is progressively weaken according to: where the symbols distinguish between the positive and negative part of the normal displacement in order to take into account the difference between tensile and compression. Considering the irreversibility of damage, d depends on the whole load history, i.e., .

In this work, a bilinear [71,72], a trilinear [38,73] and a potential model [40] are studied, which are schematized in Figure 3. In all these laws, the pan class="Disease">fracture energy corresponds to the area under the curve. However, the energy is decomposed into two parts () for the trilinear CZM, which has been attributed to micro-pan class="Chemical">cracking () and fiber-bridging () [73]. For each model, damage d is written in terms of the displacement jump w and the parameters to be identified, as summarized in Table 3.

Table 3

Damage functions and parameters of the bilinear, trilinear and potential CZM.

CZM	Trilinear [26]	Bilinear [72]	Potential [40]
	d=0, w<w_0	d=0, w≤w_0	d=minnn+1YG_c^n,1
	d=w_b(w−w_o)(1−γ)w(w_b−w_o),w_o≤w≤w_b	d=w_cw_c−w_0w−w_0w, w_0<w≤w_c	where Y=12K_0w^2
damage law	d=1−γw_b(w_c−w)w(w_c−w_b),w_b≤w≤w_c	d=1, w>w_c
	where γ=f_bw_of_tw_b,w_0=f_tK_0,	where w_0=f_tK_0,f_t=2G_cw_c
	w_b=2G_fμft,f_b=2G_fbw_c
parametersto be identified	w_c, G_fμ/G_c, f_t	w_c, K_0	K_0, n

4. CZM Characterization Using the Equivalent LEFM R-Curve

As explained in Section 1, the equivalent LEFM R-curve enables representing the influence of the FPZ development on the specimen compliance, through an elastically equivalent n class="Chemical">crack a (see Figure 2) located at some distance ahead of the initial crack (or the current stress-free crack ). In [26], a new procedure to identify the parameters of the trilinear cohesive model has been proposed. It is based on the equivalent LEFM R-curve and on a dimensional analysis in order to relate the parameters’ dependency on the geometry and material of the specimen. The main idea is to relate each model parameter with the load-displacement curve and the corresponding equivalent LEFM R-curve through numerical simulations. Among the advantages, this methodology avoids an experimental measure of the critical opening and determines the CZM parameters more efficiently than blind fitting or optimization methods such as genetic algorithms, because less numerical evaluations are needed.

In this work, FE simulations of n class="Chemical">DCB tests are performed using 1D Euler-Bernoulli beam elements, coded through an OCTAVE routine, where only one half of the specimen—due to symmetry—is modelled and discretised into 310 finite elements. Geometry and elastic properties of the specimen are according to Section 2 (due to the 1D model, only GPa is taken into account), whereas the critical energy release rate is given by the experimental test ( N/m). In the following, the procedure defined for the trilinear CZM in [26] is applied to the fracture test of Section 3, then the methodology is extended to the bilinear and potential laws.

4.1. Trilinear CZM

The impact on the equivalent LEFM R-curve of each cohesive parameter is first analysed: the critical opening , the distribution of the critical energy release rate and the tensile strength . The critical pan class="Chemical">crack extension will be then associated with these parameters and the specimen size through a dimensional analysis. In this analysis mm and n>n class="Gene">N/m are kept constant.

Influence of the critical opening: numerical simulations are carried out affecting the critical opening mm but fixing the values (-) and Pa. From Figure 6, it is observed that does not have an influence before reaching the 50% of . After that, the critical pan class="Chemical">crack extension and the critical displacement jump have a positive correlation, which means that the length of the pan class="Disease">fracture process zone increases when increases.

Figure 6

Trilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the critical opening with Pa and (-) as constants.

Influence of the pan class="Disease">fracture energy distribution: keeping constant Pa and mm, Figure 7 shows that at varying the critical energy release rate (-), the load-displacement plot and the equivalent LEFM R-curve are invariable if the dissipated energies are lower than , respectively. If the ratio tends to one, the R-curve look like a one of a brittle material and the maximal load on the load-displacement curve increases. The critical n>n class="Chemical">crack extension always remains the same.

Figure 7

Trilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the fracture energy distribution with mm and 1.15 Pa as constants.

Influence of the tensile strength: Figure 8 shows the impact at varying [ Pa] whereas (-) and mm are unchanged. It can be concluded that the n class="Disease">tensile strength impacts the response at the beginning of both curves: if increases, the behaviour becomes a brittle one. When the dissipated energy reaches the value , the n>n class="Disease">fracture response starts to be the same for any value of and the critical crack extension becomes identical.

Figure 8

Trilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the tensile strength with mm and (-) as constants.

Trilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the pan class="Disease">fracture energy distribution with mm and 1.15 Pa as constants.

From the studies carried out in Figure 6, Figure 7 and Figure 8, the influence of each cohesive parameter when other variables are constant has been known. Now, from a dimensional analysis it can be supposed that depends on the specimen size and the cohesive parameters, as proposed in [26,74]: where D is a characteristic dimension of the specimen, and are the Hillerborg’s characteristic length and a characteristic pan class="Chemical">crack opening, respectively:

Because it was observed that the ratio does not have influence on the critical pan class="Chemical">crack extension (see Figure 7b), the third argument of Equation (3) is able to be directly vanished. According to Figure 8b, neither depends on , therefore if Equation (3) is homogeneous in —i.e., depends on —it can be cancelled when factoring by [26]: where , furthermore the right equation has been multiplied by . On the other hand, a lower bound of the critical pan class="Chemical">crack extension has been previously studied in [74] for , given by:

Then, it is expected that .

The nonlinear expression relating and , Equation (5), has the ability to be solved for through the equivalent LEFM R-curve and FE computations, according to: where when .

The critical opening can be now found employing Equation (7), while the other cohesive parameter—the tensile strength —can also be studied through a dimensional analysis. Therefore, the pan class="Chemical">crack length , for a given energy release rate , is allowed to be obtained from the following general expression:

From virtual tests previously carried out for different values of and , but keeping constant, it is observed that the equivalent R-curve remains almost the same when (see Figure 6b and Figure 7b). In fact, only has an influence at the beginning of the R-curve if and are both fixed (see Figure 8b). Therefore, function is able to be considered independent of and if . For instance, calculations are here proposed with -:

As proposed in [26], an expression relating and is permitted to be obtained at multiplying Equation (9) by , then: where can be determined from Figure 8b, considering the horizontal line for and it verifies . Finally, the dimensionless functions and plotted in Figure 9 are obtained through a Hermite spline cubic interpolation.

Figure 9

Trilinear cohesive model. The dimensionless functions for model parameters as a function of the crack length: (a) function; (b) function.

4.2. Bilinear CZM

The previous methodology is here developed for a bilinear model. In this case, it has been studied the influence of the critical opening and the initial stiffness , while is always kept constant and equals to .

Influence of the critical opening: has an inverse correlation with the maximal applied load P when considering mm and keeping unchanged pan class="Gene">N/m, as observed in Figure 10a. Moreover, Figure 10b shows that the critical pan class="Chemical">crack extension is inversely proportional to the tensile strength , i.e., the interface becomes more brittle for higher values of .

Figure 10

Bilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the critical opening with N/m as constant.

Influence of the initial stiffness : from Figure 11b it is observed that pan class="Gene">N/m does not have any influence on the critical n>n class="Chemical">crack extension if is kept constant. only has an effect on the beginning of the equivalent R-curve, which is not significant on the load-displacement curve (see Figure 11a).

Figure 11

Bilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the initial stiffness with mm as constant.

Bilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the critical opening with pan class="Gene">N/m as constant.

As previously introduced in Section 4.1, the critical pan class="Chemical">crack extension can be expressed as a function which depends on the specimen size and on the cohesive parameters as: where is a characteristic stiffness defined as follows:

From Figure 11b it is assumed that does not have influence on the critical pan class="Chemical">crack extension , then the third argument disappears. Besides, the relation enables going without having the second variable , then multiplying by and considering Equation (4): where it is expected, when , i.e., in concordance with Equation (6). Finally, the critical opening is computed in the same way as for the trilinear law: where when .

To characterise the stiffness parameter , a dimensionless general expression for the pan class="Chemical">crack extension is able to be written as:

From Figure 10b it should be noticed that the effect of on the R-curve depends on the selection of . First, it is assumed that has been chosen using Equation (14); secondly, the pan class="Chemical">crack extension is searched for a given G where the impact of is significative (e.g., -), subsequently last equation becomes:

The dependence on (i.e., on ) has the chance to be vanished if factoring by , as follows:

Then, multiplying last this expression by :

Finally, the initial stiffness can be stablished employing the next expression: where it is verified that . The dimensionless functions, and , characterizing both cohesive parameters are obtained through a Hermite spline cubic interpolation and plotted in Figure 12.

Figure 12

Bilinear cohesive model. The dimensionless functions for model parameters as a function of the crack length: (a) function; (b) function.

4.3. Potential CZM

The parameters enabling identify the potential CZM are the initial stiffness and the dimensionless variable n. However, it is not possible to separately relate them to the critical pan class="Chemical">crack extension . Actually, for a fixed n, influences the whole R-curve (including ), as shown in Figure 13. Same behaviour takes place for different values of n but a fixed , as demonstrated in Figure 14. In both cases is always kept constant and equals to .

Figure 13

Potential cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the initial stiffness with (-) as constant.

Figure 14

Potential cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the parameter n with N/m as constant.

Despite this impossibility, it is feasible to find a connection between (n, ) and if paying attention to the following relation from Table 3:

More precisely, when considering , the critical displacement jump is able to be written in terms of n and as follows:

Furthermore, it is possible to obtain an expression for through a stationary point of Equation (2), then:

When examining Equation (22) for several values of n and , as plotted in Figure 15 and detailed in Table A1 in Appendix A, it is noticed that is essentially unchanged if the term (or ) remains invariable (series A and D). Nevertheless, if one parameter (n or ) is kept constant and the other is varied, increases when or n increases, respectively (series B and C). In the following, the effect of on the equivalent R-curve is studied, while is equals to .

Figure 15

Potential cohesive model. Influence of (a) n and (b) on (more details in Table A1).

Table A1

Potential cohesive model. Influence of model parameters on and .

	K_0 (N/m^3)	n (-)	n+1nK_0 (m^3/N)	w_c (m)	f_t (Pa)
	1 ·10^14	0.01	1 ·10^−12	4.43 ·10^−5	3.26 ·10^7
	5 ·10^13	0.02	1 ·10^−12	4.43 ·10^−5	3.26 ·10^7
	1 ·10^13	0.11	1 ·10^−12	4.43 ·10^−5	3.27 ·10^7
serie A	5 ·10^12	0.25	1 ·10^−12	4.43 ·10^−5	3.28 ·10^7
	3 ·10^12	0.5	1 ·10^−12	4.43 ·10^−5	3.32 ·10^7
	2 ·10^12	1	1 ·10^−12	4.43 ·10^−5	3.41 ·10^7
	1.5·10^12	2	1 ·10^−12	4.43 ·10^−5	3.56 ·10^7
	1.3·10^12	3.3	1 ·10^−12	4.43 ·10^−5	3.69 ·10^7
	1 ·10^13	2	1.5·10^−13	1.72 ·10^−5	9.18 ·10^7
	1 ·10^13	1	2·10^−13	1.98 ·10^−5	7.63 ·10^7
	1 ·10^13	0.5	3·10^−13	2.43 ·10^−5	6.07 ·10^7
serie B	1 ·10^13	0.1	1.1·10^−12	4.65 ·10^−5	3.11 ·10^7
	1 ·10^13	0.01	1·10^−11	1.41 ·10^−5	1.03 ·10^7
	1 ·10^13	0.001	1·10^−10	4.43 ·10^−4	3.26 ·10^6
	1 ·10^13	0.0001	1·10^−9	1.4 ·10^−3	1.03 ·10^6
	1 ·10^13	0.00001	1·10^−8	4.43 ·10^−3	3.26 ·10^5
	1 ·10^14	0.5	3 ·10^−14	7.68 ·10^−6	1.92 ·10^8
	5 ·10^13	0.5	6 ·10^−14	1.09 ·10^−5	1.36 ·10^8
	1 ·10^13	0.5	3 ·10^−13	2.43 ·10^−5	6.07 ·10^7
serie C	5 ·10^12	0.5	6 ·10^−13	3.43 ·10^−5	4.29 ·10^7
	3 ·10^12	0.5	1 ·10^−12	4.43 ·10^−5	3.32 ·10^7
	1 ·10^12	0.5	3 ·10^−12	7.68 ·10^−5	1.92 ·10^7
	5 ·10^11	0.5	6 ·10^−12	1.09 ·10^−4	1.36 ·10^7
	1 ·10^11	0.5	3 ·10^−11	2.43 ·10^−4	6.07 ·10^6
	1 ·10^14	0.001	1 ·10^−11	1.40 ·10^−4	1.03 ·10^7
	5 ·10^13	0.002	1 ·10^−11	1.40 ·10^−4	1.03 ·10^7
	1 ·10^13	0.0101	1 ·10^−11	1.40 ·10^−4	1.03 ·10^7
serie D	5 ·10^12	0.0204	1 ·10^−11	1.40 ·10^−4	1.03 ·10^7
	1 ·10^12	0.111	1 ·10^−11	1.40 ·10^−4	1.03 ·10^7
	5 ·10^11	0.25	1 ·10^−11	1.40 ·10^−4	1.04 ·10^7
	3 ·10^11	0.5	1 ·10^−11	1.40 ·10^−4	1.05 ·10^7
	2 ·10^11	1	1 ·10^−11	1.40 ·10^−4	1.08 ·10^7

Influence of : considering (-) (respectively n class="Gene">N/m) unaltered, it is possible to observe the influence of —or the influence of the term according to Equation (21)—at varying n>n class="Gene">N/m (respectively -). Figure 13b and Figure 14b show that the critical crack extension increases directly proportional to .

Influence of and n: when the critical opening is kept constant and equals to mm (or m/N), but modifying pan class="Gene">N/m and (-), the critical extension n>n class="Chemical">crack remains unchanged, as shown in Figure 16b. n and only have an influence in the beginning of the R-curve, while the load-displacement plot is almost the same.

Figure 16

Potential cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the cohesive parameters with m/N ( mm) as constant.

With this previous analysis in mind, the general expression for the critical pan class="Chemical">crack extension can be written as: but it is allowed to be reformulated considering Equation (21) as follows:

Due to Figure 15b, let suppose that does not have an influence on , then multiplying by Equation (24) turns into: which becomes identical to Equation (14) if solving for the critical pan class="Chemical">crack opening: where when .

Furthermore, the general expression for the pan class="Chemical">crack extension, in terms of the potential CZM variables, is:

Seeing that it is not possible to isolate and to directly choose , it is then proposed to select at first and to vanish n in Equation (27)—because n depends on for a given . Therefore, the impact of can be separated as observed in Figure 16b, especially in the first stage of the R-curve. Following that idea, the pan class="Chemical">crack extension is then looked at a specific G, for example (-), by means:

Because is supposed to be invariant when is fixed, the last relation becomes identical to Equation (16) and the initial stiffness is found in the same way as for the bilinear CZM: where it is verified that . Finally, the dimensionless value n is established employing Equation (21). Functions and , which are approached within Hermite spline cubic interpolations, are plotted in Figure 17.

Figure 17

Potential cohesive model. The dimensionless functions for model parameters as a function of the crack length: (a) function; (b) function.

5. Comparison of the CZM Characterization

The previously proposed methodology based on the equivalent LEFM R-curve is here applied to characterise the fibre-reinforced plastic tested in Section 2 under mode I pan class="Disease">fracture loading. The trilinear, bilinear and potential CZM are adjusted considering pan class="Gene">N/m and mm. At the same time, in order to verify the effectiveness of this characterization procedure, the parameters of these three laws are also looked for using an optimization algorithm. For this propose, an objective function is built from pan class="Chemical">DCB simulations considering different sets of cohesive parameters . More precisely, 36 sets of parameters are examined for each CZM, where each set enables finding one value of the objective function, considering the vertical difference between the experimental and numerical force-displacement curves, as follows: where is a point in the numerical force-displacement curve considering the set of parameters whereas is the corresponding point in the experimental one. For each , the total number of points taken into consideration was . From the 36 evaluations , the objective function is then approached through a spline interpolation leading to —in order to avoid the expensive evaluation of —and finally it is minimized using a genetic algorithm, as detailed in Figure 18. For the latter, the Scilab optimization toolbox with settings in Table 4 is employed.

Figure 18

Schema of the optimization procedure.

Table 4

Parameters for the genetic algorithm.

poblation size	5000
crossover probability	0.7
mutation probability	0.1
number of generation	300
number of couples	500
pressure	0.05

Table 5 summarizes the identified parameters for the three CZM considering both methodologies. To compare them, the value from Equation (30) is also computed. It is observed that the trilinear CZM achieves the best fitted parameters, for which the equivalent LEFM R-curve ( = 0.33) is lower than the genetic algorithm. In fact, the performance of the genetic algorithm is conditioned by the quality of which in turn depends on the amount of interpolated points , but these are expensive to obtain and thus are avoided. Furthermore, when applying the equivalent LEFM R-curve the bilinear and potential laws are not able to correctly emulate the experimental behaviour (= 16.2 and = 32.2, respectively). However, these both CZM are better adjusted with the genetic algorithm procedure ( = 1.67 and = 2.22, respectively) but in any instance they are not more favourable than the trilinear law.

Table 5

Fitted parameters for each cohesive model using equivalent LEFM and genetic algorithm.

CZM	Characterization	Number of DCB Virtual Tests	Π (ΓCZMfitt)	Δa_c (mm)	Fitted Parameters ΓCZMfitt
trilinear	eq. LEFM R-curve	10	0.33	15.0	f_t=11.51 MPa
					w_c=2.71 mm
					G_fu/G_c=0.83 (-)
	genetic algorithm	36	0.46	14.5	f_t=10.08 MPa
					w_c=2.27 mm
					G_fu/G_c=0.85 (-)
bilinear	eq. LEFM R-curve	10	16.2	15.0	w_c=2.62 mm
					K_0=2.45·10^12 N/m^3
	eq. LEFM R-curve	10	3.21	2.5	w_c=0.064 mm
					K_0=8.67·10^12 N/m^3
	genetic algorithm	36	1.67	4.96	w_c=0.32 mm
					K_0=1.151·10^12 N/m^3
potential	eq. LEFM R-curve	10	32.2	15.0	n=2.1·10^−4 (-)
					K_0=1.12·10^12 N/m^3
	eq. LEFM R-curve	10	2.98	2.5	n=0.86 (-)
					K_0=1·10^12 N/m^3
	genetic algorithm	36	2.22	4.3	n=0.1 (-)
					K_0=4.71·10^11 N/m^3

Traction-separation laws, load-displacement and R-curves are exposed in Figure 19, Figure 20 and Figure 21 for the trilinear, bilinear and potential models, respectively. It is confirmed the high-quality agreement of the trilinear CZM using both fitting methodologies, allowing reaching the critical n class="Chemical">crack extension as well as the maximal applied load P closely to the experimental values. The bilinear and potential laws are unable to follow the entire curves, indeed parameters obtained with the equivalent LEFM R-curve comply with the experimental critical crack extension but the initial response of the R-curve does not agree. On the other hand, parameters fitted with the genetic algorithm are able to follow the beginning of the load-displacement and R-curves; however, the critical crack extensions are lower than the empirical ones—33% (bilinear) and 29% (potential) lower than . Finally, and keeping in mind that the bilinear and potential laws are incapable of tuning behaviours with (i.e., a quasibrittle material), it is reasonable to apply the equivalent LEFM R-curve method taking into account a flexible way to set . For example, taking into account as the when the fracture zone process starts to grow significantly in the R-curve. In Table 5 we include the model parameters considering , subsequently a significant enhancement in the adjustment is achieved. In fact, the objective function is 80% (bilinear) and 91% (potential) lower than considering . Additionally, in Table 5 are listed the number of DCB virtual tests needed for each characterization methodology, a 72% of reduction is reached using the equivalent LEFM R-curve.

Figure 19

Trilinear cohesive model. Numerical DCB test with the fitted parameters: (a) traction-separation law; (b) load-displacement curve; (c) R-curve.

Figure 20

Bilinear cohesive model. Numerical DCB test with the fitted parameters: (a) traction-separation law; (b) load-displacement curve; (c) R-curve.

Figure 21

Potential cohesive model. Numerical DCB test with the fitted parameters: (a) traction-separation law; (b) load-displacement curve; (c) R-curve.

6. 3D Fracture Process Zone (FPZ)

This section is devoted to studying the n class="Chemical">crack front of the DCB test considering the same geometry, material properties and boundary conditions which have been defined in Section 2. The goal is to perform 3D FE simulations and to contrast them against empirical observation, exploiting the fact that damage evolution can be directly observed because specimens are made of GFRP. Actually, in the experimental setup, a camera was placed perpendicularly to the crack plane for recording propagation from the top. FE simulations are carried out using a C++ research code called “MULTI” which is based on a parallel multiscale solver [75] and where the three CZM were implemented employing the best fitted parameters found for each law in Section 5 (see Table 5). The DCB sample is modelled using a 3D mesh with 365,552 linear tetrahedron elements and over two million degrees of freedom, the CZM is treated trough interfaces elements placed on the plane of delamination (78,744 2D triangular elements) between the upper and lower arms of the double cantilever beam. The load-displacement and R-curves are schematised in Figure 22, where 1D beam simulations are also included in order to verify that neglecting the orthotropic material properties (i.e., only GPa was considered in Section 4) was a proper assumption because the geometry and boundary conditions of the problem. Two instants are chosen to compare the plane of delamination: point 1 is located near to the maximum load carrying capacity and point 2 is placed faraway from the instant where propagation begun, as marked in Figure 22.

Figure 22

Experimental and numerical DCB tests (3D and 2D simulations): (a) load-displacement curve; (b) R-curve.

Figure 23 presents a part of the delamination plane in order to show the n class="Chemical">crack tip for both instants of interest (points 1 and 2 in Figure 22). In order to monitor the n>n class="Chemical">crack growth, the experimental sample (see Figure 23a) is marked with 5 mm divisions along the delamination plane beyond of the tip of the pre-crack (that is the yellow area), but the first 5 mm are marked at 1 mm intervals. When the specimen is gradually loaded, it is possible to observe from the top view that the neighbourhood of the crack tip consistently turns a deep white. It is then supposed that this change in appearance is associated with the evolution of the FPZ and not only includes the crack propagation itself. The trilinear law has the thiner whereas the bilinear and potential have a very large , which could be attributed to the initial stiffness (trilinear), (bilinear) and (potential) N/m, respectively for each law. Because the initial stiffness is also employed for the compression behaviour in simulations, low values can induce interpenetration and to enlarge the . The damage distribution through the width is very similar in the four cases. Finally, it can be concluded that even if the load-displacement and resistence curves have a good agreement, the damage distribution over the length is not necessarily the same, especially for tiny values of d.

Figure 23

DCB delamination fronts for the experimental test and 3D simulations, according to points 1 (P1) and 2 (P2) from Figure 22: (a) experimental-P1; (b) trilinear-P1; (c) bilinear-P1; (d) potential-P1; (e) experimental-P2; (f) trilinear-P2; (g) bilinear-P2; (h) potential-P2. The initial pre-crack is indicated by a yellow area. The damage variable d ranges from 0 (blue) to 1 (red) in simulations.

7. Conclusions

In this work, a new identification methodology for bilinear and potential CZM in mode I was developed, inspired by the strategy previously introduced by [26] for a trilinear cohesive law. The main idea is to relate each model parameter with the load-displacement curve and its corresponding equivalent LEFM R-curve through dimensional analyses and numerical simulations. This is implemented mainly in two steps: (1) obtaining the experimental load-displacement test on pan class="Chemical">DCB specimens, computation of the equivalent LEFM R-curve, the critical strain energy release rate and the critical pan class="Chemical">crack extension ; (2) computation of the critical opening and the other corresponding model parameters from dimensionless functions depending on geometry and material of the specimen. Please note that the order of the data reduction for each parameter is crucial. Among the advantages, the parameters’ identification based on the equivalent LEFM R-curve only needs the experimental load-displacement curves, avoiding sophisticated experimental setups and thus it drastically reduces the costs. For validating this strategy, an optimisation algorithm for solving an inverse problem is also here implemented for comparing the identification of bilinear, trilinear and potential laws. Then, the following conclusions are drawn as a result of this research:

it is possible to characterise a bilinear and a potential CZM using a framework based on the equivalent LEFM R-curve;

for the linear, bilinear and potential CZM, the parameters’ identification based on the equivalent LEFM R-curve enables the same accuracy but reduces 72% the numerical efforts respect to a “blind fitting” which minimise the residual between experimental and numerical load-displacement curves;

when applying the equivalent LEFM R-curve framework for characterising a quasibrittle Gn class="Gene">FRP, the trilinear law achieves the best adjustment which is also proven comparing 3D simulations of the n>n class="Disease">fracture process zones. However, it is expected that a trilinear CZM fits materials with large FPZ better than bilinear and potential models. Latter will be fully exploited when characterising more brittle materials;

finally, even if optimisation techniques become popular at present due to their easy numerical implementation, strategies founded on physical models are still better solutions especially when evaluating the objective function is expensive as in mechanical problems.