Dynamics

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

The equation of motion for a typical dynamics problem has the following format: Mu¨+Cu˙+Ku=Fext\mathbf{M \ddot{u}}+\mathbf{C \dot{u}}+\mathbf{Ku}=\mathbf{F}_{ext}

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

Time integration

To solve the above equation for u\mathbf{u}, an appropriate time integration scheme needs to be chosen. Newmark (Newmark, 1959) and Hilber-Hughes-Taylor (HHT) (Hughes, 2000) 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+Δtt+\Delta t are written in terms of the displacement, velocity and acceleration at time tt and the displacement at t+Δtt+\Delta t using the NewmarkAccelAux and NewmarkVelAux Aux Kernels, respectively.

u¨(t+Δt)=u(t+Δt)−u(t)βΔt2−u˙(t)βΔt+β−0.5βu¨(t)u˙(t+Δt)=u˙(t)+(1−γ)Δtu¨(t)+γΔtu¨(t+Δt)\begin{aligned} \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{aligned}

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 (Au(t+Δt)=b\mathbf{Au}(t+\Delta t) = \mathbf{b}) from which u(t+Δt)\mathbf{u}(t+\Delta t) can be estimated.

  • For β=14\beta = \frac{1}{4} and γ=12\gamma = \frac{1}{2}, the Newmark time integration method is implicit and unconditionally stable. This is the constant average acceleration method with no numerical damping. This is recommended only when a constant timestep is used throughout the simulation. If for some reason, the simulation does not converge and the timestep is halved, this time integration method with no numerical damping can result in high frequency noise.

  • β=16\beta = \frac{1}{6} and γ=12\gamma = \frac{1}{2} results in the linear acceleration method where the acceleration is linearly varying between tt and t+Δtt+\Delta t.

  • β=0\beta = 0 and γ=12\gamma = \frac{1}{2} is identical to the central difference method.

commentnote

For 12≤γ≤2β\frac{1}{2} \le \gamma \le 2 \beta in structural dynamics problems, the Newmark method is unconditionally stable irrespective of the time-step Δt\Delta t. For γ=12\gamma = \frac{1}{2}, the Newmark method is at least second order accurate; it is first order accurate for all other values of γ\gamma.

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:

Mu¨(t+Δt)+C[(1+α)u˙(t+Δt)−αu˙(t)]+(1+α)Ku(t+Δt)−αKu(t)=Fext(t+(1+α)Δt)\mathbf{M}\mathbf{\ddot{u}}(t+\Delta t)+ \mathbf{C}[(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)

Here, α\alpha is the HHT parameter. For −13≤α≤0-\frac{1}{3} \le \alpha \le 0, β=(1−α)24\beta = \frac{(1-\alpha)^2}{4} and γ=12−α\gamma = \frac{1}{2} - \alpha, the HHT method is at least second-order accurate and unconditionally stable. For non-zero values of α\alpha, the response at high frequencies (above 12dt\frac{1}{2 dt}) are numerically damped, provided β\beta and γ\gamma are defined as above.

More details about the Newmark method and HHT method can be found in these lecture notes.

Rayleigh damping

This is the most common form of structural damping used in dynamic problems. Here, the damping matrix (C\mathbf{C}) is assumed to be a linear combination of the mass and stiffness matrices, i.e., C=ηM+ζK\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 is: Mu¨+(ηM+ζK)u˙+Ku=Fext\mathbf{M}\mathbf{\ddot{u}}+ (\eta M + \zeta K)\mathbf{\dot{u}} +\mathbf{K}\mathbf{u} = \mathbf{F_{ext}}

The degree of damping in the system depends on the coefficients ζ\zeta and η\eta as follows: ξ(ω)=η21ω+ζ2ω \xi (\omega) = \frac{\eta}{2} \frac{1}{\omega} + \frac{\zeta}{2} \omega(1)

where, ξ(ω)\xi(\omega) is the damping ratio of the system as a function of frequency ω\omega. For example, the damping ratio as a function of frequency for ζ=0.0035\zeta = 0.0035 and η=0.09\eta = 0.09 is presented in Figure 1. Note that ω\omega has units of rad/s and ff used in the figure has units of Hz.

Figure 1: Damping ratio as a function of frequency.

To model a constant damping ratio using Rayleigh damping, the aim is to find ζ\zeta and η\eta such that the ξ(f)\xi(f) is close to the target damping ratio ξt\xi_t, which is a constant value, between the frequency range [f1,f2][f_1, f_2]. This can be achieved by minimizing the difference between ξt\xi_t and ξ(f)\xi(f) for all the frequencies between f1f_1 and f2f_2, i.e., if

I=∫f1f2ξt−(η21ω+ζ2ω)dfI = \int_{f_1}^{f_2} \xi_t - \left(\frac{\eta}{2}\frac{1}{\omega} + \frac{\zeta}{2} \omega\right) df

Then, dIdη=0\frac{dI}{d \eta} = 0 and dIdζ=0\frac{dI}{d \zeta}=0 results in two equations that are linear in η\eta and ζ\zeta. Solving these two linear equations simultaneously gives:

ζ=ξt2π  3(Δf)2  (f1+f2−2f1f1Δf  lnf2f1)η=2πξt  f1f2Δf  [lnf2f1  (2+6f1f2(Δf)2)−3(f1+f2)Δf]\begin{aligned} \zeta &= \frac{\xi_t}{2 \pi} \; \frac{3}{(\Delta f)^2} \; \left(f_1 + f_2 - 2 \frac{f_1 f_1}{\Delta f} \; ln \frac{f_2}{f_1}\right) \\ \eta &= 2 \pi \xi_t \; \frac{f_1 f_2}{\Delta_f} \; \left[ln \frac{f_2}{f_1}\; \left(2 + 6 \frac{f_1 f_2}{(\Delta_f)^2}\right) - \frac{3(f_1 + f_2)}{\Delta_f}\right] \end{aligned}

A similar method can be used to estimate η\eta and ζ\zeta for non-constant damping ratios.

Implementation and Usage

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:

ρu¨(t+Δt)+ηρ[(1+α)u˙(t+Δt)−αu˙(t)]+ζ∇⋅[(1+α)ddtσ(t+Δt)−αddtσ(t)]+∇⋅[(1+α)σ(t+Δt)−ασ(t)]=Fext(t+(1+α)Δt)\rho \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}\sigma(t+\Delta t)- \alpha \frac{d}{dt}\sigma(t)] + \nabla \cdot [(1+\alpha) \sigma(t+\Delta t) - \alpha \sigma(t)] = F_{ext}(t+ (1 + \alpha) \Delta t)

Here, ρ\rho is the density of the material and σ\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 that contain ρ\rho are calculated in the InertialForce kernel. α\alpha, β\beta, γ\gamma and η\eta need to be passed as input to the InertialForce kernel.

The next two terms to the left involving σ\sigma are calculated in DynamicStressDivergenceTensors.

Note that the time derivative of σ(t+Δt)\sigma(t+\Delta t) and σ(t)\sigma(t) are approximated as follows:

ddtσ(t+Δt)=σ(t+Δt)−σ(t)Δtddtσ(t)=σ(t)−σ(t−Δt)Δt\begin{aligned} \frac{d}{dt}\sigma(t+\Delta t) &= \frac{\sigma(t+\Delta t)-\sigma(t)}{\Delta t} \\ \frac{d}{dt}\sigma(t) &= \frac{\sigma(t)-\sigma(t-\Delta t)}{\Delta t} \end{aligned}

The input file syntax for calculating the residual due to both the Inertial force and the DynamicStressDivergenceTensors is:
(contrib/moose/modules/solid_mechanics/test/tests/dynamics/rayleigh_damping/rayleigh_hht.i)

Here, ./DynamicSolidMechanics is the action that calls the DynamicStressDivergenceTensors kernel.

Finally, when using HHT time integration method, external forces like gravity and pressure also require α\alpha as input.

commentnote

For dynamic problems, it is recommended to use PresetDisplacement and PresetAcceleration for prescribing the displacement and acceleration at a boundary, respectively.

warningwarning

DynamicStressDivergenceTensors uses the undisplaced mesh by default but kernels such as InertialForce and Gravity, and boundary conditions such as Pressure use the displaced mesh by default. All calculations should be performed either using the undisplaced mesh (recommended for small strain problems) or the displaced mesh (recommended for finite strain problems). Therefore, use_displaced_mesh parameter should be appropriately set for all the kernels and boundary conditions.

Static Initialization

To initialize the system under a constant initial loading such as gravity, an initial static analysis can be conducted by turning off all the dynamics related Kernels and AuxKernels such as InertialForce, NewmarkVelAux and NewmarkAccelAux for the first time step. To turn off stiffness proportional Rayleigh damping for the first time step static_initialization flag can be set to true in DynamicSolidMechanics or DynamicStressDivergenceTensors. An example of static initialization can be found in this following test:

# One 3D element under ramped displacement loading.
#
# loading in z direction:
# time : 0.0 0.1  0.2  0.3
# disp : 0.0 0.0 -0.01 -0.01

# Gravity is applied in y direction. To equilibrate the system
# under gravity, a static analysis is run in the first time step
# by turning off the inertial terms. (see controls block and
# DynamicSolidMechanics block).

# Result: The displacement at the top node in the z direction should match
# the prescribed displacement. Also, the z acceleration should
# be two triangular pulses, one peaking at 0.1 and another peaking at
# 0.2.

# The y displacement would be offset by the gravity displacement.
# Also the y acceleration and velocity should be zero until the loading in
# the z direction starts (i.e, until 0.1s)

# Note: The time step used in the displacement data file should match
# the simulation time step (dt and dtmin in the Executioner block).

[Mesh]
  type = GeneratedMesh
  dim = 3 # Dimension of the mesh
  nx = 1 # Number of elements in the x direction
  ny = 1 # Number of elements in the y direction
  nz = 1 # Number of elements in the z direction
  xmin = 0.0
  xmax = 1
  ymin = 0.0
  ymax = 1
  zmin = 0.0
  zmax = 1
  allow_renumbering = false # So NodalVariableValue can index by id
[]

[Variables] # variables that are solved
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./disp_z]
  [../]
[]

[AuxVariables] # variables that are calculated for output
  [./accel_x]
  [../]
  [./vel_x]
  [../]
  [./accel_y]
  [../]
  [./vel_y]
  [../]
  [./accel_z]
  [../]
  [./vel_z]
  [../]
  [./stress_xx]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./strain_xx]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./strain_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./strain_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[Kernels]
  [./DynamicSolidMechanics] # zeta*K*vel + K * disp
    displacements = 'disp_x disp_y disp_z'
    stiffness_damping_coefficient = 0.000025
    static_initialization = true #turns off rayliegh damping for the first time step to stabilize system under gravity
  [../]
  [./inertia_x] # M*accel + eta*M*vel
    type = InertialForce
    variable = disp_x
    velocity = vel_x
    acceleration = accel_x
    beta = 0.25 # Newmark time integration
    gamma = 0.5 # Newmark time integration
    eta = 19.63
  [../]
  [./inertia_y]
    type = InertialForce
    variable = disp_y
    velocity = vel_y
    acceleration = accel_y
    beta = 0.25
    gamma = 0.5
    eta = 19.63
  [../]
  [./inertia_z]
    type = InertialForce
    variable = disp_z
    velocity = vel_z
    acceleration = accel_z
    beta = 0.25
    gamma = 0.5
    eta = 19.63
  [../]
  [./gravity]
    type = Gravity
    variable = disp_y
    value = -9.81
  [../]
[]

[AuxKernels]
  [./accel_x] # Calculates and stores acceleration at the end of time step
    type = NewmarkAccelAux
    variable = accel_x
    displacement = disp_x
    velocity = vel_x
    beta = 0.25
    execute_on = timestep_end
  [../]
  [./vel_x] # Calculates and stores velocity at the end of the time step
    type = NewmarkVelAux
    variable = vel_x
    acceleration = accel_x
    gamma = 0.5
    execute_on = timestep_end
  [../]
  [./accel_y]
    type = NewmarkAccelAux
    variable = accel_y
    displacement = disp_y
    velocity = vel_y
    beta = 0.25
    execute_on = timestep_end
  [../]
  [./vel_y]
    type = NewmarkVelAux
    variable = vel_y
    acceleration = accel_y
    gamma = 0.5
    execute_on = timestep_end
  [../]
  [./accel_z]
    type = NewmarkAccelAux
    variable = accel_z
    displacement = disp_z
    velocity = vel_z
    beta = 0.25
    execute_on = timestep_end
  [../]
  [./vel_z]
    type = NewmarkVelAux
    variable = vel_z
    acceleration = accel_z
    gamma = 0.5
    execute_on = timestep_end
  [../]
  [./stress_xx]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_xx
    index_i = 0
    index_j = 0
  [../]
  [./strain_xx]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = strain_xx
    index_i = 0
    index_j = 0
  [../]
  [./stress_yy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yy
    index_i = 1
    index_j = 1
  [../]
  [./strain_yy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = strain_yy
    index_i = 1
    index_j = 1
  [../]
  [./stress_zz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zz
    index_i = 2
    index_j = 2
  [../]
  [./strain_zz]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = strain_zz
    index_i = 2
    index_j = 2
  [../]
[]

[Functions]
  [./displacement_front]
    type = PiecewiseLinear
    data_file = 'displacement.csv'
    format = columns
  [../]
[]

[BCs]
  [./prescribed_displacement]
    type = PresetDisplacement
    variable = disp_z
    velocity = vel_z
    acceleration = accel_z
    beta = 0.25
    boundary = front
    function = displacement_front
  [../]
  [./anchor_x]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0.0
  [../]
  [./anchor_y]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  [../]
  [./anchor_z]
    type = DirichletBC
    variable = disp_z
    boundary = back
    value = 0.0
  [../]
[]

[Materials]
  [./elasticity_tensor]
    youngs_modulus = 325e6 #Pa
    poissons_ratio = 0.3
    type = ComputeIsotropicElasticityTensor
    block = 0
  [../]
  [./strain]
    #Computes the strain, assuming small strains
    type = ComputeSmallStrain
    block = 0
    displacements = 'disp_x disp_y disp_z'
  [../]
  [./stress]
    #Computes the stress, using linear elasticity
    type = ComputeLinearElasticStress
    block = 0
  [../]
  [./density]
    type = GenericConstantMaterial
    block = 0
    prop_names = density
    prop_values = 2000 #kg/m3
  [../]
[]

[Controls] # turns off inertial terms for the first time step
  [./period0]
    type = TimePeriod
    disable_objects = '*/vel_x */vel_y */vel_z */accel_x */accel_y */accel_z */inertia_x */inertia_y */inertia_z'
    start_time = 0.0
    end_time = 0.1 # dt used in the simulation
  [../]
[../]

[Executioner]
  type = Transient
  start_time = 0
  end_time = 3.0
  l_tol = 1e-6
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-6
  dt = 0.1
  timestep_tolerance = 1e-6
[]

[Postprocessors] # These quantites are printed to a csv file at every time step
  [./_dt]
    type = TimestepSize
  [../]
  [./accel_6x]
    type = NodalVariableValue
    nodeid = 6
    variable = accel_x
  [../]
  [./accel_6y]
    type = NodalVariableValue
    nodeid = 6
    variable = accel_y
  [../]
  [./accel_6z]
    type = NodalVariableValue
    nodeid = 6
    variable = accel_z
  [../]
  [./vel_6x]
    type = NodalVariableValue
    nodeid = 6
    variable = vel_x
  [../]
  [./vel_6y]
    type = NodalVariableValue
    nodeid = 6
    variable = vel_y
  [../]
  [./vel_6z]
    type = NodalVariableValue
    nodeid = 6
    variable = vel_z
  [../]
  [./disp_6x]
    type = NodalVariableValue
    nodeid = 6
    variable = disp_x
  [../]
  [./disp_6y]
    type = NodalVariableValue
    nodeid = 6
    variable = disp_y
  [../]
  [./disp_6z]
    type = NodalVariableValue
    nodeid = 6
    variable = disp_z
  [../]
[]

[Outputs]
  exodus = true
  perf_graph = true
[]
(contrib/moose/modules/solid_mechanics/test/tests/dynamics/prescribed_displacement/3D_QStatic_1_Ramped_Displacement_with_gravity.i)

References

  1. T. J. R Hughes. The finite element method:Linear static and dynamic finite element analysis. Dover Publications, Mineola, NY, 2000.[BibTeX]
  2. N. M. Newmark. A method of computation for structural dynamics. Journal of Engineering Mechanics, 85(EM3):67–94, 1959.[BibTeX]