Linear Viscoelasticity – Part 2 – Solve the Equations


In Part 1 of this series I derived how to calculate the stress response for any strain history using the following integral equation:

\(\sigma(t) = \displaystyle\int_0^{t} E_R(t – \tau) \dot{\varepsilon}(\tau) d\tau.\)


As I also showed in Part 1, this equation can be solved using brute force in a high level math language like julia using a function like quadgk. The problem with this approach is that is VERY slow. To illustrate just how slow, let N be the number of floating point operations (FLOP) required to calculate \(\sigma(\Delta t)\) , then it takes MN floating point operations to calculate the stress \(\sigma(M \cdot\Delta t)\). And more importantly, it takes \( \mathcal{O}(N\cdot M^2)\) floating point operations to calculate a whole stress-strain curve consisting of M points. So, if a stress-strain curve has 1,000 data points then it will take about 1 million times longer to calculate the whole curve than the first data point. That is bad! We need a better approach. In this article I will show how you can perform this calculation 100X to 1000X faster.

One-Term Prony Series

Before calculating the stress response we need to define how we want to quantify the stress relaxation response. Here I will use the following 1-term Prony series equation:

\(E_R(t) = E_0 \cdot [(1-g) + ge^{-t/\tau_0}] \equiv E_0 \cdot g_R(t).\)


In this equation, \(g=0\) corresponds to a linear elastic response, and \(g=1\) corresponds to full stress relaxation.

Integrate by Parts

It is conceptually a good idea to integrate Equation (1) by parts using the following well-known equation:

\(\displaystyle\int_0^t u(\tau) \dot{v}(\tau) d\tau = u(t) v(t) – \int_0^t \dot{u}(\tau) v(\tau) d\tau\)


This gives the following equation for the stress:

\(\sigma(t) =\displaystyle E_0 \varepsilon(t) – \frac{E_0\, g}{\tau_0} \int_0^t e^{-(t-\tau)/\tau_0} \varepsilon(\tau) d\tau = \sigma_i(t) – \sigma_v(t).\)


What is cool about this equation is that we can see that the stress consists of an instantaneous part (\(\sigma_i\)) and a viscoelastic part \((\sigma_v\)). This equation has not helped us with the run-time issue, however. We need a different approach for that.

Divide the Integral Into 2 Parts

It is now time to apply some math “magic”. Start by rewriting Equation (4) for \(\sigma(t+\Delta t)\). Note that I divided the integral into two parts:

\(\sigma(t+\Delta t) =\displaystyle E_0 \varepsilon(t+\Delta t) – \frac{E_0\, g}{\tau_0} \int_0^t e^{-(t+\Delta t-\tau)/\tau_0} \varepsilon(\tau) d\tau – \frac{E_0\, g}{\tau_0} \int_t^{t+\Delta t} e^{-(t+\Delta t-\tau)/\tau_0} \varepsilon(\tau) d\tau.\)


If you think about it for a while you should be able to convince yourself that the first integral can be simplified a lot:

\(\sigma(t+\Delta t) =\displaystyle E_0 \varepsilon(t+\Delta t) – e^{-\Delta t/\tau_0} \cdot \sigma_v(t) – \frac{E_0\, g}{\tau_0} \int_t^{t+\Delta t} e^{-(t+\Delta t-\tau)/\tau_0} \varepsilon(\tau) d\tau.\)


To calculate the final integral we need  to assume that the strain varies linearly with time during the time interval [\(t, t+\Delta t\)]:

\(\varepsilon(\tau) = \left( \varepsilon_0 – \frac{\Delta \varepsilon}{\Delta t} \, t \right) + \frac{\Delta \varepsilon}{\Delta t} \tau.\)


With this assumption it is possible to evaluate the integral in Equation (7). After 2 pages of algebra (see Figure 1), I was able to get the following final equation for the stress.

Figure 1. My pen-and-paper derivation of the stress.

\(\sigma(t+\Delta t) = E_0 \varepsilon(t+\Delta t) – \sigma_v(t+\Delta t)\)


\begin{align}\sigma_v(t+\Delta t) &= e^{-\Delta t/\tau_0} \cdot \sigma_v(t) \\&+ E_0 g\cdot \left(\varepsilon_0 – \frac{\Delta \varepsilon}{\Delta t}\, t\right) \cdot\left[ 1 – e^{-\Delta t/\tau_0}\right]\\ & + E_0 g \, \frac{\Delta \varepsilon}{\Delta t} \left[ t + \Delta t -\tau_0 + (\tau_0 – t) \cdot e^{-\Delta t/\tau_0} \right]\end{align}


Calculate the Stress

The following Julia code calculates the stress for a given strain history, and Figure 2 shows the results from running the code. This implementation runs significantly faster than the original Equation (1).

					using QuadGK, Plots
function strain(t)
    t1 = maxStrain / strainRate
    if t <= t1
        return strainRate * t
        return maxStrain - (t - t1) * strainRate
function calcStress_direct(time, strain)
    N = length(time)
    stress = zeros(N)
    stressV = zeros(N)
    for i in 2:N
        dt = time[i] - time[i-1]
        de = strain[i] - strain[i-1]
        stressE = E0 * strain[i]
        stressV1 = stressV[i-1] * exp(-dt/tau0)
        stressV2 = E0*g * (strain[i-1] - de/dt * time[i]) * (1 - exp(-dt/tau0))
        stressV3 = E0*g * de/dt * (time[i]+dt-tau0 + (tau0-time[i])*exp(-dt/tau0))
        stressV[i] = stressV1 + stressV2 + stressV3
        stress[i] = stressE - stressV[i]
    return stress
E0 = 10.0
tau0 = 4.0
g = 1.0
maxStrain = 0.1
strainRate = 0.01
N = 100
time = range(0, 2*maxStrain/strainRate, length=N)
ee = zeros(N)
for i in 1:N
    ee[i] = strain(time[i])
ss3 = calcStress_direct(time, ee)
plot!(ee, ss3, label="direct", linewidth=2)
plot!(size=(1000,800), framestyle=:box)
plot!(tickfontsize=14, guidefontsize=16, xlabel="Strain", ylabel="Stress")

Figure 2. Predicted stress-strain response.


More to explore

Leave a Comment