Linear Viscoelastic Materials
Viscoelastic materials have a time-dependent response even if the loading is constant in time. Many polymers and biological tissues exhibit this behavior. Linear viscoelasticity is a commonly used approximation where the stress depends linearly on the strain and its time derivatives (strain rate). Also, linear viscoelasticity deals with the additive decomposition of stresses and strains. It is usually assumed that the viscous part of the deformation is incompressible so that the volumetric deformation is purely elastic.
Generalized Maxwell Model
For isotropic linear elastic materials in the absence of inelastic stresses, Hooke’s law in Equation 3-11 reduces to
where the elastic strain tensor εel = ε − εinel represents the total strain minus initial and inelastic strains, such as thermal strains.
The stress tensor can be decomposed into a pressure and a deviatoric stress:
The pressure, mean stress, or volumetric stress, is given with a positive sign in compression
and the deviatoric stress is computed from the total stress minus the volumetric contribution
The elastic strain tensor εel can in the same way be decomposed into volumetric and deviatoric components
with the volumetric elastic strain given by
and the deviatoric contribution by
For isotropic linear elastic materials, the pressure is then related to the volumetric elastic deformation by the bulk modulus K
and the deviatoric stress tensor is linearly related to the deviatoric elastic strain tensor by the shear modulus G
The total stress in Equation 3-11 is then
In case of geometric nonlinearity, σ represents the second Piola-Kirchhoff stress tensor and εel the elastic Green-Lagrange strain tensor.
For viscoelastic materials, the deviatoric stress σd is not linearly related to the deviatoric strain εd but it also depends on the strain history. It is normally defined by the hereditary integral:
The function Γ(t) is called the relaxation shear modulus function (or just relaxation function) and it can be found by measuring the stress evolution in time when the material is held at a constant strain.
The relaxation function is often approximated by a Prony series:
A physical interpretation of this approach, often called the generalized Maxwell model, is shown in Figure 3-1
Figure 3-1: Generalized Maxwell model.
Hence, G is the stiffness of the main elastic branch, Gm represents the stiffness of the spring in branch m, and τm is the relaxation time constant of the spring-dashpot pair in branch m.
The auxiliary strain variable qm is introduced to represent the extension of the corresponding abstract spring, and the auxiliary variables γm = ε − qm represent the extensions in the dashpots.
The shear modulus of the elastic branch G is normally called the long-term shear modulus, or steady-state stiffness, and is often denoted with the symbol . The instantaneous shear modulus G0 is defined as the sum of the stiffness of all the branches
This is the stiffness when the external load is applied much faster than the shortest relaxation time of any viscous branch.
The relaxations time τm is normally measured in the frequency domain, so the viscosity of the dashpot is not a physical quantity but instead it is derived from stiffness and relaxation time measurements. The viscosity of each branch can be expressed in terms of the shear modulus and relaxation time as
The stress per branch can be written either in terms of the strain in the spring qm or the strain in the dashpot γm
The sum of the stresses in the viscoelastic branches is then computed from
The total stress in Hooke’s law (Equation 3-11) is then augmented by the viscoelastic stress σq
(3-17)
Computing the Stress on Each Branch
The auxiliary variable γm is a symmetric strain tensor, which has as many components as the number of strain components of the problem class. Since the stress per branch is written as
the auxiliary variables γm can be computed by solving the ODE
(3-18)
The relation between viscosity and relaxation time is
so that Equation 3-18 can equivalently be written as
(3-19)
The viscoelastic strain variables γm are treated as additional degrees of freedom. The shape functions are chosen to be one order lower than those used for the displacements because these variables add to the strains and stresses computed from displacement derivatives. The viscoelastic strain variables do not require continuity so discontinuous shape functions are used.
The viscoelastic strain variables γm are called solid.lemm1.vis1.ev, where solid is the Name of the physics interface node, and lemm1 is the name of the elastic material node.
Energy Dissipation
The dissipated energy density rate (SI unit: W/m3) in each branch m is
The rate of total dissipated energy density in the Generalized Maxwell material is then
In order to compute the dissipated energy density, the variable is integrated over time. For frequency domain studies, the dissipation of viscous forces averaged over a time period 2π/ω is computed from the shear loss modulus G’’ as
Standard Linear Solid Model
The standard linear solid model, also called SLS model, Zener model, or three-parameter model, is a simplification of the generalized Maxwell model with only one spring-dashpot branch:
Figure 3-2: Standard linear solid (SLS) model.
The stress in the single branch is computed as
where the relaxation time is related to the stiffness and relaxation time as η1 = τ1G1.
The auxiliary strain tensor γ1 is computed after solving by the ODE
and the dissipated energy density rate of the single branch is calculated from
Kelvin-Voigt Model
The Kelvin–Voigt viscoelastic model is represented by a spring connected in parallel with a damper:
Figure 3-3: Kelvin-Voigt model.
The stress tensor in the viscous branch is computed from the elastic strain rate
(3-20)
so there is no need to add the extra DOFs to compute the auxiliary strain tensor γ.
The relaxation time relates the viscosity and shear modulus by η = τG. The equivalent shear modulus is used in case of a anisotropic linear elastic material.
The dissipated energy density rate of the Kelvin-Voigt model is then computed from its rate
Burgers Model
The Burgers model, consists of a spring-dashpot branch in series with a Kelvin-Voigt branch:
Figure 3-4: Burgers model.
The strain in the fist dash-pot follows the ODE
where the shear modulus G is taken from the parent Linear Elastic material.
The strain in the second dash-pot follows the ODE
The total strain in the dampers is computed from
Combing these equations, it is possible to recover a second-order ODE for the strain tensor γ,
Note that Burgers viscoelastic material has two relaxation times related to the stiffness and viscosity in the springs and dash-pots. The relaxation times τ1 and τ2 are related to the stiffness and viscosities as τ1 = η1/G and τ2 = η2/G2.
Temperature Effects
For many polymers, the viscoelastic properties have a strong dependence on the temperature. A common assumption is that the material is thermorheologically simple (TRS). In a material of this class, a change in the temperature can be transformed directly into a change in the time scale. The reduced time is defined as
where αT(T) is a temperature-dependent shift function.
The implication is that the problem can be solved using the original material data, provided that the time is transformed into the reduced time.
Think of the shift function αT(T) as a multiplier to the viscosity in the dashpot in the Generalized Maxwell model. This shifts the relaxation time, so Equation 3-19 for a TRS material is modified to
For the SLS model, the shift applies to a single branch
and for the Kelvin-Voigt model, it applies to the viscosity in the damper.
Williams-Landel-Ferry Shift
One commonly used shift function is defined by the WLF (Williams-Landel-Ferry) equation:
where a base-10 logarithm is assumed. This shift is only valid over a certain range of temperature, typically around the glass transition temperature.
The first step to compute the shift factor αT consists of building a master curve based on experimental data. To do this, the curves of the viscoelastic properties (shear modulus, Young’s modulus, and so forth.) versus time or frequency are measured at a reference temperature T0. Then, the same properties are measured at different temperatures.
The shift value of each curve with respect to the master curve obtained at the temperature T0 defines the shift factor αT(T). The constants C1 and C2 are material dependent and are calculated after plotting log(αT) versus T − T0.
αT(T0) = 1 so that T0 is the temperature at which the master curve is given. If the temperature drops below T0 − C2, the WLF equation is no longer valid.
Since the master curve is measured at an arbitrary reference temperature T0, the shift factor αT(T) can be derived with respect to any temperature, and it is commonly taken as the shift with respect to the glass transition temperature. The values C1 = 17.4 and C2 = 51.6 K are reasonable approximations for many polymers at this reference temperature.
Arrhenius Shift
Below the Vicat softening temperature, the shift factor in polymers is normally assumed to follow an Arrhenius law. In this case, the shift factor is given by the equation
here, a base-e logarithm is assumed, Q is the activation energy (SI unit: : J/mol), and R is the universal gas constant.
Stationary Analysis
For stationary analysis it is possible to select either the long-term stiffness, in which case the stiffness of the viscoelastic branches is neglected, or the instantaneous stiffness, in which case the contribution from all branches is used.
The instantaneous shear modulus G0 is defined as the sum of the stiffness of all the branches
Frequency Domain Analysis and Damping
For frequency domain analysis, the frequency decomposition is performed as
Equation 3-17 and Equation 3-19 are then simplified to
where the shear storage modulus G and the shear loss modulus G” are defined for the generalized Maxwell model as
and
and for the SLS model as
and
and for the Kelvin-Voigt model as
and
The internal work of viscous forces averaged over a time period 2π/ω is computed as
See also the description of Viscoelasticity in the Solid Mechanics interface documentation.