Commit 2bb1f3c8 authored by Antoine Jacquey's avatar Antoine Jacquey

Stoke update in LynxADCreepModel. Fixed volumetric plastic strain incr.

parent 7a8f8b29
......@@ -39,6 +39,9 @@ public:
protected:
virtual ADReal viscousIncrement(const ADReal & pressure, const ADReal & eqv_stress, const ADReal & G);
virtual ADReal computeQpEffectiveViscosity(const ADReal & pressure);
virtual ADReal computeQpOneOnDiffViscosity(const ADReal A);
virtual ADReal computeQpOneOnDislViscosity(const ADReal A, const ADReal n, const ADReal eII);
virtual void initCreepParameters(const ADReal & pressure);
virtual ADReal iterativeResidual(const ADReal & x);
virtual ADReal iterativeResidualDerivative(const ADReal & x);
......@@ -60,6 +63,10 @@ protected:
const std::vector<Real> _E_dislocation;
const std::vector<Real> _V_dislocation;
const Real _gas_constant;
const bool _has_background_strain_rate;
const bool _has_initial_viscosity;
const std::vector<Real> _initial_viscosity;
const Real _background_strain_rate;
const std::vector<Real> _eta_min;
const std::vector<Real> _eta_max;
const unsigned int _viscous_update;
......
......@@ -38,6 +38,10 @@ defineADValidParams(
params.addParam<std::vector<Real>>("n_dislocation",
"List of power law exponents for dislocation creep.");
params.addParam<Real>("gas_constant", 8.3144621, "The universal gas constant");
params.addParam<std::vector<Real>>(
"initial_viscosity", "The vector of initial viscosity for purely viscous deformation.");
params.addParam<Real>("background_strain_rate",
"The background strain rate for purely viscous deformation.");
params.addParam<std::vector<Real>>(
"eta_min",
"The lower threshold for the effective viscosity for purely viscous deformation.");
......@@ -84,6 +88,13 @@ LynxADCreepModel<compute_stage>::LynxADCreepModel(const InputParameters & parame
? this->getLynxParam("V_dislocation")
: std::vector<Real>(_n_composition, 0.0)),
_gas_constant(getParam<Real>("gas_constant")),
_has_background_strain_rate(isParamValid("background_strain_rate")),
_has_initial_viscosity(_has_background_strain_rate ? false
: isParamValid("_has_initial_viscosity")),
_initial_viscosity(_has_initial_viscosity ? this->getLynxParam("initial_viscosity")
: std::vector<Real>(_n_composition, 0.0)),
_background_strain_rate(_has_background_strain_rate ? getParam<Real>("background_strain_rate")
: 0.0),
_eta_min(isParamValid("eta_min") ? this->getLynxParam("eta_min")
: std::vector<Real>(_n_composition, 0.0)),
_eta_max(isParamValid("eta_max") ? this->getLynxParam("eta_max")
......@@ -111,14 +122,25 @@ LynxADCreepModel<compute_stage>::creepUpdate(ADRankTwoTensor & stress_dev,
const ADReal & G,
ADRankTwoTensor & elastic_strain_incr)
{
const ADReal tau_II_tr = std::sqrt(0.5) * stress_dev.L2norm();
const ADRankTwoTensor flow_dir = (tau_II_tr != 0.0) ? stress_dev / tau_II_tr : ADRankTwoTensor();
if (G != 0.0) // Visco elastic update
{
const ADReal tau_II_tr = std::sqrt(0.5) * stress_dev.L2norm();
const ADRankTwoTensor flow_dir = (tau_II_tr != 0.0) ? stress_dev / tau_II_tr : ADRankTwoTensor();
ADReal delta_e_II = viscousIncrement(pressure, tau_II_tr, G);
ADReal delta_e_II = viscousIncrement(pressure, tau_II_tr, G);
_viscous_strain_incr[_qp] = delta_e_II * flow_dir;
stress_dev -= 2.0 * G * delta_e_II * flow_dir;
elastic_strain_incr -= _viscous_strain_incr[_qp];
}
else // Stoke
{
_viscous_strain_incr[_qp] = elastic_strain_incr;
// Compute effective viscosity
_eta_eff[_qp] = computeQpEffectiveViscosity(pressure);
_viscous_strain_incr[_qp] = delta_e_II * flow_dir;
stress_dev -= 2.0 * G * delta_e_II * flow_dir;
elastic_strain_incr -= _viscous_strain_incr[_qp];
stress_dev = 2.0 * _eta_eff[_qp] * _viscous_strain_incr[_qp].deviatoric() / _dt;
}
}
template <ComputeStage compute_stage>
......@@ -187,6 +209,53 @@ LynxADCreepModel<compute_stage>::viscousIncrement(const ADReal & pressure, const
}
}
template <ComputeStage compute_stage>
ADReal
LynxADCreepModel<compute_stage>::computeQpEffectiveViscosity(const ADReal & pressure)
{
// Map creep parameters
initCreepParameters(pressure);
ADReal strain_rate_II = std::sqrt(0.5) * _viscous_strain_incr[_qp].deviatoric().L2norm() / _dt;
if (_t_step <= 1 && strain_rate_II == 0.0 && _has_initial_viscosity)
return std::min(std::max(this->averageProperty(_initial_viscosity), this->averageProperty(_eta_min)),
this->averageProperty(_eta_max));
strain_rate_II = (_has_background_strain_rate && _t_step <= 1 && strain_rate_II == 0.0)
? _background_strain_rate
: strain_rate_II;
ADReal one_on_eta_diff = 0.0, one_on_eta_disl = 0.0;
if (_has_diffusion_creep)
one_on_eta_diff = computeQpOneOnDiffViscosity(_A_diff);
if (_has_dislocation_creep)
one_on_eta_disl = computeQpOneOnDislViscosity(_A_disl, _n_disl, strain_rate_II);
ADReal eta = 1.0 / (one_on_eta_diff + one_on_eta_disl);
return std::min(std::max(eta, this->averageProperty(_eta_min)), this->averageProperty(_eta_max));
}
template <ComputeStage compute_stage>
ADReal
LynxADCreepModel<compute_stage>::computeQpOneOnDiffViscosity(const ADReal A)
{
return 2.0 * A;
}
template <ComputeStage compute_stage>
ADReal
LynxADCreepModel<compute_stage>::computeQpOneOnDislViscosity(
const ADReal A, const ADReal n, const ADReal eII)
{
if ((eII == 0.0) && (n == 1.0))
return 2.0;
else
return 2.0 * std::pow(A, 1.0 / n) * std::pow(eII, 1.0 - 1.0 / n);
}
template <ComputeStage compute_stage>
void
LynxADCreepModel<compute_stage>::initCreepParameters(const ADReal & pressure)
......
......@@ -23,7 +23,7 @@ defineADValidParams(
params.addClassDescription("Class calculating strain and stress for an elastic rheology.");
// Elastic moduli parameters
params.addRequiredRangeCheckedParam<std::vector<Real>>("bulk_modulus",
"bulk_modulus >= 0.0",
"bulk_modulus > 0.0",
"The drained bulk modulus of the material.");
params.addRequiredRangeCheckedParam<std::vector<Real>>("shear_modulus",
"shear_modulus >= 0.0",
......
......@@ -111,7 +111,7 @@ LynxADPlasticModel<compute_stage>::plasticUpdate(ADRankTwoTensor & stress_dev,
ADReal delta_e_eqv = plasticIncrement(eqv_stress, pressure, G, K);
_plastic_strain_incr[_qp] = 1.5 * delta_e_eqv * flow_dir;
_plastic_strain_incr[_qp].addIa(_beta * delta_e_eqv);
_plastic_strain_incr[_qp].addIa(_beta * delta_e_eqv / 3.0);
stress_dev -= 3.0 * G * delta_e_eqv * flow_dir;
pressure += K * _beta * delta_e_eqv;
elastic_strain_incr -= _plastic_strain_incr[_qp];
......@@ -132,7 +132,7 @@ LynxADPlasticModel<compute_stage>::plasticIncrement(const ADReal & eqv_stress, c
initPlasticParameters(pressure, K);
// If no friction and cohesion are provided, there is no plastic update
if (_alpha == 0.0 && _k == 0.0)
if ((_alpha == 0.0 && _k == 0.0) || G == 0.0)
{
_plastic_yield_function[_qp] = -1.0;
return 0.0;
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment