## 2 Numerical Analysis of Creep Behavior
### Aging Kelvin Chain  
It is worth showing that the quasi-elastic constitutive law in Equations (8), (11), and (13) that has been just derived is fully equivalent to the solution of a rheological model consisting of a aging Kelvin chain with springs of stiffness \(E_{j}(t)\) in parallel with a dashpot with viscosities \(\tau_{j}E_{j}(t)\). To simplify the mathematical derivation, let us admit that in each time interval \(\Delta t=t_{i+1}-t_{i}\) of the loading history the spring stiffnesses is constant and equal to \(E_{j}(t^{*})\) where \(t_{i}<t^{*}<t_{i+1}\). Under this assumption the constitutive law of the aging Kelvin element in a time interval is given by the following first order differential equation  
\[\tau_{j}E_{j}(t^{*})\dot{e}_{j}+E_{j}(t^{*})\epsilon_{j}=\sigma \tag{15}\]  
An effective numerical integration of Equation (15) can be done by virtue of the so-called exponential algorithm, which makes possible the use of increasing time steps with the same accuracy and numerical stability [1]. The exponential algorithm assumes that the stress varies linearly in the time interval \(\Delta t=t_{i+1}-t_{i}\) so that the differential equation in (15) can be integrated exactly. The linear stress variation can be assumed as  
\[\sigma(t)=\sigma_{i}+(t-t_{i})\frac{\Delta\sigma}{\Delta t} \tag{16}\]  
where \(\Delta\sigma\) is the stress increment over the time step \(\Delta t\). With this stress variation the general solution of the differential equation in (15) can written as  
\[\epsilon(t)=A+B(t-t_{i})+C\left(1-e^{-\frac{t-t_{0}}{\tau}}\right) \tag{17}\]  
The integration constants \(A\) and \(B\) of the particular integral can be obtained from the substitution of Equation (17) into (15) while \(C\) is calculated from the initial condition in each time step, \(\epsilon=\epsilon(t_{i})\). This yields the strain, \(\epsilon(t_{i+1})\), at the end the time step \((t=t_{i+1}=t_{i}+\Delta t)\) as  
\[\epsilon_{j}(t_{i+1})=\epsilon_{j}(t_{i})e^{-\frac{\Delta t}{\gamma_{j}}}+ \frac{\sigma(t_{i})}{E_{j}(t^{*})}\left(1-e^{-\frac{\Delta t}{\gamma_{j}}} \right)+\frac{\Delta\sigma}{E_{j}(t^{*})}\left(1-\frac{\tau_{j}}{\Delta t}+ \frac{\tau_{j}}{\Delta t}e^{-\frac{\Delta t}{\gamma_{j}}}\right) \tag{18}\]  
So the strain increment is  
\[\Delta\epsilon_{j}=\epsilon_{j}(t_{i+1})-\epsilon_{j}(t_{i})=\left(\frac{ \sigma(t_{i})}{E_{j}(t^{*})}-\epsilon_{j}(t_{i})\right)\left(1-e^{-\frac{ \Delta t}{\gamma_{j}}}\right)+\frac{\Delta\sigma}{E_{j}(t^{*})}\left[1-\left( 1-e^{-\frac{\Delta t}{\gamma_{j}}}\right)\frac{\tau_{j}}{\Delta t}\right] \tag{19}\]Considering now a rheological model chain of \(N\) Kelvin elements in series with a spring with stiffness \(E_{0}(t^{*})\), the total strain increment is given by  
\[\Delta\epsilon=\Delta\epsilon_{0}+\sum_{j=1}^{N}\Delta\epsilon_{j} \tag{20}\]  
Substituting in this expression \(\Delta\epsilon_{j}\) from Equation (19) and \(\Delta\epsilon_{0}=\Delta\sigma/E_{0}\) one gets the stress increment in the time step from \(t_{i}\) to \(t_{i}+\Delta t\) as  
\[\Delta\sigma=E^{*}\Delta\epsilon-\sigma^{*}(t_{i}) \tag{21}\]  
where  
\[E^{*}=\left\{\frac{1}{E_{0}(t^{*})}+\sum_{j=1}^{N}\frac{1}{E_{j}(t^{*})}\left[ 1-\left(1-e^{-\frac{\Delta t}{\gamma}}\right)\frac{\tau_{j}}{\Delta t}\right] \right\}^{-1}\quad\text{and}\quad\sigma^{*}(t_{i})=E^{*}\sum_{j=1}^{N}\left(1-e ^{-\frac{\Delta t}{\gamma}}\right)\epsilon_{j}^{*}(t_{i}) \tag{22}\]  
in which \(\epsilon_{j}^{*}(t_{i})\) are the internal variables that can be updated according the recurrence formula in Equation (14). Since \(E^{*}=\widehat{E}\) and \(\sigma^{*}=\widehat{\sigma}\) the two approaches are fully equivalent. This means that the approximation of an aging compliance function through a Dirichlet (or Prony) series transforms the classical Volterra integral equation of creep (Equation (2)) into a rate-type formulation governed by an aging Kelvin chain with spring moduli \(E_{j}(t^{*})\) and viscosities \(\tau_{j}E_{j}(t^{*})\) obtained directly from the coefficient of the Dirichlet series approximation of the compliance function.