A Novel Hybrid Numerical Scheme for Solving Time-Fractional Viscoelastic Models in Structural Engineering: Application to Creep and Relaxation Behavior in Polymer Composites

This paper proposes a novel hybrid numerical scheme that augments the classical L1 finite-difference approximation of the Caputo fractional derivative of order α∈(0,1] with a selective shifted Grünwald–Letnikov correction (controlled by a shift parameter β∈[0,1)) applied only to the most recent time increment. When β=0, the scheme reduces exactly to the classical L1 scheme and retains its optimal convergence rate O(h2−α), where h denotes the uniform time-step size. For β>0 (optimally chosen as β=1−α/2), extra numerical damping is introduced at the cost of a mildly reduced convergence order O(h1−α), while long-term stability is significantly improved. The scheme is applied to the fractional Kelvin-Voigt and Standard Linear Solid models to analyze creep and relaxation responses. Numerical simulations demonstrate that the proposed hybrid scheme achieves improved accuracy, long-term stability, and computational efficiency compared to classical integer-order models and several existing fractional schemes reported in the recent literature. Results show that fractional orders capture anomalous creep behavior more accurately, aligning with experimental data from recent studies. The proposed method offers improved computational performance for real-time structural health monitoring applications.