Jint surfing

.gitignore

@@ -239,3 +239,9 @@ python/peacock/tests/postprocessor_tab/TestPostprocessorPluginManager_test_scrip
 !python/peacock/tests/**/input/*.*
 peacock_tmp_diff.exo
 *.e.diff
+
+# tmp file when computing
+.nfs*
+
+#log file from cluster computing moose
+*_1???????

doc/content/source/auxkernels/ConditionalBoundsAux.md

@@ -0,0 +1,11 @@
+# ConditionalBoundsAux
+
+!syntax description /Bounds/ConditionalBoundsAux
+
+## Example Input File Syntax
+
+!syntax parameters /Bounds/ConditionalBoundsAux
+
+!syntax inputs /Bounds/ConditionalBoundsAux
+
+!syntax children /Bounds/ConditionalBoundsAux

doc/content/source/materials/NucleationMicroForce.md

@@ -0,0 +1,11 @@
+# NucleationMicroForce
+
+!syntax description /Materials/NucleationMicroForce
+
+## Example Input File Syntax
+
+!syntax parameters /Materials/NucleationMicroForce
+
+!syntax inputs /Materials/NucleationMicroForce
+
+!syntax children /Materials/NucleationMicroForce

doc/content/source/postprocessors/PhaseFieldJIntegral.md

@@ -0,0 +1,11 @@
+# PhaseFieldJIntegral
+
+!syntax description /UserObjects/PhaseFieldJIntegral
+
+## Example Input File Syntax
+
+!syntax parameters /UserObjects/PhaseFieldJIntegral
+
+!syntax inputs /UserObjects/PhaseFieldJIntegral
+
+!syntax children /UserObjects/PhaseFieldJIntegral

doc/content/tutorials/12_surfing_boundary_problem.md

@@ -0,0 +1 @@
+

doc/menu.yml

@@ -16,6 +16,7 @@ Tutorials:
   9. Three-point bending: tutorials/09_three_point_bending.md
   10. Quenching of bi-beam: tutorials/10_quenching_bibeam.md
   11. Adaptive mesh refinement: tutorials/11_adaptivity.md
+  12. Surfing Boundary Problem: tutorials/12_surfing_boundary_problem.md
 Theory:
   Introduction:   theory/intro.md
   State variables and kinematics: theory/kinematics.md

include/auxkernels/ConditionalBoundsAux.h

@@ -0,0 +1,24 @@
+#pragma once
+
+#include "BoundsAuxBase.h"
+
+/**
+ * Provides a conditional bound of a variable using its old value or a fixed
+ * value.
+ */
+class ConditionalBoundsAux : public BoundsAuxBase
+{
+public:
+  static InputParameters validParams();
+
+  ConditionalBoundsAux(const InputParameters & parameters);
+
+protected:
+  virtual Real getBound() override;
+
+  /// The value of the fixed bound for the variable
+  Real _fixed_bound_value;
+
+  /// The threshold for conditional bound for the variable
+  Real _threshold_value;
+};

include/materials/NucleationMicroForce.h

@@ -0,0 +1,55 @@
+//* This file is part of the RACCOON application
+//* being developed at Dolbow lab at Duke University
+//* http://dolbow.pratt.duke.edu
+
+#pragma once
+
+#include "Material.h"
+#include "BaseNameInterface.h"
+
+/**
+ * The class implements the external driving force to recover a Drucker-Prager
+ * strength envelope. See Kumar et. al. https://doi.org/10.1016/j.jmps.2020.104027.
+ */
+class NucleationMicroForce : public Material, public BaseNameInterface
+{
+public:
+  static InputParameters validParams();
+
+  NucleationMicroForce(const InputParameters & parameters);
+
+protected:
+  virtual void computeQpProperties() override;
+
+  /// Name of the external driving force
+  const MaterialPropertyName _ex_driving_name;
+
+  /// The external driving force
+  ADMaterialProperty<Real> & _ex_driving;
+
+  ///@{ Phase field properties
+  /// The fracture toughness
+  const ADMaterialProperty<Real> & _Gc;
+  /// The normalization constant
+  const ADMaterialProperty<Real> & _c0;
+  /// phase field regularization length
+  const ADMaterialProperty<Real> & _L;
+  ///@}
+
+  /// Lame's first parameter
+  const ADMaterialProperty<Real> & _lambda;
+  /// The shear modulus
+  const ADMaterialProperty<Real> & _mu;
+
+  /// The critical tensile strength
+  const Real & _sigma_ts;
+
+  /// The critical compressive strength
+  const Real & _sigma_cs;
+
+  /// The regularization length dependent parameter
+  const Real & _delta;
+
+  /// The stress tensor
+  const ADMaterialProperty<RankTwoTensor> & _stress;
+};

include/postprocessors/PhaseFieldJIntegral.h

@@ -0,0 +1,31 @@
+//* This file is part of the RACCOON application
+//* being developed at Dolbow lab at Duke University
+//* http://dolbow.pratt.duke.edu
+
+#pragma once
+
+#include "SideIntegralPostprocessor.h"
+#include "RankTwoTensor.h"
+#include "BaseNameInterface.h"
+
+class PhaseFieldJIntegral : public SideIntegralPostprocessor, public BaseNameInterface
+{
+public:
+  static InputParameters validParams();
+
+  PhaseFieldJIntegral(const InputParameters & parameters);
+
+protected:
+  virtual Real computeQpIntegral() override;
+
+  /// The stress tensor
+  const ADMaterialProperty<RankTwoTensor> & _stress;
+  /// The strain energy density
+  const ADMaterialProperty<Real> & _psie;
+  /// Number of displacement variables provided
+  const unsigned int _ndisp;
+  /// Gradient of displacements
+  std::vector<const VariableGradient *> _grad_disp;
+  /// Direction of J integral
+  const RealVectorValue _t;
+};

src/auxkernels/ConditionalBoundsAux.C

@@ -0,0 +1,36 @@
+#include "ConditionalBoundsAux.h"
+
+registerMooseObject("raccoonApp", ConditionalBoundsAux);
+
+InputParameters
+ConditionalBoundsAux::validParams()
+{
+  InputParameters params = BoundsAuxBase::validParams();
+  params.addClassDescription(
+      "This class conditionally enforces a lower bound. When the variable value is below a given "
+      "threshold, a constant value is used as the bound; when the variable value is above a given "
+      "threshold, irreversibility is enforced.");
+  params.addRequiredParam<Real>("fixed_bound_value", "The value of fixed bound for the variable");
+  params.addRequiredParam<Real>("threshold_value",
+                                "The threshold for conditional history bound for the variable");
+  params.set<MooseEnum>("bound_type") = "lower";
+  params.suppressParameter<MooseEnum>("bound_type");
+  return params;
+}
+
+ConditionalBoundsAux::ConditionalBoundsAux(const InputParameters & parameters)
+  : BoundsAuxBase(parameters),
+    _fixed_bound_value(getParam<Real>("fixed_bound_value")),
+    _threshold_value(getParam<Real>("threshold_value"))
+{
+}
+
+Real
+ConditionalBoundsAux::getBound()
+{
+  Real d_old = _var.getNodalValueOld(*_current_node);
+  if (d_old >= _threshold_value)
+    return d_old;
+  else
+    return _fixed_bound_value;
+}

src/materials/NucleationMicroForce.C

@@ -0,0 +1,92 @@
+//* This file is part of the RACCOON application
+//* being developed at Dolbow lab at Duke University
+//* http://dolbow.pratt.duke.edu
+
+#include "Function.h"
+#include "NucleationMicroForce.h"
+
+registerADMooseObject("raccoonApp", NucleationMicroForce);
+
+InputParameters
+NucleationMicroForce::validParams()
+{
+  InputParameters params = Material::validParams();
+  params += BaseNameInterface::validParams();
+
+  params.addClassDescription(
+      "This class computes the external driving force for nucleation given "
+      "a Drucker-Prager strength envelope. The implementation follows Kumar et. al. (2020).");
+
+  params.addParam<MaterialPropertyName>(
+      "fracture_toughness", "Gc", "energy release rate or fracture toughness");
+  params.addParam<MaterialPropertyName>(
+      "normalization_constant", "c0", "The normalization constant $c_0$");
+  params.addParam<MaterialPropertyName>(
+      "regularization_length", "l", "the phase field regularization length");
+
+  params.addParam<MaterialPropertyName>("lambda", "lambda", "Lame's first parameter lambda");
+  params.addParam<MaterialPropertyName>("shear_modulus", "G", "shear modulus mu or G");
+
+  params.addRequiredParam<Real>(
+      "tensile_strength", "The tensile strength of the material beyond which the material fails.");
+  params.addRequiredParam<Real>(
+      "compressive_strength",
+      "The compressive strength of the material beyond which the material fails.");
+
+  params.addRequiredParam<Real>("delta", "delta");
+  params.addParam<MaterialPropertyName>(
+      "external_driving_force_name",
+      "ex_driving",
+      "Name of the material that holds the external_driving_force");
+
+  return params;
+}
+
+NucleationMicroForce::NucleationMicroForce(const InputParameters & parameters)
+  : Material(parameters),
+    BaseNameInterface(parameters),
+    _ex_driving(declareADProperty<Real>(prependBaseName("external_driving_force_name", true))),
+    _Gc(getADMaterialProperty<Real>(prependBaseName("fracture_toughness", true))),
+    _c0(getADMaterialProperty<Real>(prependBaseName("normalization_constant", true))),
+    _L(getADMaterialProperty<Real>(prependBaseName("regularization_length", true))),
+    _lambda(getADMaterialProperty<Real>(prependBaseName("lambda", true))),
+    _mu(getADMaterialProperty<Real>(prependBaseName("shear_modulus", true))),
+    _sigma_ts(getParam<Real>("tensile_strength")),
+    _sigma_cs(getParam<Real>("compressive_strength")),
+    _delta(getParam<Real>("delta")),
+    _stress(getADMaterialProperty<RankTwoTensor>(prependBaseName("stress")))
+{
+}
+
+void
+NucleationMicroForce::computeQpProperties()
+{
+  // The bulk modulus
+  ADReal K = _lambda[_qp] + 2 * _mu[_qp] / 3;
+
+  // Parameters in the strength surface
+  ADReal gamma_0 =
+      (_mu[_qp] + 3 * K) * _sigma_ts * _L[_qp] / _Gc[_qp] / 18 / _mu[_qp] / _mu[_qp] / K / K;
+  ADReal gamma_1 = (1.0 + _delta) / (2.0 * _sigma_ts * _sigma_cs);
+  ADReal gamma_2 = (8 * _mu[_qp] + 24 * K - 27 * _sigma_ts) / 144 / _mu[_qp] / K;
+
+  // The mobility
+  ADReal M = _Gc[_qp] / _L[_qp] / _c0[_qp];
+
+  // Invariants of the stress
+  ADReal I1 = _stress[_qp].trace();
+  ADRankTwoTensor stress_dev = _stress[_qp].deviatoric();
+  ADReal J2 = stress_dev.doubleContraction(stress_dev);
+
+  // Just to be extra careful... J2 is for sure non-negative.
+  mooseAssert(J2 >= 0, "Negative J2");
+
+  // Compute the external driving force required to recover the desired strength envelope.
+  ADReal beta_0 = _delta * M;
+  ADReal beta_1 = (-gamma_1 * M - gamma_2) * (_sigma_cs - _sigma_ts) -
+                  gamma_0 * (pow(_sigma_cs, 3) - pow(_sigma_ts, 3));
+  ADReal beta_2 = std::sqrt(3.0) * ((-gamma_1 * M + gamma_2) * (_sigma_cs + _sigma_ts) +
+                                    gamma_0 * (pow(_sigma_cs, 3) + pow(_sigma_ts, 3)));
+  ADReal beta_3 = _L[_qp] * _sigma_ts / _mu[_qp] / K / _Gc[_qp];
+  _ex_driving[_qp] = (beta_2 * std::sqrt(J2) + beta_1 * I1 + beta_0) / (1 + beta_3 * I1 * I1);
+}

src/postprocessors/PhaseFieldJIntegral.C

@@ -0,0 +1,47 @@
+//* This file is part of the RACCOON application
+//* being developed at Dolbow lab at Duke University
+//* http://dolbow.pratt.duke.edu
+
+#include "PhaseFieldJIntegral.h"
+
+registerMooseObject("raccoonApp", PhaseFieldJIntegral);
+
+InputParameters
+PhaseFieldJIntegral::validParams()
+{
+  InputParameters params = SideIntegralPostprocessor::validParams();
+  params += BaseNameInterface::validParams();
+  params.addClassDescription("Compute the J integral for a phase-field model of fracture");
+  params.addRequiredParam<RealVectorValue>("J_direction", "direction of J integral");
+  params.addParam<MaterialPropertyName>("strain_energy_density",
+                                        "psie"
+                                        "Name of the strain energy density");
+  params.addRequiredCoupledVar(
+      "displacements",
+      "The displacements appropriate for the simulation geometry and coordinate system");
+  return params;
+}
+
+PhaseFieldJIntegral::PhaseFieldJIntegral(const InputParameters & parameters)
+  : SideIntegralPostprocessor(parameters),
+    BaseNameInterface(parameters),
+    _stress(getADMaterialPropertyByName<RankTwoTensor>(prependBaseName("stress"))),
+    _psie(getADMaterialProperty<Real>(prependBaseName("strain_energy_density"))),
+    _ndisp(coupledComponents("displacements")),
+    _grad_disp(coupledGradients("displacements")),
+    _t(getParam<RealVectorValue>("J_direction"))
+{
+  // set unused dimensions to zero
+  for (unsigned i = _ndisp; i < 3; ++i)
+    _grad_disp.push_back(&_grad_zero);
+}
+
+Real
+PhaseFieldJIntegral::computeQpIntegral()
+{
+  RankTwoTensor H((*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]);
+  RankTwoTensor I2(RankTwoTensor::initIdentity);
+  ADRankTwoTensor Sigma = _psie[_qp] * I2 - H.transpose() * _stress[_qp];
+  RealVectorValue n = _normals[_qp];
+  return raw_value(_t * Sigma * n);
+}