Commit 6c336965 authored by Antoine Jacquey's avatar Antoine Jacquey

Capped yield working.

parent 1053a16c
......@@ -36,8 +36,10 @@ public:
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 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);
......@@ -45,6 +47,7 @@ protected:
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);
......@@ -52,6 +55,7 @@ protected:
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;
......
......@@ -55,6 +55,7 @@ protected:
virtual void deviatoricDeformation(const Real & pressure, RankTwoTensor & stress_dev);
virtual void reformStressTensor(const Real & pressure, const RankTwoTensor & stress_dev);
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);
......
......@@ -109,14 +109,32 @@ LynxDamageDeformation::initQpStatefulProperties()
LynxDeformationBase::initQpStatefulProperties();
}
void
LynxDamageDeformation::initializeQpDeformation()
{
// Need to get _xi0
initializeDamageParameters();
LynxDeformationBase::initializeQpDeformation();
}
void
LynxDamageDeformation::elasticModuli()
{
// Strain ratio
const Real xi = strainRatio(_elastic_strain_old[_qp]);
// Bulk modulus
_K[_qp] = _has_bulk_modulus ? averageProperty(_bulk_modulus) : 0.0;
_K[_qp] = _has_bulk_modulus ? averageProperty(_bulk_modulus) -
2.0 / 3.0 * _damage_old[_qp] * _damage_plasticity->_gamma *
(xi / 2.0 - _damage_plasticity->_xi0)
: 0.0;
// Shear modulus
_G[_qp] = _has_shear_modulus ? averageProperty(_shear_modulus) : 0.0;
_G[_qp] = _has_shear_modulus
? averageProperty(_shear_modulus) - _damage_old[_qp] * _damage_plasticity->_gamma *
(xi / 2.0 - _damage_plasticity->_xi0)
: 0.0;
}
void
......@@ -169,6 +187,18 @@ LynxDamageDeformation::plasticCorrection(Real & pressure, RankTwoTensor & stress
rho_0;
}
}
// Update elastic strain
_elastic_strain[_qp] = spinRotation(_elastic_strain_old[_qp].deviatoric());
_elastic_strain[_qp].addIa(_elastic_strain_old[_qp].trace() / 3.0);
_elastic_strain[_qp] +=
_strain_increment[_qp] - _viscous_strain_incr[_qp] - _plastic_strain_incr[_qp];
}
void
LynxDamageDeformation::damageCorrection()
{
}
Real
......@@ -330,6 +360,23 @@ LynxDamageDeformation::updateDamageParameters()
averageProperty(_one_on_damage_eta));
}
void
LynxDamageDeformation::initializeDamageParameters()
{
Real sin_phi = std::sin(averageProperty(_friction_angle) * libMesh::pi / 180.0);
Real xi0 = -std::sqrt(3.0) /
std::sqrt(1.0 + 1.5 * Utility::pow<2>(averageProperty(_bulk_modulus) /
averageProperty(_shear_modulus) * sin_phi));
Real porous_coeff =
averageProperty(_porous_coeff) + averageProperty(_porous_coeff_linear) * _porosity[_qp];
Real p_cr = _K[_qp] * (std::sqrt(3.0) + xi0) / (3.0 * porous_coeff);
_damage_plasticity->fill(xi0,
averageProperty(_damage_modulus),
p_cr,
averageProperty(_one_on_plastic_eta),
averageProperty(_one_on_damage_eta));
}
void
LynxDamageDeformation::updateDamageConvexParameters(const Real & pressure, const Real & eqv_stress)
{
......@@ -519,3 +566,13 @@ LynxDamageDeformation::strainRatio(const RankTwoTensor & elastic_strain)
else
return -std::sqrt(3.0);
}
RankTwoTensor
LynxDamageDeformation::rotatedElasticStrain(const RankTwoTensor & elastic_strain)
{
const Real strain_v = elastic_strain.trace();
RankTwoTensor strain_rot = spinRotation(elastic_strain.deviatoric());
strain_rot.addIa(strain_v / 3.0);
return strain_rot;
}
......@@ -340,6 +340,7 @@ LynxDeformationBase::computeQpDeformation()
reformStressTensor(pressure, stress_dev);
// Additional correction
damageCorrection();
// Update tangent operator modulus
if (_fe_problem.currentlyComputingJacobian())
......@@ -426,6 +427,11 @@ LynxDeformationBase::reformStressTensor(const Real & pressure, const RankTwoTens
_stress[_qp].addIa(-pressure);
}
void
LynxDeformationBase::damageCorrection()
{
}
void
LynxDeformationBase::tangentOperator()
{
......
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