Linear Viscoelasticity
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, but it is also possible to account for volumetric viscoelasticity.
For isotropic linear elastic materials in the absence of inelastic stresses, Hooke’s law in Equation 3-18 reduces to
where the elastic strain tensor εel = ε − εinel represents the difference between the total strain 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 defined as
and the deviatoric contribution
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-18 is then
In case of geometric nonlinearity, σ represents the second Piola–Kirchhoff stress tensor and εel the elastic Green–Lagrange strain tensor.
The Generalized Maxwell Model
For deviatoric 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 relaxation times τm are 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 auxiliary strain variable qm is introduced to represent the extension of the corresponding abstract spring, and the auxiliary variables γm = εd − 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 it is often denoted with the symbol G. The instantaneous shear modulus G0 is then defined as the long-term shear modulus plus the sum of the stiffnesses of all the viscoelastic branches
This is the equivalent stiffness when an external load is applied much faster than the shortest relaxation time of any viscoelastic branch.
Sometimes, the relaxation function Γ(t) is expressed in terms of the instantaneous stiffness and relative weights, so that the Prony series is given as
In this case, the long-term shear modulus is related to the instantaneous shear modulus by the weight w < 1
and the shear moduli in each branch are defined by the weights wm
It must be assumed then that the weights fulfill the constraint
The long-term shear modulus G is given in the parent material model (for example, Linear Elastic Material). If your material data consists of the instantaneous shear modulus G0 and the weights wm, you can convert the data using the formulas above.
In the parent Linear Elastic Material, you can enter the long-term elastic data in any form, for example in terms of Young’s modulus, E, and Poisson’s ratio, ν. It will then be converted to the corresponding shear modulus G internally.
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-18) is then augmented by the viscoelastic stress σq
(3-29)
The same derivation is valid for a volumetric viscoelastic response, see Volumetric Response.
In the case of mixed volumetric and deviatoric viscoelasticity, it is assumed that the relaxation times per branch are the same for the deviatoric and volumetric deformation.
Computing the Stress in Each Branch
For deviatoric viscoelasticity, the auxiliary variable γm is a symmetric isochoric (volume preserving) 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-30)
The relation between viscosity and relaxation time is
so that Equation 3-30 can equivalently be written as
(3-31)
The viscoelastic strain variables γm are called for instance solid.lemm1.vis1.ev1_11, where solid is the tag of the physics interface, lemm1 is the tag of the linear elastic material, and vis1 is the tag of the viscoelasticity node. The trailing numbers indicate the branch number and the strain tensor indices.
Volumetric Response
It is usually assumed that the viscous part of the deformation is incompressible so that the volumetric deformation is purely elastic, but this does not have to be the case, and the same derivation as stated in previous sections can be applied to the volumetric response.
For volumetric viscoelasticity, the bulk modulus of the elastic branch K is normally called the long-term bulk modulus, or steady-state stiffness, and it is often denoted with the symbol K. The instantaneous bulk modulus K0 is then defined as the long-term bulk modulus plus the sum of the bulk moduli of all the viscoelastic branches.
The scalar equation for the volumetric deformation per branch, em, is
(3-32)
where τm is the relaxation time per branch and εel,vol is the elastic volumetric strain.
The sum of the stresses in the viscoelastic branches is then computed from
where Km is the bulk modulus of branch m.
The volumetric viscoelastic strain variables em are called, for example, solid.lemm1.vis1.evvol1, where solid is the tag of the physics interface, lemm1 is the tag of the linear elastic material node, and vis1 is the tag of the viscoelasticity node. The trailing number represents the branch number (that is, the first branch in this case).
Equation 3-31 and Equation 3-32 are solved using a Local Time Integration algorithm. Alternatively, the viscoelastic strain variables, γm and em, can be treated as additional degrees of freedom (DOFs), in which case is it possible to select Gauss point data or Discontinuous Lagrange shape functions (see Shape Functions and Element Types for details).
Using a Suitable Number of Branches
Only branches that have a relaxation time which is of the same order of magnitude as the temporal variation of the loading are important for the viscoelastic response. The cost in terms of memory consumption and computational speed increase significantly if the model contains many branches, particularly in time domain and for eigenfrequency analyses. In these cases many auxiliary equations of the form Equation 3-31 have to be solved for.
However, if a viscoelastic branch relaxes much faster or much slower than the time scale of the excitation period of external loads, the corresponding dashpot may be considered to be either fully relaxed or rigid, depending on the relaxation time of the branch relative to the excitation frequency. These fully relaxed or rigid branches can then be removed and grouped into equivalent branches, thus reducing the computational cost.
When the frequency band of interest is bounded by two cutoff frequencies, flower and fupper, it is possible to prune viscoelastic branches that have relaxation times such as
Branches with short relaxation times τk represent the response of the viscoelastic material to high frequency excitations. When the external loads represent excitations bounded by an upper frequency fupper, it is possible to prune branches with relaxation times that fulfill the relation τkfupper < 1. These branches are grouped into an equivalent branch, with relaxation time and stiffness such as
and
Viscoelastic branches with relaxation times τk such that 1 < τkflow are grouped into an equivalent branch forming the low frequency cutoff. The stiffness and relaxation time for the equivalent branch are computed from
and
Frequency Domain
In frequency domain, the relaxation function becomes complex-valued. The complex shear modulus for the generalized Maxwell model is then defined as the sum of the shear modulus in the pure elastic branch plus the complex shear moduli in the viscoelastic branches
The storage and loss moduli are then computed as the real and imaginary parts of the complex shear modulus
and
Energy Dissipation
The dissipated energy density rate (SI unit: W/m3) in each dashpot 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
The Maxwell Model
The Maxwell model (or Maxwell element) consists of a spring and a dashpot arranged in series. Maxwell elements are the building blocks used in The Generalized Maxwell Model.
Figure 3-2: Maxwell model.
In frequency domain, the compliance of the Maxwell model is defined as the sum of the compliance modulus of the elastic part plus the compliance in the damper
The complex shear modulus then reads
The storage and loss moduli are defined as the real and imaginary parts of the complex shear modulus
and
The Generalized Kelvin–Voigt Model
The generalized Kelvin–Voigt model is used to simulate the viscoelastic deformation in a wide range of materials such as concrete, biological tissues, and glassy polymers.
Just as for the generalized Maxwell model, the deviatoric strain εd is not linearly related to the deviatoric stress σd, but it also depends on the strain history. It is normally defined by the hereditary integral:
The function ψ(t) is called the compliance function (also called creep compliance or creep function) and it can be found by measuring the strain evolution in time when the material is held at a constant stress.
The compliance function is often approximated by a Prony series:
A physical interpretation of this rheological model consists of an elastic branch plus a number of Kelvin–Voigt elements arranged in series, this approach is shown in Figure 3-3.
Figure 3-3: The generalized Kelvin–Voigt model.
Here, G is the stiffness of the main elastic branch, Gm represents the stiffness of the spring in element m, and τm is the relaxation time constant of the spring-dashpot pair in branch m. The compliance J in the pure elastic branch is related to the shear modulus by
and the compliances per branch are related to the shear moduli as
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 in each branch can be expressed in terms of the stiffness or compliance modulus and relaxation time as
The auxiliary strain variables γm are introduced to represent the extension of the corresponding Kelvin–Voigt element. Since the elements are arranged in series, the total viscoelastic strain is given as the sum of the auxiliary strains
(3-33)
The deviatoric stress σd in the pure elastic branch is the same as in all the Kelvin–Voigt elements
(3-34)
Here, the deviatoric elastic strain is defined by the difference between the total and inelastic strains.
Computing the Strain in Each Branch
The auxiliary strain variables γm represent a symmetric isochoric (volume preserving) strain tensor, which has as many components as the number of strain components of the problem class.
The stress in each element m is given by the sum of the stresses in the spring and dashpot arranged in parallel
Here, γm is the strain in the element m, and ηm = τmGm is the viscosity. The auxiliary variables γm are computed by solving the ODE
(3-35)
The compliance of the elastic branch, J = 1/G, is normally called the instantaneous compliance, and it is often denoted with the symbol J0. This gives the equivalent stiffness when the material is loaded by an abrupt load much faster than the shortest relaxation time of any branch.
The long-term compliance J is defined as the instantaneous compliance plus the sum of the compliances of the viscoelastic branches
The long-term shear modulus then reads G = 1/J.
Sometimes, the compliance function ψ(t) is expressed in terms of the instantaneous compliance J and relative weights wm, so that the Prony series reads
In this case, the compliance and shear modulus in each branch are
and
and the long-term shear compliance and long-term shear modulus read
and
Frequency Domain
In frequency domain, the compliance function becomes complex valued. The complex compliance for the generalized Kelvin–Voigt model is then defined as the sum of the compliance in the pure elastic element (parent material) plus the complex compliances in the viscoelastic elements.
The storage and loss compliances are then computed as the real and imaginary parts of the complex compliance function
and
Energy Dissipation
The dissipated energy density rate (SI unit: W/m3) in each branch m is given by the dissipation in the dashpot
The rate of total dissipated energy density in the generalized Kelvin–Voigt 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
The Kelvin–Voigt Model
The Kelvin–Voigt viscoelastic model is represented by a spring connected in parallel with a dashpot:
Figure 3-4: The Kelvin–Voigt model.
The stress tensor in the viscous branch is computed from the elastic strain rate
(3-36)
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 an anisotropic linear elastic material.
The dissipated energy density rate of the Kelvin–Voigt model is then computed from its rate
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-5: Standard linear solid (SLS) model.
The stress in the viscoelastic 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 by solving by the ODE
and the dissipated energy density rate of the single branch is calculated from
The long-term shear modulus G = G is given in the parent Linear Elastic Material, and the instantaneous stiffness is given by the sum G0 = G + G1.
The Burgers Model
The Burgers model consists of a Maxwell (spring-dashpot) branch in series with a Kelvin–Voigt branch. The rheological representation of this material model is shown in Figure 3-6:
Figure 3-6: The Burgers model.
The strain in the first dashpot follows the ODE
where the shear modulus G is taken from the parent Linear Elastic material.
The strain in the second dashpot follows the ODE
The total strain in the dashpots is computed from
Combining these equations, it is possible to recover a second-order ODE for the strain tensor γ:
Note that the Burgers material has two relaxation times related to the stiffness and viscosity in the springs and dashpots. The relaxation times τ1 and τ2 are related to the stiffness and viscosities as τ1 = η1/G and τ2 = η2/G2.
The instantaneous stiffness G0 = G is given in the parent Linear Elastic Material.
Generalized Maxwell Model with Fractional Derivatives
Typically, the rheology of viscoelasticity models consists of springs and dashpots arranged in series and in parallel. Using the framework of fractional calculus, the constitutive equations of linear viscoelasticity can be generalized with a new type of element, named spring-pot (Ref. 50). In some cases, viscoelasticity models with fractional derivatives have shown to better match experimental data.
Figure 3-7: Spring-pot element.
The basic stress-strain relation of the spring-pot element is given by a proportionality factor p and a fractional time derivative of order β
The fractional order β takes a value between 0 < β < 1. For β = 0, the parameter p plays the role of a stiffness in a spring, and for β = 1 the parameter plays the role of the viscosity in a dashpot. The SI unit for such material parameter would be Pa·sβ.
For a pair of one spring and one spring-pot element connected in series, a so-called Maxwell element with a fractional time derivative, the deviatoric stress strain relation is given by
The rheological representation of this arrangement is shown in Figure 3-8
Figure 3-8: A spring-pot element connected in series with a spring.
In frequency domain, the fractional time derivative operator is substituted by
Thus, the stress-strain relationship reads
For a standard spring-dashpot branch, the relaxation time is given by the ratio of viscosity and stiffness, τ = η/G. Equivalently, the relaxation time τ for a branch with a spring and a spring-pot element in series is derived from
Thus, p = τβ/G, and the stress-strain relationship in the spring-spring-pot system reads
The complex-valued shear modulus for this Maxwell element reads
Storage and loss moduli are defined as the real and imaginary part of the complex shear modulus
and
The loss factor is computed from
For The Generalized Maxwell Model with fractional derivatives, all dashpots are replaced with spring-pot elements. The rheological representation of this material model is shown in Figure 3-9
Figure 3-9: The generalized Maxwell model with fractional derivatives.
The complex-valued shear modulus for the generalized Maxwell model with fractional derivatives reads
where Gm, τm, and βm are the shear modulus, relaxation time, and fractional order of branch m, respectively. Storage and loss moduli are defined as the real and imaginary part of the complex shear modulus.
and
Generalized Kelvin–Voigt Model with Fractional Derivatives
The Kelvin–Voigt viscoelastic model with fractional time derivative consists of a spring connected in parallel with a spring-pot element. The stress strain relation for such arrangement reads
The rheological representation of this arrangement is shown in Figure 3-10
Figure 3-10: Spring spring-pot elements connected in parallel.
In frequency domain, substitute the fractional time derivative operator by
The stress-strain relationship for this Kelvin element then reads
The relaxation time τ for a spring and a spring-pot element in parallel is derived from
so p = τβ/G, and the stress-strain relationship in the Kelvin–Voigt element reads
The complex-valued shear modulus in a spring-spring-pot arrangement reads
Storage and loss moduli are defined as the real and imaginary part of the complex shear modulus.
and
The loss factor is computed from
The complex-valued compliance of the Kelvin–Voigt element with fractional time derivative reads
For The Generalized Kelvin–Voigt Model with fractional derivatives, all dashpots are replaced with spring-pot elements.
The rheological representation of this material model is shown in Figure 3-11.
Figure 3-11: The generalized Kelvin–Voigt model with fractional derivatives.
The complex-valued compliance reads
where Gm, τm, and βm are the shear modulus, relaxation time, and fractional order of element m, respectively. Storage and loss moduli are defined as the real and imaginary part of the complex shear modulus
and
Standard Linear Solid Model with Fractional Derivatives
In the rheological representation of the Standard Linear Solid Model with fractional derivatives, the dashpot is replaced with a spring-pot element.
Figure 3-12: Standard linear solid (SLS) model with fractional derivatives.
The deviatoric stress in the spring-spring-pot branch is computed as
here, γ1 is the strain in the spring-pot element. In the frequency domain this equation translates to
Using the relaxation time τ = (p/G1)1/β, this reads
where Gq = G1(jωτ)β/(1+(jωτ)β) is the complex shear modulus of the spring-pot branch.
Subsequently, the shear modulus is for the SLS model with fractional derivatives reads
The storage and loss moduli are defined as the real and imaginary parts of the shear modulus GSLS, respectively
and
Burgers Model with Fractional Derivatives
In the rheological representation of The Burgers Model with fractional derivatives, the dashpots are replaced with a spring-pot elements. The rheological representation of this material model is shown in Figure 3-13:
Figure 3-13: The Burgers model with fractional derivatives.
The deviatoric stress in the main branch is computed as
here, γ1 is the strain in the spring-pot element. In frequency domain, and using the relaxation times τ1 = (p1/G)1/β1, this equation translates to
The stress in the second spring-pot follows
here, γ2 is the strain in the second spring-pot element. In frequency domain, and using the relaxation time τ2 = (p2/G2)1/β2, this equation reads
The total strain in the spring-pots is computed from
which combined with the relations for the strains in the spring-pots
and
gives the stress and strain relation for Burgers model with fractional derivatives
Subsequently, the compliance for the Burgers model with fractional derivatives reads
The storage and loss moduli are defined as the real and imaginary parts of the shear modulus GB = 1/JB, respectively
and
User Defined
With the user-defined viscoelastic material model it is possible to directly enter complex-valued expressions for the bulk and shear moduli or compliances, or for the loss factor. The complex-valued bulk and shear moduli are entered in terms of storage and loss moduli or compliances, respectively. The expressions can be functions taken directly from interpolated data, or can be analytical expressions of the frequency variable ω.
When the viscoelastic strain is deviatoric, the stress is computed from the shear storage modulus and shear loss modulus as
or from the storage and loss compliances as
When the viscoelastic strain is volumetric, the pressure is computed from the bulk storage modulus and bulk loss modulus
or from the storage and loss compliances as
When loss factor damping is selected, the stress-strain relation is augmented by a complex constitutive matrix
where D is the constitutive matrix computed from the parent material, and Dc is the complex constitutive matrix used to compute the viscoelastic stresses.
The internal variables for the frequency f and angular frequency ω are named phys.freq and phys.omega. Here, phys is the name of the physics interface, for instance solid.
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 variable 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 multiplier shifts the relaxation time, so Equation 3-31 for a TRS material is modified to
for all branches.
For the SLS model, the shift applies to the single relaxation time
and for the Kelvin–Voigt model, the shift multiplies the right hand side in Equation 3-36. In general, the shift function αT will multiply all the relaxation times or viscosities in a given viscoelasticity model.
Williams–Landel–Ferry Shift
One commonly used shift function is defined by the WLF (Williams–Landel–Ferry) equation:
where a base-10 logarithm is used. This shift is only valid over a certain temperature range, typically around the glass transition temperature.
The first step to compute the shift function αT(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 Tref. 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 Tref, defines the shift factor αT(T). The constants C1 and C2 are material dependent and are calculated after plotting log(αT) versus T − Tref.
The shift factor at the reference temperature equals αT(Tref) = 1, so that Tref is the temperature at which the master curve is given. If the temperature T drops below Tref − C2, the WLF equation is no longer valid.
Since the master curve is measured at an arbitrary reference temperature Tref, the shift function α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 used, Q is the activation energy (SI unit: J/mol), and R is the universal gas constant.
Tool–Narayanaswamy–Moynihan Shift
Structural relaxation in glass can be modeled using the so-called Tool–Narayanaswamy–Moynihan shift factor which is given as
here, a base-e logarithm is used, Q is the activation energy (SI unit: J/mol), R is the universal gas constant, T is the current temperature, Tref is a reference temperature, χ is a dimensionless activation energy fraction, and Tf is the so-called fictive temperature. The fictive temperature is given as the weighted average of partial fictive temperatures.
with
Here, wi are the weights and Tfi are the partial fictive temperatures. The partial fictive temperatures are determined from a system of coupled ordinary differential equations (ODEs) which follow Tool’s equation
here, λ0i is a structural relaxation time.
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 Studies, Equation 3-29 and Equation 3-31 are simplified to
where the shear storage modulus G' and the shear loss modulus G'' are user defined expressions or are defined for the generalized Maxwell model as
and
for the SLS model as
and
for the Kelvin–Voigt model as
and
The loss tangent is then computed from the ratio between the loss and storage moduli
where δ is the phase angle between deviatoric stress and strain.
The internal work of viscous forces averaged over a time period 2π/ω is computed as
Eigenfrequency Analysis
The direct use of frequency-dependent shear or bulk moduli, as described in Frequency Domain Analysis and Damping, leads to a nonlinear eigenvalue problem. COMSOL Multiphysics solves such nonlinear problems by expanding all nonlinear expressions with quadratic polynomials around a frequency linearization point which you can specify in the Eigenvalue Solver node. The linearized eigenvalue problem then reads
where ωL is the angular frequency at the linearization point. Thus, the solution of the linearized system of equations depends on the choice of the linearization frequency ωL.
Selecting the linearization frequency closer to one of the eigenfrequencies will produce a better result for that particular frequency, but not for all eigenfrequencies in the study. Hence, the eigenfrequencies need to be computed one by one using a certain iterative process that updates the linearization frequency ωL. This iterative procedure is needed when fractional derivatives are used, and also for the User Defined viscoelasticity model.
In order to avoid solving a nonlinear eigenvalue problem for the built-in viscoelasticity models, the same procedure as for time-dependent analysis is used; that is, each viscoelastic branch is represented by additional variables. For instance, for the The Generalized Maxwell Model in an eigenfrequency analysis, Equation 3-31 reads
The viscoelastic strain variables, γm, are treated as additional degrees of freedom (DOFs), and contribute to the damping and stiffness matrices of the coupled system of equations.
This procedure results in a damped linear eigenvalue problem that can be solved using the default eigenfrequency solver in a single run, for any specified number of eigenfrequencies.
For more information see the Eigenfrequency section in the Studies and Solvers chapter in the COMSOL Multiphysics Reference Manual.
An example of a nonlinear eigenvalue problem is shown in Eigenmodes of a Viscoelastic Structural Damper: Application Library path Structural_Mechanics_Module/Dynamics_and_Vibration/viscoelastic_damper_eigenmodes
Local Time Integration
For the The Generalized Maxwell Model and the Standard Linear Solid Model it is possible to use a local time integration algorithm for time-dependent analysis. By using this algorithm, the degrees-of-freedom related to viscoelasticity are treated as internal state variables, which makes the overall solution more efficient and also leaner in terms of memory usage. The implemented algorithm is based on the method suggested in Ref. 28, and it is schematically described in the following.
Consider the time evolution of a viscoelastic branch in the Generalized Maxwell model characterized by an ODE of the type
(3-37)
where γm represents the viscoelastic strain in the m branch. The exact solution to Equation 3-37 is then
(3-38)
By assuming that the strain ε(t) varies linearly within time increments, the integral in Equation 3-38 can be computed analytically, so the viscoelastic strain at increment n + 1 reads
(3-39)
All variables in Equation 3-39 are evaluated at the increment n and are stored as internal state variables. Equation 3-39 is solved for each branch in the Generalized Maxwell model, or for the single branch in the Standard Linear Solid model.
See also the description of Viscoelasticity in the Solid Mechanics interface documentation.