Skip to main content
`Ctrl`+`K`

# System Level Matrices of Viscoelastic FEA Model

# 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