Commit 817cdaaf authored by Antoine Jacquey's avatar Antoine Jacquey

Implemented damage in two ways: fully coupled or via a MultiApp.

parent 6c336965
// /******************************************************************************/
// /* LYNX, a MOOSE-based application */
// /* */
// /* Copyright (C) 2017 by Antoine B. Jacquey and Mauro Cacace */
// /* GFZ Potsdam, German Research Centre for Geosciences */
// /* */
// /* This program is free software: you can redistribute it and/or modify */
// /* it under the terms of the GNU General Public License as published by */
// /* the Free Software Foundation, either version 3 of the License, or */
// /* (at your option) any later version. */
// /* */
// /* This program is distributed in the hope that it will be useful, */
// /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
// /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
// /* GNU General Public License for more details. */
// /* */
// /* You should have received a copy of the GNU General Public License */
// /* along with this program. If not, see <http://www.gnu.org/licenses/> */
// /******************************************************************************/
// #ifndef LYNXSTRAINRATIOAUX_H
// #define LYNXSTRAINRATIOAUX_H
// #include "LynxStrainAuxBase.h"
// class LynxStrainRatioAux;
// template <>
// InputParameters validParams<LynxStrainRatioAux>();
// class LynxStrainRatioAux : public LynxStrainAuxBase
// {
// public:
// LynxStrainRatioAux(const InputParameters & parameters);
// protected:
// virtual Real computeValue() override;
// };
// #endif // LYNXSTRAINRATIOAUX_H
\ No newline at end of file
/******************************************************************************/
/* LYNX, a MOOSE-based application */
/* */
/* Copyright (C) 2017 by Antoine B. Jacquey and Mauro Cacace */
/* GFZ Potsdam, German Research Centre for Geosciences */
/* */
/* This program is free software: you can redistribute it and/or modify */
/* it under the terms of the GNU General Public License as published by */
/* the Free Software Foundation, either version 3 of the License, or */
/* (at your option) any later version. */
/* */
/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
/* GNU General Public License for more details. */
/* */
/* You should have received a copy of the GNU General Public License */
/* along with this program. If not, see <http://www.gnu.org/licenses/> */
/******************************************************************************/
#ifndef LYNXSTRAINRATIOAUX_H
#define LYNXSTRAINRATIOAUX_H
#include "AuxKernel.h"
#include "RankTwoTensor.h"
#include "DerivativeMaterialInterface.h"
class LynxStrainRatioAux;
template <>
InputParameters validParams<LynxStrainRatioAux>();
class LynxStrainRatioAux : public DerivativeMaterialInterface<AuxKernel>
{
public:
LynxStrainRatioAux(const InputParameters & parameters);
protected:
virtual Real computeValue();
const MaterialProperty<RankTwoTensor> & _elastic_strain;
};
#endif // LYNXSTRAINRATIOAUX_H
/******************************************************************************/
/* LYNX, a MOOSE-based application */
/* */
/* Copyright (C) 2017 by Antoine B. Jacquey and Mauro Cacace */
/* GFZ Potsdam, German Research Centre for Geosciences */
/* */
/* This program is free software: you can redistribute it and/or modify */
/* it under the terms of the GNU General Public License as published by */
/* the Free Software Foundation, either version 3 of the License, or */
/* (at your option) any later version. */
/* */
/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
/* GNU General Public License for more details. */
/* */
/* You should have received a copy of the GNU General Public License */
/* along with this program. If not, see <http://www.gnu.org/licenses/> */
/******************************************************************************/
#ifndef LYNXDAMAGERATE_H
#define LYNXDAMAGERATE_H
#include "Kernel.h"
class LynxDamageRate;
template <>
InputParameters validParams<LynxDamageRate>();
class LynxDamageRate : public Kernel
{
public:
LynxDamageRate(const InputParameters & parameters);
protected:
virtual Real computeQpResidual() override;
virtual Real computeQpJacobian() override;
virtual Real computeQpOffDiagJacobian(unsigned int jvar) override;
const VariableValue & _u_old;
bool _coupled_dam;
const VariableValue & _damage_rate;
};
#endif // LYNXDAMAGERATE_H
......@@ -40,6 +40,8 @@ protected:
virtual void elasticModuli() override;
virtual void plasticCorrection(Real & pressure, RankTwoTensor & stress_dev) override;
virtual void damageCorrection() override;
virtual RankFourTensor damageTangentOperator(const RankTwoTensor & flow_direction,
const RankFourTensor & tme) override;
virtual Real computePlasticityYield(const Real & pressure, const Real & eqv_stress);
virtual Real plasticIncrement(const Real & /*pressure*/, const Real & eqv_stress);
virtual Real computeConvexPlasticityYield2(const Real & pressure, const Real & eqv_stress);
......@@ -69,6 +71,7 @@ protected:
// Plastic parameters
const std::vector<Real> _friction_angle;
const std::vector<Real> _cohesion;
const std::vector<Real> _porous_coeff;
const std::vector<Real> _porous_coeff_linear;
std::vector<Real> _one_on_plastic_eta;
......@@ -77,9 +80,21 @@ protected:
// Damage-Plasticity structure
damage_plasticity * _damage_plasticity;
// Damage-Plasticity utilities
Real _dyield_dp_tr;
Real _dyield_dq_tr;
RankTwoTensor _damaged_stress;
RankFourTensor _damaged_tensor;
// Strain properties
MaterialProperty<RankTwoTensor> & _elastic_strain;
const MaterialProperty<RankTwoTensor> & _elastic_strain_old;
// Stress properties
// MaterialProperty<RankTwoTensor> & _dstress_ddamage;
// Damage properties
MaterialProperty<Real> & _damage_rate;
};
#endif // LYNXDAMAGEDEFORMATION_H
/******************************************************************************/
/* LYNX, a MOOSE-based application */
/* */
/* Copyright (C) 2017 by Antoine B. Jacquey and Mauro Cacace */
/* GFZ Potsdam, German Research Centre for Geosciences */
/* */
/* This program is free software: you can redistribute it and/or modify */
/* it under the terms of the GNU General Public License as published by */
/* the Free Software Foundation, either version 3 of the License, or */
/* (at your option) any later version. */
/* */
/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
/* GNU General Public License for more details. */
/* */
/* You should have received a copy of the GNU General Public License */
/* along with this program. If not, see <http://www.gnu.org/licenses/> */
/******************************************************************************/
#ifndef LYNXDAMAGEDEFORMATIONNEW_H
#define LYNXDAMAGEDEFORMATIONNEW_H
#include "LynxDeformation.h"
class LynxDamageDeformationNew;
template <>
InputParameters validParams<LynxDamageDeformationNew>();
class LynxDamageDeformationNew : public LynxDeformationBase
{
public:
LynxDamageDeformationNew(const InputParameters & parameters);
virtual ~LynxDamageDeformationNew() {}
protected:
virtual void initQpStatefulProperties() override;
virtual void initializeQpDeformation() override;
// virtual void elasticModuli() override;
virtual void plasticCorrection(Real & pressure, RankTwoTensor & stress_dev) override;
virtual void damageCorrection() override;
virtual RankFourTensor damageTangentOperator(const RankTwoTensor & flow_direction,
const RankFourTensor & tme) override;
virtual Real computePlasticityYield(const Real & pressure, const Real & eqv_stress);
virtual Real plasticIncrement(const Real & /*pressure*/, const Real & eqv_stress);
virtual Real computeConvexPlasticityYield2(const Real & pressure, const Real & eqv_stress);
virtual Real computeConvexPlasticityYield2(const Real & rho);
virtual Real convexPlasticIncrement(Real & vol_plastic_incr, Real & eqv_plastic_incr);
virtual void computeDamageProperties(const Real & pressure, const Real & eqv_stress);
virtual void updateDamageParameters();
virtual void initializeDamageParameters();
virtual void updateDamageConvexParameters(const Real & pressure, const Real & eqv_stress);
virtual Real convexReferencePressure();
virtual Real dConvexPlasticYield2(const Real & rho);
virtual Real dConvexPlasticYield2_dp(const Real & pressure, const Real & eqv_stress);
virtual Real dConvexPlasticYield2_dq(const Real & pressure, const Real & eqv_stress);
virtual Real getConvexProjection(const Real & x1, const Real & x2);
virtual Real strainRatio(const RankTwoTensor & elastic_strain);
virtual RankTwoTensor rotatedElasticStrain(const RankTwoTensor & elastic_strain);
// Coupled variables
bool _coupled_dam;
const VariableValue & _damage;
const VariableValue & _damage_old;
bool _coupled_phi;
const VariableValue & _porosity;
// Elastic moduli parameters
const std::vector<Real> _damage_modulus;
// Plastic parameters
const std::vector<Real> _friction_angle;
const std::vector<Real> _cohesion;
const std::vector<Real> _porous_coeff;
const std::vector<Real> _porous_coeff_linear;
std::vector<Real> _one_on_plastic_eta;
std::vector<Real> _one_on_damage_eta;
// Damage-Plasticity structure
damage_plasticity * _damage_plasticity;
// Damage-Plasticity utilities
// Real _dyield_dp_tr;
// Real _dyield_dq_tr;
// RankTwoTensor _damaged_stress;
// RankFourTensor _damaged_tensor;
// Strain properties
MaterialProperty<RankTwoTensor> & _elastic_strain;
const MaterialProperty<RankTwoTensor> & _elastic_strain_old;
// Stress properties
// MaterialProperty<RankTwoTensor> & _dstress_ddamage;
// Damage properties
MaterialProperty<Real> & _damage_rate;
};
#endif // LYNXDAMAGEDEFORMATIONNEW_H
......@@ -57,8 +57,11 @@ protected:
virtual void plasticCorrection(Real & pressure, RankTwoTensor & stress_dev) = 0;
virtual void damageCorrection();
virtual void tangentOperator();
virtual void plasticTangentOperator(const RankTwoTensor & flow_direction,
const RankFourTensor & flow_direction_dyad);
virtual RankFourTensor viscousTangentOperator(const RankFourTensor & flow_direction_dyad);
virtual RankFourTensor plasticTangentOperator(const RankTwoTensor & flow_direction,
const RankFourTensor & flow_direction_dyad);
virtual RankFourTensor damageTangentOperator(const RankTwoTensor & /*flow_direction*/,
const RankFourTensor & /*tme*/);
virtual void finiteTangentOperator();
virtual Real computeStokeEffectiveViscosity(const Real & pressure);
virtual RankTwoTensor computeStokeEffectiveViscosityDerivative();
......
......@@ -208,12 +208,14 @@ struct plasticity
struct damage_plasticity
{
Real _xi0;
Real _gamma;
Real _p_cr;
Real _eta_p; // actually one on eta
Real _eta_d; // actually one on eta
Real _alpha0;
Real _xi0; // critical strain ratio
Real _gamma; // third elastic modulus
Real _p_cr; // critical pressure for capped yield
Real _k0; // coefficient for cohesion
Real _eta_p; // actually one on eta
Real _eta_d; // actually one on eta
Real _alpha0; // coefficient for friction
Real _p_k; // coefficient to fake cohesion for capped yield
Real _p_tr; // trial pressure
Real _q_tr; // trial eqv_stress
Real _p_r; // reference pressure for convex yield
......@@ -230,9 +232,11 @@ struct damage_plasticity
: _xi0(-std::sqrt(3.0)),
_gamma(0.0),
_p_cr(0.0),
_k0(0.0),
_eta_p(0.0),
_eta_d(0.0),
_alpha0(0.0),
_p_k(0.0),
_p_tr(0.0),
_q_tr(0.0),
_p_r(0.0),
......@@ -247,11 +251,17 @@ struct damage_plasticity
_dmu2_dxi_cr(0.0)
{
}
void fill(const Real xi0, const Real gamma, const Real p_cr, const Real eta_p, const Real eta_d)
void fill(const Real xi0,
const Real gamma,
const Real p_cr,
const Real k0,
const Real eta_p,
const Real eta_d)
{
_xi0 = xi0;
_gamma = gamma;
_p_cr = p_cr;
_k0 = k0;
_eta_p = eta_p;
_eta_d = eta_d;
}
......
// /******************************************************************************/
// /* LYNX, a MOOSE-based application */
// /* */
// /* Copyright (C) 2017 by Antoine B. Jacquey and Mauro Cacace */
// /* GFZ Potsdam, German Research Centre for Geosciences */
// /* */
// /* This program is free software: you can redistribute it and/or modify */
// /* it under the terms of the GNU General Public License as published by */
// /* the Free Software Foundation, either version 3 of the License, or */
// /* (at your option) any later version. */
// /* */
// /* This program is distributed in the hope that it will be useful, */
// /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
// /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
// /* GNU General Public License for more details. */
// /* */
// /* You should have received a copy of the GNU General Public License */
// /* along with this program. If not, see <http://www.gnu.org/licenses/> */
// /******************************************************************************/
/******************************************************************************/
/* LYNX, a MOOSE-based application */
/* */
/* Copyright (C) 2017 by Antoine B. Jacquey and Mauro Cacace */
/* GFZ Potsdam, German Research Centre for Geosciences */
/* */
/* This program is free software: you can redistribute it and/or modify */
/* it under the terms of the GNU General Public License as published by */
/* the Free Software Foundation, either version 3 of the License, or */
/* (at your option) any later version. */
/* */
/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
/* GNU General Public License for more details. */
/* */
/* You should have received a copy of the GNU General Public License */
/* along with this program. If not, see <http://www.gnu.org/licenses/> */
/******************************************************************************/
// #include "LynxStrainRatioAux.h"
#include "LynxStrainRatioAux.h"
// registerMooseObject("LynxApp", LynxStrainRatioAux);
registerMooseObject("LynxApp", LynxStrainRatioAux);
// template <>
// InputParameters
// validParams<LynxStrainRatioAux>()
// {
// InputParameters params = validParams<LynxStrainAuxBase>();
// params.addClassDescription("Access the strain ratio (volumetric strain on norm of the
// strain)."); return params;
// }
template <>
InputParameters
validParams<LynxStrainRatioAux>()
{
InputParameters params = validParams<AuxKernel>();
params.addClassDescription(
"Access the strain ratio (volumetric strain on norm of the elastic strain).");
return params;
}
// LynxStrainRatioAux::LynxStrainRatioAux(const InputParameters & parameters)
// : LynxStrainAuxBase(parameters)
// {
// }
LynxStrainRatioAux::LynxStrainRatioAux(const InputParameters & parameters)
: DerivativeMaterialInterface<AuxKernel>(parameters),
_elastic_strain(getDefaultMaterialProperty<RankTwoTensor>("elastic_strain"))
{
}
// Real
// LynxStrainRatioAux::computeValue()
// {
// Real strain_norm = (*_strain)[_qp].L2norm();
Real
LynxStrainRatioAux::computeValue()
{
Real strain_norm = _elastic_strain[_qp].L2norm();
// if (strain_norm != 0.0)
// return (*_strain)[_qp].trace() / strain_norm;
// else
// return -std::sqrt(3.0);
// }
\ No newline at end of file
if (strain_norm != 0.0)
return _elastic_strain[_qp].trace() / strain_norm;
else
return -std::sqrt(3.0);
}
/******************************************************************************/
/* LYNX, a MOOSE-based application */
/* */
/* Copyright (C) 2017 by Antoine B. Jacquey and Mauro Cacace */
/* GFZ Potsdam, German Research Centre for Geosciences */
/* */
/* This program is free software: you can redistribute it and/or modify */
/* it under the terms of the GNU General Public License as published by */
/* the Free Software Foundation, either version 3 of the License, or */
/* (at your option) any later version. */
/* */
/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
/* GNU General Public License for more details. */
/* */
/* You should have received a copy of the GNU General Public License */
/* along with this program. If not, see <http://www.gnu.org/licenses/> */
/******************************************************************************/
#include "LynxDamageRate.h"
registerMooseObject("LynxApp", LynxDamageRate);
template <>
InputParameters
validParams<LynxDamageRate>()
{
InputParameters params = validParams<Kernel>();
params.addClassDescription("Damage rate for damage rheology.");
params.addCoupledVar("damage_rate", "The damage rate auxiliary variable");
return params;
}
LynxDamageRate::LynxDamageRate(const InputParameters & parameters)
: Kernel(parameters),
_u_old(valueOld()),
_coupled_dam(isCoupled("damage_rate")),
_damage_rate(_coupled_dam ? coupledValue("damage_rate") : _zero)
{
}
/******************************************************************************/
/* RESIDUALS */
/******************************************************************************/
Real
LynxDamageRate::computeQpResidual()
{
Real rate = std::min(_damage_rate[_qp], (1.0 - _u_old[_qp]) / _dt);
return -rate * _test[_i][_qp];
}
/******************************************************************************/
/* JACOBIAN */
/******************************************************************************/
Real
LynxDamageRate::computeQpJacobian()
{
return 0.0;
}
/******************************************************************************/
/* OFF-DIAG JACOBIAN */
/******************************************************************************/
Real
LynxDamageRate::computeQpOffDiagJacobian(unsigned int jvar)
{
return 0.0;
}
......@@ -40,6 +40,9 @@ validParams<LynxDamageDeformation>()
params.addParam<std::vector<Real>>(
"friction_angle",
"The friction angle of the material for the pressure-dependant part of the yield stress.");
params.addParam<std::vector<Real>>("cohesion",
"The constant coefficient of the yield stress corresponding "
"to the cohesion of the material.");
params.addParam<std::vector<Real>>(
"porous_coeff",
"The porous coefficient controlling the capped yield in the damage rheology.");
......@@ -68,6 +71,8 @@ LynxDamageDeformation::LynxDamageDeformation(const InputParameters & parameters)
// Damage-Plasticity parameters
_friction_angle(isParamValid("friction_angle") ? getLynxParam<Real>("friction_angle")
: std::vector<Real>(_n_composition, 0.0)),
_cohesion(isParamValid("cohesion") ? getLynxParam<Real>("cohesion")
: std::vector<Real>(_n_composition, 0.0)),
_porous_coeff(isParamValid("porous_coeff") ? getLynxParam<Real>("porous_coeff")
: std::vector<Real>(_n_composition, 0.0)),
_porous_coeff_linear(isParamValid("porous_coeff_linear")
......@@ -79,7 +84,11 @@ LynxDamageDeformation::LynxDamageDeformation(const InputParameters & parameters)
: std::vector<Real>(_n_composition, 0.0)),
// Strain properties
_elastic_strain(declareProperty<RankTwoTensor>("elastic_strain")),
_elastic_strain_old(getMaterialPropertyOld<RankTwoTensor>("elastic_strain"))
_elastic_strain_old(getMaterialPropertyOld<RankTwoTensor>("elastic_strain")),
// Stress properties
// _dstress_ddamage(declareProperty<RankTwoTensor>("dstress_ddamage")),
// Damage properties
_damage_rate(declareProperty<Real>("damage_rate"))
{
_has_plasticity = isParamValid("friction_angle");
......@@ -115,6 +124,10 @@ LynxDamageDeformation::initializeQpDeformation()
// Need to get _xi0
initializeDamageParameters();
// Initialize yield derivative
_dyield_dp_tr = 0.0;
_dyield_dq_tr = 0.0;
LynxDeformationBase::initializeQpDeformation();
}
......@@ -199,12 +212,55 @@ LynxDamageDeformation::plasticCorrection(Real & pressure, RankTwoTensor & stress
void
LynxDamageDeformation::damageCorrection()
{
// Rotate old elastic strain
RankTwoTensor strain_el = spinRotation(_elastic_strain_old[_qp].deviatoric());
strain_el.addIa(_elastic_strain_old[_qp].trace() / 3.0);
Real xi = strainRatio(_elastic_strain_old[_qp]);
// Build damage correction to the elasticity tensor
RankTwoTensor e = (strain_el.L2norm() != 0.0) ? strain_el / strain_el.L2norm() : RankTwoTensor();
_damaged_tensor =
_damage_plasticity->_gamma *
(e.outerProduct(_identity_two) + _identity_two.outerProduct(e) - xi * e.outerProduct(e));
// Damage stress
_damaged_stress = _damage_plasticity->_gamma * (xi - 2.0 * _damage_plasticity->_xi0) * strain_el;
_damaged_stress.addIa(_damage_plasticity->_gamma * _elastic_strain_old[_qp].L2norm());
// _dstress_ddamage[_qp] = -damaged_stress;
// Correct stress
_stress[_qp] -=
_damage_old[_qp] * _damaged_tensor *
(_strain_increment[_qp] - _viscous_strain_incr[_qp] - _plastic_strain_incr[_qp]) +
_damaged_stress * (_damage[_qp] - _damage_old[_qp]);
// Damage rate
if ((_plastic_yield_function[_qp] > 0.0))
_damage_rate[_qp] = _damage_plasticity->_eta_d * _plastic_yield_function[_qp];
else
_damage_rate[_qp] = 0.0;
}
RankFourTensor
LynxDamageDeformation::damageTangentOperator(const RankTwoTensor & flow_direction,
const RankFourTensor & tme)
{
RankTwoTensor dyield_dstrain_tr = RankTwoTensor();
if ((_damage_rate[_qp] > 0.0) && (_damage_old[_qp] > 0.0) && (_damage[_qp] != _damage_old[_qp]))
{
dyield_dstrain_tr = 3.0 * _G[_qp] * _dyield_dq_tr * flow_direction;
dyield_dstrain_tr.addIa(-_K[_qp] * _dyield_dp_tr);
}
return _damage_old[_qp] * _damaged_tensor * tme +
_damage_plasticity->_eta_d * _damaged_stress.outerProduct(dyield_dstrain_tr) * _dt;
}
Real
LynxDamageDeformation::computePlasticityYield(const Real & pressure, const Real & eqv_stress)
{
return eqv_stress - _damage_plasticity->_alpha0 * pressure;
return eqv_stress - _damage_plasticity->_alpha0 * pressure - _damage_plasticity->_k0;
}
Real
......@@ -219,19 +275,30 @@ LynxDamageDeformation::plasticIncrement(const Real & /*pressure*/, const Real &
if (_damage_plasticity->_eta_p != 0.0)
{
vp_correction = (3.0 * _G[_qp]) * _dt * _damage_plasticity->_eta_p /
(1.0 + (3.0 * _G[_qp]) * _dt * _damage_plasticity->_eta_p);
(1.0 + 3.0 * _G[_qp] * _dt * _damage_plasticity->_eta_p);
eqv_p_strain_incr *= vp_correction;
}
// Update yield
_plastic_yield_function[_qp] -= (3.0 * _G[_qp]) * eqv_p_strain_incr;
_stress_corr_p =
(eqv_stress != 0.0) ? (eqv_stress - 3.0 * _G[_qp] * eqv_p_strain_incr) / eqv_stress : 1.0;
_dp_dp_tr_p = 1.0;
_dp_dq_tr_p = 0.0;
_dq_dp_tr_p = _damage_plasticity->_alpha0 * vp_correction;
_dq_dq_tr_p = 1.0 - vp_correction;
// Tangent operator
if (_fe_problem.currentlyComputingJacobian())
{
_stress_corr_p =
(eqv_stress != 0.0) ? (eqv_stress - 3.0 * _G[_qp] * eqv_p_strain_incr) / eqv_stress : 1.0;
_dp_dp_tr_p = 1.0;
_dp_dq_tr_p = 0.0;
_dq_dp_tr_p = _damage_plasticity->_alpha0 * vp_correction;
_dq_dq_tr_p = 1.0 - vp_correction;