# System Level Matrices of Viscoelastic FEA Model

**Special thanks to Robert Kuether, org. 01556, for preparing this section.**

In Sierra/SD, the constitutive model for an isotropic linear viscoelastic
material uses a normalized Prony series to describe the time-dependent decay
from the glassy moduli to the rubbery moduli. Following the theoretical
development of the finite element formulation in the theory manual, the
element stiffness matrices may be cast as:


\begin{eqnarray}
K_{\nu,K} &=& \left( K_g - K_\infty \right)\int B^T D_K B dV \\
K_{\nu,G} &=& \left( G_g - G_\infty \right)\int B^T D_G B dV \\
K_e &=& K_\infty \int B^T D_K B dV + G_\infty \int B^T D_G B dV
\end{eqnarray}

The matrix $B$ is the strain-displacement matrix that depends on the element
shape function, while the scalar parameters $K_\infty$, $K_g$,
$K_\infty$ and $G_g$ represent the
rubbery (subscript $\infty$) and glassy (subscript $g$) bulk and shear moduli.
Both $D_K$ and $D_G$ are the constitutive matrices for the bulk and shear terms,
respectively. These element stiffness matrices (along with the element mass
matrix) are then assembled using standard finite element techniques,
resulting in the semi-discretized equations of motion for a structure with
linear viscoelastic materials.


$$
M \ddot{x} + \int_0^t K_{\nu,K} \zeta_K(t-\tau) \dot{x}(\tau) d\tau +
             \int_0^t K_{\nu,G} \zeta_G(t-\tau) \dot{x}(\tau) d\tau +
             K_e x = f(t)
$$ (eq:viscomatrix)

This coupled integro-differential equation contains real, symmetric $N$x$N$
system-level mass ($M$), viscoelastic bulk stiffness ($K_{\nu,K}$), viscoelastic bulk
shear ($K_{\nu,K}$), and elastic stiffness ($K_e$) matrices. The $N$x1 vectors $x$ and
$f(t)$ correspond to the physical displacements and externally applied forces,
respectively, and the dot represents the time derivative. The integral
terms have a simple functional form, such that the kernel functions are
a constant matrix multiplied by a series of normalized scalar
exponential functions (Prony series).

One can extract the system level matrices ($M$, $K_{\nu,K}$ , $K_{\nu,G}$, and $K_e$) directly from
Sierra/SD by writing out the matrices of an *isotropic linear elastic* FEA model.
The mass and stiffness matrices are written to MATLAB "\*.m" files when using the "dump"
solution type in the Sierra/SD input deck. The mass matrix extraction is
straightforward since it only depends on the density; however, extracting the
individual stiffness matrices is more complicated. A method for extracting
the system-level bulk and shear stiffness matrices using the dump solution
type is given in {numref}`tbl:viscomatrix`.

:::{table} Linear elastic material parameters to output system-level stiffness matrices using the dump solution type.
:name: tbl:viscomatrix

| **Output Matrix in {math:numref}`eq:viscomatrix`** | **Input Bulk Moduli** | **Input Shear Moduli** |
| :-------------------------------------------: | :-------------------: | :--------------------: |
|                     $K_e$                     |      $K_\infty$       |       $G_\infty$       |
|                  $K_{\nu,K}$                  |   $K_g - K_\infty$    |           0            |
|                  $K_{\nu,G}$                  |           0           |    $G_g - G_\infty$    |
:::

The input below shows how to extract the $K_{\nu,K}$ stiffness matrix.

:::{input}
SOLUTION
  case 'dump matrices'
        dump
END

FILE
     geometry_file 'plate_9by9inch.exo'
END

ECHO
     mass
END

BLOCK 1
     hex20
     material 1
END

//K_g = 9.8039e6
//K_inf = 7.0e6
//G_g = 3.7594e6
//G_inf = 2.5e6

MATERIAL 1
     Isotropic
     G= 1e-4     // essentially zero
     K= 2.8039e6     // = K_g - K_inf
     density=0.00024739
END
:::
