Modeling of dynamic mode I crack growth in glass fiber-reinforced polymer composites fracture energy and failure mechanism

The mode-I dynamic fracture energy and failure mechanisms of glass fiber-reinforced polymer composites are investigated with an embedded cell model of the single-edge-notched-tension (SENT) geometry. Under an applied dynamic loading, a crack may propagate in the embed- ded microstructure, accompanied by the development of a fracture process zone in which fiber/matrix debonding, matrix cracking and ductile matrix tearing are observed. Reaching a maximum nominal strain rate of 250/s, a series of SENT tests are performed for different loading velocities and specimen sizes while the dynamic energy release rate is evaluated using the dynamic version of the 𝐽 -integral. The influence and interaction of loading rate, time-dependent material nonlinearity, structural inertia and matrix ligament bridging on the fracture toughness and failure mechanisms of composites are evaluated. It is found that with the given material parameters and studied loading rate range, the failure type is brittle with many microcracks but limited plasticity in the fracture process zone and a trend of increasing brittleness for larger strain rates is observed. The inertia effect is evident for larger strain rates but it is not dominating. An 𝑅 -curve in the average sense is found to be strain-rate independent before the fracture process zone is fully developed and afterwards a velocity–toughness mechanism is dictating the crack growth.


Introduction
Fiber reinforced polymer composites have been used in impact-resistant devices, automotives, aircraft structures due to their potential for high strength-to-weight ratios and impact energy absorption. To be able to fully exploit the potential of impact behavior of composites it is necessary to investigate dynamic crack propagation, in particular the underlying mechanisms, microstructural effects and the fracture energy.
Starting from Griffith's ideas postulated for equilibrium cracks [1] and its extension by Mott for dynamic fracture [2], dynamic fracture can be investigated on an energetic basis. The dynamic energy release rate is the energy released into the crack tip process zone per unit crack extension and must be equal to the energy required per unit extension [3]. Generally, both and are functions of crack propagation history, in particular, the crack speed . Freund [4] showed that for mode-I crack propagation in homogeneous materials under elastodynamic conditions and in plane strain state the dynamic energy release rate can be expressed in the following form: where is the Young's modulus, is the Poisson's ratio, is the mode-I dynamic stress intensity factor and is a universal function of crack speed . The dynamic stress intensity factor tends to zero as approaches the Rayleigh wave speed , which implies a limiting crack speed of in mode-I. Corresponding  Single edge notched tension levels of propagation velocity, the crack surface roughness is observed to show different features since material in the fracture process zone might experience high strain-rate plasticity, microcracks nucleation, thermomechanical interaction and other complex deformation/failure mechanisms. Upon an increase of crack speed, the crack surface first appears to be almost flat (mirror regime), next a rougher surface with conic marks forms (mist regime), and finally (micro)branching takes place (hackle regime). Phenomenologically, a relationship between fracture energy and crack speed therefore exists. The relation between the dynamic fracture energy and the crack speed for composites is determined by the rate-dependent deformation and failure process occurring across multiple length scales and time scales. More specifically, the contributing mechanisms can be roughly classified as viscous material behavior, changes in fracture mechanism, inertia effects and thermomechanical effects. Firstly, there is the role of viscosity of composite constituents (polymer, fiber and interfaces) and its interaction [5,6]. Shirinbayan et al. [5] postulated that a specific characteristic time for a local damage to occur might exist and this time scale is related to the viscoelastic behavior of the matrix or fiber/matrix interface. Fitoussi et al. [6] argued that for high rate a local strained zone around a debonded interface dissipates the strain energy and accordingly hinders the interfacial crack propagation through the matrix, which causes a delay of the damage at macroscopic scale. Secondly, there can be rate-dependency of the fracture mechanism induced by different failure processes (e.g. fiber failure with fiber pull-out, matrix damage, fiber-matrix interface failure) occurring at microscale level under different loading rates. For instance, for quasi-static tests delamination is often dominated by fiber/matrix interface failure while resin rich brittle fracture zones have been found more dominant in dynamic tests [7][8][9][10][11]. The extent of plastic deformation may decrease with increased loading rate, which represents a ductile-to-brittle transition in the process zone. Thirdly, there are inertia effects characterized as inertia resistance for rapid deformation, damage formation and crack propagation [12,13]. Due to material heterogeneity, micro-inertia effects also arise as a result of multiple wave reflection and transmission occurring at the interfaces between the constituents, which can result in complex spatiotemporal scenarios of damage and failure evolution, initiated at multiple spots [14,15]. Finally, there can be thermomechanical dissipation as a transition from isothermal to adiabatic deformation and failure process for composites is expected for increasing loading rate [16][17][18].
Computational models have been developed on the mesoscale to capture deformation and failure in composites. For such models, the composite ply is considered as a homogenized material where damage and failure can be described by continuum damage models [19,20] or extended FEM models [21] with failure-mode-based criteria and different stiffness degradation laws for the different failure processes. However, such models inevitably lead to complex constitutive and damage laws that require extensive experimental calibration and the observations obtained at these scales do not provide enough detail about the mechanical processes that explain the inelastic behavior of the material. Hence, computational micromechanical models are an appealing option for investigating the dynamic fracture energy and the interplay of different mechanisms of dynamic crack growth in composites. Microscale-based approaches can be roughly classified as: the representative volume element (RVE) based multiscale approach, the (modified) boundary layer approach and the embedded cell approach. An RVE is a characteristic sample of heterogeneous material that should be large enough to contain sufficient composite micro-heterogeneities in order to be representative, however it should also be much smaller than the macroscopic structure size [22]. The RVE-based multiscale approach assumes multiple spatial and (or) temporal scales. Solution of finer-scale problems is analyzed in an RVE and information of the finer-scale is hierarchically passed into a coarser scale by bridging laws. For a two-scale scheme, at the macroscopic level the strain localization can be represented by cohesive cracks with strong discontinuity kinematics and a proper kinematical information transfer from the macro-to-microscales [23][24][25][26]. However, the implementation of this method is not readily available in a general-purpose finite element code and the computational cost of this method can be prohibitively high. The (modified) boundary layer formulation considers a small layer of material near the crack tip with well-defined singularity displacement fields applied at the edges of the layer. Numerical solution of this problem allows a quantification of the energy dissipation under such singularity field with energy integrals. This approach has been applied to study elastic-plastic ductile cracking in a homogeneous material [27,28] and the effective fracture toughness of a heterogeneous material [29][30][31]. However, it is not clear how to apply the boundary conditions if a singularity field cannot clearly be defined. For the embedded cell approach, full details of the heterogeneous composite microstructure (including the random spatial distribution of the fibers) are explicitly resolved in the fracture region with a finer discretization. Meanwhile, the rest of the model is considered to be a homogeneous medium with simple constitutive equations (obtained a priori with any appropriate homogenization technique) and coarser discretization. The region in which the microstructure is resolved should be small so that the computational cost is affordable. However, it should be sufficiently large to include all the area in which damage occurs during the propagation of the crack, thus energy spent by the different failure micromechanisms (interface debonding, matrix cracking, matrix plastic deformation, etc.) is properly taken into account. This approach has been used in analysis of quasi-static crack propagation of in composite material and to compute the fracture toughness associated to different failure modes [32,33].
In this paper, a multiscale numerical model using the embedded cell approach is developed to evaluate the mode-I fracture energy of dynamic crack propagation in fiber-reinforced composites and to investigate the associated failure mechanisms. Specifically, the single edge notched tension (SENT) specimen is analyzed. The paper is organized as follows: in Section 2, details of the embedded cell model of the SENT specimen are given. Section 3 presents the typical deformation and failure phenomena in a series of tests on SENT specimen and the obtained relations between the dynamic fracture energy and crack speed . The failure mechanisms in the fracture process zone of the composites are discussed in Section 4.

Numerical model
To compute the mode-I fracture energy of dynamic crack growth in fiber-reinforced polymeric composites, a embedded cell model of an SENT specimen with a width of and a length of is developed. The SENT specimen is favored over other Mode-I tests, such as the double cantilever beam test, because the absence of bending deformations (with both tension and compression) in the SENT is beneficial for numerical robustness under dynamic loading condition. As it is shown in Fig. 1, an initial notch of length 0 along -direction is created on one edge of the specimen and a symmetric displacement loading is applied on the top and bottom edge of the specimen with a prescribed velocity oḟ. In the vicinity of the initial notch tip, a composite microstructure of rows and columns of repeating RVEs is embedded in a homogenized medium of the composites. The RVE has a stochastic distribution of 5 × 5 fibers with a fiber diameter = 5 and a fiber volume fraction of 60%. It is generated by a discrete element method generator called HADES, following the procedures in Liu et al. [34].
The matrix material of the microstructure is assumed to be epoxy modeled as viscoelastic-viscoplastic model as detailed in Section 2.1 while the fiber is assumed to be linear elastic. The material around the embedded microstructure is treated as a homogeneous isotropic elastic solid whose behavior is obtained by a standard computational homogenization scheme (see Appendix A in [34]) from elastic constants of the fibers and matrix in an RVE. Cracking is allowed to develop only inside the matrix and on the fiber/matrix interfaces in the embedded cell. Following Camacho and Ortiz [35], a dynamic insertion technique of cohesive elements, introduced in Section 2.2, is used to capture cracking. The whole numerical model is solved with an implicit dynamics scheme. A plane strain condition is assumed and the two-dimensional plane is considered as the transverse plane of a fiber-reinforced composite ply. The algorithm is described in detail in Section 2.3. The dynamic energy release rate for the composites is computed by utilizing the dynamic version of the -integral with its formulation shown in Section 2.4.

Polymer model
The polymer matrix of the embedded microstructure is assumed to have a constitutive behavior according to a viscoelasticviscoplastic (VE-VP) model following Rocha et al. [36]. Following the conceptual representation of the VE-VP model in Fig. 2, two contributing constitutive models are involved: a linear viscoelastic model and a viscoplastic component represented by a Perzyna-type overstress formulation with a rate-independent backbone of a pressure-dependent plasticity model.

Viscoelasticity
By assuming a linear viscoelastic model the stress for time is expressed with Boltzmann's hereditary integral of the elastic strain : in which ( ) is a time-dependent stiffness which can be expressed with a time-dependent shear stiffness ( ) and bulk stiffness ( ): where ( ) and ( ) can be further expanded as a combination of a Prony series of shear elements and bulk elements and a long-term contribution: in which ∞ and ∞ represent the long-term shear and bulk stiffness, and , , and are shear and bulk stiffness and relaxation time of the Maxwell elements, respectively. The fourth-order deviatoric and volumetric operator tensors introduced in Eq. (2) are defined as: where is the Kronecker delta. These operator tensors can also be used to decompose the elastic strain into a deviatoric part , and a hydrostatic part , : By substituting Eqs. (2) and (3) into Eq. (1), the stress can be expressed as: in which the deviatoric viscoelastic stress contribution: and the hydrostatic viscoelastic stress contribution:

Viscoplasticity
The viscoplasticity model is assumed to be a Perzyna-type model with a backbone of a pressure-dependent hardening plasticity model. The yield function of the plasticity model is defined as: ( , ) = 6 2 + 2 1 ( − ) − 2 (9) in which 1 = is the first stress invariant, 2 = 1 2 is the second invariant of the deviatoric stress , and and are the yield stress in tension and compression, respectively. The yield stress, or , is a function of the accumulated equivalent plastic strain . In an incremental form, the accumulated equivalent plastic strain is defined as: in which is the plastic Poisson's ratio. The desired contraction behavior is implemented through a non-associative flow rule which is expressed in an incremental form as: where is the incremental plastic multiplier and the parameter is: By allowing the overstress to develop beyond the yield surface, a viscous time scale is introduced in the model. A Perzyna-type of overstress formulation is adopted, which gives the evolution of the plastic multiplier : in which 0 and 0 are the yield stress values when = 0, and are viscoplastic coefficients and is the time increment.

Energy dissipation
The free energy of the VE-VP model can be expressed as: 7 Y. Liu et al. in which ℎ is the plastic hardening energy. According to the second law of thermodynamics, the Clausius-Duhem inequality for the isothermal case is imposed: Following the derivation in Rocha et al. [36], the work of energy dissipation per unit volume for viscoelasticity and viscoplasticity can be expressed as: Summing up Eqs. (16) and (17) and integration over time give the accumulated dissipation per unit volume as: To derive Eq. (18) the terṁh is neglected because the plastic hardening in the polymer also contributes to the energy for growing a macroscopic crack. In the numerical model, the total dissipated energy of the polymer matrix can be computed as the volume integral over the embedded microstructure: in which is the volume of the embedded microstructure zone.

Cohesive crack with Ortiz model
The microcracks in the embedded zone, representing fiber/matrix debonding and matrix cracking, are modeled with the cohesive zone model. Instead of inserting cohesive elements between element boundaries before the simulation starts, in this study the cohesive elements are placed on the fly following the shifted cohesive law technique described in Camacho and Ortiz [35]. A stress-based failure criterion is introduced to determine when and where the cohesive element should be inserted. The crack always starts at the middle node of edges of six-node triangle elements by splitting the nodes (see Fig. 3). Because cohesive elements are inserted on the fly, continuity of the response requires that the adopted cohesive law is an initially rigid linear softening law. As a consequence there is no initial stiffness present, of which the value could otherwise affect the overall compliance of the material or the stress development under dynamic loading conditions.

Cohesive element insertion criterion
Considering mixed-mode fracture, the adopted stress-based failure criterion reads [35]: where is the cohesive strength and the effective stress ef f is defined as: in which = ( , ) is the traction of cohesive surface along the normal direction and shear direction in the local {n,s} frame, is the displacement jump along normal and shear direction, is a shear stress factor, and is the friction coefficient.

Shifted cohesive law
To construct an initially rigid law without singularity, a shifted cohesive law is adopted [37]. As seen in Fig. 4, starting from a traction separation relation with a finite initial stiffness, a shift of this relation is applied such that the traction for zero crack opening is equal to the traction at crack initiation. This leads to the desired initially rigid behavior.
For the shifted cohesive law, the traction is computed not from the actual displacement jump, but from a translated displacement jump [[ ]]: The shift [[ ]] 0 is computed from the bulk stress at the moment of crack initiation and can be expressed as: in which 0 is the stress at crack initiation and is a dummy stiffness. The traction is updated in the local {n, s} frame as: with the effective traction defined as ef f = [[ ]] and a damage tensor is defined as: is the normal component of the effective traction ef f , is the Kronecker delta and the Macaulay bracket is define as ⟨ ⟩ = 1 2 ( + | |). The damage evolution is introduced according to a bilinear relation where the equivalent displacement jump is: and the equivalent displacement representing onset of failure 0 reads: and the equivalent displacement representing complete damage is: where is the fracture energy. The equivalent traction corresponding to onset of damage 0 is introduced in Eqs. (28) and (29) with its definition as the norm of the traction at damage initiation 0 , Note that this cohesive law formulation is generic in that it can be combined with any failure criterion, but its particular behavior is linked to the used failure criterion through Eqs. (28) and (30).

Energy dissipation in fiber/matrix interfaces
The total dissipated energy of the cohesive interfaces reads: in which is the area of the cohesive elements and the incremental dissipation density can be computed by: with the incremental of a variable

Solution scheme
Implicit dynamics analysis is carried out. The adopted semi-discretization scheme includes an implicit time integration of the Newmark-type and a spatial discretization with six-node triangular elements. The solution program flow is illustrated in Box 1.
There are a few items to be noted: (1) in step 5 of the algorithm, the dynamic system of equations is solved with a Newton-Raphson scheme. In certain circumstances, convergence cannot obtained by a large time step size. Then, an adaptive stepping algorithm is used such that the time step size is reduced and the system equation is solved with smaller time steps until convergence is reached. (2) when new cohesive elements are inserted, the mesh is updated and the same step is solved again to ensure that the final converged solution for the time step does not violate the failure criterion.

J-integral calculation
Following Anderson [38], for a fast moving crack the amount of energy flowing into the crack tip region through the contour can be calculated by the crack tip energy flux integral (see Fig. 5): and is a line segment of path , = ∫ 0̇i s the stress work density, = 1 2̇̇i s kinetic energy density, is the outward unit normal to the contour , is the stress, is the displacement. This -integral formulation is valid for time-dependent as well as history-dependent material behavior because it was derived from a generalized energy balance. In the special case of a constant crack propagation speed and steady-state crack propagation in homogeneous hyperelastic material the dynamic -integral becomes path-independent [39]. In this study, the integral contour is defined outside of the embedded microstructure similar to what was done in the embedded cell model by Herráez et al. [30].
To facilitate the application of the dynamic -integral into a FEM framework, an equivalent domain integral is introduced to replace the line integral introduced [39]. Fig. 5 shows an example of the selected path along boundaries of one ring of finite elements alongside with an extra remote path , one segment of the initial surface + and one segment of the initial surface − . A closed path = + + + − − is therefore constructed in counter-clockwise direction. In addition, a weighting function ( ), which must be continuous and differentiable and fulfills the requirements, Y. Liu et al.
Substitution of Eq. (35) into Eq. (37) gives the final expression for dynamic -integral: The conditions that the surfaces + and − are traction-free and the crack growth direction is along the -direction are used in deriving the above equation.

Fracture energy and crack speed
In this study, a series of SENT plane strain numerical specimens (see Fig. 1) is tested with different loading velocities and different specimen sizes. The considered cases of anḋare listed in Table 1. The maximum nominal strain rate investigated is 250/s, which is intermediate compared with high strain rate testing, such as in plate impact tests where the strain rates up 10 8 /s to have been reported [40]. Dimensions are normalized with respect to the length of a single RVE denoted = 0.02856 mm. For each case, the geometry of the SENT specimen satisfies that 0 = 0.05 and = 4 . The number of RVEs in the microstructure is kept fixed at RVE = 2 × 20 except where mentioned otherwise (i.e. 2 rows and 20 columns of RVEs). The initial notch tip has a distance of 1.33 from the left edge of the embedded microstructure zone. The fiber is modeled as a linear elastic material with the elasticity parameters listed in Table 2 and the material parameters for the VE-VP polymeric model are listed in Table 3. Considering that there might exist a characteristic time for a local damage to occur and this time scale is related to the relaxation times of the matrix as postulated in Shirinbayan et al. [5], it is ensured that the time steps adopted in the numerical simulations are much smaller than the relaxation times. The elastic properties corresponding to the homogenized medium outside of the embedded region determined by a computational homogenization technique are included in Table 2. Matrix cracking and fiber/matrix debonding are considered with the cohesive zone model with the shifted cohesive law described in Section 2.2. Considering that the fiber/matrix interface is generally weaker than the pure matrix, a smaller cohesive strength and fracture energy are adopted for the interface (see Table 4). The energy release rate is equal to the energy flux into the crack tip, divided by the crack speed [41,42]. The mode-I energy release rate for dynamic crack growth in composites can be computed by the dynamic -integral formulation introduced in Section 2.4 for all the considered SENT specimens. The crack speed is the time derivative of the crack length which can be determined from the numerical model. Since the dynamic energy release rate is equal to the dynamic fracture energy for a propagating crack, a relationship between dynamic fracture energy and crack speed can be established.
In this section, the crack growth process of the SENT specimen under dynamic loading for a typical case with specimen width = 600 and loading velocitẏ= 0.1 m∕s is first described. Then, the influence of the size of the embedded microstructure on the crack growth and energy release rate is discussed. Finally, the energy release rate for different crack speeds extracted from the numerical tests are presented.

Typical observations
For a typical test case with = 600 , 0 = 0.05 and = 4 , the SENT specimen is subjected to a loading velocitẏ= 0.1 m∕s. A total number of 113550 six-node triangular elements is used for the discretization of the numerical sample with a transition from a mesh size of 2 mm to 0.001 mm. Fig. 6 shows the initiation and evolution of cohesive cracks and the distribution of the normal stress of the material near the crack tip for five different time steps. It is found that the applied loading causes the typical plane-strain crack tip stress field with peanut shaped stress concentration. Inside the microstructure, an inhomogeneous stress distribution is found.
Cracks are formed in the fiber/matrix interfaces in a number of spots near the crack tip rather than the pure matrix (see Fig. 6(a)), which is due to lower cohesion strength at fiber/matrix interfaces. The spots where cracks initiate are sparsely distributed near the crack tip due to the inhomogeneous stress distribution caused by the applied dynamic loads and material inhomogeneity of the microstructure.
The material near the crack tip is experiencing complex conditions with interaction between dynamic loading, structural inertia, material nonlinearity and material failure. Importantly, the applied continuous loading generates a loading wave propagating into the structure, while the newly created crack surface unloading waves are generated. The process is also complicated due to structural inertia effects. Therefore, the material near the crack tip including the cohesive surfaces can experience several loading/unloading cycles, as visible in the change of loading/unloading state of the cracks shown in Fig. 6(a-c). Finally, a fully developed cohesive zone is formed and a dominant crack close to the mid-plane propagates in a self-similar manner (see Fig. 6(d-e)). Many cracks including both fiber/matrix debonding and matrix cracking are formed ahead of the crack tip while the cracks at the wake are unloading. Due to the inhomogeneous distribution of fibers, the growing crack is not straight but shows a certain tortuosity. The deformation of the microstructure is relatively small with a large fracture process zone.

Size of the microstructure
The adopted number of embedded RVEs in the embedded zone is RVE = 2 × 20. Justification for this choice was found in a size dependence study which is presented in this section. A study on the influence of the microstructure size on the dynamic The crack tip is defined as the appearance of the first fully damaged cohesive element with stress free surface ( = 1), as illustrated in Fig. 6e (Left). The time derivative of the crack length is the crack speed. To obtain the dynamic energy release rate of the fracture process zone, the path of the -integral is defined outside of the embedded microstructure zone. The path-dependence of the dynamic -integral is first investigated for three different prescribed paths, A, B and C shown in Fig. 7 for the case with RVE = 2 × 20. Path A is the outer boundary of the microstructure while path B is slightly further away from Path A and path C is even further than path B. Fig. 8 shows the dynamic -integral value vs. time for the three paths. It is seen that there are only very minor differences in the dynamic -integral value for the three different paths, which means that the path-independence is found for paths defined outside the embedded microstructure. It is noted that in the homogenized region where paths A, B and C are defined, the material response is modeled as elastic, which attributes to the path-independence observed here. Path A is therefore chosen as the -integral contour used in this study. Fig. 9 shows the crack extension vs. time and dynamic -integral for the four cases with different microstructure sizes. It can be observed that the crack extension curve for the four cases are not exactly the same, which is related to the fact that the location where crack occurs is not the same. However, the differences between the four cases are limited. The crack speed, i.e. the slope of the time vs. crack extension curve, seems very close among the four cases. The evolution of the dynamic for the four cases is also very close. This shows that by either increasing the number of fibers in -direction from 80 to 150 or in -direction from 10 to 20, the crack growth process is not evidently different. Therefore, the microstructure with 2 × 20 RVEs is selected in this study as it is most computationally efficient and provides size-independent responses.

Dynamic energy release rate
The dynamic energy release rate for the series of SENT tests listed in Table 1 is summarized in this section. By tracking the crack tip location during the 10 SENT tests, the crack extension is measured and shown in Fig. 10. The discrete time vs. crack extension data obtained in the numerical tests is fitted with smooth functions (e.g. exponential) using the curve fitting toolbox of MATLAB so that a smooth time vs. extension curve is obtained to compute the crack speed. It is observed that all fitted curves for the 10 cases have a 2 value larger than 0.9928, indicating that good fits are obtained. The crack speed is defined as the slope of the fitted curve, i.e. = ∕ (see Fig. 10(j)). It is noted that using numerical differentiation of the discrete time and crack extension points is not a good choice for defining the crack speed. Fig. 11 shows a comparison of the computed crack speed history with two numerical differentiation schemes and the chosen approach. The considered numerical differentiation methods are the mid-point scheme and the Euler backward scheme. For the mid-point scheme, the crack speed at discrete time ( = 1, 2, … , ) is calculated as: = ( +1 − −1 )∕( +1 − −1 ) in which is the total number of time instants. For the Euler backward scheme, the crack speed at discrete time is calculated as: = ( − −1 )∕( − −1 ). It is found that both numerical differentiation schemes give oscillating crack speeds, while the crack speed computed by the chosen approach shows a smooth crack speed history. These oscillations are not necessarily physical and may be attributed to the numerical discretization, numerical time stepping or to the microstructure. In any case, the homogenized response is of particular interest rather than the exact crack speed inside the embedded cell. Therefore, Y. Liu et al.  with the chosen approach the crack speed history is smooth and physically regarded as the average crack speed found in the SENT tests [44]. Fig. 12 shows the dynamic -integral value for different crack speeds extracted from the series of numerical tests. A number of observations are made: (1) for the first six cases of with = 70 , 100 anḋ= 0.01, 0.1, 1.0 m∕s, there appears to be a unique relation between the dynamic -integral and the crack speed , i.e.
(2) the dynamic -integral value when crack propagation starts, i.e. when crack speed > 0, among those six cases have small differences and an average value around 0.045 N/mm is identified. This value is between the fracture toughness of the matrix, 0.09 N/mm, and that of the fiber/matrix interface, 0.02 N/mm. (3) for the other four cases, the dynamic -integral value for the same crack speed is higher than that of the first six cases. (4) the case with = 600 anḋ= 0.01 m∕s has the largest dynamic -integral value. If a strain rate definition oḟ=̇∕2 is employed, the case with = 600 anḋ= 0.01 m∕s also has the lowest strain rate.  is in the same magnitude as the dynamic crack growth in carbon/epoxy composites with dynamic double cantilever beam tests in [45,46].

Discussions of mechanisms
The underlying mechanisms for the observations of the ( ) relation shown in Fig. 12 are discussed in this section, including inertia effect, ductile/brittle failure type and the R-curve effect.    Fig. 13 shows the time evolution of the dynamic -integral value before crack initiation for three cases, namely the loading velocitẏ= 0.01, 0.1, 1.0 m∕s with the same width = 100 . The dynamic -integral is gradually increasing for the three cases as a result of applied continuous displacement loading. The case with the lowest loading velocitẏshows a very smooth profile and a quadratic relation exists between the dynamic -integral and time. By increasing the loading velocity to 0.1 m/s and 1.0 m/s, the dynamic -integral clearly shows high frequency oscillations, which is due to an evident effect of system inertia activated by a larger test rate. Similar trends of the dynamic energy release rate for mode-I cracking in composites have been found in Liu et al. [10] with the interfacial thick level set (ITLS) approach.  The formulation for the dynamic -integral in Eq. (38) can be rewritten as:

Inertia effect
in which three different contributing components can be identified as , and . The is the formulation used for quasi-static loading. The represents the contribution of kinetic energy flow into the fracture process zone. The is zero in the special case of a constant crack propagation speed and steady-state crack propagation [39] and a nonzero value shows the deviation from that condition. For the case witḣ= 1.0 m∕s, the different components of the dynamic -integral are shown in Fig. 14. The dynamic -integral is close to the value of the quasi-static -integral although somewhat oscillatory. Compared with the , the and are much smaller. Of these two, gives the larger contribution. This shows that the effect of nonsteady-state crack propagation or non-constant crack speed is causing a significant inertia effect. However, the inertia effect is not dominant. A similar observation was made by Nakamura et al. [47] who reported for a three-point-bending test that the inertia effect is minor when the kinetic energy is a small fraction of the strain energy of the system. Therefore, an average fit of the ( ) data for this case as included in Fig. 12 is a reasonable representation of the material response. For a lower test rate witḣ= 0.01 m∕s, the , , and are shown in Fig. 15. The dynamic -integral is almost equal to the static -integral while the other two components and are negligible. This shows that for lower rates, the inertia effect vanishes. It is expected that for higher strain rate testing, for instance, Hopkinson Bar loaded fracture experiments, inertia effects become evident and need to be filtered out [48].

Failure type
The failure mode in the embedded microstructure zone is found to be brittle failure with limited plasticity. Fig. 16 shows the dissipation for plasticity and cohesive cracks for three cases witḣ= 0.01, 0.1, 1.0 m∕s and the same = 100 . They are calculated by Eqs. (19) and (31) and normalized with 0 , where = 0.02 N/mm is the fracture energy of fiber/matrix interface, is the RVE size and 0 = 1.0 mm is a unit thickness. It is observed that the plastic deformation in the matrix dissipates much less energy than the cohesive cracks. For instance, at the same amount of crack extension 0.1 mm, the case with loading velocitẏ= 0.1 m∕s has a plastic dissipation of 4.332×10 −4 N⋅mm, while the dissipation for cohesive crack at that point in time is 5.878×10 −4 N⋅mm. A comparison of the three cases shows that the case witḣ= 1.0 m∕s has the largest cohesive dissipation while the case witḣ = 0.01 m∕s has the largest plastic dissipation. This shows that for a lower loading velocity, plasticity is more developed while for higher loading velocity cohesive cracks dissipate more energy. This phenomenon represents the commonly referred ductile-to-brittle transition for increasing loading rate [9,49,50]. Nevertheless, failure in the fracture process zone for the case with lower loading velocity is still rather brittle. This can be seen from the distribution of the equivalent plastic strain for the case oḟ= 0.01 m∕s shown in Fig. 17. The plastic strain is limited to the area close to the crack path. The distribution of normal strain shown in Fig. 18 for the case with = 70 anḋ= 0.1 m∕s reveals that the strain near the crack tip remains small, which indicates a brittle failure.

Dynamic -curve
It is known that in laminated composites, due to fiber cross-over bridging behind the crack tip the fracture toughness can increase for a certain distance of crack growth. The increase of the apparent fracture toughness with crack extension is usually described by a function of crack growth resistance vs. crack extension, i.e. the so-called -curve [51]. Even though the crack propagation in this study takes place in the transverse plane, an -curve also exists here. This is related to the development of the fracture process zone, where microcracking at fiber/matrix interfaces and inside the pure matrix and tearing of matrix ligaments are found. In Fig. 19 the -curve is shown by plotting the evolution of the dynamic -integral as a function of crack extension. The -curve for all cases shows a rising trend.
Except for the case with = 70 anḋ= 1.0 m/and the case with = 100 anḋ= 1.0 m∕s, all cases follow approximately the same -curve. Considering that structural inertia, cohesive cracking and rate-dependent plasticity are coexisting during the development of the fracture process zone, there are minor differences of the -curve of the different cases, although, as mentioned in Section 4.2, rate-dependent plasticity is less pronounced and does not contribute much to these differences. In the cases with = 70 anḋ= 1.0 m∕s and = 100 anḋ= 1.0 m∕s, which are the two case with the highest nominal strain rates, oscillations are present in the -curve which are ascribed to inertia effects.
In Fig. 20, the distribution of normal stress in -direction is shown for two lower rate cases, one with = 70 anḋ= 0.1 m∕s and the other with = 600 anḋ= 0.01 m∕s. For both cases, two typical time instants are selected. It shows that the development of the failure zone is a gradual process with the formation of microcracks (along fiber/matrix interfaces and inside the matrix) and ductile tearing of matrix ligaments. The two cases form a very similar crack pattern when the fracture process zone is fully developed and a periodic crack pattern forms as a result of the periodicity of the embedded microstructure. The damage of cohesive cracks corresponding to Fig. 20(b) and (d) is plotted in Fig. 21 when the crack length of both cases is 0.072 mm. The case with higher loading rate only has a slightly wider spreading of cohesive cracks. The similarity of the -curve for lower rates offers an explanation for the shift observed in the plot of -integral versus crack speed in Fig. 12 for the same cases. Since there is a one-to-one relation between crack length and applied load for these cases, different crack velocities must be found for different applied loading rates.
For the higher rate cases, there are oscillations in the -curve as well as an increase in the overall fracture resistance. The differences in dynamic -integral is a numerical representation of the velocity-toughening effect that has been observed experimentally for quasi-brittle materials. Zhou et al. [52] found the failure mechanism of a PMMA plate was found to display increasingly rough crack surfaces for increasing crack propagation velocities. As seen in Figs. 10 and 12, the crack speed for the higher nominal strain rate cases is larger than that of lower nominal strain rate cases. Fig. 22 shows a comparison of the dissipation of cohesive cracks for three cases, representing the lowest loading rate and the two highest loading rates. It is observed that the higher rate cases have significantly larger cohesive dissipation, pointing at more damage in secondary microcracks.

Conclusion
In this paper, a multiscale numerical framework is established to evaluate the fracture energy of dynamic crack propagation in composites. A series of numerical simulations with different specimen sizes and different loading velocities is performed to simulate the deformation and failure process of the SENT specimen with embedded composite microstructure when subjected to continuous dynamic loading. Instead of running explicit dynamics analyses with accumulated divergence from balance of momentum, an implicit dynamics solution scheme is adopted. For each time step the dynamic version of -integral is evaluated as a measure for the dynamic fracture energy.
The introduced numerical framework allows for a quantitative evaluation of the dynamic fracture energy of composites and for analysis of how rate-dependent plasticity, distributed microcracking and inertia effects contribute to the observed fracture energy.  For all considered cases, microcracks are initially formed in the fiber/matrix interfaces in a number of locations near the crack tip. Materials near the crack tip including the newly created crack surfaces experience a complicated loading process mainly due to interaction of dynamic loading, structural inertia and material failure.
With the given material parameter set, it is seen that an increase of the applied strain rate gives rise to a trend of increasing brittleness for the failure of composites with reduced plastic energy dissipation. However, even for cases with low loading velocity and large specimen size, failure is relatively brittle with a small amount of plasticity occurring near the crack tip and in the wake of the fracture process zone. Therefore, the influence of plasticity on global rate-dependence remains limited.
The dynamic fracture energy shows an increasing trend with the crack speed but no unique ( ) relation is found [45]. For most investigated loading rates, the crack growth follows a rate-independent -curve, which is related to the fact that the amount of cohesive energy dissipation is the same for different cases with lower loading rates.  Cases with high nominal strain rate show visible inertia effects with oscillating values of . These cases also show increased microcracking which leads to a higher overall energy dissipation pointing at a velocity toughening effect. Considering that the mode-I fracture toughness is also temperature-dependent [53], incorporation of the temperature effect into the current numerical framework is further needed and can be done by following the idea of [16].

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