Student needs help with Explicit contact: Dynamics

[Rkmoffat]

Hello,
I am working on a model where a particle impacts a thin film. When running implicit, the film deforms around the particle, as I would expect. But upon switching to explicit, there's a strange gap that is created between the two surfaces. Followed by element deformation on the particle that leads to failure. The particle is very stiff and should not be deforming. Changing the primary face causes the issue to occur on the other face.

Thank you for any and all advice.

[GlobalParams] displacements = 'disp_x disp_y disp_z' []
[Mesh]
  [plate_file] type = FileMeshGenerator  file = '7_PLATE_QuartFlyerThick_50K.msh' []
  [flyer_file] type = FileMeshGenerator  file = '7_FLYER_QuartFlyerThick_50K.msh' []
  [plate]
    type = SubdomainIDGenerator
    input = plate_file
    subdomain_id = 1
  []
  [flyer]
    type = SubdomainIDGenerator
    input = flyer_file
    subdomain_id = 2
  []
  [mesh]    type   = CombinerGenerator    inputs = 'plate flyer'  []
  [flyer_surface]
    type   = SideSetsAroundSubdomainGenerator 
    input = mesh
    new_boundary = 'flyer_surface'
    normal = '0 0 -1'
    normal_tol = 1
    block = 2
  []
[]

[ExplicitDynamicsContact]
    [my_contact]
        model = frictionless_balance
    	primary = 'flyer_surface'
   	secondary = '1002'
        vel_x = 'vel_x'
        vel_y = 'vel_y'
        vel_z = 'vel_z'
    []
[]
[Variables]
  [./disp_x]  order = FIRST family = LAGRANGE [../]
  [./disp_y]  order = FIRST family = LAGRANGE [../]
  [./disp_z]  order = FIRST family = LAGRANGE [../]
[]
[Kernels]
  ##Mechanics##
    [./inertia_x]  type = InertialForce  variable = disp_x [../]
    [./inertia_y]  type = InertialForce  variable = disp_y [../]
    [./inertia_z]  type = InertialForce  variable = disp_z [../]
    [./DynamicTensorMechanics]
      displacements = 'disp_x disp_y disp_z'
      use_displaced_mesh = true
    [../]
    [gravity]
        type = Gravity
        variable = disp_z
        value = 1
	function = VeloWant
        block = '2' #Flyer
        use_displaced_mesh=true
        implicit=false
    []
[]
[Functions]
  [VeloWant]    type = ParsedFunction    expression = 'if(t=0,1,0)*-0.35/(0.0001)'  []
[]
[AuxVariables]
 #Calculating
 #Reporting
  [./min]            order = CONSTANT family = MONOMIAL [../]
  [./vel_x]          order = FIRST    family = LAGRANGE [../]
  [./accel_x]        order = FIRST    family = LAGRANGE [../]
  [./vel_y]          order = FIRST    family = LAGRANGE [../]
  [./accel_y]        order = FIRST    family = LAGRANGE [../]
  [./vel_z]          order = FIRST    family = LAGRANGE [../]
  [./accel_z]        order = FIRST    family = LAGRANGE [../]
  [./stressxx]       order = CONSTANT family = MONOMIAL [../]
  [./stressyy]       order = CONSTANT family = MONOMIAL [../]
  [./stresszz]       order = CONSTANT family = MONOMIAL [../]
  [./stressxy]       order = CONSTANT family = MONOMIAL [../]
  [./stressyz]       order = CONSTANT family = MONOMIAL [../]
  [./stressxz]       order = CONSTANT family = MONOMIAL [../]
  [./mechstrainxx]   order = CONSTANT family = MONOMIAL [../]
  [./mechstrainyy]   order = CONSTANT family = MONOMIAL [../]
  [./mechstrainzz]   order = CONSTANT family = MONOMIAL [../]
  [./mechstrainxy]   order = CONSTANT family = MONOMIAL [../]
  [./mechstrainyz]   order = CONSTANT family = MONOMIAL [../]
  [./mechstrainxz]   order = CONSTANT family = MONOMIAL [../]
  [./totstrainxx]    order = CONSTANT family = MONOMIAL [../]
  [./totstrainyy]    order = CONSTANT family = MONOMIAL [../]
  [./totstrainzz]    order = CONSTANT family = MONOMIAL [../]
  [./totstrainxy]    order = CONSTANT family = MONOMIAL [../]
  [./totstrainyz]    order = CONSTANT family = MONOMIAL [../]
  [./totstrainxz]    order = CONSTANT family = MONOMIAL [../]
  [./thstrainxx]     order = CONSTANT family = MONOMIAL [../]
  [./thstrainyy]     order = CONSTANT family = MONOMIAL [../]
  [./thstrainzz]     order = CONSTANT family = MONOMIAL [../]
  [./thstrainxy]     order = CONSTANT family = MONOMIAL [../]
  [./thstrainyz]     order = CONSTANT family = MONOMIAL [../]
  [./thstrainxz]     order = CONSTANT family = MONOMIAL [../]
  [./volstress]      order = CONSTANT family = MONOMIAL [../]
  [./volmechstrain]  order = CONSTANT family = MONOMIAL [../]
  [./vonmisesstress] order = CONSTANT family = MONOMIAL [../]
[]
[AuxKernels] 
 ##Reporting## 
  [min]            type = ElementLengthAux  variable = min    method = min    execute_on = initial  [../]
  [vel_x]          type = TestNewmarkTI     variable = vel_x        displacement = disp_x        execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'    []
  [vel_y]          type = TestNewmarkTI     variable = vel_y        displacement = disp_x        execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'    []
  [vel_z]          type = TestNewmarkTI     variable = vel_z        displacement = disp_z        execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'    []
  [accel_x]        type = TestNewmarkTI     variable = accel_x        displacement = disp_x        first = false        execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'    []
  [accel_y]        type = TestNewmarkTI     variable = accel_y        displacement = disp_y        first = false        execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'    []
  [accel_z]        type = TestNewmarkTI     variable = accel_z        displacement = disp_z        first = false        execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'    []
  [volstress]      type = RankTwoScalarAux  block = '1' variable = volstress      rank_two_tensor = stress            scalar_type = Hydrostatic      []
  [volmechstrain]  type = RankTwoScalarAux  block = '1' variable = volmechstrain  rank_two_tensor = mechanical_strain scalar_type = VolumetricStrain []
  [vonmisesstress] type = RankTwoScalarAux  block = '1' variable = vonmisesstress rank_two_tensor = stress            scalar_type = VonMisesStress   []
  [stressxx]       type = RankTwoAux        block = '1' variable = stressxx       rank_two_tensor = stress            index_j = 0    index_i = 0     []
  [stressyy]       type = RankTwoAux        block = '1' variable = stressyy       rank_two_tensor = stress            index_j = 1    index_i = 1     []
  [stresszz]       type = RankTwoAux        block = '1' variable = stresszz       rank_two_tensor = stress            index_j = 2    index_i = 2     []
  [stressxy]       type = RankTwoAux        block = '1' variable = stressxy       rank_two_tensor = stress            index_j = 1    index_i = 0     []
  [stressyz]       type = RankTwoAux        block = '1' variable = stressyz       rank_two_tensor = stress            index_j = 2    index_i = 1     []
  [stressxz]       type = RankTwoAux        block = '1' variable = stressxz       rank_two_tensor = stress            index_j = 2    index_i = 0     []
  [mechstrainxx]   type = RankTwoAux        block = '1' variable = mechstrainxx   rank_two_tensor = mechanical_strain index_j = 0    index_i = 0     []
  [mechstrainyy]   type = RankTwoAux        block = '1' variable = mechstrainyy   rank_two_tensor = mechanical_strain index_j = 1    index_i = 1     []
  [mechstrainzz]   type = RankTwoAux        block = '1' variable = mechstrainzz   rank_two_tensor = mechanical_strain index_j = 2    index_i = 2     []
  [mechstrainxy]   type = RankTwoAux        block = '1' variable = mechstrainxy   rank_two_tensor = mechanical_strain index_j = 1    index_i = 0     []
  [mechstrainyz]   type = RankTwoAux        block = '1' variable = mechstrainyz   rank_two_tensor = mechanical_strain index_j = 2    index_i = 1     []
  [mechstrainxz]   type = RankTwoAux        block = '1' variable = mechstrainxz   rank_two_tensor = mechanical_strain index_j = 2    index_i = 0     []
  [totstrainxx]    type = RankTwoAux        block = '1' variable = totstrainxx    rank_two_tensor = total_strain      index_j = 0    index_i = 0     []
  [totstrainyy]    type = RankTwoAux        block = '1' variable = totstrainyy    rank_two_tensor = total_strain      index_j = 1    index_i = 1     []
  [totstrainzz]    type = RankTwoAux        block = '1' variable = totstrainzz    rank_two_tensor = total_strain      index_j = 2    index_i = 2     []
  [totstrainxy]    type = RankTwoAux        block = '1' variable = totstrainxy    rank_two_tensor = total_strain      index_j = 1    index_i = 0     []
  [totstrainyz]    type = RankTwoAux        block = '1' variable = totstrainyz    rank_two_tensor = total_strain      index_j = 2    index_i = 1     []
  [totstrainxz]    type = RankTwoAux        block = '1' variable = totstrainxz    rank_two_tensor = total_strain      index_j = 2    index_i = 0     []
[]
################################################################################
[BCs]
  [Plate_Ends_x]
    type = DirichletBC
    variable = disp_x
    boundary = '1001'
    value = 0
  []
  [Plate_Ends_y]
    type = DirichletBC
    variable = disp_y
    boundary = '1001'
    value = 0
  []
  [Plate_Ends_z]
    type = DirichletBC
    variable = disp_z
    boundary = '1001'
    value = 0
  []
[]
[Materials]
  [./material_prop_plate] #Polystyrene
    type = GenericConstantMaterial
    block = '1'
    prop_names  = 'density gc_prop     l  visco'
    prop_values = '  1.050  0.033   0.5    100'
  [../]
  [./material_prop_flyer] #Silica
    type = GenericConstantMaterial
    block = '2'
    prop_names  = 'density'#  gc_prop   l'
    prop_values = '  2.650'#  0.002   0.7 '
  [../]
  [elasticity_tensor_plate]
    type = ComputeIsotropicElasticityTensor
    block= '1'
    youngs_modulus = 3.40#GPa
    poissons_ratio = 0.33  
  []
  [elasticity_tensor_flyer]
    type = ComputeIsotropicElasticityTensor
    block= '2'
    youngs_modulus = 74.8 #GPa
    poissons_ratio = 0.19   
  []
  [./stress_plate] type = ComputeFiniteStrainElasticStress    block = '1'  [../]
  [./stress_flyer] type = ComputeFiniteStrainElasticStress    block = '2'  [../]
  [./strain]       type = ComputeIncrementalSmallStrain  [../]
  [wave_speed]     type = WaveSpeed  []
[]
[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]
[Executioner]
  type = Transient
  start_time = 0.0
  dt      = 0.0005
  end_time = 100 #200
  [TimeIntegrator]
      type = CentralDifference #ActuallyExplicitEuler #CentralDifference #ActuallyExplicitTimeIntergrator  #CentralDifference
      solve_type = consistent#lumped
  []
[]
[Outputs]
  interval = 20
  exodus = true
  #checkpoint = true
  [Csv]    
    type = CSV
    interval = 1  
  []
[] 

[GiudGiud]

Hello

Are these results converged in the time step?
@dschwen @recuero Do you know if anyone worked on contact+explicit?

[TheGreatCid]

Ok so the issue you provided did go away, however there is an issue with a node not being released from contact, shown below as a slice:
https://github.com/user-attachments/assets/a9df3264-9f99-4873-a341-1809d378dda2
and in 3D:

This might go away with mesh refinement. It is hard to say

[TheGreatCid]

This is the input file I used
test.zip

You can use the PR by running

git fetch origin pull/28373/head:CDinSMContactPR 
git checkout CDinSMContactPR  

Then compile again

[Rkmoffat]

Alright!
I will try it with a full-sized mesh then.
Thank you for your time

[TheGreatCid]

No problem! Hopefully it works! If it doesn't let me know. I won't have time to fix it any time soon unfortunately, but It'll be good to know for when I do find time.

[TheGreatCid]

@Rkmoffat Out of curiosity, did the new method work for you after mesh refinement?