Example 6b: Displacement-controlled frictional contact problem with I-soil and using restart option for static initialization.

Model Description

This example demonstrates frictional contact between a solid block and I-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 it is given a prescribed displacement in a 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'
[]

[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
  []
  [layer_id]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
  []
[]

[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
  []
  [I_Soil<<<{"href": "../syntax/Materials/I_Soil/index.html"}>>>]
    [soil_1]
      soil_type<<<{"description": "This parameter determines the type of backbone curve used. Use 'user_defined' for a user defined backbone curve provided in a data file, 'darendeli' if the backbone curve is to be determined using Darandeli equations, 'gqh' if the backbone curve is determined using the GQ/H approach and 'thin_layer' if the soil is being used to simulate a thin-layer friction interface."}>>> = 'darendeli'
      layer_variable<<<{"description": "The auxvariable providing the soil layer identification."}>>> = layer_id
      layer_ids<<<{"description": "Vector of layer ids that map one-to-one to the rest of the soil layer parameters provided as input."}>>> = '0'
      over_consolidation_ratio<<<{"description": "The over consolidation ratio of the soil layers. Required for Darandeli backbone curve."}>>> = '1'
      plasticity_index<<<{"description": "The plasticity index of the soil layers. Required for Darandeli backbone curve."}>>> = '0'
      initial_shear_modulus<<<{"description": "The initial shear modulus of the soil layers."}>>> = '7251'
      number_of_points<<<{"description": "The total number of data points in which the backbone curve needs to be split for all soil layers (required for Darandeli or GQ/H type backbone curves)."}>>> = 10
      poissons_ratio<<<{"description": "Poissons's ratio for the soil layers. The size of the vector should be same as the size of layer_ids."}>>> = '0.3'
      block<<<{"description": "The blocks where this material model is applied."}>>> = 1
      density<<<{"description": "Vector of density values that map one-to-one with the number 'layer_ids' parameter."}>>> = '0.0001536'
      p_ref<<<{"description": "The reference pressure at which the parameters are defined for each soil layer. If 'soil_type = darendeli', then the reference pressure must be input in kilopascals."}>>> = '5'
      pressure_dependency<<<{"description": "Set to true to turn on pressure dependent stiffness and yield strength calculation."}>>> = false
    []
  []
[]

[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.2
    penalty<<<{"description": "The penalty to apply.  This can vary depending on the stiffness of your materials"}>>> = 1e4
    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.5
  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
  [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 first step to stabilize the system under the action of constant gravitational force and pressure force. The mesh shown in Figure 1 used in this model was generated externally using the mesh generation software, Cubit. MASTODON imports the mesh using the FileMesh type specified in the Mesh block.

The displacement variables are defined in the Variables block. The accelerations, velocities, normal and tangential (frictional) forces are defined as auxiliary variables in the AuxVariables block.

Figure 1: Input model in MASTODON

TensorMechanics kernel is used to set up the stress divergence kernels and Gravity kernel is used to apply the gravitational force. The ComputeIsotropicElasticityTensorBeam block is used to create the elasticity tensor of the elastic solid block, using Young's Modulus and Poisson's ratio. The stresses and strain are calculated using ComputeFiniteStrainElasticStress and ComputeFiniteStrain. The densities are assigned to the solid block and the soil using GenericConstantMaterial.

The I-soil block is used to create a nonlinear hysteretic soil model. The detailed explanation for this model can be found at I-soil.

In the soil block, the top surface is free and all other surfaces are fixed in every direction. A uniform vertical pressure of 5 psi is applied to the top surface of the solid block.

Contact is defined in MASTODON using the Contact block. This model is simulating Coulomb friction, with a coefficient of friction value equal to 0.2. The formulation implemented and demonstrated in this example uses a master/slave concept in which the nodes of the slave boundary are restricted by the nodes of the master surface boundary. Penetration is enforced using a penalty factor which is inversely related to the amount of penetration allowed (large penalty factor results in small amount of penetration and vice versa). The factor used in this example is 1e5.

[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.2
    penalty<<<{"description": "The penalty to apply.  This can vary depending on the stiffness of your materials"}>>> = 1e4
    displacements<<<{"description": "The displacements appropriate for the simulation geometry and coordinate system"}>>> = 'disp_x disp_y disp_z'
  []
[]

The use of restart in the MOOSE framework requires to specify checkpoints in the Outputs block. The checkpoint creates a set of files that be used as an initial condition for the next simulation.

[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
  [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
  []
[]

Step 2: Performing the dynamic simulation using the stabilized system

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  type = FileMesh
  file = ../ex06b/stabilize_isoil_out1_cp/LATEST #end file of previous simulation
  patch_update_strategy = iteration
  patch_size = 40
  partitioner = centroid
  centroid_partitioner_direction = y
[]

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

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

[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
  []
  [layer_id]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
  []
[]

[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
  []
  [layer_1]
    type = UniformLayerAuxKernel<<<{"description": "Computes an AuxVariable for representing a layered structure in an arbitrary direction.", "href": "../source/auxkernels/UniformLayerAuxKernel.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = layer_id
    interfaces<<<{"description": "A list of layer interface locations to apply across the domain in the specified direction."}>>> = '48'
    direction<<<{"description": "The direction to apply layering."}>>> = '0.0 0.0 1.0'
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 1
    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."}>>> = initial
  []
[]

[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
  []
  [top_x]
    type = PresetDisplacement<<<{"description": "Prescribe the displacement on a given boundary in a given direction.", "href": "../source/bcs/PresetDisplacement.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 1000
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    beta<<<{"description": "beta parameter for Newmark time integration."}>>> = 0.25
    velocity<<<{"description": "The velocity variable."}>>> = vel_x
    acceleration<<<{"description": "The acceleration variable."}>>> = accel_x
    function<<<{"description": "Function describing the displacement."}>>> = loading_bc
  []
[]

[Functions<<<{"href": "../syntax/Functions/index.html"}>>>]
  [loading_bc]
    type = PiecewiseLinear<<<{"description": "Linearly interpolates between pairs of x-y data", "href": "../source/functions/PiecewiseLinear.html"}>>>
    data_file<<<{"description": "File holding CSV data"}>>> = 'input_motion.csv'
    format<<<{"description": "Format of csv data file that is in either in columns or rows"}>>> = columns
  []
[]

[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"}>>>
    #  store_stress_old = true
    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
  []
  [I_Soil<<<{"href": "../syntax/Materials/I_Soil/index.html"}>>>]
    [soil_1]
      soil_type<<<{"description": "This parameter determines the type of backbone curve used. Use 'user_defined' for a user defined backbone curve provided in a data file, 'darendeli' if the backbone curve is to be determined using Darandeli equations, 'gqh' if the backbone curve is determined using the GQ/H approach and 'thin_layer' if the soil is being used to simulate a thin-layer friction interface."}>>> = 'darendeli'
      layer_variable<<<{"description": "The auxvariable providing the soil layer identification."}>>> = layer_id
      layer_ids<<<{"description": "Vector of layer ids that map one-to-one to the rest of the soil layer parameters provided as input."}>>> = '0'
      over_consolidation_ratio<<<{"description": "The over consolidation ratio of the soil layers. Required for Darandeli backbone curve."}>>> = '1'
      plasticity_index<<<{"description": "The plasticity index of the soil layers. Required for Darandeli backbone curve."}>>> = '0'
      initial_shear_modulus<<<{"description": "The initial shear modulus of the soil layers."}>>> = '7251'
      number_of_points<<<{"description": "The total number of data points in which the backbone curve needs to be split for all soil layers (required for Darandeli or GQ/H type backbone curves)."}>>> = 10
      poissons_ratio<<<{"description": "Poissons's ratio for the soil layers. The size of the vector should be same as the size of layer_ids."}>>> = '0.3'
      block<<<{"description": "The blocks where this material model is applied."}>>> = 1
      density<<<{"description": "Vector of density values that map one-to-one with the number 'layer_ids' parameter."}>>> = '0.0001536'
      p_ref<<<{"description": "The reference pressure at which the parameters are defined for each soil layer. If 'soil_type = darendeli', then the reference pressure must be input in kilopascals."}>>> = '5'
      pressure_dependency<<<{"description": "Set to true to turn on pressure dependent stiffness and yield strength calculation."}>>> = false
    []
  []
[]

[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.2
    penalty<<<{"description": "The penalty to apply.  This can vary depending on the stiffness of your materials"}>>> = 1e4
    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'
  start_time = 0
  end_time = 2.5
  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"}>>> = finalresult
  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
  []
[]

The final state of the previous simulation is used as an input file for the dynamic simulation using the Mesh. Problem block is used to define the advanced restart option as shown below.

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  type = FileMesh
  file = ../ex06b/stabilize_isoil_out1_cp/LATEST #end file of previous simulation
  patch_update_strategy = iteration
  patch_size = 40
  partitioner = centroid
  centroid_partitioner_direction = y
[]

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

The remaining input file is similar to the input file used in the previous step. The aux variables and the kernels required for the dynamic analysis are added and the input motion is applied to the block as shown in Figure 2.

Figure 2: Input displacement of the block

Results

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

Figure 3: Frictional force as a function of displacement for Coulomb friction model

Theoretical Solution

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

Frictional Force = coefficient of friction x Normal Force = 39.749 lbf

The results from MASTODON are in agreement with the theoretical calculations.

(examples/ex06b/stabilize_isoil.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'
[]

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

[AuxVariables]
  [nor_forc]
    order = FIRST
    family = LAGRANGE
  []
  [tang_forc_x]
    order = FIRST
    family = LAGRANGE
  []
  [layer_id]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[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
  []
  [I_Soil]
    [soil_1]
      soil_type = 'darendeli'
      layer_variable = layer_id
      layer_ids = '0'
      over_consolidation_ratio = '1'
      plasticity_index = '0'
      initial_shear_modulus = '7251'
      number_of_points = 10
      poissons_ratio = '0.3'
      block = 1
      density = '0.0001536'
      p_ref = '5'
      pressure_dependency = false
    []
  []
[]

[Contact]
  [leftright]
    secondary = 102
    primary = 103
    model = coulomb
    formulation = penalty
    normalize_penalty = true
    friction_coefficient = 0.2
    penalty = 1e4
    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.5
  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
  [console]
    type = Console
    max_rows = 1
  []
  [out1]
    type = Checkpoint
    interval = 1
    num_files = 2
  []
[]

(examples/ex06b/stabilize_isoil.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'
[]

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

[AuxVariables]
  [nor_forc]
    order = FIRST
    family = LAGRANGE
  []
  [tang_forc_x]
    order = FIRST
    family = LAGRANGE
  []
  [layer_id]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[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
  []
  [I_Soil]
    [soil_1]
      soil_type = 'darendeli'
      layer_variable = layer_id
      layer_ids = '0'
      over_consolidation_ratio = '1'
      plasticity_index = '0'
      initial_shear_modulus = '7251'
      number_of_points = 10
      poissons_ratio = '0.3'
      block = 1
      density = '0.0001536'
      p_ref = '5'
      pressure_dependency = false
    []
  []
[]

[Contact]
  [leftright]
    secondary = 102
    primary = 103
    model = coulomb
    formulation = penalty
    normalize_penalty = true
    friction_coefficient = 0.2
    penalty = 1e4
    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.5
  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
  [console]
    type = Console
    max_rows = 1
  []
  [out1]
    type = Checkpoint
    interval = 1
    num_files = 2
  []
[]

(examples/ex06b/stabilize_isoil.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'
[]

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

[AuxVariables]
  [nor_forc]
    order = FIRST
    family = LAGRANGE
  []
  [tang_forc_x]
    order = FIRST
    family = LAGRANGE
  []
  [layer_id]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[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
  []
  [I_Soil]
    [soil_1]
      soil_type = 'darendeli'
      layer_variable = layer_id
      layer_ids = '0'
      over_consolidation_ratio = '1'
      plasticity_index = '0'
      initial_shear_modulus = '7251'
      number_of_points = 10
      poissons_ratio = '0.3'
      block = 1
      density = '0.0001536'
      p_ref = '5'
      pressure_dependency = false
    []
  []
[]

[Contact]
  [leftright]
    secondary = 102
    primary = 103
    model = coulomb
    formulation = penalty
    normalize_penalty = true
    friction_coefficient = 0.2
    penalty = 1e4
    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.5
  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
  [console]
    type = Console
    max_rows = 1
  []
  [out1]
    type = Checkpoint
    interval = 1
    num_files = 2
  []
[]

(examples/ex06b/actual_simulation_isoil.i)

[Mesh]
  type = FileMesh
  file = ../ex06b/stabilize_isoil_out1_cp/LATEST #end file of previous simulation
  patch_update_strategy = iteration
  patch_size = 40
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Problem]
  restart_file_base = ../ex06b/stabilize_isoil_out1_cp/LATEST
[]

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

[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
  []
  [layer_id]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[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
  []
  [layer_1]
    type = UniformLayerAuxKernel
    variable = layer_id
    interfaces = '48'
    direction = '0.0 0.0 1.0'
    block = 1
    execute_on = initial
  []
[]

[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
  []
  [top_x]
    type = PresetDisplacement
    boundary = 1000
    variable = disp_x
    beta = 0.25
    velocity = vel_x
    acceleration = accel_x
    function = loading_bc
  []
[]

[Functions]
  [loading_bc]
    type = PiecewiseLinear
    data_file = 'input_motion.csv'
    format = columns
  []
[]

[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
    #  store_stress_old = true
    block = 2
  []
  [den_block]
    type = GenericConstantMaterial
    block = 2
    prop_names = density
    prop_values = 0.0002248 #slug/in^3
  []
  [I_Soil]
    [soil_1]
      soil_type = 'darendeli'
      layer_variable = layer_id
      layer_ids = '0'
      over_consolidation_ratio = '1'
      plasticity_index = '0'
      initial_shear_modulus = '7251'
      number_of_points = 10
      poissons_ratio = '0.3'
      block = 1
      density = '0.0001536'
      p_ref = '5'
      pressure_dependency = false
    []
  []
[]

[Contact]
  [leftright]
    secondary = 102
    primary = 103
    model = coulomb
    formulation = penalty
    normalize_penalty = true
    friction_coefficient = 0.2
    penalty = 1e4
    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'
  start_time = 0
  end_time = 2.5
  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 = finalresult
  print_linear_residuals = false
  [screen]
    type = Console
    max_rows = 1
  []
[]

(examples/ex06b/actual_simulation_isoil.i)

[Mesh]
  type = FileMesh
  file = ../ex06b/stabilize_isoil_out1_cp/LATEST #end file of previous simulation
  patch_update_strategy = iteration
  patch_size = 40
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Problem]
  restart_file_base = ../ex06b/stabilize_isoil_out1_cp/LATEST
[]

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

[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
  []
  [layer_id]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[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
  []
  [layer_1]
    type = UniformLayerAuxKernel
    variable = layer_id
    interfaces = '48'
    direction = '0.0 0.0 1.0'
    block = 1
    execute_on = initial
  []
[]

[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
  []
  [top_x]
    type = PresetDisplacement
    boundary = 1000
    variable = disp_x
    beta = 0.25
    velocity = vel_x
    acceleration = accel_x
    function = loading_bc
  []
[]

[Functions]
  [loading_bc]
    type = PiecewiseLinear
    data_file = 'input_motion.csv'
    format = columns
  []
[]

[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
    #  store_stress_old = true
    block = 2
  []
  [den_block]
    type = GenericConstantMaterial
    block = 2
    prop_names = density
    prop_values = 0.0002248 #slug/in^3
  []
  [I_Soil]
    [soil_1]
      soil_type = 'darendeli'
      layer_variable = layer_id
      layer_ids = '0'
      over_consolidation_ratio = '1'
      plasticity_index = '0'
      initial_shear_modulus = '7251'
      number_of_points = 10
      poissons_ratio = '0.3'
      block = 1
      density = '0.0001536'
      p_ref = '5'
      pressure_dependency = false
    []
  []
[]

[Contact]
  [leftright]
    secondary = 102
    primary = 103
    model = coulomb
    formulation = penalty
    normalize_penalty = true
    friction_coefficient = 0.2
    penalty = 1e4
    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'
  start_time = 0
  end_time = 2.5
  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 = finalresult
  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

Example 6b: Displacement-controlled frictional contact problem with I-soil and using restart option for static initialization.

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