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/GeneralizedExternalDrivingForce.md

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

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,30 @@
+#pragma once
+
+#include "BoundsAuxBase.h"
+
+// Forward Declarations
+class ConditionalBoundsAux;
+
+template <>
+InputParameters validParams<ConditionalBoundsAux>();
+
+/**
+ * Provides a conditional bound of a variable using its old value or fixed
+ * value.
+ */
+class ConditionalBoundsAux : public BoundsAuxBase
+{
+public:
+  static InputParameters validParams();
+
+  ConditionalBoundsAux(const InputParameters & parameters);
+
+protected:
+  virtual Real getBound() override;
+
+  /// The value of fixed bound for the variable
+  Real _fixed_bound_value;
+
+  /// The threshold for conditional bound for the variable
+  Real _threshold_value;
+};

include/materials/GeneralizedExternalDrivingForce.h

@@ -0,0 +1,46 @@
+//* 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"
+
+class GeneralizedExternalDrivingForce : public Material
+{
+public:
+  static InputParameters validParams();
+
+  GeneralizedExternalDrivingForce(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;
+
+  /// energy release rate
+  const ADMaterialProperty<Real> & _Gc;
+
+  /// phase field regularization length
+  const ADMaterialProperty<Real> & _L;
+  /// Lame's first parameter
+  const ADMaterialProperty<Real> & _Lambda;
+  /// shear modulus
+  const ADMaterialProperty<Real> & _mu;
+
+  /// critical tensile strength
+  const Real & _sigma_ts;
+
+  /// critical compressive strength
+  const Real & _sigma_cs;
+
+  /// regularization length dependent parameter
+  const Real & _delta;
+
+  /// name of the degraded stress tensor
+  const MaterialPropertyName _rank_two_tensor; // remove
+  /// the degraded stress tensor
+  const ADMaterialProperty<RankTwoTensor> & _stress;
+};

include/postprocessors/PhaseFieldJIntegral.h

@@ -0,0 +1,36 @@
+//* 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"
+
+class PhaseFieldJIntegral : public SideIntegralPostprocessor
+{
+public:
+  static InputParameters validParams();
+
+  PhaseFieldJIntegral(const InputParameters & parameters);
+
+protected:
+  virtual Real computeQpIntegral() override;
+
+  /// base name of stress
+  const std::string _base_name;
+  /// stress tensor
+  const ADMaterialProperty<RankTwoTensor> & _stress;
+  /// degraded eleastic energy
+  const ADMaterialProperty<Real> & _E_elastic;
+  /// number of displacement variables provided
+  const unsigned int _ndisp;
+  /// du_dx
+  const VariableGradient & _grad_disp_0;
+  /// du_dy
+  const VariableGradient & _grad_disp_1;
+  /// du_dz
+  const VariableGradient & _grad_disp_2;
+  /// direction of J integral
+  const RealVectorValue _t;
+};

src/auxkernels/ConditionalBoundsAux.C

@@ -0,0 +1,41 @@
+#include "ConditionalBoundsAux.h"
+
+registerMooseObject("raccoonApp", ConditionalBoundsAux);
+
+defineLegacyParams(ConditionalBoundsAux);
+
+InputParameters
+ConditionalBoundsAux::validParams()
+{
+  InputParameters params = BoundsAuxBase::validParams();
+  params.addClassDescription(
+      "Provides the bound of the variable a lower bound according to variable "
+      "value"
+      "for kumar's phase-field fracture model 2020"
+      "variable to PETSc's SNES variational inequalities solver.");
+  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");
+  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/GeneralizedExternalDrivingForce.C

@@ -0,0 +1,88 @@
+//* This file is part of the RACCOON application
+//* being developed at Dolbow lab at Duke University
+//* http://dolbow.pratt.duke.edu
+///
+// Aditya Kumar, Blaise Bourdin, Gilles A. Francfort, Oscar Lopez-Pamies,
+// Revisiting nucleation in the phase-field approach to brittle fracture,
+// Journal of the Mechanics and Physics of Solids,
+// Volume 142,
+// 2020,
+// 104027,
+// ISSN 0022-5096,
+// https://doi.org/10.1016/j.jmps.2020.104027.
+// (https://www.sciencedirect.com/science/article/pii/S0022509620302623)
+
+#include "Function.h"
+#include "GeneralizedExternalDrivingForce.h"
+
+registerADMooseObject("raccoonApp", GeneralizedExternalDrivingForce);
+
+InputParameters GeneralizedExternalDrivingForce::validParams() {
+  InputParameters params = Material::validParams();
+  params.addClassDescription("computes the Kumar coeffs"
+                             "Ce()");
+  params.addParam<MaterialPropertyName>(
+      "energy_release_rate", "Gc", "energy release rate or fracture toughness");
+  params.addParam<MaterialPropertyName>(
+      "phase_field_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",
+                                "critical tensile strength");
+  params.addRequiredParam<Real>("compressive_strength",
+                                "critical compressive strength");
+  params.addRequiredParam<Real>("delta", "delta");
+  params.addParam<MaterialPropertyName>("rank_two_tensor", "stress",
+                                        "name of the stress tensor");
+  params.addParam<MaterialPropertyName>(
+      "external_driving_force_name", "ex_driving",
+      "name of the material that holds the external_driving_force");
+  return params;
+}
+
+GeneralizedExternalDrivingForce::GeneralizedExternalDrivingForce(
+    const InputParameters &parameters)
+    : Material(parameters), _ex_driving_name(getParam<MaterialPropertyName>(
+                                "external_driving_force_name")),
+      _ex_driving(declareADProperty<Real>(_ex_driving_name)),
+      _Gc(getADMaterialProperty<Real>("energy_release_rate")),
+      _L(getADMaterialProperty<Real>("phase_field_regularization_length")),
+      _Lambda(getADMaterialProperty<Real>("Lambda")),
+      _mu(getADMaterialProperty<Real>("shear_modulus")),
+      _sigma_ts(getParam<Real>("tensile_strength")),
+      _sigma_cs(getParam<Real>("compressive_strength")),
+      _delta(getParam<Real>("delta")),
+      _rank_two_tensor(getParam<MaterialPropertyName>("rank_two_tensor")),
+      _stress(getADMaterialProperty<RankTwoTensor>(_rank_two_tensor)) {}
+
+void GeneralizedExternalDrivingForce::computeQpProperties() {
+  ADReal K = _Lambda[_qp] + 2.0 * _mu[_qp] / 3.0;
+  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;
+  ADReal _temp = _Gc[_qp] * 3.0 / _L[_qp] / 8.0;
+  ADReal I1 = _stress[_qp](0, 0) + _stress[_qp](1, 1) + _stress[_qp](2, 2);
+  ADReal J2 = (pow(_stress[_qp](0, 0) - _stress[_qp](1, 1), 2) +
+               pow(_stress[_qp](1, 1) - _stress[_qp](2, 2), 2) +
+               pow(_stress[_qp](0, 0) - _stress[_qp](2, 2), 2)) /
+                  6.0 +
+              pow(_stress[_qp](0, 1), 2) + pow(_stress[_qp](0, 2), 2) +
+              pow(_stress[_qp](2, 1), 2);
+
+  mooseAssert(J2 >= 0, "Negative J2");
+
+  ADReal _beta_0 = _delta * _temp;
+  ADReal _beta_1 = (-_gamma_1 * _temp - _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 * _temp + _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,49 @@
+//* 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.addClassDescription("Compute the J integral for a phase-field model of fracture");
+  params.addParam<std::string>("base_name",
+                               "Optional parameter that allows the user to define "
+                               "multiple mechanics material systems on the same "
+                               "block, i.e. for multiple phases");
+  params.addRequiredParam<RealVectorValue>("J_direction", "direction of J integral");
+  params.addRequiredParam<MaterialPropertyName>("elastic_energy_name",
+                                                "name of the elastic energy");
+  params.addRequiredCoupledVar(
+      "displacements",
+      "The displacements appropriate for the simulation geometry and coordinate system");
+  return params;
+}
+
+PhaseFieldJIntegral::PhaseFieldJIntegral(const InputParameters & parameters)
+  : SideIntegralPostprocessor(parameters),
+    _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
+    _stress(getADMaterialPropertyByName<RankTwoTensor>(_base_name + "stress")),
+    _E_elastic(getADMaterialProperty<Real>("elastic_energy_name")),
+    _ndisp(coupledComponents("displacements")),
+    _grad_disp_0(coupledGradient("displacements", 0)),
+    _grad_disp_1(_ndisp >= 2 ? coupledGradient("displacements", 1) : _grad_zero),
+    _grad_disp_2(_ndisp >= 3 ? coupledGradient("displacements", 2) : _grad_zero),
+    _t(getParam<RealVectorValue>("J_direction"))
+{
+}
+
+Real
+PhaseFieldJIntegral::computeQpIntegral()
+{
+  RankTwoTensor H(_grad_disp_0[_qp], _grad_disp_1[_qp], _grad_disp_2[_qp]);
+  RankTwoTensor I2(RankTwoTensor::initIdentity);
+  ADRankTwoTensor Sigma = _E_elastic[_qp] * I2 - H.transpose() * _stress[_qp];
+  RealVectorValue n = _normals[_qp];
+  ADReal value = _t * Sigma * n;
+  return value.value();
+}