12. 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:

(12.1)\[\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.

(12.2)\[ 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) \]

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 Table 12.1.

Table 12.1 Linear elastic material parameters to output system-level stiffness matrices using the dump solution type.

Output Matrix in Eq. 12.2	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.

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