# Mixed Finite Element Methods for Linear Viscoelasticity

## Introduction

One classical approach to linear isotropic elasticity is the displacement-based discretization for

Typically $u$ is discretized using continuous piecewise vector polynomials. This method yields accurate approximation for $u$. However, in practice

- The stress is usually the quantity of primary physical interest and pure displacement methods will yield stress approximations of lower order accuracy.
- The method performs poorly in the incompressible and nearly incompressible case, i.e., $\lambda \rightarrow \infty$.

An alternative approach, the stress-based formulation, addresses this problem by considering a mixed formulation. However, the main obstacle is to construct a stable pair of finite element spaces for symmetric tensors.

Instead of imposing the symmetry condition on the finite element space directly, we enforce the condition weakly by adding an additional equation and an associated Lagrange multiplier.

## Mathematical Formulation

Consider the linear viscoelastic Maxwell model

Here $\sigma$ is the stress tensor, and $A_0$, $A_1$ are fourth-order material tensors. We introduce the velocity vector

and the rotation of the velocity vector

Let $M$, $V$, and $K$ be the space of matrices, vectors, and skew symmetric matrices on $\Omega$, then the weak form for Eq. 1 is:

Find $(\sigma, v, \rho)\in H(\text{div}, \Omega; M) \times L^2(\Omega; V) \times L^2(\Omega; K)$, such that for all $(\tau, w, \eta)\in H(\text{div}, \Omega; M) \times L^2(\Omega; V) \times L^2(\Omega; K)$

Here $\Gamma_D$ is the Dirichlet boundary condition for the velocity $v$. Note that if we have traction boundary condition

The condition is part of the Dirichlet boundaries for $\sigma$.

After discretization, this leads to a DAE

where

This DAE can be solved using `TR_BDF2`

in ADCME.

## Linear Viscoelasticity Model

We consider the 2D Maxwell material model here

We will convert Eq. 2 to the form in Eq. 1. To this end, consider two linear operator

In the Voigt notation, these two linear operators have the matrix form

The following formulas can be easily derived:

Therefore, Eq. 2 can be rewritten as

which leads to

Here

In AdFem, $(a\sigma + b T\sigma, \tau)$ can be calculated using `compute_fem_bdm_mass_matrix`

.