XFEM near-tip enrichment for crack propagation

framework/include/interfaces/XFEMInterface.h

@@ -93,6 +93,7 @@ class XFEMInterface : public ConsoleStreamInterface
                       const std::vector<std::shared_ptr<NonlinearSystemBase>> & nl,
                       AuxiliarySystem & aux) = 0;
 
+  virtual bool didNearTipEnrichmentChange() = 0;
   /**
    * Initialize the solution on newly created nodes
    */

framework/src/problems/FEProblemBase.C

@@ -8276,18 +8276,27 @@ FEProblemBase::updateMeshXFEM()
           /*intermediate_change=*/false, /*contract_mesh=*/true, /*clean_refinement_flags=*/false);
 
     updated = _xfem->update(_time, _nl, *_aux);
-    if (updated)
+    if (updated) // if subcrical crack advancement, do not need to reinitialize solution
     {
       meshChanged(
           /*intermediate_change=*/false, /*contract_mesh=*/true, /*clean_refinement_flags=*/false);
       _xfem->initSolution(_nl, *_aux);
       restoreSolutions();
       _console << "\nXFEM update complete: Mesh modified" << std::endl;
     }
+  }
+  bool crack_front_advanced = false;
+  if (!updated)
+  {
+    crack_front_advanced = _xfem->didNearTipEnrichmentChange();
+    if (crack_front_advanced)
+      _console << "\nXFEM update complete: Mesh not modified but advancing crack changed near-tip "
+                  "enrichment"
+               << std::endl;
     else
       _console << "\nXFEM update complete: Mesh not modified" << std::endl;
   }
-  return updated;
+  return (updated || crack_front_advanced);
 }
 
 void

modules/solid_mechanics/include/materials/ComputeIncrementalStrainBase.h

@@ -40,4 +40,5 @@ class ComputeIncrementalStrainBase : public ComputeStrainBase
   const MaterialProperty<RankTwoTensor> & _total_strain_old;
 
   std::vector<const MaterialProperty<RankTwoTensor> *> _eigenstrains_old;
+  MaterialProperty<RankTwoTensor> & _grad_disp_rate;
 };

modules/solid_mechanics/include/userobjects/CrackFrontDefinition.h

@@ -41,6 +41,7 @@ class CrackFrontDefinition : public GeneralUserObject, public BoundaryRestrictab
   virtual void finalize() override;
   virtual void execute() override;
 
+  void updateCrackFrontPoints();
   /// used by Actions to add CrackFrontDefinitionParams
   static void includeCrackFrontDefinitionParams(InputParameters & params);
 

modules/solid_mechanics/src/materials/ComputeIncrementalStrainBase.C

@@ -26,7 +26,8 @@ ComputeIncrementalStrainBase::ComputeIncrementalStrainBase(const InputParameters
     _deformation_gradient(declareProperty<RankTwoTensor>(_base_name + "deformation_gradient")),
     _mechanical_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "mechanical_strain")),
     _total_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "total_strain")),
-    _eigenstrains_old(_eigenstrain_names.size())
+    _eigenstrains_old(_eigenstrain_names.size()),
+    _grad_disp_rate(declareProperty<RankTwoTensor>(_base_name + "grad_disp_rate"))
 {
   for (unsigned int i = 0; i < _eigenstrains_old.size(); ++i)
     _eigenstrains_old[i] = &getMaterialPropertyOld<RankTwoTensor>(_eigenstrain_names[i]);

modules/solid_mechanics/src/userobjects/CrackFrontDefinition.C

@@ -357,25 +357,26 @@ CrackFrontDefinition::initialSetup()
   }
 }
 
+void
+CrackFrontDefinition::updateCrackFrontPoints()
+{
+  _crack_front_points =
+      _crack_front_points_provider->getCrackFrontPoints(_num_points_from_provider);
+  updateCrackFrontGeometry();
+  std::size_t num_crack_front_points = getNumCrackFrontPoints();
+  if (_q_function_type == "GEOMETRY")
+    for (std::size_t i = 0; i < num_crack_front_points; ++i)
+    {
+      bool is_point_on_intersecting_boundary = isPointWithIndexOnIntersectingBoundary(i);
+      _is_point_on_intersecting_boundary.push_back(is_point_on_intersecting_boundary);
+    }
+}
+
 void
 CrackFrontDefinition::initialize()
 {
-  // Update the crack front for fracture integral calculations
-  // This is only useful for growing cracks which are currently described by the mesh
-  // cutter
   if (_use_mesh_cutter && _is_cutter_modified)
-  {
-    _crack_front_points =
-        _crack_front_points_provider->getCrackFrontPoints(_num_points_from_provider);
-    updateCrackFrontGeometry();
-    std::size_t num_crack_front_points = getNumCrackFrontPoints();
-    if (_q_function_type == "GEOMETRY")
-      for (std::size_t i = 0; i < num_crack_front_points; ++i)
-      {
-        bool is_point_on_intersecting_boundary = isPointWithIndexOnIntersectingBoundary(i);
-        _is_point_on_intersecting_boundary.push_back(is_point_on_intersecting_boundary);
-      }
-  }
+    updateCrackFrontPoints();
 }
 
 void

modules/xfem/doc/content/source/materials/ComputeCrackTipEnrichmentIncrementalStrain.md

@@ -0,0 +1,14 @@
+# ComputeCrackTipEnrichmentIncrementalStrain
+
+!syntax description /Materials/ComputeCrackTipEnrichmentIncrementalStrain
+
+# Description
+
+The displacement field contains both the standard and the enrichment solution. The
+`ComputeCrackTipEnrichmentIncrementalStrain` is inherited from '`ComputeCrackTipEnrichmentSmallStrain` to allow nonlinear material models, and is applicable only to small-strain problems.
+
+!syntax parameters /Materials/ComputeCrackTipEnrichmentIncrementalStrain
+
+!syntax inputs /Materials/ComputeCrackTipEnrichmentIncrementalStrain
+
+!syntax children /Materials/ComputeCrackTipEnrichmentIncrementalStrain

modules/xfem/include/base/XFEM.h

@@ -322,6 +322,7 @@ class XFEM : public XFEMInterface
   CutSubdomainID getCutSubdomainID(const GeometricCutUserObject * gcuo,
                                    const Elem * cut_elem,
                                    const Elem * parent_elem = nullptr) const;
+  virtual bool didNearTipEnrichmentChange() override;
 
 private:
   bool _has_secondary_cut;

modules/xfem/include/materials/ComputeCrackTipEnrichmentIncrementalStrain.h

@@ -0,0 +1,76 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+// #include "ComputeStrainBase.h"
+#include "ComputeIncrementalStrainBase.h"
+#include "Material.h"
+#include "RankTwoTensor.h"
+#include "RankFourTensor.h"
+#include "RotationTensor.h"
+#include "Assembly.h"
+#include "CrackFrontDefinition.h"
+#include "EnrichmentFunctionCalculation.h"
+
+/**
+ * ComputeIncrementalStrain defines a strain increment and rotation increment (=1), for small
+ * strains.
+ */
+class ComputeCrackTipEnrichmentIncrementalStrain : public ComputeIncrementalStrainBase,
+                                                   public EnrichmentFunctionCalculation
+{
+public:
+  static InputParameters validParams();
+
+  ComputeCrackTipEnrichmentIncrementalStrain(const InputParameters & parameters);
+  virtual ~ComputeCrackTipEnrichmentIncrementalStrain() {}
+  virtual void computeProperties() override;
+
+protected:
+  /// enrichment displacement
+  std::vector<Real> _enrich_disp;
+
+  /// gradient of enrichment displacement
+  std::vector<RealVectorValue> _grad_enrich_disp;
+  std::vector<RealVectorValue> _grad_enrich_disp_old;
+
+  /// enrichment displacement variables
+  std::vector<std::vector<MooseVariableFEBase *>> _enrich_variable;
+
+  /// the current shape functions
+  const VariablePhiValue & _phi;
+
+  /// gradient of the shape function
+  const VariablePhiGradient & _grad_phi;
+
+  const MaterialProperty<RankTwoTensor> & _mechanical_strain_old;
+  const MaterialProperty<RankTwoTensor> & _total_strain_old;
+  MaterialProperty<RankTwoTensor> & _grad_disp_tensor;
+  MaterialProperty<RankTwoTensor> & _small_strain;
+  MaterialProperty<RankTwoTensor> & _grad_enrich_disp_tensor;
+  const MaterialProperty<RankTwoTensor> & _grad_disp_tensor_old;
+  const MaterialProperty<RankTwoTensor> & _small_strain_old;
+  const MaterialProperty<RankTwoTensor> & _grad_enrich_disp_tensor_old;
+
+private:
+  /// enrichment function value
+  std::vector<Real> _B;
+  /// derivatives of enrichment function respect to global cooridnate
+  std::vector<RealVectorValue> _dBX;
+  /// derivatives of enrichment function respect to crack front cooridnate
+  std::vector<RealVectorValue> _dBx;
+  /// enrichment function at node I
+  std::vector<std::vector<Real>> _BI;
+  /// shape function
+  const std::vector<std::vector<Real>> * _fe_phi;
+  /// gradient of shape function
+  const std::vector<std::vector<RealGradient>> * _fe_dphi;
+  NonlinearSystem * _nl;
+  const NumericVector<Number> * _sln;
+};

modules/xfem/include/userobjects/GeometricCutUserObject.h

@@ -195,6 +195,8 @@ class GeometricCutUserObject : public CrackFrontPointsProvider
    */
   CutSubdomainID getCutSubdomainID(const Elem * elem) const;
 
+  virtual bool isCutterMeshChanged() const;
+
 protected:
   /// Pointer to the XFEM controller object
   std::shared_ptr<XFEM> _xfem;

modules/xfem/include/userobjects/MeshCut2DFractureUserObject.h

@@ -27,6 +27,8 @@ class MeshCut2DFractureUserObject : public MeshCut2DUserObjectBase
 
   virtual void initialize() override;
 
+  virtual bool isCutterMeshChanged() const override;
+
 protected:
   virtual void findActiveBoundaryGrowth() override;
 

modules/xfem/include/userobjects/MeshCut2DUserObjectBase.h

@@ -49,6 +49,8 @@ class MeshCut2DUserObjectBase : public MeshCutUserObjectBase
   virtual const std::vector<RealVectorValue>
   getCrackPlaneNormals(unsigned int num_crack_front_points) const override;
 
+  bool isPointOnEdgeBoundary(const Point & point, Real tolerance = 1e-10);
+
 protected:
   /// The FE solution mesh
   MooseMesh & _mesh;

modules/xfem/src/base/XFEM.C

@@ -183,6 +183,19 @@ XFEM::clearGeomMarkedElems()
   _geom_marked_elems_3d.clear();
 }
 
+bool
+XFEM::didNearTipEnrichmentChange()
+{
+  bool cutter_mesh_changed = false;
+  for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
+  {
+    cutter_mesh_changed = _geometric_cuts[i]->isCutterMeshChanged();
+    if (cutter_mesh_changed)
+      break;
+  }
+  return cutter_mesh_changed;
+}
+
 void
 XFEM::storeCrackTipOriginAndDirection()
 {
@@ -614,7 +627,6 @@ XFEM::markCutEdgesByState(Real time)
     // crack tip origin coordinates and direction
     Point crack_tip_origin(0, 0, 0);
     Point crack_tip_direction(0, 0, 0);
-
     if (isElemAtCrackTip(elem)) // crack tip element's crack intiation
     {
       orig_cut_side_id = CEMElem->getTipEdgeID();