[CM-002] Computational Modeling of Viscoelastic Materials: part 1

Constitutive Modeling of Viscoelasticity (VE)

Viscoelastic response is a mixed behavior of pure elastic and viscous deformation. Seen from different time scales, one response may be dominant over the other. Polymers are usually viscoelastic, owing to the internal macromolecular chains.

Computationally, many constitutive equations have been put forward for the modeling, such as the phenomenological models composed by springs and dashpots, denoting elastic and viscous responses, respectively, as shown in the figure below (a generalized Maxwell model). The computational model can be found at Ref [1].

In this two-part summary, here the infinitesimal strain model is demonstrated with the finite strain left to the upcoming one.

VE Constitutive Equation in Infinitesimal Strain

For linear isotropic material, the total energy can be decomposed by (volumetric and deviatoric):

$\\U^{0}(\Theta) = \frac{1}{2}K\Theta^2 \\\bar{W}^0(\mathbf{e}) = \mu \mathbf{e}:\mathbf{e}$

where the strain decomposition was also applied e = dev[ε] - tr[ε]1/3 and Θ = tr[ε] are deviatoric and volumetric strain in three-dimensional space, respectively. K and μ are bulk and shear modulus. Based on a perfect elastic response, the response of viscoelasticity can then be obtained by imposing viscous deformation on top of the initial elastic one and assume viscous deformation happens on the deviatoric part.

The total stress is decomposed into three parts: the initial elastic bulk stress, the initial elastic shear stress on the single spring element and the unsteady shear stress from all Maxwell elements:

$\mathbf{\sigma}_{n+1} = {U^{0}}^{'}(\Theta_{n+1}) \mathbf{1} + \gamma_{\infty} \mathbf{S}^{0}_{n+1} + \sum^{N}_{I=1} \gamma_I\mathbf{h}^{(I)}_{n+1}$

where γ_I is the shear modulus fraction of I-th spring, S⁰_n+1 is the initial elastic shear stress at time step n+1, and h^(I)_n+1 is the transient shear stress at time step n+1 in the I-th element, which can be calculated by:

$\\\mathbf{e}_{n+1} = \text{dev}[\mathbf{\varepsilon}_{n+1}]\\\mathbf{S}^{0}_{n+1} = \text{dev}[\partial_{\mathbf{e}} \bar{W}^0(\mathbf{e}_{n+1})]\\\mathbf{h}^{(I)}_{n+1}= \exp(-\Delta t_n/\lambda_I) \mathbf{h}^{(I)}_{n} + \exp(-\Delta t_n/2\lambda_I)(\mathbf{S}^{0}_{n+1} - \mathbf{S}^{0}_{n})$

The constitutive relation for linear loading (left panel a and c) and step loading (right panel b and d) is as follows (k=0.01; h1=0; h2=0.01; t1=0.5):

Beam Bending Example

With this constitutive equation for material response, using a numerical solver, e.g., finite element method (FEM), can solve the viscoelastic responses of structures. Here a simple cantilever viscoelastic beam is loaded with a tip force, with material parameters and geometry given below:

The duration of the applied load and the solved response of the tip displacement is:

Open-sourced code

If you are interested in how it works, please refer to the code [3] where all the “magic” are decoded explicitly.

References

[1] Juan C. Simo, and Thomas JR Hughes. Computational inelasticity. Vol. 7. Springer Science & Business Media, 2006.

[2] Strain stress decomposition [Wiki].

[3] The code