Skip to content

Commit

Permalink
Random fields spallation
Browse files Browse the repository at this point in the history
ref #98
  • Loading branch information
hugary1995 committed Jul 29, 2021
1 parent 515ecaf commit 0f2bdd3
Show file tree
Hide file tree
Showing 38 changed files with 260,564 additions and 1,250 deletions.
2 changes: 2 additions & 0 deletions include/materials/PFFComputeMultipleInelasticStress.h
Original file line number Diff line number Diff line change
Expand Up @@ -24,11 +24,13 @@ class PFFComputeMultipleInelasticStress : public ComputeMultipleInelasticStress

// @{ Strain energy density and its derivative w/r/t damage
const MaterialPropertyName _psie_name;
MaterialProperty<Real> & _psie;
MaterialProperty<Real> & _psie_active;
// @}

// @{ Interface energy density and its derivative w/r/t damage
const MaterialPropertyName _psii_name;
MaterialProperty<Real> & _psii;
MaterialProperty<Real> & _psii_active;
const MaterialProperty<Real> & _psii_active_old;
// @}
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
//* 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 "ADRankTwoTensorForward.h"
#include "ADSingleVariableReturnMappingSolution.h"
#include "BaseNameInterface.h"
#include "DerivativeMaterialPropertyNameInterface.h"

class ComputeStressWithStrengthSurface : public Material,
public ADSingleVariableReturnMappingSolution,
public BaseNameInterface,
public DerivativeMaterialPropertyNameInterface
{
public:
static InputParameters validParams();

ComputeStressWithStrengthSurface(const InputParameters & parameters);

protected:
virtual void initQpStatefulProperties() override;
virtual ADRankTwoTensor computeStress(const ADRankTwoTensor & strain);
virtual ADReal updateTrialState(const ADReal & delta_c);
virtual void computeQpProperties() override;

virtual ADReal computeResidual(const ADReal & effective_trial_stress,
const ADReal & delta_c) override;

virtual ADReal computeDerivative(const ADReal & effective_trial_stress,
const ADReal & delta_c) override;

virtual Real computeReferenceResidual(const ADReal & effective_trial_stress,
const ADReal & delta_c) override;

const ADMaterialProperty<Real> & _K;
const ADMaterialProperty<Real> & _G;

const ADMaterialProperty<RankTwoTensor> & _mechanical_strain;
ADMaterialProperty<RankTwoTensor> & _stress;
ADMaterialProperty<RankTwoTensor> & _elastic_strain;

ADMaterialProperty<Real> & _c;
const MaterialProperty<Real> & _c_old;
ADMaterialProperty<RankTwoTensor> & _structured_strain;
const MaterialProperty<RankTwoTensor> & _structured_strain_old;
ADMaterialProperty<RankTwoTensor> & _N;

const VariableName _d_name;
const MaterialPropertyName _psin_name;
ADMaterialProperty<Real> & _psin_active;
const MaterialPropertyName _psie_name;
ADMaterialProperty<Real> & _psie;
ADMaterialProperty<Real> & _psie_active;
const MaterialPropertyName _g_name;
const ADMaterialProperty<Real> & _g;
const ADMaterialProperty<Real> & _dg_dd;

const Real _gamma_1;
const Real _gamma_0;
};
14 changes: 10 additions & 4 deletions src/materials/PFFComputeMultipleInelasticStress.C
Original file line number Diff line number Diff line change
Expand Up @@ -38,10 +38,12 @@ PFFComputeMultipleInelasticStress::PFFComputeMultipleInelasticStress(
: ComputeMultipleInelasticStress(parameters),
// The strain energy density
_psie_name(_base_name + getParam<MaterialPropertyName>("strain_energy_density")),
_psie(declareProperty<Real>(_psie_name)),
_psie_active(declareProperty<Real>(_psie_name + "_active")),

// The interface energy density
_psii_name(_base_name + getParam<MaterialPropertyName>("interface_energy_density")),
_psii(declareProperty<Real>(_psii_name)),
_psii_active(declareProperty<Real>(_psii_name + "_active")),
_psii_active_old(getMaterialPropertyOldByName<Real>(_psii_name + "_active")),

Expand Down Expand Up @@ -124,7 +126,7 @@ PFFComputeMultipleInelasticStress::computeQpStress()
Real residual0 = residual;
int its = 0;
if (residual > 0)
while (abs(residual) > 1e-11 && abs(residual) > 1e-8 * residual0)
while (abs(residual) > 1e-8 && abs(residual) > 1e-6 * residual0)
{
Real jacobian = -g(_c[_qp], 2) * _psii_active_old[_qp];
Real step = -residual / jacobian;
Expand Down Expand Up @@ -172,15 +174,19 @@ PFFComputeMultipleInelasticStress::computeQpStress()
RankTwoTensor strain_op_active = strain_op_tr > 0 ? strain_op : strain_op_dev;
RankTwoTensor strain_op_inactive = strain_op - strain_op_active;

// Compute stress
_stress[_qp] = _elasticity_tensor[_qp] * (_gip[_qp] * strain_ip_active + strain_ip_inactive +
g(_c[_qp], 0) * strain_op_active + strain_op_inactive);

// Compute driving energy
_psie[_qp] = _stress[_qp].doubleContraction(strain_ip) / 2;
_psie_active[_qp] =
strain_ip_active.doubleContraction(_elasticity_tensor[_qp] * strain_ip_active);
_psii[_qp] = _stress[_qp].doubleContraction(strain_op) / 2;
_psii_active[_qp] =
strain_op_active.doubleContraction(_elasticity_tensor[_qp] * strain_op_active);

// Compute stress
_stress[_qp] = _elasticity_tensor[_qp] * (_gip[_qp] * strain_ip_active + strain_ip_inactive +
g(_c[_qp], 0) * strain_op_active + strain_op_inactive);
// Transform stress back to the global coordinates
_stress[_qp] = Q.transpose() * _stress[_qp] * Q;

// Update jacobian
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,148 @@
//* This file is part of the RACCOON application
//* being developed at Dolbow lab at Duke University
//* http://dolbow.pratt.duke.edu

#include "ComputeStressWithStrengthSurface.h"

registerMooseObject("raccoonApp", ComputeStressWithStrengthSurface);

InputParameters
ComputeStressWithStrengthSurface::validParams()
{
InputParameters params = Material::validParams();
params += ADSingleVariableReturnMappingSolution::validParams();
params += BaseNameInterface::validParams();

params.addClassDescription("Compute small deformation stress with a strength surface.");
params.addRequiredParam<MaterialPropertyName>("bulk_modulus", "The bulk modulus $K$");
params.addRequiredParam<MaterialPropertyName>("shear_modulus", "The shear modulus $G$");
params.addParam<MaterialPropertyName>(
"structure_energy_density",
"psin",
"Name of the structure energy density computed by this material model");
params.addParam<MaterialPropertyName>(
"strain_energy_density",
"psie",
"Name of the strain energy density computed by this material model");
params.addRequiredCoupledVar("phase_field", "The phase field variable");
params.addParam<MaterialPropertyName>("degradation_function", "g", "The degradation function");
params.addRequiredParam<Real>("gamma_1", "gamma 1");
params.addRequiredParam<Real>("gamma_0", "gamma 0");
return params;
}

ComputeStressWithStrengthSurface::ComputeStressWithStrengthSurface(
const InputParameters & parameters)
: Material(parameters),
ADSingleVariableReturnMappingSolution(parameters),
BaseNameInterface(parameters),
DerivativeMaterialPropertyNameInterface(),
_K(getADMaterialPropertyByName<Real>(prependBaseName("bulk_modulus", true))),
_G(getADMaterialPropertyByName<Real>(prependBaseName("shear_modulus", true))),

_mechanical_strain(getADMaterialProperty<RankTwoTensor>(prependBaseName("mechanical_strain"))),
_stress(declareADProperty<RankTwoTensor>(prependBaseName("stress"))),
_elastic_strain(declareADProperty<RankTwoTensor>(prependBaseName("elastic_strain"))),

_c(declareADProperty<Real>(prependBaseName("effective_structured_strain"))),
_c_old(getMaterialPropertyOldByName<Real>(prependBaseName("effective_structured_strain"))),
_structured_strain(declareADProperty<RankTwoTensor>(prependBaseName("structured_strain"))),
_structured_strain_old(
getMaterialPropertyOldByName<RankTwoTensor>(prependBaseName("structured_strain"))),
_N(declareADProperty<RankTwoTensor>(prependBaseName("flow_direction"))),

_d_name(getVar("phase_field", 0)->name()),
_psin_name(prependBaseName("structure_energy_density", true)),
_psin_active(declareADProperty<Real>(_psin_name + "_active")),
_psie_name(prependBaseName("strain_energy_density", true)),
_psie(declareADProperty<Real>(_psie_name)),
_psie_active(declareADProperty<Real>(_psie_name + "_active")),
_g_name(prependBaseName("degradation_function", true)),
_g(getADMaterialProperty<Real>(_g_name)),
_dg_dd(getADMaterialProperty<Real>(derivativePropertyName(_g_name, {_d_name}))),

_gamma_1(getParam<Real>("gamma_1")),
_gamma_0(getParam<Real>("gamma_0"))
{
}

void
ComputeStressWithStrengthSurface::initQpStatefulProperties()
{
_c[_qp] = 0;
_structured_strain[_qp].zero();
}

ADRankTwoTensor
ComputeStressWithStrengthSurface::computeStress(const ADRankTwoTensor & strain)
{
const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity);
return _g[_qp] * (_K[_qp] * strain.trace() * I2 + 2 * _G[_qp] * strain.deviatoric());
}

ADReal
ComputeStressWithStrengthSurface::updateTrialState(const ADReal & delta_c)
{
const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity);

// Trial strain
ADRankTwoTensor trial_strain = _mechanical_strain[_qp] - _structured_strain_old[_qp];

// Flow direction
ADReal n2 = trial_strain.doubleContraction(trial_strain);
if (MooseUtils::absoluteFuzzyEqual(n2, 0))
n2.value() = libMesh::TOLERANCE * libMesh::TOLERANCE;
ADReal n = std::sqrt(n2);
_N[_qp] = 2 / (2 - delta_c) * (trial_strain.deviatoric() / n / std::sqrt(2) + _gamma_1 * I2);

// Effective trial stress
return computeStress(trial_strain).doubleContraction(_N[_qp]);
}

void
ComputeStressWithStrengthSurface::computeQpProperties()
{
ADReal delta_c = 0;
ADReal effective_trial_stress = updateTrialState(delta_c);

// Return mapping
ADReal phi = computeResidual(effective_trial_stress, delta_c);
if (phi > 0)
returnMappingSolve(effective_trial_stress, delta_c, _console);

// Update internal variables
_c[_qp] = _c_old[_qp] + delta_c;
_structured_strain[_qp] = _structured_strain_old[_qp] + delta_c * _N[_qp];

// Update stress and strain
_elastic_strain[_qp] = _mechanical_strain[_qp] - _structured_strain[_qp];
_stress[_qp] = computeStress(_elastic_strain[_qp]);
_psin_active[_qp] = _gamma_0 * _c[_qp];
_psie[_qp] = _stress[_qp].doubleContraction(_elastic_strain[_qp]) / 2;
_psie_active[_qp] = _psie[_qp] / _g[_qp];
}

Real
ComputeStressWithStrengthSurface::computeReferenceResidual(
const ADReal & /*effective_trial_stress*/, const ADReal & delta_c)
{
ADReal sigma = updateTrialState(delta_c);
return raw_value(sigma - computeStress(delta_c * _N[_qp]).doubleContraction(_N[_qp]));
}

ADReal
ComputeStressWithStrengthSurface::computeResidual(const ADReal & /*effective_trial_stress*/,
const ADReal & delta_c)
{
ADReal sigma = updateTrialState(delta_c);
return sigma - computeStress(delta_c * _N[_qp]).doubleContraction(_N[_qp]) - _g[_qp] * _gamma_0;
}

ADReal
ComputeStressWithStrengthSurface::computeDerivative(const ADReal & /*effective_trial_stress*/,
const ADReal & delta_c)
{
ADReal sigma = updateTrialState(delta_c);
return -computeStress(_N[_qp]).doubleContraction(_N[_qp]) -
computeStress(delta_c * _N[_qp]).doubleContraction(_N[_qp]) / (2 - delta_c);
}
Loading

0 comments on commit 0f2bdd3

Please sign in to comment.