Fixing MeshCut2DFractureUserObject to use just stress based or consta…

…nt Kcrit or kcrit from vpp crack growth. Added error checking and testing for these options.
idaholab · Dec 20, 2024 · 8da5e53 · 8da5e53
1 parent 3740102
commit 8da5e53
Show file tree

Hide file tree

Showing 11 changed files with 77 additions and 95 deletions.
diff --git a/modules/xfem/include/userobjects/MeshCut2DFractureUserObject.h b/modules/xfem/include/userobjects/MeshCut2DFractureUserObject.h
@@ -31,15 +31,6 @@ class MeshCut2DFractureUserObject : public MeshCut2DUserObjectBase
   virtual void findActiveBoundaryGrowth() override;
 
 private:
-  /**
-   * Compute all of the maximum hoop stress fracture integrals for all crack trips from the fracture
-   * integral vector post processors
-   * @param k1 fracture integrals from KI vector postprocessors
-   * @param k2 fracture integrals from KII vector postprocessors
-   * @return computed fracture integral squared
-   */
-  std::vector<Real> getKSquared(const std::vector<Real> & k1, const std::vector<Real> & k2) const;
-
   /// amount to grow crack by for each xfem update step
   const Real & _growth_increment;
 

diff --git a/modules/xfem/src/userobjects/MeshCut2DFractureUserObject.C b/modules/xfem/src/userobjects/MeshCut2DFractureUserObject.C
@@ -56,7 +56,8 @@ MeshCut2DFractureUserObject::MeshCut2DFractureUserObject(const InputParameters &
     _growth_increment(getParam<Real>("growth_increment")),
     _use_k(isParamValid("k_critical") || isParamValid("k_critical_vectorpostprocessor")),
     _use_stress(isParamValid("stress_threshold")),
-    _k_critical(_use_k ? getParam<Real>("k_critical") : std::numeric_limits<Real>::max()),
+    _k_critical(isParamValid("k_critical") ? getParam<Real>("k_critical")
+                                           : std::numeric_limits<Real>::max()),
     _stress_threshold(_use_stress ? getParam<Real>("stress_threshold")
                                   : std::numeric_limits<Real>::max()),
     _ki_vpp(_use_k ? &getVectorPostprocessorValue(
@@ -81,6 +82,11 @@ MeshCut2DFractureUserObject::MeshCut2DFractureUserObject(const InputParameters &
     paramError("k_critical",
                "Must set crack extension criterion with k_critical, k_critical_vectorpostprocessor "
                "or stress_threshold.");
+
+  if (isParamValid("k_critical") && isParamValid("k_critical_vectorpostprocessor"))
+    paramError("k_critical",
+               "Fracture toughness cannot be specified by both k_critical and "
+               "k_critical_vectorpostprocessor.");
 }
 
 void
@@ -133,34 +139,36 @@ MeshCut2DFractureUserObject::findActiveBoundaryGrowth()
         "\n  cracktips in MeshCut2DFractureUserObject = ",
         _original_and_current_front_node_ids.size());
 
-  std::vector<Real> k_squared = getKSquared(*_ki_vpp, *_kii_vpp);
   _active_front_node_growth_vectors.clear();
   for (unsigned int i = 0; i < _original_and_current_front_node_ids.size(); ++i)
   {
-    Real k_crit = _k_critical;
-    if (_k_critical_vpp)
-      k_crit = std::min(_k_critical_vpp->at(i), _k_critical);
-
-    if (_use_k && k_squared[i] > (k_crit * k_crit) && _ki_vpp->at(i) > 0)
+    if (_use_k)
     {
-      // growth direction in crack front coord (cfc) system based on the  max hoop stress
-      // criterion
-      Real ki = _ki_vpp->at(i);
-      Real kii = _kii_vpp->at(i);
-      Real sqrt_k = std::sqrt(ki * ki + kii * kii);
-      Real theta = 2 * std::atan((ki - sqrt_k) / (4 * kii));
-      RealVectorValue dir_cfc;
-      dir_cfc(0) = std::cos(theta);
-      dir_cfc(1) = std::sin(theta);
-      dir_cfc(2) = 0;
-
-      // growth direction in global coord system based on the max hoop stress criterion
-      RealVectorValue dir_global =
-          _crack_front_definition->rotateFromCrackFrontCoordsToGlobal(dir_cfc, i);
-      Point dir_global_pt(dir_global(0), dir_global(1), dir_global(2));
-      Point nodal_offset = dir_global_pt * _growth_increment;
-      _active_front_node_growth_vectors.push_back(
-          std::make_pair(_original_and_current_front_node_ids[i].second, nodal_offset));
+      Real k_crit = _k_critical;
+      if (_k_critical_vpp)
+        k_crit = std::min(_k_critical_vpp->at(i), _k_critical);
+      Real k_squared = _ki_vpp->at(i) * _ki_vpp->at(i) + _kii_vpp->at(i) * _kii_vpp->at(i);
+      if (k_squared > (k_crit * k_crit) && _ki_vpp->at(i) > 0)
+      {
+        // growth direction in crack front coord (cfc) system based on the  max hoop stress
+        // criterion
+        Real ki = _ki_vpp->at(i);
+        Real kii = _kii_vpp->at(i);
+        Real sqrt_k = std::sqrt(ki * ki + kii * kii);
+        Real theta = 2 * std::atan((ki - sqrt_k) / (4 * kii));
+        RealVectorValue dir_cfc;
+        dir_cfc(0) = std::cos(theta);
+        dir_cfc(1) = std::sin(theta);
+        dir_cfc(2) = 0;
+
+        // growth direction in global coord system based on the max hoop stress criterion
+        RealVectorValue dir_global =
+            _crack_front_definition->rotateFromCrackFrontCoordsToGlobal(dir_cfc, i);
+        Point dir_global_pt(dir_global(0), dir_global(1), dir_global(2));
+        Point nodal_offset = dir_global_pt * _growth_increment;
+        _active_front_node_growth_vectors.push_back(
+            std::make_pair(_original_and_current_front_node_ids[i].second, nodal_offset));
+      }
     }
     else if (_use_stress && _stress_vpp->at(i) > _stress_threshold)
     {
@@ -175,14 +183,3 @@ MeshCut2DFractureUserObject::findActiveBoundaryGrowth()
     }
   }
 }
-
-std::vector<Real>
-MeshCut2DFractureUserObject::getKSquared(const std::vector<Real> & k1,
-                                         const std::vector<Real> & k2) const
-{
-  std::vector<Real> k_squared(k1.size());
-  for (unsigned int i = 0; i < k_squared.size(); ++i)
-    k_squared[i] = k1[i] * k1[i] + k2[i] * k2[i];
-
-  return k_squared;
-}
diff --git a/modules/xfem/test/tests/mesh_cut_2D_fracture/edge_crack_2d_propagation.i b/modules/xfem/test/tests/mesh_cut_2D_fracture/edge_crack_2d_propagation.i
@@ -46,39 +46,17 @@
   used_by_xfem_to_grow_crack = true
 []
 
-[VectorPostprocessors]
-  [CrackFrontNonlocalStressVpp]
-    type = CrackFrontNonlocalStress
-    material_name = stress
-    crack_front_definition = crackFrontDefinition
-    box_length = 0.05
-    box_height = 0.1
-    execute_on = NONLINEAR
-  []
-  [CrackFrontNonlocalScalarVpp]
-    type = CrackFrontNonlocalScalar
-    material_name = k_crit_mat
-    crack_front_definition = crackFrontDefinition
-    box_length = 0.05
-    box_height = 0.1
-    execute_on = NONLINEAR
-  []
-[]
-[UserObjects]
-  [cut_mesh2]
-    type = MeshCut2DFractureUserObject
-    mesh_file = make_edge_crack_in.e
-    growth_increment = 0.05
-    ki_vectorpostprocessor = "II_KI_1"
-    kii_vectorpostprocessor = "II_KII_1"
-    k_critical = 100
-    k_critical_vectorpostprocessor=CrackFrontNonlocalScalarVpp
-    k_critical_vector_name = "crack_tip_k_crit_mat"
-    # stress_vectorpostprocessor = "CrackFrontNonlocalStressVpp"
-    # stress_vector_name = "crack_tip_stress"
-    # stress_threshold = 120
-  []
-[]
+# MeshCut2DFractureUserObject are included in seperate input files
+# [UserObjects]
+#   [cut_mesh2]
+#     type = MeshCut2DFractureUserObject
+#     mesh_file = make_edge_crack_in.e
+#     growth_increment = 0.05
+#     ki_vectorpostprocessor = "II_KI_1"
+#     kii_vectorpostprocessor = "II_KII_1"
+#     k_critical = 100
+#   []
+# []
 
 [Physics/SolidMechanics/QuasiStatic]
   [all]
@@ -127,12 +105,6 @@
   [stress]
     type = ComputeFiniteStrainElasticStress
   []
-  [k_critical]
-    type = ParsedMaterial
-    property_name = k_crit_mat
-    extra_symbols = 'x'
-    expression = 'if(x > -0.4,100,90)'
-  []
 []
 
 [Executioner]
@@ -141,6 +113,7 @@
   petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -pc_factor_shift_type -pc_factor_shift_amount'
   petsc_options_value = ' lu       superlu_dist                 NONZERO               1e-20'
   line_search = 'none'
+  nl_abs_tol = 1e-7
   [Predictor]
     type = SimplePredictor
     scale = 1.0
@@ -153,7 +126,6 @@
 []
 
 [Outputs]
-  exodus = true
   [xfemcutter]
     type = XFEMCutMeshOutput
     xfem_cutter_uo = cut_mesh2

diff --git a/modules/xfem/test/tests/mesh_cut_2D_fracture/gold/edge_crack_2d_propagation_out.e-s002 b/modules/xfem/test/tests/mesh_cut_2D_fracture/gold/edge_crack_2d_propagation_out.e-s002
diff --git a/modules/xfem/test/tests/mesh_cut_2D_fracture/gold/k_stress_growth.e-s003 b/modules/xfem/test/tests/mesh_cut_2D_fracture/gold/k_stress_growth.e-s003
diff --git a/modules/xfem/test/tests/mesh_cut_2D_fracture/gold/k_stress_growth_XFEMCutMeshOutput.e-s0004 b/modules/xfem/test/tests/mesh_cut_2D_fracture/gold/k_stress_growth_XFEMCutMeshOutput.e-s0004
diff --git a/...e_crack_2d_propagation_xfemcutter.e-s0002 → ...t_uo_xfemcutter_XFEMCutMeshOutput.e-s0004 b/...e_crack_2d_propagation_xfemcutter.e-s0002 → ...t_uo_xfemcutter_XFEMCutMeshOutput.e-s0004
diff --git a/...h_cut_2D_fracture/gold/kcrit_stress_based_meshCut_uo_xfemcutter_XFEMCutMeshOutput.e-s0004 b/...h_cut_2D_fracture/gold/kcrit_stress_based_meshCut_uo_xfemcutter_XFEMCutMeshOutput.e-s0004
diff --git a/...ests/mesh_cut_2D_fracture/gold/kvpp_based_meshCut_uo_xfemcutter_XFEMCutMeshOutput.e-s0004 b/...ests/mesh_cut_2D_fracture/gold/kvpp_based_meshCut_uo_xfemcutter_XFEMCutMeshOutput.e-s0004
diff --git a/...tion_xfemcutter_XFEMCutMeshOutput.e-s0004 → ...t_uo_xfemcutter_XFEMCutMeshOutput.e-s0004 b/...tion_xfemcutter_XFEMCutMeshOutput.e-s0004 → ...t_uo_xfemcutter_XFEMCutMeshOutput.e-s0004
diff --git a/modules/xfem/test/tests/mesh_cut_2D_fracture/tests b/modules/xfem/test/tests/mesh_cut_2D_fracture/tests
@@ -6,27 +6,49 @@
     'where the crack is defined using a topologically defined surface cutting mesh by MeshCut2DFunctionUserObject'
     [k_growth]
       type = Exodiff
-      input = edge_crack_2d_propagation.i
-      exodiff = 'edge_crack_2d_propagation_xfemcutter_XFEMCutMeshOutput.e-s0004 edge_crack_2d_propagation_out.e-s002'
+      input = 'edge_crack_2d_propagation.i kcrit_based_meshCut_uo.i'
+      exodiff = 'kcrit_based_meshCut_uo_xfemcutter_XFEMCutMeshOutput.e-s0004'
       abs_zero = 1e-8
       map = false
       # XFEM requires --enable-unique-ids in libmesh
       unique_id = true
-      detail = ' and crack growth and direction is computed by a fracture integral.'
+      detail = ' and crack growth and direction is computed by a fracture integral with constant fracture toughness.'
+    []
+    [kvpp_growth]
+      type = Exodiff
+      input = 'edge_crack_2d_propagation.i kvpp_based_meshCut_uo.i'
+      exodiff = 'kvpp_based_meshCut_uo_xfemcutter_XFEMCutMeshOutput.e-s0004'
+      abs_zero = 1e-8
+      map = false
+      # XFEM requires --enable-unique-ids in libmesh
+      unique_id = true
+      detail = ' and crack growth and direction is computed by a fracture integral with fracture toughness sampled from a material at the crack tip.'
+    []
+    [stress_growth]
+      type = Exodiff
+      input = 'edge_crack_2d_propagation.i stress_based_meshCut_uo.i'
+      exodiff = 'stress_based_meshCut_uo_xfemcutter_XFEMCutMeshOutput.e-s0004'
+      abs_zero = 1e-8
+      map = false
+      # XFEM requires --enable-unique-ids in libmesh
+      unique_id = true
+      detail = ' and crack growth is a maximum stress criterion and the growth direction is in the original crack direction.'
     []
     [k_stress_growth]
       type = Exodiff
-      input = edge_crack_2d_propagation.i
-      cli_args = 'Outputs/file_base=k_stress_growth '
-                 'UserObjects/cut_mesh2/stress_vectorpostprocessor=CrackFrontNonlocalStressVpp '
-                 'UserObjects/cut_mesh2/stress_vector_name="crack_tip_stress" '
-                 'UserObjects/cut_mesh2/stress_threshold=120'
-      exodiff = 'k_stress_growth.e-s003 k_stress_growth_XFEMCutMeshOutput.e-s0004'
+      input = 'edge_crack_2d_propagation.i kcrit_stress_based_meshCut_uo.i'
+      exodiff = 'kcrit_stress_based_meshCut_uo_xfemcutter_XFEMCutMeshOutput.e-s0004'
       abs_zero = 1e-8
       map = false
       # XFEM requires --enable-unique-ids in libmesh
       unique_id = true
-      detail = ' and crack growth and direction is computed by a fracture integral or maximum stress criterion.'
+      detail = ' and crack growth and direction is computed by a fracture integral with constant fracture toughness or maximum stress criterion.'
+    []
+    [k_error]
+      type = 'RunException'
+      input = 'edge_crack_2d_propagation.i kvpp_based_meshCut_uo.i'
+      cli_args = "UserObjects/cut_mesh2/k_critical=100"
+      expect_err = 'Fracture toughness cannot be specified by both k_critical and k_critical_vectorpostprocessor.'
     []
     [stress_check]
       type = CSVDiff