Rate-Dependent Viscoelastic–Damage Constitutive Modeling and Numerical Simulation of Engineered Cementitious Composites Under Dynamic Compression

Concrete, the most extensively used material in building structures, is subjected not only to normal service loads but also to severe dynamic loads like impacts and explosions during a structure's service life [1]. In military protective structures, these loads can rise to 100 MPa within 1 ms, demanding material data at strain rates of 10$^2$–10$^3$ s$^{-1}$ [2]. Under explosion loads, three main types of local damage can occur in concrete structures [3], [4], [5]. As for crater damage, compressive forces crush the concrete on the blast-exposed face, creating a crater-like depression. In terms of spall damage, when the compressive wave from the explosion reaches the back face, it reflects into a tensile wave. If the resulting tensile stress exceeds the concrete's tensile strength, spalling occurs. In high-energy explosions, the spalled fragments can be ejected at velocities capable of causing injury. In terms of breach damage, if the shock wave from the explosion is strong enough to overcome the structure's resistance, a breach occurs. This results in high-energy fragment ejection and potential injury. To boost concrete's blast-resistant capacity, it's key to enhance its compressive and tensile strength and toughness. The compressive properties are crucial for resisting penetration and perforation in blast-resistant applications [6]. They determine whether a material can withstand impact or explosive loads without failure.

Existing studies have shown that fiber-reinforced concrete, such as Engineered cementitious composites (ECCs), can significantly improve these properties. ECCs are a class of ultra-ductile, fiber-reinforced cementitious composites that exhibit strain-hardening behavior and multiple micro-cracking under tension, resulting in a tensile strain capacity of 3–7%—two orders of magnitude higher than ordinary concrete [7]. Unlike conventional fiber-reinforced concrete, which typically shows tension-softening after first cracking, ECCs maintain load-carrying capacity through closely spaced micro-cracks ($<$100 $\mu m$), providing superior energy absorption and damage. Currently, ECCs are applied to seismic-resistant structures [8], shotcrete reinforcement [9], etc., and have shown good engineering performance [10].

Kai et al. [11] conducted dynamic mechanical experiments on ECCs. It was found that when the strain rate increased from 100 s$^{-1}$ to 190 s$^{-1}$, the compressive strength of ECCs doubled, and the peak strain increased from 0.55% to 0.81%. Shi et al. [12] conducted dynamic compression tests on polypropylene-ECC specimens with different polypropylene fiber volume fractions using a Split Hopkinson Pressure Bar (SHPB) device. The results indicated that the specimen with a 1.5% fiber volume fraction exhibited the highest dynamic compressive strength across all strain rates. However, the specimen with a 1.8% fiber volume fraction demonstrated a lower dynamic compressive strength than the matrix material, indicating that increasing the fiber volume fraction does not always enhance the material's dynamic compressive strength. Liu et al. [13] conducted SHPB tests on polyethylene-ECCs. It was found that as the strain rate increases, the material's peak strength and elastic modulus both rise, suggesting that the excellent energy absorption of polyethylene-ECCs comes from its strain hardening and the long post-peak descending part in the stress-strain curve.

Xu et al. [14] calibrated the Holmquist-Johnson-Cook (HJC) constitutive model based on the SHPB dynamic compression test data of polyvinyl alcohol-ECCs, identifying 21 HJC parameters and validating the simulation's accuracy. By analyzing dynamic compression stress-strain curves at five strain rates, the strain rate effect on the dynamic enhancement factor of peak stress was examined. Numerical simulations using Livermore Software for DYnamic Analysis (LS-DYNA) software explored the relationship between the material's failure process, failure mode, and strain rate. Li et al. [15] also demonstrated via numerical simulation and experimental studies that ECC plates have excellent explosion-resistant protective properties. The dynamic compressive strength of the material significantly influences the blast-resistant performance of the structure.

Although previous studies have examined the impact resistance and energy absorption of ECCs, their rate-dependent compressive response has not yet been adequately represented within a unified constitutive framework. In particular, the concurrent effects of nonlinear viscoelastic deformation and progressive material damage on the ascending, peak, and post-peak stages of the dynamic stress–strain response require further investigation. This study therefore examined the dynamic compressive behavior of ECC containing 2.0% polyvinyl alcohol fibers through SHPB tests conducted at average strain rates ranging from approximately 15 to 200 s$^{-1}$. A rate-dependent constitutive model combining the Zhu–Wang–Tang nonlinear viscoelastic formulation with a Weibull statistical damage function was developed and calibrated against the experimental results. The model was subsequently implemented in LS-DYNA to reproduce the stress–strain response and failure characteristics observed at different strain-rate levels. The study aims to clarify the interaction between rate-dependent deformation and damage evolution in ECC and to provide a constitutive basis for the numerical analysis of ECC components subjected to impact-type loading.

This study adopts the ECC matrix mix proportion proposed by Yang and Li [16] (Table 1) because it has been extensively validated under both static and low-rate dynamic loading, ensuring reproducibility and consistency for high-rate SHPB tests. It was made using ordinary Portland cement (P.O. 42.5), Class I fly ash (low-calcium), extra-fine river sand (average grain size $<$ 1.3 mm), and a high-efficiency water reducer (polycarboxylate). The fiber used was polyvinyl alcohol fiber. The physical parameters of the polyvinyl alcohol fiber are presented in Table 2. Specimens were made in polyvinyl chloride molds (outer diameter 50 mm, height 30 mm). After 28 days of standard curing, the specimen ends were ground on a flat-bed grinder. The ground specimens had a diameter of 46 mm and a height of 25 mm. Specimen fabrication and grinding are shown in Figure 1.

This test employed a $ \Phi$50 mm conical variable-cross-section SHPB device with a steel striker rod (density 7850 kg/m$^3$, wave velocity 5172 m/s, and elastic modulus 210 GPa). Data were collected using a DH5960 multi-channel ultra-dynamic signal parallel synchronous testing system and processed via the two-wave method [17]. A total of 16 specimens was prepared, including 12 specimens used for the formal SHPB tests and four spare specimens. The 12 test specimens were numbered C1–C12 and divided into four loading groups, S1–S4, with three specimens in each group. Figure 2 shows the schematic diagram of the SHPB apparatus. Prior to formal testing, preliminary tests were conducted to determine the impact velocities corresponding to four strain rates. The lowest strain rate (S1) corresponded to the impact velocity at which visible cracks initiated in all specimens, while the highest strain rate (S4) represented the impact velocity reaching the load-bearing capacity limit of the pressure bar. The intermediate strain rates S2 and S3 were then established by approximately trisecting this range. Specimen grouping configurations are detailed in Table 3.

The SHPB dynamic impact compression test results for ECC specimens are summarized in Table 4. It is evident that the average strain rate of the specimens increases with the impact velocity, and there are significant differences in strain rates across different loading rates. This indicates the successful implementation of the planned variable strain rate loading. Under these conditions, the peak stress and peak strain of the ECC specimens all increase with strain rate, showing a clear strain rate effect.

To compare the impact of varying strain rates on mechanical properties, this study uses the mechanical property values obtained at the S1 impact level as a reference and analyzes the relative changes at S2–S4 levels. Compared to S1, the peak stress and peak strain at S2–S4 levels increase by 19.7%–41.0% and 28.4%–98.6%, respectively. Peak strain is the most strain rate-sensitive parameter. As shown in Figure 3, the relationship between property changes and strain rate follows an initially rapid and subsequently slower pattern. This suggests that the mechanical properties of ECCs are more sensitive to strain rate variations in the low-to-moderate range. From the relative trends of the curves, peak strain is most affected by strain rate (uppermost curve with the steepest slope), and then peak stress (lowermost curve with the gentlest slope). The peak stress curve also shows a gradual transition from steep to gentle with increasing strain rate, indicating a diminishing influence of strain rate on peak stress beyond a certain threshold.

For the constitutive relations of concrete-like materials under impact loading, nonlinear viscoelastic constitutive models are commonly employed. This study adopts the representative Zhu–Wang–Tang (ZWT) nonlinear viscoelastic model [18], with its constitutive equation given by:

where, $\sigma$, $\varepsilon$, and $t$ represent stress, strain, and time, respectively; $E_0$, $\alpha$, and $\beta$ denote elastic coefficients; $E_1$ and $\theta_1$ correspond to the equivalent modulus and relaxation time for viscoelastic response at low strain rates, while $E_2$ and $\theta_2$ represent those for high strain rates. Based on the ZWT constitutive model [19], Umar et al. [20] posited that, under impact loading, material elements lack sufficient time for low-frequency relaxation processes. Consequently, the low-frequency response can be neglected, and the ZWT constitutive model can be simplified as follows:

where, $E=E_0+E_1$ represents the equivalent comprehensive elastic modulus.

To more accurately characterize the influence of internal material degradation on constitutive behavior, this study integrates a damage variable D based on the Weibull distribution into the constitutive framework to describe damage evolution, as given by Eq. (3) [20]:

where, $\varepsilon_c$ denotes the elastic strain limit of the cementitious material, conventionally taken as $\varepsilon_c=0.7 \varepsilon_p$, where $\varepsilon_p$ is the peak strain. $\lambda$ (scale parameter) and $m$ (shape parameter) are Weibull distribution parameters characterizing damage evolution behavior:

Based on the stress-strain data obtained from SHPB impact tests, the constitutive relationship was calibrated using a FORTRAN computational program. The calibrated parameters for each strain-rate level are listed in Table 5, and the model predictions are compared with the experimental stress–strain curves in Figure 4. The constitutive model reproduced the ascending branch and peak region of the experimental stress–strain curves with reasonable accuracy, whereas larger discrepancies were observed in the post-peak softening stage, particularly at the intermediate and high strain-rate levels. The blue “Segment 1” (indicating the linear ascending phase) experimental data aligns well with the fitting curve. This indicates that during the elastic deformation stage, the material's impact response can be accurately characterized by the parameters $E$, $E_2$, and $\theta_2$. In this phase, the S1-S4 data was fitted to $E \approx 43.2 \mathrm{GPa}, E_2 \approx 13.5 \mathrm{GPa}$, and $\theta_2 \approx 1.0 \mathrm{~S}$. Here, $E$ and $E_2$ are modulus parameters with units of GPa, whereas $\theta_2$ is the relaxation time of the high-frequency viscoelastic element and therefore has a unit of seconds. The nearly constant fitted values of these parameters across the investigated loading levels suggest that no pronounced systematic variation in the identified elastic and short-term viscoelastic parameters was observed within the present dataset. The peak stress increases with the strain rate, and the peak strain slightly shifts to the right. This is consistent with the positive gain of the rate-amplification term $E_2 \theta_2 \dot{\varepsilon}\left[ 1-\exp \left(-\frac{\varepsilon}{\dot{\varepsilon} \theta_2}\right)\right]$ near $\varepsilon_c$ in Eq. (4).

The experimental data and fitting data for the yellow “Segment 2” (representing the nonlinear rising and softening stage) both exhibit the characteristics of sustained residual bearing capacity and slow decay in the post-peak stage from S2 to S4. Compared with the experimental data, the fitting data shows a faster decay rate in the softening stage. This is because the single exponential softening form $\exp \left[-\frac{\left(\varepsilon-\varepsilon_c\right)^m}{\lambda}\right]$ in the constitutive model underestimates the material properties in the softening stage. On the parameter level, the softening shape index is set as $m=2.72$ at S1, and tends to $m \approx 2.0$ at medium and high strain rates. This indicates that as the strain rate increases, the internal cracks of the material penetrate more rapidly, the damage localization is more severe, and the failure mode becomes more dominated by brittleness. The softening scale factor $\lambda$ spans approximately two orders of magnitude with the increase in strain rate (from $1.18 \times 10^{-6}$ to approximately $10^{-4}$). This reflects that as the strain rate increases, the threshold for the initiation of ECC damage is raised, and the material's resistance to damage in the softening stage is enhanced.

In terms of the accuracy of fitting parameters, the root mean square error and standard error of each parameter increase with the rise in strain rate. Specifically, the overall fitting root mean square error increased from 0.98 at S1 to 8.19 at S4. Meanwhile, the standard error of $\lambda$ increased substantially from $5.30 \times 10^{-8}$ to $2.51 \times 10^{-4}$, indicating greater uncertainty in the identification of $\lambda$ at higher loading-rate levels. The sources of error can mainly be attributed to three categories. The first category comprises factors related to the SHPB system and the experimental process, such as dynamic imbalance during the initial loading phase, as well as high-frequency noise caused by incomplete alignment and unevenness of the specimen end faces. The second category encompasses factors associated with data processing and identification, such as approximating the time-varying strain rate with a single $\dot{\varepsilon}$, where the determination of $\varepsilon_{\mathrm{c}}$ is influenced by smoothing and sampling, as well as parameter correlations resulting in an ill-posed problem. The third category pertains to the simplification in the model mechanism. Specifically, in Eq. (4), the post-peak behavior is characterized solely by a single exponential degradation form, which is incapable of accurately representing the complex softening characteristics associated with fragmentation, fiber bridging, and densification. As a result, this leads to an underestimation of the softening and residual load-bearing capacity within the S2-S4 range. Within the investigated strain-rate range, the model captured the pre-peak stiffness and the rate-dependent variation in peak response using a limited number of physically interpretable parameters. Its predictive accuracy decreased in the post-peak region because the single exponential damage function could not fully represent the combined effects of matrix fragmentation, fiber bridging, and material compaction. The model should therefore be used with caution when the post-peak response or residual load-carrying capacity is of primary concern.

Numerical simulation, as a complement to direct experimentation, is crucial for capturing elusive details and phenomena in tests. To refine the experimental outcomes, this study employs LS-DYNA software to conduct a comprehensive analysis of the SHPB dynamic compression test.

The choice of material model and parameters is crucial for the accuracy of numerical analysis. In this numerical simulation, an elastic model, defined by the keyword “MAT_ELASTIC” and based on the actual parameters of the striker rod, was selected for the striker rod material. Due to the unique properties of ECCs, conventional concrete models are inadequate for simulating its behavior. Consequently, an LS-DYNA subroutine was developed using FORTRAN to define a novel ECC material model. The subroutine flow chart is illustrated in Figure 5. This model, featuring three independent strength surfaces, effectively captures the hardening and softening processes of ECCs.

The subroutine execution logic and the corresponding strength surfaces are depicted in Figure 5. To provide a clear description of the computational pipeline, the key variables and parameters featured in the flowchart are defined as follows: $P$ represents the hydrostatic pressure ($P<0$ denotes dynamic tension, whereas $P \geq 0$ denotes compression). $\sigma_{\text {trial }}$ is the trial stress component evaluated before checking the yield criterion, and $J_2$ is the second invariant of the stress deviator tensor. $DIF_t$ and $DIF_c$ denote the Dynamic Increase Factors for dynamic tensile strength and compressive strength, respectively, accounting for rate-dependent strain-rate enhancement. $D_t$ signifies the dynamic tensile damage variable; $D_s$ and $D_h$ represent the shear damage variable and hydrostatic crushing damage variable under dynamic compression, respectively. $D$ denotes the overall scalar damage variable updated at each time step. $\sigma_m$, $\sigma_y$, and $\sigma_r$ correspond to the maximum strength surface, yield strength surface, and residual strength surface of the ECCs, respectively. $\sigma$ and $\varepsilon$ represent the updated stress tensor and total strain tensor transferred to the LS-DYNA main solver at the end of the load step.

Eq. (5) is employed to characterize the strain-rate effect of ECCs [21]. $W_x$ is the strain-rate threshold coefficient that governs the onset of enhancement and is fitted to 0.8 from test data; $W_y$ is the strength saturation coefficient (capping the dynamic increase factor) taken as 1.05 from the study by Gebbeken and Ruppert [21]; $F_m$ is the static reference strength based on the quasi-static test results.

In LS-DYNA, the equation of state defined by the keyword *EOS_TABULATED is used to describe the relationship between volumetric strain and hydrostatic pressure for concrete (Eq. (6)). In the equation, $C(\mu)$ denotes the pressure coefficient, $\theta(\mu)$ represents the temperature parameter, $\gamma_0$ designates the specific heat capacity, and $E_0$ corresponds to the volumetric energy density.

Under uniaxial dynamic compression, the uniaxial stress-strain curve is projected into a “hydrostatic pressure - volumetric strain” relationship using Eqs. (7) and (8).

where, $v_{eff}$ is the evolving effective Poisson's ratio, calculated using the following empirical formula [22]:

where, $\eta=\sigma / f_c$, with $f_c$ and $\varepsilon_c$ as the dynamic peak stress and the corresponding peak strain measured in the experiment, respectively; $\eta_1$ is taken as 0.75, $\eta_2$ as 1.0 , and $v_0$ as 0.2 [23]. When $\varepsilon>\varepsilon_c$, it indicates the material has entered the post-peak strain-softening phase, meaning ECCs exhibit volumetric expansion after reaching the dynamic peak stress. In this phase, $v_{eff}$ is jointly defined by $\alpha_d=0.02-0.08$ (taken as 0.05 in this study) and the ultimate strain $\varepsilon_d$.

To address the inconsistency between numerical simulations and experiments caused by bullet impact on the incident rod failing to generate a triangular incident wave, this study directly applies the incident wave at the rod end, aligning the loading conditions in both experiments and simulations. Additionally, penalty function erosion contact is used between specimens to prevent rod penetration, and mesh refinement is applied at contact points to avoid mesh distortion during calculations, as shown in Figure 6.

The total simulation time was set to 1000 $\mu s$ with a time step of 1 $\mu s$. Considering the rod length and wave speed in the specimen, it takes about 38 $\mu s$ for the stress wave to reflect three times within the specimen. The load rise time was set to 120 $\mu s$ to ensure stress uniformity. Following the research by Shuai et al. [24], failure was defined by maximum principal strain.

Figure 7 shows the stress wave cloud diagram in the rod under the S1 strain rate. Figure 7a illustrates the stress state induced by the incident wave propagating through the incident bar at $t$ = 278 $\mu s$. The maximum stress of 0.847 MPa, corresponding to the red contour region, is observed at the wavefront location. Symmetrical stress variations, represented by the green contour regions, occur on both sides of the peak stress zone due to the passage of the stress transient (either post-wavefront or pre-wavefront). This stress distribution pattern aligns with characteristic elastic wave propagation behavior. Figure 7b depicts the stress state at $t$ = 423 $\mu s$ as the incident wave propagates into the specimen. At this stage, the specimen exhibits a relatively uniform von Mises stress of approximately 0.614 MPa. Crucially, the stress distributions near the bar-specimen interfaces at both ends closely match, satisfying the fundamental premise of dynamic stress equilibrium required by one-dimensional stress wave theory for valid SHPB analysis. Figure 7c depicts the stress state at $t$ = 533 $\mu s$, where the impedance mismatch at the ECC specimen interface causes partial reflection of the incident wave back into the incident bar (left side), generating a reflected wave with a peak stress of approximately 0.335 MPa. Concurrently, the transmitted wave propagates into the output bar (right side), exhibiting a peak stress of about 0.534 MPa. The significant contrast in magnitude between the reflected wave (approximately 0.335 MPa) and the transmitted wave (approximately 0.534 MPa) directly characterizes the dynamic stress-strain response and energy dissipation mechanisms of the ECC material under high-strain-rate loading. The sequential visualization of stress wave propagation, clear demonstration of dynamic stress equilibrium within the specimen, and physically consistent reflection/transmission ratio collectively validate the reliability of the present numerical model in replicating the fundamental mechanical behavior of ECC materials under SHPB testing conditions.

Figure 8 illustrates the cyclic reflection of stress waves between the specimen's end faces. It takes 533 $\mu s$ for stress waves to traverse the specimen and reach the transmission bar strain gauge. By this time, the specimen undergoes four distinct stages of stress wave evolution internally. In Stage 1 (initial wave entry), the incident compressive stress wave propagates from the incident bar into the specimen. A steep stress gradient forms near the input interface as the wavefront advances into the material. In Stage 2 (first reflection & wave superposition), upon reaching the opposite (output bar) interface, the initial wave is partially reflected back into the specimen due to the impedance mismatch. This reflected tensile/compressive wave (depending on the impedance relationship) travels back towards the input interface. Simultaneously, the initial wave continues its propagation into the output bar (initiating transmission). Crucially, the incident wave still entering the specimen overlaps (superimposes) with this first reflected wave within the specimen volume. In Stage 3 (multiple reflections and homogenization), the wave reflected from the output interface reaches the input interface, where it is partially reflected back into the specimen again. This process leads to multiple internal reflections of the stress wave between the two end faces. The continual superposition of successive incident, reflected, and re-reflected wave components progressively distributes the stress state throughout the specimen volume. In Stage 4 (dynamic stress equilibrium), after sufficient cycles of internal reflection and superposition (typically within three wave transits, depending on specimen length and wave speed), the stress field within the specimen becomes spatially uniform. The amplitude of the oscillating superimposed waves attenuates towards a constant mean value. This state of dynamic stress equilibrium is achieved when the stresses at both specimen-bar interfaces are equal and constant over time, satisfying the fundamental assumption for valid SHPB data analysis using the one-dimensional wave theory.

The incident and transmitted bar stress-time curves were extracted from the simulation. Using the two-wave method, the stress-strain curves for ECC specimens with 2.0% fiber volume fraction were obtained across S1–S4 strain rates, as shown in Figure 9. Numerical simulations accurately replicate the key experimental trends, with the increases in peak stress and peak strain from S1 to S4 deviating from experimental measurements by a maximum of 8.2% and 6.5%, respectively. Critical mechanical parameters exhibit high-fidelity correspondence, with maximum errors of 6.8% in elastic modulus. Crucially, the model captures strain-rate-dependent damage progression through two distinct regimes. At lower strain rates (S1–S2), the fiber-governed strain-softening behavior exhibits a minimal 7% deviation in post-peak slope ($R^2$ = 0.94) – a direct correlation with microscale fiber pull-out mechanisms observed experimentally. Dynamic regimes (S3–S4) manifest matrix compaction-induced hardening, where tangent modulus predictions maintain an error within 9%. Overall, the numerical simulation results from the ECC dynamic damage constitutive model match the experimental data well. The error sources are mainly mesh size effects, where denser meshes improve accuracy and the deletion of distorted elements due to erosion criteria.

Figure 10 shows the compression failure modes of ECC specimens with 2.0% fiber volume fraction under S1–S4 strain rates, comparing experimental and numerical results. The numerical model accurately replicates the transition from localized macro-cracking at lower strain rates (S1–S2) to distributed micro-cracking and fragmentation at higher strain rates (S3–S4). This progression aligns with experimental observations where specimens evolve from exhibiting dominant shear planes to complete structural disintegration. Both experimental and simulated results demonstrate ECCs’ characteristic strain-rate dependence. As for S1–S2, diagonal shear failure predominates, and the numerical simulation successfully replicates shear band formations consistent with experimental observations, reasonably capturing both the fracture plane inclination and localized strain patterns. As for S3–S4, progressive fragmentation occurs, and the constitutive model reproduces the key features of the fragmentation failure mode in ECC specimens. The close correspondence between simulated and experimental failure morphology and fragmentation characteristics validates the fidelity of the constitutive model in capturing strain-rate-dependent damage evolution of ECCs under dynamic compression.

In this study, the dynamic compressive properties of ECCs with 2% fiber content were studied through SHPB experiments, constitutive relationship derivation, and LS-DYNA simulations. The main conclusions are summarized as follows:

(i) Both the dynamic peak stress and peak strain of ECCs showed significant strain-rate-dependent behavior. The maximum increase in strength reached 41%, while the peak strain rose by up to 98.6%. The strain-rate-dependent behavior diminished progressively as the strain rate increased. Peak strain was more significantly affected by the strain rate level than the peak stress. The impact of strain rate on peak stress diminished gradually once a certain threshold was surpassed.

(ii) The constitutive relationship, which combines a nonlinear viscoelastic model based on the ZWT framework with a Weibull distribution damage model, was obtained by fitting relevant parameters using experimental data. The results showed that it characterized effectively and safely the material's damage evolution. It achieved high simulation accuracy for the ascending pre-peak linear stage at different strain levels but lower accuracy for the softening stage.

(iii) The LS-DYNA software was used to successfully simulate the entire process of the SHPB dynamic compression test. The simulation results aligned well with the experimental data, especially at high strain rates, indicating that the simulation can guide future research on ECCs' explosion-resistant properties.

However, this study has limitations. Although the constitutive model used in this study provides a relatively conservative estimate of the softening stage, which is safe from an engineering design perspective, it still falls short in describing the characteristics of ECCs' softening stage. Further investigation into the mechanical properties of the softening stage could be conducted based on more accurate models. The findings of this study provide valuable insights for the design of protective structures subjected to blast and impact loads. The calibrated numerical model can be directly integrated into engineering design codes to predict the dynamic response of ECCs under various loading conditions.