Fractional-Order Viscoelastic Models

• Fractional-order viscoelastic models are defined by non-integer derivatives that capture power-law relaxation and anomalous creep in materials like biological tissues and polymers.
• They generalize classical models such as Kelvin–Voigt and Maxwell by incorporating springpot elements, which reduce parameter complexity while enhancing predictive accuracy.
• Efficient computational techniques like diffusive representations enable scalable numerical simulations for applications in biomechanics and material science.

Fractional-order viscoelastic models are a class of constitutive frameworks for soft materials in which stress–strain relationships are encoded through derivatives of non-integer (fractional) order in time. Unlike classical integer-order models based on combinations of springs and dashpots, fractional-order models efficiently capture observed power-law memory effects, broad relaxation spectra, and anomalous dissipation with a minimal number of physically interpretable parameters. Fractional viscoelasticity is of particular importance in modeling biological tissues, polymers, engineered soft solids, and complex fluids, where experimental characterizations reveal sub-exponential (non-Debye) stress relaxation, anomalous creep, and frequency-dependent scaling of both storage and loss moduli.

1. Mathematical Foundations and Canonical Models

The foundational element of fractional viscoelasticity is the Scott–Blair element or “springpot,” which linearly relates stress and strain through a fractional derivative of order $\alpha$ ($0<\alpha<1$):

$$\sigma(t) = c\,{}_{0}D_{t}^{\alpha} \varepsilon(t),$$

where ${}_{0}D_{t}^{\alpha}$ denotes the Caputo or Riemann–Liouville fractional derivative (commonly chosen based on initial condition requirements), and $c$ is a material stiffness coefficient carrying units $\mathrm{kPa}\cdot\mathrm{s}^{\alpha}$ (Pasupathy et al., 2023, Mainardi et al., 2011).

The two classical rheological models extended by fractional calculus are the Kelvin–Voigt and Maxwell models:

• Fractional Kelvin–Voigt (KV):

$$\sigma(t) = E\,\varepsilon(t) + \eta_{\alpha}\,{}^C D_t^{\alpha}\varepsilon(t)$$

The parallel arrangement of a spring and a springpot produces a stress–relaxation behaviour interpolating between pure elasticity and power-law creep.

• Fractional Maxwell (FM):

$$\sigma(t) + \tau^{\alpha}\,{}^C D_t^{\alpha} \sigma(t) = E\,\tau^{\alpha}\,{}^C D_t^{\alpha} \varepsilon(t)$$

This series configuration provides a power-law (algebraic) decay in relaxation modulus $G(t)$ and a transition from solid-like to fluid-like creep at long times.

Mitigating the limitations of purely springpot models, the fractional Zener or standard linear solid (SLS) model introduces augmented topologies to flexibly reproduce both instantaneous elastic and long-time dissipative regimes (Tolasa et al., 29 Nov 2025, Nguedjio et al., 16 Jun 2025). More general models, such as the “distributed-order Maxwell” and fractional Burgers/Jeffreys models, systematically assemble multiple fractional elements, allowing tailored fading-memory spectra and accurate fits across decades of time and frequency (Ferrás et al., 2022, Mainardi et al., 2011, Ruderman, 2020).

2. Constitutive Law Construction and Tensorial Generalizations

A general viscoelastic law using fractional derivatives in an isotropic or anisotropic continuum can be expressed using hereditary integrals:

$$\sigma_{ij}(t) = \frac{1}{\Gamma(1-\alpha)} \int_{0}^{t} \left[ (K_\alpha - \tfrac{2}{3}G_\alpha)\delta_{ij}\,\dot\varepsilon_{kk}(\tau) + 2G_\alpha\,\dot\varepsilon_{ij}(\tau) \right] (t-\tau)^{-\alpha} d\tau,$$

where $K_\alpha$ and $G_\alpha$ are fractional bulk and shear moduli, respectively, and the splitting into volumetric and deviatoric parts allows natural extension to three-dimensional tensorial settings (Pasupathy et al., 2023).

For inhomogeneous microstructures (e.g., axons in white matter), homogenized orthotropic fractional parameters $\{E_{ij}, \alpha_{ij}\}$ are extracted via numerical relaxation tests on finite element RVEs, enforcing periodic boundary conditions and solving an inverse problem to minimize deviation between simulated and calculated stress responses. Such upscaling bridges micro- to macro-scale viscoelasticity while fully representing the anisotropy and memory effects inherent to the composite (Pasupathy et al., 2023).

3. Parameter Identification and Model Fitting

Fractional-order parameters are typically fitted to experimental relaxation, creep, or oscillatory shear data by minimizing an objective function involving the storage ($G'$) and loss ($G''$) moduli in frequency space or direct time-domain kernel matching. For the springpot:

$$G(\omega) = c\,(i\omega)^{\alpha}, \quad G'(w) = c\,\omega^{\alpha} \cos \tfrac{\pi\alpha}{2}, \quad G''(w) = c\,\omega^{\alpha} \sin \tfrac{\pi\alpha}{2}$$

The cost function for parameter fitting can be formulated as: $$E = \frac{1}{2m} \sum_{k=1}^m \left[ (\ln c + \alpha \ln \omega_k - \Re G^*(\omega_k))^2 + \left(\frac{\pi}{2}\alpha - \Im G^*(\omega_k)\right)^2 \right]$$ Optimization over $(c, \alpha)$ yields concise, interpretable representations of tissue properties, often with fewer parameters and greater fidelity than integer-order models (Pasupathy et al., 2023, Nguedjio et al., 16 Jun 2025).

Other approaches incorporate human qualitative feedback (HiL Bayesian optimization) to identify perceptually optimal parameters, aggregating across subjects via committee machine fusion for population-level parameters relevant to biomechanics and haptic rendering (Tolasa et al., 29 Nov 2025).

4. Distributed-Order and Infinite-Order Formulations

Distributed-order fractional models generalize the classical springpot by introducing a weight function $c(\alpha)$ over $\alpha \in [0,1]$:

$$G(t) = \int_0^1 \frac{c(\alpha)}{\Gamma(1-\alpha)}\,t^{-\alpha} d\alpha$$

A delta function $c(\alpha) = \delta(\alpha-\beta)$ recovers the single-order case; a broad $c(\alpha)$ yields a continuum of relaxation mechanisms, physically representing broad or fractal memory spectra (Ferrás et al., 2022, Konjik et al., 2017). Series or parallel combinations (as in the generalised Maxwell model) with distributed-order elements allow flexible design of material responses, including effective approximations to Prony-series networks with infinitely many modes.

Infinite-order operators (e.g., Bessel models) use entire function generating functions for the differential operator, encoding an infinite discrete spectrum of relaxation times and yielding Mittag–Leffler or more elaborate time-domain responses. These models are thermodynamically admissible and can match observed viscoelastic responses without resorting to artificially large integer-order networks (Giusti, 2017).

5. Computational Methods and Numerical Implementation

Naïve discretizations of fractional derivatives have $O(N^2)$ time and memory cost due to nonlocal history summation. Modern fast solution methods represent the fractional kernel as a superposition of exponentials (“diffusive representation”), introducing a finite set of internal variables governed by ODEs for each quadrature node:

$${}^C D_t^\alpha y(t) \approx \sum_{j=1}^M w_j [-\partial_t z_j(t) + \text{other terms}],$$

with

$$\frac{d z_j}{d t} = y'(t) - \lambda_j z_j, \quad \frac{d Z_j}{d t} = y(t) - \lambda_j Z_j$$

This reduces complexity to $O(NM)$ or $O(N)$ for $M \ll N$. Truncating the history or employing short-memory principles yields further efficiency, enabling large-scale finite element simulations and real-time computation for applications such as haptic rendering and viscoelastic imaging (Diethelm, 2021, Jazia et al., 2013, Pasupathy et al., 2023). These approaches maintain accuracy while controlling computational resource requirements.

6. Thermodynamic and Frame-Indifference Considerations

Objectivity (material frame-indifference) is enforced by careful definition of the fractional operator acting on kinematic quantities (e.g., rate-of-deformation tensors). The models of Berjamin & Destrade (2024) provide three admissible classes for incompressible materials: upper-convected, lower-convected, and Capilnasiu–Palade types. These models are constructed to satisfy both objectivity under superposed rigid body motions and thermodynamic consistency (non-negative dissipation via explicit free-energy constructions) (Berjamin et al., 2023).

Criteria for admissibility can be systematically stated:

• The relaxation modulus $G(t)$ must be completely monotonic.
• The creep compliance $J(t)$ must be a Bernstein function.
• The Laplace transforms $G(s)/s$ and $sJ(s)$ are Stieltjes functions (analytic in $\mathbb{C} \setminus (-\infty,0]$, mapping the right half-plane to the left).

For complex-order derivatives, admissibility requires conjugate pairs and suitable weighting of terms to ensure energy dissipation for all sinusoidal inputs (Atanacković et al., 2014).

7. Applications, Extensions, and Practical Recommendations

Fractional-order models successfully fit soft tissue, brain white matter, wood, gels, and engineered soft polymers—any system in which empirical relaxation moduli exhibit extended power-law scaling spanning multiple decades in time or frequency (Pasupathy et al., 2023, Bonfanti et al., 2020, Nguedjio et al., 16 Jun 2025). Model selection and parameter tuning should consider the strain regime (linear or nonlinear), presence of anisotropy/orthotropy, and the measurement type (creep, relaxation, oscillatory).

Extensions to coupled visco-elasto-plasticity, memory-dependent damage (via Lemaitre-like evolution), and distributed-order models enable robust treatment of strong history effects, plasticity, and complex multi-scale relaxation phenomena (Suzuki et al., 2022, Suzuki et al., 2019). Fractional models are particularly powerful in simulating nonlinear dynamics of beams, microstructured solids, and as acousto-elastic calibrations for wave propagation (Suzuki et al., 2020, Konjik et al., 2017).

Numerically, fractional Zener and SLS models—often with $\alpha$ in $0.1-0.4$ for biological or polymeric solids—provide compact, robust, and interpretable alternatives to long Prony-series models, outperforming integer-order models in predictive fidelity for both viscoelastic and recovery regimes (Nguedjio et al., 16 Jun 2025, Mainardi et al., 2011). When higher computational efficiency is required, diffusive approximation with a small number of internal variables is recommended.

---

References

• (Pasupathy et al., 2023) A Fractional Viscoelastic Model Of The Axon In Brain White Matter
• (Ferrás et al., 2022) A generalized distributed-order Maxwell model
• (Tolasa et al., 29 Nov 2025) Active Learning of Fractional-Order Viscoelastic Model Parameters for Realistic Haptic Rendering
• (Berjamin et al., 2023) Models of fractional viscous stresses for incompressible materials
• (Nguedjio et al., 16 Jun 2025) Fractional order derivative approach of viscoelastic behavior of tropical wood
• (Giusti, 2017) On infinite order differential operators in fractional viscoelasticity
• (Diethelm, 2021) Fast Solution Methods for Fractional Differential Equations in the Modeling of Viscoelastic Materials
• (Suzuki et al., 2022) A General Return-Mapping Framework for Fractional Visco-Elasto-Plasticity
• (Suzuki et al., 2019) A Thermodynamically Consistent Fractional Visco-Elasto-Plastic Model with Memory-Dependent Damage for Anomalous Materials
• (Atanacković et al., 2014) Complex order fractional derivatives in viscoelasticity
• (Mainardi et al., 2011) Creep, Relaxation and Viscosity Properties for Basic Fractional Models in Rheology