Introduction
This is part 6 of my series on Linear Viscoelasticity. In this article I will show how the integral equation form of linear viscoelasticity (which I introduced in Part 1) is actually identical to a rheological model with parallel networks of linear springs and dashpots. There are 2 reasons why this discovery is important: (1) for most people it is easier to visualize how a spring-dashpot system behaves than an integral equation; (2) It provides a way to easily generalized linear viscoelasticity to non-linear viscoelasticity.
Maxwell Model Theory
To show that the integral form of linear viscoelasticity is the same a linear rheological model, I will start by analyzing a Maxwell model consisting of a linear spring in series with a linear dashpot. The total strain is given by: \(\varepsilon = \varepsilon_1 + \varepsilon_2\).
Recall that the definition of a linear spring is: \(\sigma = E \varepsilon_1\), and the definition of a linear dashpot is \(\dot{\varepsilon}_2 = \sigma / \eta\).
From these equations we can formulate force equilibrium in different ways. Here are 2 options (one in terms of the strain \(\varepsilon_2\), and one in terms of the stress):
\[\dot{\varepsilon}_2(t) = \displaystyle \frac{E}{\eta}\left[ \varepsilon(t) – \varepsilon_2 \right]\]
\[\dot{\sigma} + \displaystyle\frac{\sigma}{\eta/E} = E \dot{\varepsilon}\]
Strain Jump Response
We can now calculate the stress response from a Maxwell model loaded with a strain jump at time=0: \(\varepsilon(t) = \varepsilon_0\). By inserting this into the equilibrium equation based on stress we directly get:
\[\sigma(t) = \sigma_0 \exp\displaystyle\left[ \frac{-t}{\eta/E} \right]\]
From which we directly get the normlized relaxation modulus:
\[E_R(t) = E \exp\displaystyle\left[ \frac{-t}{\eta/E} \right]\]
This is identical to the expression that I derived from the integral formulation back in Part 1. Hence, a Maxwell spring-dashpot model behaves exactly like a one-term Prony series linear viscoelastic model! And a multiple Prony series model is equivalent to multiple parallel Maxwell elements.
Stress Calculation
In Part 2 of this series I showed how to calculate the stress response from the integral form of linear viscoelasticity. In this article I have shown how the stress equations can also be as a first-order ordinary differential equation. In the next two subsections i will examine forward and backward Euler numerical solution schemes for the stress.
Forward Euler Solution
The forward Euler solution to the stress is given by:
\[\sigma(t+\Delta t) = \displaystyle\left[ 1 – \frac{\Delta t}{\eta/E} \right] \sigma(t) + E \Delta\varepsilon\]
The figure below shows the load-unload stress response of a Maxwell model with E=1 MPa, \(\eta\)=1 MPa s, max strain = 0.1, and strain rate=0.01/s. In the figure I plotted the response for three cases that used 15, 25, and 55 time increments to solve the problem. As expected, the forward Euler method can become unstable if not enough time increments are used. The Julia code that I used to generate the image is shown below the figure.
using Plots
function calc_strain(time, maxStrain, strainRate)
t1 = maxStrain / strainRate
if time < t1
return strainRate * time
else
return maxStrain - (time - t1) * strainRate
end
end
function calc_stress(E, eta, maxStrain, strainRate, N)
time = range(0, 2*maxStrain/strainRate, length=N)
strain = zeros(N)
for i in 1:N
strain[i] = calc_strain(time[i], maxStrain, strainRate)
end
stress = zeros(N)
for i in 2:N
de = strain[i] - strain[i-1]
dt = time[i] - time[i-1]
stress[i] = (1 - dt/(eta/E)) * stress[i-1] + E * de
end
return time, strain, stress
end
time, strain, stress = calc_stress(1, 1, 0.1, 0.01, 55)
plot(strain, stress, label="Forward Euler (N=55)", linewidth=2)
time, strain, stress = calc_stress(1, 1, 0.1, 0.01, 25)
plot!(strain, stress, label="Forward Euler (N=25)", linewidth=2)
time, strain, stress = calc_stress(1, 1, 0.1, 0.01, 15)
plot!(strain, stress, label="Forward Euler (N=15)", linewidth=2)
plot!(size=(1000,800), framestyle=:box)
plot!(tickfontsize=14, guidefontsize=16, xlabel="Strain", ylabel="Stress")
savefig("Forward_Euler.png")
Backward Euler Solution
The backward Euler solution is given by:
\[\sigma(t+\Delta t) = \displaystyle\frac{1}{1 + \frac{\Delta t}{\eta/E}} \left[ \sigma(t) + E \Delta\varepsilon \right]\]
The figure below shows the load-unload stress response of a Maxwell model with the same properties as above. The figure shows the response for three cases that used 15, 25, and 55 time increments to solve the problem. The backward Euler method is more stable than the Forward Euler approach, but the number of time increments still have to be relatively large in order to get an accurate solution. The Julia code that I used to generate the image is shown below the figure.
using Plots
function calc_strain(time, maxStrain, strainRate)
t1 = maxStrain / strainRate
if time < t1
return strainRate * time
else
return maxStrain - (time - t1) * strainRate
end
end
function calc_stress(E, eta, maxStrain, strainRate, N)
time = range(0, 2*maxStrain/strainRate, length=N)
strain = zeros(N)
for i in 1:N
strain[i] = calc_strain(time[i], maxStrain, strainRate)
end
stress = zeros(N)
for i in 2:N
de = strain[i] - strain[i-1]
dt = time[i] - time[i-1]
stress[i] = (stress[i-1] + E*de) / (1 + dt/(eta/E))
end
return time, strain, stress
end
time, strain, stress = calc_stress(1, 1, 0.1, 0.01, 55)
plot(strain, stress, label="Backward Euler (N=55)", linewidth=2)
time, strain, stress = calc_stress(1, 1, 0.1, 0.01, 25)
plot!(strain, stress, label="Backward Euler (N=25)", linewidth=2)
time, strain, stress = calc_stress(1, 1, 0.1, 0.01, 15)
plot!(strain, stress, label="Backward Euler (N=15)", linewidth=2)
plot!(size=(1000,800), framestyle=:box)
plot!(tickfontsize=14, guidefontsize=16, xlabel="Strain", ylabel="Stress")
savefig("Backward_Euler.png")
Non-Linear Viscoelasticity
I started this series on linear viscoelasticity with the Bolzmann’s super position principle. I then showed how that directly leads to a integral equation formulation for the stress. One problem with this approach is that is it not easy to generalize it to a non-linear viscoelastic theory. One of the benefits of the rheological approach discussed in this article is that it is easy to extend that to a non-linear viscoelastic theory. One way to do this, for example, is to replace the linear dashpot with a non-linear dashpot. A common example is to use the following flow equation for the dashpot:
\[\dot{\varepsilon}_2 = \displaystyle\left(\frac{\sigma}{\tau_b}\right)^m\]
This flow equation is a simplified version of the Bergstrom-Boyce model that I developed for predicting the non-linear viscoelastic response of rubber-like materials.
Summary
- Linear viscoelasticity is often derived in terms of an integral equation, but is also equivalent to a spring-dashpot model.
- Non-linear viscoelasticity can be developed from a spring-dashpot rheological model
