Contact Verification

Overview

Finite element contact addresses the interaction of two surfaces that apply loads to one another.

Hertz Contact of Cylinder and Rigid Surface

The problem is an infinite cylinder in contact with a rigid plane. The problem is taken from Wang and Zhu (2013). See also Johnson (1985).

The analytic expressions for this problem rely on $var element = document.getElementById("moose-equation-35dd49c8-72de-4486-bba8-335628ee94a5");katex.render("R_e", element, {displayMode:false,throwOnError:false});$ , the radius of a so-called equivalent cylinder.

$var element = document.getElementById("moose-equation-7e3715b3-658d-483f-ba4c-5ec3414f2a20");katex.render(" R_e = \\frac{R_1+R_2}{R_1 R_2}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-9d5712d0-f266-4838-bdb0-5fb7dd47f906");katex.render("R_1", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-7ec49588-3e91-48e0-8624-036c9b1212a6");katex.render("R_2", element, {displayMode:false,throwOnError:false});$ are the radii of two cylinders in contact that have parallel axes. In our case, $var element = document.getElementById("moose-equation-4c720333-be37-40e5-aba8-122469cfa8e3");katex.render("R_1", element, {displayMode:false,throwOnError:false});$ = 3, and $var element = document.getElementById("moose-equation-b55e8300-5ed9-4ed1-b279-bdeca3e04138");katex.render("R_2", element, {displayMode:false,throwOnError:false});$ is infinity. This gives $var element = document.getElementById("moose-equation-fc392038-cd95-4750-ab25-011866b0934f");katex.render("R_e", element, {displayMode:false,throwOnError:false});$ = 3.

We also require $var element = document.getElementById("moose-equation-224dc0ef-a3f2-419a-a1e3-27dad9ea4b97");katex.render("E^*", element, {displayMode:false,throwOnError:false});$ : $var element = document.getElementById("moose-equation-4e682829-77f6-4d5c-b47d-f070a206967e");katex.render(" \\frac{1}{E^*} = \\frac{1-\\nu_1^2}{E_1} + \\frac{1-\\nu_2^2}{E_2}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-0b449d83-a95a-4ff1-97a0-20fe6b6ab406");katex.render("E", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-43b5e8bc-cbee-4932-9abf-81d039cd6b00");katex.render("\\nu", element, {displayMode:false,throwOnError:false});$ are Young's modulus and Poisson's ratio with subscripts representing each of the two bodies. In our case, $var element = document.getElementById("moose-equation-eb4e1af0-27e2-4bae-aab1-0f45c590c6f8");katex.render("E_1", element, {displayMode:false,throwOnError:false});$ = 1E6, $var element = document.getElementById("moose-equation-865c6f5f-38d2-448a-bbe0-88891594455b");katex.render("\\nu_1", element, {displayMode:false,throwOnError:false});$ = 0.3, $var element = document.getElementById("moose-equation-f135e064-aab1-447e-be82-cdd3a055782b");katex.render("E_2", element, {displayMode:false,throwOnError:false});$ = infinity, and $var element = document.getElementById("moose-equation-be55f376-3496-4399-856e-2ad20c6a2b35");katex.render("\\nu_2", element, {displayMode:false,throwOnError:false});$ = 0. This gives $var element = document.getElementById("moose-equation-e72be4f9-ae54-4ede-877b-bb027520da6e");katex.render("1/E^*", element, {displayMode:false,throwOnError:false});$ = 9.1E-7.

$var element = document.getElementById("moose-equation-597e7aed-3bf2-4ae0-ae22-8c0afe7ca3b7");katex.render("P", element, {displayMode:false,throwOnError:false});$ is the line load per unit length. We take this to be a pressure applied to the bottom half of the cylinder, 2E3, times $var element = document.getElementById("moose-equation-ef2b0eb6-6250-44c2-ac4f-25a2171d1398");katex.render("2R_1", element, {displayMode:false,throwOnError:false});$ . $var element = document.getElementById("moose-equation-16d74225-3e88-4e19-8a8c-ac89465f6a6c");katex.render("P", element, {displayMode:false,throwOnError:false});$ = 1.2E4.

The peak load, which occurs at the centerline of the cylinder, is $var element = document.getElementById("moose-equation-00f3643b-131c-4a57-ad71-6fdb5e5b0ed7");katex.render(" P_{\\text{max}} = \\frac{2P}{\\pi a}", element, {displayMode:true,throwOnError:false});$ with $var element = document.getElementById("moose-equation-7dcae8e8-8906-4940-94ed-172c5a806574");katex.render(" a = \\sqrt{\\frac{4 P R_e}{\\pi E^*}}", element, {displayMode:true,throwOnError:false});$

This gives $var element = document.getElementById("moose-equation-28371e55-1e64-4c97-bfc6-fff15b399c5e");katex.render("a", element, {displayMode:false,throwOnError:false});$ = 0.204 and $var element = document.getElementById("moose-equation-79b8fea2-db0c-45da-a685-7e0264226219");katex.render("P_{\\text{max}}", element, {displayMode:false,throwOnError:false});$ = 37405. Bison computes $var element = document.getElementById("moose-equation-031fb103-8e4b-4ad6-af3f-2c30707471be");katex.render("P_{\\text{max}}", element, {displayMode:false,throwOnError:false});$ as ~3.80E4, giving an error of 1.52%.

Figure 1: Contact traction for Hertz cylinder contact problem.

Input Syntax

#
# Hertz Contact of Cylinder and Rigid Surface
#   See https://doi.org/10.1007/978-0-387-92897-5_491
#
# 1/E* = (1-v1^2)/E1 + (1-v2^2)/E2
# v1 = 0.3
# E1 = 1e6
# v2 = 0
# E2 = infinity
# 1/E* = (1-0.09)/1e6 = 9.1e-7
# EstarInv = 9.1e-7
# 1/Re = (R1+R2)/R1/R2
# R1 = 3
# R2 = infinity
# 1/Re = 1/3
# Re = 3.0
# P = pressure at 0.1 seconds * 2*R1
# P = 2e3 * 6
# P = 1.2e4
#
# a is half width of contacting area
# a = sqrt(4 * P * Re / pi * EstarInv)
# a = 0.204
#
# Pmax = 2 * P / pi / a
# Pmax = 37405
#

[GlobalParams<<<{"href": "../../syntax/GlobalParams/index.html"}>>>]
  displacements = 'disp_x disp_y'
[]

[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
  [input_file]
    type = FileMeshGenerator<<<{"description": "Read a mesh from a file.", "href": "../../source/meshgenerators/FileMeshGenerator.html"}>>>
    file<<<{"description": "The filename to read."}>>> = hertz_cyl_coarser.e
  []
  [secondary]
    type = LowerDBlockFromSidesetGenerator<<<{"description": "Adds lower dimensional elements on the specified sidesets.", "href": "../../source/meshgenerators/LowerDBlockFromSidesetGenerator.html"}>>>
    new_block_id<<<{"description": "The lower dimensional block id to create"}>>> = 10001
    new_block_name<<<{"description": "The lower dimensional block name to create (optional)"}>>> = 'secondary_lower'
    sidesets<<<{"description": "The sidesets from which to create the new block"}>>> = '3'
    input<<<{"description": "The mesh we want to modify"}>>> = input_file
  []
  [primary]
    type = LowerDBlockFromSidesetGenerator<<<{"description": "Adds lower dimensional elements on the specified sidesets.", "href": "../../source/meshgenerators/LowerDBlockFromSidesetGenerator.html"}>>>
    new_block_id<<<{"description": "The lower dimensional block id to create"}>>> = 10000
    sidesets<<<{"description": "The sidesets from which to create the new block"}>>> = '2'
    new_block_name<<<{"description": "The lower dimensional block name to create (optional)"}>>> = 'primary_lower'
    input<<<{"description": "The mesh we want to modify"}>>> = secondary
  []
[]

[Problem<<<{"href": "../../syntax/Problem/index.html"}>>>]
  type = ReferenceResidualProblem
  extra_tag_vectors = 'ref'
  reference_vector = 'ref'
[]

[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
  [disp_x]
  []
  [disp_y]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = -0.01
  []
  [frictionless_normal_lm]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = FIRST
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = LAGRANGE
    block = 'secondary_lower'
    use_dual = true
  []
[]

[AuxVariables<<<{"href": "../../syntax/AuxVariables/index.html"}>>>]
  [saved_x]
  []
  [saved_y]
  []
[]

[Functions<<<{"href": "../../syntax/Functions/index.html"}>>>]
  [pressure_function]
    type = PiecewiseLinear<<<{"description": "Linearly interpolates between pairs of x-y data", "href": "../../source/functions/PiecewiseLinear.html"}>>>
    x<<<{"description": "The abscissa values"}>>> = '0. 1. 3.5'
    y<<<{"description": "The ordinate values"}>>> = '0. 2e4 2e4'
  []
  [analytic]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../source/functions/MooseParsedFunction.html"}>>>
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'EstarInv Re width P'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = '9.1e-7 3.0 6 pressure_function'
    expression<<<{"description": "The user defined function."}>>> = '2 * P*width / pi / sqrt(4 * P*width * Re / pi * EstarInv)'
  []
[]

[Physics<<<{"href": "../../syntax/Physics/index.html"}>>>/SolidMechanics<<<{"href": "../../syntax/Physics/SolidMechanics/index.html"}>>>/QuasiStatic<<<{"href": "../../syntax/Physics/SolidMechanics/QuasiStatic/index.html"}>>>]
  [fred]
    block<<<{"description": "The list of ids of the blocks (subdomain) that the stress divergence kernels will be applied to"}>>> = '1 2 3 4 5 6 7'
    strain<<<{"description": "Strain formulation"}>>> = finite
    save_in<<<{"description": "The displacement residuals"}>>> = 'saved_x saved_y'
    extra_vector_tags<<<{"description": "The tag names for extra vectors that residual data should be saved into"}>>> = 'ref'
    generate_output<<<{"description": "Add scalar quantity output for stress and/or strain"}>>> = 'stress_xx stress_yy stress_xy'
  []
[]

[Postprocessors<<<{"href": "../../syntax/Postprocessors/index.html"}>>>]
  [traction_applied]
    type = NodalSum<<<{"description": "Computes the sum of all of the nodal values of the specified variable. Note: This object sets the default \"unique_node_execute\" flag to true to avoid double counting nodes between shared blocks.", "href": "../../source/postprocessors/NodalSum.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = saved_y
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 4
  []
  [_dt]
    type = TimestepSize<<<{"description": "Reports the timestep size", "href": "../../source/postprocessors/TimestepSize.html"}>>>
  []
  [bison_peak_traction]
    type = NodalExtremeValue<<<{"description": "Finds either the min or max elemental value of a variable over the domain.", "href": "../../source/postprocessors/NodalExtremeValue.html"}>>>
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 3
    variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = frictionless_normal_lm
  []
  [analytical_peak_traction]
    type = FunctionValuePostprocessor<<<{"description": "Computes the value of a supplied function at a single point (scalable)", "href": "../../source/postprocessors/FunctionValuePostprocessor.html"}>>>
    function<<<{"description": "The function which supplies the postprocessor value."}>>> = analytic
  []
  [error]
    type = ParsedPostprocessor<<<{"description": "Computes a parsed expression with post-processors", "href": "../../source/postprocessors/ParsedPostprocessor.html"}>>>
    pp_names<<<{"description": "Post-processors arguments"}>>> = 'bison_peak_traction analytical_peak_traction'
    expression<<<{"description": "function expression"}>>> = '(bison_peak_traction-analytical_peak_traction)/analytical_peak_traction*100'
  []
[]

[BCs<<<{"href": "../../syntax/NuclearMaterials/BCs/index.html"}>>>]
  [base_x]
    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_x
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = '1 2'
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [base_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"}>>> = '1 2'
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [Pressure<<<{"href": "../../syntax/BCs/Pressure/index.html"}>>>]
    [pressure]
      boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 4
      function<<<{"description": "The function that describes the pressure"}>>> = pressure_function
      use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used. For small strain formulations pressure should be applied to the undisplaced mesh to obtain agreement with analytical benchmark solutions."}>>> = false
    []
  []
[]

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [elas_tens]
    type = ComputeIsotropicElasticityTensor<<<{"description": "Compute a constant isotropic elasticity tensor.", "href": "../../source/materials/ComputeIsotropicElasticityTensor.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = '1 2 3 4 5 6 7'
    youngs_modulus<<<{"description": "Young's modulus of the material."}>>> = 1e6
    poissons_ratio<<<{"description": "Poisson's ratio for the material."}>>> = 0.3
  []
  [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"}>>> = '1 2 3 4 5 6 7'
  []
[]

[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -pc_factor_shift_type -pc_factor_shift_amount -mat_mffd_err'
  petsc_options_value = 'lu       superlu_dist               NONZERO               1e-15                   1e-5'

  line_search = 'none'

  nl_abs_tol = 1e-7

  start_time = 0.0
  end_time = 0.1
  l_tol = 1e-4
  dt = 0.1
  dtmin = 0.1

  [Predictor<<<{"href": "../../syntax/Executioner/Predictor/index.html"}>>>]
    type = SimplePredictor
    scale = 1.0
  []
[]

[Preconditioning<<<{"href": "../../syntax/Preconditioning/index.html"}>>>]
  [SMP]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
  []
[]

[VectorPostprocessors<<<{"href": "../../syntax/VectorPostprocessors/index.html"}>>>]
  [cont_press]
    type = NodalValueSampler<<<{"description": "Samples values of nodal variable(s).", "href": "../../source/vectorpostprocessors/NodalValueSampler.html"}>>>
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = frictionless_normal_lm
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = '3'
    sort_by<<<{"description": "What to sort the samples by"}>>> = id
  []
[]

[Outputs<<<{"href": "../../syntax/Outputs/index.html"}>>>]
  print_linear_residuals<<<{"description": "Enable printing of linear residuals to the screen (Console)"}>>> = true
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
  [console]
    type = Console<<<{"description": "Object for screen output.", "href": "../../source/outputs/Console.html"}>>>
    max_rows<<<{"description": "The maximum number of postprocessor/scalar values displayed on screen during a timestep (set to 0 for unlimited)"}>>> = 5
  []
  [csv1]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../source/outputs/CSV.html"}>>>
    show<<<{"description": "A list of the variables and postprocessors that should be output to the Exodus file (may include Variables, ScalarVariables, and Postprocessor names)."}>>> = 'cont_press'
    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."}>>> = cylinder_frictionless
  []
[]

[UserObjects<<<{"href": "../../syntax/UserObjects/index.html"}>>>]
  [weighted_vel_uo]
    type = LMWeightedGapUserObject<<<{"description": "Provides the mortar normal Lagrange multiplier for constraint enforcement.", "href": "../../source/userobjects/LMWeightedGapUserObject.html"}>>>
    primary_boundary<<<{"description": "The name of the primary boundary sideset."}>>> = 2
    secondary_boundary<<<{"description": "The name of the secondary boundary sideset."}>>> = 3
    primary_subdomain<<<{"description": "The name of the primary subdomain."}>>> = 10000
    secondary_subdomain<<<{"description": "The name of the secondary subdomain."}>>> = 10001
    disp_x<<<{"description": "The x displacement variable"}>>> = disp_x
    disp_y<<<{"description": "The y displacement variable"}>>> = disp_y
    lm_variable<<<{"description": "The Lagrange multiplier variable representing the contact pressure."}>>> = frictionless_normal_lm
  []
[]

[Constraints<<<{"href": "../../syntax/Constraints/index.html"}>>>]
  [weighted_gap_lm]
    type = ComputeWeightedGapLMMechanicalContact<<<{"description": "Computes the weighted gap that will later be used to enforce the zero-penetration mechanical contact conditions", "href": "../../source/constraints/ComputeWeightedGapLMMechanicalContact.html"}>>>
    primary_boundary<<<{"description": "The name of the primary boundary sideset."}>>> = 2
    secondary_boundary<<<{"description": "The name of the secondary boundary sideset."}>>> = 3
    primary_subdomain<<<{"description": "The name of the primary subdomain."}>>> = 10000
    secondary_subdomain<<<{"description": "The name of the secondary subdomain."}>>> = 10001
    variable<<<{"description": "The name of the lagrange multiplier variable that this constraint is applied to. This parameter may not be supplied in the case of using penalty methods for example"}>>> = frictionless_normal_lm
    disp_x<<<{"description": "The x displacement variable"}>>> = disp_x
    disp_y<<<{"description": "The y displacement variable"}>>> = disp_y
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation.  Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = true
    c<<<{"description": "Parameter for balancing the size of the gap and contact pressure"}>>> = 1.0e6
    weighted_gap_uo<<<{"description": "The weighted gap user object"}>>> = weighted_vel_uo
  []
  [x]
    type = NormalMortarMechanicalContact<<<{"description": "This class is used to apply normal contact forces using lagrange multipliers", "href": "../../source/constraints/NormalMortarMechanicalContact.html"}>>>
    primary_boundary<<<{"description": "The name of the primary boundary sideset."}>>> = '2'
    secondary_boundary<<<{"description": "The name of the secondary boundary sideset."}>>> = '3'
    primary_subdomain<<<{"description": "The name of the primary subdomain."}>>> = '10000'
    secondary_subdomain<<<{"description": "The name of the secondary subdomain."}>>> = '10001'
    variable<<<{"description": "The name of the lagrange multiplier variable that this constraint is applied to. This parameter may not be supplied in the case of using penalty methods for example"}>>> = frictionless_normal_lm
    secondary_variable<<<{"description": "Primal variable on secondary surface."}>>> = disp_x
    component<<<{"description": "The force component constraint that this object is supplying"}>>> = x
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation.  Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = true
    compute_lm_residuals<<<{"description": "Whether to compute Lagrange Multiplier residuals"}>>> = false
    weighted_gap_uo<<<{"description": "The weighted gap user object."}>>> = weighted_vel_uo
  []
  [y]
    type = NormalMortarMechanicalContact<<<{"description": "This class is used to apply normal contact forces using lagrange multipliers", "href": "../../source/constraints/NormalMortarMechanicalContact.html"}>>>
    primary_boundary<<<{"description": "The name of the primary boundary sideset."}>>> = '2'
    secondary_boundary<<<{"description": "The name of the secondary boundary sideset."}>>> = '3'
    primary_subdomain<<<{"description": "The name of the primary subdomain."}>>> = '10000'
    secondary_subdomain<<<{"description": "The name of the secondary subdomain."}>>> = '10001'
    variable<<<{"description": "The name of the lagrange multiplier variable that this constraint is applied to. This parameter may not be supplied in the case of using penalty methods for example"}>>> = frictionless_normal_lm
    secondary_variable<<<{"description": "Primal variable on secondary surface."}>>> = disp_y
    component<<<{"description": "The force component constraint that this object is supplying"}>>> = y
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation.  Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = true
    compute_lm_residuals<<<{"description": "Whether to compute Lagrange Multiplier residuals"}>>> = false
    weighted_gap_uo<<<{"description": "The weighted gap user object."}>>> = weighted_vel_uo
  []
[]

References

K. L. Johnson. Contact Mechanics. Cambridge University Press, 1985. doi:10.1017/CBO9781139171731.

@book{Johnson1985,
    author = "Johnson, K. L.",
    place = "Cambridge",
    title = "Contact Mechanics",
    DOI = "10.1017/CBO9781139171731",
    publisher = "Cambridge University Press",
    year = "1985"
}

Q. Jane Wang and Dong Zhu. Hertz Theory: Contact of Cylindrical Surfaces, pages 1639–1647. Springer US, Boston, MA, 2013. doi:10.1007/978-0-387-92897-5_491.

@Inbook{Wang2013,
    author = "Wang, Q. Jane and Zhu, Dong",
    editor = "Wang, Q. Jane and Chung, Yip-Wah",
    title = "Hertz Theory: Contact of Cylindrical Surfaces",
    bookTitle = "Encyclopedia of Tribology",
    year = "2013",
    publisher = "Springer US",
    address = "Boston, MA",
    pages = "1639--1647",
    isbn = "978-0-387-92897-5",
    doi = "10.1007/978-0-387-92897-5\_491"
}

(assessment/verification/Contact/hertz/cylinder_frictionless_pressure.i)

#
# Hertz Contact of Cylinder and Rigid Surface
#   See https://doi.org/10.1007/978-0-387-92897-5_491
#
# 1/E* = (1-v1^2)/E1 + (1-v2^2)/E2
# v1 = 0.3
# E1 = 1e6
# v2 = 0
# E2 = infinity
# 1/E* = (1-0.09)/1e6 = 9.1e-7
# EstarInv = 9.1e-7
# 1/Re = (R1+R2)/R1/R2
# R1 = 3
# R2 = infinity
# 1/Re = 1/3
# Re = 3.0
# P = pressure at 0.1 seconds * 2*R1
# P = 2e3 * 6
# P = 1.2e4
#
# a is half width of contacting area
# a = sqrt(4 * P * Re / pi * EstarInv)
# a = 0.204
#
# Pmax = 2 * P / pi / a
# Pmax = 37405
#

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [input_file]
    type = FileMeshGenerator
    file = hertz_cyl_coarser.e
  []
  [secondary]
    type = LowerDBlockFromSidesetGenerator
    new_block_id = 10001
    new_block_name = 'secondary_lower'
    sidesets = '3'
    input = input_file
  []
  [primary]
    type = LowerDBlockFromSidesetGenerator
    new_block_id = 10000
    sidesets = '2'
    new_block_name = 'primary_lower'
    input = secondary
  []
[]

[Problem]
  type = ReferenceResidualProblem
  extra_tag_vectors = 'ref'
  reference_vector = 'ref'
[]

[Variables]
  [disp_x]
  []
  [disp_y]
    initial_condition = -0.01
  []
  [frictionless_normal_lm]
    order = FIRST
    family = LAGRANGE
    block = 'secondary_lower'
    use_dual = true
  []
[]

[AuxVariables]
  [saved_x]
  []
  [saved_y]
  []
[]

[Functions]
  [pressure_function]
    type = PiecewiseLinear
    x = '0. 1. 3.5'
    y = '0. 2e4 2e4'
  []
  [analytic]
    type = ParsedFunction
    symbol_names = 'EstarInv Re width P'
    symbol_values = '9.1e-7 3.0 6 pressure_function'
    expression = '2 * P*width / pi / sqrt(4 * P*width * Re / pi * EstarInv)'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fred]
    block = '1 2 3 4 5 6 7'
    strain = finite
    save_in = 'saved_x saved_y'
    extra_vector_tags = 'ref'
    generate_output = 'stress_xx stress_yy stress_xy'
  []
[]

[Postprocessors]
  [traction_applied]
    type = NodalSum
    variable = saved_y
    boundary = 4
  []
  [_dt]
    type = TimestepSize
  []
  [bison_peak_traction]
    type = NodalExtremeValue
    boundary = 3
    variable = frictionless_normal_lm
  []
  [analytical_peak_traction]
    type = FunctionValuePostprocessor
    function = analytic
  []
  [error]
    type = ParsedPostprocessor
    pp_names = 'bison_peak_traction analytical_peak_traction'
    expression = '(bison_peak_traction-analytical_peak_traction)/analytical_peak_traction*100'
  []
[]

[BCs]
  [base_x]
    type = DirichletBC
    variable = disp_x
    boundary = '1 2'
    value = 0.0
  []
  [base_y]
    type = DirichletBC
    variable = disp_y
    boundary = '1 2'
    value = 0.0
  []
  [Pressure]
    [pressure]
      boundary = 4
      function = pressure_function
      use_displaced_mesh = false
    []
  []
[]

[Materials]
  [elas_tens]
    type = ComputeIsotropicElasticityTensor
    block = '1 2 3 4 5 6 7'
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
    block = '1 2 3 4 5 6 7'
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -pc_factor_shift_type -pc_factor_shift_amount -mat_mffd_err'
  petsc_options_value = 'lu       superlu_dist               NONZERO               1e-15                   1e-5'

  line_search = 'none'

  nl_abs_tol = 1e-7

  start_time = 0.0
  end_time = 0.1
  l_tol = 1e-4
  dt = 0.1
  dtmin = 0.1

  [Predictor]
    type = SimplePredictor
    scale = 1.0
  []
[]

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

[VectorPostprocessors]
  [cont_press]
    type = NodalValueSampler
    variable = frictionless_normal_lm
    boundary = '3'
    sort_by = id
  []
[]

[Outputs]
  print_linear_residuals = true
  exodus = true
  csv = true
  [console]
    type = Console
    max_rows = 5
  []
  [csv1]
    type = CSV
    show = 'cont_press'
    file_base = cylinder_frictionless
  []
[]

[UserObjects]
  [weighted_vel_uo]
    type = LMWeightedGapUserObject
    primary_boundary = 2
    secondary_boundary = 3
    primary_subdomain = 10000
    secondary_subdomain = 10001
    disp_x = disp_x
    disp_y = disp_y
    lm_variable = frictionless_normal_lm
  []
[]

[Constraints]
  [weighted_gap_lm]
    type = ComputeWeightedGapLMMechanicalContact
    primary_boundary = 2
    secondary_boundary = 3
    primary_subdomain = 10000
    secondary_subdomain = 10001
    variable = frictionless_normal_lm
    disp_x = disp_x
    disp_y = disp_y
    use_displaced_mesh = true
    c = 1.0e6
    weighted_gap_uo = weighted_vel_uo
  []
  [x]
    type = NormalMortarMechanicalContact
    primary_boundary = '2'
    secondary_boundary = '3'
    primary_subdomain = '10000'
    secondary_subdomain = '10001'
    variable = frictionless_normal_lm
    secondary_variable = disp_x
    component = x
    use_displaced_mesh = true
    compute_lm_residuals = false
    weighted_gap_uo = weighted_vel_uo
  []
  [y]
    type = NormalMortarMechanicalContact
    primary_boundary = '2'
    secondary_boundary = '3'
    primary_subdomain = '10000'
    secondary_subdomain = '10001'
    variable = frictionless_normal_lm
    secondary_variable = disp_y
    component = y
    use_displaced_mesh = true
    compute_lm_residuals = false
    weighted_gap_uo = weighted_vel_uo
  []
[]

Overview
Input Syntax
References