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.
\(\sigma(t+\Delta t) = E_0 \varepsilon(t+\Delta t) – \sigma_v(t+\Delta t)\)
(8)
(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).
using QuadGK, Plots
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.
