Volumetric Locking Correction

Volumetric locking is the over stiffening of elements when the material is close to being incompressible (Poisson's ratio $var element = document.getElementById("moose-equation-22701582-278c-4557-a7ed-d5a87f8da740");katex.render("\\nu", element, {displayMode:false,throwOnError:false});$ nearing 0.5). This stiffening happens when a fully integrated element (such as Hex8 elements with 8 quadrature points or Quad4 elements with 4 quadrature points) is used. This is a numerical artifact introduced because shape functions used in finite element analysis cannot properly approximate the incompressibility condition throughout the element. To avoid this locking of elements, B-bar correction Hughes (1987) is implemented in MOOSE.

Theory

In this method, both the strain ( $var element = document.getElementById("moose-equation-76d51b1a-4645-4d4d-a11d-1415278784be");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ ) and the virtual strain ( $var element = document.getElementById("moose-equation-d670c883-f2c8-4ff5-8a65-242b1094cdd4");katex.render("\\delta \\epsilon", element, {displayMode:false,throwOnError:false});$ ) in an element are separated into volumetric and deviatoric components. The volumetric component is then replaced with an element averaged volumetric strain. This ensures that the volumetric strain remains constant throughout the element.

For example, in the case of small strain linear elasticity, the equation of motion is: $var element = document.getElementById("moose-equation-2e7b3469-881e-4388-93fd-1b86d09c20c7");katex.render("\\begin{aligned} \\int_V \\sigma(\\epsilon)\\delta \\epsilon dV - \\int_V b \\delta v dV - \\int_{\\partial V} t \\delta v dA = 0 \\end{aligned}", element, {displayMode:true,throwOnError:false});$ The element averaged volumetric strain (assuming small strain formulation) is: $var element = document.getElementById("moose-equation-5bfdcd38-5742-4990-b50a-64295be63679");katex.render("\\begin{aligned} w=\\frac{1}{V_e} \\int_{V_e} \\frac{1}{3} tr(\\epsilon) dV, \\end{aligned}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-5887cc54-c2e5-4ee6-9e09-3d2d746d6305");katex.render("V_e", element, {displayMode:false,throwOnError:false});$ is the volume of the element and tr(.) is the trace of the matrix.

The strain in each element is replaced by the approximation: $var element = document.getElementById("moose-equation-91b6a7e5-097e-4bcd-a2c1-369ef289a289");katex.render("\\begin{aligned} \\bar{\\epsilon} = \\epsilon +(w - \\frac{tr(\\epsilon)}{3})I \\end{aligned}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-951e2b50-b752-48b6-b1d3-c93f21832d3c");katex.render("I", element, {displayMode:false,throwOnError:false});$ is the $var element = document.getElementById("moose-equation-340d2e83-49dd-4905-ac91-668ae42e4dcb");katex.render("3 \\times 3", element, {displayMode:false,throwOnError:false});$ identity matrix. Similarly, the virtual strain is also approximated by: $var element = document.getElementById("moose-equation-3c5a59d9-f6bb-456a-9b63-7759cb1175bc");katex.render("\\begin{aligned} \\bar{\\delta \\epsilon} = \\delta \\epsilon + (\\delta w - \\frac{tr(\\delta \\epsilon)}{3})I \\end{aligned}", element, {displayMode:true,throwOnError:false});$ The modified equation of motion is: $var element = document.getElementById("moose-equation-b21b5c1a-1f29-4445-8408-d58cf70433ef");katex.render("\\begin{aligned} \\int_V \\sigma(\\bar{\\epsilon})\\bar{\\delta \\epsilon} dV - \\int_V b \\delta v dV - \\int_{\\partial V} t \\delta v dA = 0 \\end{aligned}", element, {displayMode:true,throwOnError:false});$ More details about this method can be found in Section 8.6 of Bower (2009).

When finite strain formulation is used, the volumetric component of the strain is separated using the determinant of the deformation matrix.

Usage

Volumetric locking correction is set to false by default in tensor mechanics. When dealing with problems involving plasticity or incompressible materials, it can be turned on by setting volumetric_locking_correction=true in both the stress divergence kernel and the strain calculator or in the Tensor Mechanics master action.

When volumetric locking correction is turned on, using a SMP preconditioner with coupled displacement variables may help with convergence. For a 3-D problem with only displacement as unknown variables, the following pre-conditioner block may be used:

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  [../]
[]

Verification of locking correction on Cook's membrane

A 2D trapezoidal membrane Figure 1 with Poisson's ratio of 0.4999 is fixed at the left edge and sheared at the right edge. Locking behavior of the membrane is observed with the use of first order (Quad4) elements with no volumetric locking correction. This locking results in much lower vertical displacement at Point A than that observed in other scenarios, see Figure 2 . Locking can be avoided with the use of Quad4 elements along with volumetric locking correction or with the use of second order elements (Quad8) with or without volumetric locking correction. These results match with that presented in Fig. 6 of Nakshatrala et al. (2008).

Figure 1: 2D problem to demonstrate volumetric locking.

Figure 2: Vertical displacement at Point A for different element types and mesh density. Locking behavior is observed when Quad4 elements with no volumetric locking correction are used.

Note that at least 20 nodes per edge are required to converge to the correct solution even when volumetric locking correction is used for this problem. Also, second order elements do not display locking behavior, so volumetric locking correction is not required for second order elements. Using volumetric locking correction with second order elements causes a different convergence pattern for the displacements Figure 2 and it is also shown to result in noisy zigzag pattern in stress or strain profiles.

References

A. F. Bower. Applied Mechanics of Solids. CRC press, 2009.

T. J. R Hughes. The Finite Element Method, Linear Static and Dynamic Finite Element Analysis. New Jersey:Prentice-Hall, 1987.

@book{hughes1987finite,
    author = "Hughes, T. J. R",
    title = "The Finite Element Method, Linear Static and Dynamic Finite Element Analysis",
    publisher = "New Jersey:Prentice-Hall",
    year = "1987"
}

K. B. Nakshatrala, A. Masud, and K. D. Hjelmstad. On finite element formulations for nearly incompressible linear elasticity. Computational Mechanics, 41:547–561, 2008.

@article{nakshatrala2008fem,
    author = "Nakshatrala, K. B. and Masud, A. and Hjelmstad, K. D.",
    journal = "Computational Mechanics",
    title = "On finite element formulations for nearly incompressible linear elasticity",
    volume = "41",
    pages = "547-561",
    year = "2008"
}

(moose/modules/tensor_mechanics/test/tests/finite_strain_elastic/elastic_rotation_test.i)

#
# Rotation Test
#
# This test is designed to compute stress based on uniaxial strain
# and then follow that stress as the mesh is rotated 90 degrees.
#
# The mesh is composed of one block with a single element.  The nodal
# displacements in the three directions are prescribed.  Poisson's
# ratio is 1/3, and Young's modulus is 1e6.
#
# This test is mentioned in
# K. Kamojjala, R. Brannon, A. Sadeghirad, and J. Guilkey, "Verification
#   tests in solid mechanics," Engineering with Computers, Vol. 31, 2015.
#   DOI: 10.1007/s00366-013-0342-x
#
[Mesh]
  type = FileMesh
  file = rotation_test.e
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[Functions]
  [./x_200]
    type = ParsedFunction
    symbol_names = 'delta t0'
    symbol_values = '1e-6 1.0'
    expression = 'if(t<=1.0, delta*t, (1.0+delta)*cos(pi/2*(t-t0)) - 1.0)'
  [../]
  [./y_200]
    type = ParsedFunction
    symbol_names = 'delta t0'
    symbol_values = '1e-6 1.0'
    expression = 'if(t<=1.0, 0.0, (1.0+delta)*sin(pi/2*(t-t0)))'
  [../]
  [./x_300]
    type = ParsedFunction
    symbol_names = 'delta t0'
    symbol_values = '1e-6 1.0'
    expression = 'if(t<=1.0, delta*t, (1.0+delta)*cos(pi/2.0*(t-t0)) - sin(pi/2.0*(t-t0)) - 1.0)'
  [../]
  [./y_300]
    type = ParsedFunction
    symbol_names = 'delta t0'
    symbol_values = '1e-6 1.0'
    expression = 'if(t<=1.0, 0.0, cos(pi/2.0*(t-t0)) + (1+delta)*sin(pi/2.0*(t-t0)) - 1.0)'
  [../]
  [./x_400]
    type = ParsedFunction
    symbol_names = 'delta t0'
    symbol_values = '1e-6 1.0'
    expression = 'if(t<=1.0, 0.0, -sin(pi/2.0*(t-t0)))'
  [../]
  [./y_400]
    type = ParsedFunction
    symbol_names = 'delta t0'
    symbol_values = '1e-6 1.0'
    expression = 'if(t<=1.0, 0.0, cos(pi/2.0*(t-t0)) - 1.0)'
  [../]
[]

[Modules]
  [./TensorMechanics]
    [./Master]
      [./all]
        strain = FINITE
        add_variables = true
        generate_output = 'stress_xx stress_yy stress_zz stress_xy stress_yz stress_zx'
      [../]
    [../]
  [../]
[]

[BCs]
  [./no_x]
    type = DirichletBC
    variable = disp_x
    boundary = 100
    value = 0.0
  [../]
  [./no_y]
    type = DirichletBC
    variable = disp_y
    boundary = 100
    value = 0.0
  [../]
  [./x_200]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = 200
    function = x_200
  [../]
  [./y_200]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = 200
    function = y_200
  [../]
  [./x_300]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = 300
    function = x_300
  [../]
  [./y_300]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = 300
    function = y_300
  [../]
  [./x_400]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = 400
    function = x_400
  [../]
  [./y_400]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = 400
    function = y_400
  [../]
  [./no_z]
    type = DirichletBC
    variable = disp_z
    boundary = '100 200 300 400'
    value = 0.0
  [../]
[]

[Materials]
  [./stress]
    type = ComputeFiniteStrainElasticStress
  [../]
  [./elasticity_tensor]
    type = ComputeElasticityTensor
    fill_method = symmetric9
    C_ijkl = '1.5e6 0.75e6 0.75e6 1.5e6 0.75e6 1.5e6 0.375e6 0.375e6 0.375e6'
  [../]
[]

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type '
  petsc_options_value = lu
  nl_rel_tol = 1e-30
  nl_abs_tol = 1e-20
  l_max_its = 20
  start_time = 0.0
  dt = 0.01
  end_time = 2.0
[]

[Outputs]
  exodus = true
[]