# A Modified Strain-Rate-Dependent Spring-Mass Model for Response Prediction of Composite Laminates Subjected to Low-Velocity Impact

Authors: Qiao Zhu; Zhefeng Yu
Journal: International Journal of Aerospace Engineering
 (2023)

---

## Abstract

A modified strain-rate dependent (SRD) spring-mass model was first developed to capture the strain rate effect on response prediction of fiber-reinforced polymer (FRP) laminates subjected to low-velocity impact (LVI). The constitutive relations of FRP material were modified with the impact-induced strain rate using Yen-Caiazzo’s function in space and time dimensions simultaneously. The stiffness coefficients of the spring-mass model were updated step-by-step during solving the Bubnov-Galerkin equation by the variational method, allowing the SRD contact history to be obtained via recursive integration. The SRD expressions of stiffness coefficients under four typical boundary conditions were presented. Drop-weight tests on FRP laminates and corresponding VUMAT-based LVI simulation with the finite element method (FEM) were provided to prove the validity of the proposed SRD model in evaluating the strain rate effect on the LVI response of composites. Further parameter studies were carried out to investigate the reliability of the proposed SRD model for evaluating the influence of reinforced fibers with different strain rate dependency on impact response. The proposed lumped parameter model has been proven to be more efficient than the traditional FEM, which can be combined with some existing damage models to accurately analyze the delamination evolution process of FRP laminates under LVI in further studies.


---

## Body

## 1. Introduction

As fiber-reinforced polymer (FRP) laminates are increasingly used in aeronautical engineering, their relatively weak resistance to impact has attracted substantial attention in recent years [1–3]. The analysis of impact on FRP laminates mainly evaluates the contact behavior and the structural response. The coupling effect between these two aspects can be described using specific impact models from simplified analytical approaches to more precise finite element method (FEM) [4, 5]. However, it has been a challenge to develop models that provide adequate accuracy with reasonable computation time and a wide selection of abstraction scales using numerical FEM. As such, the derivation of the equivalent reduced analytical impact models becomes more attractive to conveniently design FRP laminate structures [6].The study of foreign object impact on FRP laminate structures mainly involves the prediction of contact behavior and structural response [7], recognizing the onset and assessing the extent of impact-induced damage [8], and evaluating the residual properties of the structure [9]. The first step of impact analysis is the most basic and important which requires the evaluation of the motion of the impactor, the response of the structure, and the indentation of the contact area. The development of analytical impact model has generated a great deal of interest and high level of research activities as evidenced by the numerous articles [4, 6], especially on low-velocity impact (LVI) by projectile with relatively “large mass” common in engineering field. The spring-mass model, first presented by Shivakumar et al. [10], is a simple and accurate model, which provides the contact history for most LVI scenarios. The effects of contact, bending shear, and membrane on structural response are represented by the corresponding stiffness in this model. For impacts with negligible indentation but large structural deflection and significant membrane stiffening, the spring-mass model can be simplified to a nonlinear single-degree-of-freedom (SDOF) version [6, 10]. The numerical solution of the motion equation yields accurate estimations of contact history and the impact response of the structure.Many existing researches on the impact behavior of FRP laminate structures have been conducted under the condition of quasi-static deformation. However, even at a low initial velocity, the impact would easily induce a rapid increase in strain at the surrounding area of the impact point, especially in the early stage of contact [11]. The matrix-dominated properties of laminates have been proven to be strain-rate sensitive due to the viscoplastic characteristic of the resin [12], but no consistency is observed for fibers [13]. Glass fiber [14, 15] and Kevlar fiber [16] are found to be strain-rate dependent (SRD). In contrast, the elastic tensile modulus of polyvinyl-alcohol (PVA) fiber [17] remains unchanged with the strain rate, but its mechanical strength properties, such as ultimate stress and failure strain, are SRD. For carbon fiber, the strain rate exhibits no effect on both elastic and plastic mechanical properties [18–20].As the modulus and strength of the unidirectional lamina increase with the impact-induced strain rate, the mechanical properties of FRP laminates in dynamic and static loading conditions differ significantly [21, 22]. Dynamic constitutive relations are crucial for the design of FRP laminate structures under impact loadings to prevent the introduction of non-negligible errors in the structural response prediction [23], which are often overlooked in previous studies. Therefore, it is essential to consider the strain rate effect when developing an analytical model for the accurate prediction of impact response and consequent damage evolution. The constitutive relations between strain rate and material mechanical properties can be modeled using SRD functions [24, 25], such as Cowper-Symonds’ power function [26, 27], Johnson-Cook’s (J-C) logarithmic function [28, 29], and Yen-Caiazzo’s (Y-C) logarithmic function [30], which were all experimental data-based empirical models. These models were all built through multiplying static properties by a dynamic increase function.Due to the strain rate effect on material mechanical properties vary with time and distance from the impact position, which cannot be directly represented by the explicit expressions given by Shivakumar [10]. These expressions were applied directly by later studies without considering the variation of elastic modulus caused by the strain rate effect [6, 31–33]. Most existing studies considering the strain rate effect on impact response of laminate are FEM-based [34–37], indicating that there is still a lack of analytical solutions for versatility and efficiency. This study mainly focused on the influence of strain rate on the LVI response of FRP laminates. The spring-mass model was modified with the dynamic stiffness coefficients obtained by considering the strain rate effect on material mechanical properties. The derivation of the bending and membrane stiffness was retraced, and the expressions were presented in terms of integrals in this study. Then the dynamic stiffness can be obtained by numerical recursive integrations. The strain rate effect on mechanical properties of the laminate and impact response of the laminate was evaluated with time and space simultaneously.The remainder of this paper is organized as follows: the basic formulations and implementation flow of the proposed SRD model are presented on a circular orthotropic FRP laminate subjected to LVI in Section2. Setup of the experiment and establishment of the VUMAT-based finite element model for the comparative validation of the proposed SRD model are introduced in Section 3 and Section 4, respectively. In section 5, convergence analysis of the circumferential integration step and the time step is performed to determine the appropriate value of these two uncertain parameters of the proposed SRD model. The proposed model is verified by experiments with three different impact energies, and the efficiency compared to the FEM is investigated. Further, parametric studies for strain rate dependency of the reinforced fiber on LVI response prediction of FRP laminates are performed. The main conclusions and prospects of future research are summarized in the last section.


## 2. Theory

### 2.1. Strain-Rate Independent (SRI) “Large Mass” Impact Model

Consider a circular orthotropic thin FRP laminate of dimensionsRp×hp subjected to vertical impact at the center by an elastic sphere with an initial velocity V0. For convenience of expression, the cylindrical coordinate system orφz was set coincident with the midplane of the laminate, as shown in Figure 1. It has been experimentally proven that when the impactor/structure mass ratio is larger than 2, the contact duration of such “large mass” impacts would be significantly longer than the time that the flexural waves take to reach the structure boundary and reflect. The response of the laminate is dominated by quasi-static deformation under a concentrated load [38].Figure 1
Circular FRP laminate impacted by an elastic sphere at the center.As the load distribution, structural deflection, and boundary conditions are all symmetric about the center of the axisymmetric bending problem, the deflection equations of the circular laminate are independent of the angleφ. In the absence of material damping and surface friction, the deflection of the point at a distance r from the center on the midplane of the laminate wr are given by [39].
(1a)wr=3FcRp21−μe24πEehp31−r2Rp2+2r2Rp2lnrRpclamped,(1b)wr=3FcRp21−μe24πEehp33+μe1+μe1−r2Rp2+21+μe3+μer2Rp2lnrRpsimplesupported,where Fc is the contact force between the impactor and laminate. Ee and μe are averaged homogenized modulus and Poisson ratio, respectively, which are denoted as [32]
(2)Ee=12π∫02πEφdφμe=12π∫02πμφdφ,where Eφ and μφ are the equivalent elastic modulus and Poisson ratios of the laminate at angle φ. These material property coefficients can be obtained by applying the material and coordinate transformation matrices to the constitutive relations of the FRP material [40]
(3)Q¯φk=CTφMTθ−1QkMTθ−TCTφT,where Q¯φk is the transformed reduced stiffness matrix of the kth lamina at angle φ in cylindrical coordinates. Qk is the reduced stiffness matrix of the kth lamina in material coordinate system [41]. As shown in Figure 2, MTθ is the material transformation matrix from lamina to laminate in Cartesian coordinate system oxyz. θ is the angle from the x-axis counterclockwise to the fiber direction. CTφ is the coordinate system transformation matrix from Cartesian coordinate system to cylindrical coordinate system. φ is the angle from the x-axis counterclockwise to the r-axis.Figure 2
Transformations between coordinate systems and lamina fiber direction.The extensional stiffness coefficients of the laminate at angleφ can be obtained by stacking the transformed reduced stiffness matrix along the laminate thickness
(4)Aijφ=∑k=1NQ¯ijkφhk,where hk is the thickness of the lamina.The equivalent elastic modulus of the laminate in cylindrical coordinates at angleφ can be obtained using
(5)Eφ=1hpα11φ,where αφ=Aφ−1.Based on the deflection equations as shown in Equation (1a), the Bubnov-Galerkin variational equation derived from the principle of virtual displacement is formulated to establish a relationship between the contact force and the central deflection of the laminate [42]
(6)∫0RpEehp3121−μe2ddr∇2wr−Ψ−hprdΦdrdwrdrdηidrrdr=0,where Ψ and Φ are the load and stress functions, respectively. ηi represents the radius-related term in the series expansion of the deflection equation, such as the part in square brackets of Equation (1a). ∇2 is the rotation-angle independent Laplacian operator in the cylindrical coordinates, which can be expressed as
(7)∇2=1rddrrddr.Take the peripheral edge-moveable clamped (EMC) circular laminate acted with a concentrated load as an example. It can be seen from Equation (1a), when r=0, the maximum deflection of the laminate (at the loading point) can be represented by the terms outside the square bracket on the right-hand side of the equation
(8)wmax=w0=3FcRp21−μe24πEehp3.As defined in [42], Ψ, Φ, and ηi can be denoted as
(9a)η1=1−r2Rp2+2r2Rp2lnrRp,(9b)Ψ=Fc2πr,(9c)ddr∇2Φ=−Ee2rdwrdr2.The derivative ofΦ on the right side of Equation (6) can be obtained by applying the following integral formula to Equation (9a).
(10)∫ρmlnnρdρ=ρm+1lnnρm+1−nm+1∫ρmlnn−1ρdρ.The result of integration by parts is given as [42]
(11)dΦdr=−4Eewmax2Rpr34Rp3ln2rRp−3r38Rp3lnrRp+7r332Rp3−7r32Rp.Substituting Equations. (1a), (9a), and (11) into Equation. (6), the Bubnov-Galerkin variation equation is expressed as a function of the contact force Fc and the laminate deflection at the loading point w0 as
(12)w3016hpEeRp6∫0Rpr5Rp2ln4rRp−3r52Rp2ln3rRp+7r58Rp2ln2rRp−7r38ln2rRpdr+w08hp3Ee3Rp41−μe2∫0RprlnrRpdr=Fc2πRp2∫0RprlnrRpdr.Without considering the strain-rate effect, the averaged homogenized properties in Equation (12) remain constant throughout the entire contact process at all points of the structure. Equation (12) can be simplified though continuous integration on the coefficients of w30, w0, and Fc as
(13)191πEehp648Rp2w30+4πEehp331−μe2Rp2w0=Fc.Equation (13) provides the specific relationship between Fc and w0. From this, the spring-mass model was developed by Shivakumar et al. [10] to predict the contact history and estimate the structural response of the “large mass” impact. The schematic representation is shown in Figure 3, where the superscripts i and p of mass m and displacement w represent the impactor and FRP laminate, respectively. The indentation α can be denoted as the difference between the displacement of the impactor and laminate, where α=wi−wp.Figure 3
Schematic representation of the spring-mass model: (a) physical quantities and effects in the impact problem. (b) TDOF impact model, and (c) SDOF impact modelDue to the local indentation is generally negligible (α≈0) compared to the structural deflection [6] and the total plate mass is significantly smaller than the impactor mass (mp≪mi), the originally developed two-degree-of freedom (TDOF) impact model (Figure 3(b)) can be reduced to an SDOF system (Figure 3(c)) to improve efficiency while ensuring the accuracy. The motion equation of such model is given as
(14)mi∂2wp∂t2+kbswp+kmwp2=0,where kbs is denoted by the combination of bending stiffness kb and shear stiffness ks as kbs−1=kb−1+ks−1. km is the nonlinear membrane stiffness. For a thin plate structure (large span to thickness ratio), the relatively low shear stiffness ks can be neglected [43] reducing the combination spring constant kbs to kb.By comparing the coefficients of terms with the same power on both sides, the specific expressions of stiffness in Equation (14) can be obtained from Equation (13) correspondingly. The expressions of these stiffness coefficients under simply supported boundary are also given by Shivakumar et al. [10] The displacement of the FRP laminate at the impact point wp corresponds to the maximum deflection of the laminate w0. The contact force Fc is denoted by the product of the impactor mass and its acceleration as:
(15)Fc=mi∂2wi∂t2=mi∂2wp∂t2.


### 2.2. Strain-Rate Dependent (SRD) “Large Mass” Impact Model

The expressions of bending and membrane stiffness for circular composite plates given by Shivakumar et al. [10] generally remain constant during the entire contact process, which makes the impact response different from the realistic situation due to the sensitivity of laminate material properties to impact-induced strain rate variations. In dynamic conditions, the dependence of laminate material properties on the strain rate is nonnegligible. In addition, the strain rate varies with the distance from the laminate center during the impact process. Therefore, the expressions of the stiffness coefficients in Equation (14) should be modified with the variation of strain rate in space and time dimensions simultaneously.The bending strain of the laminate is a function of the thickness, thus the strain component of thekth lamina at the nth time step can be expressed as
(16)εkr,tn=−zk∂2wr∂r2,where zk is the coordinate of the middle surface of the kth lamina.Due to the nonlinearity of Equation (14), it can only be solved using step-by-step recursive numerical method. On the time dimension, the strain rate of the kth lamina ε̇k can be obtained by taking the difference between the strain values at corresponding time step and before
(17)ε̇kr,tn=εkr,tn−εkr,tn−1Δt,where Δt is the length of time step.In this study, without considering the slight effect of strain hardening and thermal softening, only the strain rate effects were considered. The Y-C model [30], which is the simplified version of the general J-C model [28], was adopted to modify the material properties of the laminate under the dynamic stress state, which was given as
(18)Ek∗r,tn=Emelogε̇kr,tnε̇0+1,where Ek∗ is the SRD modulus of the kth lamina. E is the quasi-static modulus obtained from the standard test. me is the strain rate coefficient obtained from the empirical fitting of the experimental data. Δε0 is the reference strain rate, which was set to 10−4 for a quasi-static stress state [23].After updating the material properties in time dimension, the derivation of the stiffness coefficients in Equation (14) should be rededuced through integration along the entire radius due to the variation in strain rate on the spatial scale. In the numerical integration simplification of Equation (12) to Equation (13), Ee and μe are no longer the constant. The averaged SRD modulus Ee∗ can be obtained by substituting Equation (18) into Equations (2)-(5). Because no explicit expressions of these averaged effective properties can be obtained, the bending and membrane stiffness can only be expressed in terms of integrals and calculated numerically as
(19a)kbr,tn=∫0Rp8Eer,tnhp3/3Rp41−μe2Vbdr∫0Rp2/πRp2VFdr,(19b)kmr,tn=∫0Rp16Eer,tnhp/Rp5Vmdr∫0Rp2/πRp2VFdr,where Vb, Vm, and VF are the radius correlation coefficients which depend on the boundary conditions of the laminate. The expressions of these coefficients are presented in Table 1 along with the other three typical boundary conditions, namely edge-immoveable clamped (EIC), edge-moveable simple-supported (EMSS), and edge-immoveable simple-supported (EISS).Table 1
Expressions ofVb, Vm, and VF with different boundary conditions.
Boundary conditionsVbVmVFEdge-moveable clampedrlnr/Rpr5/Rp3ln4r/Rp−3r5/2Rp3ln3r/Rp+7r5/8Rp3ln2r/Rp−7r3/8Rpln2r/Rprlnr/RpEdge-immoveable clampedrlnr/Rpr5/Rp3ln4r/Rp−3r5/2Rp3ln3r/Rp+7r5/8Rp3ln2r/Rp−r3/8Rp9−7μ/1−μln2r/Rprlnr/RpEdge-moveable simple-supportedr1+μ/3+μ21+μlnr/Rp−113+μ41+μlnrRp−12r5Rp3−r3Rp−21+μr5Rp3lnrRp−34r5Rp3−34r3Rp+1+μ2r5Rp3ln2rRp−32r5Rp3lnrRp+78r5Rp3−78r3Rpr/3+μ1+μlnr/Rp−1Edge-immoveable simple-supportedr1+μ/3+μ21+μlnr/Rp−113+μ41+μlnrRp−12r5Rp3−3−μ1−μr3Rp−21+μr5Rp3lnrRp−34r5Rp3−5−3μ41−μr3Rp+1+μ2r5Rp3ln2rRp−32r5Rp3lnrRp+78r5Rp3−9−7μ81−μr3Rpr/3+μ1+μlnr/Rp−1For rectangular laminate of dimensionsLx×Ly, equivalent radius Rp∗ in Equations (19a) and (19b) can be determined from the circular plate with equal effective area: Rp∗=LxLy/π [44].By substituting Equation (19) into Equation (14), the modified SDOF impact model is established, which can be solved using the Runge-Kutta method. The contact history and structural response at the impact point can be obtained using point-by-point numerical integration with step-by-step recursive integration. The one-cycle flow chart of the developed SRD “large mass” impact model is shown in Figure 4.Figure 4
Flow chart of the developed SRD “large mass” impact model.


## 2.1. Strain-Rate Independent (SRI) “Large Mass” Impact Model

Consider a circular orthotropic thin FRP laminate of dimensionsRp×hp subjected to vertical impact at the center by an elastic sphere with an initial velocity V0. For convenience of expression, the cylindrical coordinate system orφz was set coincident with the midplane of the laminate, as shown in Figure 1. It has been experimentally proven that when the impactor/structure mass ratio is larger than 2, the contact duration of such “large mass” impacts would be significantly longer than the time that the flexural waves take to reach the structure boundary and reflect. The response of the laminate is dominated by quasi-static deformation under a concentrated load [38].Figure 1
Circular FRP laminate impacted by an elastic sphere at the center.As the load distribution, structural deflection, and boundary conditions are all symmetric about the center of the axisymmetric bending problem, the deflection equations of the circular laminate are independent of the angleφ. In the absence of material damping and surface friction, the deflection of the point at a distance r from the center on the midplane of the laminate wr are given by [39].
(1a)wr=3FcRp21−μe24πEehp31−r2Rp2+2r2Rp2lnrRpclamped,(1b)wr=3FcRp21−μe24πEehp33+μe1+μe1−r2Rp2+21+μe3+μer2Rp2lnrRpsimplesupported,where Fc is the contact force between the impactor and laminate. Ee and μe are averaged homogenized modulus and Poisson ratio, respectively, which are denoted as [32]
(2)Ee=12π∫02πEφdφμe=12π∫02πμφdφ,where Eφ and μφ are the equivalent elastic modulus and Poisson ratios of the laminate at angle φ. These material property coefficients can be obtained by applying the material and coordinate transformation matrices to the constitutive relations of the FRP material [40]
(3)Q¯φk=CTφMTθ−1QkMTθ−TCTφT,where Q¯φk is the transformed reduced stiffness matrix of the kth lamina at angle φ in cylindrical coordinates. Qk is the reduced stiffness matrix of the kth lamina in material coordinate system [41]. As shown in Figure 2, MTθ is the material transformation matrix from lamina to laminate in Cartesian coordinate system oxyz. θ is the angle from the x-axis counterclockwise to the fiber direction. CTφ is the coordinate system transformation matrix from Cartesian coordinate system to cylindrical coordinate system. φ is the angle from the x-axis counterclockwise to the r-axis.Figure 2
Transformations between coordinate systems and lamina fiber direction.The extensional stiffness coefficients of the laminate at angleφ can be obtained by stacking the transformed reduced stiffness matrix along the laminate thickness
(4)Aijφ=∑k=1NQ¯ijkφhk,where hk is the thickness of the lamina.The equivalent elastic modulus of the laminate in cylindrical coordinates at angleφ can be obtained using
(5)Eφ=1hpα11φ,where αφ=Aφ−1.Based on the deflection equations as shown in Equation (1a), the Bubnov-Galerkin variational equation derived from the principle of virtual displacement is formulated to establish a relationship between the contact force and the central deflection of the laminate [42]
(6)∫0RpEehp3121−μe2ddr∇2wr−Ψ−hprdΦdrdwrdrdηidrrdr=0,where Ψ and Φ are the load and stress functions, respectively. ηi represents the radius-related term in the series expansion of the deflection equation, such as the part in square brackets of Equation (1a). ∇2 is the rotation-angle independent Laplacian operator in the cylindrical coordinates, which can be expressed as
(7)∇2=1rddrrddr.Take the peripheral edge-moveable clamped (EMC) circular laminate acted with a concentrated load as an example. It can be seen from Equation (1a), when r=0, the maximum deflection of the laminate (at the loading point) can be represented by the terms outside the square bracket on the right-hand side of the equation
(8)wmax=w0=3FcRp21−μe24πEehp3.As defined in [42], Ψ, Φ, and ηi can be denoted as
(9a)η1=1−r2Rp2+2r2Rp2lnrRp,(9b)Ψ=Fc2πr,(9c)ddr∇2Φ=−Ee2rdwrdr2.The derivative ofΦ on the right side of Equation (6) can be obtained by applying the following integral formula to Equation (9a).
(10)∫ρmlnnρdρ=ρm+1lnnρm+1−nm+1∫ρmlnn−1ρdρ.The result of integration by parts is given as [42]
(11)dΦdr=−4Eewmax2Rpr34Rp3ln2rRp−3r38Rp3lnrRp+7r332Rp3−7r32Rp.Substituting Equations. (1a), (9a), and (11) into Equation. (6), the Bubnov-Galerkin variation equation is expressed as a function of the contact force Fc and the laminate deflection at the loading point w0 as
(12)w3016hpEeRp6∫0Rpr5Rp2ln4rRp−3r52Rp2ln3rRp+7r58Rp2ln2rRp−7r38ln2rRpdr+w08hp3Ee3Rp41−μe2∫0RprlnrRpdr=Fc2πRp2∫0RprlnrRpdr.Without considering the strain-rate effect, the averaged homogenized properties in Equation (12) remain constant throughout the entire contact process at all points of the structure. Equation (12) can be simplified though continuous integration on the coefficients of w30, w0, and Fc as
(13)191πEehp648Rp2w30+4πEehp331−μe2Rp2w0=Fc.Equation (13) provides the specific relationship between Fc and w0. From this, the spring-mass model was developed by Shivakumar et al. [10] to predict the contact history and estimate the structural response of the “large mass” impact. The schematic representation is shown in Figure 3, where the superscripts i and p of mass m and displacement w represent the impactor and FRP laminate, respectively. The indentation α can be denoted as the difference between the displacement of the impactor and laminate, where α=wi−wp.Figure 3
Schematic representation of the spring-mass model: (a) physical quantities and effects in the impact problem. (b) TDOF impact model, and (c) SDOF impact modelDue to the local indentation is generally negligible (α≈0) compared to the structural deflection [6] and the total plate mass is significantly smaller than the impactor mass (mp≪mi), the originally developed two-degree-of freedom (TDOF) impact model (Figure 3(b)) can be reduced to an SDOF system (Figure 3(c)) to improve efficiency while ensuring the accuracy. The motion equation of such model is given as
(14)mi∂2wp∂t2+kbswp+kmwp2=0,where kbs is denoted by the combination of bending stiffness kb and shear stiffness ks as kbs−1=kb−1+ks−1. km is the nonlinear membrane stiffness. For a thin plate structure (large span to thickness ratio), the relatively low shear stiffness ks can be neglected [43] reducing the combination spring constant kbs to kb.By comparing the coefficients of terms with the same power on both sides, the specific expressions of stiffness in Equation (14) can be obtained from Equation (13) correspondingly. The expressions of these stiffness coefficients under simply supported boundary are also given by Shivakumar et al. [10] The displacement of the FRP laminate at the impact point wp corresponds to the maximum deflection of the laminate w0. The contact force Fc is denoted by the product of the impactor mass and its acceleration as:
(15)Fc=mi∂2wi∂t2=mi∂2wp∂t2.


## 2.2. Strain-Rate Dependent (SRD) “Large Mass” Impact Model

The expressions of bending and membrane stiffness for circular composite plates given by Shivakumar et al. [10] generally remain constant during the entire contact process, which makes the impact response different from the realistic situation due to the sensitivity of laminate material properties to impact-induced strain rate variations. In dynamic conditions, the dependence of laminate material properties on the strain rate is nonnegligible. In addition, the strain rate varies with the distance from the laminate center during the impact process. Therefore, the expressions of the stiffness coefficients in Equation (14) should be modified with the variation of strain rate in space and time dimensions simultaneously.The bending strain of the laminate is a function of the thickness, thus the strain component of thekth lamina at the nth time step can be expressed as
(16)εkr,tn=−zk∂2wr∂r2,where zk is the coordinate of the middle surface of the kth lamina.Due to the nonlinearity of Equation (14), it can only be solved using step-by-step recursive numerical method. On the time dimension, the strain rate of the kth lamina ε̇k can be obtained by taking the difference between the strain values at corresponding time step and before
(17)ε̇kr,tn=εkr,tn−εkr,tn−1Δt,where Δt is the length of time step.In this study, without considering the slight effect of strain hardening and thermal softening, only the strain rate effects were considered. The Y-C model [30], which is the simplified version of the general J-C model [28], was adopted to modify the material properties of the laminate under the dynamic stress state, which was given as
(18)Ek∗r,tn=Emelogε̇kr,tnε̇0+1,where Ek∗ is the SRD modulus of the kth lamina. E is the quasi-static modulus obtained from the standard test. me is the strain rate coefficient obtained from the empirical fitting of the experimental data. Δε0 is the reference strain rate, which was set to 10−4 for a quasi-static stress state [23].After updating the material properties in time dimension, the derivation of the stiffness coefficients in Equation (14) should be rededuced through integration along the entire radius due to the variation in strain rate on the spatial scale. In the numerical integration simplification of Equation (12) to Equation (13), Ee and μe are no longer the constant. The averaged SRD modulus Ee∗ can be obtained by substituting Equation (18) into Equations (2)-(5). Because no explicit expressions of these averaged effective properties can be obtained, the bending and membrane stiffness can only be expressed in terms of integrals and calculated numerically as
(19a)kbr,tn=∫0Rp8Eer,tnhp3/3Rp41−μe2Vbdr∫0Rp2/πRp2VFdr,(19b)kmr,tn=∫0Rp16Eer,tnhp/Rp5Vmdr∫0Rp2/πRp2VFdr,where Vb, Vm, and VF are the radius correlation coefficients which depend on the boundary conditions of the laminate. The expressions of these coefficients are presented in Table 1 along with the other three typical boundary conditions, namely edge-immoveable clamped (EIC), edge-moveable simple-supported (EMSS), and edge-immoveable simple-supported (EISS).Table 1
Expressions ofVb, Vm, and VF with different boundary conditions.
Boundary conditionsVbVmVFEdge-moveable clampedrlnr/Rpr5/Rp3ln4r/Rp−3r5/2Rp3ln3r/Rp+7r5/8Rp3ln2r/Rp−7r3/8Rpln2r/Rprlnr/RpEdge-immoveable clampedrlnr/Rpr5/Rp3ln4r/Rp−3r5/2Rp3ln3r/Rp+7r5/8Rp3ln2r/Rp−r3/8Rp9−7μ/1−μln2r/Rprlnr/RpEdge-moveable simple-supportedr1+μ/3+μ21+μlnr/Rp−113+μ41+μlnrRp−12r5Rp3−r3Rp−21+μr5Rp3lnrRp−34r5Rp3−34r3Rp+1+μ2r5Rp3ln2rRp−32r5Rp3lnrRp+78r5Rp3−78r3Rpr/3+μ1+μlnr/Rp−1Edge-immoveable simple-supportedr1+μ/3+μ21+μlnr/Rp−113+μ41+μlnrRp−12r5Rp3−3−μ1−μr3Rp−21+μr5Rp3lnrRp−34r5Rp3−5−3μ41−μr3Rp+1+μ2r5Rp3ln2rRp−32r5Rp3lnrRp+78r5Rp3−9−7μ81−μr3Rpr/3+μ1+μlnr/Rp−1For rectangular laminate of dimensionsLx×Ly, equivalent radius Rp∗ in Equations (19a) and (19b) can be determined from the circular plate with equal effective area: Rp∗=LxLy/π [44].By substituting Equation (19) into Equation (14), the modified SDOF impact model is established, which can be solved using the Runge-Kutta method. The contact history and structural response at the impact point can be obtained using point-by-point numerical integration with step-by-step recursive integration. The one-cycle flow chart of the developed SRD “large mass” impact model is shown in Figure 4.Figure 4
Flow chart of the developed SRD “large mass” impact model.


## 3. Experiment

The proposed SRD model described above was validated in this section by comparing the contact history predicted with the results from drop-weight impact tests and corresponding VUMAT-based LVI FEM simulation. As the object of this study is to evaluate the strain rate effect on LVI response of composite independently and the modified SRD spring-mass model was built regardless of impact damage, the contact force during the drop-weight impact tests was controlled not to exceed the delamination threshold load (DTL) of the laminate.


### 3.1. Specimen

The square FRP laminate specimen with 125 mm on each side employed in the following experiments and simulations was fabricated by Guangwei Composites Material Co., Ltd. using the unidirectional carbon fiber prepreg USN20000 with a thickness of 0.2 mm. The ply configuration of the laminate was [(0/90)4]s. Table 2 gives the mechanical properties of the lamina. me was obtained by fitting test data with different impact velocities following the Y-C logarithmic function in Equation (18).Table 2
Lamina properties of the laminate specimen.
E11E22=E33G12=G13V12=V13V23ρhme120 GPa7.8 GPa5.4 GPa0.280.451530 kg/m30.2 mm0.052


### 3.2. Experimental Setup

The LVI experiments were carried out on the INSTRON CEAST 9350 drop weight impact testing machine, as shown in Figure5(a). A square steel fixture with a circular hole of radius 50 mm was customized, as shown in Figure 5(c). The laminate specimen was fixed on the steel fixture by toggle clamps pressing the corner, as shown in Figure 5(b). It has been proven through natural frequency analysis that the boundary condition approximate is simply supported. In order to get a relatively high initial impact velocity with lower impact energy, the lightweight tip holder matching the testing machine was chosen. The total mass of the drop weight, including the hemispherical tip with a diameter of 16 mm, is 2.277 kg.Figure 5
Drop weight impact testing machine and customized steel fixture (a) The INSTRON CEAST 9350, (b) experimental state of the specimen during the test, (c) steel fixture with a coaxial circular hole.To ensure that no delamination was occurring during the validation tests, the DTL of the laminate was first measured according to the ASTM D7136/D7136M-12 [45]. The LVI test was conducted with energy of 10 J. The contact force history curve of the test is shown in Figure 6, where the DTL of 3226 N was identified at the point where the sudden drop first occurred. The LVI tests with three different levels of energy were conducted, and the corresponding initial impact velocities are listed in Table 3. When the impact energy was increased to about 5 J, it could be seen from the C-scan images that delamination occurred.Figure 6
DTL identification from contact force histories curve.Table 3
Impact energy and corresponding initial velocity of LVI tests.
Impact energy2.14 J3.14 J4.24 JInitial velocity1.37 m/s1.66 m/s1.93 m/s


## 3.1. Specimen

The square FRP laminate specimen with 125 mm on each side employed in the following experiments and simulations was fabricated by Guangwei Composites Material Co., Ltd. using the unidirectional carbon fiber prepreg USN20000 with a thickness of 0.2 mm. The ply configuration of the laminate was [(0/90)4]s. Table 2 gives the mechanical properties of the lamina. me was obtained by fitting test data with different impact velocities following the Y-C logarithmic function in Equation (18).Table 2
Lamina properties of the laminate specimen.
E11E22=E33G12=G13V12=V13V23ρhme120 GPa7.8 GPa5.4 GPa0.280.451530 kg/m30.2 mm0.052


## 3.2. Experimental Setup

The LVI experiments were carried out on the INSTRON CEAST 9350 drop weight impact testing machine, as shown in Figure5(a). A square steel fixture with a circular hole of radius 50 mm was customized, as shown in Figure 5(c). The laminate specimen was fixed on the steel fixture by toggle clamps pressing the corner, as shown in Figure 5(b). It has been proven through natural frequency analysis that the boundary condition approximate is simply supported. In order to get a relatively high initial impact velocity with lower impact energy, the lightweight tip holder matching the testing machine was chosen. The total mass of the drop weight, including the hemispherical tip with a diameter of 16 mm, is 2.277 kg.Figure 5
Drop weight impact testing machine and customized steel fixture (a) The INSTRON CEAST 9350, (b) experimental state of the specimen during the test, (c) steel fixture with a coaxial circular hole.To ensure that no delamination was occurring during the validation tests, the DTL of the laminate was first measured according to the ASTM D7136/D7136M-12 [45]. The LVI test was conducted with energy of 10 J. The contact force history curve of the test is shown in Figure 6, where the DTL of 3226 N was identified at the point where the sudden drop first occurred. The LVI tests with three different levels of energy were conducted, and the corresponding initial impact velocities are listed in Table 3. When the impact energy was increased to about 5 J, it could be seen from the C-scan images that delamination occurred.Figure 6
DTL identification from contact force histories curve.Table 3
Impact energy and corresponding initial velocity of LVI tests.
Impact energy2.14 J3.14 J4.24 JInitial velocity1.37 m/s1.66 m/s1.93 m/s


## 4. Finite Element Model

To further verify the effectiveness of the proposed SRD model, a VUMAT-based FEM simulation was also conducted using the commercial software ABAQUS/Explicit as a comparative assessment [46, 47]. The finite element model of the circular laminate impacted by a hemispherical rigid impactor at the center was developed under the boundary condition of EMSS, as shown in Figure 7. The lamina was modeled using C3D8R elements, which means 8-node cubic elements with reduced integration. To save computing time, various mesh sizes in different regions of the laminate were introduced in the finite element model. The mesh grid became sparser as the distance from the laminate center increased, from 0.2 mm of the central contact area to 1 mm of the boundary area. The general contact algorithm with a penalty enforcement formulation of friction was applied to the tangential interaction between the impactor and the laminate plate. The friction coefficient was set to 0.3 based on several experimental summaries [48].Figure 7
Finite element model of impact on laminate plate.The SRD updating of mechanical properties was implemented through the secondary development of the user material subroutine VUMAT. The flowchart of this secondary development subroutine of FEM is shown in Figure8. The strain rate at each grid point of the laminate was calculated using the invoked strain increment (“strainInc”) at every time step (“stepTime”). Then, the material properties (“props”) matrix was modified and reassigned correspondingly for the next time step. The strain components for the next time step (“stateNew”) were updated by adding the strain increment to the previous old strain components (“stateOld”). The same update operation was applied to the stress components (“stressNew”) by accumulating the product of the modified stiffness matrix and new strain components to the previous old stress components (“stressOld”). The updated components of this time step were taken as the old for the next time step. The entire contact history of stress and strain distribution was obtained by repeating these operations of the cycle step-by-step.Figure 8
Flowchart of the VUMAT-based FEM simulations.


## 5. Results and Discussion

### 5.1. Convergence Analysis

There are two integral increments in the proposed SRD model: the circumferential stepΔφ in Equation (2) and the time step Δt in Equation (17). For the balance of prediction accuracy and time efficiency, the convergence analysis was conducted to determine the appropriate value of these two uncertain parameters.First, the averaged homogenized modulus of the specimen above was obtained with different values ofΔφ according to Equation (2). The results are shown in Figure 9. It can be seen that the value of Ee converges with Δφ decreases. When Δφ is less than 15°, the difference between Ee and its convergence value is negligible. Hence, Δφ=10° was adopted in the subsequent numerical examples.Figure 9
Averaged homogenized modulus with differentΔφ.The same procedure was taken for the time stepΔt. The model for the circular FRP laminate with a diameter of 100 mm impacted by a hemispherical rigid impactor at the center was developed, as shown in Figure 1. The impactor had a diameter of 16 mm and a mass of 3 kg. The contact force history at V0=1m/s with various Δt was obtained by recursive integration of Equation (14). The results under EIC boundary are shown in Figure 10, and the corresponding time consumption on Intel® Core™ i7-9750H, with 12 M cache and frequency up to 4.5 GHz, is listed in Table 4. It can be found that variation of the time step Δt has no effect on the SRI contact force history results. As Δt decreases, the impact period becomes shorter and the extreme value of contact force increases. Figure 10 and Table 4 show that decrease of Δt from 1e−6 s to 5e−7 s only has a slight effect on the SRD history but cause a dramatic increase of time consumption. Hence, after deliberation, Δt=1e−6 s was selected for model calculation in the following section.Figure 10
Contact force histories with variousΔt.Table 4
Time consumption of differentΔt.
Δt (s)5e−71e−65e−61e−5Time consumption (s)6881.473021.03633.38303.46


### 5.2. Model Validation

Figures11–13 give the comparison of contact force history curves obtained by the SRD/SRI model with the experiment and FEM simulation under three different levels of impact energy. The FEM curves are obtained by filtering the original data. Since the elastic modulus of carbon fiber has been proven to be insensitive to the strain rate, fiber-dominated mechanical properties (E11) of the lamina do not vary with the strain rate during the model calculation and FEM.Figure 11
Comparison of contact force history from model calculation with FEM simulation and experimental results (2.14 J).Figure 12
Comparison of contact force history from model calculation with FEM simulation and experimental results (3.14 J).Figure 13
Comparison of contact force history from model calculation with FEM simulation and experimental results (4.24 J).It can be found that the contact force histories from the spring-mass model agree well with the FEM results in both SRD and SRI cases. But there are some discrepancies in the impact period and the extreme value of contact force between the SRI model and experimental results. These discrepancies have been decrease greatly by considering the strain rate effect through the proposed SRD model. The contact force history curves from the SRD spring-mass model are consistent with the experimental results well. The SRD FEM results also coincides with that from the experiment. It shows that the strain rate effect on “large mass” impact response has been accurately evaluated and the validity of the proposed SRD model has been verified, also with the VUMAT-based FEM. The time consumption of the proposed SRD model was only about one hour, showing a high efficiency compared with the FEM which takes more than ten hours.It is also worth noting in Figure13 that when the impact energy is 4.24 J, the experimental contact force is slightly lower than the SRD model results near the peak of the history curve, due to that the extreme value of the impact force is approaching to DTL of the specimen with the increase of the impact energy. Under this impact energy, even if delamination damage does not occur, local crushing of the resin matrix in the contact area would result in a slight decrease of contact force extreme value.


### 5.3. Parametric Analysis

The experimental verification described in the previous section was conducted on the laminates reinforced by carbon fiber, which has previously been reported to be SRI. The performance of the laminate which is dominated by the carbon fiber is virtually unaffected by the strain rate variation. However, for laminates whose reinforcing fiber is strain-rate sensitive material, the fiber-dominated properties of the lamina vary with the strain rate during the contact process, which affects the impact response of the laminate to a greater extent compared with the laminates with strain-rate insensitive fibers.To further investigate the reliability and practicability of the proposed SRD model, parametric analysis was performed to test its capacities for accurately evaluating the strain rate effect on the LVI response of laminates reinforced by fiber with different strain rate dependency. The contact histories calculated from the proposed SRD model were compared with the results from the VUMAT-based FEM simulation in three different cases: SRI and SRD excludeE11 and SRD include E11.Parameters of the model calculation and FEM simulation were the same as in the above experiments except for the strain rate dependency of fiber-dominated mechanical properties (E11) of the lamina. The history curves of contact force, impactor displacement, and velocity with impact energy of 2.12 J and 3.14 J are shown in Figures 14 and 15. The SRI curves were presented as the basis for comparison. The latter two cases involved modifying the constitutive relations of the lamina with and without the tensile modulus E11, corresponding to the laminate composed of strain-rate-sensitive (such as glass and Kevlar) and strain-rate-insensitive (such as PVA and carbon) fibers, respectively.Figure 14
Contact histories with impact energy of 2.12 J: (a) contact force, (b) displacement, and (c) velocity of the impactor.
(a)(b)(c)Figure 15
Contact histories with impact energy of 3.14 J: (a) contact force, (b) displacement, and (c) velocity of the impactor.
(a)(b)(c)Comparing the curves of two SRD cases with the SDI results, the strain rate effect on impact response is more pronounced for the laminate with strain-rate-sensitive fibers, which is attributed to the dominant role of the fiber on the macroscopic mechanical properties of the laminate. For the laminate with strain-rate-insensitive fibers, the strain rate effect on the impact response, which was mainly controlled by the viscoplasticity of the matrix, was less noticeable than in the case of strain-rate sensitive fibers.Good consistency between the proposed SRD model and FEM results was observed in the variation of period and amplitude of the curves. The varying trends in SRI an SRD curves proved the validity and reliability of the modified SRD spring-mass model in predicting the “large mass” impact response of the laminate. This conclusion was also supported by the displacement and velocity curves of the impactor in Figures14(b) and 14(c) and Figures 15(b) and 15(c) where the model calculation curves showed similar trends to the FEM curves in each case. The applicability of the proposed model to the simulation of key impact quantities shows its relatively higher value.


## 5.1. Convergence Analysis

There are two integral increments in the proposed SRD model: the circumferential stepΔφ in Equation (2) and the time step Δt in Equation (17). For the balance of prediction accuracy and time efficiency, the convergence analysis was conducted to determine the appropriate value of these two uncertain parameters.First, the averaged homogenized modulus of the specimen above was obtained with different values ofΔφ according to Equation (2). The results are shown in Figure 9. It can be seen that the value of Ee converges with Δφ decreases. When Δφ is less than 15°, the difference between Ee and its convergence value is negligible. Hence, Δφ=10° was adopted in the subsequent numerical examples.Figure 9
Averaged homogenized modulus with differentΔφ.The same procedure was taken for the time stepΔt. The model for the circular FRP laminate with a diameter of 100 mm impacted by a hemispherical rigid impactor at the center was developed, as shown in Figure 1. The impactor had a diameter of 16 mm and a mass of 3 kg. The contact force history at V0=1m/s with various Δt was obtained by recursive integration of Equation (14). The results under EIC boundary are shown in Figure 10, and the corresponding time consumption on Intel® Core™ i7-9750H, with 12 M cache and frequency up to 4.5 GHz, is listed in Table 4. It can be found that variation of the time step Δt has no effect on the SRI contact force history results. As Δt decreases, the impact period becomes shorter and the extreme value of contact force increases. Figure 10 and Table 4 show that decrease of Δt from 1e−6 s to 5e−7 s only has a slight effect on the SRD history but cause a dramatic increase of time consumption. Hence, after deliberation, Δt=1e−6 s was selected for model calculation in the following section.Figure 10
Contact force histories with variousΔt.Table 4
Time consumption of differentΔt.
Δt (s)5e−71e−65e−61e−5Time consumption (s)6881.473021.03633.38303.46


## 5.2. Model Validation

Figures11–13 give the comparison of contact force history curves obtained by the SRD/SRI model with the experiment and FEM simulation under three different levels of impact energy. The FEM curves are obtained by filtering the original data. Since the elastic modulus of carbon fiber has been proven to be insensitive to the strain rate, fiber-dominated mechanical properties (E11) of the lamina do not vary with the strain rate during the model calculation and FEM.Figure 11
Comparison of contact force history from model calculation with FEM simulation and experimental results (2.14 J).Figure 12
Comparison of contact force history from model calculation with FEM simulation and experimental results (3.14 J).Figure 13
Comparison of contact force history from model calculation with FEM simulation and experimental results (4.24 J).It can be found that the contact force histories from the spring-mass model agree well with the FEM results in both SRD and SRI cases. But there are some discrepancies in the impact period and the extreme value of contact force between the SRI model and experimental results. These discrepancies have been decrease greatly by considering the strain rate effect through the proposed SRD model. The contact force history curves from the SRD spring-mass model are consistent with the experimental results well. The SRD FEM results also coincides with that from the experiment. It shows that the strain rate effect on “large mass” impact response has been accurately evaluated and the validity of the proposed SRD model has been verified, also with the VUMAT-based FEM. The time consumption of the proposed SRD model was only about one hour, showing a high efficiency compared with the FEM which takes more than ten hours.It is also worth noting in Figure13 that when the impact energy is 4.24 J, the experimental contact force is slightly lower than the SRD model results near the peak of the history curve, due to that the extreme value of the impact force is approaching to DTL of the specimen with the increase of the impact energy. Under this impact energy, even if delamination damage does not occur, local crushing of the resin matrix in the contact area would result in a slight decrease of contact force extreme value.


## 5.3. Parametric Analysis

The experimental verification described in the previous section was conducted on the laminates reinforced by carbon fiber, which has previously been reported to be SRI. The performance of the laminate which is dominated by the carbon fiber is virtually unaffected by the strain rate variation. However, for laminates whose reinforcing fiber is strain-rate sensitive material, the fiber-dominated properties of the lamina vary with the strain rate during the contact process, which affects the impact response of the laminate to a greater extent compared with the laminates with strain-rate insensitive fibers.To further investigate the reliability and practicability of the proposed SRD model, parametric analysis was performed to test its capacities for accurately evaluating the strain rate effect on the LVI response of laminates reinforced by fiber with different strain rate dependency. The contact histories calculated from the proposed SRD model were compared with the results from the VUMAT-based FEM simulation in three different cases: SRI and SRD excludeE11 and SRD include E11.Parameters of the model calculation and FEM simulation were the same as in the above experiments except for the strain rate dependency of fiber-dominated mechanical properties (E11) of the lamina. The history curves of contact force, impactor displacement, and velocity with impact energy of 2.12 J and 3.14 J are shown in Figures 14 and 15. The SRI curves were presented as the basis for comparison. The latter two cases involved modifying the constitutive relations of the lamina with and without the tensile modulus E11, corresponding to the laminate composed of strain-rate-sensitive (such as glass and Kevlar) and strain-rate-insensitive (such as PVA and carbon) fibers, respectively.Figure 14
Contact histories with impact energy of 2.12 J: (a) contact force, (b) displacement, and (c) velocity of the impactor.
(a)(b)(c)Figure 15
Contact histories with impact energy of 3.14 J: (a) contact force, (b) displacement, and (c) velocity of the impactor.
(a)(b)(c)Comparing the curves of two SRD cases with the SDI results, the strain rate effect on impact response is more pronounced for the laminate with strain-rate-sensitive fibers, which is attributed to the dominant role of the fiber on the macroscopic mechanical properties of the laminate. For the laminate with strain-rate-insensitive fibers, the strain rate effect on the impact response, which was mainly controlled by the viscoplasticity of the matrix, was less noticeable than in the case of strain-rate sensitive fibers.Good consistency between the proposed SRD model and FEM results was observed in the variation of period and amplitude of the curves. The varying trends in SRI an SRD curves proved the validity and reliability of the modified SRD spring-mass model in predicting the “large mass” impact response of the laminate. This conclusion was also supported by the displacement and velocity curves of the impactor in Figures14(b) and 14(c) and Figures 15(b) and 15(c) where the model calculation curves showed similar trends to the FEM curves in each case. The applicability of the proposed model to the simulation of key impact quantities shows its relatively higher value.


## 6. Conclusions

This study presented a modified SRD spring-mass model to more accurately predict the contact history of an orthotropic FRP laminate subjected to out-of-plane LVI. Impact response prediction of FRP laminates considering the strain rate effect was theoretically established for the first time, although it has been previously realized by numerical FEM simulations. The constitutive relation of the laminate needs to be constantly updated with the strain rate variation in space and time dimensions simultaneously. The contact history and structural response of the FRP laminate can be obtained by solving the modified SDOF spring-mass model using a step-by-step recursive integration. The drop-weight tests and VUMAT-based FEM simulations were conducted to validate the model. The results show that the strain rate effect reflected in the contact force during “large mass” impact can be accurately predicted. Also, the proposed SRD model was proven to be more efficient than FEM simulations, as better results were achieved with less time consumption. Further, parametric analysis of strain rate effect evaluation on the LVI response of laminates reinforced by fiber with different strain rate dependency was performed. The good consistency between the analytical solution and FEM results before and after containing different strain rate dependencies showed the reliability and practicability of the proposed SRD model in evaluating various degrees of the strain rate effect on the LVI response. For future studies, some developed damage models can be combined with this modified SRD spring-mass model to further analyze the delamination evolution processes of FRP laminates under LVI more accurately.



---

*Source: 2121209-2023-03-02.xml*