Next Article in Journal
Conversion of Pectin-Containing By-Products to Pectinases by Bacillus amyloliquefaciens and Its Applications on Hydrolyzing Banana Peels for Prebiotics Production
Next Article in Special Issue
Constant Temperature Approach for the Assessment of Injection Molding Parameter Influence on the Fatigue Behavior of Short Glass Fiber Reinforced Polyamide 6
Previous Article in Journal
Patchy Micelles with a Crystalline Core: Self-Assembly Concepts, Properties, and Applications
Previous Article in Special Issue
Fabrication and Thermo-Electro and Mechanical Properties Evaluation of Helical Multiwall Carbon Nanotube-Carbon Fiber/Epoxy Composite Laminates
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Magíster en Ciencias de la Ingeniería c/m Ingeniería Mecánica, Facultad de Ingeniería, Campus Curicó, Universidad de Talca, Curicó 3340000, Chile
2
Departamento de Tecnologías Industriales, Campus Curicó, Universidad de Talca, Curicó 3340000, Chile
3
Departamento de Ingeniería en Obras Civiles, Facultad de Ingeniería, Universidad de Santiago de Chile (USACH), Santiago 9170124, Chile
*
Author to whom correspondence should be addressed.
Polymers 2021, 13(9), 1482; https://doi.org/10.3390/polym13091482
Submission received: 10 February 2021 / Revised: 19 April 2021 / Accepted: 24 April 2021 / Published: 4 May 2021
(This article belongs to the Special Issue Mechanics of Polymer and Polymer Composite Materials and Structures)

Abstract

:
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.

Graphical Abstract

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 prediction 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 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 G i , according respectively to the three fracture modes ( G I , G I I and G I I I ) 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 G i c (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 a 0 , 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]:
G = P 2 2 B d C ( a ) d a
where C ( a ) = δ / P is the compliance or displacement δ to applied load P ratio.
The double cantilever beam (DCB) test shown in Figure 1a and used in this work is the most widely used procedure for the measurement of mode I delamination fracture toughness. For the load-displacement curve (see Figure 1b) and considering a rigid foundation at the crack ending, the compliance can be given by the beam theory as C ( a 0 ) = 2 a 0 3 / 3 E I , where I = B h 3 / 12 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].
The R-curve shown in Figure 1c illustrates the variation of the material crack resistance with respect to the crack propagation length Δ a = a a 0 . 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 G I remains constant and equal to G I c (for the sake of simplicity subindex ( . ) I 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 Δ a c , 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 l f p z 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 C ( a ) and the crack growth resistance is then determined from Equation (1) [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 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 f ( w ) 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 K 0 and a maximum critical traction level f t . Upon reaching f t , the softening function is described by a negative tangent stiffness until a critical displacement jump w c is achieved. At this point, the system presents no more load-bearing capacity, i.e., f ( w c ) = 0 . 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 G c , 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 f ( w ) 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 P δ 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 f ( w ) . 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).

2. Experimental Test

DCB specimens were manufactured through a vacuum infusion process, considering a fibre/resin weight ratio of 0.47 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 company Fibratec (Santiago, Chile) was employed, while the reinforcement is a unidirectional glass fibre from the German company P-D Interglas Technologies GmbH (now acquired by Porcher Industries, Eclose Badinieres, France), currently named as UD 220 g/m 2 (i.e., with 207 g/m 2 in the warp direction and 13 g/m 2 in the fill direction). Then, the composite, previously characterised in [70], has the following elastic properties: E 1 = 32.1 GPa, E 2 = 12.6 GPa, ν 12 = 0.1 (-), ν 23 = 0.24 (-) and G 12 e x p = 3 GPa, where 1 direction is in the fibre direction and aligned through the specimen’s length.
The 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.

3. Cohesive Zone Models

Traction-separation laws f ( w ) 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 f ( w ) is a vector-valued function, while the initial stiffness K 0 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 K 0 is progressively weaken according to:
f ( w ) = ( 1 d ( w + ) ) K 0 w + + K 0 w
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., d | t = max τ t ( d | τ ) .
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 fracture energy G c corresponds to the area under the curve. However, the energy is decomposed into two parts ( G c = G f μ + G f b ) for the trilinear CZM, which has been attributed to micro-cracking ( G f μ ) and fiber-bridging ( G f b ) [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.

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 crack a (see Figure 2) located at some distance ahead of the initial crack a 0 (or the current stress-free crack a s f ). 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 w c 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 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 E 1 = 32.1 GPa is taken into account), whereas the critical energy release rate is given by the experimental test ( G c e x p = 982.2 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 w c , the distribution of the critical energy release rate G f μ / G c and the tensile strength f t . The critical crack extension Δ a c will be then associated with these parameters and the specimen size through a dimensional analysis. In this analysis w 0 = 10 7 mm and G c = 982.2 N/m are kept constant.
  • Influence of the critical opening: numerical simulations are carried out affecting the critical opening w c { 1 , 2 , 4 , 5 , 7 } mm but fixing the values G f μ / G c = 0.5 (-) and f t = 1.15 · 10 7 Pa. From Figure 6, it is observed that w c does not have an influence before reaching the 50% of G c . After that, the critical crack extension Δ a c and the critical displacement jump w c have a positive correlation, which means that the length of the fracture process zone l f p z increases when w c increases.
  • Influence of the fracture energy distribution: keeping constant f t = 1.15 · 10 7 Pa and w c = 2.72 mm, Figure 7 shows that at varying the critical energy release rate G f μ / G c { 0.5 , 0.65 , 0.75 } (-), the load-displacement plot and the equivalent LEFM R-curve are invariable if the dissipated energies G f μ / G c are lower than { 0.5 , 0.65 , 0.75 } , respectively. If the ratio G f μ / G c 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 crack extension Δ a c always remains the same.
  • Influence of the tensile strength: Figure 8 shows the impact at varying f t { 3 , 5 , 7 , 10 , 15 } [ 10 7 Pa] whereas G f μ / G c = 0.5 (-) and w c = 2.72 mm are unchanged. It can be concluded that the tensile strength impacts the response at the beginning of both curves: if f t increases, the behaviour becomes a brittle one. When the dissipated energy reaches the value G f μ , the fracture response starts to be the same for any value of f t and the critical crack extension Δ a c becomes identical.
Figure 6. Trilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the critical opening w c with f t = 1.15 · 10 7 Pa and G f μ / G c = 0.5 (-) as constants.
Figure 6. Trilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the critical opening w c with f t = 1.15 · 10 7 Pa and G f μ / G c = 0.5 (-) as constants.
Polymers 13 01482 g006
Figure 7. Trilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the fracture energy distribution G f μ / G c with w c = 2.72 mm and f t = 1.15 · 10 7 Pa as constants.
Figure 7. Trilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the fracture energy distribution G f μ / G c with w c = 2.72 mm and f t = 1.15 · 10 7 Pa as constants.
Polymers 13 01482 g007
Figure 8. Trilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the tensile strength f t with w c = 2.72 mm and G f μ / G c = 0.5 (-) as constants.
Figure 8. Trilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the tensile strength f t with w c = 2.72 mm and G f μ / G c = 0.5 (-) as constants.
Polymers 13 01482 g008
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 Δ a c depends on the specimen size and the cohesive parameters, as proposed in [26,74]:
Δ a c l c h = ϕ 1 D l c h , w c w c h , G f μ G c
where D is a characteristic dimension of the specimen, l c h and w c h are the Hillerborg’s characteristic length and a characteristic crack opening, respectively:
l c h = E G c f t 2 , w c h = G c f t
Because it was observed that the ratio G f μ / G c does not have influence on the critical crack extension Δ a c (see Figure 7b), the third argument of Equation (3) is able to be directly vanished. According to Figure 8b, Δ a c neither depends on f t , therefore if Equation (3) is homogeneous in f t 2 —i.e., depends on w c 2 —it can be cancelled when factoring by w c 2 w c h 2 [26]:
Δ a c l c h = w c 2 w c h 2 ϕ 2 D w c h 2 w c 2 l c h , 1 or Δ a c D = w c 2 E G c D ϕ w c 2 E G c D
where ϕ ( · ) = ϕ 2 ( 1 / · , 1 ) , furthermore the right equation has been multiplied by 1 / D . On the other hand, a lower bound of the critical crack extension has been previously studied in [74] for D , given by:
lim D Δ a c π 32 w c 2 E G c
Then, it is expected that ϕ ( 0 ) π / 32 .
The nonlinear expression relating Δ a c and w c , Equation (5), has the ability to be solved for w c through the equivalent LEFM R-curve and FE computations, according to:
w c = 32 Δ a c G c π E ψ Δ a c D
where ψ ( 0 ) 1 when D .
The critical opening w c can be now found employing Equation (7), while the other cohesive parameter—the tensile strength f t —can also be studied through a dimensional analysis. Therefore, the crack length Δ a , for a given energy release rate G < G c , is allowed to be obtained from the following general expression:
Δ a l c h = ζ 1 G G c , D l c h , w c w c h , G f μ G c
From virtual tests previously carried out for different values of w c and G f μ / G c , but keeping f t constant, it is observed that the equivalent R-curve remains almost the same when G < 0.5 G c (see Figure 6b and Figure 7b). In fact, f t only has an influence at the beginning of the R-curve if w c and G f μ / G c are both fixed (see Figure 8b). Therefore, function ζ 1 ( · ) is able to be considered independent of w c / w c h and G f μ / G c if G < 0.5 G c . For instance, calculations are here proposed with G / G c = 0.2 -:
Δ a 0.2 l c h = ζ 2 0.2 , D l c h = ζ 3 D l c h
As proposed in [26], an expression relating Δ a 0.2 / D and D / l c h is permitted to be obtained at multiplying Equation (9) by l c h / D , then:
f t = E G c D χ D Δ a 0.2
where χ ( · ) can be determined from Figure 8b, considering the horizontal line for G = 0.2 G c and it verifies χ ( 0 ) 0 . Finally, the dimensionless functions ψ ( · ) and χ ( · ) plotted in Figure 9 are obtained through a Hermite spline cubic interpolation.

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 w c and the initial stiffness K 0 , while G c is always kept constant and equals to G c e x p .
  • Influence of the critical opening: w c has an inverse correlation with the maximal applied load P when considering w c { 0.1 , 0.2 , 0.4 , 1 , 4 } mm and keeping unchanged K 0 = 10 13 N/m 3 , as observed in Figure 10a. Moreover, Figure 10b shows that the critical crack extension Δ a c is inversely proportional to the tensile strength f t , i.e., the interface becomes more brittle for higher values of f t .
  • Influence of the initial stiffness K 0 : from Figure 11b it is observed that K 0 { 5 , 10 , 22 , 10 2 , 10 4 } · 10 11 N/m 3 does not have any influence on the critical crack extension Δ a c if f t is kept constant. K 0 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 10. Bilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the critical opening w c with K 0 = 10 13 N/m 3 as constant.
Figure 10. Bilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the critical opening w c with K 0 = 10 13 N/m 3 as constant.
Polymers 13 01482 g010
Figure 11. Bilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the initial stiffness K 0 with w c = 0.13 mm as constant.
Figure 11. Bilinear cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the initial stiffness K 0 with w c = 0.13 mm as constant.
Polymers 13 01482 g011
As previously introduced in Section 4.1, the critical crack extension Δ a c can be expressed as a function which depends on the specimen size and on the cohesive parameters as:
Δ a c l c h = γ 1 D l c h , w c w c h , K 0 K c h
where K c h is a characteristic stiffness defined as follows:
K c h = E 4 G c f t 2
From Figure 11b it is assumed that K 0 does not have influence on the critical crack extension Δ a c , then the third argument disappears. Besides, the relation f t = 2 G c / w c enables going without having the second variable w c / w c h = 0.5 , then multiplying by 1 / D and considering Equation (4):
Δ a c D l c h = 1 D γ 2 D l c h or Δ a c D = w c 2 E 4 G c D γ 2 4 G c D w c 2 E or Δ a c D = w c 2 E G c D γ w c 2 E G c D
where it is expected, when D , i.e., γ ( 0 ) π / 32 in concordance with Equation (6). Finally, the critical opening is computed in the same way as for the trilinear law:
w c = 32 Δ a c G c π E η Δ a c D
where η ( 0 ) 1 when D .
To characterise the stiffness parameter K 0 , a dimensionless general expression for the crack extension is able to be written as:
Δ a l c h = ν 1 G G c , D l c h , w c w c h , K 0 K c h
From Figure 10b it should be noticed that the effect of K 0 on the R-curve depends on the selection of w c . First, it is assumed that w c has been chosen using Equation (14); secondly, the crack extension is searched for a given G where the impact of K 0 is significative (e.g., G / G c = 0.01 -), subsequently last equation becomes:
Δ a 0.01 l c h = ν 2 D l c h , K 0 K c h
The dependence on f t (i.e., on w c ) has the chance to be vanished if factoring by E 4 K 0 G c f t 2 , as follows:
Δ a 0.01 l c h = ν 2 D f t 2 E G c , K 0 G c f t 2 E 4 = K 0 G c f t 2 E 4 ν 3 D E 3 G c 2 K 0 , 1
Then, multiplying last this expression by 1 / D :
Δ a 0.01 D = K 0 G c 2 D E 3 ν 3 D E 3 G c 2 K 0 , 1 = K 0 G c 2 D E 3 ν 4 D E 3 G c 2 K 0 = ν K 0 G c 2 D E 3
Finally, the initial stiffness can be stablished employing the next expression:
K 0 = D E 3 G c 2 ξ D Δ a 0.01
where it is verified that K 0 ( ) . The dimensionless functions, η ( · ) and ξ ( · ) , characterizing both cohesive parameters are obtained through a Hermite spline cubic interpolation and plotted in Figure 12.

4.3. Potential CZM

The parameters enabling identify the potential CZM are the initial stiffness K 0 and the dimensionless variable n. However, it is not possible to separately relate them to the critical crack extension Δ a c . Actually, for a fixed n, K 0 influences the whole R-curve (including Δ a c ), as shown in Figure 13. Same behaviour takes place for different values of n but a fixed K 0 , as demonstrated in Figure 14. In both cases G c is always kept constant and equals to G c e x p .
Despite this impossibility, it is feasible to find a connection between (n, K 0 ) and Δ a c if paying attention to the following relation from Table 3:
d ( w ) = ( n + 1 ) K 0 w 2 2 n G c n
More precisely, when considering d = 1 , the critical displacement jump is able to be written in terms of n and K 0 as follows:
w c 2 = 2 G c n + 1 n K 0
Furthermore, it is possible to obtain an expression for f t through a stationary point of Equation (2), then:
f t = 2 1 / 2 n n K 0 ( n + 1 2 ) n + 1 2 n 2 G c n + 1 n K 0 1 / 2
When examining Equation (22) for several values of n and K 0 , as plotted in Figure 15 and detailed in Table A1 in Appendix A, it is noticed that f t is essentially unchanged if the term ( n + 1 ) / ( n K 0 ) (or w c ) remains invariable (series A and D). Nevertheless, if one parameter (n or K 0 ) is kept constant and the other is varied, f t increases when K 0 or n increases, respectively (series B and C). In the following, the effect of ( n + 1 ) / ( n K 0 ) on the equivalent R-curve is studied, while G c is equals to G c e x p .
  • Influence of w c : considering n = 0.7 (-) (respectively K 0 = 10 13 N/m 3 ) unaltered, it is possible to observe the influence of w c —or the influence of the term ( n + 1 ) / ( n K 0 ) according to Equation (21)—at varying K 0 { 1.12 , 5.6 , 11.2 , 22.4 , 33.6 } · 10 11 N/m 3 (respectively n { 10 5 , 10 4 , 10 3 , 10 2 , 10 0 } -). Figure 13b and Figure 14b show that the critical crack extension Δ a c increases directly proportional to w c .
  • Influence of K 0 and n: when the critical opening is kept constant and equals to w c = 2.9 mm (or ( n + 1 ) / ( n K 0 ) = 4.25 · 10 9 m 3 /N), but modifying K 0 { 2 , 8 , 10 , 1120 } · 10 9 N/m 3 and n { 1.13 , 0.3 , 3.8 , 0.24 , 0.0021 } · 10 1 (-), the critical extension crack remains unchanged, as shown in Figure 16b. n and K 0 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 n + 1 n K 0 = 4.25 · 10 9 m 3 /N ( w c = 2.9 mm) as constant.
Figure 16. Potential cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the cohesive parameters with n + 1 n K 0 = 4.25 · 10 9 m 3 /N ( w c = 2.9 mm) as constant.
Polymers 13 01482 g016
With this previous analysis in mind, the general expression for the critical crack extension can be written as:
Δ a c l c h = μ 1 D l c h , n , K 0 K c h
but it is allowed to be reformulated considering Equation (21) as follows:
Δ a c l c h = μ 2 D l c h , w c 2 w c h 2
Due to Figure 15b, let suppose that f t does not have an influence on Δ a c , then multiplying by 1 / D Equation (24) turns into:
Δ a c l c h = w c 2 w c h 2 μ 3 D w c h 2 w c 2 l c h , 1 or Δ a c D = w c 2 E G c D μ E w c 2 G c D
which becomes identical to Equation (14) if solving for the critical crack opening:
w c = 32 Δ a c G c π E ς Δ a c D
where ς ( 0 ) 1 when D .
Furthermore, the general expression for the crack extension, in terms of the potential CZM variables, is:
Δ a l c h = β G G c , D l c h , n , K 0 K c h
Seeing that it is not possible to isolate and to directly choose K 0 , it is then proposed to select w c at first and to vanish n in Equation (27)—because n depends on K 0 for a given w c . Therefore, the impact of K 0 can be separated as observed in Figure 16b, especially in the first stage of the R-curve. Following that idea, the crack extension is then looked at a specific G, for example G / G c = 0.01 (-), by means:
Δ a 0.01 l c h = β 1 D l c h , K 0 K c h
Because f t is supposed to be invariant when w c 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:
K 0 = D E 3 G c 2 ϰ D Δ a 0.01
where it is verified that K 0 ( ) . 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.

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 fracture loading. The trilinear, bilinear and potential CZM are adjusted considering G c e x p = 982.2 N/m and Δ a c e x p = 15 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 DCB simulations considering different sets of cohesive parameters Γ j . More precisely, 36 sets of parameters Γ j 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:
Π ( Γ j ) = Π j = i = 1 N 1 N | P n , i j P e , i | 2
where P n , i j is a point in the numerical force-displacement curve considering the set of parameters Γ j whereas P e , i is the corresponding point in the experimental one. For each Π j , the total number of points taken into consideration was N = 20 . From the 36 evaluations Π j , the objective function Π is then approached through a spline interpolation leading to Π appr —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.
Table 5 summarizes the identified parameters Γ CZM fitt for the three CZM considering both methodologies. To compare them, the value Π ( Γ CZM fitt ) 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 28.3 % lower than the genetic algorithm. In fact, the performance of the genetic algorithm is conditioned by the quality of Π appr which in turn depends on the amount of interpolated points Π j , 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.
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 crack extension Δ a c 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 Δ a c e x p . Finally, and keeping in mind that the bilinear and potential laws are incapable of tuning behaviours with Δ a c > > 0 (i.e., a quasibrittle material), it is reasonable to apply the equivalent LEFM R-curve method taking into account a flexible way to set Δ a c . For example, taking into account Δ a c as the Δ a when the fracture zone process starts to grow significantly in the R-curve. In Table 5 we include the model parameters considering Δ a 0.56 , subsequently a significant enhancement in the adjustment is achieved. In fact, the objective function is 80% (bilinear) and 91% (potential) lower than considering Δ a c e x p . 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.

6. 3D Fracture Process Zone (FPZ)

This section is devoted to studying the 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 E 1 = 32.1 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 23 presents a part of the delamination plane in order to show the crack tip for both instants of interest (points 1 and 2 in Figure 22). In order to monitor the 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 l f p z whereas the bilinear and potential have a very large l f p z , which could be attributed to the initial stiffness K 0 10 18 (trilinear), 10 12 (bilinear) and 10 11 (potential) N/m 3 , respectively for each law. Because the initial stiffness K 0 is also employed for the compression behaviour in simulations, low values can induce interpenetration and to enlarge the l f p z . 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.

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 DCB specimens, computation of the equivalent LEFM R-curve, the critical strain energy release rate G c and the critical crack extension Δ a c ; (2) computation of the critical opening w c 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 GFRP, the trilinear law achieves the best adjustment which is also proven comparing 3D simulations of the 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.

Author Contributions

Conceptualization, K.S. and J.C.P.; Data curation, L.T.; Formal analysis, L.T.; Investigation, L.T.; Methodology, K.S.; Resources, K.S. and G.P.; Supervision, K.S. and G.P.; Visualization, L.T.; Writing—original draft, L.T. and K.S.; Writing—review & editing, K.S. and J.C.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by CONICYT through the FONDECYT Initiation projects N 11160633 and N 11130623.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is contained within the article.

Acknowledgments

The authors wish to thank LMT-Cachan (ENS Paris-Saclay/CNRS/Université Paris-Saclay) for enabling using the research code “MULTI” and to thank Juan Pablo Muñoz for providing us part of the 1D FE code.

Conflicts of Interest

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.

Appendix A

Table A1 summaries the study carried out to understand the influence of the cohesive parameters on the critical traction f t and on the critical displacement jump w c in the case of the potential law.
Table A1. Potential cohesive model. Influence of model parameters on f t and w c .
Table A1. Potential cohesive model. Influence of model parameters on f t and w c .
K 0 (N/m 3 )n (-) n + 1 nK 0 (m 3 /N) w c (m) f t (Pa)
1 · 10 14 0.011 · 10 12 4.43 · 10 5 3.26 · 10 7
5 · 10 13 0.021 · 10 12 4.43 · 10 5 3.26 · 10 7
1 · 10 13 0.111 · 10 12 4.43 · 10 5 3.27 · 10 7
serie A5 · 10 12 0.251 · 10 12 4.43 · 10 5 3.28 · 10 7
3 · 10 12 0.51 · 10 12 4.43 · 10 5 3.32 · 10 7
2 · 10 12 11 · 10 12 4.43 · 10 5 3.41 · 10 7
1.5 · 10 12 21 · 10 12 4.43 · 10 5 3.56 · 10 7
1.3 · 10 12 3.31 · 10 12 4.43 · 10 5 3.69 · 10 7
1 · 10 13 21.5 · 10 13 1.72 · 10 5 9.18 · 10 7
1 · 10 13 12 · 10 13 1.98 · 10 5 7.63 · 10 7
1 · 10 13 0.53 · 10 13 2.43 · 10 5 6.07 · 10 7
serie B1 · 10 13 0.11.1 · 10 12 4.65 · 10 5 3.11 · 10 7
1 · 10 13 0.011 · 10 11 1.41 · 10 5 1.03 · 10 7
1 · 10 13 0.0011 · 10 10 4.43 · 10 4 3.26 · 10 6
1 · 10 13 0.00011 · 10 9 1.4 · 10 3 1.03 · 10 6
1 · 10 13 0.000011 · 10 8 4.43 · 10 3 3.26 · 10 5
1 · 10 14 0.53 · 10 14 7.68 · 10 6 1.92 · 10 8
5 · 10 13 0.56 · 10 14 1.09 · 10 5 1.36 · 10 8
1 · 10 13 0.53 · 10 13 2.43 · 10 5 6.07 · 10 7
serie C5 · 10 12 0.56 · 10 13 3.43 · 10 5 4.29 · 10 7
3 · 10 12 0.51 · 10 12 4.43 · 10 5 3.32 · 10 7
1 · 10 12 0.53 · 10 12 7.68 · 10 5 1.92 · 10 7
5 · 10 11 0.56 · 10 12 1.09 · 10 4 1.36 · 10 7
1 · 10 11 0.53 · 10 11 2.43 · 10 4 6.07 · 10 6
1 · 10 14 0.0011 · 10 11 1.40 · 10 4 1.03 · 10 7
5 · 10 13 0.0021 · 10 11 1.40 · 10 4 1.03 · 10 7
1 · 10 13 0.01011 · 10 11 1.40 · 10 4 1.03 · 10 7
serie D5 · 10 12 0.02041 · 10 11 1.40 · 10 4 1.03 · 10 7
1 · 10 12 0.1111 · 10 11 1.40 · 10 4 1.03 · 10 7
5 · 10 11 0.251 · 10 11 1.40 · 10 4 1.04 · 10 7
3 · 10 11 0.51 · 10 11 1.40 · 10 4 1.05 · 10 7
2 · 10 11 11 · 10 11 1.40 · 10 4 1.08 · 10 7

References

  1. González, C.; Vilatela, J.; Molina-Aldareguía, J.; Lopes, C.; LLorca, J. Structural composites for multifunctional applications: Current challenges and future trends. Prog. Mater. Sci. 2017, 89, 194–251. [Google Scholar] [CrossRef] [Green Version]
  2. Tabiei, A.; Zhang, W. Composite Laminate Delamination Simulation and Experiment: A Review of Recent. Appl. Mech. Rev. 2018, 70, 030801. [Google Scholar] [CrossRef]
  3. Fish, J.E. Multiscale Methods: Bridging the Scales in Science and Engineering; OUP Oxford: Oxford, UK, 2009. [Google Scholar]
  4. Geers, M.; Kouznetsova, V.; Brekelmans, W. Multi-scale computational homogenization: Trends and challenges. J. Comput. Appl. Math. 2010, 234, 2175–2182. [Google Scholar] [CrossRef]
  5. LLorca, J.; González, C.; Molina-Aldareguía, J.M.; Lópes, C.S. Multiscale Modeling of Composites: Toward Virtual Testing … and Beyond. JOM 2013, 65, 215–225. [Google Scholar] [CrossRef] [Green Version]
  6. Okereke, M.; Akpoyomare, A.; Bingley, M. Virtual testing of advanced composites, cellular materials and biomaterials: A review. Compos. Part B Eng. 2014, 60, 637–662. [Google Scholar] [CrossRef]
  7. Bouvet, C.; Castanié, B.; Bizeul, M.; Barrau, J. Low velocity impact modelling in laminate composite panels with discrete interface elements. Int. J. Solids Struct. 2009, 46, 2809–2821. [Google Scholar] [CrossRef] [Green Version]
  8. Pakdel, H.; Mohammadi, B. Stiffness degradation of composite laminates due to matrix cracking and induced delamination during tension-tension fatigue. Eng. Fract. Mech. 2019, 216, 106489. [Google Scholar] [CrossRef]
  9. Green, B.; Wisnom, M.; Hallett, S. An experimental investigation into the tensile strength scaling of notched composites. Compos. Part A Appl. Sci. Manuf. 2007, 38, 867–878. [Google Scholar] [CrossRef]
  10. Sridharan, S. (Ed.) Delamination Behaviour of Composites; Series in Composites Science and Engineering; Woodhead Publishing: Cambridge, UK, 2008. [Google Scholar]
  11. Wisnom, M.R. The role of delamination in failure of fibre-reinforced composites. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 2012, 370, 1850–1870. [Google Scholar] [CrossRef] [PubMed]
  12. Jäger, S.; Pickett, A.; Middendorf, P. A Discrete Model for Simulation of Composites Plate Impact Including Coupled Intra- and Inter-ply Failure. Appl. Compos. Mater. 2016, 23, 179–195. [Google Scholar] [CrossRef]
  13. Saavedra, K.; Allix, O.; Gosselet, P.; Hinojosa, J.; Viard, A. An enhanced nonlinear multi-scale strategy for the simulation of buckling and delamination on 3D composite plates. Comput. Methods Appl. Mech. Eng. 2017, 317, 952–969. [Google Scholar] [CrossRef]
  14. ASTM-D5528-94a. Standard Test Method for Mode I Inter-laminar Fracture Toughness of Unidirectional Continuous Fiber Rein- forced Polymer Matrix Composites; ASTM: Philadelphia, PA, USA, 1994. [Google Scholar]
  15. ASTM-D6671. Standard Test Method for Mixed Mode I–Mode II Interlaminar Fracture Toughness of Unidirectional Fiber Reinforced Polymer Matrix Composites; ASTM International: West Conshohocken, PA, USA, 2013. [Google Scholar]
  16. ASTM-D9705-14. Standard Test Method for Determination of the Mode II Interlaminar Fracture Toughness of Unidirectional Fiber-Reinforced Polymer Matrix Composites; Protocols for Interlaminar Fracture Testing of Composites; European Structural Integrity Society (ESIS): Delft, The Netherlands, 1993. [Google Scholar]
  17. ISO-15024. Fiber-reinforced Plastic Composites—Determination of Mode I Interlaminar Fracture Toughness, GIc, for Unidirectionally Reinforced Materials; International Organization for Standardization: Geneva Switzerland, 2001. [Google Scholar]
  18. Irwin, G.R.; Kies, J.A. Critical Energy Rate Analysis of Fracture Strength. Weld. J. Res. Suppl. 1954, 33, 193–198. [Google Scholar]
  19. Kanninen, M.F. An augmented double cantilever beam model for studying crack propagation and arrest. Int. J. Fract. 1973, 9, 83–92. [Google Scholar]
  20. Allix, O.; Ladevéze, P.; Corigliano, A. Damage analysis of interlaminar fracture specimens. Compos. Struct. 1995, 31, 61–74. [Google Scholar] [CrossRef]
  21. Suo, Z.; Bao, G.; Fan, B. Delamination R-curve phenomena due to damage. J. Mech. Phys. Solids 1992, 40, 1–16. [Google Scholar] [CrossRef]
  22. Frossard, G.; Cugnoni, J.; Gmür, T.; Botsis, J. Mode I interlaminar fracture of carbon epoxy laminates: Effects of ply thickness. Compos. Part A Appl. Sci. Manuf. 2016, 91, 1–8. [Google Scholar] [CrossRef]
  23. Canal, L.P.; Alfano, M.; Botsis, J. A multi-scale based cohesive zone model for the analysis of thickness scaling effect in fiber bridging. Compos. Sci. Technol. 2017, 139, 90–98. [Google Scholar] [CrossRef]
  24. Cameselle-Molares, A.; Vassilopoulos, A.P.; Renart, J.; Turon, A.; Keller, T. Numerical simulation of two-dimensional in-plane crack propagation in FRP laminates. Compos. Struct. 2018, 200, 396–407. [Google Scholar] [CrossRef]
  25. Bazǎnt, Z. Concrete fracture models: Testing and practice. Eng. Fract. Mech. 2002, 69, 165–205. [Google Scholar] [CrossRef]
  26. Morel, S.; Lespine, C.; Coureau, J.L.; Planas, J.; Dourado, N. Bilinear softening parameters and equivalent LEFM R-curve in quasibrittle failure. Int. J. Solids Struct. 2010, 47, 837–850. [Google Scholar] [CrossRef] [Green Version]
  27. Rybicki, E.F.; Kanninen, M.F. A finite element calculation of stress intensity factors by a modified crack closure integral. Eng. Fract. Mech. 1977, 9, 931–938. [Google Scholar] [CrossRef]
  28. Krueger, R. The Virtual Crack Closure Technique for modeling interlaminar failure and delamination in advanced composite materials. In Numerical Modelling of Failure in Advanced Composite Materials; Camanho, P.P., Hallett, S.R., Eds.; Woodhead Publishing Series in Composites Science and Engineering; Woodhead Publishing: Cambridge, UK, 2015; pp. 3–53. [Google Scholar]
  29. Marjanović, M.; Meschke, G.; Vuksanović, D. A finite element model for propagating delamination in laminated composite plates based on the Virtual Crack Closure method. Compos. Struct. 2016, 150, 8–19. [Google Scholar] [CrossRef]
  30. Belytschko, T.; Black, T. Elastic crack growth in finite elements with minimal remeshing. Int. J. Numer. Methods Eng. 1999, 45, 601–620. [Google Scholar] [CrossRef]
  31. Moës, N.; Dolbow, J.; Belytschko, T. A finite element method for crack growth without remeshing. Int. J. Numer. Methods Eng. 1999, 46, 131–150. [Google Scholar] [CrossRef]
  32. Sosa, J.C.; Karapurath, N. Delamination modelling of GLARE using the extended finite element method. Compos. Sci. Technol. 2012, 72, 788–791. [Google Scholar] [CrossRef]
  33. Yazdani, S.; Rust, W.J.; Wriggers, P. An XFEM approach for modelling delamination in composite laminates. Compos. Struct. 2016, 135, 353–364. [Google Scholar] [CrossRef]
  34. Dugdale, D.S. Yielding of steel sheets containing slits. J. Mech. Phys. Solids 1960, 8, 100–104. [Google Scholar] [CrossRef]
  35. Barenblatt, G. The Mathematical Theory of Equilibrium Cracks in Brittle Fracture. In Advances in Applied Mechanics; Elsevier: Amsterdam, The Netherlands, 1962; Volume 7, pp. 55–129. [Google Scholar]
  36. Hillerborg, A.; Modéer, M.; Petersson, P.E. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem. Concr. Res. 1976, 6, 773–781. [Google Scholar] [CrossRef]
  37. Heidari-Rarani, M.; Sayedain, M. Finite element modeling strategies for 2D and 3D delamination propagation in composite DCB specimens using VCCT, CZM and XFEM approaches. Theor. Appl. Fract. Mech. 2019, 103. [Google Scholar] [CrossRef]
  38. Petersson, P. Crack Growth and Development of Fracture Zones in Plain Concrete and Similar Materials. Ph.D. Thesis, Lund University, Lund, Sweden, 1981. [Google Scholar]
  39. Mai, Y.; Lawn, B. Crack-Interface Grain Bridging as a Fracture Resistance Mechanism in Ceramics: II, Theoretical Fracture Mechanics Model. J. Am. Ceram. Soc. 1984, 70, 289–294. [Google Scholar] [CrossRef]
  40. Allix, O.; Lévêque, D.; Perret, L. Identification and forecast of delamination in composite laminates by an interlaminar interface model. Compos. Sci. Technol. 1998, 58, 671–678. [Google Scholar] [CrossRef]
  41. Tvergaard, V.; Hutchinson, J. Effect of T-Stress on mode I crack growth resistance in a ductile solid. Int. J. Solids Struct. 1994, 31, 823–833. [Google Scholar] [CrossRef]
  42. Rakin, M.; Medjo, B.; Gubeljak, N.; Sedmak, A. Micromechanical assessment of mismatch effects on fracture of high-strength low alloyed steel welded joints. Eng. Fract. Mech. 2013, 109, 221–235. [Google Scholar] [CrossRef]
  43. Tijssens, M.; van der Giessen, E.; Sluys, L. Modeling of crazing using a cohesive surface methodology. Mech. Mater. 2000, 32, 19–35. [Google Scholar] [CrossRef] [Green Version]
  44. Dourado, N.; Morel, S.; de Moura, M.; Valentin, G.; Morais, J. Comparison of fracture properties of two wood species through cohesive crack simulations. Compos. Part A Appl. Sci. Manuf. 2008, 39, 415–427. [Google Scholar] [CrossRef]
  45. de Morais, A.; Pereira, A. Application of the effective crack method to mode I and mode II interlaminar fracture of carbon/epoxy unidirectional laminates. Compos. Part A Appl. Sci. Manuf. 2007, 38, 785–794. [Google Scholar] [CrossRef]
  46. Su, Z.; Tay, T.; Ridha, M.; Chen, B. Progressive damage modeling of open-hole composite laminates under compression. Compos. Struct. 2015, 122, 507–517. [Google Scholar] [CrossRef]
  47. Flores, E.S.; Saavedra, K.; Hinojosa, J.; Chandra, Y.; Das, R. Multi-scale modelling of rolling shear failure in cross-laminated timber structures by homogenisation and cohesive zone models. Int. J. Solids Struct. 2016, 81, 219–232. [Google Scholar] [CrossRef]
  48. Yang, Y.; Liu, X.; Wang, Y.Q.; Gao, H.; Li, R.; Bao, Y. A progressive damage model for predicting damage evolution of laminated composites subjected to three-point bending. Compos. Sci. Technol. 2017, 151, 85–93. [Google Scholar] [CrossRef]
  49. Koloor, S.R.; Ayatollahi, M.; Tamin, M. Elastic-damage deformation response of fiber-reinforced polymer composite laminates with lamina interfaces. J. Reinf. Plast. Compos. 2017, 36, 832–849. [Google Scholar] [CrossRef]
  50. Koloor, S.R.; Tamin, M. Mode-II interlaminar fracture and crack-jump phenomenon in CFRP composite laminate materials. Compos. Struct. 2018, 204, 594–606. [Google Scholar] [CrossRef]
  51. Confalonieri, F.; Perego, U. A new framework for the formulation and validation of cohesive mixed-mode delamination models. Int. J. Solids Struct. 2019, 164, 168–190. [Google Scholar] [CrossRef]
  52. Park, K.; Paulino, G.H. Cohesive Zone Models: A Critical Review of Traction-Separation Relationships Across Fracture Surfaces. Appl. Mech. Rev. 2013, 64, 060802. [Google Scholar] [CrossRef]
  53. Sorensen, L.; Botsis, J.; Gmür, T.; Humbert, L. Bridging tractions in mode I delamination: Measurements and simulations. Compos. Sci. Technol. 2008, 68, 2350–2358. [Google Scholar] [CrossRef]
  54. Stutz, S.; Cugnoni, J.; Botsis, J. Crack–fiber sensor interaction and characterization of the bridging tractions in mode I delamination. Eng. Fract. Mech. 2011, 78, 890–900. [Google Scholar] [CrossRef]
  55. Pappas, G.; Canal, L.; Botsis, J. Characterization of intralaminar mode I fracture of AS4/PPS composite using inverse identification and micromechanics. Compos. Part A Appl. Sci. Manuf. 2016, 91, 117–126. [Google Scholar] [CrossRef]
  56. Alfano, M.; Lubineau, G.; Paulino, G.H. Global sensitivity analysis in the identification of cohesive models using full-field kinematic data. Int. J. Solids Struct. 2015, 55, 66–78. [Google Scholar] [CrossRef]
  57. Ostapska, K.; Malo, K.A. Crack path tracking using DIC and XFEM modelling of mixed-mode fracture in wood. Theor. Appl. Fract. Mech. 2021, 112, 102896. [Google Scholar] [CrossRef]
  58. Sorensen, B.F.; Jacobsen, T.K. Determination of cohesive laws by the J integral approach. Eng. Fract. Mech. 2003, 70, 1841–1858. [Google Scholar] [CrossRef]
  59. Frossard, G.; Cugnoni, J.; Gmür, T.; Botsis, J. An efficient method for fiber bridging traction identification based on the R-curve: Formulation and experimental validation. Compos. Struct. 2017, 175, 135–144. [Google Scholar] [CrossRef]
  60. Kharratzadeh, M.; Shokrieh, M.; Salamat-talab, M. Effect of interface fiber angle on the mode I delamination growth of plain woven glass fiber-reinforced composites. Theor. Appl. Fract. Mech. 2018, 98, 1–12. [Google Scholar] [CrossRef]
  61. Gheibi, M.R.; Shojaeefard, M.H.; Googarchin, H.S. Direct determination of a new mode-dependent cohesive zone model to simulate metal-to-metal adhesive joints. J. Adhes. 2019, 95, 943–970. [Google Scholar] [CrossRef]
  62. Montenegro, D.; Pappas, G.; Botsis, J.; Zogg, M.; Wegener, K. 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. [Google Scholar] [CrossRef]
  63. Kottner, R.; Hynek, R.; Kroupa, T. Identification of parameters of cohesive elements for modeling of adhesively bonded joints of epoxy composites. Appl. Comput. Mech. 2013, 7, 137–144. [Google Scholar]
  64. Xu, Y.; Li, X.; Wang, X.; Liang, L. Inverse parameter identification of cohesive zone model for simulating mixed-mode crack propagation. Int. J. Solids Struct. 2014, 51, 2400–2410. [Google Scholar] [CrossRef] [Green Version]
  65. de Morais, A.; Pereira, A.; de Moura, M.; Silva, F.; Dourado, N. Bilinear approximations to the mixed-mode I–II delamination cohesive law using an inverse method. Compos. Struct. 2015, 122, 361–366. [Google Scholar] [CrossRef]
  66. Pincheira, G.; Ferrada, N.; Hinojosa, J.; Montecino, G.; Torres, L.; Saavedra, K. A study of interlaminar properties for a unidirectional glass fiber reinforced epoxy composite. Proc. Inst. Mech. Eng. Part L J. Mater. Des. Appl. 2019, 233, 348–357. [Google Scholar] [CrossRef]
  67. Joki, R.; Grytten, F.; Hayman, B.; Sørensen, B. Determination of a cohesive law for delamination modelling – Accounting for variation in crack opening and stress state across the test specimen width. Compos. Sci. Technol. 2016, 128, 49–57. [Google Scholar] [CrossRef] [Green Version]
  68. Abdel Monsef, S.; Ortega, A.; Turon, A.; Maimí, P.; Renart, J. An efficient method to extract a mode I cohesive law for bonded joints using the double cantilever beam test. Compos. Part B Eng. 2019, 178, 107424. [Google Scholar] [CrossRef]
  69. Skec, L. Identification of parameters of a bi-linear cohesive-zone model using analytical solutions for mode-I delamination. Eng. Fract. Mech. 2019, 214, 558–577. [Google Scholar] [CrossRef]
  70. Avalos, F.; Pincheira, G.; Inostroza, J.; Flores, P. Material parameter identification for vacuum infusion manufactured components. Int. J. Mater. Form. 2010, 3, 579–582. [Google Scholar] [CrossRef]
  71. Camacho, G.; Ortiz, M. Computational modelling of impact damage in brittle materials. Int. J. Solids Struct. 1996, 33, 2899–2938. [Google Scholar] [CrossRef]
  72. Alfano, G.; Crisfield, M.A. Finite element interface models for the delamination analysis of laminated composites: Mechanical and computational issues. Int. J. Numer. Methods Eng. 2001, 50, 1701–1736. [Google Scholar] [CrossRef]
  73. Stanzl-Tschegg, S.; Tan, D.M.; Tschegg, E. New splitting method for wood fracture characterization. Wood Sci. Technol. 1995, 29, 31–50. [Google Scholar] [CrossRef]
  74. Elices, M.; Planas, J. The equivalent elastic crack: 1. Load-Y equivalences. Int. J. Fract. 1993, 61, 159–172. [Google Scholar] [CrossRef]
  75. Allix, O.; Gosselet, P.; Kerfriden, P.; Saavedra, K. Virtual delamination testing through non-linear multi-scale computational methods: Some recent progress. Comput. Mater. Contin. 2012, 32, 107–132. [Google Scholar]
Figure 1. Double cantilever beam test.
Figure 1. Double cantilever beam test.
Polymers 13 01482 g001
Figure 2. Fracture process zone and equivalent LEFM: l f p z is defined as the distance between the tip of the stress-free crack a s f and the point along the potential crack path where damage begins (schema inspired from [26]).
Figure 2. Fracture process zone and equivalent LEFM: l f p z is defined as the distance between the tip of the stress-free crack a s f and the point along the potential crack path where damage begins (schema inspired from [26]).
Polymers 13 01482 g002
Figure 3. Traction-separation laws: (a) trilinear; (b) bilinear; and (c) potential CZM.
Figure 3. Traction-separation laws: (a) trilinear; (b) bilinear; and (c) potential CZM.
Polymers 13 01482 g003
Figure 4. Experimental DCB tests: layout according to [17].
Figure 4. Experimental DCB tests: layout according to [17].
Polymers 13 01482 g004
Figure 5. Experimental DCB tests: (a) load-displacement curve; (b) resistance curve.
Figure 5. Experimental DCB tests: (a) load-displacement curve; (b) resistance curve.
Polymers 13 01482 g005
Figure 9. Trilinear cohesive model. The dimensionless functions for model parameters as a function of the crack length: (a) w c function; (b) f t function.
Figure 9. Trilinear cohesive model. The dimensionless functions for model parameters as a function of the crack length: (a) w c function; (b) f t function.
Polymers 13 01482 g009
Figure 12. Bilinear cohesive model. The dimensionless functions for model parameters as a function of the crack length: (a) w c function; (b) K 0 function.
Figure 12. Bilinear cohesive model. The dimensionless functions for model parameters as a function of the crack length: (a) w c function; (b) K 0 function.
Polymers 13 01482 g012
Figure 13. Potential cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the initial stiffness K 0 with n = 0.7 (-) as constant.
Figure 13. Potential cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the initial stiffness K 0 with n = 0.7 (-) as constant.
Polymers 13 01482 g013
Figure 14. Potential cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the parameter n with K 0 = 1 · 10 13 N/m 3 as constant.
Figure 14. Potential cohesive model: (a) load-displacement curve; (b) equivalent R-curve. Influence of the parameter n with K 0 = 1 · 10 13 N/m 3 as constant.
Polymers 13 01482 g014
Figure 15. Potential cohesive model. Influence of (a) n and (b) K 0 on f t (more details in Table A1).
Figure 15. Potential cohesive model. Influence of (a) n and (b) K 0 on f t (more details in Table A1).
Polymers 13 01482 g015
Figure 17. Potential cohesive model. The dimensionless functions for model parameters as a function of the crack length: (a) w c function; (b) K 0 function.
Figure 17. Potential cohesive model. The dimensionless functions for model parameters as a function of the crack length: (a) w c function; (b) K 0 function.
Polymers 13 01482 g017
Figure 18. Schema of the optimization procedure.
Figure 18. Schema of the optimization procedure.
Polymers 13 01482 g018
Figure 19. Trilinear cohesive model. Numerical DCB test with the fitted parameters: (a) traction-separation law; (b) load-displacement curve; (c) 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.
Polymers 13 01482 g019
Figure 20. Bilinear 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.
Polymers 13 01482 g020
Figure 21. Potential 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.
Polymers 13 01482 g021
Figure 22. Experimental and numerical DCB tests (3D and 2D simulations): (a) load-displacement curve; (b) R-curve.
Figure 22. Experimental and numerical DCB tests (3D and 2D simulations): (a) load-displacement curve; (b) R-curve.
Polymers 13 01482 g022
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.
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.
Polymers 13 01482 g023
Table 1. DCB specimen geometry according to [17].
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.0720.4 ± 0.08124.68 ± 0.5247
Table 2. Experimental DCB tests: critical energy release rate and maximal applied load.
Table 2. Experimental DCB tests: critical energy release rate and maximal applied load.
Specimen G c exp (N/m) P max exp (N)
1981.6531.08
2926.6032.16
31062.8029.72
4964.9027.98
5974.830.70
average value982.1530.33
standard deviation49.851.58
Table 3. Damage functions d ( w ) and parameters of the bilinear, trilinear and potential CZM.
Table 3. Damage functions d ( w ) and parameters of the bilinear, trilinear and potential CZM.
CZMTrilinear [26]Bilinear [72]Potential [40]
d = 0 , w < w 0 d = 0 , w w 0 d = m i n n n + 1 Y G c n , 1
d = w b ( w w o ) ( 1 γ ) w ( w b w o ) , w o w w b d = w c w c w 0 w w 0 w , w 0 < w w c where Y = 1 2 K 0 w 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 b w o f t w b , w 0 = f t K 0 , where w 0 = f t K 0 , f t = 2 G c w c
w b = 2 G f μ f t , f b = 2 G f b w c
parametersto be identified w c , G f μ / G c , f t w c , K 0 K 0 , n
Table 4. Parameters for the genetic algorithm.
Table 4. Parameters for the genetic algorithm.
poblation size5000
crossover probability0.7
mutation probability0.1
number of generation300
number of couples500
pressure0.05
Table 5. Fitted parameters for each cohesive model using equivalent LEFM and genetic algorithm.
Table 5. Fitted parameters for each cohesive model using equivalent LEFM and genetic algorithm.
CZMCharacterizationNumber of DCB Virtual Tests Π ( Γ CZM fitt ) Δ a c (mm)Fitted Parameters Γ CZM fitt
trilineareq. LEFM R-curve100.3315.0 f t = 11.51 MPa
w c = 2.71 mm
G f u / G c = 0.83 (-)
genetic algorithm360.4614.5 f t = 10.08 MPa
w c = 2.27 mm
G f u / G c = 0.85 (-)
bilineareq. LEFM R-curve1016.215.0 w c = 2.62 mm
K 0 = 2.45 · 10 12 N/m 3
eq. LEFM R-curve103.212.5 w c = 0.064 mm
K 0 = 8.67 · 10 12 N/m 3
genetic algorithm361.674.96 w c = 0.32 mm
K 0 = 1.151 · 10 12 N/m 3
potentialeq. LEFM R-curve1032.215.0 n = 2.1 · 10 4 (-)
K 0 = 1.12 · 10 12 N/m 3
eq. LEFM R-curve102.982.5 n = 0.86 (-)
K 0 = 1 · 10 12 N/m 3
genetic algorithm362.224.3 n = 0.1 (-)
K 0 = 4.71 · 10 11 N/m 3
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Torres, L.; Saavedra, K.; Pincheira, G.; Pina, J.C. 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. Polymers 2021, 13, 1482. https://doi.org/10.3390/polym13091482

AMA Style

Torres L, Saavedra K, Pincheira G, Pina JC. 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. Polymers. 2021; 13(9):1482. https://doi.org/10.3390/polym13091482

Chicago/Turabian Style

Torres, Luis, Karin Saavedra, Gonzalo Pincheira, and Juan Carlos Pina. 2021. "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" Polymers 13, no. 9: 1482. https://doi.org/10.3390/polym13091482

APA Style

Torres, L., Saavedra, K., Pincheira, G., & Pina, J. C. (2021). 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. Polymers, 13(9), 1482. https://doi.org/10.3390/polym13091482

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop