Discrepancy between Mortar CZM and CZM Action

[mpai60]

I'm working building a polycrystalline model with CZMs for grain boundaries. I am having trouble with the CZM action as the interfaces do not seem to separate when I use this CZM Action input block specifically. Here are the results compared to mortar CZM with a toy geometry to illustrate the issue:

Both cases have a sliding constraint on the right and bottom faces (Dirichlet BC). The mortar CZM (left) provides the expected results. All the CZM parameters are the same for both cases. I cannot figure out why the CZMs on the right are separating differently from the ones on the left.

Input files, mesh, and Euler angle files are here: https://github.com/mpai60/czm_testing
Both inputs use the same mesh and Euler angle files.

[GiudGiud]

hello

a couple things I noticed:

• the solid mechanics and other actions need to be block restricted to the full-dimensional blocks. e.g. you need to set block = '1 2 3' to these physics
• the property file reader is set up differently in both inputs, it gives you a warning for that though results are not affected

the big difference is the material properties. They differ by orders of magnitude between the two inputs. why is that?
I changed this, and got closer results

  [czm]
    type = BiLinearMixedModeTraction
    boundary = czm_boundary # '561 562 564' - This is one sideset from each interface, but using this doesn't seem to affect the results.
    penalty_stiffness = 1e6
    GI_c = 1e5
    GII_c =1e2
    normal_strength = 1e6
    shear_strength = 1e3
    displacements = 'disp_x disp_y'
    eta = 2.2
    viscosity = 1e-3
    # base_name = 'czm_test'
    # lag_mode_mixity = true
    # lag_displacement_jump = true
  []

[mpai60]

Understood, thank you for the clarification! Yes, I am using BreakMeshByBlockGenerator to create the opposing interfaces, and I have the split_interface set to false. This seems to work for node-face contact, and I can use the CZM action to call that boundary.

However, I can set split_interface to true, meaning there's a pair of sidesets for each contact interface. I have added an exodus file to the original github repo (binder_in.e) showing the setup with split_interface=true.

If I'm understanding you correctly, the automatic pairing with mortar contact is not fully implemented? I would be very interested in that feature. I will play with the ContactAction pairing ability for now, please let me know if there's anything I can send over to help.

[lindsayad]

I've opened #33129. It will need to wait on a couple of libMesh pull requests to make their way into MOOSE first before the MOOSE PR can be merged. This puts mortar on par with node-face in terms of ability to perform automatic boundary pairing (and it automatically constructs the lower-D subdomains/LM variables etc.)

[mpai60]

Thank you!

[mpai60]

@lindsayad just wanted to check one other thing, I see that multiple contact pairs are not currently usable with the mortar formulation in the contact action, would the above PR add that functionality as well or would that need to be done separately?

[lindsayad]

Yes that PR is adding support for multiple mortar pairs within a single contact action as well

[hugary1995]

The CZM action issue appears to be identical to this one: #33062 (comment) The fix is discussed in that thread as well.

[mpai60]

I believe I have attempted this already. In the following input, I use just one sideset from each interface in the breakmesh. This method of interface boundary selection is commented out in the czm_action.i file in the original repo I linked. It seems to give me the same results as the original. Am I implementing it incorrectly here?

Input file

[Mesh]

  construct_side_list_from_node_list=true

  [tess]
    type = FileMeshGenerator
    file = 'test_3_refine.msh'
    #file=  'oriented.msh'
  []

  [right_boundary]
    type = BoundingBoxNodeSetGenerator
    input = tess
    new_boundary = 560
    bottom_left = '0.9999999 -0.2 0'
    top_right = '1.2 1.2 0'
  []

  [bottom_boundary]
    type = BoundingBoxNodeSetGenerator
    input = right_boundary
    new_boundary = 559
    bottom_left = '-0.2 -0.2 0'
    top_right = '1.2 0.00015 0'
  []

  [breakmesh]
    type = BreakMeshByBlockGenerator
    input = bottom_boundary
    #split_interface = false
    #interface_name = czm_boundary
  []
[]

[UserObjects]
  [euler_angles]
    type = PropertyReadFile
    nprop = 3
    prop_file_name = 'euler_3.txt'
    read_type = 'block'
    nblock = 3
    use_zero_based_block_indexing = false
  []

[]


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

[Variables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
  []
  [disp_y]
    order = FIRST
    family = LAGRANGE
  []
[]

[AuxVariables]
  [temp]
    order = CONSTANT
    family = MONOMIAL
  []
  [irr]
    order = CONSTANT
    family = MONOMIAL
  []
  [m_ec]
    order = CONSTANT
    family = MONOMIAL
  []
  [m_ea]
    order = CONSTANT
    family = MONOMIAL
  []
  [m_vc]
    order = CONSTANT
    family = MONOMIAL
  []
  [CTE_0]
    order = CONSTANT
    family = MONOMIAL
  []
[]

###################################################################################
[Functions]

  [temp_def]
    type = ConstantFunction
    value = 800
  []
  [irr_def]
    type = ConstantFunction
    value = 2
  []
  [m_ec_def]
    type = ConstantFunction
    value = -0.01
  []
  [m_ea_def]
    type = ConstantFunction
    value = 0.05
  []
  [m_vc_def]
    type = ConstantFunction
    value = -0.01
  []
  [CTE_0_def]
    type = ConstantFunction
    value = 1.3e-5
  []
[]



###################################################################################
[Physics]

  [SolidMechanics]

    [QuasiStatic]
      [all]
        block = '1 2 3' 
        eigenstrain_names = 'thermal_strain irr_strain'
        add_variables = true
        generate_output = 'vonmises_stress'
      []
    []
    
    [CohesiveZone]
      [all]
        boundary = '561 562 564' 
        strain = SMALL
        #verbose = true
        #base_name = czm_test
      []
    []  
  []
[]

###################################################################################
[AuxKernels]
  [tempfuncaux]
    type = FunctionAux
    variable = temp
    function = temp_def
    use_displaced_mesh = false
  []
  [irrfuncaux]
    type = FunctionAux
    variable = irr
    function = irr_def
    use_displaced_mesh = false
  []
  [cte0funcaux]
    type = FunctionAux
    variable = CTE_0
    function = CTE_0_def
    use_displaced_mesh = false
  []
  [mvcfuncaux]
    type = FunctionAux
    variable = m_vc
    function = m_vc_def
    use_displaced_mesh = false
  []
  [mecfuncaux]
    type = FunctionAux
    variable = m_ec
    function = m_ec_def
    use_displaced_mesh = false
  []
  [meafuncaux]
    type = FunctionAux
    variable = m_ea
    function = m_ea_def
    use_displaced_mesh = false
  []
[]

###################################################################################
[BCs]
  [right]
    type = DirichletBC
    boundary = 560
    variable = disp_x
    value = 0.
  []

  [bottom]
    type = DirichletBC
    boundary = 559
    variable = disp_y
    value = 0.
  []  
[]



###################################################################################
[Materials]
  [therm_prefactor]
    type = DerivativeParsedMaterial
    coupled_variables = 'CTE_0 temp'
    property_name = therm_prefactor
    constant_names = 'm_cte T'
    constant_expressions = '2.6e-5 298' #'2.65e-5 298'
    expression = 'CTE_0*(temp-T)' #'(CTE_0*((1+m_cte)*irr)*(temp-T))'
  []
  [thermal_strain]
    type = ComputeVariableEigenstrain
    eigen_base = '-0.0577 0 0 0 1 0 0 0 1'
    args = 'temp'
    prefactor = therm_prefactor
    eigenstrain_name = thermal_strain
  []
  [irr_prefactor]
    type = DerivativeParsedMaterial
    #block = '100'
    coupled_variables = 'irr'
    property_name = irr_prefactor
    constant_names = 'm'
    constant_expressions = '1.185'
    expression = '((m*irr)/100)'
  []
  [irr_strain]
    type = ComputeVariableEigenstrain
    eigen_base = '-0.31 0 0 0 1 0 0 0 1'
    args = 'irr'
    prefactor = irr_prefactor
    eigenstrain_name = irr_strain
  []
  [elasticity_tensor]
    type = ComputeElasticityTensorCP
    block = '1 2 3' 
    C_ijkl = '1.095e12 3.65e10 1.095e12 2.8568e8 9.549e6 9.549e6 0.01 0.01 0.3 0.3 0.01 0.01' 
    fill_method = orthotropic
    read_prop_user_object = euler_angles
  []

  [czm]
    type = BiLinearMixedModeTraction
    boundary = '561 562 564'
    penalty_stiffness = 1e13
    GI_c = 1e5
    GII_c = 1e2
    normal_strength = 1e6
    shear_strength = 1e3
    displacements = 'disp_x disp_y'
    eta = 2.2
    viscosity = 1e-3
    # base_name = 'czm_test'
    # lag_mode_mixity = true
    # lag_displacement_jump = true
  []

  [stress]
    type = ComputeLinearElasticStress
  []
[]

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


[Executioner]

  type = Transient

  solve_type = 'NEWTON'
  line_search = none

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  automatic_scaling = true
  dt = 1
  end_time = 1

[]

[Outputs]
  exodus = true
  csv = true
[]

[mpai60]

Just wanted to follow up on this, is my implementation here incorrect or would I have to use mortar CZMs to capture the interaction of surfaces sliding parallel to each other?