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.