Hilber-Hughes-Taylor (HHT) Time Integration Verification

Problem Statement and Results

In this one element test, the base and top of the 2D element is fixed in y direction. An acceleration is prescribed at the base of the element in the x direction.

HHT time integration is used to obtain the displacement, velocity and acceleration at the top of the element. The equation of motion is

$var element = document.getElementById("moose-equation-d74ef801-b53e-40fa-b3da-9079a3f8fa9c");katex.render("\\mathbf{M} a_x(t + \\delta t) + \\mathbf{K}((1+\\alpha)d_x(t+\\delta t) - \\alpha d_x(t)) = 0,", element, {displayMode:true,throwOnError:false});$ the $var element = document.getElementById("moose-equation-bc2d3e06-41ec-43ee-b6a1-9658e81121d7");katex.render("a_x", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-8b406a92-c557-49fe-a434-9cbebf9f3cf4");katex.render("d_x", element, {displayMode:false,throwOnError:false});$ are column vectors that contain the acceleration and displacement in the x-direction of the top and bottom node, respectively. The row of the mass matrix ( $var element = document.getElementById("moose-equation-4a5f61ba-21c7-43d8-bc29-e53f1657f638");katex.render("\\textbf{M}", element, {displayMode:false,throwOnError:false});$ ) that corresponds to the top node is computed as $var element = document.getElementById("moose-equation-49d4df76-dd7b-4c9c-a25e-98ae788c2a53");katex.render("\\rho \\cdot A \\cdot l \\cdot [1/3\\, 1/6]", element, {displayMode:false,throwOnError:false});$ , where $var element = document.getElementById("moose-equation-f72d416e-b9b0-4c25-ad20-1a053ca2de34");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$ is the density, $var element = document.getElementById("moose-equation-8b1d315c-1381-495e-830d-57225b0fa754");katex.render("A", element, {displayMode:false,throwOnError:false});$ is the area, and $var element = document.getElementById("moose-equation-108f0c0a-6939-4386-a204-b9f6af511ab7");katex.render("l", element, {displayMode:false,throwOnError:false});$ is the length. The row of the stiffness matrix that corresponds to the top node is $var element = document.getElementById("moose-equation-6f83f43c-b074-43ef-ba3d-43b6137d3645");katex.render("E\\cdot A\\cdot l \\cdot [1\\,-1]", element, {displayMode:false,throwOnError:false});$ , where $var element = document.getElementById("moose-equation-ef331aed-91fd-41a5-b49d-7d00737d0085");katex.render("E", element, {displayMode:false,throwOnError:false});$ is Young's modulus.

Solving the above equation analytically for $var element = document.getElementById("moose-equation-de00df89-cecb-4ea2-9e2d-e353421ab23a");katex.render("\\alpha = -0.3", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-08c5f9dd-b864-462f-96af-b5ef0d86b032");katex.render("\\beta = 0.4225", element, {displayMode:false,throwOnError:false});$ , and $var element = document.getElementById("moose-equation-b454d1ab-c62b-4f0e-93e6-6ed80dde043c");katex.render("\\gamma = 0.8", element, {displayMode:false,throwOnError:false});$ gives the results in Table 1. The interactive results shown in Figure 1 compare the analytical result ("previous") with the result computed from MASTODON for this problem ("current").

Table 1: Analytical displacement, acceleration, and velocity in the x-direction.

t	Displacement	Velocity	Acceleration
0	0	0	0
1	0.25	0.473	0.591
2	0.871\| 0.785	0.241
3	1.467	0.439	-0.493

Figure 1: Acceleration, velocity and displacement values computed using HHT integration in MASTODON

Complete Input File

[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 2
  nx = 1
  ny = 1
  xmin = 0.0
  xmax = 1.0
  ymin = 0.0
  ymax = 1.0
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
[]

[AuxVariables<<<{"href": "../../../../syntax/AuxVariables/index.html"}>>>]
  [./vel_x]
  [../]
  [./accel_x]
  [../]
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [./DynamicTensorMechanics<<<{"href": "../../../../syntax/Kernels/DynamicTensorMechanics/index.html"}>>>]
     displacements<<<{"description": "The nonlinear displacement variables for the problem"}>>> = 'disp_x disp_y'
     hht_alpha<<<{"description": "alpha parameter for mass dependent numerical damping induced by HHT time integration scheme"}>>> = -0.3
  [../]
  [./interia_x]
    type = InertialForce<<<{"description": "Calculates the residual for the inertial force ($M \\cdot acceleration$) and the contribution of mass dependent Rayleigh damping and HHT time  integration scheme ($\\eta \\cdot M \\cdot ((1+\\alpha)velq2-\\alpha \\cdot vel-old) $)", "href": "../../../../source/kernels/InertialForce.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    velocity<<<{"description": "velocity variable"}>>> = vel_x
    acceleration<<<{"description": "acceleration variable"}>>> = accel_x
    beta<<<{"description": "beta parameter for Newmark Time integration"}>>> = 0.4225
    gamma<<<{"description": "gamma parameter for Newmark Time integration"}>>> = 0.8
  [../]
[]

[AuxKernels<<<{"href": "../../../../syntax/AuxKernels/index.html"}>>>]
  [./accel_x]
    type = NewmarkAccelAux<<<{"description": "Computes the current acceleration using the Newmark method.", "href": "../../../../source/auxkernels/NewmarkAccelAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = accel_x
    displacement<<<{"description": "displacement variable"}>>> = disp_x
    velocity<<<{"description": "velocity variable"}>>> = vel_x
    beta<<<{"description": "beta parameter for Newmark method"}>>> = 0.4225
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  [../]
  [./vel_x]
    type = NewmarkVelAux<<<{"description": "Calculates the current velocity using Newmark method.", "href": "../../../../source/auxkernels/NewmarkVelAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = vel_x
    acceleration<<<{"description": "acceleration variable"}>>> = accel_x
    gamma<<<{"description": "gamma parameter for Newmark method"}>>> = 0.8
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  [../]
[]

[BCs<<<{"href": "../../../../syntax/BCs/index.html"}>>>]
  [./fixed_base_x]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    value<<<{"description": "Value of the BC"}>>> = 0.0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = bottom
  [../]
  [./fixed_base_y]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    value<<<{"description": "Value of the BC"}>>> = 0.0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = bottom
  [../]
  [./top_x]
    type = PresetAcceleration<<<{"description": "Prescribe acceleration on a given boundary in a given direction", "href": "../../../../source/bcs/PresetAcceleration.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = top
    function<<<{"description": "Function describing the velocity."}>>> = acceleration_top
    beta<<<{"description": "beta parameter for Newmark time integration."}>>> = 0.25
    acceleration<<<{"description": "The acceleration variable."}>>> = accel_x
    velocity<<<{"description": "The velocity variable."}>>> = vel_x
  [../]
  [./top_y]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = top
    value<<<{"description": "Value of the BC"}>>> = 0.0
  [../]
[]

[Functions<<<{"href": "../../../../syntax/Functions/index.html"}>>>]
  [./acceleration_top]
    type = PiecewiseLinear<<<{"description": "Linearly interpolates between pairs of x-y data", "href": "../../../../source/functions/PiecewiseLinear.html"}>>>
    x<<<{"description": "The abscissa values"}>>> = '0.0 1.0 2.0 3.0'
    y<<<{"description": "The ordinate values"}>>> = '0.0 1.0 0.0 -1.0'
  [../]
[../]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor<<<{"description": "Compute a constant isotropic elasticity tensor.", "href": "../../../../source/materials/ComputeIsotropicElasticityTensor.html"}>>>
    poissons_ratio<<<{"description": "Poisson's ratio for the material."}>>> = 0.3
    youngs_modulus<<<{"description": "Young's modulus of the material."}>>> = 1e9
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 0
  [../]
  [./strain]
    type = ComputeIncrementalSmallStrain<<<{"description": "Compute a strain increment and rotation increment for small strains.", "href": "../../../../source/materials/ComputeIncrementalStrain.html"}>>>
    displacements<<<{"description": "The displacements appropriate for the simulation geometry and coordinate system"}>>> = 'disp_x disp_y'
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 0
  [../]
  [./stress]
    type = ComputeFiniteStrainElasticStress<<<{"description": "Compute stress using elasticity for finite strains", "href": "../../../../source/materials/ComputeFiniteStrainElasticStress.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 0
  [../]
  [./density]
    type = GenericConstantMaterial<<<{"description": "Declares material properties based on names and values prescribed by input parameters.", "href": "../../../../source/materials/GenericConstantMaterial.html"}>>>
    prop_names<<<{"description": "The names of the properties this material will have"}>>> = 'density'
    prop_values<<<{"description": "The values associated with the named properties"}>>> = '2000.0'
  [../]
[../]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  start_time = 0.0
  end_time = 3.0
  dt = 1.0
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [./disp_x]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '1.0 1.0 0.0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = disp_x
  [../]
  [./vel_x]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '1.0 1.0 0.0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = vel_x
  [../]
  [./accel_x]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '1.0 1.0 0.0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = accel_x
  [../]
[]

[Outputs<<<{"href": "../../../../syntax/Mastodon/Outputs/index.html"}>>>]
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
[]

(test/tests/kernels/time_integration/hht.i)

# @requirement F1.2
# HHT time integration test

# In this one element test, the base and top of the 2D element is fixed in y
# direction. An acceleration is prescribed at the base of the element in the x direction.

# HHT time integration is used to obtain the displacement, velocity and acceleration
# at the top of the element.

# The equation of motion is

# [M] accel_x(t+ delta t) + [K]((1+alpha)disp_x(t+deltat) - alpha disp_x(t)) = 0

# The row of the mass matrix that corresponds to the top node is density * area
# * length * [1/3 1/6]. The row of the stiffness matrix that corresponds to the
# top node is Youngs's modulus * Area * l * [1 -1]. The accel_x and disp_x in
# above equation are column vectors which contain the x acceleration and
# displacement of the top and bottom node, respectively.

# Solving the above equation analytically for alpha = -0.3, beta = 0.4225 and
# gamma = 0.8 gives:

#     t        disp_x     vel_x       accel_x
#     0          0          0           0
#     1          0.25       0.473       0.591
#     2          0.871      0.785       0.241
#     3          1.467      0.439      -0.493

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 1
  ny = 1
  xmin = 0.0
  xmax = 1.0
  ymin = 0.0
  ymax = 1.0
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
[]

[AuxVariables]
  [./vel_x]
  [../]
  [./accel_x]
  [../]
[]

[Kernels]
  [./DynamicTensorMechanics]
     displacements = 'disp_x disp_y'
     hht_alpha = -0.3
  [../]
  [./interia_x]
    type = InertialForce
    variable = disp_x
    velocity = vel_x
    acceleration = accel_x
    beta = 0.4225
    gamma = 0.8
  [../]
[]

[AuxKernels]
  [./accel_x]
    type = NewmarkAccelAux
    variable = accel_x
    displacement = disp_x
    velocity = vel_x
    beta = 0.4225
    execute_on = timestep_end
  [../]
  [./vel_x]
    type = NewmarkVelAux
    variable = vel_x
    acceleration = accel_x
    gamma = 0.8
    execute_on = timestep_end
  [../]
[]

[BCs]
  [./fixed_base_x]
    type = DirichletBC
    value = 0.0
    variable = disp_x
    boundary = bottom
  [../]
  [./fixed_base_y]
    type = DirichletBC
    value = 0.0
    variable = disp_y
    boundary = bottom
  [../]
  [./top_x]
    type = PresetAcceleration
    variable = disp_x
    boundary = top
    function = acceleration_top
    beta = 0.25
    acceleration = accel_x
    velocity = vel_x
  [../]
  [./top_y]
    type = DirichletBC
    variable = disp_y
    boundary = top
    value = 0.0
  [../]
[]

[Functions]
  [./acceleration_top]
    type = PiecewiseLinear
    x = '0.0 1.0 2.0 3.0'
    y = '0.0 1.0 0.0 -1.0'
  [../]
[../]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = 0.3
    youngs_modulus = 1e9
    block = 0
  [../]
  [./strain]
    type = ComputeIncrementalSmallStrain
    displacements = 'disp_x disp_y'
    block = 0
  [../]
  [./stress]
    type = ComputeFiniteStrainElasticStress
    block = 0
  [../]
  [./density]
    type = GenericConstantMaterial
    prop_names = 'density'
    prop_values = '2000.0'
  [../]
[../]

[Executioner]
  type = Transient
  start_time = 0.0
  end_time = 3.0
  dt = 1.0
[]

[Postprocessors]
  [./disp_x]
    type = PointValue
    point = '1.0 1.0 0.0'
    variable = disp_x
  [../]
  [./vel_x]
    type = PointValue
    point = '1.0 1.0 0.0'
    variable = vel_x
  [../]
  [./accel_x]
    type = PointValue
    point = '1.0 1.0 0.0'
    variable = accel_x
  [../]
[]

[Outputs]
  csv = true
  exodus = true
[]

Problem Statement and Results
Complete Input File

Usage

Development

Demonstration

Hilber-Hughes-Taylor (HHT) Time Integration Verification

Problem Statement and Results

Complete Input File