Capturing full-process creep in cemented backfill: a modified generalized Kelvin model with fractional derivative and hardening function

1 Introduction

Cemented backfill mining is an important construction process for green mining, in which water, cementing material and coal gangue are mixed and cured to form a cemented backfill to replace coal pillars, thus realizing the goal of effectively supporting the overlying strata [1]. Recent advancements in geotechnical and mining engineering have highlighted the importance of innovative control strategies and material technologies. These approaches are critical for addressing complex engineering challenges. For instance, [2] proposed a three-dimensional fuzzy control system combined with an Avoiding Peak Filling Valley strategy, achieving a 20.69% reduction in drainage costs for coal mines by dynamically optimizing water level management. Concurrently, studies on rock mechanics under extreme conditions, such as [3] thermal-mechanical triaxial compression tests on granite, revealed how temperature and pressure govern deformation patterns and thermal degradation mechanisms, providing insights for deep underground engineering. Further, [4] highlighted the importance of geocell-soil interactions in subgrade reinforcement, showing how material properties and pocket geometries influence load-bearing capacity. [5] advanced topology optimization for lightweight tensegrity structures, while [6] developed accelerated creep models for polymer-alloy geocells, offering new methodologies for structural efficiency and long-term performance prediction. In addition, the time-dependent behavior of deep underground structures has attracted increasing attention. [7] investigated the time-dependent behavior of deep roadway surrounding rock by developing a coupled elastoplastic-damage model that accounts for excavation and mining disturbances. [8] further proposed a new elasto-visco-plastic damage model for predicting the long-term deformation of deep tunnels, and validated its effectiveness through numerical simulations. These studies collectively underscore the necessity of integrating multi-disciplinary approaches—from dynamic control systems to advanced material modeling—to enhance safety, sustainability, and cost-effectiveness in engineering practices.

Currently, research on the creep properties of cemented backfill materials remains limited. [9] conducted graded loading creep tests on cemented backfill under different water pressures. The results showed that the long-term strength of the cemented backfill decreased with increasing water pressure, exhibiting a nonlinear relationship with the water pressure. [10] further carried out triaxial compression creep tests on cemented backfill under seepage-stress coupling and constructed a creep constitutive equation that accounts for seepage pressure and stress level. [11] investigated the effects of physicochemical properties and mineral components of sulfur-containing tailings on the creep strength of cemented paste materials. [12] examined the rheological properties of cemented paste materials during pipeline transport and characterized the rheological behavior of tailings during deep thickening using compressive yield stress, gel point, and interference settling coefficient. [13] conducted graded loading creep tests on gangue cemented backfill under high stress intensity ratios. Their results indicated that the creep strength of the backfill was slightly higher than its uniaxial compressive strength, with the ratio ranging from 108.9% to 116.5%, and that the gangue cemented backfill exhibited obvious creep hardening characteristics. Similarly, [14] investigated the time-dependent behavior and creep-damage model of basalt-fiber-reinforced high-water material under uniaxial compression, revealing the coupled creep-damage evolution mechanism in cementitious materials. However, it remains challenging to reproduce the creep hardening-damage evolution law of cemented backfill materials, as most existing studies only explore the deformation behavior from a single perspective—either creep hardening or damage characteristics. Additionally, previous investigators [15] developed integer-order constitutive models to describe the creep properties of cemented backfill materials. Nevertheless, these models struggle to characterize the deformation behavior during the accelerated creep stage. The application of fractional calculus has recently expanded beyond traditional fields, demonstrating significant advantages in modeling complex systems with memory-dependent and non-local behaviors. Recent studies highlight its versatility in diverse areas such as communicable disease models [16], time dilation [17], water treatment [18], and fractional systems [19–21]. Particularly in civil and mining engineering, fractional-order models have proven effective in capturing the time-dependent mechanical responses of geomaterials due to their ability to interpolate between solid-like and fluid-like states with fewer empirical parameters [22, 23]. These advancements underscore the potential of fractional operators to overcome the limitations of integer-order models in representing hardening-damage coupling phenomena. Therefore, considering the coupled creep hardening-damage characteristics of cemented backfill materials, it is imperative to thoroughly investigate a fractional-order creep constitutive model that can capture the full creep process.

2 Preliminaries

The development of fractional derivatives has led mathematicians to propose diverse definitions from various theoretical perspectives. While each definition offers unique advantages, most are constrained to specific applications and exhibit inherent limitations. Currently, no universally accepted standard exists for fractional derivative definitions. Nevertheless, extensive research has revealed a unifying characteristic among these operators: they can all be expressed as convolution integrals with power-law kernels. Below, we briefly introduce the most widely used definitions of fractional derivatives.

• Grünwald-Letnikov fractional derivative

Liouville [24] first proposed the general form of Grünwald-Letnikov’s definition, which was later elaborated more specifically by [25] and [26] in 1867 and 1868, respectively. The Grünwald Letnikov definition is often used in numerical calculations and can be seen as an extension of the limit form of integer order calculus difference definitions [27].

G L α D t α f x = lim h → 0 + 1 h α ∑ r = 0 t ‐ α h ‐ 1 r α r f t ‐ r h

G L t D b α f x = lim h → 0 + 1 h α ∑ r = 0 b ‐ t h ‐ 1 r α r f t + r h

• Riemann–Liouville fractional derivative

The Grünwald-Letnikov fractional derivative is defined as a difference scheme, when f(t) has an order of m+1 and m takes ⌊ α ⌋ = n − 1, the Grünwald-Letnikov fractional derivative is equivalent to the Riemann–Liouville fractional derivative. The Riemann Liouville fractional derivative adopt the form of differential integral, which has a wider range of applications and plays an important role in mathematical theory [28].

D a + α f x = 1 Γ n − α d n d x n ∫ a x x − ξ n − α − 1 f ξ d ξ , x ≥ a , n − 1 < α < n R L

D b − α f x = − 1 n Γ n − α d n d x n ∫ x b ξ − x n − α − 1 f ξ d ξ , x ≤ b , n − 1 < α < n R L

• Caputo fractional derivative

Italian geophysicist Caputo proposed a fractional order differential definition in the late 1960s to solve the fractional order initial value problem in Riemann Liouville’s fractional definition. The advantage of Caputo type derivatives [29] is that the initial conditions that need to be satisfied in modeling applications and integral transformations are given in the form of integer order calculus.

D a + α f x = 1 Γ n − α ∫ a x x − ξ n − α − 1 f n ξ d ξ , x ≥ a , n − 1 < α < n c

D b − α f x = − 1 n Γ n − α ∫ x b ξ − x n − α − 1 f n ξ d ξ , x ≤ b , n − 1 < α < n c

The Riemann–Liouville and Caputo fractional derivatives represent refined formulations derived from the Grünwald-Letnikov definition. Their fundamental distinction lies in the sequence of operations: the Riemann–Liouville derivative applies integration prior to differentiation, whereas the Caputo operator reverses this order. Mathematically, the Caputo definition imposes stricter requirements, necessitating the existence of n-order derivatives for the function. This formulation offers a significant advantage in engineering applications, as its Laplace transform maintains consistency with integer-order calculus by employing standard initial conditions. In contrast, the Riemann–Liouville transform requires fractional-order initial values of f(t), which are often physically non-intuitive. Consequently, the Caputo derivative is generally preferred for modeling physical systems and engineering phenomena.

3 The fractional hardening-damage constitutive model of cemented backfill material

The creep process of cemented backfill material is the result of the continuous adjustment of internal stress, hardening and damage effects and the joint action. In the step loading creep test of cemented backfill material, the hardening effect of the loading process is stronger due to the increase of deformation modulus, resulting in obvious hardening phenomenon; by the influence of instantaneous loading, the viscosity coefficient shows fluctuation, but the overall damage softening as the dominant.

Building upon these theoretical foundations, we enhance the classical generalized Kelvin model by incorporating both a hardening function and a fractional dashpot element. This modification enables the development of a novel fractional-order creep constitutive model that effectively captures the coupled creep-hardening-damage behavior characteristic of cemented backfill materials.

3.1 Classical generalized Kelvin model

More inclusive linear viscoelastic-plastic constitutive models can be created by combining the basic components, this subsection introduces the generalised Kelvin model, which consists of a spring and a Kelvin-Voigt model combined in series (see in Figure 1), and its mathematical representation is given by Equation 1. It has been found that describing the dynamic viscoelasticity of a material using the generalised Kelvin model is more flexible than the Kelvin-Voigt model [28].

The constitutive model of generalized Kelvin model is shown as follows:where E 1 , E 2 represent the Yang’s modulus (determined by the properties of the material), η denotes the Viscous coefficient, εand σ represent the strain (constantly positive number) and stress.

E 1 E 2 ε t + E 1 η ε ′ t = E 1 + E 2 σ t + η σ ′ t (1)

The creep deformation of the generalized Kelvin model under constant load consists of the sum of the two components, the spring and the Kelvin-Voigt model. As for the spring, it’s transient deformation equals to σ 0 E 1, and for the Kelvin-Voigt model, the creep equation is ε = σ 0 E 1 − e − E η t. Then the creep deformation of the generalized Kelvin model under constant load is:

ε t = σ 0 E 1 + σ 0 E 2 · 1 − e − E 2 η t

3.2 Modified generalized Kelvin model

Hardening of the deformation modulus occurs during instantaneous loading under step loading conditions, therefore, for the springs in the first part of the classical generalized Kelvin model, the hardening function is introduced [30]:where r and α are material constants, and 0 < r < 1,0<α < 1. We have:that is: ε 1 = σ E 1 σ , t = σ 1 − r E 1 t 1 − α.

E 1 σ , t = E 1 σ r t 1 − α

σ = E 1 σ , t ε 1

For deformation modulus E 1 σ , t with hardening function, assign values to some of the parameters, and observe the change of strain with time through the image, as shown in Figure 2.

The creep hardening phenomenon in cemented backfill materials manifests as a time-dependent increase in deformation modulus, with two key characteristics emerging from our analysis: First, the creep strain rate exhibits progressive deceleration over time, indicating material stiffening. Second, parametric studies reveal that: (i) For fixed values of other parameters, the creep deformation decreases monotonically with increasing hardening factor r (0 < r < 1); and (ii) The strain rate demonstrates a non-monotonic dependence on α, initially increasing before reaching a maximum and subsequently decreasing as α continues to increase.

Meanwhile, for the second part, we introduce the fractional order derivatives into the dashpot to construct Caputo fractional dashpot, then the modified generalized Kelvin model is shown in Figure 3.

The modified generalized Kelvin model is characterized by the following constitutive equation.where is the Caputo fractional dashpot, β is the fractional order satisfied 0<β < 1.

σ 1 = E 1 σ , t ε 1 t σ 2 = E 2 ε 2 t σ 3 = η · D β C ε 2 t σ = σ 1 = σ 2 + σ 3 ε = ε 1 + ε 2 (2)

Substituting Equation 2 gives:

D β C ε 2 t = σ η − E 2 η ε 2 t (3)

The Laplace transform of Caputo fractional derivative is obtained by applying both sides of Equation 3:where is the Laplace transformation function of.

s β ε ∼ 2 s = σ η s − 1 − E 2 η ε ∼ 2 s (4)

From Equation 4 we have:

ε ∼ 2 s = σ η s − 1 s β + E 2 η (5)

Combined with the properties of the Mittag–Leffler function [31], Equation 5 leads to:where the Mittag–Leffler function M a , b t = ∑ k = 0 ∞ t k Γ k a + b , a , b > 0.

ε 2 t = σ η · t β · M β , 1 + β − E 2 η t β

In summary, it can be obtained that the constitutive equation of the modified generalized Kelvin model is:

ε = ε 1 t + ε 2 t = σ 1 − r E 1 t 1 − α + σ η · t β · M β , 1 + β − E 2 η t β = σ 1 − r E 1 t 1 − α + σ η · t β · ∑ k = 0 ∞ − E 2 η t β k Γ 1 + β + k β (6)

Assigning values to some of the parameters in the model, a three-dimensional plot of the creep strain of the cemented backfill material as a function of time t and the order β of the fractional derivative can be obtained:

It can be seen from Figure 4 that when the parameter assignment is appropriate, the modified generalized Kelvin model can describe the harden-damage characteristics of cemented backfill materials in the creep process. When the loading stress level is low, the creep rate of the material continues to decrease, reflecting the hardening characteristics in the creep process, while under the action of high-level stress, the material will gradually fail and the creep rate will increase instantaneously.

3.3 Model verification and parameter determination

The creep test data under three typical stress levels from the C-8 creep test curve of cemented backfill material in literature [32] were fitted, and parameter identification was performed to validate the correctness and rationality of the proposed fractional hardening-damage creep constitutive model. The stress levels selected were 3.06, 4.08 and 5.10MPa, respectively. The ant colony algorithm was used for curve fitting and parameter identification of Equation 6, and the creep parameters obtained were shown in Table 1. The comparison between the fitting curve and the test curve was shown in Figure 5.

TABLE 1. Parameter fitting results of different confining pressures.

From the comparison curves in Figure 5, the curves are in good agreement. Therefore, the modified generalized Kelvin model established in this paper can better reflect the hardening-damage characteristics of the whole process of graded loading creep of cemented filling materials.

To further validate the superiority of the proposed model in capturing accelerated creep behavior, we compared its performance with the classical fractional Nishihara model under identical stress conditions (3.06, 4.08, and 5.10 MPa). As shown in Figure 6, while both models adequately describe primary and secondary creep stages, the modified generalized Kelvin model demonstrates significant advantage in simulating the accelerated creep phase under high stress (5.10 MPa). Specifically, the Nishihara model underestimates strain accumulation by 12%–18% in the accelerated regime due to its limited ability to coupling hardening and damage effects. cn contrast, our model incorporates a hardening function and fractional dashpot to explicitly represent the transition from hardening-dominated to damage-dominated behavior, resulting in a 94% correlation with experimental data during acceleration. This confirms the proposed model’s enhanced capability in predicting long-term failure thresholds in backfill materials.

3.4 Sensitivity analysis of model parameters

To evaluate the stability of the proposed model under parameter variations, we conducted numerical simulations by perturbing fractional order β, hardening factor r and α. The strain responses were compared to assess the model’s robustness.

The sensitivity analysis of key model parameters (r, α, β) under low-stress conditions (3.06 MPa) reveals distinct influences on creep strain evolution. Figure 7a shows a non-monotonic relationship with α, where intermediate values (α = 0.5) maximize strain accumulation (2.1%) before declining, suggesting an optimal range for damage-healing balance. Figure 7b demonstrates that increasing the hardening factor r (0.1→0.9) systematically reduces total creep strain by 42%–67%, confirming its role in suppressing time-dependent deformation through material stiffening. Figure 7c indicates that smaller fractional orders β (0.1–0.3) produce flatter creep curves characteristic of solid-like behavior, while higher β values (0.7–0.9) approach viscous fluid response with 2.8 increased strain rates. Collectively, these results validate the model’s ability to capture hardening-dominated regimes through parameter-dependent strain modulation, with r providing primary deformation control while α and β govern nonlinear temporal evolution. The absence of accelerated creep at this stress level confirms the threshold-dependent damage activation observed experimentally. The results in Figure 7d demonstrate that under high-stress conditions (5.10 MPa), increasing the fractional order β from 0.1 to 0.9 leads to a transition from primary creep to accelerated creep, with β ≥ 0.7 triggering a 3.2-fold increase in strain rate, clearly capturing the damage-induced instability characteristic of backfill materials near failure thresholds.

Physically, β governs the viscoelastic transition: lower β (0.1–0.3) produces solid-like behavior with minimal creep (Figure 7c), while higher β (0.7–0.9) approaches fluid-like response with accelerated deformation. Under high stress, increasing β triggers the transition from primary to tertiary creep (Figure 7d), capturing damage-induced instability. The hardening factor r suppresses creep strain by 42%–67% as it increases from 0.1 to 0.9 (Figure 7b), reflecting stress-induced stiffening. Parameter α exhibits a non-monotonic effect (Figure 7a), with intermediate values producing peak strain due to the competition between hardening and damage mechanisms.

The three-dimensional sensitivity analysis in Figure 8 comprehensively illustrates the coupled influence of parameters α and β on the creep strain evolution. For parameter α (Figure 8a), the surface plot reveals a pronounced nonlinear relationship, where intermediate α values (0.3–0.5) yield peak strain accumulation before declining at higher values, indicating an optimal range for damage-healing equilibrium. The β analysis (Figure 7b) demonstrates how increasing fractional order values progressively enhance strain rates and total deformation, with the steepest transitions occurring between β = 0.5–0.7. Notably, the curvature variations across both parameter spaces confirm their interdependent effects on viscoelastic response - while α primarily modulates the strain magnitude, β dominates the temporal evolution pattern. These 3D visualizations effectively capture the complex, nonlinear interactions governing the material’s time-dependent behavior under the combined hardening-damage mechanism.

3.5 Criteria for fractional order and parameter selection

• Physical Constraints: β∈(0,1) ensures the model interpolates between solid β = 0 and fluid β = 1 states, with values outside this range deemed physically inadmissible for viscoelastic materials.
• Numerical Stability: Truncation of the Mittag–Leffler series at 50 terms ensures convergence within typical time scales (t ∈ 0 , T).
• Experimental Calibration: Parameters are optimized via the ant colony algorithm to minimize error between predicted and experimental creep curves, with (β) prioritized to capture the transition from primary to accelerated creep.
• Sensitivity Guidance: Figures 6, 7 provide practical guidance on how α , β and r jointly influence strain evolution, aiding users in selecting values aligned with material responses.

4 Conclusion

The creep hardening behavior of cemented backfill materials plays a critical role in practical engineering applications. In this study, we developed a fractional-order modified generalized Kelvin model to comprehensively characterize the coupled creep hardening-damage effects in these materials. The proposed model integrates a hardening function and a fractional dashpot, enabling the analysis of time-dependent evolution in material strength under combined hardening and damage mechanisms.

Key contributions of this work include:

• Formulation of axial creep equations that account for elastic modulus variation under hardening-damage interactions, incorporating time-dependent effects.
• Experimental validation demonstrating the model’s accuracy in simulating creep behavior across both low- and high-stress regimes.
• Parametric studies revealing distinct hardening-dominated (decelerating creep) and damage-driven (accelerated creep) phases, aligned with observed material responses.

The model’s effectiveness was confirmed through curve fitting and parameter sensitivity analysis, supporting its applicability in predicting long-term stability of backfill structures.

From an engineering perspective, the proposed model offers a practical tool for predicting the long-term stability of cemented backfill structures. By identifying the transition from decelerating to accelerating creep through the evolution of β and the hardening factor r, the model can be used to estimate failure thresholds under high-stress conditions. This capability is particularly valuable for optimizing support system design and scheduling maintenance interventions in deep mining environments, where time-dependent deformation poses a significant risk to operational safety. The model thus provides a theoretical basis for more reliable stability assessments in cemented backfill engineering.

Future work could extend this framework to incorporate multiphysics interactions (e.g., hydration or seepage effects) for broader engineering relevance.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Shuang Shu, Hohai University, China

Reviewed by: Zhen-Chao Teng, Northeast Petroleum University, China

Rongbin Hou, North China University of Water Resources and Electric Power, China