# Linear Viscoelasticity – Part 2 – Solve the Equations

## Introduction

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.$$

(1)

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).$$

(2)

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$$

(3)

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).$$

(4)

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.$$

(5)

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.$$

(6)

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.$$

(7)

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)$$

(8)

\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}

(9)

## 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).

				
function strain(t)
t1 = maxStrain / strainRate
if t <= t1
return strainRate * t
else
return maxStrain - (t - t1) * strainRate
end
end
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]
end
return stress
end
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])
end
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.