A stochastic damage model for concrete considering strain rate effect

Abstract

To represent the rate-dependent mechanical behavior of concrete, a stochastic damage model is proposed that integrates Langevin dynamics within the Micro-Meso Stochastic Fracture framework. In this model, the fracturing process of each micro-spring is described as a barrier-crossing event where the effective barrier height depends on the external loading rate. Evolution of the reaction coordinate follows the Langevin dynamics. Solving the corresponding Fokker–Planck–Kolmogorov equation yields the failure probability of each micro-spring under dynamic loading, from which the dynamic fracture strain is derived. Numerical results demonstrate that the model naturally captures the characteristic two-stage strength enhancement of concrete: a gradual increase of dynamic increase factor at low strain rates, followed by a rapid rise beyond a critical threshold. By linking the static fracture strain of each micro-spring to its unique initial energy barrier, the model assigns distinct rate sensitivities to different material constituents. This approach enables the model to reflect different rate sensitivities in tension and compression without separate empirical adjustments. Finally, simulations of reinforced concrete beams under various loading velocities, including high-rate impact, validate the model's capability to predict structural dynamic responses in practical engineering scenarios.

Keywords

Concrete damage model strain rate effect Langevin dynamics stochastic damage

Introduction

Concrete mechanical properties exhibit pronounced sensitivity to loading rate, a phenomenon commonly termed the strain rate effect. Multiple experimental evidence confirms that both concrete tensile and compressive strength increase with increasing strain rate (Cowell, 1966; Grote et al., 2001; Klepaczko and Brara, 2001; Ross et al., 1996; Schuler et al., 2006; Tedesco and Ross, 1993). This enhancement is typically quantified by the dynamic increase factor (DIF), defined as the ratio of dynamic strength at a given strain rate to the corresponding static strength. Generally, the evolution of DIF follows a two-stage pattern: at low strain rates, DIF increases gradually and steadily; beyond a critical threshold, DIF rises sharply (Bischoff and Perry, 1991; Malvar and Ross, 1998).

A variety of damage models have been developed to capture the above behavior. Direct empirical approaches fit expressions to DIF data and embed them within static damage formulations to achieve apparent rate dependency (Holmquist et al., 1993; Riedel et al., 1999). While these models reproduce the observed strength increase, they provide limited insight into the underlying mechanisms and often lack robustness and generality beyond the calibration range. Motivated by physical understanding, alternative explanations exist. For instance, experimental findings by Reinhard (Reinhardt and Weerheijm, 1991) indicate that wet concrete has a higher rate sensitivity than dry concrete, which suggests that pore water has an influence on the strain rate effect. This phenomenon is often referred to as the Stefan effect. Guided by this observation, viscous theories have been incorporated into rate-dependent damage models. Two principal viscous frameworks have emerged: the Perzyna type (Perzyna, 1966; Ren and Li, 2013) and the Duvaut–Lions type (Duvant and Lions, 2012; Wei and Ren, 2023). The primary distinction between the two types of models lies in their respective definitions of overstress. On the other hand, experiments reported by Cadoni et al. (2001) and Brara and Klepaczko (2006) show similar rate sensitivities for dry and wet concrete under high strain rates, implying water has minimal influence on strength enhancement at these extremes. Bažant and Caner (2014) attributed high-strain-rate strength enhancement primarily to inertial effects. Building upon this theoretical framework, subsequent researchers have developed several damage constitutive models that explicitly account for inertial contributions (Feng and Hai, 2025; Hai et al., 2020; Häußler-Combe and Kitzig, 2009).

Despite viscous and inertial perspectives, reaction rate theory has been applied to dynamic fracture phenomena (Benichou and Givli, 2015; Lea and Jardine, 2018; Yang et al., 2020). This perspective holds that material failure happens as a result of a thermally activated bond-breaking process, where the system overcomes a free energy barrier at a rate that is dependent on both temperature and barrier height. Reaction rate theory fundamentally requires the barrier height to be substantially larger than $k_{B} T$ , where T is the temperature and $k_{B}$ is the Boltzmann constant. However, under high-speed loading, this high-energy barrier assumption may not hold at the time of failure, thus limiting the applicability of the theory (Hänggi et al., 1990). Moreover, when the external loading timescale approaches or falls below that of thermal fluctuations, reaction rate predictions become invalid (Kramers, 1940). To address these limitations, Feng et al. (2022) recently formulated a unified nonequilibrium statistical framework based on Langevin dynamics and the Fokker–Planck–Kolmogorov (FPK) equation, which captures the abrupt transition in concrete strength between low and high loading rates.

Notably, empirical DIF measurements consistently indicate a higher rate of sensitivity under tension compared to compression. This observed difference poses a challenge for many existing approaches, which frequently resort to independent parameter adjustments for tension and compression within a given modeling framework. Such adjustments fail to derive the distinct sensitivities from a unified physical basis, highlighting the need for a physics-driven model inherently capable of reflecting both tensile and compressive rate effects without separate calibrations.

Motivated by this need, the present study aims to extend the established Micro-Meso Stochastic Fracture (MMSF) model (Li and Ren, 2009; Li and Guo, 2023) to dynamic loading cases. It is important to note that while previous rate-dependent MMSF models exist (Guo and Li, 2024; Ren and Li, 2013; Ren et al., 2015), they typically apply a uniform DIF to the entire fracture strain random field. This approach results in an isotropic amplification of the static response, where every micro-spring shares an identical rate sensitivity.

In contrast, the present work introduces a micro-spring-level heterogeneous rate sensitivity. By linking each micro-spring's static fracture strain to its unique energy barrier height, we enable different micro-springs to exhibit distinct dynamic responses. This approach physically reflects the realistic scenario in concrete, where different microstructural components possess not only different static strengths but also distinct rate sensitivities. Through this mechanism, the static fracture strain of each micro-spring evolves into a rate-dependent quantity, leading to a fully dynamic stochastic damage formulation. Crucially, this heterogeneous coupling allows the model to inherently reproduce the distinct rate sensitivities of concrete in tension and compression.

Furthermore, we address the critical challenge of computational feasibility. While Feng et al. (2022) rigorously derived the Langevin-FPK formulation, directly solving the FPK equation at every integration point for each time step is computationally prohibitive for structural-scale simulations. This study overcomes this barrier by developing an efficient fitted expression for the FPK solutions, thus creating a practical pathway from the physical principles of Langevin dynamics to macroscopic continuum damage modeling suitable for finite-element analysis.

Continuous elastoplastic damage framework

Bi-scalar damage variable

Generally, the elastoplastic damage framework is expressed as

σ = (I - D) : E : (ε - ε^{p})

(1)

where

σ

and

ε

are the Cauchy stress and strain tensor, respectively.

ε^{p}

denotes the plastic strain tensor;

I

is the fourth-order identity tensor;

D

is the fourth-order damage tensor; and

E

is the fourth-order elastic stiffness tensor.

According to the strain equivalent assumption (Lemaitre, 1996), the damage portion can be eliminated from the effective stress space. Thus, the effective stress tensor $\bar{σ}$ can be expressed as

\bar{σ} = E : (ε - ε^{p})

(2)

Since the tensile properties of concrete are distinct from compressive properties, two scalar damage variables, $d^{\pm}$ , are introduced for the tensile stress state and the compressive stress state. Therefore, the damage tensor can be expressed in terms of bi-scalar damage variables as

D = d^{+} P^{+} + d^{-} P^{-}

(3)

where

P^{\pm}

are the projection tensors, represented as

{\begin{matrix} P^{+} = \sum_{i = 1}^{3} H ({\bar{σ}}_{i}) n^{(i)} \otimes n^{(i)} \otimes n^{(i)} \otimes n^{(i)} \\ P^{-} = I - P^{+} \end{matrix}

(4)

where

{\bar{σ}}_{i}

and

n^{(i)}

represent the ith eigenvalues and eigenvectors of the effective stress tensor, and

H (\cdot)

denotes the Heaviside function.

Correspondingly, the effective stress tensor is decomposed into the positive part and the negative part:

\bar{σ} = {\bar{σ}}^{+} + {\bar{σ}}^{-}, {\bar{σ}}^{\pm} = P^{\pm}

(5)

According to irreversible thermodynamics, the damage energy release rate is defined as the work conjugate to the damage variable. Following Wu et al. (2006), the tensile and compressive damage energy release rates are expressed as

{\begin{matrix} Y^{+} = \sqrt{E_{0} ({\bar{σ}}^{+} : E^{- 1} : \bar{σ})} \\ Y^{-} = α {\bar{I}}_{1} + \sqrt{3 {\bar{J}}_{2}} \end{matrix}

(6)

where

E_{0}

is the elastic modulus,

α

denotes the biaxial strength increase factor, usually

α = 0.1212

{\bar{I}}_{1}

and

{\bar{J}}_{2}

are the first and second invariants of the effective stress tensor.

The damage evolution law is defined as the function of the damage energy release rate:

d^{\pm} = g (Y^{\pm})

(7)

In this work, a simple empirical expression for the plastic strain is adopted (Ren et al., 2015). It is assumed that the plastic strain increment ${\dot{ε}}^{p \pm}$ is proportional to the elastic strain increment ${\dot{ε}}^{e \pm}$ :

{\dot{ε}}^{p \pm} = H (d^{\pm}) f^{\pm} (d^{\pm}) {\dot{ε}}^{e \pm}

(8)

where the proportionality function is expressed as

f^{\pm} (d^{\pm}) = ξ_{p}^{\pm} (d^{\pm})^{n_{p}^{\pm}}

(9)

with

ξ_{p}^{\pm}

and

n_{p}^{\pm}

the fitting parameters.

Micro-meso stochastic fracture model

To complete the theory, one needs to provide the damage evolution law. In this work, we adopted the stochastic damage evolution law from the MMSF (Li and Ren, 2009). According to random media assumption (Li et al., 2014; Liu et al., 2018), the representative volume element of concrete is modeled as a series of parallel micro-springs with rigid bars at both ends, as shown in Figure 1. Once the energy equivalent strain exceeds the fracture strain of the micro-spring, it breaks, and the damage evolves.

Figure 1.

Micro-meso stochastic fracture model.

To reflect the random property of concrete, the fracture strain of the micro-spring is assumed as a one-dimensional random field. It is important to clarify that the spatial variability of material strength is captured through a statistically equivalent random field, avoiding the prohibitive cost of explicit three-phase mesoscale models (Peng et al., 2023). The stochastic damage evolution law is therefore derived as

D^{\pm} = \int_{0}^{1} H (ε^{eq \pm} - Δ^{\pm} (x)) d x

(10)

where

Δ^{\pm} (x)

is the 1D random field that denotes the fracture strain of the micro-spring along the x direction.

ε^{eq \pm}

is the energy equivalent strain deduced from the damage energy release rates:

{\begin{matrix} ε^{eq +} = \frac{Y^{+}}{E_{0}} \\ ε^{eq -} = \frac{Y^{-}}{(1 - α) E_{0}} \end{matrix}

(11)

D^{\pm}

denotes the stochastic damage variables. In this study, the capital letter D is employed to represent stochastic damage, thereby distinguishing it from the previously used deterministic damage variable d.

The 1-D probability distribution of $Δ^{\pm} (x)$ is assumed to be a log-normal distribution, consistent with nano- and micro-indentation test results on concrete's constituent phases (Liu et al., 2018). Taking tensile stress state as an example, let $Z (x) = \ln Δ (x)$ , $λ$ and $ζ$ are the mean and standard deviation of $Z (x)$ , respectively. The auto-coefficient function of $Z (x)$ is assumed as

ρ_{Z} (y) = \exp (- ω | y |)

(12)

where

ω

is the correlation parameter, and

y = x_{1} - x_{2}

represents the relative distance between two micro-springs.

Therefore, the mean and standard deviation of damage evolution can be solved as

μ_{D} = Φ_{1} (\frac{\ln ε^{eq} - λ}{ζ})

(13)

V_{D}^{2} = 2 \int_{0}^{1} (1 - y) Φ_{2} (\frac{\ln ε^{eq} - λ}{ζ}, \frac{\ln ε^{eq} - λ}{ζ} | ρ_{Z} (y)) d y - μ_{D}^{2}

(14)

where

Φ_{1} (\cdot)

and

Φ_{2} (\cdot)

are the 1D and 2D cumulative probability functions of the normal distribution, respectively.

While the MMSF model effectively captures concrete's nonlinearity and stochastic nature under quasi-static conditions, its original formulation cannot account for strain rate effects. The present study addresses this gap by extending the MMSF framework to dynamic loading through the incorporation of a rate-dependent fracture mechanism within the micro-springs.

Rate-dependent fracture mechanism of the micro-spring

Langevin dynamics

Material failure is conceptualized as a rate-dependent energy barrier crossing process within reaction rate theory. Following Feng et al. (2022), this work adopts an interval-N-type energy to describe the micro-spring fracture process, as illustrated in Figure 2.

Figure 2.

The interval-N energy barrier.

The internal energy U of this micro-spring is expressed as

U (ϕ, δ) = A ϕ + \frac{1}{2} (1 - ϕ) k δ^{2}, 0 \leq ϕ \leq 1

(15)

where

ϕ

is the reaction coordinate, denoting the damage state of the micro-spring. A is the initial energy barrier height; k represents the micro-spring's stiffness; and

δ

represents its extension.

The effective energy barrier height $E_{b}$ decreases gradually with the increase of external loading.

E_{b} = A - \frac{1}{2} k δ^{2}

(16)

Obviously, when the extension reaches $δ_{c} = \sqrt{2 A / k}$ , the height of the energy barrier decreases to 0.

The evolution of the reaction coordinate is assumed to follow the Langevin equation (Feng et al., 2022):

η_{ϕ} \dot{ϕ} (t) = - U^{'} (ϕ, δ) + W_{ϕ} (t)

(17)

where

η_{ϕ}

is the damping coefficient, and the prime represents the derivative with respect to

ϕ

$W_{ϕ} (t)$ is the white noise that is satisfying:

E [W_{ϕ} (t)] = 0, E [W_{ϕ} (t) W_{ϕ} (τ)] = 2 η_{ϕ} k_{B} T δ_{D} (t - τ)

(18)

where

δ_{D}

is the Dirac delta function;

k_{B}

and T denote the Boltzmann constant and absolute temperature, respectively.

Then the probability density evolution of the reaction coordinate follows the FPK equation (Jiang et al., 2025; Risken and Risken, 1996):

\frac{\partial p (ϕ, t)}{\partial t} = - \frac{\partial}{\partial ϕ} J (ϕ, t), J (ϕ, t) = - [\frac{1}{η_{ϕ}} U^{'} (ϕ, δ) + \frac{k_{B} T}{η_{ϕ}} \frac{\partial}{\partial ϕ}] p (ϕ, t)

(19)

The boundary conditions of equation (19) are as follows:

{\begin{matrix} J (ϕ = 0, t) = 0 \\ p (ϕ = 1, t) = 0 \end{matrix}

(20)

Equation (20) indicates that a reflection boundary at $ϕ = 0$ and an absorbing boundary at $ϕ = 1$ .

For ease of analysis, the dimensionless operation is performed for equations (15), (16), and (19).

{\begin{matrix} \hat{U} (ϕ, \hat{δ}) = \hat{A} [ϕ + (1 - ϕ) {\hat{δ}}^{2} (\hat{t})] \\ {\hat{E}}_{b} (\hat{δ}) = \hat{A} (1 - {\hat{δ}}^{2} (\hat{t})) \\ \frac{\partial p (ϕ, \hat{t})}{\partial \hat{t}} = \frac{\partial}{\partial ϕ} [\hat{U}^{'} (ϕ, \hat{δ}) p (ϕ, \hat{t})] + \frac{\partial^{2} p (ϕ, \hat{t})}{\partial ϕ^{2}} \end{matrix}

(21)

where

\hat{A} = A / k_{B} T

{\hat{E}}_{b} = E_{b} / k_{B} T

\hat{U} = U / k_{B} T

\hat{t} = t / t_{ϕ}

t_{ϕ} = η_{ϕ} / k_{B} T

\hat{δ} = δ / δ_{c}

. Equation (21) uses

k_{B} T

t_{ϕ}

, and

δ_{c}

as the reference physical quantities for energy, time, and length. Likewise, the reference loading speed is defined as

v_{0} = δ_{c} / t_{ϕ}

. Thus, the normalization loading speed

\hat{v} = v / v_{0}

As the energy barrier decreases, the probability density function of the reaction coordinate is continuously absorbed by the absorbing boundary. Therefore, the failure probability of the micro-spring can be defined as

F_{R} (t) = 1 - \int_{0}^{1} p (ϕ, t) d ϕ

(22)

And the corresponding survival probability is expressed as

S_{R} (t) = 1 - F_{R} (t) = \int_{0}^{1} p (ϕ, t) d ϕ

(23)

The micro-spring is considered failed when its survival probability reaches zero. Its dynamic fracture properties will be analyzed in the following subsection.

Dynamic property of the micro-spring

Equation (21) is numerically solved via the path integral approach (Risken and Risken, 1996), and the relevant factors influencing the computation results at various loading rates can be examined.

Assuming $\hat{A} = 50$ and $\hat{v} = 1$ , Figure 3 displays the evolution of residual probability density function (RDF) during the solution process. As loading progresses, the potential barrier height progressively decreases, causing the response coordinate's probability density to shift rightward toward the absorbing boundary. The survival probability of the micro-spring and the evolution of the failure probability can be further obtained by integrating the RDF.

Figure 3.

Numerical solution of the FPK equation. (a) Evolution process of RDF, (b) Evolution process of failure probability.

The evolution of the failure probability is significantly affected by the loading rate, as shown in Figure 4. It can be found that, for the same tensile length, the micro-spring's failure probability gradually decreases with increasing loading rate. In other words, the fracture length of the micro-spring increases with increasing loading rate.

Figure 4.

Failure probability under different loading rates.

In addition to the loading rate, the initial energy barrier height is also an important factor affecting the critical micro-spring fracture length. Figure 5 demonstrates the relationship between fracture length and loading rate across varying $\hat{A}$ values, revealing a distinctive two-stage response: fracture length exhibits low-rate sensitivity below a critical loading threshold, transitioning to dramatically enhanced sensitivity beyond this limit. Furthermore, fracture length sensitivity inversely scales with $\hat{A}$ magnitude, indicating the energy barrier's moderating effect on rate-dependence.

Figure 5.

Fracture length with different initial energy barrier heights.

The characteristic two-stage evolution of the DIF observed in Figure 5 can be explained by the competing physical mechanisms inherent to the Langevin dynamics framework. At low and moderate loading rates, the failure process is dominated by thermal fluctuations overcoming the energy barrier. Mathematically, the failure probability is governed by the first eigenvalue of the FPK operator, which characterizes slow, thermally assisted barrier crossing. As shown in Figure 6, with the increase of loading rate, the influence of high-order eigenvalues becomes greater, and the first-order eigenvalues approximation on FPK solutions deteriorates. This leads to a gradual, steady increase in DIF observed at low-to-moderate strain rates.

Figure 6.

Strength increase mechanism under low loading rate.

Beyond a critical loading rate, the random thermal fluctuations become negligible compared to the deterministic external loading. The reaction coordinate evolution is no longer dominated by thermal noise but rather by the viscous term. As presented in Figure 7, with the increase of loading rate, the viscous term approximation on the FPK solution becomes better. Mathematically, this regime is characterized by contributions from higher-order eigenvalues of the FPK operator, which correspond to faster, non-exponential barrier-crossing dynamics. The result is a sharp, nonlinear increase in DIF at high strain rates.

Figure 7.

Strength increase mechanism under high loading rate.

Fitting the dynamic fracture strain

Within the MMSF framework, micro-spring fracture is governed by energy equivalent strain, while Langevin dynamics determines failure through elongation. To reconcile these approaches, this study establishes the quasi-static condition at a normalization loading rate ${\hat{v}}_{norm} = 5 \times 10^{- 5}$ , which corresponds to a strain rate of $\dot{ε} = 1 \times 10^{- 5} s^{- 1}$ . The DIF $Δ_{DIF}$ for any loading rate is defined as:

Δ_{DIF} = \frac{δ_{v}}{δ_{v_{norm}}}

(24)

where

δ_{v}

represents elongation at the target strain rate and

δ_{v_{norm}}

denotes elongation at the normalization rate.

This coordinate mapping transforms dynamic elongation (vertical axis) into fracture strain DIF while converting normalized loading rate (horizontal axis) to strain rate, with Figure 8 presenting the resulting strain-rate-dependent enhancement relationships.

Figure 8.

Dynamic fracture strain conversion. (a) Dynamic elongation of the micro-spring, (b) Dynamic increase factor of the micro-spring.

In the numerical implementation of the MMSF framework, the fracture strain random field needs to be discretized by the stochastic harmonic function method (Chen et al., 2013). Recognizing that micro-springs with distinct static fracture strains exhibit differential rate sensitivities, initial energy barrier heights must be uniquely assigned. This study establishes a linear relationship between the static fracture strain and initial potential barrier height. Defining a reference micro-spring with standard fracture strain $Δ_{1}$ and corresponding barrier height ${\hat{A}}_{1}$ , the barrier height for any micro-spring i is given by:

{\hat{A}}_{i} = \frac{Δ_{i}}{Δ_{1}} {\hat{A}}_{1}

(25)

where

Δ_{i}

denotes the static fracture strain of micro-spring i.

This assumption is conceptually grounded in the fact that concrete is a composite material where different phases possess distinct strengths. It is physically intuitive that a stronger material region, which requires more energy to fail and thus has a higher static fracture strain, would also present a higher activation energy barrier against dynamic failure. The linear form is adopted as the simplest, first-order approximation of this positive correlation, avoiding the introduction of additional parameters with unclear physical meaning.

This physically motivated assumption is the key to capturing the tension-compression asymmetry in rate sensitivity. Since concrete's tensile fracture strain is much lower than its compressive fracture strain, the corresponding energy barrier for tensile failure is also significantly lower. According to the principles of Langevin dynamics, materials with lower energy barriers exhibit higher rate sensitivity because the transition to the viscosity-dominated regime occurs at a lower strain rate. This mechanism allows the model to naturally reproduce the distinct rate sensitivities of tension and compression.

In this way, the dynamic fracture strain of the micro-spring can be expressed as

Δ_{d} (x) = Δ_{DIF} \times Δ_{s} (x)

(26)

where the subscripts s and d denote static and dynamic loading, respectively.

Substituting the dynamic fracture strain derived in equation (26) into (10) yields a rate-dependent stochastic damage evolution law, thereby enabling explicit representation of concrete's strain rate effects. Obviously, the suggested model requires a significant computational cost for obtaining the relevant dynamic fracture strain increase factor, which requires solving the corresponding FPK equations for each micro-spring at each incremental step. To enhance computational efficiency, this study develops a generalized fitting function for the micro-spring's dynamic strength factor across all relevant loading rates and energy barrier heights. The fitting expression takes the form:

Δ_{DIF} = α_{DIF} (\hat{A}) ({\dot{ε}}^{β_{DIF} (\hat{A})} - {\dot{ε}}_{s}^{β_{DIF} (\hat{A})}) + 1

(27)

where

α_{DIF}

and

β_{DIF}

are the fitting parameters, respectively, which are functions of the energy barrier height

\hat{A}

\dot{ε}

denotes the arbitrary strain rate, and

{\dot{ε}}_{s}

is the quasi-static strain rate; here

{\dot{ε}}_{s} = 1 \times 10^{- 5} s^{- 1}

. Figure 9 illustrates the typical fitting results with different energy barrier heights.

Figure 9.

Dynamic increase factor fitting of fracture strain. (a) $\hat{A} = 20$ , (b) $\hat{A} = 50$ , (c) $\hat{A} = 200$ , (d) $\hat{A} = 1000$ .

Given that the parameters $α_{DIF}$ and $β_{DIF}$ in equation (27) are directly related to the barrier height, it is necessary to further establish a fitted expression between these two parameters and the energy barrier height. The specific expressions are as follows:

α_{DIF} (\hat{A}) or β_{DIF} (\hat{A}) = a_{α / β} {\hat{A}}^{b_{α / β}} + c_{α / β} {\hat{A}}^{d_{α / β}}

(28)

where the fitting parameters

a_{α} = 9.572 \times 10^{10}

b_{α} = - 9.986

c_{α} = 9.778

d_{α} = - 0.683

a_{β} = 0.298

b_{β} = 0.040

c_{β} = - 4.085

d_{β} = - 2.422

. The subscripts of the fitted parameters indicate which variable is being targeted. The final fit is shown in Figure 10.

Figure 10.

Fitting results for parameters $α_{DIF}$ and $β_{DIF}$ .

To clarify the role of the various parameters used in the model, they can be divided into two distinct categories. The first category consists of material parameters (e.g. $λ^{\pm}$ , $ζ^{\pm}$ , $ω^{\pm}$ ), which are calibrated from standard quasi-static material tests following established procedures in the MMSF literature (Chen and Li, 2023; Li and Ren, 2009). The second category includes numerical coefficients (e.g. $α_{DIF}$ , $β_{DIF}$ , and the fitting parameters in equation (28)). These coefficients are purely mathematical artifacts derived from fitting the numerical solutions of the FPK equation. Their sole purpose is to construct a computationally efficient surrogate model that avoids solving the FPK equation at every time step during simulations. These coefficients carry no direct physical meaning and, in principle, would be unnecessary if computational resources were unlimited.

It is important to note that using the fitting function (equation (27)) neglects the history dependence inherent to the full FPK equation. This approximation is most accurate for monotonic, high-rate loading scenarios where the strain rate does not vary drastically. However, for complex loading paths involving significant variations or reversals in strain rate, such as cyclic dynamic loading, this fitting approach would no longer be suitable, and the full FPK equation would need to be solved.

In the proposed model, the elastic modulus exhibits typical rate sensitivity in addition to the fracture strain of the micro-spring. For ease of simplicity, the empirical equations for the dynamic elastic modulus from the Euro-International Committee for Concrete (CEB) recommendations (Bischoff and Perry, 1991) are directly adopted here:

\frac{E_{0, d}}{E_{0, s}} = {(\frac{{\dot{ε}}_{d}}{{\dot{ε}}_{s}})}^{0.0 26}

(29)

where the subscripts s and d denote static and dynamic loading, respectively.

Case studies

Numerical scheme

An explicit numerical scheme is developed to implement the damage model. Based on the operator split method (Simo and Hughes, 1998), the stress update process can be decomposed into three steps: elastic prediction, plastic correction, and damage correction. In the first step, the trial effective stress is calculated by freezing the damage and plastic strain:

{\begin{matrix} {\bar{σ}}_{n + 1}^{trial} = E_{0} : (ε_{n + 1} - ε_{n + 1}^{p}) \\ ε_{n + 1} = ε_{n} + Δ ε_{n + 1} \\ ε_{n + 1}^{p \pm} = ε_{n}^{p \pm} \\ D_{n + 1}^{\pm} = D_{n}^{\pm} \end{matrix}

(30)

where the subscripts n and n + 1 denote the physical quantities of incremental step n and incremental step n + 1, respectively.

We can update the stress directly if $Y_{n + 1}^{trial \pm} \leq r_{n}^{\pm}$ , where $r_{n}^{\pm} = max {r_{0}^{\pm}, max_{τ \in [0, n]} Y_{τ}^{\pm}}$ are the damage threshold. When $Y_{n + 1}^{trial \pm} > r_{n}^{\pm}$ , plasticity correction and damage correction steps are required. The plastic stress is defined as $σ^{p \pm} = E_{0} : ε^{p \pm}$ . Rewriting this into incremental form as

σ_{n + 1}^{p \pm} = E_{0} : (ε_{n}^{p \pm} + Δ ε_{n + 1}^{p \pm}) = σ_{n}^{p \pm} + η_{n + 1}^{\pm} ({\bar{σ}}_{n + 1}^{\pm} - {\bar{σ}}_{n}^{\pm})

(31)

where the plastic factor

η_{n + 1}^{\pm} = ξ_{p}^{\pm} D_{n + 1}^{n_{p}^{\pm}}

The effective stress after plastic correction can be expressed as

{\bar{σ}}_{n + 1}^{\pm} = E_{0} : (ε_{n + 1}^{\pm} - ε_{n + 1}^{p \pm}) = {\bar{σ}}_{n + 1}^{trial} + σ_{n}^{p \pm} - σ_{n + 1}^{p \pm}

(32)

Substituting equation (32) into equation (31), the plastic stress can be solved as

σ_{n + 1}^{p \pm} = σ_{n}^{p \pm} + \frac{η_{n + 1}^{\pm}}{η_{n + 1}^{\pm} + 1} ({\bar{σ}}_{n + 1}^{trial \pm} - {\bar{σ}}_{n}^{\pm})

(33)

After the plasticity correction, the damage energy release rate needs to be recalculated, and the dynamic strength increase factor of the micro-spring is calculated based on the loaded strain rate at this increment. Finally, the damage variables and stresses are updated. The numerical scheme flow chart is presented in Figure 11.

Figure 11.

Numerical scheme flow chart.

Uniaxial test

Uniaxial tensile simulations were conducted to validate the model's description of dynamic concrete behavior and its representation of material stochasticity. Three distinct strain rates were examined, with 50 samples generated per strain rate using the harmonic function approach (Chen et al., 2013) to capture inherent variability. Figure 12 displays a comparison between the mean and standard deviation determined from the simulated findings and the experimental results (Gao et al., 2021). The simulation parameters are as follows: $E_{0, s} = 3.4 \times 10^{4} N / m m^{2}$ , $λ^{+} = 5.2$ , $ζ^{+} = 0.51$ , $ω^{+} = 40$ , $Δ_{1} = 1 \times 1 0^{- 5}$ , ${\hat{A}}_{1} = 4$ .

Figure 12.

Uniaxial tensile simulation. (a) ${\dot{ε}}_{1} = 1 \times 10^{- 5} / s$ , (b) ${\dot{ε}}_{2} = 1 \times 10^{- 4} / s$ , (c) ${\dot{ε}}_{3} = 1 \times 10^{- 2} / s$ .

Similarly, for the uniaxial compression case, three different strain rates were calculated. The statistical characteristics are compared to the experimental result (Zeng et al., 2013), as shown in Figure 13. The simulation parameters are as follows: $E_{0, s} = 3.0 \times 10^{4} N / m m^{2}$ , $λ^{-} = 7.7$ , $ζ^{-} = 0.18$ , $ω^{-} = 40$ , $ξ_{p}^{-} = 1.4$ , $n_{p}^{-} = 0.2$ , $Δ_{1} = 1 \times 1 0^{- 5}$ , ${\hat{A}}_{1} = 4$ .

Figure 13.

Uniaxial compressive simulation. (a) ${\dot{ε}}_{1} = 1 \times 10^{- 5} / s$ , (b) ${\dot{ε}}_{2} = 1 \times 10^{- 4} / s$ , (c) ${\dot{ε}}_{3} = 3.5 \times 10^{- 2} / s$ .

Critically, Figures 12 and 13 show that simulated responses not only accurately capture the experimental mean behavior across strain rates but also successfully reproduce the statistical variability observed in concrete tests. This demonstrates the model's dual capability in representing both rate-dependent mechanics and intrinsic material randomness.

The DIF of concrete strength was computed across an extensive strain rate range to assess model performance. Analysis of 20 representative numerical samples confirms the proposed formulation accurately reproduces both the characteristic biphasic DIF evolution and inherent material stochasticity, with most experimental data points residing within ±2 standard deviations of the simulated mean values, as presented in Figure 14. These DIF experimental data were collected from Dilger et al. (1984), Guo et al. (2017), Levi Hevroni et al. (2018), Ma et al. (2020), Tang et al. (2020), Wu et al. (2005), Yan et al. (2006), Zhang et al. (2009), etc. This dual validation establishes the model’s capability to characterize concrete strength enhancement under dynamic loading while capturing intrinsic randomness. Notably, identical parameter sets successfully describe tension-compression strength increments without requiring separate calibrations, demonstrating unified predictive capacity across loading modes.

Figure 14.

DIF under tension and compression. (a) Uniaxial tensile DIF, (b) Uniaxial compressive DIF.

Four-point bending test of RC beam

Numerical simulations of the four-point bending test of reinforced concrete beams without webs by Kulkarni and Shah (1998) were conducted using the proposed dynamic damage model. The relevant finite-element model is created using the reinforcing details and beam geometry shown in Figure 15. Two beam series with different shear-span ratios were selected for simulation and analysis to verify model applicability, with dimensional and loading parameters detailed in Table 1.

Figure 15.

Four-point bending reinforced concrete beam dimensions.

Table 1.

Loading information of the reinforced concrete beam.

Set	Specimen number	a/d	z (mm)	$f_{c}$ (MPa)	V (mm/s)	Failure mode
1	B3DE03-S	5.0	152	46.2	0.0071	Flexure
1	B3DE03-H	5.0	152	46.2	380	Flexure
2	B3NO15-S	4.0	152	43.0	0.0071	Shear
2	B3NO15-H	4.0	152	43.0	380	Flexure

To address the mesh sensitivity problem of the local damage model, this study adopts the viewpoint of the crack band model (Bažant and Oh, 1983). The principle is to ensure that the fracture energy dissipated per unit crack area remains constant regardless of element size. In the MMSF model, where the strain energy density of the softening material is proportional to the square of the fracture strain, this regularization is achieved by scaling the fracture strain random field at each integration point. Specifically, the fracture strain $Δ_{i}$ sampled from the random field is adjusted as: $Δ_{i}^{reg} = Δ_{i} \cdot \sqrt{h_{ref} / h}$ , where $h_{ref}$ is the reference mesh size used for calibration and h is the characteristic length of the current element. Following Zhang et al. (2020), the tensile stiffening effect from reinforcement in concrete beams is incorporated by assigning twice the fracture energy of plain concrete to reinforced concrete. Numerical simulations employ 10-mm triangular plane stress elements for concrete and 80-mm truss elements for reinforcement. The bond slip between the steel reinforcement and the concrete is not considered, and only the mean value of the stochastic damage evolution of the concrete is considered in the simulation. The concrete damage model parameters are as follows: $E_{0, s} = 3.4 \times 10^{4} N / m m^{2}$ , $λ^{+} = 5.2 5$ _, $ζ^{+} = 0.5$ , $ω^{+} = 40$ , $ξ_{p}^{+} = 40$ , $n_{p}^{+} = 2$ , $λ^{-} = 7.9$ , $ζ^{-} = 0.25$ , $ω^{-} = 40$ , $ξ_{p}^{-} = 30$ , $n_{p}^{-} = 1$ , $Δ_{1} = 1 \times 1 0^{- 5}$ , ${\hat{A}}_{1} = 4$ .

Reinforcement follows a linear hardening model with elastic modulus $E_{steel} = 2.2 \times 10^{5} N / m m^{2}$ , yield strength $f_{y} = 490 MPa$ , and plastic hardening modulus $E_{p} = 0.02 E_{steel}$ under static conditions, with $f_{y}$ increased to 550 MPa for dynamic cases to account for strain-rate sensitivity.

Figures 16 to 21 display tensile damage contours and midspan load–displacement curves for two beam series. The simulation results successfully capture distinct failure modes across shear-span ratios and accurately represent nonlinear development under varying loading rates. The load–displacement curves obtained from the simulation are in good agreement with the experimental results. The B3DE03 series beams exhibit flexural failure under both static and dynamic loading conditions. However, under dynamic loading, the load-carrying capacity and ductility of the members were increased due to the increase in material strength. As for the B3NO15 series beams, failure transitions from shear mode under static loading to flexural mode under dynamic loading, with simulations accurately reflecting this rate-dependent failure mode shift.

Figure 16.

Tensile damage contour of the B3DE03-S specimen.

Figure 17.

Tensile damage contour of the B3DE03-H specimen.

Figure 18.

B3DE03 series load–displacement curve.

Figure 19.

Tensile damage contour of the B3NO15-S specimen.

Figure 20.

Tensile damage contour of the B3NO15-H specimen.

Figure 21.

B3NO15 series load–displacement curve.

RC beam under drop-weight impact

To further validate the proposed model, the mechanical response of a reinforced concrete beam under drop-weight impact (Leppänen et al., 2020) was simulated. The dimensions and reinforcement details of the beam are shown in Figure 22. The drop hammer was released from a height of 5.5 m with a weight of 10.1 kg, and the average impact velocity during the test was 10.35 m/s. The experimentally obtained tensile strength of the concrete was 3.28 MPa, compressive strength 45.5 MPa, and yield strength of the reinforcement 500 MPa. The parameters used in the simulation are as follows: $E_{0, s} = 3.2 \times 10^{4} N / m m^{2}$ , $λ^{+} = 5.2 5$ _, $ζ^{+} = 0.5$ , $ω^{+} = 40$ , $ξ_{p}^{+} = 10$ , $n_{p}^{+} = 0.5$ , $λ^{-} = 7.6$ , $ζ^{-} = 0.25$ , $ω^{-} = 40$ , $ξ_{p}^{-} = 30$ , $n_{p}^{-} = 1$ , $Δ_{1} = 1 \times 1 0^{- 5}$ , ${\hat{A}}_{1} = 4$ . Penalty contact was applied between the supports and the beam with a friction coefficient of 0.15. Figure 23 presents the development of tensile damage contour during the impact of one sample. The overall response of the beam can be considered as a fixed-end beam with an effective span l, which continuously shifts toward the supports as loading progresses until the entire beam is affected. At the moment of impact, a conical punching zone forms at the mid-span of the beam. Subsequently, as the effective span increases, the upper part of the beam becomes a negative bending moment region, forming vertical cracks under tensile stress. When the effective span reaches the support, tensile cracking occurs at the bottom of the beam. Ultimately, due to the absence of vertical constraints at the supports, the beam rebounds upward and detaches from the supports.

Figure 22.

Dimensions and reinforcement details of the beam.

Furthermore, a mesh sensitivity analysis was performed for this impact case. As shown in Figure 24, the displacement–time curves obtained from both a coarse and a fine mesh are in excellent agreement, confirming the mesh objectivity of the proposed model. Finally, a stochastic analysis was conducted by simulating ten numerical samples. The results are presented in Figure 25. The mean response bounded by ±2 standard deviation shows good agreement with the experimental data. This not only validates the model for high-rate impact but also highlights its capability to capture the stochastic nature of structural responses.

Figure 23.

Tensile damage contour of RC beam under impact.

Figure 24.

Mesh-dependent analysis. (a) Coarse mesh and fine mesh in FE simulation, (b) Deflection in the midpoint of the beams.

Figure 25.

Stochastic analysis for the impact case.

Concluding remarks

This study developed a rate-dependent stochastic damage model for concrete by coupling Langevin dynamics with the MMSF framework. Within this formulation, micro-spring fracture is conceptualized as a loading-rate-modulated energy barrier crossing process. By linking initial energy barrier heights to static fracture strains of the micro-springs, a heterogeneous rate-sensitivity mechanism was established. This allows the model to intrinsically replicate concrete's characteristic two-stage strength enhancement and capture the distinct rate sensitivities under tensile and compressive loading without separate empirical adjustments. Numerical examples, including uniaxial tests and reinforced concrete beam simulations under various loading velocities up to impact conditions, validate the model's predictive capability and demonstrate good agreement with reported dynamic strength data. This proposed approach provides enhanced physical insight into rate-dependent damage mechanisms.

Several limitations of the current study should be acknowledged. The use of a fitted function for the DIF, while computationally efficient, is an approximation that neglects the history dependence inherent to the full FPK equation. This simplification makes the current model most suitable for monotonic loading paths but less accurate for complex histories involving significant strain rate variations. Furthermore, while the model has been validated for impact loading, it does not yet account for the effects of high confinement pressure, which are critical in many blast and penetration scenarios. Future work will focus on developing more sophisticated computational strategies to address history dependence and incorporating confinement effects to extend the model's applicability to a wider range of extreme loading conditions.

Footnotes

ORCID iD

Jie Li

Author contributions

Chenggong Guo: conceptualization, methodology, software, visualization, writing—original draft, writing—review & editing, funding acquisition. Ye Feng: conceptualization, methodology, visualization, writing—original draft, writing—review & editing. Jie Li: conceptualization, methodology, writing—review & editing, funding acquisition.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (Grant Number: 51538010) and Innovative projects of Key Disciplines of Civil Engineering of Changsha University of Science and Technology (Grant Number: 25ZDXK04).

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data availability

Data will be made available on request.

Appendix

References

Bažant

Caner

(2014) Impact comminution of solids due to local kinetic energy of high shear strain rate: I. Continuum theory and turbulence analogy. Journal of the Mechanics and Physics of Solids 64: 223–235.

Bažant

(1983) Crack band theory for fracture of concrete. Matériaux et Construction 16: 155–177.

Benichou

Givli

(2015) Rate dependent response of nanoscale structures having a multiwell energy landscape. Physical Review Letters 114(9): 095504.

Bischoff

Perry

(1991) Compressive behaviour of concrete at high strain rates. Materials and Structures 24: 425–450.

Brara

Klepaczko

(2006) Experimental characterization of concrete in dynamic tension. Mechanics of Materials 38(3): 253–267.

Cadoni

Labibes

Albertini

, et al. (2001) Strain-rate effect on the tensile behaviour of concrete at different relative humidity levels. Materials and Structures 34(1): 21–26.

Chen

Sun

, et al. (2013) Stochastic harmonic function representation of stochastic processes. Journal of Applied Mechanics 80(1): 011001.

Chen

(2023) Identification of probabilistic distribution parameters for the mesoscopic stochastic fracture model. Probabilistic Engineering Mechanics 71: 103415.

Cowell

(1966) Dynamic properties of plain Portland cement concrete. California: US Naval Civil Engineering Laboratory.

10.

Dilger

Koch

Kowalczyk

(1984) Ductility of plain and confined concrete under different strain rates. Journal Proceedings 81(1): 73–81.

11.

Duvant

Lions

(2012) Inequalities in mechanics and physics. New York: Springer Science & Business Media.

12.

Feng

Fan

Tadmor

(2022) A rigorous universal model for the dynamic strength of materials across loading rates. Journal of the Mechanics and Physics of Solids 159: 104715.

13.

Feng

Hai

(2025) 3D phase-field cohesive fracture: Unifying energy, driving force, and stress criteria for crack nucleation and propagation direction. Journal of the Mechanics and Physics of Solids 196: 106036.

14.

Gao

Zhou

Ren

, et al. (2021) Rate effect on the stress–strain behavior of concrete under uniaxial tensile stress. Structural Concrete 22: E815–E830.

15.

Grote

Park

Zhou

(2001) Dynamic behavior of concrete at high strain rates and pressures: I. Experimental characterization. International Journal of Impact Engineering 25(9): 869–886.

16.

Guo

(2024) A unified stochastic damage model for concrete based on multi-scale energy dissipation analysis. Science China Technological Sciences 67(3): 863–877.

17.

Guo

Gao

Jing

, et al. (2017) Response of high-strength concrete to dynamic compressive loading. International Journal of Impact Engineering 108: 114–135.

18.

Hai

(2020) A phase-field damage model with micro inertia effect for the dynamic fracture of quasi-brittle solids. Engineering Fracture Mechanics 225: 106821.

19.

Hänggi

Talkner

Borkovec

(1990) Reaction-rate theory: Fifty years after Kramers. Reviews of Modern Physics 62(2): 51.

20.

Häußler-Combe

Kitzig

(2009) Modeling of concrete behavior under high strain rates with inertially retarded damage. International Journal of Impact Engineering 36(9): 1106–1115.

21.

Holmquist

Johnson

Cook

(1993) A computational constitutive model for concrete subjected to large strains, high strain rates, and high pressures. In: 14th International Symposium on Ballistics, pp.1–10. Quebec, Canada.

22.

Jiang

Zhang

Beer

, et al. (2025) First excursion probability of dynamical systems: A review on computational methods. Mechanical Systems and Signal Processing 232: 112751.

23.

Klepaczko

Brara

(2001) An experimental method for dynamic tensile testing of concrete by spalling. International Journal of Impact Engineering 25(4): 387–409.

24.

Kramers

(1940) Brownian motion in a field of force and the diffusion model of chemical reactions. Physica 7(4): 284–304.

25.

Kulkarni

Shah

(1998) Response of reinforced concrete beams at high strain rates. Structural Journal 95(6): 705–715.

26.

Lea

Jardine

(2018) Characterisation of high rate plasticity in the uniaxial deformation of high purity copper at elevated temperatures. International Journal of Plasticity 102: 41–52.

27.

Lemaitre

(1996) A course on damage mechanics. Berlin, Heidelberg: Springer.

28.

Leppänen

Johansson

Grassl

(2020) On the dynamic response of reinforced concrete beams subjected to drop weight impact. Finite Elements in Analysis and Design 180: 103438.

29.

Levi-Hevroni

Kochavi

Kofman

, et al. (2018) Experimental and numerical investigation on the dynamic increase factor of tensile strength in concrete. International Journal of Impact Engineering 114: 93–104.

30.

Guo

(2023) A unified stochastic damage model for concrete under monotonic and fatigue loading. International Journal of Fatigue 175: 107766.

31.

Ren

(2009) Stochastic damage model for concrete based on energy equivalent strain. International Journal of Solids and Structures 46(11-12): 2407–2419.

32.

Chen

(2014) Stochastic damage mechanics of concrete structures. Beijing: Science Press.

33.

Liu

Ren

(2018) Indentation tests based multi-scale random media modeling of concrete. Construction and Building Materials 168: 209–220.

34.

, et al. (2020) Experimental and three-dimensional mesoscopic investigation of coral aggregate concrete under dynamic splitting-tensile loading. Materials and Structures 53: 1–20.

35.

Malvar

Ross

(1998) Review of strain rate effects for concrete in tension. Materials Journal 95(6): 735–739.

36.

Peng

Jia

, et al. (2023) Dynamic behavior of UHPC member under lateral low-velocity impact: Mesoscale analysis. International Journal of Impact Engineering 172: 104418.

37.

Perzyna

(1966) Fundamental problems in viscoplasticity. Advances in Applied Mechanics 9: 243–377.

38.

Reinhardt

Weerheijm

(1991) Tensile fracture of concrete at high loading rates taking account of inertia and crack velocity effects. International Journal of Fracture 51(1): 31–42.

39.

Ren

(2013) A unified dynamic model for concrete considering viscoplasticity and rate-dependent damage. International Journal of Damage Mechanics 22(4): 530–555.

40.

Ren

Zeng

(2015) A rate-dependent stochastic damage–plasticity model for quasi-brittle materials. Computational Mechanics 55: 267–285.

41.

Riedel

Thoma

Hiermaier

, et al. (1999) Penetration of reinforced concrete by BETA-B-500 numerical analysis using a new macroscopic concrete model for hydrocodes. In: 9th International Symposium on the Effects of Munitions with Structures, pp.1–8. Berlin, Germany.

42.

Risken

(1996) Fokker-planck equation. Berlin, Heidelberg: Springer.

43.

Ross

Jerome

Tedesco

, et al. (1996) Moisture and strain rate effects on concrete strength. Materials Journal 93(3): 293–300.

44.

Schuler

Mayrhofer

Thoma

(2006) Spall experiments for the measurement of the tensile strength and fracture energy of concrete at high strain rates. International Journal of Impact Engineering 32(10): 1635–1650.

45.

Simo

Hughes

TJR

(1998) Computational immortality. New York: Springer.

46.

Tang

Tam

VWY

, et al. (2020) Investigation on dynamic mechanical properties of fly ash/slag-based geopolymeric recycled aggregate concrete. Composites Part B: Engineering 185: 107776.

47.

Tedesco

Ross

(1993) Experimental and numerical analysis of high strain rate splitting-tensile tests. Materials Journal 90(2): 162–169.

48.

Wei

Ren

(2023) Coupled viscosity-damage model for concrete under high strain rate. Engineering Fracture Mechanics 277: 108985.

49.

Zhang

Huang

, et al. (2005) Experimental and numerical investigation on the dynamic tensile strength of concrete. International Journal of Impact Engineering 32(1-4): 605–617.

50.

Faria

(2006) An energy release rate-based plastic-damage model for concrete. International Journal of Solids and Structures 43(3-4): 583–612.

51.

Yan

Lin

(2006) Dynamic properties of concrete in direct tension. Cement and Concrete Research 36(7): 1371–1378.

52.

Yang

Liechti

Huang

(2020) A multiscale cohesive zone model for rate-dependent fracture of interfaces. Journal of the Mechanics and Physics of Solids 145: 104142.

53.

Zeng

Ren

(2013) Triaxial behavior of concrete subjected to dynamic compression. Journal of Structural Engineering 139(9): 1582–1592.

54.

Zhang

, et al. (2009) Further investigation on the dynamic compressive strength enhancement of concrete-like materials based on split Hopkinson pressure bar tests. Part I: Experiments. International Journal of Impact Engineering 36(12): 1327–1334.

55.

Zhang

Cheng

Wang

, et al. (2020) Modeling the hysteretic responses of RC shear walls under cyclic loading by an energy-based plastic-damage model for concrete. International Journal of Damage Mechanics 29(1): 184–200.