Dynamic Phase-field Fracture in RACCOON

[zuckarm27]

hi all,
i am trying to solve some benchmark problems mentioned in the following literature: https://arxiv.org/abs/2411.16393
i am seeing errors in almost all the problems, i am using the tutorial input files and modifying them to include DynamicStressDivergenceTensors for dynamic problems.
the phase-field variable is getting solved only on a particular boundary where the displacement/traction is applied, even after providing refinement in mesh to facilitate crack propagation.
am i missing something? these problems can be build upon the tutorials in RACCOON, right?

also, i wanted to apply IC on strain history field in the vicinity of the crack. is it possible to do so in RACCOON?

[TheGreatCid]

All the results for the paper you cited were generated in RACCOON ha
I need more details about the type of errors you are seeing

What do you mean by an IC on a strain history field? If you mean initializing a state variable, you might need to go into the source code and code something in.

[zuckarm27]

okay, thank you, i will try this out

in mechanical and fracture file, where we are specifying ICs for d it seems to not work if i provide a mesh file using FileMeshGenerator, it works only when i use the MeshGenerators in MOOSE. is this correct or am i missing something?

i am trying to provide a Gmsh file, its not necessary for the current problem, but for more complex geometries i would prefer to use Gmsh

[TheGreatCid]

It should work. Could you provide more details?

[zuckarm27]

i just changed the Mesh block to provide a Gmsh file.

E = 32e3 #MPa
nu = 0.2
Gc = 0.003 # N/mm
l = 0.375 # mm
rho = 2.45e3 #micrograms/mm^3 

K = '${fparse E/3/(1-2*nu)}'
G = '${fparse E/2/(1+nu)}'

Lambda = '${fparse E*nu/(1+nu)/(1-2*nu)}'
sigma_ts = 5 #MPa
sigma_cs = 18 #MPa
sigma_hs = '${fparse 2/3*sigma_ts*sigma_cs/(sigma_cs - sigma_ts)}'

[MultiApps]
  [phase_field_fracture]
    type = TransientMultiApp
    input_files = 'edge_notch_test_github_fracture.i'
    cli_args = 'E=${E};K=${K};G=${G};Lambda=${Lambda};Gc=${Gc};l=${l};sigma_ts=${sigma_ts};sigma_hs=${sigma_hs}'
    execute_on = 'TIMESTEP_END' 
  []
[]

[Transfers]
  [from_sub_app]
    type = MultiAppCopyTransfer
    from_multi_app = 'phase_field_fracture'
    variable = d
    source_variable = d
    execute_on = 'TIMESTEP_END'
  []

  [to_main_app]
    type = MultiAppCopyTransfer
    to_multi_app = 'phase_field_fracture'
    variable = 'disp_x disp_y strain_zz psie_active'
    source_variable = 'disp_x disp_y strain_zz psie_active'
    # execute_on = MULTIAPP_FIXED_POINT_BEGIN 
  []
 []

[GlobalParams]
    displacements = 'disp_x disp_y'
[]

[Mesh]
  [gen]
    type = FileMeshGenerator
    file = edge_notch_test_dynamic.msh
  []
[]
    

[Variables]
  [disp_x]
  []

  [disp_y]
  []

  [strain_zz]
  []
[]

[AuxVariables]
  [reac_x]
  []

  [reac_y]
  []

  [d]
    [InitialCondition]
      type = FunctionIC
      function = 'if(y=0&x>=0&x<=20,1,0)'
    []  
  []

  [vel_x]
  []

  [accel_x]
  []

  [vel_y]
  []

  [accel_y]
  []
[]

[Kernels]
  [dynamic_x]
    type = ADDynamicStressDivergenceTensors
    variable = disp_x
    component = 0
    save_in = reac_x
    alpha = -0.25
  []

  [dynamic_y]
    type = ADDynamicStressDivergenceTensors
    variable = disp_y
    component = 1
    save_in = reac_y
    alpha = -0.25
  []

  [inertia_x]
    type = ADInertialForce
    variable = disp_x
    beta = 0.39
    gamma = 0.75
    velocity = vel_x
    acceleration = accel_x
    density = ad_density
  []

  [inertia_y]
    type = ADInertialForce
    variable = disp_y
    beta = 0.39
    gamma = 0.75
    velocity = vel_y
    acceleration = accel_y
        density = ad_density

  []

  [plane_stress]
    type = ADWeakPlaneStress
    variable = 'strain_zz'
    displacements = 'disp_x disp_y'
  []
[]

[AuxKernels]
  [accel_x]
    type = TestNewmarkTI
    variable = accel_x
    displacement = disp_x
    first = false
  []

  [vel_x]
    type = TestNewmarkTI
    variable = vel_x
    displacement = disp_x
  []

  [accel_y]
    type = TestNewmarkTI
    variable = accel_y
    displacement = disp_y
    first = false
  []

  [vel_y]
    type = TestNewmarkTI
    variable = vel_y
    displacement = disp_y
  []
[]

[Functions]
	[func]
		type = PiecewiseLinear
    x = '0  20  40'
    y = '0  0.0015   0.0015'
	[]
  [bottom_func] # added bottom func
		type = PiecewiseLinear
    x = '0  20  40'
    y = '0  -0.0015   -0.0015'
	[]
[]

[BCs]
	[top_disp]
		type = FunctionDirichletBC
		variable = disp_y
		boundary = top
		function = func
	[]

  [bottom_disp] # using bottom func
		type = FunctionDirichletBC
		variable = disp_y
		boundary = bottom
		function = bottom_func
	[]
[]

[Materials]
  [bulk_mat_props]
    type = ADGenericConstantMaterial
    prop_names = 'E K G lambda Gc l'
    prop_values = '${E} ${K} ${G} ${Lambda} ${Gc} ${l}'
  []

  [density]
    type = GenericConstantMaterial
    prop_names = 'density'
    prop_values = '${rho}'
  []

  [ad_density]
    type = ADGenericConstantMaterial
    prop_names = 'ad_density'
    prop_values = '${rho}'
  []

  [degradation_func]
    type = PowerDegradationFunction
    property_name = g
    expression = (1-d)^p*(1-eta)+eta
    phase_field = d
    parameter_names = 'p eta'
    parameter_values = '2 0'
  []

  [strain_calc]
    type = ADComputePlaneSmallStrain
    displacements = 'disp_x disp_y'
    out_of_plane_strain = strain_zz
  []

  [elasticity]
    type = SmallDeformationIsotropicElasticity
    bulk_modulus = K
    shear_modulus = G
    phase_field = d
    degradation_function = g
    decomposition = NONE
    output_properties = 'elastic_strain	psie_active psie'
    outputs = exodus
  []

  [stress_calc]
    type = ComputeSmallDeformationStress
    elasticity_model = elasticity
    output_properties = 'stress'
    outputs = exodus
  []
[]

[Postprocessors]
  [Kinetic_Energy]
    type = KineticEnergy
    displacements = 'disp_x disp_y'
    density = ad_density
  []

  [External_Work]
    type = ExternalWork
    displacements = 'disp_x disp_y'
    forces = 'reac_x reac_y'
  []

  [Strain_Energy]
    type = ADElementIntegralMaterialProperty
    mat_prop = psie
  []

  [Dissipated_Energy]
    type = ParsedPostprocessor
    expression = 'External_Work-Kinetic_Energy-Strain_Energy'
    pp_names = 'External_Work Kinetic_Energy Strain_Energy'
  []

  [Fy]
    type = NodalSum
    variable = reac_y
  []
[]

[Executioner]
  type = Transient
  start_time = 0
  end_time = 50 #microseconds
  dt = 0.1 #microseconds
  dtmin = 1e-15
  # timestep_tolerance = 1e-8

  [TimeIntegrator]
    type = NewmarkBeta
  []
  fixed_point_min_its = 2
  fixed_point_max_its = 500

  fixed_point_rel_tol = 1e-6
  fixed_point_abs_tol = 1e-8

  nl_max_its = 500
  l_max_its = 500

  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-10
  l_tol = 1e-8
  l_abs_tol = 1e-10
[]

[Outputs]
  exodus = true
  csv = false
  file_base = outputs/github/out
  perf_graph = true
  checkpoint = false
  print_linear_residuals = false
[]

[Mesh]
  [gen]
    type = FileMeshGenerator
    file = edge_notch_test_dynamic.msh
  []
[]


[Variables]
  [d]
    [InitialCondition]
      type = FunctionIC
      function = 'if(y=0&x>=0&x<=20,1,0)'
    []  
  []
[]

[AuxVariables]
  [bounds_dummy]
  []

  [disp_x]
  []

  [disp_y]
  []

  [strain_zz]
  []

  [psie_active]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Bounds]
  [irreversibility]
    type = VariableOldValueBounds
    variable = bounds_dummy
    bounded_variable = d
    bound_type = lower
  []
  
  [upper]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = d
    bound_type = upper
    bound_value = 1
  []
[]

[Kernels]
  [diff]
    type = ADPFFDiffusion
    variable = d
    fracture_toughness = Gc_delta
    regularization_length = l
    normalization_constant = c0
  []

  [source]
    type = ADPFFSource
    variable = d
    free_energy = psi
  []

  [nucleation_force]
    type = ADCoefMatSource
    variable = d
    prop_names = 'ce'
    coefficient = 1.0
  []
[]

[Materials]
  [fracture_properties]
    type = ADGenericConstantMaterial
    prop_names = 'E K G lambda Gc l sigma_ts sigma_hs'
    prop_values = '${E} ${K} ${G} ${Lambda} ${Gc} ${l} ${sigma_ts} ${sigma_hs}'
  []

  [degradation]
    type = PowerDegradationFunction
    property_name = g
    expression = (1-d)^p*(1-eta)+eta
    phase_field = d
    parameter_names = 'p eta '
    parameter_values = '2 0'
  []

  [crack_geometric]
    type = CrackGeometricFunction
    property_name = alpha
    expression = 'd'
    phase_field = d
  []

  [psi]
    type = ADDerivativeParsedMaterial
    property_name = psi
    expression = 'alpha*(Gc*delta/c0/l)+g*psie_active'
    coupled_variables = 'd psie_active'
    material_property_names = 'delta alpha(d) g(d) Gc c0 l'
    derivative_order = 1
    output_properties = 'psi'
    outputs = exodus
  []

  [nucleation_micro_force]
    type = LDLNucleationMicroForce
    phase_field = d
    degradation_function = g
    regularization_length = l
    normalization_constant = c0
    fracture_toughness = Gc
    tensile_strength = sigma_ts
    hydrostatic_strength = sigma_hs
    delta = delta
    h_correction = true
    external_driving_force_name = ce
  []

  [strain]
    type = ADComputePlaneSmallStrain
    displacements = 'disp_x disp_y'
    out_of_plane_strain = strain_zz
  []

  [elasticity]
    type = SmallDeformationIsotropicElasticity
    bulk_modulus = K
    shear_modulus = G
    phase_field = d
    degradation_function = g
    decomposition = NONE
  []

  [stress]
    type = ComputeSmallDeformationStress
    elasticity_model = elasticity
  []

  [Gc_delta]
    type = ADParsedMaterial
    property_name = Gc_delta
    expression = 'Gc*delta'
    material_property_names = 'Gc delta'
  []
[]

[Executioner]
  type = Transient

  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -snes_type'
  petsc_options_value = 'lu       superlu_dist                  vinewtonrsls'
  automatic_scaling = true

  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-10
[]

[Postprocessors]
  [strain_energy]
    type = ElementAverageValue
    variable = psi
  []
[]

[Outputs]
  exodus = true
  csv = true
  print_linear_residuals = false
[]

this is the phase-field at timestep 5

when i do --mesh-only i get the foll output

Mesh Information:
  elem_dimensions()={2}
  elem_default_orders()={FIRST}
  supported_nodal_order()=1
  spatial_dimension()=2
  n_nodes()=42945
    n_local_nodes()=2833
  n_elem()=42784
    n_local_elem()=2716
    n_active_elem()=42784
  n_subdomains()=1
  n_elemsets()=0
  n_partitions()=16
  n_processors()=16
  n_threads()=1
  processor_id()=0
  is_prepared()=true
  is_replicated()=true

 Mesh Bounding Box:
  Minimum: (x,y,z)=(       0,      -20,        0)
  Maximum: (x,y,z)=(      40,       20,        0)
  Delta:   (x,y,z)=(      40,       40,        0)

 Mesh Element Type(s):
  QUAD4

 Mesh Nodesets:
  Nodeset 0, 4 nodes
   Bounding box minimum: (x,y,z)=(       0,      -20,        0)
   Bounding box maximum: (x,y,z)=(      40,       20,        0)
   Bounding box delta: (x,y,z)=(      40,       40,        0)
  Nodeset 1, 41 nodes
   Bounding box minimum: (x,y,z)=(       0,      -20,        0)
   Bounding box maximum: (x,y,z)=(      40,      -20,        0)
   Bounding box delta: (x,y,z)=(      40,        0,        0)
  Nodeset 2, 121 nodes
   Bounding box minimum: (x,y,z)=(      40,      -20,        0)
   Bounding box maximum: (x,y,z)=(      40,       20,        0)
   Bounding box delta: (x,y,z)=(       0,       40,        0)
  Nodeset 3, 41 nodes
   Bounding box minimum: (x,y,z)=(       0,       20,        0)
   Bounding box maximum: (x,y,z)=(      40,       20,        0)
   Bounding box delta: (x,y,z)=(      40,        0,        0)
  Nodeset 4, 121 nodes
   Bounding box minimum: (x,y,z)=(       0,      -20,        0)
   Bounding box maximum: (x,y,z)=(       0,       20,        0)
   Bounding box delta: (x,y,z)=(       0,       40,        0)

 Mesh Sidesets:
  Sideset 1 (bottom), 40 sides (EDGE2), 40 elems (QUAD4), 41 nodes
   Side volume: 40
   Bounding box minimum: (x,y,z)=(       0,      -20,        0)
   Bounding box maximum: (x,y,z)=(      40,      -20,        0)
   Bounding box delta: (x,y,z)=(      40,        0,        0)
  Sideset 2 (right), 120 sides (EDGE2), 120 elems (QUAD4), 121 nodes
   Side volume: 40
   Bounding box minimum: (x,y,z)=(      40,      -20,        0)
   Bounding box maximum: (x,y,z)=(      40,       20,        0)
   Bounding box delta: (x,y,z)=(       0,       40,        0)
  Sideset 3 (top), 40 sides (EDGE2), 40 elems (QUAD4), 41 nodes
   Side volume: 40
   Bounding box minimum: (x,y,z)=(       0,       20,        0)
   Bounding box maximum: (x,y,z)=(      40,       20,        0)
   Bounding box delta: (x,y,z)=(      40,        0,        0)
  Sideset 4 (left), 120 sides (EDGE2), 120 elems (QUAD4), 121 nodes
   Side volume: 40
   Bounding box minimum: (x,y,z)=(       0,      -20,        0)
   Bounding box maximum: (x,y,z)=(       0,       20,        0)
   Bounding box delta: (x,y,z)=(       0,       40,        0)

 Mesh Edgesets:
  None

 Mesh Subdomains:
  Subdomain 0: 42784 elems (QUAD4, 42784 active), 42945 active nodes
   Volume: 1600.02
   Bounding box minimum: (x,y,z)=(       0,      -20,        0)
   Bounding box maximum: (x,y,z)=(      40,       20,        0)
   Bounding box delta: (x,y,z)=(      40,       40,        0)
  Global mesh volume = 1600.02

edge_notch_test_dynamic.zip

[TheGreatCid]

I think what is happening is that you are defining a intial condition over a line y==0. MOOSE by default generates a structured mesh which has nodes along that line. Your gmsh mesh might not have a set of nodes along y==0 (they might be at y==0.1 for example). You can check this in paraview by inspecting the nodes.
Try defining a bounding box instead
Increase the y thickness until you see an initial damage band.

[zuckarm27]

thanks