Rate-Dependent Cohesive Models for Dynamic Mode I Interfacial Propagation and Failure of Unidirectional Composite Laminates

: With the increasing application of composite materials in anti-impact structure, the development of reliable rate-dependent interlaminar constitutive model becomes necessary. This study aims to assess and evaluate the applicability of three types of rate-dependent cohesive models (logarithmic, exponential and power) in numerical delamination simulation, through comparison with dynamic test results of double cantilever beam (DCB) specimens made from T700/MTM28-1 composite laminate. Crack propagation length history proﬁles are extracted to calibrate the numerical models. Crack propagation contours and fracture toughness data are predicted, extracted and compared to investigate the difference of the three different rate-dependent cohesive models. The variation of cohesive zone length and force proﬁles with the implemented models is also investigated. The results suggest that the crack propagation length can be better predicted by logarithmic and power models. Although crack propagation length proﬁles are well predicted, the numerical calculated dynamic fracture toughness tends to be higher than that of experimental measured results. The three models also show differences in the prediction of maximum loading forces. The results of this work provide useful guidance for the development of more efﬁcient cohesive models and more reliable interface failure simulation of impact problems.


Introduction
Delamination is one of the main damage modes of composite structures subjected to impact loads, which will result in significant stiffness degradation of the structure [1][2][3][4]. The fracture toughness of mode I delamination is generally lower than that of mode II, so the mode I delamination is more likely to occur and has been extensively studied by double cantilever beam (DCB) experiments [5]. For the simulation of delamination of composite materials, the cohesive model is the most widely used method [6]. The concept of cohesive zone was first proposed by Dugdale [7] and Barenblatt [8], in the study of crack-tip yielding and fracture of brittle materials. Initial studies of cohesive model were mainly on the crack propagation in brittle materials [9,10], which was gradually applied to the fracture simulation of other materials and structures, including interface failure of adhesively bonded joints [11] and delamination in composite laminates [12].
The original cohesive model is rate-independent and shows good correlation with quasi-static test results. However, composite structures often undergo dynamic loading during its service life. Additionally, experimental results also elaborated the sensitivity of A few studies were conducted to understand the feasibility of existing cohesive models in delamination studies of composite interfaces [29][30][31]. May et al. [29] identified the necessary of considering simultaneously the rate effect of interface strength and fracture toughness, through the comparison study of four different cohesive models for the impact simulation of composites. Fernandes et al. [30] evaluated the effects of bilinear, linearexponential and trapezoidal cohesive zone model (CZM) laws on the simulation of bonded composite joints, and concluded that the shape of the CZM has no significant effect on the force-displacement curve of the brittle bonding interface, but has a certain effect on the ductile interface. Salih et al. [31] compared a variety of rate-dependent trapezoidal cohesive models and found that the shape of the cohesive model has certain influences on the simulation results of composite interfaces, and proposed a new rate-dependent cohesive model. The modeling results of Salih et al. [31] agree well with the experimental results; it is found that the initial stress has a positive increasing tendency with strain rate until reaching a critical value, following which the fracture energy starts to increase with the increase of strain rate.
Although some researchers have analyzed the impact of different cohesive models, the influence of rate-dependent cohesive models with different formulae based on bilinear cohesive model is not discussed in the previous research. Therefore, the effects of different types of rate-dependent cohesive models on the delamination of composite materials have not been systematically studied and need to be explored further. The influence of different types of rate-dependent cohesive models on the DCB test of carbon fiber composite laminate is analyzed in this paper. Firstly, the materials and experiments are described in Section 2. Then, different rate-dependent cohesive models are introduced in Section 3 and implemented in commercial finite element software Abaqus ® (6.14-1, SIMULIA, Providence, RI, USA) using the VUMAT subroutine. The finite element model is introduced in Section 4. Then, the simulation results of different rate-dependent cohesive models are discussed in Section 5. Conclusions are explored in the last section.

Materials and Tests
The Cytec TM T700/MTM28-1 (Woodland Park, CO, USA) carbon/epoxy unidirectional laminates with fiber volume fraction of 61% are employed as specimens in the DCB experiment. The basic properties of the epoxy resin are ρ = 1.21 g/cm 3 , E = 8.4 GPa, S = 55.4 MPa and ν = 0.33. A non-adhesive polytetrafluoroethylene film (i.e., Teflon) with the thickness of 13 µm is put in the mid-plane of the specimen during the lay-up to introduce the pre-crack. Therefore, the stacking sequence of all DCB specimens is [0 24 /Teflon/0 24 ]. The loading blocks are made by Aluminum and glued to the arms of the DCB laminate using an epoxy adhesive. The temperature in the laboratory is room temperature and the relative humidity is within the normal range.
The simulation in this paper is based on the dynamic double cantilever beam (DCB) experiment in the literature [32]. The dynamic DCB test was carried out by an electromagnetic Hopkinson bar apparatus which can guarantee pure mode I loading, shown in Figure 1. The schematic diagram of the electromagnetic Hopkinson bar is shown in Figure 1a and the testing bench is shown in Figure 1b. As the Hopkinson bar is very long, only the testing bench near the DCB specimen is shown. The Hopkinson tension bar test system and the data acquisition system are contained in this experiment equipment. The Hopkinson tension bar test system contains the incident bars, active coils, inductive coils and electromagnetic control circuit. The active coils are connected with the electromagnetic control circuit. After the electromagnetic induction and interaction between active coil and inductive coil, the tension stress waves can be synchronously generated in the incident bars. Then, the specimen can be loaded symmetrically to get the pure mode I fracture.
There are two pairs of strain gages mounting on each incident bar to record the incident and reflected stress waves. Then, the loading displacement, u 1 and u 2 in Figure 2, can be obtained by one-dimensional stress wave theory according to the incident and reflected strain signals. To ensure accurateness, the crack propagation was monitored by a high-speed There are two pairs of strain gages mounting on each incident bar to record the incident and reflected stress waves. Then, the loading displacement, u1 and u2 in Figure 2, can be obtained by one-dimensional stress wave theory according to the incident and reflected strain signals. To ensure accurateness, the crack propagation was monitored by a highspeed photography (Kirana ® , Pitstone, UK) and two crack-propagation gauges. The maximum frame rate of the camera is 5,000,000 fps, and the frame rate is set to 250,000 fps in this experiment. The results in [32] show that the crack propagation is correctly captured. There are also two strain gages attaching on the two sides of the DCB specimen to prove the symmetrical loading. All the strain signals and the high-speed photography are recorded in data acquisition system.

Theories of Cohesive Models
The rate-independent cohesive model and three rate-dependent cohesive models are introduced in this section. Bilinear cohesive model is chosen as the baseline rate-independent cohesive model because of the simplicity and widely usage in the delamination simulation of composites. The rate-dependent cohesive models are taken from literature [12,21,24] as discussed in introduction section and numerically implemented in Abaqus ® (Dassault Systèmes, Vélizy-Villacoublay, France) as user subroutines (VUMAT).

Bilinear Cohesive Model
For pure mode I or II loading, the bilinear cohesive law is defined by the point of damage initiation δ , σ and complete failure δ , 0 . The separation at damage initiation is defined as: where i means mode I or mode II facture, δ is the separation at damage initiation, is the stress at damage initiation and E is the stiffness of the element.
According to the triangular shape of bilinear cohesive law, the separation at complete failure, δ , can be calculated as: where Gic is the fracture toughness and is defined as the area under the traction-separation curve.

Theories of Cohesive Models
The rate-independent cohesive model and three rate-dependent cohesive models are introduced in this section. Bilinear cohesive model is chosen as the baseline rateindependent cohesive model because of the simplicity and widely usage in the delamination simulation of composites. The rate-dependent cohesive models are taken from literature [12,21,24] as discussed in introduction section and numerically implemented in Abaqus ® (Dassault Systèmes, Vélizy-Villacoublay, France) as user subroutines (VUMAT).

Bilinear Cohesive Model
For pure mode I or II loading, the bilinear cohesive law is defined by the point of damage initiation δ 0 , σ 0 and complete failure (δ c , 0). The separation at damage initiation is defined as: where i means mode I or mode II facture, δ 0 i is the separation at damage initiation, σ 0 i is the stress at damage initiation and E is the stiffness of the element.
According to the triangular shape of bilinear cohesive law, the separation at complete failure, δ c i , can be calculated as: where G ic is the fracture toughness and is defined as the area under the traction-separation curve.
Coatings 2021, 11, 191 5 of 16 In mixed-mode facture, the coupling relationship between mode I and mode II facture is shown in Figure 2a. The quadratic stress failure criterion is adopted as the mixed-mode damage initiation criterion [21]: where σ I and σ II are the current stresses under mode I and mode II loading. The separation under mixed-mode loading, δ m , can be expressed as δ m = δ 2 I + δ 2 II . The mixed-mode initiation separation, δ 0 m , can be obtained by substituting this equation and Equation (1) into Equation (3) as follows [21]: where β is a coefficient representing the degree of mixed-mode and is defined as β = δ II /δ I [21]. The damage propagation under mixed-mode is governed by the mixed-mode fracture toughness, G c . Additionally, G c is given by B-K criterion proposed by Benzeggagh and Kenane [33] as follows: where η is an experience coefficient determined from experiments. Similarly, critical mixed-mode separation, δ c m , can be derived from Equation (5) [21]: In mixed-mode facture, the coupling relationship between mode I and mode II facture is shown in Figure 2a. The quadratic stress failure criterion is adopted as the mixedmode damage initiation criterion [21]: where σI and σII are the current stresses under mode I and mode II loading. The separation under mixed-mode loading, δm, can be expressed as δ = δ δ . The mixed-mode initiation separation, δ , can be obtained by substituting this equation and Equation (1) into Equation (3) as follows [21]: where β is a coefficient representing the degree of mixed-mode and is defined as β = δ δ ⁄ [21].
The damage propagation under mixed-mode is governed by the mixed-mode fracture toughness, Gc. Additionally, Gc is given by B-K criterion proposed by Benzeggagh and Kenane [33] as follows: where η is an experience coefficient determined from experiments. Similarly, critical mixed-mode separation, δ , can be derived from Equation (5) [21]:

Rate-Dependent Cohesive Models
The rate-dependent cohesive model is established based on the rate-independent bilinear cohesive model introduced in last section. There are three types of rate-dependent cohesive models to be discussed in this paper, including logarithmic model, exponential model, and power model. The logarithmic rate-dependent cohesive model is derived from the Johnson-Cook model which is used in characterizing metal failure and is popular used in adhesive joint and composite laminates. The formula of logarithmic model is given as follows [12]: where . δ i is the separation opening rate at the crack tip, defined as .
δ is the rate-dependent initiation stress, σ 0 i,0 is the quasi-static initiation stress and c i is the rate-dependent parameter of stress; G ic . δ is the rate-dependent fracture toughness; G ic,0 is the quasi-static fracture toughness; and m i is the rate-dependent parameter of fracture toughness.
The exponential rate-dependent cohesive model is proposed for composite laminates firstly. The exponential model is shown as follows [21]: where δ 0 i . δ is the rate-dependent separation at damage initiation and δ c i . δ is the ratedependent critical separation. The subscript of "0" represents quasi-static loading conditions and the subscript of "∞" represents limiting loading conditions. The limiting loading condition is an ideal condition where the crack propagation speed is infinite and is impossible to appear. The limiting values can be obtained by amplifying the quasi-static values by the rate-dependent parameter k. The value of the rate-dependent parameter changes between the quasi-static value and the limiting loading value.
The power model is first used in a brittle material, i.e., PMMA plate. The power model is given as follows [24]: Coatings 2021, 11, 191 7 of 16 where n is the rate-dependent parameter. δ c i,0 can be obtained by Equation (2). The initiation stress remains constant in this model. Therefore, the Equation (14) is used in programming, shown as follows: The pure mode I or II traction-separation law of rate-dependent cohesive models are depicted in Figure 2b-d. The difference between these models can be seen in Figure 2. With the increasing of the separation opening rate, the initiation stress and the fracture toughness is both increasing in logarithmic model and exponential model, while only fracture toughness is increasing in power model. The parallelism is observed in the growth of exponential model, while the increase of initiation stress and fracture toughness is relatively independent in logarithmic model.

Model Description
The finite element (FE) model was established in Abaqus ® 6.14-1, shown in Figure 3. The FE model consists of loading block and DCB specimen. The size of the specimen is 80 mm × 15 mm × 6 mm, and the prefabricated crack tip is 30 mm away from the loading side. The mesh size in the width direction is 0.5 mm and the thickness direction is 0.75 mm. To ensure convergence, the mesh size in the longitudinal direction is set as small as 0.2 mm. The zerothickness cohesive interface elements are generated in the middle face of the laminate with the cohesive law presented in Section 3. The element type of laminate is three-dimensional eight-node reduction integral element (C3D8R) with number of 96,000, and the element type of cohesive interface is three-dimensional eight-node cohesive element (COH3D8) with number of 7500. Due to its low impact on the result, the mesh of the loading block is automatically generated and the overall mesh size is 0.5 mm. The computation time for each simulation takes about one hour using two CPUs of i5-3230M @ 2.60 GHz. (Intel, Santa Clara, CA, USA).
To consist with the experimental loading condition, the loading blocks and the external face of DCB specimen are constrained by tie constraint. The four outermost centers of round hole are set as reference points. The reference points and the inside surfaces of the hole are connected by coupling constraints. Only the freedom of rotation in Z direction of the inner surface of the hole is released. The loading displacement, u 1 and u 2 , are loaded on the reference point. There are no additional constraints on the DCB specimen. Coatings 2021, 11, x FOR PEER REVIEW 8 of 16

Material Parameters
There is no damage in the intralaminar of the laminate during DCB test. Therefore, the material property of the laminate is set as orthotropic linear elasticity with engineering constants in Table 1 [32] without damage criterion. The material of loading block is aluminum, so the material property of the loading block is set as isotropic linear elasticity with typical values of E = 70 GPa, ν = 0.32 and ρ = 2700 kg/m 3 .
The parameters need to be obtained can be divided into two categories: rate-independent parameters and rate-dependent parameters. Considering that the crack propagation in DCB test is mode I dominated, the parameters of mode II are assumed to be equal to that of mode I except that the quasi-static fracture toughness of mode I and II were obtained by experiments from literature [32]. The value of E should be as high as possible to make sure that the artificial deformation and artificial strain energy are minimization [34]. Therefore, the values of stiffness of mode I, II and mixed-mode are considered to be rate-independent and both taken as 5 × 10 5 N/mm 3 . Typical values are employed for η [12]. Quasi-static mode I and mode II strength, σ , and σ , are appropriate adjusted based on the typical value of literature [12] to gain a better result.
The reference separation opening rate of logarithmic model adopts the value in the literature [12]. The reference separation opening rate of exponential model is suggested to be 70,000 mm/s in the literature [21], however, this value is too high for the DCB simulation in this paper. After several trial simulations, the reference separation opening rate of exponential model is determined to be 4500 mm/s in this work. The value of reference separation opening rate of power model used in the literature [24] is 5000 mm/s which is still not suitable for this model, so 20 mm/s in the literature [12] is choose for power model. The other rate-dependent parameters, i.e., ci, mi, k and n are calibrated based on the values in [12,21,24] to make the simulation results best match with the curves of crack length vs. time of experiments [32]. The parameters of rate-dependent cohesive interface element are summarized in Table 2.

Material Parameters
There is no damage in the intralaminar of the laminate during DCB test. Therefore, the material property of the laminate is set as orthotropic linear elasticity with engineering constants in Table 1 [32] without damage criterion. The material of loading block is aluminum, so the material property of the loading block is set as isotropic linear elasticity with typical values of E = 70 GPa, ν = 0.32 and ρ = 2700 kg/m 3 .
The parameters need to be obtained can be divided into two categories: rate-independent parameters and rate-dependent parameters. Considering that the crack propagation in DCB test is mode I dominated, the parameters of mode II are assumed to be equal to that of mode I except that the quasi-static fracture toughness of mode I and II were obtained by experiments from literature [32]. The value of E should be as high as possible to make sure that the artificial deformation and artificial strain energy are minimization [34]. Therefore, the values of stiffness of mode I, II and mixed-mode are considered to be rate-independent and both taken as 5 × 10 5 N/mm 3 . Typical values are employed for η [12]. Quasi-static mode I and mode II strength, σ 0 I,0 and σ 0 II,0 are appropriate adjusted based on the typical value of literature [12] to gain a better result.
The reference separation opening rate of logarithmic model adopts the value in the literature [12]. The reference separation opening rate of exponential model is suggested to be 70,000 mm/s in the literature [21], however, this value is too high for the DCB simulation in this paper. After several trial simulations, the reference separation opening rate of exponential model is determined to be 4500 mm/s in this work. The value of reference separation opening rate of power model used in the literature [24] is 5000 mm/s which is still not suitable for this model, so 20 mm/s in the literature [12] is choose for power model. The other rate-dependent parameters, i.e., c i , m i , k and n are calibrated based on the values in [12,21,24] to make the simulation results best match with the curves of crack length vs. time of experiments [32]. The parameters of rate-dependent cohesive interface element are summarized in Table 2.

Model Correlation and Prediction of Crack Propagation Length
The crack propagation length-time curves simulated by different rate-dependent cohesive models are shown in Figure 4. The average speeds of different curves are calculated by linear fitting and the error between experiment and different models are shown in Table 3. Overall, the crack propagation curves of all three models are close to the experiment because the parameters are calibrated by these curves. The crack propagation length curves of logarithmic model and power model are very similar and close to the experimental curves, while the curves of exponential model are relatively far from the experimental curves. However, the crack initial time of exponential model is close to that of experiment while the crack initial time of logarithmic model and power model are always later than that of experiment.
The crack speed of logarithmic model is higher than that of experiment initially, and becomes lower gradually in 108 m/s condition, while the average speed of the logarithmic model is close to that of experiment with error of 2.78%. The same trend of crack speed of logarithmic model can be seen in 175 and 253 m/s conditions, though the trend is less and less obvious and the curve of 253 m/s condition is nearly linear. The average speed errors of logarithmic model are relatively larger in 175 and 253 m/s conditions, with errors of 18.86% and 9.09%, respectively, than that shown by the 108 m/s condition. The crack propagation length curves of exponential model are all nearly linear in three conditions. However, the average speed errors of exponential model are higher with 19.44% and −16.60% in 108 and 253 m/s conditions, respectively, than that of logarithmic model. Although the average speed error of exponential model in 175 m/s condition is very small of −3.43%, the initial crack speed is smaller than that of experiment which result the curve away from experiment. The behavior of crack propagation length curve of power model is similar with that of logarithmic model and the average speed errors is 8.33%, 18.29% and 7.51% in 108, 175 and 253 m/s conditions, respectively.
The propagation contours of cohesive element during the dynamic tests are simulated and compared in right side of Figure 4. The total cracking time decreases with the increase of crack-propagation velocity. Due to the approximate plane stress state of the side and the approximate plane strain state of the middle of the specimen, the crack tip is concave in x-z plane and the crack is always initiated firstly at the middle section [35]. To be consistent and comparable with the experiment, the numerical predicted crack-propagation length is counted by erosion of cohesive element along the edges. The predicted cohesive-zone length (length of red-color zone) of power model is obviously larger than that of logarithmic and exponential models which will be further discussed in Section 5.3.

Prediction of Crack Propagation Toughness
The crack propagation is governed by the dynamic fracture toughness. For the crack propagation velocity much smaller than the Rayleigh wave speed, the energy release rate can be used to express the fracture toughness. The energy history of the whole model can be obtained in numerical simulation and the energy release rate can be calculated by the following equation based on the energy conservation:

Prediction of Crack Propagation Toughness
The crack propagation is governed by the dynamic fracture toughness. For the crack propagation velocity much smaller than the Rayleigh wave speed, the energy release rate can be used to express the fracture toughness. The energy history of the whole model can be obtained in numerical simulation and the energy release rate can be calculated by the following equation based on the energy conservation: where the U E , U S and U K are the external work, total strain energy and kinetic energy, respectively. The B is the thickness of the specimen and the a is the crack propagation length. The energy release rate-time curves of three rate-dependent cohesive models are depicted in Figure 5. Overall, since the crack propagation curves of logarithmic model and power model are similar and crack propagation length is governed by fracture toughness, the energy release rate curves of these two models are also very similar. Note that the initial stress of power model remains constant, which means that the change of initial stress is not very important in this loading condition. With the average crack propagation speed increasing, the general energy release rates of the three models increase in which the exponential model is the most obvious. Valleys can be seen in the final stage of all the curves.
The energy release rate curves of logarithmic model and power model keep oscillation during the crack propagation. The energy release rate of logarithmic model is a little higher than that of power model after 0.45 ms in the 108 m/s condition, which means that the crack propagation speed of logarithmic model is lower than that of power model as can be seen in Figure 4a. The first point of the power model in 175 and 253 m/s condition is extremely high because the first point of da is very small, which should be neglected.
where the UE, US and UK are the external work, total strain energy and kinetic energy, respectively. The B is the thickness of the specimen and the a is the crack propagation length. The energy release rate-time curves of three rate-dependent cohesive models are depicted in Figure 5. Overall, since the crack propagation curves of logarithmic model and power model are similar and crack propagation length is governed by fracture toughness, the energy release rate curves of these two models are also very similar. Note that the initial stress of power model remains constant, which means that the change of initial stress is not very important in this loading condition. With the average crack propagation speed increasing, the general energy release rates of the three models increase in which the exponential model is the most obvious. Valleys can be seen in the final stage of all the curves.
The energy release rate curves of logarithmic model and power model keep oscillation during the crack propagation. The energy release rate of logarithmic model is a little higher than that of power model after 0.45 ms in the 108 m/s condition, which means that the crack propagation speed of logarithmic model is lower than that of power model as can be seen in Figure 4a. The first point of the power model in 175 and 253 m/s condition is extremely high because the first point of da is very small, which should be neglected. The curves of exponential model increase initially then decrease gradually in all conditions, while the variation of logarithmic and power models is little, especially in 253 m/s condition. The overall energy release rates of exponential model are lower than that of logarithmic and power model in 108 m/s condition which results that the average crack speed is the highest in three models. Additionally, the overall energy release rates of exponential model are higher than that of logarithmic and power model in 253 m/s condition which results that the average crack speed is the lowest in three models. The curves of exponential model increase initially then decrease gradually in all conditions, while the variation of logarithmic and power models is little, especially in 253 m/s condition. The overall energy release rates of exponential model are lower than that of logarithmic and power model in 108 m/s condition which results that the average crack speed is the highest in three models. Additionally, the overall energy release rates of exponential model are higher than that of logarithmic and power model in 253 m/s condition which results that the average crack speed is the lowest in three models.
The dynamic fracture toughness varies from 0.1 to 1.1 N/mm in experiment and may be totally different between two tests with the same average crack velocity. According to the experiment, the fracture toughness is governed by the crack tip opening rate and the crack propagation speed. However, within the phenomenology cohesive model discussed in this paper, the dynamic fracture toughness is only the function of crack tip opening rate. Therefore, the energy release rates in numerical results is 2-3 times than the highest energy release rate in experiment when the crack propagation length curves of numerical results are close to that of experiment. The higher fracture toughness is a compromise to compensate some physical mechanisms which is the disadvantage of phenomenology model. Although this error is relatively large in this DCB simulation, the resulting error of absorbed energy can be neglected in impact simulation of big structures. Therefore, the phenomenology models are still popular used in many studies.

Prediction of Cohesive-Zone Length
The stability of finite element simulation of cohesive zone is strongly dependent on cohesive zone length. The cohesive zone length is the length of the failure process zone ahead of the crack tip where the cohesive forces are acting. Within the finite element method framework, the cohesive zone length can be defined as: (16) where N e is the number of damaged cohesive elements and d e is the size of the cohesive elements which is 0.2 mm in this simulation. The cohesive zone lengths of three models in all conditions are counted and shown in Figure 6 and the N e is expressed as units of the Y-axis. According to [11], the cohesive zone length in an infinite body can be expressed by the function of strength and fracture toughness: where r = 2 for the infinite body. While in other situations, r may be other values, e.g., r = 0.8~0.9 in the literature [11].
The phenomenon depicted in Figure 6 can be explained according to Equation (17). The most appeared cohesive zone lengths of logarithmic model, exponential model and power model are 4, 6 and 8, respectively. For the initial stress remaining unchanged in power model, the initial stress is the smallest in three models which results that the cohesive zone length is the largest. The cohesive zone length of logarithmic model is smaller than that of exponential model which means that the initial stress of logarithmic model is larger than that of exponential model when the fracture toughness is the same.
For power model, the fracture toughness becomes larger with the increase of crack speed which results that the cohesive zone length becomes larger. For exponential model, the fracture toughness changes as the square of stress. Therefore, if the r = 2 in this simulation, the cohesive zone length of exponential model should remain constant. However, the trend of cohesive zone length of exponential model is similar with that of power model. Therefore, the r is smaller than 2 in this simulation. For logarithmic model, the trend is contrary to that of power model which means that the change of initial stress is larger than that of fracture toughness with the change of crack speed. Coatings 2021, 11, x FOR PEER REVIEW 13 of 16

Sensitivity of Loading Force
The load forces are calculated by add the reaction forces of the two reference points of top loading block. The load force-time curves of three models are shown in Figure 7. Totally speaking, the forces of the elastic loading section between different models and conditions are the same which is easy to understand. After the delamination beginning, the trend of propagation force curves is gradually decrease because of the energy release which is similar to that of the quasi-static condition. After the laminates totally cracking, the forces become very oscillating and meaningless. Before the laminates totally cracking, the forces in 108 and 175 m/s conditions is normal and higher than 0 N. While in the 253 m/s condition, the forces are lower than 0 N at 0.35 ms in all models, which may be due to the higher loading rate.
The abscissa of vertical lines in Figure 7 means the initial time of the crack. In the typical loading force-time curves of quasi-static, the delamination begins at the peak force. While in dynamic loading condition, the time of delamination initiation is earlier than that of the peak force except the logarithmic and power models in 108 m/s condition. This may be due to the inertial effect under dynamic loading condition. Additionally, in dynamic loading conditions, the peak forces are almost 10 times than that in quasi-static loading conditions which is also a typical phenomenon under dynamic loading.

Sensitivity of Loading Force
The load forces are calculated by add the reaction forces of the two reference points of top loading block. The load force-time curves of three models are shown in Figure 7. Totally speaking, the forces of the elastic loading section between different models and conditions are the same which is easy to understand. After the delamination beginning, the trend of propagation force curves is gradually decrease because of the energy release which is similar to that of the quasi-static condition. After the laminates totally cracking, the forces become very oscillating and meaningless. Before the laminates totally cracking, the forces in 108 and 175 m/s conditions is normal and higher than 0 N. While in the 253 m/s condition, the forces are lower than 0 N at 0.35 ms in all models, which may be due to the higher loading rate.
The abscissa of vertical lines in Figure 7 means the initial time of the crack. In the typical loading force-time curves of quasi-static, the delamination begins at the peak force. While in dynamic loading condition, the time of delamination initiation is earlier than that of the peak force except the logarithmic and power models in 108 m/s condition. This may be due to the inertial effect under dynamic loading condition. Additionally, in dynamic loading conditions, the peak forces are almost 10 times than that in quasi-static loading conditions which is also a typical phenomenon under dynamic loading. Similar to the crack length and fracture toughness, the forces of logarithmic and power models are very close, which again shows that the interface strength parameter does not affect the results much in this dynamic DCB simulation. Interestingly, the peak forces of exponential model in all conditions are lower than that of logarithmic and power models, while the propagation forces of exponential model in 175 and 253 m/s conditions is higher than that of logarithmic and power models.

Conclusions
This paper studies and examines the applicability of different rate-dependent cohesive models in the simulating of model I delamination of unidirectional composite. The capabilities of logarithmic, exponential and power-type rate-dependent cohesive models in predicting crack propagation length, energy release rate, cohesive zone length and reaction force are systematically discussed and compared. The main conclusions of this study are listed as follows:

•
All three models predict well the crack length profiles with average velocity errors within 20%. The exponential model shows relatively larger error compared with that of the logarithmic and power models.

•
The predicted crack propagation toughness values of all three cohesive models are much higher than that of experimental results. This indicates that physical mechanism is sacrificed by the phenomenological model in purpose of predicting accurately the cracking behavior.

•
In the presented numerical simulation results, the correlation of cohesive-zone length with crack length of DCB specimen are consistent with the results and theory predictions made in the literature [11].

•
The peak load of the exponential model is always lower than that of the logarithmic and power models, while in propagation stage, the load force of exponential model Similar to the crack length and fracture toughness, the forces of logarithmic and power models are very close, which again shows that the interface strength parameter does not affect the results much in this dynamic DCB simulation. Interestingly, the peak forces of exponential model in all conditions are lower than that of logarithmic and power models, while the propagation forces of exponential model in 175 and 253 m/s conditions is higher than that of logarithmic and power models.

Conclusions
This paper studies and examines the applicability of different rate-dependent cohesive models in the simulating of model I delamination of unidirectional composite. The capabilities of logarithmic, exponential and power-type rate-dependent cohesive models in predicting crack propagation length, energy release rate, cohesive zone length and reaction force are systematically discussed and compared. The main conclusions of this study are listed as follows:

•
All three models predict well the crack length profiles with average velocity errors within 20%. The exponential model shows relatively larger error compared with that of the logarithmic and power models.

•
The predicted crack propagation toughness values of all three cohesive models are much higher than that of experimental results. This indicates that physical mechanism is sacrificed by the phenomenological model in purpose of predicting accurately the cracking behavior.

•
In the presented numerical simulation results, the correlation of cohesive-zone length with crack length of DCB specimen are consistent with the results and theory predictions made in the literature [11].

•
The peak load of the exponential model is always lower than that of the logarithmic and power models, while in propagation stage, the load force of exponential model is higher than that of the logarithmic and power models especially under high loading rate.