Delft University of Technology
Cohesive zone and interfacial thick level set modeling of the dynamic double cantilever
beam test of composite laminate
Liu, Y.; van der Meer, F. P.; Sluys, L. J. DOI
10.1016/j.tafmec.2018.07.004 Publication date
2018
Document Version
Accepted author manuscript Published in
Theoretical and Applied Fracture Mechanics
Citation (APA)
Liu, Y., van der Meer, F. P., & Sluys, L. J. (2018). Cohesive zone and interfacial thick level set modeling of the dynamic double cantilever beam test of composite laminate. Theoretical and Applied Fracture
Mechanics, 96, 617-630. https://doi.org/10.1016/j.tafmec.2018.07.004 Important note
To cite this publication, please use the final published version (if applicable). Please check the document version above.
Copyright
Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy
Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.
Cohesive zone and interfacial thick level set modeling of the dynamic double
cantilever beam test of composite laminate
Y. Liu, F.P. van der Meer and L.J. Sluys
Faculty of Civil Engineering and Geosciences, Delft University of Technology PO Box 5048, 2600 GA Delft, The Netherlands
Email: y.liu-7@tudelft.nl
Abstract
The mode-I interlaminar fracture toughness of composite laminates under different loading rates can be measured by the double cantilever beam (DCB) test. It is observed from the DCB test of a unidirectional PEEK/carbon composite laminate that as the loading rate increases from quasi-static to dynamic range: (1) delamination crack growth exhibits a transition from stable to unstable (“stick/slip”) and back to a stable type; (2) the interlaminar fracture toughness is not constant as the loading rate increases. In this paper, two numerical approaches are used to reproduce the experimental observations: a cohesive zone model (CZM) and the interfacial thick level set (ITLS) model. CZM simulations with rate-independent and rate-dependent cohesive laws are carried out. A new version of the ITLS is introduced with a phenomenological relation between crack speed and energy release rate. The simulation results of the CZM and the ITLS model are compared with the real DCB test data to evaluate the capability of these two types of models. It is found that the used CZM can reproduce rate-dependence of the fracture energy, but not the stick/slip behavior. The ITLS can capture the stick/slip behavior, but needs different parameter sets for different loading rates.
Keywords: Double cantilever beam, Rate dependency, Crack arrest, Cohesive zone model, Thick level set
1. Introduction
Delamination is one of the crucial degradation mechanisms of composite laminates [1]. Engineering composite laminates can be subjected to complex working load conditions including quasi-static and dynamic loading (e.g. low velocity impact). In order to predict the extent of delamination under given load conditions, it is important to quantify the interlaminar fracture toughness of composite laminates for both quasi-static and dynamic loading.
The double cantilever beam (DCB) test is one of the most commonly used experimental methods for determining the mode-I interlaminar fracture toughness [2]. Studies on measurement of the mode-II, mode-III interlaminar fracture toughness can be found in [3–8]. Based on published results on the DCB test with a smallest test rate of 1.67×10−7
m/s in [9] and a largest test rate of 15 m/s in [10, 11], the following observations can be made. Firstly, depending on the investigated test rate and composite system, the crack propagation in the DCB test could be either stable or unstable (“stick/slip”) [3, 9, 10, 12–14]. In [9], where a carbon/epoxy composite material, T300/2500, was tested within the crosshead speed range of 1.67×10−7 to 8.33×10−3 m/s, the delamination crack growth was unstable for test rates lower than 8.33 ×10−6 m/s and became stable for higher test rates. In [12], a unidirectional and a woven carbon/epoxy composite laminate with 19% inclusion of transversal E-glass fibres was tested with different crosshead velocities ranging from 8.3 × 10−5m/s to 0.19 m/s. For both materials and all loading rates, the crack propagated in an unstable fashion. Secondly, the loading rate may influence the fracture toughness although there is no universal trend in increased loading rate effects on the toughness. Aliyu and Daniel [15] investigated an AS-4/3501-6 carbon/epoxy system with DCB tests at a crosshead displacement rate up to 8.5× 10−3m/s. It was found that the mode-I interlaminar
fracture toughness increased 28% over three orders of magnitude of loading rate. Smiley and Pipes [13] found that the mode-I interlaminar fracture toughness for APC-2 carbon/PEEK composite remains constant over four decades of low loading rates while further increase in loading rate caused a decrease up to 70% over the next decade of loading rate. This inconsistency in the trends might be attributed to differences in material constituent, specimen geometry, data-reduction scheme, measurement technique and definition of the rate parameter [16, 17].
© 2018 Manuscript version made available under CC-BY-NC-ND 4.0 license
https://creativecommons.org/licenses/by-nc-nd/4.0/
In terms of numerical modeling, the cohesive zone model (CZM) is widely used for modeling delamination. Mo-tivated by the possible rate-sensitivity of composite laminates (or polymeric adhesives), rate-dependent cohesive laws have been proposed. They can be roughly divided into four categories, namely, the dynamic increase factor (DIF) models, the damage-delay models, the viscoplasticity models and the viscoelasticity models. The first category of models assumes that some characteristic parameters (e.g. the fracture energy) of a chosen cohesive law are functions of the separation rate of the cohesive surface [18–24]. The second category modifies the classical damage evolution formulation into a rate form, thus limiting the damage rate to a certain maximum [7, 25, 26]. The third category intro-duces an overstress function to the rate-independent softening plasticity law [27–29]. The fourth category combines viscoelastic components (e.g. Maxwell element) with a damage-type of rate-independent cohesive law [30–32].
Rather than treating damage along the interface as a result of a local displacement jump as it is done for the CZM model, the damage can be driven by a non-local energy release rate. This approach is possible with the interfacial thick level set (ITLS) model. The thick level set (TLS) model was first introduced by M¨oes et al. [33] for the modeling of damage growth in a continuum under quasi-static loading conditions. Contrary to conventional continuum damage mechanics, the damage d in the TLS is not a direct function of the local strain but rather a function of a level set field φ, whose definition is the signed distance to a moving damage front. This damage band has a predefined fixed length acting as an intrinsic length scale to avoid spurious localization and the associated pathological mesh-dependency observed in local damage models. Recent developments and applications of the TLS method include an improved numerical implementation scheme [34], 3D crack simulation for quasi-brittle materials [35], shear failure simulation of sandwich structures [36] and simulation of fragmentation under impact [37, 38]. The TLS model was translated to interface element by Latifi et al. [39] to provide an alternative to the CZM for interfacial cracking. The damage in this interfacial thick level set (ITLS) model is a function of a level set field defined on the interface. The ITLS allows for non-local evaluation of the energy release rate. Consequently, the damage growth criterion can be related to a fracture mechanics based crack growth criterion. The ITLS model, in this sense, is advantageous since it provides a robust numerical tool for straightforward implementation of energy release rate based crack growth relations. The application of the ITLS model for simulating fatigue crack growth, where the Paris relation gives such a relation, was reported in [40, 41].
This paper is organized as follows: in Section 2, a rate-independent and a rate-dependent cohesive law are in-troduced and their ability to reproduce experimental observations is evaluated. The formulation, solution scheme and relevant details of the ITLS model are presented in Section 3. The capability of this method is demonstrated by modeling the dynamic DCB test with a parametric study in Section 4. The final section gives a short comparison and summary about the CZM model and the ITLS model.
2. Cohesive zone modeling
In this section, both rate-independent and rate-dependent cohesive modeling are used to simulate the DCB tests on unidirectional PEEK/carbon composite laminate reported by Blackman et al. [10, 11]. Starting point is the cohesive law by Turon et al. [42]. Rate-dependency is introduced for the cohesive strength and fracture energy following May [21, 22]. By setting the rate sensitivity parameters in the rate-dependent cohesive law to be zero, the cohesive law introduced by Turon et al. [42] is recovered.
2.1. Rate-independent cohesive law
For the sake of completeness, the rate-independent cohesive law is first reviewed. The studied DCB test is assumed to be a plane strain problem, so only a two-dimensional formulation is introduced below. The following relation between traction tiand displacement jump ~uiis used:
ti= (1 − d) Kδi j~uj− dKh−~u1iδ1i, i= 1, 2; j = 1, 2 (1)
where index 1 represents the normal direction, 2 is the shear direction, d is a scalar damage variable, K is a penalty stiffness, δi jis the Kronecker delta and the MacAuley bracket is defined as hxi=12(x+ |x|).
The displacement jump when damage initiates, ~u0i, is given as,
~u0i =
σi
~u2 t ~u1 (~u0 2, σ2) (~u2f, 0) (~u01, σ1) (~u1f, 0) (∆f, 0) (∆0, σm) Mode-II Mode-I
Fig. 1. Mixed-mode cohesive law
in which σiis the cohesive strength. The displacement jump when full decohesion occurs, ~u f i, reads ~uif = 2GIc σi , i= 1 2GIIc σi , i= 2 (3)
where GIcis the mode-I fracture energy and GIIcis the mode-II fracture energy.
The coupling between independent modes is demonstrated in Fig. 1. For a given mixed-mode ratio defined as
β = |~u2|
h~u1i+ |~u2|
(4)
The mixed-mode crack initiation criterion is introduced as∆ = ∆0, in which the equivalent displacement jump∆ is
∆ = q(h~u1i)2+ (|~u2|)2 (5) and, ∆0= r ~u01 2 + ~u02 2 −~u0 1 2 Bη (6)
where η is a mode interaction coefficient and
B= β
2
1+ 2β2− 2β (7)
The crack propagation criterion is imposed by∆ = ∆f where
∆f = ~u01~u f 1+ ~u02~u f 2− ~u 0 1~u f 1 Bη ∆0 (8)
and the mixed-mode fracture energy Gcis calculated as
The damage at time T is defined as, d(T )= max τ≤T 0, ∆ ≤ ∆0 ∆f(∆−∆0) ∆(∆f−∆0), ∆0 < ∆ < ∆f 1, ∆ ≥ ∆f (10)
Detailed derivations on the above relations can be found in [42].
2.2. Rate-dependency relation
Rate dependency is introduced for both the cohesive strength and the fracture energy of independent modes with a Johnson-Cook law similar to May [21, 22]. For pure mode-I or mode-II, the following relation is introduced,
σi( ˙∆) = σ0 i 1+ ciln ˙ ∆i ˙ ∆ref i , ˙∆i≥ ˙∆refi σ0 i, ∆˙i< ˙∆refi , i= 1, 2 (11)
where σiis the rate-dependent cohesive strength, σ0i is the quasi-static cohesive strength, ciis a rate-sensitivity
con-stant, ˙∆i is the displacement jump rate defined as ˙∆1 = ∂h~u1i/∂T and ˙∆2 = ∂|[u]2|/∂T , and ˙∆refi is a reference
displacement jump rate.
The rate dependency of the fracture toughness is assumed to be,
GIc( ˙∆) = G0 Ic, ∆˙1< ˙∆ ref 1 G0Ic 1+ m1ln ˙ ∆1 ˙ ∆ref 1 , ˙∆ref 1 ≤ ˙∆1≤ ˙∆ inf 1 GinfIc, ∆˙1> ˙∆inf1 (12) GIIc( ˙∆) = G0 IIc, ∆˙2< ˙∆ ref 2 G0IIc 1+ m2ln ˙ ∆2 ˙ ∆ref 2 , ˙∆ref 2 ≤ ˙∆2≤ ˙∆ inf 2
GinfIIc, ∆˙2> ˙∆inf2
(13)
where G0 Ic and G
0
IIc are the mode-I and mode-II fracture energy under quasi-static loading, m1 and m2 are
rate-sensitivity constants, ˙∆inf 1 and ˙∆
inf
2 are the reference displacement rates for defining constants G inf Ic and G
inf IIc that
are introduced to bound the fracture energy when negative rate-dependency parameters are used such that Ginf Ic = G0 Ic 1+ m1ln ∆˙inf 1 ˙ ∆ref 1 and Ginf IIc = G 0 IIc 1+ m2ln ∆˙inf 2 ˙ ∆ref 2
. It is further assumed that when rate dependency is activated, the coupling between independent modes introduced in the previous section is still valid. An advantage of the phe-nomenological DIF formulation in Eq. (12) and (13) is that it can accommodate both a positive and negative influence of displacement jump rate on either cohesive strength or fracture energy while the existing rate-dependent CZMs de-veloped with damage-delay [7, 25, 26] or viscosity formulation [27–32] can only capture a positive influence. This versatility provides more possibilities in describing different failure mechanisms in different composite material sys-tem and loading conditions. The measurements by Blackman et al. [10, 11] which show a decrease in fracture energy for increasing loading rate are exemplary of this advantage of the chosen cohesive formulation.
2.3. Model description
Numerical models have been created to simulate both dynamic and quasi-static DCB tests. Depending on the type of loading, i.e. dynamic or quasi-static, the numerical model set-ups have differences.
The geometry and boundary conditions of the numerical model used for simulating dynamic DCB tests are de-picted in Fig. 2. Load is applied at a point on the top surface, the velocity of which is raised to a constant test rate
˙
wwithin a short time and maintained until the end of the simulation. The vertical degree-of-freedom of the node symmetric to the loading point with the crack plane is constrained. The bulk material of the beam is represented by plane strain triangular elements with an orthotropic linear elastic constitutive model. The initially intact section of the middle surface between the two arms is modeled by zero-thickness cohesive zone elements with the cohesive
˙ w
L h
a0 b
Fig. 2. Geometry and boundary conditions of the double cantilever beam test under dynamic loading
law introduced in the previous section. The precracked section of middle surface is considered by zero-thickness interface elements with penalty stiffness in case of a negative ~u1to enforce the contact condition. Implicit dynamic
analysis is performed with a Newmark time integration scheme. When the rate-dependent cohesive law is used, the displacement jump rate in Eqs. (11), (12) and (13) is calculated with the Euler-backward scheme, which means that the displacement jump rate is implicitly updated at the end of each load step. The derivation of the consistent tangent is given in Appendix A.
For quasi-static DCB tests, the numerical model differs from the model used for dynamic loading at three points. Firstly, the horizontal degree-of-freedom of the loading point and supported point is constrained to eliminate the rigid body motion of the system. Secondly, the loading velocity is instantly set to the constant test speed ˙wat the start of the simulation. Thirdly, no time integration scheme is needed as inertia is not considered.
2.4. Material parameters
The experiments from [10, 11] were preformed on a unidirectional PEEK/carbon composite laminate (APC-2) with a nominal fiber volume fraction of 65%. The bulk material is assumed to be orthotropic linear elastic with parameters as defined in Table 1. Of these, the longitudinal Young’s modulus E1, density ρ and a major Poisson’s
ratio ν12are taken from the measured value in [10] while the other values which are taken from typical values in [43]
have limited influence on the response.
The parameters needed to be determined for the rate-dependent cohesive law can be divided into: parameters related to the rate-independent cohesive law and parameters related to the rate-dependency. Considering that the simulated DCB test is mode-I dominated, the material properties for mode-I and mode-II are assumed to be the same. The quasi-static fracture energies, G0
Icand G 0
IIc, are determined to be 1.5 N/mm from the range of measurement values
shown in Fig. 12 in [10]. Typical values are used for σ0 1, σ
0
2and η. A high value is prescribed for the penalty stiffness
K.
Concerning rate-dependency, the quasi-static reference displacement jump rates, ˙∆ref 1 and ˙∆
ref
2 , determine the
threshold value for the rate effect comes into play. Several trial simulations of a quasi-static DCB test, with a loading rate of 3.3 × 10−5m/s, using different reference displacement jump rates are performed and the results are compared
with the rate-independent cohesive model. The smallest value adopted in the rate-dependent CZM model that results in the same response as the rate-independent CZM model is determined to be the quasi-static reference displacement jump rates. Finally, the limit jump rate ˙∆inf
1 and ˙∆ inf
2 are set sufficiently high not to be reached in the presented
simula-tions. They are included in the formulation only to exclude the possibility of negative fracture energy from the model. The calibrations of c1, c2, m1and m2are done by trial simulations to find the best match with the time vs. crack length
of the studied four DCB tests. A summary of CZM related parameters is shown in Table 2. The rate-dependency parameters for cohesive strength are found to be positive. This means that the cohesive strength increases with the displacement jump rate, which could be due to the viscosity of the PEEK polymer matrix [13, 44, 45]. The rate-dependency parameters for fracture energy are negative, meaning that a decrease of fracture energy with increasing displacement jump rate is found. This increased brittleness is possibly due to lack of time for plasticity to develop. The strain rate dependency of a brittle crack results in a smaller dissipated energy with larger strain rate [45].
2.5. Simulation results
The numerical models are used to simulate four DCB tests of a unidirectional PEEK/carbon composite laminate with different test rates reported in [10, 11]: 3.3 × 10−5, 1.0, 6.5 and 10.0 m/s. Both the independent and
rate-Table 1: Bulk material parameters: E1, ρ and ν12from [10], G12, E2and ν23from [43]
elastic modulus E1= 115000 N/mm2 E2= 8000 N/mm2
Poission’s ratio ν12= 0.28 ν23= 0.43
shear modulus G12= 5000 N/mm2
density ρ = 1540 Kg/m3
Table 2: CZM related properties
Quasi-static cohesive strength (mode-I/II) σ01= σ02= 60 N/mm2 Quasi-static fracture energy (mode-I/II) G0Ic= G0IIc= 1.5 N/mm
Interface penalty stiffness K= 100000 N/mm3 Mode interaction coefficient η = 2.284 Quasi-static reference displacement jump rate (mode-I/II) ∆˙ref
1 = ˙∆ ref
2 = 0.02 m/s
Strength coefficient (mode-I/II) c1= c2= 0.22
Fracture energy coefficient (mode-I/II) m1= m2= −0.0783
Maximum reference displacement jump rate (mode-I/II) ∆˙inf 1 = ˙∆
inf
2 = 100.0 m/s
Limit fracture energy (mode-I/II) Ginf Ic = G
inf
IIc= 0.5 N/mm
dependent cohesive law are considered for each analysis. The dimensions for the first three test cases are the same, i.e. L= 200.0 mm, h = 1.5 mm, a0= 35.0 mm, b = 7.1 mm (Fig. 2), while the dimensions for the last test case only differ
in its initial crack length, i.e. a0= 30 mm. The time vs. crack length curve for the four simulated cases are plotted in
Fig. 3a-3d. It can be observed that both the rate-independent (denoted in the figure with c= m = 0) and rate-dependent cohesive model predict a continuous crack propagation process for all the analysis. The model fails to capture the crack arrest phenomenon that was observed experimentally at the test rates of 1.0 and 6.5 m/s. It can be seen from Fig. 3a that the rate-independent model and rate-dependent model produce almost the same time vs. crack length curve for the lowest loading rate. They both match reasonably well with the measured curve. For the 1.0 m/s case (see Fig. 3b), the rate-dependent model produces a time vs. crack length curve which is almost coincident with the left boundary of the test curve. The average crack speed from the calibrated rate-dependent model matches with the measurement and the only difference is that in the test the crack does not propagate continuously. The rate-independent cohesive model matches the right boundary of the test curve and the overall crack propagation is slower than the measurement. For the 6.5 m/s case (see Fig. 3c), the time vs. crack length curve obtained with the rate-dependent cohesive model follows the measurement well except that it can not capture the crack arrest event for a crack length of around 81 mm. For this loading rate, the time vs. crack length curve produced by the rate-independent cohesive model has a relatively large deviation from the experimental curve. The model predicts a later crack initiation and the crack speed is slower. For the test rate of 10.0 m/s (Fig. 3d), the time vs. crack length curve obtained with the rate-dependent model still matches reasonably well with the measurement in terms of crack speed, while the rate-independent model significantly underpredicts the crack speed. In conclusion, the rate-dependent model can reproduce the measurements except for the crack arrests which are not captured. The mismatch between measurements and rate-independent model result indicates the significance of rate dependency for this series of experiments.
3. Interfacial thick level set modeling
As an alternative to the CZM, the ITLS model proposed by Latifi et al. [39] is introduced here to simulate the dynamic DCB tests. Similar to the CZM, a traction-separation relation with damage is assumed as the constitutive model of the interface. However, unlike the CZM, damage in the ITLS is not a function of the displacement jump across the surface but a function of a level set field φ on the interface. The level set function is defined as a signed distance function from a damage front. In 2D the damage front collapses into a point and φ gives the signed distance along the line of the interface to this point. The damage band has a fixed length lcover which damage d increases
from 0 to 1 as the distance to the front grows from 0 to lc. This band separates uncracked material with d = 0 from
0 200 400 600 800 1,000 1,200 40 60 80 100 Time (s) Crack length (mm) Experiment ci= 0, mi= 0 ci= 0.22, mi= −0.0783 (a) ˙w= 3.3 × 10−5m/s 0 10 20 30 40 40 60 80 100 Time (ms) Crack length (mm) Experiment ci= 0, mi= 0 ci= 0.22, mi= −0.0783 (b) ˙w= 1.0 m/s 0 1 2 3 4 5 6 40 60 80 100 Time (ms) Crack length (mm) Experiment ci= 0, mi= 0 ci= 0.22, mi= −0.0783 (c) ˙w= 6.5 m/s 0 0.5 1 1.5 2 2.5 3 40 60 80 100 Time (ms) Crack length (mm) Experiment ci= 0, mi= 0 ci= 0.22, mi= −0.0783 (d) ˙w= 10.0 m/s
Fig. 3. Comparison between experimental and cohesive zone modeling result for test rate of (a) ˙w= 3.3 × 10−5m/s; (b) ˙w = 1.0 m/s; (c) ˙w = 6.5 m/s; (d) ˙w = 10.0 m/s.
be computed. Moreover, a crack growth rate can be imposed as the front velocity. This makes the method suitable for implementation of crack growth relations that relate crack growth rate to energy release rate, such as the Paris law for fatigue crack growth [40, 41]. In this section, the formulation and solution scheme for the ITLS method are described.
3.1. Local governing equation
Following Latifi et al. [39], the free energy per unit surface of interface is defined as:
ψ(~u, d) = (1 − d) ψ0(~u) +
1
2dKh−~u1i
2 (14)
where ~u = (~u1, ~u2)Tis the displacement jump vector between the two facets of the interface, ~u1is the normal
displacement jump, ~u2is the shear displacement jump and d is a scalar damage variable. The variable ψ0is defined
as:
ψ0(~u) =
1
2~uiKδi j~uj (15)
where K is a penalty stiffness and δi jis the Kronecker delta.
The relation between the traction t and the displacement jump ~u across the surface is obtained by:
ti=
∂ψ ∂~ui
= (1 − d) K~ui− dKh−~u1iδ1i, i= 1, 2 (16)
The local driving force for damage growth is obtained by differentiating the free energy function with respect to the damage variable:
Y= −∂ψ ∂d =ψ0(~u) − 1 2Kh−~u1i 2 (17) 3.2. Damage definition
Fig. 4 shows the damage band, the level set field φ defined on the interface and its associated damage distribution. The damage function proposed by Voormeeren et al. [41] is adopted.
d(φ)= ( 0, φ < 0 arctanγlφ c arctan−1(γ), 0 < φ < l c 1, φ > lc (18)
where γ is a shape constant. For the current study a value of γ= 15 is used following Voormeeren et al. [41]. The adopted damage function is plotted in Fig. 5.
3.3. Energy release rate
Because the damage distribution inside the band is given, the energy release rate for crack growth or movement of the band can be calculated by integrating the local driving force over the damage band. In 2D, the integration takes place along a line:
G= − Z lc 0 ∂ψ ∂φdφ = − Z lc 0 ∂ψ ∂d ∂d ∂φdφ = Z lc 0 Y∂d ∂φdφ (19)
A two-point Gauss integration scheme is used for elementwise numerical integration of Eq. (19). However, during the movement of the damage band, its front (φ= 0) and wake (φ = lc) can intersect the interface element causing
an unsmooth damage distribution inside the element. In order to improve the accuracy of numerical integration for Eq. (19) with discontinuous ∂d∂φ, the ITLS element crossed by the damage band front or wake is subdivided into two sub-elements. In each of the two elements, two-point Gauss integration is applied (Fig. 6).
lc φ(x) x y φ ≤ 0 0 ≤ φ ≤ lc φ ≥ lc damage band uncracked cracked
damage front crack tip
φ = 0 φ = lc
d
x d= 0
0 ≤ d ≤ 1 d= 1
Fig. 4. Initial level set field and damage distribution: The level set field at one point (red point) is defined as the signed distance to the damage band front (blue point).
−0.50 0 0.5 1 1.5 0.2 0.4 0.6 0.8 1 1.2 φ/lc d
Fig. 5. Damage function d(φ) with γ= 15
φ ≤ 0, d = 0 undamaged zone
0 ≤ φ ≤ lc, 0 ≤ d ≤ 1
damaged zone
φ ≥ lc, d= 1
fully damaged zone
× × × ×
subdivision
φ = 0
× × × × φ = lc
zero-thickness ITLS element × Gauss integration point
Ga Gi Va Vi Vm Arrest Initiation Propagation
Energy release rate G
Crack
speed
V
Fig. 7. Relation between crack velocity V and energy release rate G (the dashed line schematically represents the shape for larger crack speed range [18, 50])
3.4. Crack speed
Completing the ITLS formulation requires a relation between the crack speed V and the energy release rate G. This can be an advantage. For instance, for fatigue crack growth experimental observations can be described well with the Paris law, which gives such a relation. In [40, 41], the ITLS method has been validated for the calculation of energy release rate and prediction of fatigue crack growth. In this paper, a function between energy release rate G and crack speed V for dynamic crack growth with possible stick/slip behavior is introduced (see Fig. 7). This function is inspired by the relation between the dynamic stress intensity factor K and crack speed V proposed by Ravi-Chandar [46] for describing dynamic crack growth. Similar to Ravi-Chandar [46], we explicitly differentiate three states for dynamic crack growth, namely, crack initiation, propagation and arrest. The crack starts to grow when the energy release rate G reaches the crack initiation toughness Gi. Crack growth starts at a nonzero crack speed Vi. During crack
growth, the crack speed V is related to the energy release rate G according to a phenomenological relation which possesses three features: (1) the crack speed has an asympototic maximum value Vm[46, 47]; (2) the V(G) curve has
a positive slope, representing the influence of increased microcracking on the fracture toughness [7, 46, 48]; (3) the crack speed jumps from a finite value Vato zero when the energy release rate drops below the crack arrest toughness
Ga[49]. The crack initiation point is not on the curve characterizing the dynamic crack growth criterion is possibly
due to bluntness of the initial crack, intrinsic rate dependence of the material, or inertia effects [46]. This difference between crack initiation and growth toughness causes the crack growth rate to jump to a finite speed immediately upon initiation. This behaviour has been observed experimentally for crack growth in polymers [49].
For the studied DCB tests, the Rayleigh wave speed cRis the theoretical crack speed limit Vm[47]. In the
exper-iments the crack speed is smaller than 85.0 m/s which is far from the theoretic limit (cR = 2523 m/s). Therefore, we
introduce a simple linear relation for the studied range of tests as (see Fig. 7)
V= Vi+
Vi− Va
Gi− Ga
(G − Gi) (20)
3.5. Solution scheme
The numerical model is created to simulate the delamination crack in dynamic DCB tests. The model is the same as mentioned in Section 2.2 except that the constitutive behavior of the cohesive elements inserted along the middle surface of the beam is now described with the ITLS.
The program flow is illustrated in Box 1. There are a few items to be noted: (1) a damage band of length lcis
predefined ahead of the initial crack tip before the first load increment is applied; (2) the Newmark time integration scheme is used for computing displacements for a given damage distribution; (3) a staggered solution scheme is
1. Set the time step number n= 1;
2. Define a initial damage band with length lcahead of the crack tip and calculate the level set field φ(x);
3. Initialize damage distribution based on the level set field φ(x) using Eq. (18); 4. Set the initial time step δT1;
5. Apply the nth load/displacement increment;
6. Solve dynamic equilibrium with the Newton-Raphson scheme;
7. Calculate the energy release rate of the nth time step, Gn, using Eq. (19);
8. Check the crack state at (n − 1)th time step, stationary or propagating?
• Stationary—The crack speed at the end of the nth time step, Vn, is calculated with Eq. (20) if the energy
release rate Gnis larger than the initiation toughness Gi, otherwise it is zero.
• Propagating—The Vnis calculated with Eq. (20) if Gnis larger than the arrest toughness Ga, otherwise it is
zero.
9. Set time step size for the (n+1)th time step: δTn= min(δT1, 0.6l
ie/Vn), where lieis the minimum interface element
size;
10. Move the level set field in the crack growth direction by a distance of δsnwith δsn = Vn−1·δTn(see Fig. 8);
11. Update the damage distribution with the updated level set field using Eq. (18); 12. n= n + 1, go to step 5;
Box 1: Staggered solution algorithm for the ITLS model
X Y Z φ < 0, d = 0 φ > lc, d= 1 0 ≤ φ ≤ lc, 0 ≤ d ≤ 1 T = Tn T = Tn−1 δsn
crack growth direction
Fig. 8. Movement of damage band within the nth time step (from time Tn−1to Tn): triangular meshes are used for discretization.
used for the damage update and the displacement field because of its robustness and simplicity. The translational distance of the damage band at the end of time step n is determined by the speed at the time step n − 1. Therefore, when solving the dynamic equilibrium equation with the Newton-Raphson scheme, the damage is only updated when global convergence is achieved through the update of the level set field.
4. Simulation of dynamic DCB tests
The ITLS model is used to simulate the same series of DCB tests discussed in Section 2.4 except the quasi-static case with a test rate of 3.3 × 10−5m/s, since the V(G) relation introduced in Section 3.4 is a dynamic crack growth
criterion assumed for delamination crack growth in the context of dynamic loading. A V(G) relation with artificial viscoscity for simulating delamination crack generated by quasi-static loading with the ITLS has been introduced in [39].
In this section, the accuracy of the ITLS model for calculating the energy release rate is first validated. Afterwards, the influence of parameters used for the dynamic crack growth criterion is discussed and then the ITLS model is used to calibrate the material parameters for the dynamic DCB tests.
4.1. Verification of energy release rate computation
The accuracy of the energy release rate calculation in the ITLS model is demonstrated by modeling a quasi-static DCB specimen in which the deflection of the beam w is fixed while the crack length is gradually increased. Bulk
30 40 50 60 70 80 90 100 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 ITLS A ITLS B MBT Ener gy release rate G (N /mm) Crack length a (mm)
Fig. 9. Evolution of energy release rate with crack length for ITLS A (ITLS model without element subdivision used), ITLS B (ITLS model with subdivision applied) and MBT (modified beam theory solution).
material parameters listed in Table 1 are used. The value of 15.0 is adopted for the shape constant γ of the damage function. The damage band width lcis chosen as 0.9 mm in the current study, which is a typical cohesive zone length
found in the cohesive modeling study shown in Section 2. The interface element size is around 0.125 mm, which ensures 8 elements inside the damage band and a high accuracy in calculating the energy release rate. The interface penalty stiffness K = 100000 N/mm3. The energy release rate from Eq. (19) is compared with the value from modified
beam theory (MBT) introduced in [51]. According to the MBT, the energy release rate for the quasi-static DCB test can be calculated by
G= 3 4
E1h3w2
(a+ χh)4 (21)
where the correction factor χ is estimated by
χ = E1 11G12 ! 3 − 2 Γ 1+ Γ !2 1/2 (22) in whichΓ = 1.18E1E2 G12 1/2
. Fig. 9 shows the crack length vs. energy release rate G calculated by the ITLS model with element subdivision and without element subdivision technique introduced in Section 3.3 along with the solution given by the MBT. It is shown that the ITLS model with element subdivision removes, to a large extent, the oscillations that are present for the model without element subdivision. The result by the ITLS model with element subdivision has a good match with the MBT solution, which demonstrates that the ITLS model is capable of calculating the energy release rate accurately.
In order to investigate the influence of the integration scheme on the accuracy of calculating G, a relative error is introduced as
=Z 55
35
|G(a) − G(a)|
G(a) da (23)
in which G(a) is the calculated energy release rate G at corresponding crack length a and G(a) is a reference exact solution. Three different mesh sizes are used for the ITLS model, changing from the coarsest mesh with ITLS interface element size of 0.125 mm, to medium mesh with a size of 0.0625 mm, and to the finest mesh with a size of 0.03125 mm. For each mesh size, four simulations, i.e., ITLS A with 2-point Gauss, ITLS A with 3-point Gauss, ITLS B with 2-point Gauss and ITLS B with 3-point Gauss, are performed and for each simulation we can calculate the relative
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 0.2 0.4 0.6 0.8
ITLS element size
Relati v e error Gauss2 (ITLS A) Gauss3 (ITLS A) Gauss2 (ITLS B) Gauss3 (ITLS B)
Fig. 10. Error estimation of the integration scheme
error . In this case, the numerical solution obtained from the ITLS model with interface element size of 0.03125 mm and with 3-point Gauss integration scheme and element-subdivision (ITLS B) is chosen as the reference exact solution G(a) used in Eq. (23). This result is plotted in Fig. 10 with the horizontal axis showing the number of the elements and the vertical axis showing the relative error . It can be found that for each mesh size, the ITLS B model with 3-point Gauss integration scheme has the smallest error while the ITLS A with 2-point Gauss integration scheme has the biggest error. The use of a 3-point Gauss integration slightly reduces the error for ITLS A compared with 2-point integration. However, with ITLS B the influence of higher order integration is limited. The difference of ITLS B with 2-point and 3-point Gauss integration scheme is minor. It is concluded that 2-point Gauss integration with element subdivision (ITLS B) is optimal.
4.2. Parameter sensitivity study
The influence of the parameters in the dynamic crack growth criterion introduced in Section 3.4 on crack growth is clarified by simulating a dynamic DCB test with a loading velocity of 6.5 m/s. An appropriate parameter set is found by trial and error and used as baseline for this parameter study. The baseline values are Gi= 1.4 N/mm, Ga = 0.67
N/mm, Vi = 40.0 m/s and Va = 14.5 m/s. The bulk material parameters, mesh discretization, interface parameters,
damage function shape parameters γ and damage band length lc are the same as used in the previous section. The
boundary conditions are the same as described in Section 2.2.
4.2.1. Influence of Gi
Five numerical simulations are performed with the initiation toughness Giincreasing from 1.0 to 1.8 N/mm while
keeping the other three input parameters (Ga, Vi, Va) the same as the baseline values. Fig. 11a shows the evolution of
crack length with time for all the simulations along with the experimental test result. The enclosed subplot illustrates how the V(G) relation changes by varying Gi. Fig. 11b depicts the development of the energy release for each
simulation. The time needed for crack initiation to occur is the longest for the case of Gi= 1.8. The accumulation of
energy release rate with time before crack initiation is independent of Gi(see Fig. 11b). Therefore it needs more time
to reach a larger initiation toughness. This means that Gican be calibrated from the initiation time in experimental
measurements.
Crack arrest is observed for all five simulations. The energy release rate tends to decrease during crack growth and when it reaches Ga crack growth is halted. For the two smallest Gicases, more than one crack arrest occurred
while the other three cases only show one crack arrest event. The first crack arrest happens later for larger Gi. Since
the crack speed for the five cases is in the same range [Va, Vi], the time needed for the energy release rate to drop from
0 1 2 3 4 5 6 7 40 60 80 100 0 1 2 3 4 5 6 7 0 1 1.2 1.4 1.6 1.8 2 Gi= 1.0 1.2 1.4 1.6 1.8 Experiment G V Time (ms) Crack length a (mm) Time (ms) Ener gy release rate G (N /mm) Gi= 1.0 1.2 1.4 1.6 1.8 Ga (a) (b)
Fig. 11. (a) Comparison of the time vs. crack length for simulation and experimental test and (b) the evolution of energy release with time for simulation with different Gi.
4.2.2. Influence of Ga
Next, five numerical simulations in which Gais varied from 0.47 to 0.87 N/mm are performed. Results are plotted
in Fig. 12. It is observed from Fig. 12a that for all five cases the crack initiation time is the same. More than one crack arrest event appears for the two cases with highest arrest toughness Gawhile the first three cases only have one
crack arrest. The time at which crack arrest first occurs decreases for increasing Ga. This is a natural outcome of the
dynamic crack growth criterion. Since the crack speed range is unchanged, a smaller arrest toughness takes longer to be reached. Another influence of the arrest toughness is that smaller Ga with all other parameters constant leads to
lower curvature time vs. crack length curve. This observation helps in calibration of the parameters, as a comparison of the average crack speed between simulation and experimental curve can determine the way of adjusting the value for Ga. The case with Ga= 0.67 N/mm gives the best match with the experimental result.
4.2.3. Influence of Vi
The influence of initiation velocity Vi is studied by changing the value from 30 to 50 m/s in five numerical
simulations. Fig. 13 shows the results. Again, the crack initiation time is the same for the five cases. For smaller initiation velocity Vi, the average crack velocity is also smaller as it is shown in the subplot. The shape of the time
vs. crack length curve of the five cases is similar. This means that Vidoes not have a big influence on the evolution
of the energy release rate after crack initiation. The most notable influence is that increasing Vi decreases the time
it takes before the first arrest occurs. This is due to the fact that smaller initiation velocity means smaller averaged crack growth velocity which in turn implies that the energy release rate G decreases slower. The case with Vi = 40
m/s shows the best match with the experimental curve. 4.2.4. Influence of Va
Five simulations with different Va, ranging from 9.5 m/s to 19.5 m/s, are carried out to study the influence of Va.
Results are displayed in Fig. 14. The case with Va = 9.5 m/s shows no crack arrest which can be attributed to the fact
that the averaged crack speed is the smallest. For the other four cases, crack arrest is observed to occur earlier with larger Va. This is because the energy release rate G decreases faster for higher Va. When compared with the influence
of Vishown in previous section, the time at which crack arrest occurs is more sensitive to Va. It is shown that the case
0 1 2 3 4 5 6 7 40 60 80 100 0 1 2 3 4 5 6 7 0 0.47 0.57 0.67 0.77 0.87 1.5 Ga= 0.47 0.57 0.67 0.77 0.87 Experiment G V Time (ms) Crack length a (mm) Time (ms) Ener gy release rate G (N /mm) Ga= 0.47 0.57 0.67 0.77 0.87 Gi (a) (b)
Fig. 12. (a) Comparison of the time vs. crack length for simulation and experimental test and (b) the evolution of energy release with time for simulation with different Ga.
0 1 2 3 4 5 6 7 40 60 80 100 0 1 2 3 4 5 6 7 0 0.5 1 1.5 Vi= 30 35 40 45 50 Experiment G V Time (ms) Crack length a (mm) Time (ms) Ener gy release rate G (N /mm) Vi= 30 35 40 45 50 Gi Ga (a) (b)
Fig. 13. (a) Comparison of the time vs. crack length for simulation and experimental test and (b) the evolution of energy release with time for simulation with different Vi.
0 1 2 3 4 5 6 7 40 60 80 100 0 1 2 3 4 5 6 7 0 0.5 1 1.5 Va= 9.5 12.0 14.5 17.0 19.5 Experiment G V Time (ms) Crack length a (mm) Time (ms) Ener gy release rate G (N /mm) Vi= 9.5 12.0 14.5 17.0 19.5 Gi Ga (a) (b)
Fig. 14. (a) Comparison of the time vs. crack length for simulation and experimental test and (b) the evolution of energy release with time for simulation with different Va.
4.3. Calibration results
Based on the parametric sensitivity study, the ITLS model is independently calibrated on the dynamic DCB tests for each of the test rates of 1.0, 6.5 and 10.0 m/s. For each case, the initiation toughness Gi, the arrest toughness
Ga, the initiation velocity Viand the arrest velocity Va, are determined by comparison with the test curve. All other
simulation inputs are the same as those used for the previous section.
4.3.1. Case I:w˙ = 1 m/s
For the DCB test at ˙w = 1.0 m/s, the calibrated parameters are Gi = 1.7 N/mm, Ga = 0.7 N/mm, Vi = 20 m/s
and Va = 4.0 m/s. The time vs. crack length curve for the experiment and the ITLS simulation are displayed in Fig.
15a. The experimental test result shows three major crack arrests events, which are well reproduced by the numerical simulation. The crack length at arrest is also very close to that in the test data. It is also observed that the crack speed, i.e. the slope of the time vs. crack length curve, has a good match between experimental test and numerical simulation.
The numerical model can also show the evolution of energy release rate during the loading process (see Fig. 15b). Based on data of the crack length at each time from the ITLS model, the evolution of the energy release rate with time can be calculated by invoking the dynamic version of the MBT [11, 51], i.e. when the crack speed ˙a is zero,
G=3 4 E1h3( ˙wT)2 (a+ χh)4 − 33 140 E1h( ˙w)2 η2 (24)
and when the crack speed ˙a is larger than zero,
G=3 4 E1h3( ˙wT)2 (a+ χh)4 − 111 280 E1h( ˙w)2 η2 (25)
where η= pE1/ρ and T is the physical time. A comparison between the energy release rate evaluated in the ITLS
model with Eq. (19) and the MBT as plotted in Fig. 15b reveals that the MBT is still a reasonable approximation for this low loading rate. Hence, the MBT can be a useful tool to gain insight into the energy release rate evolution of the ITLS model. Before crack initiation, the energy release rate keeps increasing until the initiation toughness Gi
0 10 20 30 40 40 60 80 100 0 10 20 30 40 0 0.5 1 1.5 ITLS Experiment ITLS MBT Ga Gi Time (ms) Time (ms) Crack length a (mm) Ener gy release rate G (N /mm) (a) (b)
Fig. 15. (a) Comparison of the time vs. crack length curve for simulation and experimental test at 1.0 m/s and (b) evolution of energy release rate with time for the simulation (solid line) and comparison with the MBT (dash dot line).
initiation. This effect can be understood by Eq. (25), from which two observations can be made. Firstly, an increase of the crack length a by crack growth has an effect of reducing the energy release rate G. Secondly, the external loading prescribed by the loading velocity ˙wcauses G to increase with time. Considering that ˙wfor this case is relatively small, the increase of G caused by external loading can not compensate for the loss of G triggered by a fast growing crack. Therefore, G gradually decreases after crack initiation, although a certain amount of fluctuation is present due to dynamic vibration of the structure. When the arrest toughness Gais reached, the crack is arrested and the crack
speed drops to zero from a speed of Va = 4 m/s instantly. It is also observed in Fig. 15 that the duration time for the
third plateau is the longest while the duration for the first plateau is the shortest. This means that the time needed for the arrested crack to re-initiate becomes longer, which is also in agreement with Eq. (25). The crack length at arrest for a later event is longer, hence, it takes more time to increase the energy release rate from Gato Gito trigger crack
initiation. This feature is visible in experimental measurements as well as in the numerical results. However, it is observed that the reinitiation occurs earlier in the test than in the numerical simulation. This might be due to more complex physics at crack initiation and initial crack tip bluntness.
4.3.2. Case II:w˙ = 6.5 m/s
The numerical results for the test at ˙w= 6.5 m/s have been discussed already in detail in Section 4.2. The results obtained with calibrated parameters Gi = 1.4 N/mm, Ga = 0.67 N/mm, Vi = 40 m/s and Va = 14.5 m/s are shown
in Fig. 16. The experimental curve shows one crack arrest at a crack length of 84 mm, which is very close to what is obtained in the numerical simulation. The crack velocity is well captured by the numerical model as the simulation curve follows the path of the test curve. It should be noted that the calibrated arrest toughness Ga = 0.67 N/mm differs
from the arrest value of around 1.0 N/mm that was extracted from these experimental measurements [10, 11]. The reason for the difference is that in [11] a different definition of the arrest toughness was used for this loading rate, where Gawas determined as a fit through the complete curve, and therefore rather the averaged G during propagation.
The evolution of the energy release rate G during the whole loading process is shown in Fig. 16b. Similar to the previous case, the energy release rate G tends to increase when the crack is stationary and to decrease when the crack
0 1 2 3 4 5 6 7 40 60 80 100 0 1 2 3 4 5 6 7 0 0.5 1 1.5 ITLS Experiment Ga Gi Time (ms) Time (ms) Crack length a (mm) Ener gy release rate G (N /mm) (a) (b)
Fig. 16. (a) Comparison of the time vs. crack length curve for simulation and experimental test at 6.5 m/s and (b) evolution of energy release rate with time for the simulation.
is growing. However, compared to case I the energy release rate curve has more oscillations especially after the crack is arrested. This oscillation is due to an evident effect of system inertia activated by a larger test rate. When the crack is arrested from a finite crack speed, the system inertia has the tendency to resist the sudden change from a propagation crack state to a stationary crack state. Afterwards, the external loading also has to overcome system inertia to cause the stress at the crack tip to build up and to propagate the arrested crack again. For a larger test rate, the interaction between external loading (which promotes crack propagation) and system inertia (which resists the change of system state) leads to more evident oscillating energy release rate evolution.
4.3.3. Case III:w˙ = 10.0 m/s
The calibrated material parameters for the test rate of 10.0 m/s are Gi= 1.3 N/mm, Ga= 0.3 N/mm, Vi= 85 m/s
and Va = 20 m/s. Fig. 17a shows the evolution of crack length with time for experiment and simulation. No clear
crack arrest is observed in the experiment, which means that the crack propagation pattern transitions from unstable to stable crack growth as the test rate increases from 1.0 m/s to 10.0 m/s. The numerical simulation also shows no crack arrest. The crack propagation speed is also reasonably well captured. In absence of arrest events in the test data, calibration does not lead a unique set of values for Gaand Va.
The evolution of energy release rate G with time is shown in Fig. 17b. Before crack initiation, the energy release rate already shows visible oscillations compared to case I and II because higher frequency components are activated for higher loading speed. With the high loading rate, G drops more slowly during crack growth than in the previous cases. This offers a possible explanation for the transition from stick/slip to stable crack growth for increasing loading rate.
4.4. Discussion
Fig. 18 shows the calibrated relation between crack speed V and energy release rate G for the three dynamic DCB tests we have studied. It is shown that the calibrated relation differs significantly for the different test rates. Within
0 0.5 1 1.5 2 2.5 3 40 60 80 100 0 0.5 1 1.5 2 2.5 3 0 0.4 0.8 1.2 ITLS Experiment Ga Gi Time (ms) Time (ms) Crack length a (mm) Ener gy release rate G (N /mm) (a) (b)
Fig. 17. (a) Comparison of the time vs. crack length curve for simulation and experimental test at 10.0 m/s and (b) evolution of energy release rate with time for the simulation.
0 0.5 1 1.5 2 0 20 40 60 80 100 6.5 1 ˙ w= 10
Energy release rate G (N/mm)
Crack
speed
V
(m
/s)
Fig. 18. Calibrated relation of crack speed V and energy release rate G at ˙w= 1.0, 6.5 and 10.0 m/s
the studied test rate range, the initiation toughness Gifor test rate of 1.0, 6.5 and 10.0 m/s is 1.7, 1.4 and 1.3 N/mm
respectively. The arrest toughness Ga for the three cases is 0.7, 0.67 and 0.3 N/mm correspondingly although the
arrest toughness Gaand arrest velocity Vafor 10.0 m/s case are postulated as the crack propagates in a stable manner.
It appears that both Giand Gadecrease with the test rate where Giis more sensitive to the loading rate than Ga. The
initiation velocity Vifor the three cases is 20, 40 and 85 m/s. The arrest velocity Vafor the three cases is 4, 14.5 and
20 m/s. This shows that both Vi and Va increase with the test rate and therefore the average crack velocity is also
increasing. Although the proposed model performs well in reproducing crack arrest, a single set of parameters does not exist to describe the studied dynamic DCB tests at different rates. It can be concluded that in this series of tests, the crack speed is not only a function of the energy release rate at the crack tip. As suggested by Ravi-Chandar [46], strain rate and temperature may influence this relation.
5. Conclusion
In this paper, a rate-dependent cohesive zone model (CZM) and the interfacial thick level set (ITLS) model are used to simulate a series of DCB tests carried out over a test rate range of 3.3 × 10−5to 10.0 m/s on an unidirectional
PEEK/carbon composite laminate reported in [10, 11]. Within the studied test rate range and for this composite material, the following conclusions are drawn.
The rate-independent cohesive zone model fails to capture the inherent loading rate dependency of this composite material and therefore cannot capture the crack propagation process well. Unlike other available rate-dependent cohesive laws which can not capture a decrease in fracture energy for increasing loading rate, the presented rate-dependent cohesive zone model is capable of capturing the average crack growth speed. However, it fails to reproduce crack arrests that do occur in some of the experiments.
Using the interfacial thick level set method, an explicit relation between the crack growth speed V and energy release rate G can be implemented to capture both stable crack propagation and unstable crack propagation in DCB tests. With carefully calibrated parameters for this relation, the ITLS model could well reproduce the crack growth process for the different rates. In case of unstable crack growth, the ITLS approach can capture the number of crack arrest events as well as the timing and duration of these events. This indicates that the adopted dynamic fracture cri-terion with V(G) relation, initiation toughness and arrest toughness is related to the physics of the problem. However, no single V(G) relation could be identified that works for all rates, which means that some physics is still missing in the formulation.
A comparison between the presented rate-dependent CZM and ITLS model shows that: (1) The ITLS method can readily implement a relation of the type V(G), which makes it a good tool for implementing a dynamic crack growth
criterion in fracture mechanics. (2) The ITLS model is more advantageous than the rate-dependent cohesive zone model in terms of the capability of capturing unstable crack growth. (3) Instead of adjusting the parameters for the V(G) relation used in ITLS model, the rate-dependent cohesive zone model requires only a single set of parameters to capture the average response of the DCB at different loading rates. Research is thus needed to establish a more general dynamic crack growth criterion that could give the ITLS approach the ability to be more predictive.
The ITLS is for the first time embedded in a dynamics solution scheme which extends the application of this method from quasi-static and cyclic loading to dynamic loading. It is shown that the V(G) relation, inspired by Ravi-Chandar’s V(K) relation [46], can provide a very good match with experimentally observed arrests and reinitiation phenomena. The good match that has been observed for individual loading rates points out that the V(G) relation is a useful ingredient for describing dynamic crack growth.
Acknowledgement
Financial support from the China Scholarship Council is gratefully acknowledged by YL.
Appendix A. Consistent tangent
Unlike the original implementation of the rate-dependent cohesive law by May [21] for explicit finite element analysis, the cohesive law in this work was implemented in an implicit framework. Therefore the constitutive tangent matrix was derived for this work. The traction vs. displacement jump relation in Eq. (1) can be rewritten in matrix notation as t= K(I − dP)~u (26) in which, I= " 1 0 0 1 # , P= " h~u1i ~u1 0 0 1 # (27)
During damage growth, the consistent tangent Dconsfor a fixed time step size of Tcis defined as
Dcons= ∂t ∂~u = K I − d P − P~u ∂d ∂~u !T (28) with ∂d ∂~u = ∂d ∂∆ ∂∆ ∂~u + ∂d ∂∆0 ∂∆0 ∂~u + ∂d ∂∆f ∂∆f ∂~u (29)
The Eq. (29) is expanded as
∂d ∂∆ = ∆f∆0 ∆2(∆ f−∆0) (30) ∂∆ ∂~u = 1 ∆( h~u1i, ~u2)T (31) ∂d ∂∆0 = − ∆f(∆f−∆) ∆(∆f −∆0)2 (32) ∂∆0 ∂~u = ∂∆0 ∂B ∂B ∂~u+ ∂∆0 ∂~u0 1 ∂~u0 1 ∂σ1 ∂σ1 ∂~u+ ∂∆0 ∂~u0 2 ∂~u0 2 ∂σ2 ∂σ2 ∂~u (33) ∂d ∂∆f = − ∆0(∆ − ∆0) ∆(∆f−∆0)2 (34) ∂∆f ∂~u = ∂∆f ∂B ∂B ∂~u + ∂∆f ∂∆0 ∂∆0 ∂~u+ ∂∆f ∂~u0 1 ∂~u0 1 ∂σ1 ∂σ1 ∂~u+ ∂∆f ∂~u0 2 ∂~u0 2 ∂σ2 ∂σ2 ∂~u+ ∂∆f ∂~uf 1 ∂~uf 1 ∂~u + ∂∆f ∂~uf 2 ∂~uf 2 ∂~u (35)
with ∂∆0 ∂B = ηBη−1 2∆0 (~u02) 2 − (~u01) 2 (36) ∂B ∂~u = ∂B ∂~u1 , ∂B ∂~u2 !T = −2h~u1i(~u2) 2 ∆4 , 2(h~u1i)2~u2 ∆4 !T (37) ∂∆0 ∂~u0 1 =(1 − B∆ η) 0 ~u 0 1, ∂∆0 ∂~u0 2 =∆Bη 0~u 0 2 (38) ∂~u0 1 ∂σ1 = ∂~u0 2 ∂σ2 = 1 K (39) ∂σ1 ∂~u = c1σ01 ˙ ∆1Tc h~u 1i ~u1 , 0 T , ˙∆1≥ ˙∆ref1 0, ∆˙1< ˙∆ref1 , ∂σ2 ∂~u = c2σ02 ˙ ∆2Tc 0, |~u2| ~u2 T , ˙∆2≥ ˙∆ref2 0, ∆˙2< ˙∆ref2 (40) and ∂∆f ∂B = ηBη−1 ∆0 ~u02~u f 2− ~u 0 1~u f 1 (41) ∂∆f ∂∆0 = −∆f ∆0 (42) ∂∆f ∂~u0 1 =1 − B∆ η 0 ~u f 1, ∂∆f ∂~u0 2 =∆Bη 0~u f 2 (43) ∂∆f ∂~uf 1 =(1 − B η)~u0 1 ∆0 , ∂∆f ∂~uf 2 = B η ~u02 ∆0 (44) ∂~uf 1 ∂~u = ∂~uf 1 ∂GIC ∂GIC ∂~u + ∂~uf 1 ∂σ1 ∂σ1 ∂~u (45) ∂~uf 2 ∂~u = ∂~uf 2 ∂GIIC ∂GIIC ∂~u + ∂~uf 2 ∂σ2 ∂σ2 ∂~u (46) with ∂~uf 1 ∂GIC =σ2 1 , ∂~u f 2 ∂GIIC =σ2 2 (47) ∂~uf 1 ∂σ1 = −2GIC (σ1)2 , ∂~u f 2 ∂σ2 = −2GIIC (σ2)2 (48) ∂GIC ∂~u = m1G0IC ˙ ∆1Tc h~u 1i ~u1 , 0 T , ˙∆ref 1 ≤ ˙∆1≤ ˙∆ inf 1 0, otherwise , ∂GIIC ∂~u = m2G0IIC ˙ ∆2Tc 0, |~u2| ~u2 T , ˙∆ref 2 ≤ ˙∆2≤ ˙∆ inf 2 0, otherwise (49) References
[1] M. R. Wisnom, “The role of delamination in failure of fibre-reinforced composites,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 370, no. 1965, pp. 1850–1870, 2012.
[2] M. May, “Measuring the rate-dependent mode I fracture toughness of composites A review,” Composites Part A: Applied Science and Manufacturing, vol. 81, pp. 1–12, 2016.
[3] S. Hashemi, A. Kinloch, and J. Williams, “The Effects of Geometry, Rate and Temperature on the Mode I, Mode II and Mixed-Mode I/II Interlaminar Fracture of Carbon-Fibre/Poly(ether-ether ketone) Composites,” Journal of Composite Materials, vol. 24, no. 9, pp. 918–956, 1990.
[4] H. Zabala, L. Aretxabaleta, G. Castillo, and J. Aurrekoetxea, “Dynamic 4 ENF test for a strain rate dependent mode II interlaminar fracture toughness characterization of unidirectional carbon fibre epoxy composites,” Polymer Testing, vol. 55, pp. 212–218, 2016.
[5] A. Smiley and R. Pipes, “Rate sensitivity of mode II interlaminar fracture toughness in graphite/epoxy and graphite/PEEK composite mate-rials,” Composites Science and Technology, vol. 29, no. 1, pp. 1–15, 1987.
[6] M. Colin de Verdiere, A. Skordos, A. Walton, and M. May, “Influence of loading rate on the delamination response of untufted and tufted carbon epoxy non-crimp fabric composites/Mode II,” Engineering Fracture Mechanics, vol. 96, pp. 1–10, dec 2012.
[7] J.-M. Guimard, O. Allix, N. Pechnik, and P. Th´evenet, “Characterization and modeling of rate effects in the dynamic propagation of mode-II delamination in composite laminates,” International Journal of Fracture, vol. 160, no. 1, pp. 55–71, 2009.
[8] D. Pennas, W. Cantwell, and P. Compston, “The Influence of Strain Rate on the Mode III Interlaminar Fracture of Composite Materials,” Journal of Composite Materials, vol. 41, no. 21, pp. 2595–2614, 2007.
[9] T. Kusaka, M. Hojo, Y.-W. Mai, T. Kurokawa, T. Nojima, and S. Ochiai, “Rate dependence of mode I fracture behaviour in carbon-fibre/epoxy composite laminates,” Composites Science and Technology, vol. 58, no. 3-4, pp. 591–602, 1998.
[10] B. R. K. Blackman, J. P. Dear, A. J. Kinloch, H. Macgillivray, Y. Wang, J. G. Williams, and P. Yayla, “The failure of fibre composites and adhesively bonded fibre composites under high rates of test,” Journal of Materials Science, vol. 30, no. 23, pp. 5885–5900, 1995.
[11] B. R. K. Blackman, A. J. Kinloch, Y. Wang, and J. G. Williams, “The failure of fibre composites and adhesively bonded fibre composites under high rates of test,” Journal of Materials Science, vol. 31, no. 17, pp. 4451–4466, 1996.
[12] H. Zabala, L. Aretxabaleta, G. Castillo, and J. Aurrekoetxea, “Loading rate dependency on mode I interlaminar fracture toughness of unidi-rectional and woven carbon fibre epoxy composites,” Composite Structures, vol. 121, pp. 75–82, 2015.
[13] A. Smiley and R. Pipes, “Rate Effects on Mode I Interlaminar Fracture Toughness in Composite Materials,” Journal of Composite Materials, vol. 21, no. 7, pp. 670–687, 1987.
[14] J. Gillespie, L. Carlsson, and A. Smiley, “Rate-dependent mode I interlaminar crack growth mechanisms in graphite/epoxy and graphite/PEEK,” Composites Science and Technology, vol. 28, no. 1, pp. 1–15, 1987.
[15] A. Aliyu and I. Daniel, “Effects of Strain Rate on Delamination Fracture Toughness of Graphite/Epoxy,” in Delamination and Debonding of Materials, pp. 336–336–13, 100 Barr Harbor Drive, PO Box C700, West Conshohocken, PA 19428-2959: ASTM International, 1987. [16] G. C. Jacob, J. M. Starbuck, J. F. Fellers, S. Simunovic, and R. G. Boeman, “The effect of loading rate on the fracture toughness of fiber
reinforced polymer composites,” Journal of Applied Polymer Science, vol. 96, no. 3, pp. 899–904, 2005.
[17] W. J. Cantwell and M. Blyton, “Influence of Loading Rate on the Interlaminar Fracture Properties of High Performance Composites - A Review,” Applied Mechanics Reviews, vol. 52, no. 6, p. 199, 1999.
[18] F. Zhou, J.-F. Molinari, and T. Shioya, “A rate-dependent cohesive model for simulating dynamic crack propagation in brittle materials,” Engineering Fracture Mechanics, vol. 72, no. 9, pp. 1383–1410, 2005.
[19] A. Corigliano, S. Mariani, and A. Pandolfi, “Numerical analysis of rate-dependent dynamic composite delamination,” Composites Science and Technology, vol. 66, no. 6, pp. 766–775, 2006.
[20] A. Corigliano, S. Mariani, and A. Pandolfi, “Numerical modeling of rate-dependent debonding processes in composites,” Composite Struc-tures, vol. 61, pp. 39–50, jul 2003.
[21] M. May, “Numerical evaluation of cohesive zone models for modeling impact induced delamination in composite materials,” Composite Structures, vol. 133, pp. 16–21, 2015.
[22] M. May, O. Hesebeck, S. Marzi, W. B¨ohme, J. Lienhard, S. Kilchert, M. Brede, and S. Hiermaier, “Rate dependent behavior of crash-optimized adhesives Experimental characterization, model development, and simulation,” Engineering Fracture Mechanics, vol. 133, pp. 112–137, jan 2015.
[23] X. Wei, A. de Vaucorbeil, P. Tran, and H. D. Espinosa, “A new rate-dependent unidirectional composite model - Application to panels subjected to underwater blast,” Journal of the Mechanics and Physics of Solids, vol. 61, no. 6, pp. 1305–1318, 2013.
[24] B. Gozluklu and D. Coker, “Modeling of dynamic crack propagation using rate dependent interface model,” Theoretical and Applied Fracture Mechanics, vol. 85, pp. 191–206, oct 2016.
[25] O. Allix and J.-F. De¨u, “Delayed-damage modelling for fracture prediction of laminated composites under dynamic loading,” Engineering Transactions, 1997.
[26] O. Allix, P. Feissel, and P. Th´evenet, “A delay damage mesomodel of laminates under dynamic loading: basic aspects and identification issues,” Computers& Structures, vol. 81, no. 12, pp. 1177–1191, 2003.
[27] A. Corigliano and M. Ricci, “Rate-dependent interface models: formulation and numerical applications,” International Journal of Solids and Structures, vol. 38, pp. 547–576, jan 2001.
[28] C. M. Landis, T. Pardoen, and J. W. Hutchinson, “Crack velocity dependent toughness in rate dependent materials,” Mechanics of Materials, vol. 32, pp. 663–678, nov 2000.
[29] G. Giambanco and G. Fileccia Scimemi, “Mixed mode failure analysis of bonded joints with rate-dependent interface models,” International Journal for Numerical Methods in Engineering, vol. 67, pp. 1160–1192, aug 2006.
[30] J. Wang, Q. Qin, Y. Kang, X. Li, and Q. Rong, “Viscoelastic adhesive interfacial model and experimental characterization for interfacial parameters,” Mechanics of Materials, vol. 42, pp. 537–547, may 2010.
[31] M. Musto and G. Alfano, “A novel rate-dependent cohesive-zone model combining damage and visco-elasticity,” Computers& Structures, vol. 118, pp. 126–133, mar 2013.
[32] C. Xu, T. Siegmund, and K. Ramani, “Rate-dependent crack growth in adhesives: I. modeling approach,” International Journal of Adhesion and Adhesives, vol. 23, no. 1, pp. 9 – 13, 2003.
[33] N. Mo¨es, C. Stolz, P.-E. Bernard, and N. Chevaugeon, “A level set based model for damage growth: The thick level set approach,” Interna-tional Journal for Numerical Methods in Engineering, vol. 86, no. 3, pp. 358–380, 2011.
[34] P. Bernard, N. Mo¨es, and N. Chevaugeon, “Damage growth modeling using the Thick Level Set (TLS) approach: Efficient discretization for quasi-static loadings,” Computer Methods in Applied Mechanics and Engineering, vol. 233-236, pp. 11–27, 2012.
[35] A. Salzman, N. Mo¨es, and N. Chevaugeon, “On use of the thick level set method in 3D quasi-static crack simulation of quasi-brittle material,” International Journal of Fracture, vol. 202, no. 1, pp. 21–49, 2016.
[36] F. P. van der Meer and L. J. Sluys, “The Thick Level Set method: Sliding deformations and damage initiation,” Computer Methods in Applied Mechanics and Engineering, vol. 285, pp. 64–82, 2015.
[37] K. Moreau, N. Mo¨es, D. Picart, and L. Stainier, “Explicit dynamics with a non-local damage model using the thick level set approach,” International Journal for Numerical Methods in Engineering, vol. 102, no. 3-4, pp. 808–838, 2015.
[38] A. J. Stershic, J. E. Dolbow, and N. Mo¨es, “The Thick Level-Set model for dynamic fragmentation,” Engineering Fracture Mechanics, vol. 172, pp. 39–60, 2017.
[39] M. Latifi, F. P. van der Meer, and L. J. Sluys, “An interface thick level set model for simulating delamination in composites,” International Journal for Numerical Methods in Engineering, vol. 111, no. 4, pp. 303–324, 2017.
[40] M. Latifi, F. P. van der Meer, and L. J. Sluys, “Fatigue modeling in composites with the thick level set interface method,” Composites Part A: Applied Science and Manufacturing, vol. 101, pp. 72–80, 2017.
[41] L. O. Voormeeren, F. P. van der Meer, J. Maljaars, and L. J. Sluys, “A new method for fatigue life prediction based on the Thick Level Set approach,” Engineering Fracture Mechanics, vol. 182, no. Supplement C, pp. 449–466, 2017.
[42] A. Turon, P. Camanho, J. Costa, and C. D´avila, “A damage model for the simulation of delamination in advanced composites under variable-mode loading,” Mechanics of Materials, vol. 38, no. 11, pp. 1072–1089, 2006.
[43] I. M. Daniel and O. Ishai, Engineering mechanics of composite materials. Oxford University Press, 2006.
[44] K. Friedrich, R. Walter, L. A. Carlsson, A. J. Smiley, and J. W. Gillespie, “Mechanisms for rate effects on interlaminar fracture toughness of carbon/epoxy and carbon/peek composites,” Journal of Materials Science, vol. 24, pp. 3387–3398, Sep 1989.
[45] M. A. Meyers, Dynamic behavior of materials. John wiley & sons, 1994.
[46] K. Ravi-Chandar, “Energy Balance and Fracture Criteria,” in Dynamic Fracture (K. Ravi-Chandar, ed.), pp. 71–79, Oxford: Elsevier, 2004. [47] L. B. Freund, Energy concepts in dynamic fracture, pp. 221–295. Cambridge Monographs on Mechanics, Cambridge University Press, 1990. [48] T. Crump, G. Fert´e, A. Jivkov, P. Mummery, and V.-X. Tran, “Dynamic fracture analysis by explicit solid dynamics and implicit crack
propagation,” International Journal of Solids and Structures, vol. 110-111, pp. 113–126, apr 2017.
[49] K. Ravi-Chandar and W. G. Knauss, “An experimental investigation into dynamic fracture: I. Crack initiation and arrest,” International Journal of Fracture, vol. 25, no. 4, pp. 247–262, 1984.
[50] L. Lamberson, V. Eliasson, and A. J. Rosakis, “In Situ Optical Investigations of Hypervelocity Impact Induced Dynamic Fracture,” Experi-mental Mechanics, vol. 52, pp. 161–170, feb 2012.
[51] J. G. Williams, “The fracture mechanics of delamination tests,” The Journal of Strain Analysis for Engineering Design, vol. 24, no. 4, pp. 207–214, 1989.
