Commit 478519da authored by Antoine Jacquey's avatar Antoine Jacquey

Damage model with AD.

parent 2890cf40
......@@ -6,7 +6,7 @@
"type": "cppdbg",
"request": "launch",
"program": "${workspaceFolder}/lynx-dbg",
"args": ["-i", "test/tests/thermo_mech/thermal_expansion.i"],
"args": ["-i", "test/tests/damage/damage_shear_no_dam.i"],
"stopAtEntry": false,
"cwd": "${workspaceFolder}",
"environment": [],
......
......@@ -21,8 +21,10 @@ class LynxADDamageDeformation : public LynxADElasticDeformation
public:
static InputParameters validParams();
LynxADDamageDeformation(const InputParameters & parameters);
void initialSetup() override;
protected:
virtual void initQpStatefulProperties() override;
virtual void initializeQpDeformation() override;
virtual void computeQpStress() override;
......@@ -34,5 +36,5 @@ protected:
const MaterialProperty<RankTwoTensor> & _elastic_strain_old;
// Elasticity tensor
ADRankFourTensor _Cijkl;
RankFourTensor _Cijkl;
};
\ No newline at end of file
/******************************************************************************/
/* This file is part of */
/* LYNX, a MOOSE-based application */
/* Lithosphere dYnamic Numerical toolboX */
/* */
/* Copyright (C) 2017 by Antoine B. Jacquey and Mauro Cacace */
/* GFZ Potsdam, German Research Centre for Geosciences */
/* */
/* Licensed under GNU General Public License 3, */
/* please see LICENSE for details */
/* or http://www.gnu.org/licenses/gpl.html */
/******************************************************************************/
#pragma once
#include "LynxADMaterialBase.h"
class LynxADDamageModelBase : public LynxADMaterialBase
{
public:
static InputParameters validParams();
LynxADDamageModelBase(const InputParameters & parameters);
void setQp(unsigned int qp);
virtual void elasticGuess(ADRankTwoTensor & stress,
const ADRankTwoTensor & stress_old,
const RankFourTensor & Cijkl,
const ADRankTwoTensor & elastic_strain_old,
const ADRankTwoTensor & elastic_strain_incr) = 0;
virtual void damageUpdate(ADRankTwoTensor & stress, ADRankTwoTensor & elastic_strain_incr) = 0;
void resetQpProperties() final {}
void resetProperties() final {}
protected:
virtual void initQpStatefulProperties() override;
const Real _abs_tol;
const Real _rel_tol;
const unsigned int _max_its;
const std::vector<Real> _damage0;
// Damage properties
ADMaterialProperty<Real> & _damage;
const MaterialProperty<Real> & _damage_old;
ADMaterialProperty<Real> & _damage_drive;
ADMaterialProperty<RankTwoTensor> & _plastic_strain_incr;
ADMaterialProperty<Real> & _damage_incr;
ADMaterialProperty<Real> & _yield_function;
};
\ No newline at end of file
/******************************************************************************/
/* This file is part of */
/* LYNX, a MOOSE-based application */
/* Lithosphere dYnamic Numerical toolboX */
/* */
/* Copyright (C) 2017 by Antoine B. Jacquey and Mauro Cacace */
/* GFZ Potsdam, German Research Centre for Geosciences */
/* */
/* Licensed under GNU General Public License 3, */
/* please see LICENSE for details */
/* or http://www.gnu.org/licenses/gpl.html */
/******************************************************************************/
#pragma once
#include "LynxADDamageModelBase.h"
class LynxADLyakhovskyDamage : public LynxADDamageModelBase
{
public:
static InputParameters validParams();
LynxADLyakhovskyDamage(const InputParameters & parameters);
virtual void elasticGuess(ADRankTwoTensor & stress,
const ADRankTwoTensor & stress_old,
const RankFourTensor & Cijkl,
const ADRankTwoTensor & elastic_strain_old,
const ADRankTwoTensor & elastic_strain_incr) override;
virtual void damageUpdate(ADRankTwoTensor & stress,
ADRankTwoTensor & elastic_strain_incr) override;
protected:
virtual ADRankFourTensor damageStiffness(const RankFourTensor & Cijkl);
virtual void returnMap(ADReal & gamma_s, ADReal & gamma_d);
virtual void
residual(const ADReal & gamma_vs, const ADReal & gamma_d, ADReal & ress, ADReal & resd);
virtual void jacobian(const ADReal & gamma_s,
const ADReal & gamma_d,
ADReal & jacss,
ADReal & jacdd,
ADReal & jacsd,
ADReal & jacds);
virtual void
overStress(const ADReal & gamma_s, const ADReal & gamma_d, ADReal & over_s, ADReal & over_d);
virtual void overStressDerivS(const ADReal & gamma_s,
const ADReal & gamma_d,
ADReal & over_s_s,
ADReal & over_d_s);
virtual void overStressDerivD(const ADReal & gamma_s,
const ADReal & gamma_d,
ADReal & over_s_d,
ADReal & over_d_d);
virtual void updateYieldParameters(const ADReal & gamma_s, const ADReal & gamma_d);
virtual void updateYieldParametersDerivS(const ADReal & gamma_s, const ADReal & gamma_d);
virtual void updateYieldParametersDerivD(const ADReal & gamma_s, const ADReal & gamma_d);
virtual ADReal flowRate();
virtual ADReal flowDirectionS();
virtual ADReal flowDirectionD();
virtual ADReal dFlowRatedS();
virtual ADReal dFlowRatedD();
virtual ADReal dFlowDirectionSdS();
virtual ADReal dFlowDirectionSdD();
virtual ADReal dFlowDirectionDdS();
virtual ADReal dFlowDirectionDdD();
virtual ADReal yieldFunction();
virtual ADReal dYieldFunctiondDev();
virtual ADReal dYieldFunctiondPres();
virtual ADReal dYieldFunctiondDam();
virtual ADReal dYieldFunctiondMua();
virtual ADReal d2YieldFunctiondDev2();
virtual ADReal d2YieldFunctiondDevdPres();
virtual ADReal d2YieldFunctiondDevdMua();
virtual ADReal d2YieldFunctiondDamdDev();
virtual ADReal d2YieldFunctiondDamdPres();
virtual ADReal d2YieldFunctiondDamdMua();
virtual ADReal d2YieldFunctiondDam2();
virtual ADReal strainRatio(const ADRankTwoTensor & elastic_strain);
virtual ADRankTwoTensor reformPlasticStrainTensor(const ADReal & gamma_s);
// Damage parameters
const std::vector<Real> _damage_modulus;
const std::vector<Real> _friction_angle;
const std::vector<Real> _cohesion;
const std::vector<Real> _plastic_viscosity;
// Damage utilities
Real _Gam0;
Real _eta_p;
Real _K;
Real _G;
Real _xi0;
Real _k;
ADReal _e_norm;
ADReal _xi;
ADRankTwoTensor _e;
Real _mu0;
ADRankFourTensor _Dijkl;
// Return map utilities
// Trial states
ADRankTwoTensor _stress_tr;
ADReal _pressure_tr;
ADReal _dev_stress_tr;
ADReal _damage_force_tr;
// Forces
ADReal _pressure;
ADReal _dev_stress;
ADReal _damage_force;
ADReal _mua;
// Forces derivatives
ADReal _dpressure_dS;
ADReal _dpressure_dD;
ADReal _ddev_stress_dS;
ADReal _ddev_stress_dD;
ADReal _ddamage_force_dS;
ADReal _ddamage_force_dD;
ADReal _dmua_dS;
ADReal _dmua_dD;
};
\ No newline at end of file
//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
#pragma once
namespace libMesh
{
template <typename>
class VectorValue;
typedef VectorValue<Real> RealGradient;
}
using libMesh::RealGradient;
namespace ElasticityTensorTools
{
/**
* Get the shear modulus for an isotropic elasticity tensor
* param elasticity_tensor the tensor (must be isotropic, but not checked for efficiency)
*/
template <typename T>
T
getIsotropicShearModulus(const RankFourTensorTempl<T> & elasticity_tensor)
{
return elasticity_tensor(0, 1, 0, 1);
}
/**
* Get the bulk modulus for an isotropic elasticity tensor
* param elasticity_tensor the tensor (must be isotropic, but not checked for efficiency)
*/
template <typename T>
T
getIsotropicBulkModulus(const RankFourTensorTempl<T> & elasticity_tensor)
{
const T shear_modulus = getIsotropicShearModulus(elasticity_tensor);
// dilatational modulus is defined as lambda plus two mu
const T dilatational_modulus = elasticity_tensor(0, 0, 0, 0);
const T lambda = dilatational_modulus - 2.0 * shear_modulus;
const T bulk_modulus = lambda + 2.0 * shear_modulus / 3.0;
return bulk_modulus;
}
/**
* Get the Young's modulus for an isotropic elasticity tensor
* param elasticity_tensor the tensor (must be isotropic, but not checked for efficiency)
*/
template <typename T>
T
getIsotropicYoungsModulus(const RankFourTensorTempl<T> & elasticity_tensor)
{
const T shear_modulus = getIsotropicShearModulus(elasticity_tensor);
// dilatational modulus is defined as lambda plus two mu
const T dilatational_modulus = elasticity_tensor(0, 0, 0, 0);
const T lambda = dilatational_modulus - 2.0 * shear_modulus;
const T youngs_modulus =
shear_modulus * (3.0 * lambda + 2.0 * shear_modulus) / (lambda + shear_modulus);
return youngs_modulus;
}
/**
* Get the Poisson's modulus for an isotropic elasticity tensor
* param elasticity_tensor the tensor (must be isotropic, but not checked for efficiency)
*/
template <typename T>
T
getIsotropicPoissonsRatio(const RankFourTensorTempl<T> & elasticity_tensor)
{
const T poissons_ratio = elasticity_tensor(1, 1, 0, 0) /
(elasticity_tensor(1, 1, 1, 1) + elasticity_tensor(1, 1, 0, 0));
return poissons_ratio;
}
/**
* Construct the elasticity tensor based on the bulk and shear modulus
*/
template <typename T>
RankFourTensorTempl<T>
elasticityTensorKandG(const T & K, const T & G)
{
std::vector<T> iso_const(2);
iso_const[0] = K - 2.0 / 3.0 * G;
iso_const[1] = G;
return RankFourTensorTempl<T>(iso_const, RankFourTensorTempl<T>::symmetric_isotropic);
}
}
\ No newline at end of file
......@@ -13,8 +13,7 @@
#include "LynxADDamageDeformation.h"
#include "ElasticityTensorTools.h"
#include "LynxADCreepModel.h"
#include "LynxADPlasticModel.h"
#include "LynxADDamageModelBase.h"
registerMooseObject("LynxApp", LynxADDamageDeformation);
......@@ -23,7 +22,7 @@ LynxADDamageDeformation::validParams()
{
InputParameters params = LynxADElasticDeformation::validParams();
params.addClassDescription("Class calculating the deformation of a material following a "
"damage viscoelasto-(visco)plastic rheology.");
"damage elasto-viscoplastic rheology.");
params.addRequiredParam<MaterialName>(
"damage_model",
"The material object to use for the damage rheology.");
......@@ -32,16 +31,32 @@ LynxADDamageDeformation::validParams()
LynxADDamageDeformation::LynxADDamageDeformation(const InputParameters & parameters)
: LynxADElasticDeformation(parameters),
// Damage model
// _damage_model(&getMaterialByName("damage_model")),
// Strain properties
_elastic_strain(declareADProperty<RankTwoTensor>("elastic_strain")),
_elastic_strain_old(getMaterialPropertyOld<RankTwoTensor>("elastic_strain"))
{
}
void
LynxADDamageDeformation::initialSetup()
{
LynxADDeformationBase::initialSetup();
MaterialName damage_model = getParam<MaterialName>("damage_model");
LynxADDamageModelBase * damage_r =
dynamic_cast<LynxADDamageModelBase *>(&this->getMaterialByName("damage_model"));
dynamic_cast<LynxADDamageModelBase *>(&this->getMaterialByName(damage_model));
_damage_model = damage_r;
_damage_model = damage_r;
}
void
LynxADDamageDeformation::initQpStatefulProperties()
{
LynxADElasticDeformation::initQpStatefulProperties();
_elastic_strain[_qp].zero();
_elastic_strain[_qp].addIa(-_plith[_qp] / (3.0 * averageProperty(_bulk_modulus)));
}
void
......@@ -70,8 +85,8 @@ LynxADDamageDeformation::computeQpStress()
_stress[_qp].addIa(_plith_old[_qp] - _plith[_qp]);
// Damage - plastic correction
_damage_model->damageUpdate(_stress[_qp], elastic_strain_old);
_damage_model->damageUpdate(_stress[_qp], _elastic_strain_incr[_qp]);
_elastic_strain[_qp] = elastic_strain_old + _elastic_strain_incr[_qp];
}
// void
......
/******************************************************************************/
/* This file is part of */
/* LYNX, a MOOSE-based application */
/* Lithosphere dYnamic Numerical toolboX */
/* */
/* Copyright (C) 2017 by Antoine B. Jacquey and Mauro Cacace */
/* GFZ Potsdam, German Research Centre for Geosciences */
/* */
/* Licensed under GNU General Public License 3, */
/* please see LICENSE for details */
/* or http://www.gnu.org/licenses/gpl.html */
/******************************************************************************/
#include "LynxADDamageModelBase.h"
InputParameters
LynxADDamageModelBase::validParams()
{
InputParameters params = LynxADMaterialBase::validParams();
params.addClassDescription("Base class for the damage rheology.");
params.set<bool>("compute") = false;
params.suppressParameter<bool>("compute");
params.addRangeCheckedParam<Real>("abs_tolerance",
1.0e-10,
"abs_tolerance > 0.0",
"The absolute tolerance for the iterative update.");
params.addRangeCheckedParam<Real>("rel_tolerance",
1.0e-10,
"rel_tolerance > 0.0",
"The relative tolerance for the iterative update.");
params.addRangeCheckedParam<unsigned int>(
"max_iterations",
200,
"max_iterations >= 1",
"The maximum number of iterations for the iterative update");
params.addParam<std::vector<Real>>("initial_damage", "The initial damage value.");
return params;
}
LynxADDamageModelBase::LynxADDamageModelBase(const InputParameters & parameters)
: LynxADMaterialBase(parameters),
_abs_tol(getParam<Real>("abs_tolerance")),
_rel_tol(getParam<Real>("rel_tolerance")),
_max_its(getParam<unsigned int>("max_iterations")),
_damage0(isParamValid("initial_damage") ? getLynxParam<Real>("initial_damage")
: std::vector<Real>(_n_composition, 0.0)),
_damage(declareADProperty<Real>("damage")),
_damage_old(getMaterialPropertyOld<Real>("damage")),
_damage_drive(declareADProperty<Real>("damage_force")),
_plastic_strain_incr(declareADProperty<RankTwoTensor>("plastic_strain_increment")),
_damage_incr(declareADProperty<Real>("damage_increment")),
_yield_function(declareADProperty<Real>("plastic_yield_function"))
{
}
void
LynxADDamageModelBase::setQp(unsigned int qp)
{
_qp = qp;
}
void
LynxADDamageModelBase::initQpStatefulProperties()
{
_damage[_qp] = averageProperty(_damage0);
}
\ No newline at end of file
This diff is collapsed.
[Mesh]
type = GeneratedMesh
dim = 2
nx = 5
ny = 5
xmin = 0
xmax = 1
ymin = 0
ymax = 1
[]
[Variables]
[./disp_x]
order = FIRST
family = LAGRANGE
[../]
[./disp_y]
order = FIRST
family = LAGRANGE
[../]
[]
[Kernels]
[./mech_x]
type = LynxADSolidMomentum
variable = disp_x
component = 0
[../]
[./mech_y]
type = LynxADSolidMomentum
variable = disp_y
component = 1
[../]
[]
[AuxVariables]
[./p_lith]
order = FIRST
family = LAGRANGE
initial_condition = 1.0e+01
[../]
[./mises_stress]
order = CONSTANT
family = MONOMIAL
[../]
[./vol_strain]
order = CONSTANT
family = MONOMIAL
[../]
[./eqv_strain]
order = CONSTANT
family = MONOMIAL
[../]
[./eqv_in_strain]
order = CONSTANT
family = MONOMIAL
[../]
[./eqv_in_strain_rate]
order = CONSTANT
family = MONOMIAL
[../]
[./yield]
order = CONSTANT
family = MONOMIAL
[../]
[./damage]
order = CONSTANT
family = MONOMIAL
[../]
[./damage_force]
order = CONSTANT
family = MONOMIAL
[../]
[./pressure]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./p_lith_aux]
type = ConstantAux
variable = p_lith
value = 1.0e+01
[../]
[./mises_stress_aux]
type = LynxADVonMisesStressAux
variable = mises_stress
execute_on = 'TIMESTEP_END'
[../]
[./vol_strain_aux]
type = LynxADVolStrainAux
variable = vol_strain
execute_on = 'TIMESTEP_END'
[../]
[./eqv_strain_aux]
type = LynxADEqvStrainAux
variable = eqv_strain
execute_on = 'TIMESTEP_END'
[../]
[./eqv_in_strain_aux]
type = LynxADEqvStrainAux
variable = eqv_in_strain
strain_type = plastic
execute_on = 'TIMESTEP_END'
[../]
[./eqv_in_strain_rate_aux]
type = LynxADEqvStrainRateAux
variable = eqv_in_strain_rate
strain_type = plastic
execute_on = 'TIMESTEP_END'
[../]
[./yield_aux]
type = ADMaterialRealAux
variable = yield
property = plastic_yield_function
execute_on = 'TIMESTEP_END'
[../]
[./damage_aux]
type = ADMaterialRealAux
variable = damage
property = damage
execute_on = 'TIMESTEP_END'
[../]
[./damage_force_aux]
type = ADMaterialRealAux
variable = damage_force
property = damage_force
execute_on = 'TIMESTEP_END'
[../]
[./pressure_aux]
type = LynxADEffectivePressureAux
variable = pressure
execute_on = 'TIMESTEP_END'
[../]
[]
[BCs]
[./no_ux]
type = DirichletBC
variable = disp_x
boundary = left
value = 0.0
preset = true
[../]
[./ux_right]
type = LynxVelocityBC
variable = disp_x
boundary = right
value = -1.0
[../]
[./uy_top]
type = LynxVelocityBC
variable = disp_y
boundary = top
value = 1.0
[../]
[./no_uy]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
preset = true
[../]
[]
[Materials]
[./elastic_mat]
type = LynxADDamageDeformation
displacements = 'disp_x disp_y'
lithostatic_pressure = 'p_lith'
bulk_modulus = 1.0e+04
shear_modulus = 1.0e+04
damage_model = 'damage_vladi'
[../]
[./damage_vladi]
type = LynxADLyakhovskyDamage
damage_modulus = 1.0e+04
friction_angle = 30
plastic_viscosity = 1.0e-01
initial_damage = 1.0e-03
[../]
[]
[Postprocessors]
[./D]
type = ElementAverageValue
variable = damage
outputs = csv
[../]
[./Se]
type = ElementAverageValue
variable = mises_stress
outputs = csv
[../]
[./E_vol]
type = ElementAverageValue
variable = vol_strain
outputs = csv
[../]
[./E_eqv]
type = ElementAverageValue
variable = eqv_strain
outputs = csv
[../]
[./E_in_eqv]
type = ElementAverageValue
variable = eqv_in_strain
outputs = csv
[../]
[./E_dot_in_eqv]
type = ElementAverageValue
variable = eqv_in_strain_rate
outputs = csv
[../]
[./F]
type = ElementAverageValue
variable = yield
outputs = csv
[../]
[./Dam_force]
type = ElementAverageValue
variable = damage_force
outputs = csv
[../]
[]
[Preconditioning]
[./precond]
type = SMP
full = true
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it
-pc_type -pc_hypre_type -pc_hypre_boomeramg_strong_threshold
-pc_hypre_boomeramg_agg_nl -pc_hypre_boomeramg_agg_num_paths
-pc_hypre_boomeramg_max_levels
-pc_hypre_boomeramg_coarsen_type -pc_hypre_boomeramg_interp_type
-pc_hypre_boomeramg_P_max -pc_hypre_boomeramg_truncfactor
-snes_linesearch_type'
petsc_options_value = 'fgmres 1.0e-06 1.0e-10 1000 100
hypre boomeramg 0.7
4 5
25
HMIS ext+i
2 0.3
basic'
[../]
[]
[Executioner]
type = Transient
solve_type = NEWTON
start_time = 0.0
end_time = 3.0e-03
num_steps = 50
[]
[Outputs]
execute_on = 'timestep_end'
print_linear_residuals = false
perf_graph = true
exodus = true
csv = true
[]
[Mesh]
type = GeneratedMesh
dim = 2
nx = 5
ny = 5
xmin = 0
xmax = 1
ymin = 0
ymax = 1
[]
[Variables]
[./disp_x]