Prony Series Fitting: A Step-by-Step Workflow
The Complete Pipeline
The workflow from raw experimental data to a validated Prony series follows six steps:
- Isothermal relaxation tests at 9–12 temperatures (e.g., −10 to 80°C), each covering 3–4 decades of log(t).
- TTS shifting to construct a master curve at reference temperature T0 (typically 20°C). Automated shifting via the CFS algorithm is preferred for objectivity.
- WLF equation fit to the shift factors → material constants C1, C2.
- Prony series fit to the master curve → coefficients {Gi, τi} and equilibrium modulus G∞.
- Validation: check r², compare predicted dynamic properties against DMA data.
- FEM input: enter Prony + WLF parameters into the structural model.
The Fitting Algorithm
Pre-select relaxation times
The relaxation times τi are not fitted — they are pre-selected as logarithmically spaced values spanning the master curve time range. For example, with N = 12 terms covering 10−4 s to 108 s: one term per decade. This ensures stable, unique solutions.
Linear least squares for the weights
With τi fixed, the problem is linear. Given M data points (tj, Gj) from the master curve:
Build matrix A: A[j, i] = exp(−tj / τi)
Target vector b: b[j] = Gj − G∞
Solve: min ‖A·x − b‖² subject to xi ≥ 0
Where xi = Gi (the Prony weights).
The non-negativity constraint is essential: negative weights would imply the material stiffens during relaxation at certain time scales — thermodynamically impossible and numerically unstable in FEM solvers.
Log-log space fitting
Fitting should be performed in log-log space (log G vs log t), not linear space. The modulus spans several orders of magnitude (0.1–1000 MPa), so linear fitting would weight the high-modulus glassy region disproportionately. Log-log fitting gives equal weight to each decade of both modulus and time.
Choosing the Number of Terms
Centelles et al. (2021) used 12–13 Prony terms and achieved r² > 0.997 for all seven materials tested. Ferry (1980, Ch. 3, p. 59) provides the theoretical justification:
“Any experimentally observed stress relaxation curve which decreases monotonically can in principle be fitted with any desired degree of accuracy to a series of terms [...] by taking n sufficiently large.”
| Application | N terms | Expected r² | Use case |
|---|---|---|---|
| Quick estimate | 2–5 | > 0.95 | Preliminary sizing, feasibility checks |
| Standard engineering | 8–12 | > 0.99 | EN 16612 design checks, building control submissions |
| Research / FEM | 12–20 | > 0.999 | Abaqus viscoelastic models, academic publications |
Quality Metrics and Validation
r² (coefficient of determination)
Target: r² > 0.99 for engineering, r² > 0.999 for research. All seven materials in Centelles et al. achieved r² > 0.997 with 12–13 terms.
Visual inspection
Always plot the master curve data alongside the Prony fit on a log-log scale. Check that the fit captures the glassy plateau, the transition slope, and the rubbery plateau. Residuals should show no systematic pattern.
Interconversion validation
The strongest validation: compute G′(ω) and G″(ω) from the Prony coefficients and compare with independent DMA frequency sweep data. Centelles et al. (2021) confirmed excellent agreement for all seven materials, proving the Prony series accurately predicts both time-domain and frequency-domain behaviour.
Practical Tips
Use all decimal places
Centelles et al. (2021, p. 5) explicitly warn: “All presented decimals of Prony series should be used in order to successfully input the data in the FEM software.” Rounding creates cumulative errors across the fit because the terms are interdependent.
Verify boundary values
After fitting, check: G(0) from the Prony series ≈ G0 from the data (instantaneous modulus). G(∞) from the Prony series ≈ G∞ from the data. G(t) must be monotonically decreasing for all t > 0.
Handle negative weights
If negative weights appear despite the non-negativity constraint, this indicates ill-conditioning: too many terms, poorly spaced τi, or insufficient data in certain time regions. Reduce N or adjust the τi spacing.
From Prony Coefficients to FEM Input
Normalised form (Abaqus, ANSYS)
Most FEM codes accept the normalised Prony series:
G(t) = G0 · [1 − ∑ gi(1 − e−t/τi)]
Input: G0 (instantaneous modulus) and pairs (gi, τi) for each term.
Temperature dependence
For temperature-dependent analysis, also provide the WLF constants C1, C2, and T0. The FEM solver computes the shift factor aT internally at each temperature using the TRS (Thermo-Rheological Simplicity) material model.
Galic et al. (2022) confirmed that Abaqus with Prony + WLF input produces results “in nearly 100% agreement with the experiment” for laminated glass beam deflections across the full temperature range −20°C to 50°C.
Automate the Workflow
Our Prony Calculator handles the master curve, WLF shifting, and least-squares fitting automatically. Select a material, choose the number of terms, and download the coefficients.
Launch Prony Calculator