Step 1 - First Contact

Continuing from the final step in the solid mechanics tutorial we add a first attempt at mechanical contact.

#
# A first attempt at mechanical contact
# https://mooseframework.inl.gov/modules/contact/tutorials/introduction/step01.html
#

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

[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  [generated1]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = 5
    ny<<<{"description": "Number of elements in the Y direction"}>>> = 15
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = -0.6
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = -0.1
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = 5
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = 0.9
    boundary_name_prefix<<<{"description": "If provided, prefix the built in boundary names with this string"}>>> = pillar1
  []
  [generated2]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = 5
    ny<<<{"description": "Number of elements in the Y direction"}>>> = 15
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = 0.1
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = 0.6
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = 5
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = 0.9
    boundary_name_prefix<<<{"description": "If provided, prefix the built in boundary names with this string"}>>> = pillar2
    boundary_id_offset<<<{"description": "This offset is added to the generated boundary IDs"}>>> = 4
  []
  [collect_meshes]
    type = MeshCollectionGenerator<<<{"description": "Collects multiple meshes into a single (unconnected) mesh.", "href": "../../../../source/meshgenerators/MeshCollectionGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'generated1 generated2'
  []

  patch_update_strategy = iteration
[]

[Physics<<<{"href": "../../../../syntax/Physics/index.html"}>>>/SolidMechanics<<<{"href": "../../../../syntax/Physics/SolidMechanics/index.html"}>>>/QuasiStatic<<<{"href": "../../../../syntax/Physics/SolidMechanics/QuasiStatic/index.html"}>>>]
  [all]
    add_variables<<<{"description": "Add the displacement variables"}>>> = true
    strain<<<{"description": "Strain formulation"}>>> = FINITE
    generate_output<<<{"description": "Add scalar quantity output for stress and/or strain"}>>> = 'vonmises_stress'
  []
[]

[Contact<<<{"href": "../../../../syntax/Contact/index.html"}>>>]
  [pillars]
    primary<<<{"description": "The list of boundary IDs referring to primary sidesets"}>>> = pillar1_right
    secondary<<<{"description": "The list of boundary IDs referring to secondary sidesets"}>>> = pillar2_left
    model<<<{"description": "The contact model to use"}>>> = frictionless
    formulation<<<{"description": "The contact formulation"}>>> = penalty
    penalty<<<{"description": "The penalty to apply.  This can vary depending on the stiffness of your materials"}>>> = 1e9
    normalize_penalty<<<{"description": "Whether to normalize the penalty parameter with the nodal area."}>>> = true
  []
[]

[BCs<<<{"href": "../../../../syntax/BCs/index.html"}>>>]
  [bottom_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"}>>> = 'pillar1_bottom pillar2_bottom'
    value<<<{"description": "Value of the BC"}>>> = 0
  []
  [bottom_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"}>>> = 'pillar1_bottom pillar2_bottom'
    value<<<{"description": "Value of the BC"}>>> = 0
  []
  [Pressure<<<{"href": "../../../../syntax/BCs/Pressure/index.html"}>>>]
    [sides]
      boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 'pillar1_left pillar2_right'
      # we square time here to get a more progressive loading curve
      # (more pressure later on once contact is established)
      function<<<{"description": "The function that describes the pressure"}>>> = 1e4*t^2
    []
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [elasticity]
    type = ComputeIsotropicElasticityTensor<<<{"description": "Compute a constant isotropic elasticity tensor.", "href": "../../../../source/materials/ComputeIsotropicElasticityTensor.html"}>>>
    youngs_modulus<<<{"description": "Young's modulus of the material."}>>> = 1e9
    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"}>>>
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = NEWTON
  line_search = none
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  end_time = 5
  dt = 0.5
  [Predictor<<<{"href": "../../../../syntax/Executioner/Predictor/index.html"}>>>]
    type = SimplePredictor
    scale = 1
  []
[]

[Outputs<<<{"href": "../../../../syntax/Mastodon/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
  print_linear_residuals<<<{"description": "Enable printing of linear residuals to the screen (Console)"}>>> = false
  perf_graph<<<{"description": "Enable printing of the performance graph to the screen (Console)"}>>> = true
[]
(moose/modules/contact/tutorials/introduction/step01.i)

If you recall, in the final step of the mechanics tutorial we modeled two cantilevers that were bent towards each other and ended up occupying the same space without interacting. We're about to change that.

Input file

Mesh

We add the "patch_update_strategy" parameter and set it to iteration. MOOSE actually recommends this setting when you run the model without it. The parameter configures the geometric search that is required for modeling contact.

Contact

This is the major functional addition to the previous input. The [Contact] action block is doing the heavy lifting for us in setting up all objects required to enforce mechanical contact in one of the numerous ways supported by MOOSE (penalty, kinematic, mortar - frictionless, glued, frictional).

We use the "primary" and "secondary" parameters to specify the two interfaces we want to interact through mechanical contact. The "model" parameter is set to select frictionless contact and as the "formulation" we chose the penalty method. Frictionless penalty based contact is a good initial choice for contact modeling as it exhibits benign convergence properties and works in 2D as well as 3D.

We select a "penalty" factor of 1e9. The choice of the penalty should be guided by the stiffness of your system. A factor close to the Young's modulus of the materials involved in contact, but not exceeding it, is recommended. Penalty contact allows for a non-zero interpenetration of the contact surfaces. Both the kinematic and mortar formulation will enforce penetrations as low as the solution convergence allows.

We utilize the "normalize_penalty" option here to compensate for mesh size effects.

BCs

Note that we've made a small change to the pressure BC, changing the applied pressure from 1e4*t to 1e4*t^2. This is simply for demonstration purposes to apply a stronger pressure later in the simulation to force the cantilevers to bend more and establish contact over a larger area.

Tasks and questions

First run the example and look at the output. You should see some new fields that can be visualized in Paraview (or the Exodus viewer of your choice). In particular please pay attention to contact_force and penetration.

Penetration

Visualize the penetration variable. You may have to rescale the visualization to start the color scale at 0. Negative penetrations are not of interest here (they effectively are the gap width)

You should see a maximum interpenetration of about 6e-4, which is about half a percent of the element width. It is application dependent whether this amount is acceptable.

Once you've answered the questions and run this example we will move on to Step 2 which introduces mortar based contact.