ComputeDynamicWeightedGapLMMechanicalContact

Computes the normal contact mortar constraints for dynamic simulations

This object is virtually analogous to ComputeWeightedGapLMMechanicalContact, but it is used for dynamic simulations.

When the mortar mechanical contact constraints are used in dynamic simulations, the normal contact constraints need to be stabilized by ensuring that the normal gap time derivative is included to guarantee a dynamic contact. This is means of guaranteeing the 'persistency' condition, i.e. not only do we enforce the constraints instantaneously, but also that it will remain in contact. This approximate contact constraint stabilization is performed in ComputeDynamicWeightedGapLMMechanicalContact. Once nodal contact is established, the constraint enforcement is switched to the persistency equation:

$var element = document.getElementById("moose-equation-4e7d3d0a-0116-440d-9e03-c750bb6dcbb9");katex.render("\\begin{aligned} (\\tilde{g}_n)_j \\coloneqq \\dot{\\tilde{g}}_{nj} \\cdot \\Delta t \\end{aligned}", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

where $var element = document.getElementById("moose-equation-1623dd53-87b0-4504-b4ba-dca9aacd46ef");katex.render("\\,\\dot{}\\,", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ denotes the first time derivative, $var element = document.getElementById("moose-equation-02974bbd-8601-49ad-adb6-64049b206944");katex.render("\\Delta t", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the time step increment, $var element = document.getElementById("moose-equation-8b214aaf-cb48-4608-a86b-d47d74dde9c5");katex.render("j", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ refers to an arbitrary node.

The capture_tolerance is an optional contact parameter used in dynamic contact constraints to determine when to impose the persistency condition for normal contact. For relevant, general equations, see (Tal and Hager, 2018).

Input Parameters

disp_xThe x displacement variable
C++ Type:std::vector<VariableName>
Controllable:No
Description:The x displacement variable
disp_yThe y displacement variable
C++ Type:std::vector<VariableName>
Controllable:No
Description:The y displacement variable
newmark_betaBeta parameter for the Newmark time integrator
C++ Type:double
Controllable:No
Description:Beta parameter for the Newmark time integrator
newmark_gammaGamma parameter for the Newmark time integrator
C++ Type:double
Controllable:No
Description:Gamma parameter for the Newmark time integrator
primary_boundaryThe name of the primary boundary sideset.
C++ Type:BoundaryName
Controllable:No
Description:The name of the primary boundary sideset.
primary_subdomainThe name of the primary subdomain.
C++ Type:SubdomainName
Controllable:No
Description:The name of the primary subdomain.
secondary_boundaryThe name of the secondary boundary sideset.
C++ Type:BoundaryName
Controllable:No
Description:The name of the secondary boundary sideset.
secondary_subdomainThe name of the secondary subdomain.
C++ Type:SubdomainName
Controllable:No
Description:The name of the secondary subdomain.

Required Parameters

c1e+06Parameter for balancing the size of the gap and contact pressure
Default:1e+06
C++ Type:double
Controllable:No
Description:Parameter for balancing the size of the gap and contact pressure
capture_tolerance1e-05Parameter describing a gap threshold for the application of the persistency constraint in dynamic simulations.
Default:1e-05
C++ Type:double
Controllable:No
Description:Parameter describing a gap threshold for the application of the persistency constraint in dynamic simulations.
compute_lm_residualsTrueWhether to compute Lagrange Multiplier residuals
Default:True
C++ Type:bool
Controllable:No
Description:Whether to compute Lagrange Multiplier residuals
compute_primal_residualsTrueWhether to compute residuals for the primal variable.
Default:True
C++ Type:bool
Controllable:No
Description:Whether to compute residuals for the primal variable.
correct_edge_droppingFalseWhether to enable correct edge dropping treatment for mortar constraints. When disabled any Lagrange Multiplier degree of freedom on a secondary element without full primary contributions will be set (strongly) to 0.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to enable correct edge dropping treatment for mortar constraints. When disabled any Lagrange Multiplier degree of freedom on a secondary element without full primary contributions will be set (strongly) to 0.
debug_meshFalseWhether this constraint is going to enable mortar segment mesh debug information. An exodusfile will be generated if the user sets this flag to true
Default:False
C++ Type:bool
Controllable:No
Description:Whether this constraint is going to enable mortar segment mesh debug information. An exodusfile will be generated if the user sets this flag to true
disp_zThe z displacement variable
C++ Type:std::vector<VariableName>
Controllable:No
Description:The z displacement variable
ghost_point_neighborsFalseWhether we should ghost point neighbors of secondary face elements, and consequently also their mortar interface couples.
Default:False
C++ Type:bool
Controllable:No
Description:Whether we should ghost point neighbors of secondary face elements, and consequently also their mortar interface couples.
interpolate_normalsTrueWhether to interpolate the nodal normals (e.g. classic idea of evaluating field at quadrature points). If this is set to false, then non-interpolated nodal normals will be used, and then the _normals member should be indexed with _i instead of _qp
Default:True
C++ Type:bool
Controllable:No
Description:Whether to interpolate the nodal normals (e.g. classic idea of evaluating field at quadrature points). If this is set to false, then non-interpolated nodal normals will be used, and then the _normals member should be indexed with _i instead of _qp
normalize_cFalseWhether to normalize c by weighting function norm. When unnormalized the value of c effectively depends on element size since in the constraint we compare nodal Lagrange Multiplier values to integrated gap values (LM nodal value is independent of element size, where integrated values are dependent on element size).
Default:False
C++ Type:bool
Controllable:No
Description:Whether to normalize c by weighting function norm. When unnormalized the value of c effectively depends on element size since in the constraint we compare nodal Lagrange Multiplier values to integrated gap values (LM nodal value is independent of element size, where integrated values are dependent on element size).
periodicFalseWhether this constraint is going to be used to enforce a periodic condition. This has the effect of changing the normals vector for projection from outward to inward facing
Default:False
C++ Type:bool
Controllable:No
Description:Whether this constraint is going to be used to enforce a periodic condition. This has the effect of changing the normals vector for projection from outward to inward facing
prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
quadratureDEFAULTQuadrature rule to use on mortar segments. For 2D mortar DEFAULT is recommended. For 3D mortar, QUAD meshes are integrated using triangle mortar segments. While DEFAULT quadrature order is typically sufficiently accurate, exact integration of QUAD mortar faces requires SECOND order quadrature for FIRST variables and FOURTH order quadrature for SECOND order variables.
Default:DEFAULT
C++ Type:MooseEnum
Options:DEFAULT, FIRST, SECOND, THIRD, FOURTH
Controllable:No
Description:Quadrature rule to use on mortar segments. For 2D mortar DEFAULT is recommended. For 3D mortar, QUAD meshes are integrated using triangle mortar segments. While DEFAULT quadrature order is typically sufficiently accurate, exact integration of QUAD mortar faces requires SECOND order quadrature for FIRST variables and FOURTH order quadrature for SECOND order variables.
variableThe name of the lagrange multiplier variable that this constraint is applied to. This parameter may not be supplied in the case of using penalty methods for example
C++ Type:NonlinearVariableName
Controllable:No
Description:The name of the lagrange multiplier variable that this constraint is applied to. This parameter may not be supplied in the case of using penalty methods for example
wear_depthThe name of the mortar auxiliary variable that is used to modify the weighted gap definition
C++ Type:std::vector<VariableName>
Controllable:No
Description:The name of the mortar auxiliary variable that is used to modify the weighted gap definition

Optional Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
use_displaced_meshTrueWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:True
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Options:nontime, time
Controllable:No
Description:The tag for the vectors this Kernel should fill

Tagging Parameters

Input Files

(modules/contact/test/tests/mortar_aux_kernels/block-dynamics-aux-wear.i)
(modules/contact/test/tests/mortar_aux_kernels/block-dynamics-aux-vel.i)
(modules/contact/test/tests/mortar_dynamics/block-dynamics.i)

Child Objects

(modules/contact/include/constraints/ComputeDynamicFrictionalForceLMMechanicalContact.h)

References

Yuval Tal and Bradford H Hager. Dynamic mortar finite element method for modeling of shear rupture on frictional rough surfaces. Computational Mechanics, 61(6):699–716, 2018.

@article{tal2018dynamic,
    author = "Tal, Yuval and Hager, Bradford H",
    title = "Dynamic mortar finite element method for modeling of shear rupture on frictional rough surfaces",
    journal = "Computational Mechanics",
    volume = "61",
    number = "6",
    pages = "699--716",
    year = "2018",
    publisher = "Springer"
}

(modules/contact/test/tests/mortar_aux_kernels/block-dynamics-aux-wear.i)

starting_point = 0.5e-1
offset = -0.05

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  file = long-bottom-block-1elem-blocks.e
[]

[Variables]
  [disp_x]
    block = '1 2'
  []
  [disp_y]
    block = '1 2'
  []
  [normal_lm]
    block = 3
    use_dual = true
  []
[]

[ICs]
  [disp_y]
    block = 2
    variable = disp_y
    value = '${fparse starting_point + offset}'
    type = ConstantIC
  []
[]

[Kernels]
  [DynamicTensorMechanics]
    displacements = 'disp_x disp_y'
    generate_output = 'stress_xx stress_yy'
    strain = FINITE
    block = '1 2'
    zeta = 1.0
    hht_alpha = 0.0
  []
  [inertia_x]
    type = InertialForce
    variable = disp_x
    velocity = vel_x
    acceleration = accel_x
    beta = 0.25
    gamma = 0.5
    alpha = 0
    eta = 0.0
    block = '1 2'
  []
  [inertia_y]
    type = InertialForce
    variable = disp_y
    velocity = vel_y
    acceleration = accel_y
    beta = 0.25
    gamma = 0.5
    alpha = 0
    eta = 0.0
    block = '1 2'
  []
[]

[Materials]
  [elasticity_2]
    type = ComputeIsotropicElasticityTensor
    block = '2'
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  []
  [elasticity_1]
    type = ComputeIsotropicElasticityTensor
    block = '1'
    youngs_modulus = 1e8
    poissons_ratio = 0.3
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
    block = '1 2'
  []
  [strain]
    type = ComputeFiniteStrain
    block = '1 2'
  []
  [density]
    type = GenericConstantMaterial
    block = '1 2'
    prop_names = 'density'
    prop_values = '7750'
  []
[]

[AuxVariables]
  [worn_depth]
    block = '3'
  []
  [vel_x]
    block = '1 2'
  []
  [accel_x]
    block = '1 2'
  []
  [vel_y]
    block = '1 2'
  []
  [accel_y]
    block = '1 2'
  []
  [vel_z]
    block = '1 2'
  []
  [accel_z]
    block = '1 2'
  []
[]

[AuxKernels]
  [worn_depth]
    type = MortarArchardsLawAux
    variable = worn_depth
    primary_boundary = 20
    secondary_boundary = 10
    primary_subdomain = 4
    secondary_subdomain = 3
    displacements = 'disp_x disp_y'
    friction_coefficient = 0.5
    energy_wear_coefficient = 1.0
    normal_pressure = normal_lm
  []
  [accel_x]
    type = NewmarkAccelAux
    variable = accel_x
    displacement = disp_x
    velocity = vel_x
    beta = 0.25
    execute_on = 'linear timestep_end'
  []
  [vel_x]
    type = NewmarkVelAux
    variable = vel_x
    acceleration = accel_x
    gamma = 0.5
    execute_on = 'linear timestep_end'
  []
  [accel_y]
    type = NewmarkAccelAux
    variable = accel_y
    displacement = disp_y
    velocity = vel_y
    beta = 0.25
    execute_on = 'linear timestep_end'
  []
  [vel_y]
    type = NewmarkVelAux
    variable = vel_y
    acceleration = accel_y
    gamma = 0.5
    execute_on = 'linear timestep_end'
  []
[]

[Constraints]
  [weighted_gap_lm]
    type = ComputeDynamicWeightedGapLMMechanicalContact
    primary_boundary = 20
    secondary_boundary = 10
    primary_subdomain = 4
    secondary_subdomain = 3
    variable = normal_lm
    disp_x = disp_x
    disp_y = disp_y
    use_displaced_mesh = true
    c = 1e4
    interpolate_normals = false
    capture_tolerance = 1.0e-5
    newmark_beta = 0.25
    newmark_gamma = 0.5
  []
  [normal_x]
    type = NormalMortarMechanicalContact
    primary_boundary = 20
    secondary_boundary = 10
    primary_subdomain = 4
    secondary_subdomain = 3
    variable = normal_lm
    secondary_variable = disp_x
    component = x
    use_displaced_mesh = true
    compute_lm_residuals = false
    interpolate_normals = false
  []
  [normal_y]
    type = NormalMortarMechanicalContact
    primary_boundary = 20
    secondary_boundary = 10
    primary_subdomain = 4
    secondary_subdomain = 3
    variable = normal_lm
    secondary_variable = disp_y
    component = y
    use_displaced_mesh = true
    compute_lm_residuals = false
    interpolate_normals = false
  []
[]

[BCs]
  [botx]
    type = DirichletBC
    variable = disp_x
    boundary = 40
    value = 0.0
  []
  [boty]
    type = DirichletBC
    variable = disp_y
    boundary = 40
    value = 0.0
  []
  [topy]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = 30
    function = '${starting_point} * cos(8.0 * pi / 4 * t) + ${offset}'
  []
  [leftx]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = 50
    function = '1e-2 * t'
  []
[]

[Executioner]
  type = Transient
  end_time = 0.675
  dt = 0.075
  dtmin = .075
  solve_type = 'NEWTON'
  petsc_options = '-snes_converged_reason -ksp_converged_reason -pc_svd_monitor '
                  '-snes_linesearch_monitor -snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -pc_factor_shift_type -pc_factor_shift_amount '
  petsc_options_value = 'lu       superlu_dist                  NONZERO               1e-15'
  nl_max_its = 30
  line_search = 'l2'
  snesmf_reuse_base = false

  [TimeIntegrator]
    type = NewmarkBeta
    beta = 0.25
    gamma = 0.5
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[Outputs]
  exodus = true
  checkpoint = true
[]

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

[Postprocessors]
  active = 'num_nl cumulative contact'
  [num_nl]
    type = NumNonlinearIterations
  []
  [cumulative]
    type = CumulativeValuePostprocessor
    postprocessor = num_nl
  []
  [contact]
    type = ContactDOFSetSize
    variable = normal_lm
    subdomain = '3'
    execute_on = 'nonlinear timestep_end'
  []
[]

(modules/contact/test/tests/mortar_aux_kernels/block-dynamics-aux-vel.i)

starting_point = 2e-1
offset = -0.19

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  file = long-bottom-block-1elem-blocks.e
[]

[Variables]
  [disp_x]
    block = '1 2'
  []
  [disp_y]
    block = '1 2'
  []
  [normal_lm]
    block = 3
    use_dual = true
  []
[]

[ICs]
  [disp_y]
    block = 2
    variable = disp_y
    value = '${fparse starting_point + offset}'
    type = ConstantIC
  []
[]

[Kernels]
  [DynamicTensorMechanics]
    displacements = 'disp_x disp_y'
    generate_output = 'stress_xx stress_yy'
    strain = FINITE
    block = '1 2'
    zeta = 1.0
    alpha = 0.0
  []
  [inertia_x]
    type = InertialForce
    variable = disp_x
    velocity = vel_x
    acceleration = accel_x
    beta = 0.25
    gamma = 0.5
    alpha = 0
    eta = 0.0
    block = '1 2'
  []
  [inertia_y]
    type = InertialForce
    variable = disp_y
    velocity = vel_y
    acceleration = accel_y
    beta = 0.25
    gamma = 0.5
    alpha = 0
    eta = 0.0
    block = '1 2'
  []
[]

[Materials]
  [elasticity_2]
    type = ComputeIsotropicElasticityTensor
    block = '2'
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  []
  [elasticity_1]
    type = ComputeIsotropicElasticityTensor
    block = '1'
    youngs_modulus = 1e8
    poissons_ratio = 0.3
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
    block = '1 2'
  []
  [strain]
    type = ComputeFiniteStrain
    block = '1 2'
  []
  [density]
    type = GenericConstantMaterial
    block = '1 2'
    prop_names = 'density'
    prop_values = '7750'
  []
[]

[AuxVariables]
  [gap_vel]
    block = '3'
  []
  [vel_x]
    block = '1 2'
  []
  [accel_x]
    block = '1 2'
  []
  [vel_y]
    block = '1 2'
  []
  [accel_y]
    block = '1 2'
  []
  [vel_z]
    block = '1 2'
  []
  [accel_z]
    block = '1 2'
  []
[]

[AuxKernels]
  [gap_vel]
    type = WeightedGapVelAux
    variable = gap_vel
    primary_boundary = 20
    secondary_boundary = 10
    primary_subdomain = 4
    secondary_subdomain = 3
    disp_x = disp_x
    disp_y = disp_y
  []
  [accel_x]
    type = NewmarkAccelAux
    variable = accel_x
    displacement = disp_x
    velocity = vel_x
    beta = 0.25
    execute_on = 'linear timestep_end'
  []
  [vel_x]
    type = NewmarkVelAux
    variable = vel_x
    acceleration = accel_x
    gamma = 0.5
    execute_on = 'linear timestep_end'
  []
  [accel_y]
    type = NewmarkAccelAux
    variable = accel_y
    displacement = disp_y
    velocity = vel_y
    beta = 0.25
    execute_on = 'linear timestep_end'
  []
  [vel_y]
    type = NewmarkVelAux
    variable = vel_y
    acceleration = accel_y
    gamma = 0.5
    execute_on = 'linear timestep_end'
  []
[]

[Constraints]
  [weighted_gap_lm]
    type = ComputeDynamicWeightedGapLMMechanicalContact
    primary_boundary = 20
    secondary_boundary = 10
    primary_subdomain = 4
    secondary_subdomain = 3
    variable = normal_lm
    disp_x = disp_x
    disp_y = disp_y
    use_displaced_mesh = true
    c = 1e4
    interpolate_normals = false
    capture_tolerance = 1.0e-5
    newmark_beta = 0.25
    newmark_gamma = 0.5
  []
  [normal_x]
    type = NormalMortarMechanicalContact
    primary_boundary = 20
    secondary_boundary = 10
    primary_subdomain = 4
    secondary_subdomain = 3
    variable = normal_lm
    secondary_variable = disp_x
    component = x
    use_displaced_mesh = true
    compute_lm_residuals = false
    interpolate_normals = false
  []
  [normal_y]
    type = NormalMortarMechanicalContact
    primary_boundary = 20
    secondary_boundary = 10
    primary_subdomain = 4
    secondary_subdomain = 3
    variable = normal_lm
    secondary_variable = disp_y
    component = y
    use_displaced_mesh = true
    compute_lm_residuals = false
    interpolate_normals = false
  []
[]

[BCs]
  [botx]
    type = DirichletBC
    variable = disp_x
    boundary = 40
    value = 0.0
  []
  [boty]
    type = DirichletBC
    variable = disp_y
    boundary = 40
    value = 0.0
  []
  [topy]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = 30
    function = '${starting_point} * cos(2 * pi / 4 * t) + ${offset}'
  []
  [leftx]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = 50
    function = '1e-2 * t'
  []
[]

[Executioner]
  type = Transient
  end_time = 0.1
  dt = 0.05
  dtmin = 0.05
  solve_type = 'PJFNK'
  petsc_options = '-snes_converged_reason -ksp_converged_reason -pc_svd_monitor '
                  '-snes_linesearch_monitor -snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -pc_factor_shift_type -pc_factor_shift_amount -mat_mffd_err '
  petsc_options_value = 'lu       superlu_dist                  NONZERO               1e-15                   1e-5'
  nl_max_its = 20
  line_search = 'none'
  snesmf_reuse_base = false

  [TimeIntegrator]
    type = NewmarkBeta
    beta = 0.25
    gamma = 0.5
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[Outputs]
  exodus = true
  checkpoint = true
[]

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

[Postprocessors]
  active = 'num_nl cumulative contact'
  [num_nl]
    type = NumNonlinearIterations
  []
  [cumulative]
    type = CumulativeValuePostprocessor
    postprocessor = num_nl
  []
  [contact]
    type = ContactDOFSetSize
    variable = normal_lm
    subdomain = '3'
    execute_on = 'nonlinear timestep_end'
  []
[]

(modules/contact/test/tests/mortar_dynamics/block-dynamics.i)

starting_point = 2e-1
offset = -0.19

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  file = long-bottom-block-1elem-blocks.e
[]

[Variables]
  [disp_x]
    block = '1 2'
  []
  [disp_y]
    block = '1 2'
  []
  [normal_lm]
    block = 3
    use_dual = true
  []
[]

[ICs]
  [disp_y]
    block = 2
    variable = disp_y
    value = '${fparse starting_point + offset}'
    type = ConstantIC
  []
[]

[Kernels]
  [DynamicTensorMechanics]
    displacements = 'disp_x disp_y'
    generate_output = 'stress_xx stress_yy'
    strain = FINITE
    block = '1 2'
    zeta = 1.0
    alpha = 0.0
  []
  [inertia_x]
    type = InertialForce
    variable = disp_x
    velocity = vel_x
    acceleration = accel_x
    beta = 0.25
    gamma = 0.5
    alpha = 0
    eta = 0.0
    block = '1 2'
  []
  [inertia_y]
    type = InertialForce
    variable = disp_y
    velocity = vel_y
    acceleration = accel_y
    beta = 0.25
    gamma = 0.5
    alpha = 0
    eta = 0.0
    block = '1 2'
  []
[]

[Materials]
  [elasticity_2]
    type = ComputeIsotropicElasticityTensor
    block = '2'
    youngs_modulus = 1e6
    poissons_ratio = 0.3
  []
  [elasticity_1]
    type = ComputeIsotropicElasticityTensor
    block = '1'
    youngs_modulus = 1e8
    poissons_ratio = 0.3
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
    block = '1 2'
  []
  [strain]
    type = ComputeFiniteStrain
    block = '1 2'
  []
  [density]
    type = GenericConstantMaterial
    block = '1 2'
    prop_names = 'density'
    prop_values = '7750'
  []
[]

[AuxVariables]
  [vel_x]
    block = '1 2'
  []
  [accel_x]
    block = '1 2'
  []
  [vel_y]
    block = '1 2'
  []
  [accel_y]
    block = '1 2'
  []
  [vel_z]
    block = '1 2'
  []
  [accel_z]
    block = '1 2'
  []
[]

[AuxKernels]
  [accel_x]
    type = NewmarkAccelAux
    variable = accel_x
    displacement = disp_x
    velocity = vel_x
    beta = 0.25
    execute_on = 'LINEAR timestep_end'
  []
  [vel_x]
    type = NewmarkVelAux
    variable = vel_x
    acceleration = accel_x
    gamma = 0.5
    execute_on = 'LINEAR timestep_end'
  []
  [accel_y]
    type = NewmarkAccelAux
    variable = accel_y
    displacement = disp_y
    velocity = vel_y
    beta = 0.25
    execute_on = 'LINEAR timestep_end'
  []
  [vel_y]
    type = NewmarkVelAux
    variable = vel_y
    acceleration = accel_y
    gamma = 0.5
    execute_on = 'LINEAR timestep_end'
  []
[]

[Constraints]
  [weighted_gap_lm]
    type = ComputeDynamicWeightedGapLMMechanicalContact
    primary_boundary = 20
    secondary_boundary = 10
    primary_subdomain = 4
    secondary_subdomain = 3
    variable = normal_lm
    disp_x = disp_x
    disp_y = disp_y
    use_displaced_mesh = true
    c = 1e4
    interpolate_normals = false
    capture_tolerance = 1.0e-5
    newmark_beta = 0.25
    newmark_gamma = 0.5
  []
  [normal_x]
    type = NormalMortarMechanicalContact
    primary_boundary = 20
    secondary_boundary = 10
    primary_subdomain = 4
    secondary_subdomain = 3
    variable = normal_lm
    secondary_variable = disp_x
    component = x
    use_displaced_mesh = true
    compute_lm_residuals = false
    interpolate_normals = false
  []
  [normal_y]
    type = NormalMortarMechanicalContact
    primary_boundary = 20
    secondary_boundary = 10
    primary_subdomain = 4
    secondary_subdomain = 3
    variable = normal_lm
    secondary_variable = disp_y
    component = y
    use_displaced_mesh = true
    compute_lm_residuals = false
    interpolate_normals = false
  []
[]

[BCs]
  [botx]
    type = DirichletBC
    variable = disp_x
    boundary = 40
    value = 0.0
  []
  [boty]
    type = DirichletBC
    variable = disp_y
    boundary = 40
    value = 0.0
  []
  [topy]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = 30
    function = '${starting_point} * cos(2 * pi / 4 * t) + ${offset}'
  []
  [leftx]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = 50
    function = '1e-2 * t'
  []
[]

[Executioner]
  type = Transient
  end_time = 75
  dt = 0.05
  dtmin = .05
  solve_type = 'PJFNK'
  petsc_options = '-snes_converged_reason -ksp_converged_reason -pc_svd_monitor '
                  '-snes_linesearch_monitor -snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount -mat_mffd_err '
  petsc_options_value = 'lu       NONZERO               1e-15                   1e-5'
  nl_max_its = 20
  line_search = 'none'
  snesmf_reuse_base = false

  [TimeIntegrator]
    type = NewmarkBeta
    beta = 0.25
    gamma = 0.5
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[Outputs]
  exodus = true
  checkpoint = true
[]

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

[Postprocessors]
  active = 'num_nl cumulative contact'
  [num_nl]
    type = NumNonlinearIterations
  []
  [cumulative]
    type = CumulativeValuePostprocessor
    postprocessor = num_nl
  []
  [contact]
    type = ContactDOFSetSize
    variable = normal_lm
    subdomain = '3'
    execute_on = 'nonlinear timestep_end'
  []
[]

(modules/contact/include/constraints/ComputeDynamicFrictionalForceLMMechanicalContact.h)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "ComputeDynamicWeightedGapLMMechanicalContact.h"

/**
 * Computes the mortar tangential frictional forces for dynamic simulations.
 */
class ComputeDynamicFrictionalForceLMMechanicalContact
  : public ComputeDynamicWeightedGapLMMechanicalContact
{
public:
  static InputParameters validParams();

  ComputeDynamicFrictionalForceLMMechanicalContact(const InputParameters & parameters);
  void residualSetup() override;
  void post() override;

  /**
   * Copy of the post routine but that skips assembling inactive nodes
   */
  void
  incorrectEdgeDroppingPost(const std::unordered_set<const Node *> & inactive_lm_nodes) override;

protected:
  /**
   * Computes properties that are functions only of the current quadrature point (\p _qp), e.g.
   * indepedent of shape functions
   */
  virtual void computeQpProperties() override;

  /**
   * Computes properties that are functions both of \p _qp and \p _i, for example the weighted gap
   */
  virtual void computeQpIProperties() override;

  /**
   * Method called from \p post(). Used to enforce node-associated constraints. E.g. for the base \p
   * ComputeFrictionalForceLMMechanicalContact we enforce PDASS frictional constraints. This is also
   * where we actually feed the node-based constraint information into the system residual and
   * Jacobian
   */
  virtual void enforceConstraintOnDof(const DofObject * const dof) override;

  /**
   * Method called from \p post(). Used to enforce node-associated constraints. E.g. for the base \p
   * ComputeFrictionalForceLMMechanicalContact we enforce PDASS frictional constraints. This is also
   * where we actually feed the node-based constraint information into the system residual and
   * Jacobian. Method used to enforce frictional constraints in three-dimensional simulations.
   */
  virtual void enforceConstraintOnDof3d(const DofObject * const dof);

  /**
   * Communicate weighted velocities to the owning process
   */
  void communicateVelocities();

  /**
   * Apply constant or function-based friction coefficient
   */
  ADReal computeFrictionValue(const ADReal & contact_pressure,
                              const ADReal & tangential_vel,
                              const ADReal & tangential_vel_dir);

  /// A map from node to two tangential velocities
  std::unordered_map<const DofObject *, std::array<ADReal, 2>> _dof_to_weighted_tangential_velocity;

  /// An array of two pointers to avoid copies
  std::array<const ADReal *, 2> _tangential_vel_ptr = {{nullptr, nullptr}};

  /// The value of the tangential velocity values at the current quadrature point
  std::array<ADReal, 2> _qp_tangential_velocity;

  /// The value of the tangential velocity vectors at the current node
  ADRealVectorValue _qp_tangential_velocity_nodal;

  /// Numerical factor used in the tangential constraints for convergence purposes
  const Real _c_t;

  /// Frictional Lagrange's multiplier variable pointers
  std::vector<MooseVariable *> _friction_vars;

  /// x-velocity on the secondary face
  const ADVariableValue & _secondary_x_dot;

  /// x-velocity on the primary face
  const ADVariableValue & _primary_x_dot;

  /// y-velocity on the secondary face
  const ADVariableValue & _secondary_y_dot;

  /// y-velocity on the primary face
  const ADVariableValue & _primary_y_dot;

  /// z-velocity on the secondary face
  const ADVariableValue * const _secondary_z_dot;

  /// z-velocity on the primary face
  const ADVariableValue * const _primary_z_dot;

  /// Small contact pressure value to trigger computation of frictional forces
  const Real _epsilon;

  /// Friction coefficient
  const Real _mu;

  /// input function
  const Function * const _function_friction;

  /// Boolean to determine whether the friction coefficient is taken from a function
  const bool _has_friction_function;

  /// Automatic flag to determine whether we are doing three-dimensional work
  bool _3d;
};

Input Parameters
Input Files
Child Objects
References

Install MOOSE

New Users

Examples and Tutorials

Application Usage

Physics and Syntax

Application Development

Framework Development

MOOSEDocs

Infrastructure

Questions

Information and Tools

INL Applications and Remote Access