Example 6c: Force-controlled frictional contact problem.

Model Description

This example demonstrates frictional contact between a solid block and elastic soil in MASTODON. A 6 inch cubic block is centered on top of a 4 foot cubic block of soil. The top surface of the soil is free and the remaining surfaces are fixed in all directions. A uniform normal pressure of 5 psi is applied to the top of the smaller solid block and force is applied to the block in direction parallel to the contact surface. The resulting normal and frictional forces at the interface of the two materials are then obtained and compared with analytical results.

Modeling in MASTODON

The simulation is carried out in two steps.

Step 1: Stabilizing the system under the gravity and pressure forces

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  type = FileMesh
  file = contact_modified.e
  patch_update_strategy = iteration
  patch_size = 40
  partitioner = centroid
  centroid_partitioner_direction = y
[]

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

[Variables<<<{"href": "../syntax/Variables/index.html"}>>>]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [nor_forc]
    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
  []
  [tang_forc_x]
    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
  []
[]

[Kernels<<<{"href": "../syntax/Kernels/index.html"}>>>]
  [TensorMechanics<<<{"href": "../syntax/Kernels/TensorMechanics/index.html"}>>>]
    use_displaced_mesh<<<{"description": "Whether to use displaced mesh in the kernels"}>>> = true
    displacements<<<{"description": "The nonlinear displacement variables for the problem"}>>> = 'disp_x disp_y disp_z'
  []
  [gravity]
    type = Gravity<<<{"description": "Apply gravity. Value is in units of acceleration.", "href": "../source/kernels/Gravity.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_z
    value<<<{"description": "Value multiplied against the residual, e.g. gravitational acceleration"}>>> = -386.09 #in/s2
  []
[]

[AuxKernels<<<{"href": "../syntax/AuxKernels/index.html"}>>>]
  [nor_forc]
    type = PenetrationAux<<<{"description": "Auxiliary Kernel for computing several geometry related quantities between two contacting bodies.", "href": "../source/auxkernels/PenetrationAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = nor_forc
    quantity<<<{"description": "The quantity to recover from the available penetration information"}>>> = normal_force_magnitude
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 102
    paired_boundary<<<{"description": "The boundary to be penetrated"}>>> = 103
  []
  [tang_forc_x]
    type = PenetrationAux<<<{"description": "Auxiliary Kernel for computing several geometry related quantities between two contacting bodies.", "href": "../source/auxkernels/PenetrationAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = tang_forc_x
    quantity<<<{"description": "The quantity to recover from the available penetration information"}>>> = tangential_force_x
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 102
    paired_boundary<<<{"description": "The boundary to be penetrated"}>>> = 103
  []
[]

[BCs<<<{"href": "../syntax/BCs/index.html"}>>>]
  [fix_x_soil]
    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"}>>> = 100
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [fix_y_soil]
    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"}>>> = 100
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [fix_z_soil]
    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_z
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 100
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [concrete_pressure]
    type = Pressure<<<{"description": "Applies a pressure on a given boundary in a given direction", "href": "../source/bcs/Pressure.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 101
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_z
    component<<<{"description": "The component for the pressure"}>>> = 2
    factor<<<{"description": "The magnitude to use in computing the pressure"}>>> = 5 #psi
  []
[]

[Materials<<<{"href": "../syntax/Materials/index.html"}>>>]
  [elasticity_tensor_block]
    youngs_modulus<<<{"description": "Young's modulus of the material."}>>> = 4e6 #psi
    poissons_ratio<<<{"description": "Poisson's ratio for the material."}>>> = 0.25
    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"}>>> = 2
  []
  [strain_block]
    type = ComputeFiniteStrain<<<{"description": "Compute a strain increment and rotation increment for finite strains.", "href": "../source/materials/ComputeFiniteStrain.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 2
    displacements<<<{"description": "The displacements appropriate for the simulation geometry and coordinate system"}>>> = 'disp_x disp_y disp_z'
  []
  [stress_block]
    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"}>>> = 2
  []
  [den_block]
    type = GenericConstantMaterial<<<{"description": "Declares material properties based on names and values prescribed by input parameters.", "href": "../source/materials/GenericConstantMaterial.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 2
    prop_names<<<{"description": "The names of the properties this material will have"}>>> = density
    prop_values<<<{"description": "The values associated with the named properties"}>>> = 0.0002248 #slug/in^3
  []
  [elasticity_tensor_soil]
    youngs_modulus<<<{"description": "Young's modulus of the material."}>>> = 1.3983e+05 #psi
    poissons_ratio<<<{"description": "Poisson's ratio for the material."}>>> = 0.3
    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
  []
  [strain_soil]
    type = ComputeFiniteStrain<<<{"description": "Compute a strain increment and rotation increment for finite strains.", "href": "../source/materials/ComputeFiniteStrain.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 1
    displacements<<<{"description": "The displacements appropriate for the simulation geometry and coordinate system"}>>> = 'disp_x disp_y disp_z'
  []
  [stress_soil]
    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
  []
  [den_soil]
    type = GenericConstantMaterial<<<{"description": "Declares material properties based on names and values prescribed by input parameters.", "href": "../source/materials/GenericConstantMaterial.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 1
    prop_names<<<{"description": "The names of the properties this material will have"}>>> = density
    prop_values<<<{"description": "The values associated with the named properties"}>>> = 0.0001356 #slug/in^3
  []
[]

[Contact<<<{"href": "../syntax/Contact/index.html"}>>>]
  [leftright]
    secondary<<<{"description": "The list of boundary IDs referring to secondary sidesets"}>>> = 102
    primary<<<{"description": "The list of boundary IDs referring to primary sidesets"}>>> = 103
    model<<<{"description": "The contact model to use"}>>> = coulomb
    formulation<<<{"description": "The contact formulation"}>>> = penalty
    normalize_penalty<<<{"description": "Whether to normalize the penalty parameter with the nodal area."}>>> = true
    friction_coefficient<<<{"description": "The friction coefficient"}>>> = 0.7
    penalty<<<{"description": "The penalty to apply.  This can vary depending on the stiffness of your materials"}>>> = 1e5
    displacements<<<{"description": "The displacements appropriate for the simulation geometry and coordinate system"}>>> = 'disp_x disp_y disp_z'
  []
[]

[Preconditioning<<<{"href": "../syntax/Preconditioning/index.html"}>>>]
  [andy]
    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
  []
[]

[Postprocessors<<<{"href": "../syntax/Postprocessors/index.html"}>>>]
  [nor_forc]
    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"}>>> = nor_forc
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 102
  []
  [tang_forc_x]
    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"}>>> = tang_forc_x
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 102
  []
  [dispx]
    type = NodalExtremeValue<<<{"description": "Finds either the min or max elemental value of a variable over the domain.", "href": "../source/postprocessors/NodalExtremeValue.html"}>>>
    value_type<<<{"description": "Type of extreme value to return. 'max' returns the maximum value. 'min' returns the minimum value. 'max_abs' returns the maximum of the absolute value."}>>> = max
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = '2'
    variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = disp_x
  []
[]

[Executioner<<<{"href": "../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu     superlu_dist'
  line_search = 'none'
  end_time = 2.0
  dt = 0.005
  dtmin = 0.001
  nl_abs_tol = 1e-1
  nl_rel_tol = 1e-4
  l_tol = 1e-2
  l_max_its = 20
  timestep_tolerance = 1e-3
[]

[Outputs<<<{"href": "../syntax/Mastodon/Outputs/index.html"}>>>]
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
  perf_graph<<<{"description": "Enable printing of the performance graph to the screen (Console)"}>>> = 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)"}>>> = 1
  []
  [out1]
    type = Checkpoint<<<{"description": "Output for MOOSE recovery checkpoint files.", "href": "../source/outputs/Checkpoint.html"}>>>
    interval = 1
    num_files<<<{"description": "Number of the restart files to save"}>>> = 2
  []
[]

A static analysis is carried out in the mesh shown in Figure 1 in the first step to stabilize the system under the action of constant gravitational force and pressure force as explained in Example 6b

Figure 1: Input model in MASTODON.

Step 2: Performing the dynamic simulation using the stabilized system

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  type = FileMesh
  file = ../ex06c/stabilize_lc_out1_cp/LATEST
  patch_update_strategy = iteration
  patch_size = 40
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Problem<<<{"href": "../syntax/Problem/index.html"}>>>]
  restart_file_base = ../ex06c/stabilize_lc_out1_cp/LATEST
[]

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

[Variables<<<{"href": "../syntax/Variables/index.html"}>>>]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [vel_x]
  []
  [accel_x]
  []
  [accel_y]
  []
  [vel_y]
  []
  [vel_z]
  []
  [accel_z]
  []
  [nor_forc]
    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
  []
  [tang_forc_x]
    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
  []
[]

[Kernels<<<{"href": "../syntax/Kernels/index.html"}>>>]
  [TensorMechanics<<<{"href": "../syntax/Kernels/TensorMechanics/index.html"}>>>]
    use_displaced_mesh<<<{"description": "Whether to use displaced mesh in the kernels"}>>> = true
    displacements<<<{"description": "The nonlinear displacement variables for the problem"}>>> = 'disp_x disp_y disp_z'
  []
  [gravity]
    type = Gravity<<<{"description": "Apply gravity. Value is in units of acceleration.", "href": "../source/kernels/Gravity.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_z
    value<<<{"description": "Value multiplied against the residual, e.g. gravitational acceleration"}>>> = -386.09 #in/s2
  []
  [inertia_x]
    type = InertialForce<<<{"description": "Calculates the residual for the inertial force ($M \\cdot acceleration$) and the contribution of mass dependent Rayleigh damping and HHT time  integration scheme ($\\eta \\cdot M \\cdot ((1+\\alpha)velq2-\\alpha \\cdot vel-old) $)", "href": "../source/kernels/InertialForce.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    velocity<<<{"description": "velocity variable"}>>> = vel_x
    acceleration<<<{"description": "acceleration variable"}>>> = accel_x
    beta<<<{"description": "beta parameter for Newmark Time integration"}>>> = 0.25
    gamma<<<{"description": "gamma parameter for Newmark Time integration"}>>> = 0.5
    eta<<<{"description": "Name of material property or a constant real number defining the eta parameter for the Rayleigh damping."}>>> = 0.0
  []
  [inertia_y]
    type = InertialForce<<<{"description": "Calculates the residual for the inertial force ($M \\cdot acceleration$) and the contribution of mass dependent Rayleigh damping and HHT time  integration scheme ($\\eta \\cdot M \\cdot ((1+\\alpha)velq2-\\alpha \\cdot vel-old) $)", "href": "../source/kernels/InertialForce.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    velocity<<<{"description": "velocity variable"}>>> = vel_y
    acceleration<<<{"description": "acceleration variable"}>>> = accel_y
    beta<<<{"description": "beta parameter for Newmark Time integration"}>>> = 0.25
    gamma<<<{"description": "gamma parameter for Newmark Time integration"}>>> = 0.5
    eta<<<{"description": "Name of material property or a constant real number defining the eta parameter for the Rayleigh damping."}>>> = 0.0
  []
  [inertia_z]
    type = InertialForce<<<{"description": "Calculates the residual for the inertial force ($M \\cdot acceleration$) and the contribution of mass dependent Rayleigh damping and HHT time  integration scheme ($\\eta \\cdot M \\cdot ((1+\\alpha)velq2-\\alpha \\cdot vel-old) $)", "href": "../source/kernels/InertialForce.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_z
    velocity<<<{"description": "velocity variable"}>>> = vel_z
    acceleration<<<{"description": "acceleration variable"}>>> = accel_z
    beta<<<{"description": "beta parameter for Newmark Time integration"}>>> = 0.25
    gamma<<<{"description": "gamma parameter for Newmark Time integration"}>>> = 0.5
    eta<<<{"description": "Name of material property or a constant real number defining the eta parameter for the Rayleigh damping."}>>> = 0.0
  []
[]

[AuxKernels<<<{"href": "../syntax/AuxKernels/index.html"}>>>]
  [nor_forc]
    type = PenetrationAux<<<{"description": "Auxiliary Kernel for computing several geometry related quantities between two contacting bodies.", "href": "../source/auxkernels/PenetrationAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = nor_forc
    quantity<<<{"description": "The quantity to recover from the available penetration information"}>>> = normal_force_magnitude
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 102
    paired_boundary<<<{"description": "The boundary to be penetrated"}>>> = 103
  []
  [tang_forc_x]
    type = PenetrationAux<<<{"description": "Auxiliary Kernel for computing several geometry related quantities between two contacting bodies.", "href": "../source/auxkernels/PenetrationAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = tang_forc_x
    quantity<<<{"description": "The quantity to recover from the available penetration information"}>>> = tangential_force_x
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 102
    paired_boundary<<<{"description": "The boundary to be penetrated"}>>> = 103
  []
  [accel_x]
    type = NewmarkAccelAux<<<{"description": "Computes the current acceleration using the Newmark method.", "href": "../source/auxkernels/NewmarkAccelAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = accel_x
    displacement<<<{"description": "displacement variable"}>>> = disp_x
    velocity<<<{"description": "velocity variable"}>>> = vel_x
    beta<<<{"description": "beta parameter for Newmark method"}>>> = 0.25
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [vel_x]
    type = NewmarkVelAux<<<{"description": "Calculates the current velocity using Newmark method.", "href": "../source/auxkernels/NewmarkVelAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = vel_x
    acceleration<<<{"description": "acceleration variable"}>>> = accel_x
    gamma<<<{"description": "gamma parameter for Newmark method"}>>> = 0.5
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [accel_y]
    type = NewmarkAccelAux<<<{"description": "Computes the current acceleration using the Newmark method.", "href": "../source/auxkernels/NewmarkAccelAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = accel_y
    displacement<<<{"description": "displacement variable"}>>> = disp_y
    velocity<<<{"description": "velocity variable"}>>> = vel_y
    beta<<<{"description": "beta parameter for Newmark method"}>>> = 0.25
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [vel_y]
    type = NewmarkVelAux<<<{"description": "Calculates the current velocity using Newmark method.", "href": "../source/auxkernels/NewmarkVelAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = vel_y
    acceleration<<<{"description": "acceleration variable"}>>> = accel_y
    gamma<<<{"description": "gamma parameter for Newmark method"}>>> = 0.5
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [accel_z]
    type = NewmarkAccelAux<<<{"description": "Computes the current acceleration using the Newmark method.", "href": "../source/auxkernels/NewmarkAccelAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = accel_z
    displacement<<<{"description": "displacement variable"}>>> = disp_z
    velocity<<<{"description": "velocity variable"}>>> = vel_z
    beta<<<{"description": "beta parameter for Newmark method"}>>> = 0.25
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [vel_z]
    type = NewmarkVelAux<<<{"description": "Calculates the current velocity using Newmark method.", "href": "../source/auxkernels/NewmarkVelAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = vel_z
    acceleration<<<{"description": "acceleration variable"}>>> = accel_z
    gamma<<<{"description": "gamma parameter for Newmark method"}>>> = 0.5
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
[]

[BCs<<<{"href": "../syntax/BCs/index.html"}>>>]
  [fix_x_soil]
    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"}>>> = 100
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [fix_y_soil]
    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"}>>> = 100
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [fix_z_soil]
    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_z
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 100
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [concrete_pressure]
    type = Pressure<<<{"description": "Applies a pressure on a given boundary in a given direction", "href": "../source/bcs/Pressure.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 101
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_z
    component<<<{"description": "The component for the pressure"}>>> = 2
    factor<<<{"description": "The magnitude to use in computing the pressure"}>>> = 5 #psi
  []
[]

[NodalKernels<<<{"href": "../syntax/NodalKernels/index.html"}>>>]
  [force_x]
    type = UserForcingFunctionNodalKernel<<<{"description": "Residual contribution to an ODE from a source function acting at nodes.", "href": "../source/nodalkernels/UserForcingFunctionNodalKernel.html"}>>>
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 1000
    function<<<{"description": "The forcing function"}>>> = force
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
  []
[]

[Functions<<<{"href": "../syntax/Functions/index.html"}>>>]
  [force]
    type = PiecewiseLinear<<<{"description": "Linearly interpolates between pairs of x-y data", "href": "../source/functions/PiecewiseLinear.html"}>>>
    x<<<{"description": "The abscissa values"}>>> = '0 1 2'
    y<<<{"description": "The ordinate values"}>>> = '0  0.1953125 0.390625' #number of nodes = 512, force 0,100,200
  []
[]

[Materials<<<{"href": "../syntax/Materials/index.html"}>>>]
  [elasticity_tensor_block]
    youngs_modulus<<<{"description": "Young's modulus of the material."}>>> = 4e6 #psi
    poissons_ratio<<<{"description": "Poisson's ratio for the material."}>>> = 0.25
    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"}>>> = 2
  []
  [strain_block]
    type = ComputeFiniteStrain<<<{"description": "Compute a strain increment and rotation increment for finite strains.", "href": "../source/materials/ComputeFiniteStrain.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 2
    displacements<<<{"description": "The displacements appropriate for the simulation geometry and coordinate system"}>>> = 'disp_x disp_y disp_z'
  []
  [stress_block]
    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"}>>> = 2
  []
  [den_block]
    type = GenericConstantMaterial<<<{"description": "Declares material properties based on names and values prescribed by input parameters.", "href": "../source/materials/GenericConstantMaterial.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 2
    prop_names<<<{"description": "The names of the properties this material will have"}>>> = density
    prop_values<<<{"description": "The values associated with the named properties"}>>> = 0.0002248 #slug/in^3
  []
  [elasticity_tensor_soil]
    youngs_modulus<<<{"description": "Young's modulus of the material."}>>> = 1.3983e+05 #psi
    poissons_ratio<<<{"description": "Poisson's ratio for the material."}>>> = 0.3
    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
  []
  [strain_soil]
    type = ComputeFiniteStrain<<<{"description": "Compute a strain increment and rotation increment for finite strains.", "href": "../source/materials/ComputeFiniteStrain.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 1
    displacements<<<{"description": "The displacements appropriate for the simulation geometry and coordinate system"}>>> = 'disp_x disp_y disp_z'
  []
  [stress_soil]
    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
  []
  [den_soil]
    type = GenericConstantMaterial<<<{"description": "Declares material properties based on names and values prescribed by input parameters.", "href": "../source/materials/GenericConstantMaterial.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 1
    prop_names<<<{"description": "The names of the properties this material will have"}>>> = density
    prop_values<<<{"description": "The values associated with the named properties"}>>> = 0.0001356 #slug/in^3
  []
[]

[Contact<<<{"href": "../syntax/Contact/index.html"}>>>]
  [leftright]
    secondary<<<{"description": "The list of boundary IDs referring to secondary sidesets"}>>> = 102
    primary<<<{"description": "The list of boundary IDs referring to primary sidesets"}>>> = 103
    model<<<{"description": "The contact model to use"}>>> = coulomb
    formulation<<<{"description": "The contact formulation"}>>> = penalty
    normalize_penalty<<<{"description": "Whether to normalize the penalty parameter with the nodal area."}>>> = true
    friction_coefficient<<<{"description": "The friction coefficient"}>>> = 0.7
    penalty<<<{"description": "The penalty to apply.  This can vary depending on the stiffness of your materials"}>>> = 1e5
    displacements<<<{"description": "The displacements appropriate for the simulation geometry and coordinate system"}>>> = 'disp_x disp_y disp_z'
  []
[]

[Preconditioning<<<{"href": "../syntax/Preconditioning/index.html"}>>>]
  [andy]
    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
  []
[]

[Postprocessors<<<{"href": "../syntax/Postprocessors/index.html"}>>>]
  [nor_forc]
    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"}>>> = nor_forc
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 102
  []
  [tang_forc_x]
    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"}>>> = tang_forc_x
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 102
  []
  [dispx]
    type = NodalExtremeValue<<<{"description": "Finds either the min or max elemental value of a variable over the domain.", "href": "../source/postprocessors/NodalExtremeValue.html"}>>>
    value_type<<<{"description": "Type of extreme value to return. 'max' returns the maximum value. 'min' returns the minimum value. 'max_abs' returns the maximum of the absolute value."}>>> = max
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = '2'
    variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = disp_x
  []
  [force2x]
    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."}>>> = force
    scale_factor<<<{"description": "A scale factor to be applied to the function"}>>> = 512 #Number of nodes
  []
[]

[Executioner<<<{"href": "../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu     superlu_dist'
  line_search = 'none'
  start_time = 0
  end_time = 1.55
  dt = 0.005
  dtmin = 0.001
  nl_abs_tol = 1e-1
  nl_rel_tol = 1e-4
  l_tol = 1e-2
  l_max_its = 20
  timestep_tolerance = 1e-3
[]

[Outputs<<<{"href": "../syntax/Mastodon/Outputs/index.html"}>>>]
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = results55
  perf_graph<<<{"description": "Enable printing of the performance graph to the screen (Console)"}>>> = true
  print_linear_residuals<<<{"description": "Enable printing of linear residuals to the screen (Console)"}>>> = false
  [screen]
    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)"}>>> = 1
  []
[]

Advanced restart is used to continue the simulation as discussed in Example 6b. A linearly ramping force from 0 lbf to 155 lbf is applied to the solid block.

Results

The graph for the frictional force as a function of displacement obtained from MASTODON is shown in Figure 2.

Figure 2: Graph showing the frictional force as a function of applied force for Columb Friction Model

Theoretical Solution

Normal Force = density x volume x gravity + pressure x area = 198.7473 lbf

Maximum Frictional Force = coefficient of friction x Normal Force = 139.123 lbf

The frictional force increases linearly from 0 to 139.123 and remains constant after the block starts sliding. The same behavior is seen in the Figure 2 obtained from MASTODON.

(examples/ex06c/stabilize_lc.i)

[Mesh]
  type = FileMesh
  file = contact_modified.e
  patch_update_strategy = iteration
  patch_size = 40
  partitioner = centroid
  centroid_partitioner_direction = y
[]

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

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

[AuxVariables]
  [nor_forc]
    order = FIRST
    family = LAGRANGE
  []
  [tang_forc_x]
    order = FIRST
    family = LAGRANGE
  []
[]

[Kernels]
  [TensorMechanics]
    use_displaced_mesh = true
    displacements = 'disp_x disp_y disp_z'
  []
  [gravity]
    type = Gravity
    variable = disp_z
    value = -386.09 #in/s2
  []
[]

[AuxKernels]
  [nor_forc]
    type = PenetrationAux
    variable = nor_forc
    quantity = normal_force_magnitude
    boundary = 102
    paired_boundary = 103
  []
  [tang_forc_x]
    type = PenetrationAux
    variable = tang_forc_x
    quantity = tangential_force_x
    boundary = 102
    paired_boundary = 103
  []
[]

[BCs]
  [fix_x_soil]
    type = DirichletBC
    variable = disp_x
    boundary = 100
    value = 0.0
  []
  [fix_y_soil]
    type = DirichletBC
    variable = disp_y
    boundary = 100
    value = 0.0
  []
  [fix_z_soil]
    type = DirichletBC
    variable = disp_z
    boundary = 100
    value = 0.0
  []
  [concrete_pressure]
    type = Pressure
    boundary = 101
    variable = disp_z
    component = 2
    factor = 5 #psi
  []
[]

[Materials]
  [elasticity_tensor_block]
    youngs_modulus = 4e6 #psi
    poissons_ratio = 0.25
    type = ComputeIsotropicElasticityTensor
    block = 2
  []
  [strain_block]
    type = ComputeFiniteStrain
    block = 2
    displacements = 'disp_x disp_y disp_z'
  []
  [stress_block]
    type = ComputeFiniteStrainElasticStress
    block = 2
  []
  [den_block]
    type = GenericConstantMaterial
    block = 2
    prop_names = density
    prop_values = 0.0002248 #slug/in^3
  []
  [elasticity_tensor_soil]
    youngs_modulus = 1.3983e+05 #psi
    poissons_ratio = 0.3
    type = ComputeIsotropicElasticityTensor
    block = 1
  []
  [strain_soil]
    type = ComputeFiniteStrain
    block = 1
    displacements = 'disp_x disp_y disp_z'
  []
  [stress_soil]
    type = ComputeFiniteStrainElasticStress
    block = 1
  []
  [den_soil]
    type = GenericConstantMaterial
    block = 1
    prop_names = density
    prop_values = 0.0001356 #slug/in^3
  []
[]

[Contact]
  [leftright]
    secondary = 102
    primary = 103
    model = coulomb
    formulation = penalty
    normalize_penalty = true
    friction_coefficient = 0.7
    penalty = 1e5
    displacements = 'disp_x disp_y disp_z'
  []
[]

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

[Postprocessors]
  [nor_forc]
    type = NodalSum
    variable = nor_forc
    boundary = 102
  []
  [tang_forc_x]
    type = NodalSum
    variable = tang_forc_x
    boundary = 102
  []
  [dispx]
    type = NodalExtremeValue
    value_type = max
    block = '2'
    variable = disp_x
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu     superlu_dist'
  line_search = 'none'
  end_time = 2.0
  dt = 0.005
  dtmin = 0.001
  nl_abs_tol = 1e-1
  nl_rel_tol = 1e-4
  l_tol = 1e-2
  l_max_its = 20
  timestep_tolerance = 1e-3
[]

[Outputs]
  csv = true
  exodus = true
  perf_graph = true
  [console]
    type = Console
    max_rows = 1
  []
  [out1]
    type = Checkpoint
    interval = 1
    num_files = 2
  []
[]

(examples/ex06c/actual_simulation_lc.i)

[Mesh]
  type = FileMesh
  file = ../ex06c/stabilize_lc_out1_cp/LATEST
  patch_update_strategy = iteration
  patch_size = 40
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Problem]
  restart_file_base = ../ex06c/stabilize_lc_out1_cp/LATEST
[]

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

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

[AuxVariables]
  [vel_x]
  []
  [accel_x]
  []
  [accel_y]
  []
  [vel_y]
  []
  [vel_z]
  []
  [accel_z]
  []
  [nor_forc]
    order = FIRST
    family = LAGRANGE
  []
  [tang_forc_x]
    order = FIRST
    family = LAGRANGE
  []
[]

[Kernels]
  [TensorMechanics]
    use_displaced_mesh = true
    displacements = 'disp_x disp_y disp_z'
  []
  [gravity]
    type = Gravity
    variable = disp_z
    value = -386.09 #in/s2
  []
  [inertia_x]
    type = InertialForce
    variable = disp_x
    velocity = vel_x
    acceleration = accel_x
    beta = 0.25
    gamma = 0.5
    eta = 0.0
  []
  [inertia_y]
    type = InertialForce
    variable = disp_y
    velocity = vel_y
    acceleration = accel_y
    beta = 0.25
    gamma = 0.5
    eta = 0.0
  []
  [inertia_z]
    type = InertialForce
    variable = disp_z
    velocity = vel_z
    acceleration = accel_z
    beta = 0.25
    gamma = 0.5
    eta = 0.0
  []
[]

[AuxKernels]
  [nor_forc]
    type = PenetrationAux
    variable = nor_forc
    quantity = normal_force_magnitude
    boundary = 102
    paired_boundary = 103
  []
  [tang_forc_x]
    type = PenetrationAux
    variable = tang_forc_x
    quantity = tangential_force_x
    boundary = 102
    paired_boundary = 103
  []
  [accel_x]
    type = NewmarkAccelAux
    variable = accel_x
    displacement = disp_x
    velocity = vel_x
    beta = 0.25
    execute_on = timestep_end
  []
  [vel_x]
    type = NewmarkVelAux
    variable = vel_x
    acceleration = accel_x
    gamma = 0.5
    execute_on = timestep_end
  []
  [accel_y]
    type = NewmarkAccelAux
    variable = accel_y
    displacement = disp_y
    velocity = vel_y
    beta = 0.25
    execute_on = timestep_end
  []
  [vel_y]
    type = NewmarkVelAux
    variable = vel_y
    acceleration = accel_y
    gamma = 0.5
    execute_on = timestep_end
  []
  [accel_z]
    type = NewmarkAccelAux
    variable = accel_z
    displacement = disp_z
    velocity = vel_z
    beta = 0.25
    execute_on = timestep_end
  []
  [vel_z]
    type = NewmarkVelAux
    variable = vel_z
    acceleration = accel_z
    gamma = 0.5
    execute_on = timestep_end
  []
[]

[BCs]
  [fix_x_soil]
    type = DirichletBC
    variable = disp_x
    boundary = 100
    value = 0.0
  []
  [fix_y_soil]
    type = DirichletBC
    variable = disp_y
    boundary = 100
    value = 0.0
  []
  [fix_z_soil]
    type = DirichletBC
    variable = disp_z
    boundary = 100
    value = 0.0
  []
  [concrete_pressure]
    type = Pressure
    boundary = 101
    variable = disp_z
    component = 2
    factor = 5 #psi
  []
[]

[NodalKernels]
  [force_x]
    type = UserForcingFunctionNodalKernel
    boundary = 1000
    function = force
    variable = disp_x
  []
[]

[Functions]
  [force]
    type = PiecewiseLinear
    x = '0 1 2'
    y = '0  0.1953125 0.390625' #number of nodes = 512, force 0,100,200
  []
[]

[Materials]
  [elasticity_tensor_block]
    youngs_modulus = 4e6 #psi
    poissons_ratio = 0.25
    type = ComputeIsotropicElasticityTensor
    block = 2
  []
  [strain_block]
    type = ComputeFiniteStrain
    block = 2
    displacements = 'disp_x disp_y disp_z'
  []
  [stress_block]
    type = ComputeFiniteStrainElasticStress
    block = 2
  []
  [den_block]
    type = GenericConstantMaterial
    block = 2
    prop_names = density
    prop_values = 0.0002248 #slug/in^3
  []
  [elasticity_tensor_soil]
    youngs_modulus = 1.3983e+05 #psi
    poissons_ratio = 0.3
    type = ComputeIsotropicElasticityTensor
    block = 1
  []
  [strain_soil]
    type = ComputeFiniteStrain
    block = 1
    displacements = 'disp_x disp_y disp_z'
  []
  [stress_soil]
    type = ComputeFiniteStrainElasticStress
    block = 1
  []
  [den_soil]
    type = GenericConstantMaterial
    block = 1
    prop_names = density
    prop_values = 0.0001356 #slug/in^3
  []
[]

[Contact]
  [leftright]
    secondary = 102
    primary = 103
    model = coulomb
    formulation = penalty
    normalize_penalty = true
    friction_coefficient = 0.7
    penalty = 1e5
    displacements = 'disp_x disp_y disp_z'
  []
[]

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

[Postprocessors]
  [nor_forc]
    type = NodalSum
    variable = nor_forc
    boundary = 102
  []
  [tang_forc_x]
    type = NodalSum
    variable = tang_forc_x
    boundary = 102
  []
  [dispx]
    type = NodalExtremeValue
    value_type = max
    block = '2'
    variable = disp_x
  []
  [force2x]
    type = FunctionValuePostprocessor
    function = force
    scale_factor = 512 #Number of nodes
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu     superlu_dist'
  line_search = 'none'
  start_time = 0
  end_time = 1.55
  dt = 0.005
  dtmin = 0.001
  nl_abs_tol = 1e-1
  nl_rel_tol = 1e-4
  l_tol = 1e-2
  l_max_its = 20
  timestep_tolerance = 1e-3
[]

[Outputs]
  csv = true
  exodus = true
  file_base = results55
  perf_graph = true
  print_linear_residuals = false
  [screen]
    type = Console
    max_rows = 1
  []
[]

Model Description
Modeling in MASTODON
Step 1 : Stabilizing the system under the gravity and pressure forces
Step 2 : Performing the dynamic simulation using the stabilized system
Results
Theoretical Solution

Usage

Development

Demonstration