Commit 12f04b84 authored by Antoine Jacquey's avatar Antoine Jacquey

Damage rheology working and tested.

parent 817cdaaf
/******************************************************************************/
/* 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 LYNXDAMAGEAUX_H
#define LYNXDAMAGEAUX_H
#include "AuxKernel.h"
#include "DerivativeMaterialInterface.h"
class LynxDamageAux;
template <>
InputParameters validParams<LynxDamageAux>();
class LynxDamageAux : public DerivativeMaterialInterface<AuxKernel>
{
public:
LynxDamageAux(const InputParameters & parameters);
protected:
virtual Real computeValue() override;
const MaterialProperty<Real> & _damage_rate;
};
#endif // LYNXDAMAGEAUX_H
\ No newline at end of file
......@@ -22,13 +22,14 @@
#define LYNXDAMAGERATE_H
#include "Kernel.h"
#include "DerivativeMaterialInterface.h"
class LynxDamageRate;
template <>
InputParameters validParams<LynxDamageRate>();
class LynxDamageRate : public Kernel
class LynxDamageRate : public DerivativeMaterialInterface<Kernel>
{
public:
LynxDamageRate(const InputParameters & parameters);
......@@ -40,7 +41,12 @@ protected:
const VariableValue & _u_old;
bool _coupled_dam;
bool _coupled_disp;
unsigned int _ndisp;
std::vector<unsigned> _disp_var;
const VariableValue & _damage_rate;
const MaterialProperty<Real> & _damage_rate_mat;
const MaterialProperty<RankTwoTensor> & _ddamage_rate_dstrain;
};
#endif // LYNXDAMAGERATE_H
......@@ -58,6 +58,7 @@ protected:
const VariableValue & _pf;
bool _coupled_plith;
bool _coupled_pdyn;
bool _coupled_dam;
const unsigned int _component;
bool _vol_locking_correction;
const MaterialProperty<RankTwoTensor> & _stress;
......@@ -68,6 +69,7 @@ protected:
const MaterialProperty<Real> & _rho_b;
const MaterialProperty<Real> & _drho_dtemp;
const MaterialProperty<Real> & _drho_dev;
const MaterialProperty<RankTwoTensor> & _dstress_ddamage;
// Gradient of test function averaged over the element. Used in volumetric locking correction
// calculation.
std::vector<std::vector<Real>> _avg_grad_test;
......@@ -79,6 +81,7 @@ protected:
unsigned int _pf_var;
unsigned int _plith_var;
unsigned int _pdyn_var;
unsigned int _damage_var;
};
#endif // LYNXSOLIDMOMENTUM_H
......@@ -35,13 +35,10 @@ public:
virtual ~LynxDamageDeformation() {}
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 RankFourTensor damageTangentOperator(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);
......@@ -83,15 +80,10 @@ protected:
// 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;
MaterialProperty<RankTwoTensor> & _dstress_ddamage;
MaterialProperty<RankTwoTensor> & _ddamage_rate_dstrain;
// Damage properties
MaterialProperty<Real> & _damage_rate;
......
/******************************************************************************/
/* 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
......@@ -51,8 +51,8 @@ protected:
virtual void computeQpDeformation();
virtual void initializeQpDeformation();
virtual void elasticModuli();
virtual void volumetricDeformation(Real & pressure);
virtual void deviatoricDeformation(const Real & pressure, RankTwoTensor & stress_dev);
virtual Real volumetricDeformation();
virtual RankTwoTensor deviatoricDeformation(const Real & pressure);
virtual void reformStressTensor(const Real & pressure, const RankTwoTensor & stress_dev);
virtual void plasticCorrection(Real & pressure, RankTwoTensor & stress_dev) = 0;
virtual void damageCorrection();
......@@ -60,8 +60,7 @@ protected:
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 RankFourTensor damageTangentOperator(const RankFourTensor & /*tme*/);
virtual void finiteTangentOperator();
virtual Real computeStokeEffectiveViscosity(const Real & pressure);
virtual RankTwoTensor computeStokeEffectiveViscosityDerivative();
......@@ -150,6 +149,8 @@ protected:
Real _dq_dq_tr_p;
// Strain properties
MaterialProperty<RankTwoTensor> & _elastic_strain;
const MaterialProperty<RankTwoTensor> & _elastic_strain_old;
MaterialProperty<RankTwoTensor> & _strain_increment;
MaterialProperty<RankTwoTensor> & _spin_tensor;
const MaterialProperty<RankTwoTensor> & _thermal_strain_incr;
......@@ -164,7 +165,6 @@ protected:
// Stress properties
MaterialProperty<RankTwoTensor> & _stress;
const MaterialProperty<RankTwoTensor> & _stress_old;
MaterialProperty<Real> & _K;
MaterialProperty<Real> & _G;
MaterialProperty<RankFourTensor> & _tangent_modulus;
......
......@@ -38,7 +38,6 @@ protected:
virtual Real computeDT() override;
private:
const Real _log_dt;
const Real _first_dt;
const Real _max_dt;
......
/******************************************************************************/
/* 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 "LynxDamageAux.h"
registerMooseObject("LynxApp", LynxDamageAux);
template <>
InputParameters
validParams<LynxDamageAux>()
{
InputParameters params = validParams<AuxKernel>();
params.addClassDescription("Updating the damage auxiliary variable.");
return params;
}
LynxDamageAux::LynxDamageAux(const InputParameters & parameters)
: DerivativeMaterialInterface<AuxKernel>(parameters),
_damage_rate(getDefaultMaterialProperty<Real>("damage_rate"))
{
}
Real
LynxDamageAux::computeValue()
{
Real damage_rate = _damage_rate[_qp];
if (damage_rate > (1.0 - _u_old[_qp]) / _dt)
damage_rate = (1.0 - _u_old[_qp]) / _dt;
return _u_old[_qp] + damage_rate * _dt;
}
\ No newline at end of file
......@@ -27,8 +27,7 @@ InputParameters
validParams<LynxEqvStrainAux>()
{
InputParameters params = validParams<LynxStrainAuxBase>();
params.addClassDescription(
"Calculates the equivalent strain of the given tensor.");
params.addClassDescription("Calculates the equivalent strain of the given tensor.");
return params;
}
......
......@@ -27,8 +27,7 @@ InputParameters
validParams<LynxHeatFluxAux>()
{
InputParameters params = validParams<AuxKernel>();
params.addClassDescription(
"Calculates the heat flux in each element for the given direction.");
params.addClassDescription("Calculates the heat flux in each element for the given direction.");
params.addRequiredCoupledVar("temperature", "The temperature variable.");
params.addRequiredRangeCheckedParam<unsigned int>(
"component",
......@@ -54,5 +53,5 @@ LynxHeatFluxAux::computeValue()
else
_cond_bulk = _thermal_diff[_qp] * _rhoC_b[_qp];
return -1*_cond_bulk*_grad_T[_qp](_component);
return -1 * _cond_bulk * _grad_T[_qp](_component);
}
......@@ -27,15 +27,13 @@ InputParameters
validParams<LynxTemperatureCelsiusAux>()
{
InputParameters params = validParams<AuxKernel>();
params.addRequiredCoupledVar("temperature",
"The temperature variable to convert to Celsius.");
params.addRequiredCoupledVar("temperature", "The temperature variable to convert to Celsius.");
params.addClassDescription("Converts Kelvin to degree Celsius.");
return params;
}
LynxTemperatureCelsiusAux::LynxTemperatureCelsiusAux(const InputParameters & parameters)
: AuxKernel(parameters),
_T_K(coupledValue("temperature"))
: AuxKernel(parameters), _T_K(coupledValue("temperature"))
{
}
......
......@@ -30,7 +30,9 @@ validParams<LynxVariableRateAux>()
params.addClassDescription("Computes the rate of change of a coupled scalar variable.");
params.addRequiredCoupledVar("coupled_variable", "The coupled variable.");
params.addParam<Real>("time_scale_factor", 1.0, "Multiply the time step size by this factor.");
params.addParam<bool>("output_relative", false,
params.addParam<bool>(
"output_relative",
false,
"If true, write gradient relative to the absolute value of coupled_variable.");
return params;
}
......@@ -48,7 +50,7 @@ Real
LynxVariableRateAux::computeValue()
{
if (_relative)
return (_cvar[_qp] - _cvar_old[_qp])/_cvar_old[_qp]/_dt/_tscale;
return (_cvar[_qp] - _cvar_old[_qp]) / _cvar_old[_qp] / _dt / _tscale;
else
return (_cvar[_qp] - _cvar_old[_qp])/_dt/_tscale;
return (_cvar[_qp] - _cvar_old[_qp]) / _dt / _tscale;
}
......@@ -28,11 +28,13 @@ InputParameters
validParams<LynxWinklerBC>()
{
InputParameters params = validParams<IntegratedBC>();
params.addClassDescription("Applies a Winkler foundation BC on a given boundary in a given direction.");
params.addClassDescription(
"Applies a Winkler foundation BC on a given boundary in a given direction.");
params.addRequiredCoupledVar(
"displacements",
"The displacements appropriate for the simulation geometry and coordinate system.");
params.addRequiredParam<unsigned int>("component", "The component of the displacement orthogonal to the boundary.");
params.addRequiredParam<unsigned int>(
"component", "The component of the displacement orthogonal to the boundary.");
params.addParam<Real>("value", 0.0, "Value of the reference pressure.");
params.addParam<FunctionName>("function", "The function that describes the reference pressure.");
params.addParam<Real>("external_density", 0.0, "The density of the external material.");
......
......@@ -27,17 +27,34 @@ InputParameters
validParams<LynxDamageRate>()
{
InputParameters params = validParams<Kernel>();
params.addClassDescription("Damage rate for damage rheology.");
params.addClassDescription(
"Damage rate for damage rheology. Works fully coupled way or in a MultiApp");
params.addCoupledVar("damage_rate", "The damage rate auxiliary variable");
params.addCoupledVar("displacements", "The displacements variables.");
return params;
}
LynxDamageRate::LynxDamageRate(const InputParameters & parameters)
: Kernel(parameters),
: DerivativeMaterialInterface<Kernel>(parameters),
_u_old(valueOld()),
_coupled_dam(isCoupled("damage_rate")),
_damage_rate(_coupled_dam ? coupledValue("damage_rate") : _zero)
_coupled_disp(isCoupled("displacements")),
_ndisp(coupledComponents("displacements")),
_disp_var(_ndisp),
_damage_rate(_coupled_dam ? coupledValue("damage_rate") : _zero),
_damage_rate_mat(getDefaultMaterialProperty<Real>("damage_rate")),
_ddamage_rate_dstrain(getDefaultMaterialProperty<RankTwoTensor>("ddamage_rate_dstrain"))
{
if (_coupled_dam && _coupled_disp)
mooseError("LynxDamageRate: you provided both damage_rate and displacements. Choose either "
"damage_rate if running in a subApp or displacements if running a fully coupled "
"simulation.");
else if (!_coupled_dam && !_coupled_disp)
mooseError("LynxDamageRate: you need to provide either damage_rate for running in a subApp or "
"displacements for fully coupled!");
for (unsigned i = 0; i < _ndisp; ++i)
_disp_var[i] = _coupled_disp ? coupled("displacements", i) : -1;
}
/******************************************************************************/
......@@ -47,7 +64,12 @@ LynxDamageRate::LynxDamageRate(const InputParameters & parameters)
Real
LynxDamageRate::computeQpResidual()
{
Real rate = std::min(_damage_rate[_qp], (1.0 - _u_old[_qp]) / _dt);
Real rate = (1.0 - _u_old[_qp]) / _dt;
if (_coupled_dam)
rate = std::min(_damage_rate[_qp], rate);
else if (_coupled_disp)
rate = std::min(_damage_rate_mat[_qp], rate);
return -rate * _test[_i][_qp];
}
......@@ -68,5 +90,20 @@ LynxDamageRate::computeQpJacobian()
Real
LynxDamageRate::computeQpOffDiagJacobian(unsigned int jvar)
{
return 0.0;
Real jac = 0.0;
for (unsigned int coupled_component = 0; coupled_component < _ndisp; ++coupled_component)
if (_coupled_disp && (jvar == _disp_var[coupled_component]))
{
Real rate = (1.0 - _u_old[_qp]) / _dt;
if (_coupled_dam)
rate = std::min(_damage_rate[_qp], rate);
else if (_coupled_disp)
rate = std::min(_damage_rate_mat[_qp], rate);
if (rate != ((1.0 - _u_old[_qp]) / _dt))
jac += -_ddamage_rate_dstrain[_qp].row(coupled_component) * _grad_phi[_j][_qp] *
_test[_i][_qp];
}
return jac;
}
......@@ -39,6 +39,7 @@ validParams<LynxSolidMomentum>()
params.addCoupledVar("fluid_pressure", "The fluid pressure variable.");
params.addCoupledVar("lithostatic_pressure", "The lithostatic pressure variable.");
params.addCoupledVar("dynamic_pressure", "The dynamic pressure variable.");
params.addCoupledVar("damage", "The damage intensity variable.");
params.set<bool>("use_displaced_mesh") = false;
params.addRequiredParam<unsigned int>("component",
"An integer corresponding to the direction "
......@@ -57,6 +58,7 @@ LynxSolidMomentum::LynxSolidMomentum(const InputParameters & parameters)
_pf(_coupled_pf ? coupledValue("fluid_pressure") : _zero),
_coupled_plith(isCoupled("lithostatic_pressure")),
_coupled_pdyn(isCoupled("dynamic_pressure")),
_coupled_dam(isCoupled("damage")),
_component(getParam<unsigned int>("component")),
_vol_locking_correction(getParam<bool>("volumetric_locking_correction")),
_stress(getMaterialProperty<RankTwoTensor>("stress")),
......@@ -67,13 +69,15 @@ LynxSolidMomentum::LynxSolidMomentum(const InputParameters & parameters)
_rho_b(getDefaultMaterialProperty<Real>("bulk_density")),
_drho_dtemp(getDefaultMaterialProperty<Real>("drho_dtemp")),
_drho_dev(getDefaultMaterialProperty<Real>("drho_dev")),
_dstress_ddamage(getDefaultMaterialProperty<RankTwoTensor>("dstress_ddamage")),
_avg_grad_test(_test.size(), std::vector<Real>(3, 0.0)),
_avg_grad_phi(_phi.size(), std::vector<Real>(3, 0.0)),
_disp_var(_ndisp),
_temp_var(_coupled_temp ? coupled("temperature") : 0),
_pf_var(_coupled_pf ? coupled("fluid_pressure") : 0),
_plith_var(_coupled_plith ? coupled("lithostatic_pressure") : 0),
_pdyn_var(_coupled_pdyn ? coupled("dynamic_pressure") : 0)
_pdyn_var(_coupled_pdyn ? coupled("dynamic_pressure") : 0),
_damage_var(_coupled_dam ? coupled("damage") : 0)
{
// Checking for consistency between mesh size and length of the provided displacements vector
if (_ndisp != _mesh.dimension())
......@@ -274,6 +278,9 @@ LynxSolidMomentum::computeQpOffDiagJacobian(unsigned int jvar)
if (_coupled_pdyn && jvar == _pdyn_var)
return -_phi[_j][_qp] * _grad_test[_i][_qp](_component);
if (_coupled_dam && jvar == _damage_var)
return _dstress_ddamage[_qp].row(_component) * _phi[_j][_qp] * _grad_test[_i][_qp];