Linear Elasticity Problem

Summary

Solves a 3D linear elasticity problem for the deformation of a multi-material cantilever beam. This example is based on MFEM Example 2.

Description

This problem solves a linear elasticity problem for small deformations with strong form:

$var element = document.getElementById("moose-equation-262a59b9-e148-427b-81b1-862f2bb53911");katex.render("\\begin{split} -\\partial_j \\sigma_{ij} =0 \\,\\,\\,&\\mathrm{on}\\,\\, \\Omega \\\\ u_i = 0 \\,\\,\\, &\\mathrm{on}\\,\\, \\Gamma_1 \\\\ \\sigma_{ij} \\hat n_j= f_i \\,\\,\\, &\\mathrm{on}\\,\\, \\Gamma_2 \\end{split}", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

where Einstein notation for summation over repeated indices has been used, and the pull-down force $var element = document.getElementById("moose-equation-557ead6e-9854-464c-8a87-96c3e56f1275");katex.render("\\vec f = f_i \\hat e_i", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is given by

$var element = document.getElementById("moose-equation-0f95d148-245f-4574-9717-2ce25ce39783");katex.render("\\vec f = \\begin{pmatrix} 0.0\\\\ 0.0\\\\ -0.01 \\end{pmatrix}", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

In this example, the stress/strain relation is taken to be isotropic, such that

$var element = document.getElementById("moose-equation-8b7b3644-5fba-43da-929f-8adb31961a57");katex.render("\\begin{split} \\sigma_{ij} \\equiv C_{ijkl} \\varepsilon_{kl} = \\lambda \\delta_{ij} \\varepsilon_{kk} + 2 \\mu \\varepsilon_{ij} \\\\ \\varepsilon_{ij} = \\frac{1}{2} \\left(\\partial_i u_j + \\partial_j u_i\\right) \\end{split}", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

In this example, we solve this using the weak form

$var element = document.getElementById("moose-equation-d7be24c1-0e04-4aa8-805f-7c6523a24ee3");katex.render("(\\sigma_{ij}, \\partial_j v_i)_\\Omega - (\\sigma_{ij} \\hat n_j, v_i)_{\\partial \\Omega} = 0 \\,\\,\\,\\, \\forall v_i \\in V", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

where

$var element = document.getElementById("moose-equation-dc55369d-169a-4bb9-9d17-67d5b093da0a");katex.render("v_i \\in H^1(\\Omega) : v_i = 0 \\,\\,\\, \\mathrm{on} \\,\\, \\Gamma_1", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

Material Parameters

In this example, the cantilever beam is comprised of two materials; a stiffer material on the block with domain attribute 1, and a more flexible material defined on the block with domain attribute 2. These materials are parametrised with different values for the Lamé parameters, $var element = document.getElementById("moose-equation-ab63e596-d87f-43e6-bc01-4827314a40de");katex.render("\\lambda", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ and $var element = document.getElementById("moose-equation-9ab3e7ec-e87a-489c-ab64-9c19c696e10f");katex.render("\\mu", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , on the two domains.

The two Lamé parameters can be expressed in terms of the Young's modulus $var element = document.getElementById("moose-equation-63da8b6c-8338-4f5b-a87f-026f84390ca3");katex.render("E", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ and the Poisson ratio $var element = document.getElementById("moose-equation-bb1e0bc4-b509-4fda-bed8-3191f3ba9f34");katex.render("\\nu", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ in each material, using

$var element = document.getElementById("moose-equation-d3ef11a5-92eb-4db5-9ac5-fc7774f0c5c2");katex.render("\\begin{split} \\lambda &= \\frac{E\\nu}{(1-2\\nu)(1+\\nu)} \\\\ \\mu &= \\frac{E}{2(1+\\nu)} \\end{split}", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

Example File

[Mesh<<<{"href": "../Mesh/index.html"}>>>]
  type = MFEMMesh
  file = ../mesh/beam-tet.mesh
  dim = 3
  uniform_refine<<<{"description": "Specify the level of uniform refinement applied to the initial mesh"}>>> = 2
  displacement = "displacement"
[]

[Problem<<<{"href": "../Problem/index.html"}>>>]
  type = MFEMProblem
[]

[FESpaces<<<{"href": "../FESpaces/index.html"}>>>]
  [H1FESpace]
    type = MFEMVectorFESpace<<<{"description": "Convenience class to construct vector finite element spaces, abstracting away some of the mathematical complexity of specifying the dimensions.", "href": "../../source/mfem/fespaces/MFEMVectorFESpace.html"}>>>
    fec_type<<<{"description": "Specifies the family of FE shape functions."}>>> = H1
    fec_order<<<{"description": "Order of the FE shape function to use."}>>> = FIRST
    range_dim<<<{"description": "The number of components of the vectors in reference space. Zero (the default) means it will be the same as the problem dimension. Note that MFEM does not currently support 2D vectors in 1D space for ND and RT elements."}>>> = 3
    ordering<<<{"description": "Ordering style to use for vector DoFs."}>>> = "vdim"
  []
[]

[Variables<<<{"href": "../Variables/index.html"}>>>]
  [displacement]
    type = MFEMVariable<<<{"description": "Class for adding MFEM variables to the problem (`mfem::ParGridFunction`s).", "href": "../../source/mfem/variables/MFEMVariable.html"}>>>
    fespace<<<{"description": "The finite element space this variable is defined on."}>>> = H1FESpace
  []
[]

[BCs<<<{"href": "../BCs/index.html"}>>>]
  [dirichlet]
    type = MFEMVectorDirichletBC<<<{"description": "Applies a Dirichlet condition to all components of a vector variable.", "href": "../../source/mfem/bcs/MFEMVectorDirichletBC.html"}>>>
    variable<<<{"description": "Variable on which to apply the boundary condition"}>>> = displacement
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies. Defaults to all boundaries."}>>> = '1'
  []
  [pull_down]
    type = MFEMVectorBoundaryIntegratedBC<<<{"description": "Adds the boundary integrator to an MFEM problem for the linear form $(\\vec f, \\vec v)_{\\partial\\Omega}$", "href": "../../source/mfem/bcs/MFEMVectorBoundaryIntegratedBC.html"}>>>
    variable<<<{"description": "Variable on which to apply the boundary condition"}>>> = displacement
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies. Defaults to all boundaries."}>>> = '2'
    vector_coefficient<<<{"description": "Vector coefficient used in the boundary integrator. A functor is any of the following: a variable, an MFEM material property, a function, a postprocessor or a numeric vector value (enclosed in curly braces)."}>>> = '0.0 0.0 -0.01'
  []
[]

[FunctorMaterials<<<{"href": "../FunctorMaterials/index.html"}>>>]
  [Rigidium]
    type = MFEMGenericFunctorMaterial<<<{"description": "Declares material scalar properties based on names and coefficients prescribed by input parameters.", "href": "../../source/mfem/functormaterials/MFEMGenericFunctorMaterial.html"}>>>
    prop_names<<<{"description": "The names of the properties this material will have"}>>> = 'lambda mu'
    prop_values<<<{"description": "The corresponding names of coefficients associated with the named properties. A functor is any of the following: a variable, an MFEM material property, a function, a postprocessor or a number."}>>> = '50.0 50.0'
    block<<<{"description": "The list of subdomains (names or ids) that this object will be restricted to. Leave empty to apply to all subdomains."}>>> = 1
  []
  [Bendium]
    type = MFEMGenericFunctorMaterial<<<{"description": "Declares material scalar properties based on names and coefficients prescribed by input parameters.", "href": "../../source/mfem/functormaterials/MFEMGenericFunctorMaterial.html"}>>>
    prop_names<<<{"description": "The names of the properties this material will have"}>>> = 'lambda mu'
    prop_values<<<{"description": "The corresponding names of coefficients associated with the named properties. A functor is any of the following: a variable, an MFEM material property, a function, a postprocessor or a number."}>>> = '1.0 1.0'
    block<<<{"description": "The list of subdomains (names or ids) that this object will be restricted to. Leave empty to apply to all subdomains."}>>> = 2
  []
[]

[Kernels<<<{"href": "../Kernels/index.html"}>>>]
  [diff]
    type = MFEMLinearElasticityKernel<<<{"description": "The isotropic linear elasticity operator with weak form $(c_{ikjl} \\nabla u_j, \\nabla v_i)$, to be added to an MFEM problem, where $c_{ikjl}$ is the isotropic elasticity tensor, $c_{ikjl} = \\lambda \\delta_{ik} \\delta_{jl} + \\mu \\left( \\delta_{ij} \\delta_{kl} + \\delta_{il} \\delta_{jk} \\right)$, $\\lambda$ is the first Lame parameter, $\\lambda = \\frac{E\\nu}{(1-2\\nu)(1+\\nu)}$, $\\mu$ is the second Lame parameter, $\\mu = \\frac{E}{2(1+\\nu)}$, where $E$ is Young's modulus and $\\nu$ is Poisson's ratio.", "href": "../../source/mfem/kernels/MFEMLinearElasticityKernel.html"}>>>
    variable<<<{"description": "Variable labelling the weak form this kernel is added to"}>>> = displacement
    lambda<<<{"description": "Name of MFEM Lame constant lambda to multiply the div(u)*I term by. A functor is any of the following: a variable, an MFEM material property, a function, a postprocessor or a number."}>>> = lambda
    mu<<<{"description": "Name of MFEM Lame constant mu to multiply the gradients term by. A functor is any of the following: a variable, an MFEM material property, a function, a postprocessor or a number."}>>> = mu
  []
[]

[Preconditioner<<<{"href": "../Preconditioner/index.html"}>>>]
  [boomeramg]
    type = MFEMHypreBoomerAMG<<<{"description": "Hypre BoomerAMG solver and preconditioner for the iterative solution of MFEM equation systems.", "href": "../../source/mfem/solvers/MFEMHypreBoomerAMG.html"}>>>
    l_max_its<<<{"description": "Set the maximum number of iterations."}>>> = 500
    l_tol<<<{"description": "Set the relative tolerance."}>>> = 1e-8
    print_level<<<{"description": "Set the solver verbosity."}>>> = 2
  []
[]

[Solver<<<{"href": "../Solver/index.html"}>>>]
  type = MFEMHyprePCG
  #preconditioner = boomeramg
  l_max_its = 5000
  l_tol = 1e-8
  l_abs_tol = 0.0
  print_level = 2
[]

[Executioner<<<{"href": "../Executioner/index.html"}>>>]
  type = MFEMSteady
  device = "cpu"
[]

[Outputs<<<{"href": "../Outputs/index.html"}>>>]
  [ParaViewDataCollection]
    type = MFEMParaViewDataCollection<<<{"description": "Output for controlling export to an mfem::ParaViewDataCollection.", "href": "../../source/mfem/outputs/MFEMParaViewDataCollection.html"}>>>
    file_base<<<{"description": "The desired solution output name without an extension. If not provided, MOOSE sets it with Outputs/file_base when available. Otherwise, MOOSE uses input file name and this object name for a master input or uses master file_base, the subapp name and this object name for a subapp input to set it."}>>> = OutputData/LinearElasticity
    vtk_format<<<{"description": "Select VTK data format to use, choosing between BINARY, BINARY32, and ASCII."}>>> = ASCII
  []
[]

(test/tests/mfem/kernels/linearelasticity.i)

[Mesh]
  type = MFEMMesh
  file = ../mesh/beam-tet.mesh
  dim = 3
  uniform_refine = 2
  displacement = "displacement"
[]

[Problem]
  type = MFEMProblem
[]

[FESpaces]
  [H1FESpace]
    type = MFEMVectorFESpace
    fec_type = H1
    fec_order = FIRST
    range_dim = 3
    ordering = "vdim"
  []
[]

[Variables]
  [displacement]
    type = MFEMVariable
    fespace = H1FESpace
  []
[]

[BCs]
  [dirichlet]
    type = MFEMVectorDirichletBC
    variable = displacement
    boundary = '1'
  []
  [pull_down]
    type = MFEMVectorBoundaryIntegratedBC
    variable = displacement
    boundary = '2'
    vector_coefficient = '0.0 0.0 -0.01'
  []
[]

[FunctorMaterials]
  [Rigidium]
    type = MFEMGenericFunctorMaterial
    prop_names = 'lambda mu'
    prop_values = '50.0 50.0'
    block = 1
  []
  [Bendium]
    type = MFEMGenericFunctorMaterial
    prop_names = 'lambda mu'
    prop_values = '1.0 1.0'
    block = 2
  []
[]

[Kernels]
  [diff]
    type = MFEMLinearElasticityKernel
    variable = displacement
    lambda = lambda
    mu = mu
  []
[]

[Preconditioner]
  [boomeramg]
    type = MFEMHypreBoomerAMG
    l_max_its = 500
    l_tol = 1e-8
    print_level = 2
  []
[]

[Solver]
  type = MFEMHyprePCG
  #preconditioner = boomeramg
  l_max_its = 5000
  l_tol = 1e-8
  l_abs_tol = 0.0
  print_level = 2
[]

[Executioner]
  type = MFEMSteady
  device = "cpu"
[]

[Outputs]
  [ParaViewDataCollection]
    type = MFEMParaViewDataCollection
    file_base = OutputData/LinearElasticity
    vtk_format = ASCII
  []
[]

Summary
Description
Material Parameters
Example File