Example 6a: Displacement-controlled contact problem with elastic soil

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

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  type = FileMesh
  file = contact_modified.e
  patch_update_strategy = iteration #For contact
  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"}>>>]
  [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
  []
  [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
  []
  [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"}>>>
    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 commented before I got it
  []
[]

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

[Controls<<<{"href": "../syntax/Controls/index.html"}>>>]
  [inertia_switch]
    type = TimePeriod<<<{"description": "Control the enabled/disabled state of objects with time.", "href": "../source/controls/TimePeriod.html"}>>>
    start_time<<<{"description": "The time at which the objects are to be enabled/disabled."}>>> = 0.9
    end_time<<<{"description": "The time at which the objects are to be enable/disabled."}>>> = 1
    disable_objects<<<{"description": "A list of object tags to disable."}>>> = '*/inertia_x */inertia_y */inertia_z */vel_x */vel_y */vel_z */accel_x */accel_y */accel_z'
    set_sync_times<<<{"description": "Set the start and end time as execute sync times."}>>> = true
  []
[]

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

The mesh used in this model was generated externally using the mesh generation software, Cubit. MASTODON imports the mesh shown in Figure 1 using the FileMesh type specified in the Mesh block. Details regarding supported mesh file types can be found at GettingStarted.

Figure 1: Input model in MASTODON.

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.

TensorMechanics and InertialForce kernels are used to model the dynamics of the 3D solid elements. The Gravity block is specified in Kernels in order to apply gravitational force to the entire model. The Newmark time integration parameters used in this problem correspond to the Newmark's average acceleration method, i.e. beta = 0.25 and gamma = 0.5 and damping is ignored.

The ComputeIsotropicElasticityTensorBeam block is used to create the elasticity tensor of the solid block and the elastic soil, 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.

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 and it is given the prescribed displacement history shown in Figure 2.

Figure 2: Input displacement of the 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.7. 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.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'
  []
[]

Static analysis is performed by turning off inertial terms for a certain time period in the Controls block. It stabilizes the system under the gravitational loading and pressure to perform dynamic analysis. The Executioner and the Postprocessors blocks are provided to solve the nonlinear system and extract the required data.

Results

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

Figure 3: Graph showing the frictional force as a function of displacement for Columb Friction Model

The graph for the normal force as a function of time obtained from MASTODON is shown in Figure 4.

Figure 4: Graph showing the normal force as a function of time for Columb 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 = 139.123 lbf

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

(examples/ex06a/psf_grav.i)

[Mesh]
  type = FileMesh
  file = contact_modified.e
  patch_update_strategy = iteration #For contact
  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]
  [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
  []
  [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
  []
  [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
    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 commented before I got it
  []
[]

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

[Controls]
  [inertia_switch]
    type = TimePeriod
    start_time = 0.9
    end_time = 1
    disable_objects = '*/inertia_x */inertia_y */inertia_z */vel_x */vel_y */vel_z */accel_x */accel_y */accel_z'
    set_sync_times = true
  []
[]

[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.9
  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
  perf_graph = true
  print_linear_residuals = false
  [screen]
    type = Console
    max_rows = 1
  []
[]

(examples/ex06a/psf_grav.i)

[Mesh]
  type = FileMesh
  file = contact_modified.e
  patch_update_strategy = iteration #For contact
  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]
  [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
  []
  [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
  []
  [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
    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 commented before I got it
  []
[]

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

[Controls]
  [inertia_switch]
    type = TimePeriod
    start_time = 0.9
    end_time = 1
    disable_objects = '*/inertia_x */inertia_y */inertia_z */vel_x */vel_y */vel_z */accel_x */accel_y */accel_z'
    set_sync_times = true
  []
[]

[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.9
  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
  perf_graph = true
  print_linear_residuals = false
  [screen]
    type = Console
    max_rows = 1
  []
[]

Model Description
Modeling in MASTODON
Results
Theoretical Solution

Usage

Development

Demonstration

Example 6a: Displacement-controlled contact problem with elastic soil

Model Description

Modeling in MASTODON

Results

Theoretical Solution