Dynamic problems like response of a multi degree of freedom structure to external forcing and wave propagation in a medium can also be solved using the Tensor Mechanics module.

The equation of motion for a typical dynamics problem has the following format: $$$$\begin{eqnarray} \mathbf{M}\mathbf{\ddot{u}}+ \mathbf{C}\mathbf{\dot{u}} +\mathbf{K}\mathbf{u} = \mathbf{F_{ext}} \end{eqnarray}$$$$

Here, $$$\mathbf{M}$$$ is the mass matrix, $$$\mathbf{C}$$$ is the damping matrix, $$$\mathbf{K}$$$ is the stiffness matrix and $$$\mathbf{F_{ext}}$$$ is the vector of external forces acting at the nodes. $$$\mathbf{u}$$$, $$$\mathbf{\dot{u}}$$$ and $$$\mathbf{\ddot{u}}$$$ are the vector of displacement, velocity and acceleration at the nodes, respectively.

Rayleigh damping

This is the most common form of damping used in dynamic problems. Here, the damping matrix ($$$\mathbf{C}$$$) is assumed to be a linear combination of the mass and stiffness matrices, i.e., $$$\mathbf{C} = \eta \mathbf{M} +\zeta\mathbf{K}$$$. Here, $$$\eta$$$ and $$$\zeta$$$ are the mass and stiffness dependent Rayleigh damping parameters, respectively.

The equation of motion in the presence of Rayleigh damping becomes: $$$$\begin{eqnarray} \mathbf{M}\mathbf{\ddot{u}}+ (\eta M + \zeta K)\mathbf{\dot{u}} +\mathbf{K}\mathbf{u} = \mathbf{F_{ext}} \end{eqnarray}$$$$

Time Integration

To solve the above equation for $$$\mathbf{u}$$$, an appropriate time integration scheme needs to be chosen. Newmark and Hilber-Hughes-Taylor (HHT) time integration schemes are two of the commonly used methods in dynamics.

Newmark time integration

In Newmark time integration, the acceleration and velocity at $$$t+\Delta t$$$ are written in terms of the displacement, velocity and acceleration at time $$$t$$$ and the displacement at $$$t+\Delta t$$$.

$$$$\begin{eqnarray} \mathbf{\ddot{u}}(t+\Delta t) &=& \frac{\mathbf{u}(t+\Delta t)-\mathbf{u}(t)}{\beta \Delta t^2}- \frac{\mathbf{\dot{u}}(t)}{\beta \Delta t}+\frac{\beta -0.5}{\beta}\mathbf{\ddot{u}}(t) \\ \mathbf{\dot{u}}(t+ \Delta t) &=& \mathbf{\dot{u}}(t)+ (1-\gamma)\Delta t \mathbf{\ddot{u}}(t) + \gamma \Delta t \mathbf{\ddot{u}}(t+\Delta t) \end{eqnarray}$$$$

In the above equations, $$$\beta$$$ and $$$\gamma$$$ are Newmark time integration parameters. Substituting the above two equations into the equation of motion will result in a linear system of equations ($$$\mathbf{Au}(t+\Delta t) = \mathbf{b}$$$) from which $$$\mathbf{u}(t+\Delta t)$$$ can be estimated.

Hilber-Hughes-Taylor time integration

The HHT time integration scheme is built upon Newmark time integration method. Here, in addition to the Newmark equations, the equation of motion is also altered resulting in:

$$$$\begin{eqnarray} \mathbf{M}\mathbf{\ddot{u}}(t+\Delta t)+ (\eta\mathbf{M}+\zeta\mathbf{K})[(1+\alpha)\mathbf{\dot{u}}(t+\Delta t)-\alpha\mathbf{\dot{u}}(t)] +(1+\alpha)\mathbf{K}\mathbf{u}(t+\Delta t) - \alpha K \mathbf{u}(t) = \mathbf{F_{ext}}(t+ (1 + \alpha) \Delta t) \end{eqnarray}$$$$

Here, $$$\alpha$$$ is the HHT parameter.

Implementation in Moose

In the moose framework, the mass and stiffness matrices are not explicitly calculated. Only the residuals are calculated. To get the residual, the equation of motion (with Rayleigh damping and HHT time integration) can be written as: $$$$\begin{eqnarray} \rho\mathbf{\ddot{u}}(t+\Delta t) + \eta \rho [(1+\alpha)\dot{u}(t + \Delta t)-\alpha \dot{u}(t)] + \zeta \nabla \cdot [(1+\alpha)\frac{d}{dt}\boldsymbol{\sigma}(t+\Delta t)- \alpha \frac{d}{dt}\boldsymbol{\sigma}(t)] + \nabla \cdot [(1+\alpha) \boldsymbol{\sigma}(t+\Delta t) - \alpha \boldsymbol{\sigma}(t)] = \mathbf{F_{ext}}(t+ (1 + \alpha) \Delta t) \end{eqnarray}$$$$

Here, $$$\rho$$$ is the density of the material and $$$\boldsymbol{\sigma}$$$ is the stress tensor. The weak form of the above equation is used to get the residuals.

The first two terms to the left ($$$\rho\mathbf{\ddot{u}}(t+\Delta t) + \eta \rho [(1+\alpha)\dot{u}(t + \Delta t)-\alpha \dot{u}(t)]$$$)are addressed in InertialForce.C. $$$\alpha$$$, $$$\beta$$$, $$$\gamma$$$ and $$$\eta$$$ need to be passed as input to the InertialForce kernel. The input file syntax for calculating the residual due to Inertial force is:

     type = InertialForce
     variable = disp_x 
     velocity = vel_x
     acceleration = accel_x
     alpha = 0.11
     beta = 0.25
     gamma = 0.5
     eta = 0.1

The default values for $$$\alpha$$$, $$$\eta$$$ and $$$\zeta$$$ are zero.

The next two terms to the left involving $$$\boldsymbol{\sigma}$$$ are addressed in DynamicStressDivergenceTensors.C.

Note that the time derivative of $$$\boldsymbol{\sigma}(t+\Delta t)$$$ and $$$\boldsymbol{\sigma}(t)$$$ are approximated as follows: $$$$\begin{eqnarray} \frac{d}{dt}\boldsymbol{\sigma}(t+\Delta t) &=& \frac{\boldsymbol{\sigma}(t+\Delta t)-\boldsymbol{\sigma}(t)}{\Delta t} \\ \frac{d}{dt}\boldsymbol{\sigma}(t) &=& \frac{\boldsymbol{\sigma}(t)-\boldsymbol{\sigma}(t-\Delta t)}{\Delta t} \end{eqnarray}$$$$

The input file syntax for calculating the residual due to terms containing $$$\boldsymbol{\sigma}$$$ is:

     displacement = 'disp_x disp_y disp_z'
     zeta = 0.1
     alpha = 0.11  

Here, ./DynamicTensorMechanics is the action that calls DynamicStressDivergence.C

For calculating the time derivative of the stress, we require old and older states of stress, i.e., $$$\boldsymbol{\sigma}(t)$$$ and $$$\boldsymbol{\sigma}(t-\Delta t)$$$, respectively. These values need to be stored only when DynamicStressDivergence.C is used. Therefore, we can request ComputeStressBase.C to store the old states of stress through the boolean flag store_stress_old when DynamicStressDivergence.C is used. This boolean flag can be set to true in the input file through one of the child classes of ComputeStressBase.C

The input file syntax for setting store_stress_old to true is as follows:

      type = ComputeLinearElasticStress
      store_stress_old = true
      block = 0

Finally, external forces like gravity and pressure also require $$$\alpha$$$ as input. For example, the input file syntax for a Pressure boundary condition is as follows:

          boundary = bottom
          function = pressure_function
          disp_x = disp_x
          disp_y = disp_y
          disp_z = disp_z
          factor = 1
          alpha = 0.11

When the default values for $$$\eta$$$ and $$$\zeta$$$ are used, then Rayleigh damping is not included in the problem. Similarly, for the default value of $$$\alpha$$$, Newmark time integration method is used.

For examples on using Rayleigh damping and Time integration schemes, please check tensor_mechanics/tests/dynamics/rayleigh_damping and tensor_mechanics/tests/dynamics/time_integration, respectively.

For 1D wave propagation in a bar element, please check tensor_mechanics/tests/wave_1D


  1. Newmark time integration: Newmark, N. M. (1959) A method of computation for structural dynamics. Journal of Engineering Mechanics, ASCE, 85 (EM3) 67-94.

  2. HHT time integration: Hughes, T. J. R. (2000) The finite element method: Linear static and dynamic finite element analysis, Dover Publications, Mineola, NY.