Generalized Maxwell Viscoelastic Models#

The most general form of linear viscoelastic model is the generalized Maxwell model, which can be represented by a spring in parallel with a number of Maxwell models (see Fig. 5). With this model it is possible to represent a number of simpler viscoelastic models. For example, a simple Maxwell model corresponds to setting the elastic constants of all springs to zero, with the exception of the spring contained in the first Maxwell model (\(\mathit{\mu}_{1}\)). Similarly, the Kelvin-Voigt model corresponds to setting the elastic constants \(\mathit{\mu}_{2} = \mathit{\mu}_{3} = 0\), and setting \(\mathit{\mu}_{1} = \infty\) (or a very large number).

TODO

Add more information about using the generalized Maxwell model to implement various viscoelastic models (Burgers material).

We follow formulations similar to those used by [Zienkiewicz and Taylor, 2000] and [Taylor, 2003]. In this formulation, we specify the total shear modulus of the model (\(\mu_{tot}\)) and the bulk modulus (\(K\)). We then provide the fractional shear modulus for each Maxwell element spring in the model. It is not necessary to specify the fractional modulus for \(\mu_{0}\), because we can find it by subtracting the sum of the other ratios from 1. Note that the sum of all these fractions must equal 1. We use a similar formulation for our linear Maxwell viscoelastic model, but in that case \(\mu_{0}\) is always zero and we only use a single Maxwell model. The parameters defining the both materials are listed in Table 16.

As for all our viscoelastic models, the volumetric strain is completely elastic, and the viscoelastic deformation may be expressed purely in terms of the deviatoric components:

(52)#\[\boldsymbol{\sigma}^{dev} = 2\mu_{tot}\left[\mu_{0}\boldsymbol{\epsilon}^{dev} + \sum_{i=1}^{N}\mu_{i}\boldsymbol{q}^{i} - \boldsymbol{\epsilon}^{refdev}\right] + \boldsymbol{\sigma}^{refdev}; P = 3K(\theta - \theta^{ref}) + P^{ref},\]

where \(K\) is the bulk modulus, \(N\) is the number of Maxwell models, the terms with \(ref\) superscripts refer to reference states, and the variable \(q_{i}\) follows the evolution equations

(53)#\[\boldsymbol{\dot{q}}^{i} + \frac{1}{\tau_{i}}\boldsymbol{q}^{i} = \boldsymbol{\dot{\epsilon}}^{dev}.\]

The \(\tau_{i}\) are the relaxation times for each Maxwell model:

(54)#\[\tau_{i} = \frac{\eta_{i}}{\mu_{tot}\mu_{i}}.\]

An alternative to the differential equation form above is an integral equation form expressed in terms of the relaxation modulus function. This function is defined in terms of an idealized experiment in which, at time zero (\(t = 0\)), a specimen is subjected to a constant strain, \(\boldsymbol{\epsilon}^{dev}_{0}\), and the stress response, \(\boldsymbol{\sigma}^{dev}(t)\), is measured. For a linear material we obtain:

(55)#\[\boldsymbol{\sigma}(t) = 2\mu(t)(\boldsymbol{\epsilon_{0}^{dev} - \boldsymbol{\epsilon}^{refdev}}) + \boldsymbol{\sigma}^{refdev},\]

where \(\mu(t)\) is the shear relaxation modulus function. Using linearity and superposition for an arbitrary state of strain yields an integral equation:

(56)#\[\boldsymbol{\sigma}^{dev}(t) = \intop_{-\infty}^{t}\mu(t - T)\boldsymbol{\dot{\epsilon}}^{dev} dT.\]

Writing the modulus function in Prony series form we obtain

(57)#\[\mu(t) = \mu_{tot}\left(\mu_{0} + \sum_{i=1}^{N}\mu_{i}\exp\frac{-t}{\tau_{i}}\right),\]

where

(58)#\[\mu_{0} + \sum_{i=1}^{N}\mu_{i} = 1.\]

With the form in (57), the integral equation form is identical to the differential equation form.

If we assume the material is undisturbed until a strain is suddenly applied at time zero, we can divide the integral into

(59)#\[\intop_{-\infty}^{t} (\cdot) dT = \intop_{-\infty}^{0^{-}} (\cdot) dT + \intop_{0^{-}}^{0^{+}} (\cdot) dT + \intop_{0^{+}}^{t} (\cdot) dT .\]

The first term is zero, the second term includes a jump term associated with \(\boldsymbol{\epsilon}^{dev}_{0}\) at time zero, and the last term covers the subsequent history of strain. Applying this separation to (56),

(60)#\[\boldsymbol{\sigma}^{dev}(t) = 2\mu(t)(\boldsymbol{\epsilon}^{dev}_{0} - \boldsymbol{\epsilon}^{refdev}) + \boldsymbol{\sigma}^{refdev} + 2\intop_{0}^{t}\mu(t - T) \boldsymbol{\dot{\epsilon}}^{dev}(T) dT,\]

where we have left the sign off of the lower limit on the integral.

Substituting (57) into (60), we obtain

(61)#\[\boldsymbol{\sigma}^{dev}(t) = 2\mu_{tot}\left\{ \mu_{0}\boldsymbol{\epsilon}^{dev}(t) + \sum_{i=1}^{N}\left[\mu_{i}\exp\frac{-t}{\tau_{i}}\left(\boldsymbol{\epsilon}_{0}^{dev} + \intop_{0}^{t}\exp\frac{t}{\tau_{i}} \boldsymbol{\dot{\epsilon}}^dev(T) dT\right)\right] - \boldsymbol{\epsilon^{refdev}}\right\} + \boldsymbol{\sigma}^{refdev}.\]

We then split each integral into two ranges: from \(0\) to \(t_{n}\), and from \(t_{n}\) to \(t\), and define each integral as

(62)#\[\boldsymbol{i}_{i}^{1}(t) = \intop_{0}^{t}\exp\frac{T}{\tau_{i}}\boldsymbol{\dot{\epsilon}}^{dev}(T) dT.\]

The integral then becomes

(63)#\[\boldsymbol{i}_{i}^{1}(t) = \boldsymbol{i}_{i}^{1}(t_{n}) + \intop_{t_{n}}^{t}\exp\frac{T}{\tau_{i}}\boldsymbol{\dot{\epsilon}}^{dev}(T) dT.\]

Including the negative exponential multiplier, we have

(64)#\[\boldsymbol{h}_{i}^{1}(t) = \exp\frac{-t}{\tau_{i}}\boldsymbol{i}_{i}^{1}.\]

Then

(65)#\[\boldsymbol{h}_{i}^{1}(t) = \exp\frac{-\Delta t}{\tau_{i}}\boldsymbol{h}_{i}^{1}(t_{n}) + \Delta \boldsymbol{h}_{i}.\]

where

(66)#\[\Delta \boldsymbol{h}_{i} = \exp\frac{-t}{\tau_{i}}\intop_{t_{n}}^{t}\exp\frac{T}{\tau_{i}}\boldsymbol{\dot{\epsilon}}^{dev}(T) dT.\]

Approximating the strain rate as constant over each time step, the solution may be found as

(67)#\[\Delta \boldsymbol{h}_{i} = \frac{\tau_{i}}{\Delta t}\left(1 - \exp\frac{-\Delta t}{\tau_{i}}\right)\left(\boldsymbol{\epsilon}^{dev} - \boldsymbol{\epsilon}_{n}^{dev}\right) = \Delta h_{i}\left(\boldsymbol{\epsilon}^{dev} - \boldsymbol{\epsilon}_{n}^{dev}\right).\]

The approximation is singular for zero time steps, but a series expansion may be used for small time-step sizes:

(68)#\[\Delta h_{i}\approx1-\frac{1}{2}\left(\frac{\Delta t}{\tau_{i}}\right)+\frac{1}{3!}\left(\frac{\Delta t}{\tau_{i}}\right)^{2}-\frac{1}{4!}\left(\frac{\Delta t}{\tau_{i}}\right)^{3}+\cdots\,.\]

This converges with only a few terms. With this formulation, the constitutive relation now has the simple form:

(69)#\[\boldsymbol{\sigma}^{dev}(t) = 2\mu_{tot}\left( \mu_{0}\boldsymbol{\epsilon}^{dev}(t) + \sum_{i=1}^{N}\mu_{i}\boldsymbol{h}_{i}^{1}(t) - \boldsymbol{\epsilon^{refdev}}\right) + \boldsymbol{\sigma}^{refdev}.\]

We need to compute the tangent constitutive matrix when forming the Jacobian matrix. In addition to the volumetric contribution to the tangent constitutive matrix, we require the deviatoric part:

(70)#\[\frac{\partial\sigma_{ij}^{dev}}{\partial\epsilon_{kl}} = \frac{\partial\sigma_{ij}^{dev}}{\partial\epsilon_{mn}^{dev}}\frac{\partial\epsilon_{mn}^{dev}}{\partial\epsilon_{kl}},\]

where the second derivative on the right may be easily deduced from the definition of deviatoric strain in Mathematical notation for viscoelastic formulations.. The other derivative is given by

(71)#\[\frac{\partial\sigma_{ij}^{dev}}{\partial\epsilon_{mn}^{dev}} = 2\mu_{tot} \left[\mu_{0}\delta_{im}\delta_{jn} + \sum_{l=1}^{N}\mu_{l}\frac{\partial h_{lij}^{1}}{\partial \epsilon_{mn}^{dev}}\right].\]

From (65) and (67), the derivative inside the brackets is

(72)#\[\frac{\partial h_{lij}^{1}}{\partial\epsilon_{mn}^{dev}} = \Delta h_{l}(\Delta t)\delta_{im}\delta_{jn}.\]

The complete deviatoric tangent relation is then

(73)#\[\frac{\partial\sigma_{ij}^{dev}}{\partial\epsilon_{mn}} = 2\mu_{tot} \left[\mu_{0} + \sum_{l=1}^{N}\mu_{l}\Delta h_{l}(\Delta t)\right]\frac{\partial \epsilon_{ij}^{dev}}{\partial \epsilon_{mn}}.\]

We use this formulation for both our Maxwell and generalized Maxwell viscoelastic models. For the Maxwell model, \(\mu_{0} = 0\) and \(N = 1\). For the generalized Maxwell model, \(N = 3\). The stable time step is equal to 1/5 of the minimum relaxation time for all of the Maxwell models (54).