Viscoelastic Material

Constitutive Laws

This MPM material has separate constitutive laws for deviatoric stress and pressure.

Deviatoric Constitutive Law

The deviatoric constitutive law is always a small-strain, linear viscoelastic material with time-dependent shear modulus, G(t), which is given by a sum of n exponentials:

      [math]\displaystyle{ G(t) = G_0 + \sum_{i=1}^n G_i e^{-t/\tau_i} }[/math]

Here G₀ is the long-time shear modulus and the short-time shear modulus is the sum:

      [math]\displaystyle{ G(0) = \sum_{i=0}^n G_i }[/math]

The updates for components of the deviatoric stress become

      [math]\displaystyle{ ds_{ij} = 2\left( G(0) de_{ij} - \sum_{k=1}^n G_k d\alpha_{ij,k} \right) }[/math]

where α_ij,k are a series of internal variables that are tracked in history variables on each particle.

Pressure Constitutive Law

The pressure constitutive law has two options. The first in to use a small strain linear viscoelastic law time-dependent bulk modulus of

      [math]\displaystyle{ K(t) = K_0 + \sum_{i=1}^n K_i e^{-t/\tau_i} }[/math]

Here K₀ is the long-time bulk modulus and the short-time shear modulus is the sum:

      [math]\displaystyle{ K(0) = \sum_{i=0}^n K_i }[/math]

The pressure update becomes

      [math]\displaystyle{ dP = - K(0)(de_{ii} - 3d\varepsilon_{res}) + \sum_{k=1}^n K_k d\alpha_{V,k} }[/math]

where [math]\displaystyle{ de_{ii} }[/math] is increment volumetric strain and [math]\displaystyle{ d\varepsilon_{res} }[/math] accounts for thermal and solvent expansion effects (the thermal and solvent expansion coefficients are assumed to be independent of time). The [math]\displaystyle{ \alpha_{V,k} }[/math] are a series of internal variables that are tracked on each particle. To use this law, which is the default pressure law, set pressureLaw to 0 and enter all bulk moduli properties. For time-independent bulk modulus, enter K₀ but no relaxation times. The bulk modulus relaxation times can differ from the shear modulus relaxation times.

The second option is to use then Mie-Grüneisen equation of state (MGEOS). To use this law, set pressureLaw to 1 and enter the MGEOS properties. The pressure law models a non-linear elastic bulk modulus but cannot model time dependence of the bulk modulus.

Plane Stress Analysis

This material can be used in plane stress analysis, but only if it uses the linear pressure law (pressureLaw=0) and it does not add artificial viscosity. Support for plane stress in other conditions may be provided soon.

Effective Time Implementation

The above constitutive law assumes G(t) depends only on time, but in real materials, it will depend strongly on temperature. When modeling diffusion with solvent concentration, it may depend on solvent concentration as well. These dependencies are modeling by assuming the material obeys time-temperature and time-solvent superposition whereby G(t) is given by:

      [math]\displaystyle{ G(t) = G_0 + \sum_{i=1}^n G_i e^{-t/(a\tau_i)} }[/math]

where a is a shift factor. In other words, the coefficients of G(t) remain constant, but the relaxation times shift by a factor a. Furthermore, all relaxation times are assumed to shift by the same factor. We are led to define a reduced time:

      [math]\displaystyle{ t_r = \int_0^t \frac{dt}{a} \quad{\rm and}\quad \Delta t_r = \int_t^{t+\Delta t} \frac{dt}{a} }[/math]

is the effect time increment corresponding to actual time step [math]\displaystyle{ \Delta t }[/math]. The above constitutive law can then be extended to variable environments by converting convolution integrals to integrals over reduced or effective time. To implement this modeling, the material methods need input for a as a function of environment. When both temperature and solvent vary, a is replaced by a_Ta_m or product of thermal and moisture (i.e., solvent) shift factors. These shift factors are explained in the next two sections.

Variable Temperature

For temperature variations, a_T is the thermal shift factor. In polymer materials, this shift factor is approximated by the WLF equation:

      [math]\displaystyle{ \log a_T = \log \frac{\tau(T)}{\tau(T_{ref})} = - \frac{C_1(T-T_{ref})}{C_2+T-T_{ref}} }[/math]

This shift factor is implemented in simulations by entering Tref, C1, and C2. In polymer materials, if T_ref=T_g, or the glass transition temperature, C₁=17.44 and C₂=51.6 are average values over a range of polymers (the values for one specific polymer, however, may vary significantly). If T_ref differs from T_g, it can be used in a shifted WLF equation using C₃ and C₄ defined from C₁ and C₂ at T_g by

      [math]\displaystyle{ \log a_T = \log \frac{\tau(T)}{\tau(T_{ref})} = - \frac{C_3(T-T_{ref})}{C_4+T-T_{ref}} }[/math]

where

      [math]\displaystyle{ C_3 = \frac{C_1C_2}{C_2+T_{ref}-T_g} \quad{\rm and}\quad C_4 = C_2+T_{ref}-T_g }[/math]

Notice that [math]\displaystyle{ \log a_T\to\infty }[/math] as [math]\displaystyle{ T\to T_{ref}-C_2 }[/math] (e.g., if T reaches T_g-51.6 for an average polymer). This condition corresponds to infinite relaxation time and this material will respond as an elastic material for any temperature below this limit. Real materials may still have viscoelasticity effects, but those effects are not well modeled by extrapolating from T_g to far below T_g using the WLF equation.

The WLF equation works near and above T_g. An option for simulations well below T_g is to replace the non-Arrhenius response in the WLF equation with Arrhenius activation energy. The net result is an alternate shift factor:

      [math]\displaystyle{ \log a_T = \frac{\Delta H_a}{R\ln 10}\left(\frac{1}{T}-\frac{1}{T_{ref}}\right) }[/math]

where [math]\displaystyle{ \Delta H_a }[/math] is an apparent activation energy and R is the gas constant. This alternate shift factor is selected Tref and then setting C₁ to a negative number given by

      [math]\displaystyle{ C_1 = -\frac{\Delta H_a}{R} }[/math]

in units of degrees K.

Variable Solvent Concentration

For concentration variations, a_m is the solvent shift factor. In absence of a WLF equation, the solvent shift is modeled using a WLF-style equation based on assumption that solvent expansion leads to free volume that promotes relaxation:

      [math]\displaystyle{ \log a_m = \log \frac{\tau(m)}{\tau(m_{ref})} = - \frac{C_{m1}(m-m_{ref})}{C_{m2}+m-m_{ref}} }[/math]

This shift factor is implemented in simulations by entering mref, Cm1, and Cm2.

If [math]\displaystyle{ C_{m1} }[/math] and [math]\displaystyle{ C_{m1} }[/math] are known at one [math]\displaystyle{ m_{ref} }[/math], the shift relative to a new reference concentration, [math]\displaystyle{ m_{ref}' }[/math], would be:

      [math]\displaystyle{ \log a_m = \log \frac{\tau(m)}{\tau(m_{ref}')} = - \frac{C_{m1}'(m-m_{ref}')}{C_{m2}'+m-m_{ref}'} }[/math]

where

      [math]\displaystyle{ \quad{\rm where} \quad C_{m1}' = \frac{C_{m1}C_{m2}}{C_{m2}+m_{ref}'-m_{ref}} \quad{\rm and} \quad C_{m2}' = C_{m2} + m_{ref}'-m_{ref} }[/math]

For example, if [math]\displaystyle{ m_{ref}=c_{sat} }[/math], which is likely the condition with the shortest relaxation time, the shift factor would vary from [math]\displaystyle{ C_{m1}c_{sat}/(C_{m2}-c_{sat}) }[/math] to 0 between zero and saturation solvent conditions. But, if one switched to [math]\displaystyle{ m_{ref}'=0 }[/math], one would find [math]\displaystyle{ C_{m1}'= C_{m1}C_{m2}/(C_{m2}-c_{sat}) }[/math] and [math]\displaystyle{ C_{m2}'= C_{m2}-c_{sat} }[/math]. The shift factor would vary from 0 to [math]\displaystyle{ -C_{m1}c_{sat}/(C_{m2}-c_{sat}) }[/math] between zero and saturation solvent conditions (i.e., same relative change as when using [math]\displaystyle{ m_{ref}=c_{sat} }[/math]).

Simulations to include solvent effects on relaxations times must activate Diffusion Calculations, and enter material properties for saturation concentration, solvent expansion coefficient, and diffusion constants. Even constant-concentration simulations must activate diffusion calculations. To use constant concentration, set all particles to the same concentration (such that no diffusion occurs).

Elastic Mechanical Properties

The elastic properties of a viscoelastic material are the effective elastic properties at zero time. For example, the elastic shear modulus is

      [math]\displaystyle{ G(0) = G_e = G_0 + \sum_{i=1}^n G_i }[/math]

When simulations include a change in temperature or moisture, it is possible that both relaxations times and initial elastic properties will change. In time-temperature-superposition theory for time plots, these effects correspond to horizontal shifting (for changes in relaxation times) and vertical shifting (for changes in elastic properties). Vertical shifting is modeled in this material by setting all mechanical properties to

      [math]\displaystyle{ P(T,m) = b_Tb_mP(T_0,m_0) }[/math]

where [math]\displaystyle{ b_T }[/math] and [math]\displaystyle{ b_m }[/math] relative shifts of mechanical property [math]\displaystyle{ P(T,m) }[/math] as a function of temperature or moisture. The vertical shift factors are determined by entering any number of [math]\displaystyle{ (T,b_T) }[/math] and [math]\displaystyle{ (m,b_m) }[/math] pairs with monotonically increasing [math]\displaystyle{ T }[/math] and [math]\displaystyle{ m }[/math]. If simulation value is outside the entered range, it is found by linear extrapolation for first two or last two entered points. The reference temperature and moisture are arbitrary and equal to temperature where vertical shifts factors equal one. Note that vertical shifting can be used even when horizontal shifting is not used. Any material that [math]\displaystyle{ (T,b_T) }[/math] or [math]\displaystyle{ (m,b_m) }[/math] pairs will calculate vertical shift factors.

Material Properties

The unusual task for this material is to use multiple Gk and tauk properties (all with the same property name) to enter a material with multiple relaxation times.

Property	Description	Units	Default
G0	The long term (or fully-relaxed) shear modulus	pressure units	0
ntaus	The number of shear modulus relaxation times. This property is only needed in XML files and must come before any Gk or tauk properties. In scripted files, the number is automatically determined from the number of relaxation times you provide.	none	none
Gk	The shear modulus for the next relaxation time. Enter multiple Gk properties to have multiple relaxation times.	pressure units	none
tauk	The next shear modulus relaxation time. Enter multiple tauk properties to have multiple relaxation times.	time units	none
pressureLaw	Picks the constitutive law used for time independent pressure. The options are 0 to use linear viscoelastic law and 1 to use MGEOS equation of state.	none	0
K	Time-independent bulk modulus (when using linear viscoelastic law)	pressure units	none
ntausK	The number of bulk modulus relaxation times. This property is only needed in XML files and must come before any Kk or tauKk properties. In scripted files, the number is automatically determined from the number of relaxation times you provide.	none	none
Kk	The bulk modulus for the next relaxation time. Enter multiple Kk properties to have multiple relaxation times.	pressure units	none
tauKk	The next bulk modulus relaxation time. Enter multiple tauKk properties to have multiple relaxation times.	time units	none
(MGEOS)	Enter MGEOS properties C0, S1, S2, S3, gamma, and Kmax. The UJOption is fixed at 1.	varies	varies
Tref	Reference temperature to shift relaxation times using the WLF equation. If Tref<0, then no shifting is done and relaxation times will be independent of temperature.	K	-1
C1	Coefficient in WLF equation used when Tref≥0 to shift relaxation times. If the entered value is negative, the thermal shift switches to Arrhenious activation energy with ΔHa = -RC1,	none	17.44
C2	Coefficient in WLF equation used when Tref≥0 to shift relaxation times (not used if C1 is negative).	none	51.6
mref	Reference concentration to shift relaxation times using the WLF-style equation. If mref<0, then no shifting is done and relaxation times will be independent of concentration.	K	-1
Cm1	Coefficient in WLF-style equation used when mref≥0 to shift relaxation times	none	10
Cm2	Coefficient in WLF-style equation used when mref≥0 to shift relaxation times	none	0.0625
bTemp bTValue	Enter (bTemp,bTValue) pairs for piecewise interpolation of vertical temperature shift of elastic properties. Enter any number of pairs with bTemp values monotonically increasing	(degrees,none)	none
bConc bCValue	Enter (bConc,bCValue) pairs for piecewise interpolation of vertical concentration (or moisture) shift of elastic properties. Enter any number of pairs with bConc values monotonically increasing	none	none
alpha	Thermal expansion coefficient (ignored when using MGEOS law)	ppm/K	40
(other)	Properties common to all materials	varies	varies

The total number of Gk and tauk properies must be equal. In XML files, that total number must match the supplied ntaus property.

The default value for Kmax is -1, which means to not limit the bulk modulus. This mode is almost always stable, but simulations with high compression should always add the AdjustTimeStep Custom Task to keep calculation stable under high tangent bulk modulus conditions.

Viscoelastic Solids and Liquids

If G0 is not zero, the material is a viscoelastic solid, which means the shear stress at infinite time reamains a finite number. Viscoelastic solids are used to model materials such as elastomers that do not show long time flow due to their cross links or have a plateau shear modulus equal to G0.

If G0 is zero, the material is a viscoelastic liquid that will flow like a liquid if you wait long enough. For example, to emulate a liquid (i.e., similar to a Tait Liquid Material), set G0 to zero, use a single relaxation time with tauk short (on time scale of the simulation), and set the one Gk modulus to:

      [math]\displaystyle{ G_1 = {\eta\over 2\tau_1} }[/math]

where [math]\displaystyle{ \eta }[/math] is desired viscosity and [math]\displaystyle{ \tau_1 }[/math] is the single relaxation time.

History Variables

This material tracks internal history variables (one for each relaxation time and each component of stress) for implementation of linear viscoelastic properties, but currently none of these internal variables are available for archiving.

This material also tracks J (total relative volume change) and J_res (volume change of free expansion state) as history variables 1 and 2. Note that J_res is only needed, and therefore only tracked, when using MGEOS for pressure constitutive law (when pressureLaw is 1). If not tracked, it is always 1.