Two-dimensional spherical indenter (mortar)

A two-dimensional problem with RZ symmetry is used to model the penetration of a spherical indenter into an inelastic base material.

Background

Indentation tests are often used to characterize the behavior of materials at small scales. In this example, we use a spherical indenter driven by a prescribed displacement as a boundary condition. Frictionless contact with a lower-dimensional enforcement (mortar) formulation is employed to drive base material deformation. As a result, a load displacement curve can be obtained.

Creating contact input

Mechanical contact can be enforced on lower-dimensional domains in a weak sense. This type of approach is usually referred to as mortar. To employ this approach, the user can manually build the lower-dimensional subdomains. primary and secondary subdomains are created from mesh sidesets.

[Mesh]
  patch_update_strategy = auto
  patch_size = 2
  partitioner = centroid
  centroid_partitioner_direction = y

  [simple_mesh]
    type = FileMeshGenerator
    file = indenter_rz_fine_bigsideset.e
  []
  # For NodalVariableValue to work with distributed mesh
  allow_renumbering = false
[]

[Functions]
  [disp_y]
    type = PiecewiseLinear
    x = '0.  1.0     2.0    2.6   3.0'
    y = '0.  -4.5   -5.7   -5.7  -4.0'
  []
[]

[Variables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
    block = '1 2'
  []

  [disp_y]
    order = FIRST
    family = LAGRANGE
    block = '1 2'
  []
[]

[AuxVariables]
  [saved_x]
  []
  [saved_y]
  []
[]

[Modules/TensorMechanics/Master]
  [all]
    add_variables = true
    strain = FINITE
    block = '1 2'
    use_automatic_differentiation = false
    generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_zz'
    save_in = 'saved_x saved_y'
  []
[]

[BCs]
  # Symmetries of the Problem
  [symm_x_indenter]
    type = DirichletBC
    variable = disp_x
    boundary = 5
    value = 0.0
  []

  [symm_x_material]
    type = DirichletBC
    variable = disp_x
    boundary = 9
    value = 0.0
  []

  # Material should not fly away
  [material_base_y]
    type = DirichletBC
    variable = disp_y
    boundary = 8
    value = 0.0
  []

  # Drive indenter motion
  [disp_y]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = 1
    function = disp_y
  []
[]

[Contact]
  [contact]
    secondary = 4
    primary = 6
    model = frictionless
    # Investigate von Mises stress at the edge
    correct_edge_dropping = true
    formulation = mortar
    c_normal = 1e+2
  []
[]

[UserObjects]
  [slip_rate_gss]
    type = CrystalPlasticitySlipRateGSS
    variable_size = 48
    slip_sys_file_name = input_slip_sys_bcc48.txt
    num_slip_sys_flowrate_props = 2
    flowprops = '1 48 0.0001 0.01'
    uo_state_var_name = state_var_gss
    slip_incr_tol = 10.0
    block = 2
  []
  [slip_resistance_gss]
    type = CrystalPlasticitySlipResistanceGSS
    variable_size = 48
    uo_state_var_name = state_var_gss
    block = 2
  []
  [state_var_gss]
    type = CrystalPlasticityStateVariable
    variable_size = 48
    groups = '0 24 48'
    group_values = '900 1000' #120
    uo_state_var_evol_rate_comp_name = state_var_evol_rate_comp_gss
    scale_factor = 1.0
    block = 2
  []
  [state_var_evol_rate_comp_gss]
    type = CrystalPlasticityStateVarRateComponentGSS
    variable_size = 48
    hprops = '1.4 1000 1200 2.5'
    uo_slip_rate_name = slip_rate_gss
    uo_state_var_name = state_var_gss
    block = 2
  []
[]

[Materials]
  [tensor]
    type = ComputeIsotropicElasticityTensor
    block = '1'
    youngs_modulus = 1.0e7
    poissons_ratio = 0.25
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
    block = '1'
  []

  [crysp]
    type = FiniteStrainUObasedCP
    block = 2
    stol = 1e-2
    tan_mod_type = exact
    uo_slip_rates = 'slip_rate_gss'
    uo_slip_resistances = 'slip_resistance_gss'
    uo_state_vars = 'state_var_gss'
    uo_state_var_evol_rate_comps = 'state_var_evol_rate_comp_gss'
    maximum_substep_iteration = 20
  []
  [elasticity_tensor]
    type = ComputeElasticityTensorCP
    block = 2
    C_ijkl = '265190 113650 113650 265190 113650 265190 75769 75769 75760'
    fill_method = symmetric9
  []
[]

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

[Executioner]
  type = Transient
  solve_type = 'PJFNK'
  petsc_options = '-snes_ksp_ew'

  petsc_options_iname = '-pc_type -snes_linesearch_type -pc_factor_shift_type '
                        '-pc_factor_shift_amount'
  petsc_options_value = 'lu       basic                 NONZERO               1e-15'
  line_search = 'none'
  automatic_scaling = true
  nl_abs_tol = 2.0e-07
  nl_rel_tol = 2.0e-07
  l_max_its = 40
  l_abs_tol = 1e-08
  l_tol = 1e-08
  start_time = 0.0
  dt = 0.01
  end_time = 3.0 # Executioner
[]

[Postprocessors]
  [maxdisp]
    type = NodalVariableValue
    nodeid = 39
    variable = disp_y
  []
  [resid_y]
    type = NodalSum
    variable = saved_y
    boundary = 1
  []
[]

[Outputs]
  [out]
    type = Exodus
    elemental_as_nodal = true
  []
  perf_graph = true
  csv = true
[]
(modules/contact/examples/2d_indenter/indenter_rz_fine.i)
commentnote

Mortar-based mechanical contact can be defined through the contact action. Here, a more manual, user-driven definition is used.

Fig. 1: Spherical indenter.

For frictionless contact in two dimensions, three blocks need to be defined. First, the NormalNodalLMMechanicalContact constraint is used to enforce the Karush-Kuhn-Tucker contact conditions. Then, NormalMortarMechanicalContact enforces contact constaints in an integral or weak sense in both problem dimensions.

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

[Problem]
  coord_type = RZ
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
[]

[Mesh]
  patch_update_strategy = auto
  patch_size = 2
  partitioner = centroid
  centroid_partitioner_direction = y

  [simple_mesh]
    type = FileMeshGenerator
    file = indenter_rz_fine_bigsideset.e
  []
  # For NodalVariableValue to work with distributed mesh
  allow_renumbering = false
[]

[Functions]
  [disp_y]
    type = PiecewiseLinear
    x = '0.  1.0     2.0    2.6   3.0'
    y = '0.  -4.5   -5.7   -5.7  -4.0'
  []
[]

[Variables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
    block = '1 2'
  []

  [disp_y]
    order = FIRST
    family = LAGRANGE
    block = '1 2'
  []
[]

[AuxVariables]
  [saved_x]
  []
  [saved_y]
  []
[]

[Modules/TensorMechanics/Master]
  [all]
    add_variables = true
    strain = FINITE
    block = '1 2'
    use_automatic_differentiation = false
    generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_zz'
    save_in = 'saved_x saved_y'
  []
[]

[BCs]
  # Symmetries of the Problem
  [symm_x_indenter]
    type = DirichletBC
    variable = disp_x
    boundary = 5
    value = 0.0
  []

  [symm_x_material]
    type = DirichletBC
    variable = disp_x
    boundary = 9
    value = 0.0
  []

  # Material should not fly away
  [material_base_y]
    type = DirichletBC
    variable = disp_y
    boundary = 8
    value = 0.0
  []

  # Drive indenter motion
  [disp_y]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = 1
    function = disp_y
  []
[]

[Contact]
  [contact]
    secondary = 4
    primary = 6
    model = frictionless
    # Investigate von Mises stress at the edge
    correct_edge_dropping = true
    formulation = mortar
    c_normal = 1e+2
  []
[]

[UserObjects]
  [slip_rate_gss]
    type = CrystalPlasticitySlipRateGSS
    variable_size = 48
    slip_sys_file_name = input_slip_sys_bcc48.txt
    num_slip_sys_flowrate_props = 2
    flowprops = '1 48 0.0001 0.01'
    uo_state_var_name = state_var_gss
    slip_incr_tol = 10.0
    block = 2
  []
  [slip_resistance_gss]
    type = CrystalPlasticitySlipResistanceGSS
    variable_size = 48
    uo_state_var_name = state_var_gss
    block = 2
  []
  [state_var_gss]
    type = CrystalPlasticityStateVariable
    variable_size = 48
    groups = '0 24 48'
    group_values = '900 1000' #120
    uo_state_var_evol_rate_comp_name = state_var_evol_rate_comp_gss
    scale_factor = 1.0
    block = 2
  []
  [state_var_evol_rate_comp_gss]
    type = CrystalPlasticityStateVarRateComponentGSS
    variable_size = 48
    hprops = '1.4 1000 1200 2.5'
    uo_slip_rate_name = slip_rate_gss
    uo_state_var_name = state_var_gss
    block = 2
  []
[]

[Materials]
  [tensor]
    type = ComputeIsotropicElasticityTensor
    block = '1'
    youngs_modulus = 1.0e7
    poissons_ratio = 0.25
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
    block = '1'
  []

  [crysp]
    type = FiniteStrainUObasedCP
    block = 2
    stol = 1e-2
    tan_mod_type = exact
    uo_slip_rates = 'slip_rate_gss'
    uo_slip_resistances = 'slip_resistance_gss'
    uo_state_vars = 'state_var_gss'
    uo_state_var_evol_rate_comps = 'state_var_evol_rate_comp_gss'
    maximum_substep_iteration = 20
  []
  [elasticity_tensor]
    type = ComputeElasticityTensorCP
    block = 2
    C_ijkl = '265190 113650 113650 265190 113650 265190 75769 75769 75760'
    fill_method = symmetric9
  []
[]

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

[Executioner]
  type = Transient
  solve_type = 'PJFNK'
  petsc_options = '-snes_ksp_ew'

  petsc_options_iname = '-pc_type -snes_linesearch_type -pc_factor_shift_type '
                        '-pc_factor_shift_amount'
  petsc_options_value = 'lu       basic                 NONZERO               1e-15'
  line_search = 'none'
  automatic_scaling = true
  nl_abs_tol = 2.0e-07
  nl_rel_tol = 2.0e-07
  l_max_its = 40
  l_abs_tol = 1e-08
  l_tol = 1e-08
  start_time = 0.0
  dt = 0.01
  end_time = 3.0 # Executioner
[]

[Postprocessors]
  [maxdisp]
    type = NodalVariableValue
    nodeid = 39
    variable = disp_y
  []
  [resid_y]
    type = NodalSum
    variable = saved_y
    boundary = 1
  []
[]

[Outputs]
  [out]
    type = Exodus
    elemental_as_nodal = true
  []
  perf_graph = true
  csv = true
[]
(modules/contact/examples/2d_indenter/indenter_rz_fine.i)

Note that the subdomain blocks had been created in the mesh input using LowerDBlockFromSidesetGenerator.

commentnote

Mortar enforcement is only available for two-dimensional contact.

Other input

The problem is axisymmetric Compute Axisymmetric RZ Finite Strain and symmetric boundary conditions are used.

Numerical results

The resulting force exerted as material resistance on the indenter may be plotted against the vertical displacement. In this problem, the base material is a monocrystal with body-centered cubic (bcc) unit cell with arbitrary parameters. Platic deformation causes the piling up of the base material's contact surface, as shown in the animation in Fig. 1.

Fig. 2: Load-displacement curve.

Crystal plasticity parameters can be calibrated to match a given experimental nano-indentation test. For this example, the load-displacement curve is shown in Fig. 2.

Notes:

  • Friction may alter results

  • Indenter geometry does not reproduce that of a real problem's. Its geometry can be modified in the journal file.

  • Element distortions may become large