Introduction
This is part 1 in my series on linear viscoelasticity. It turns out that the theory of linear viscoelasticity (LVE) is very simply and easy to derive from first principles. By presenting the thinking behind the theory it will be much easier to understand when linear viscoelasticity works, and when it will not be a good approach. In this first part I will focus on small strain uniaxial loading. In later parts I will extend the theory to other more general situations.
Derivation - 1D Loading
Linear viscoelasticity is a clean way to combine linear elasticity with viscoelasticity. It is based on the assumption shown in Figure1.
Figure 1. Ludwig Boltzmann (1844-1906) assumed that each loading step independently contributes to the final state.
Let’s start by considering the following thought experiment. Let’s apply a jump in strain (at time=0) and then hold that strain constant. At the same time measure the stress as a function of time. For most real polymers the response will look something like what is shown in Figure 2. We can then define the stress relaxation modulus from: \( E_r(t) = \sigma(t) / \varepsilon_0 \). As I will show below, the stress relaxation modulus is all that is needed to fully determine the stress for any strain history.
Figure 2. The stress will gradually relax after a jump in strain.
To make this more useful we need to be able to analyze general strain histories. The first step to do this is to write the applied strain as a sum of strain jumps:
\( \varepsilon(t) =\displaystyle \sum_{i=1}^\infty \Delta \varepsilon_i H(t – \tau_i) \)
Figure 3. The applied strain history can be approximated by a number of strain jumps.
By using Boltzmann’s superposition idea, the total stress simply becomes the sum of the stresses from each strain jump:
\( \sigma(t) =\displaystyle \sum_{i=1}^\infty \Delta \varepsilon_i E_R(t – \tau_i). \)
This equation can be converted to an integral since we are summing over an infinite number of small strain jumps:
\( \sigma(t) =\displaystyle \int_{0}^{t} E_R(t – \tau) d\varepsilon(\tau). \)
This equation can be converted to a time integral by dividing by and multiplying by \( d\tau \):
\( \sigma(t) =\displaystyle \int_0^{t} E_R(t – \tau) \dot{\varepsilon}(\tau) d\tau. \)
(1)
This master equation shows us how to calculate the stress from any uniaxial strain history. Pretty cool!
Example 1: Stress Relaxation
Let’s test the master linear viscoelasticity Equation (1) for a case with a strain jump at time 0 (as is shown to the left of Figure 2). To calculate the stress response we also need to have an equation for the stress relaxation modulus as a function of time. Assume the following form: \(E_r(t) = E_0 \exp(-t/\tau_0)\). This equation is particularly convenient since the integral of an exponential function is also an exponential function. If you perform the calculation you will get the following equation of the stress: \( \sigma(t) = \varepsilon_0 E_0 \exp(-t / \tau_0)\). The results are plotted in Figure 4.
Figure 4. Predicted stress relaxation response for different values of \(\tau_0\).
Example 2: Monotonic Tension
As a second example consider a linear viscoelastic material that is loaded with a constant strain rate in tension: \( \varepsilon(t) = \dot{\varepsilon}_0 t\). If I insert this into the master linear viscoelastic equation (1) I get the following stress expression:
\( \sigma =\displaystyle E_0 \dot{\varepsilon}_0 \tau_0 \left[ 1 – \exp\left( \frac{-\varepsilon}{\dot{\varepsilon} \tau_0} \right) \right]\).
(2)
At really large strains the stress reaches a steady state value of \(E_0 \dot{\varepsilon} \tau_0\), and the initial Young’s modulus is \(E_0\).
Figure 5. Predicted stress in monotonic tension.
Example 3: Arbitrary Strain History
The stress response due to an arbitrary strain history can be calculated using Equation (1). In a general case the stress cannot be expressed in closed form, but must be numerically calculated. One way to do this is to use the Julia code shown below. In this case I used the same stress relaxation equation as in Examples 1 and 2, but I apply a load-unload cycle. The strain rate was 0.01/s, and the max strain was 0.1. Note that the code presented here is perhaps the most inefficient way to calculate the stress. I simply solved Equation (1) numerically using the quadgk() function in julia. In an upcoming article I will show a significantly more efficient way to calculate the stress response.
Figure 6. Predicted load-unload response using the Julia code.
using QuadGK
using Plots
function strain(t)
maxStrain = 0.1
strainRate = 0.01
t1 = maxStrain / strainRate
if t <= t1
return strainRate * t
else
return maxStrain - (t - t1) * strainRate
end
end
function strainRate(t1)
dt = 1.0e-4
return (strain(t1+dt) - strain(t1)) / dt
end
function integrand(tau)
global t
E0 = 10
tau0 = 4.0
ER = E0 * exp(-(t-tau)/tau0)
return ER * strainRate(tau)
end
N = 100
time = range(0, 20, length=N)
ee = zeros(N)
ss = zeros(N)
for (i,tt) in enumerate(time)
global t
t = tt
ee[i] = strain(t)
val, err = quadgk(integrand, 0, t, rtol=1.0e-3)
ss[i] = val
end
plot(ee, ss, label="LVE Prediction", linewidth=2)
plot!(size=(1000,800), framestyle=:box)
plot!(tickfontsize=14, guidefontsize=16, xlabel="Strain", ylabel="Stress")
savefig("LVE1.png")
It is, of course, also possible to calculate the stress response for this specific case using MCalibration. Figure 7 shows that the MCalibration results are very similar to the results from Julia.
Figure 7. MCalibration stress prediction.
Summary
- It is easy to derive a 1D version of linear viscoelasticity from first principles.
- In the following articles in this series I will cover:
- Multiaxial loading
- Prony series
- A better computer implementation
- Cyclic loading (E’ and E’’)
- Differential equation formulation
