Skip to content

Commit

Permalink
Cleanup ref idaholab#184
Browse files Browse the repository at this point in the history
  • Loading branch information
@vprithiv
vprithiv committed Feb 14, 2023
1 parent 694b323 commit 7c51277
Show file tree
Hide file tree
Showing 5 changed files with 140 additions and 34 deletions.
4 changes: 3 additions & 1 deletion doc/content/source/materials/ComputeMultipleInelasticDamageStress.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,9 @@

## Description

`ComputeMultipleInelasticDamageStress` computes the damage stress.
This class `ComputeMultipleInelasticDamageStress` computes the stress with the damage obtained from
the previous time step. This is done mainly to avoid convergence problems

This ComputeMultipleInelasticStress is to be used with (/DamagePlasticityStressUpdate.md).
## Example Input Files

Expand Down
1 change: 0 additions & 1 deletion include/materials/ComputeMultipleInelasticDamageStress.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,6 @@ class ComputeMultipleInelasticDamageStress : public ComputeMultipleInelasticStre
const MaterialProperty<Real> & _D;
const MaterialProperty<Real> & _D_old;


virtual void computeQpJacobianMult() override;

virtual void computeAdmissibleState(unsigned model_number,
Expand Down
150 changes: 127 additions & 23 deletions include/materials/DamagePlasticityStressUpdate.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,31 +22,41 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate
static InputParameters validParams();
DamagePlasticityStressUpdate(const InputParameters & parameters);

/**
* Does the model require the elasticity tensor to be isotropic?
*/
bool requiresIsotropicTensor() override { return true; }

protected:
virtual void initQpStatefulProperties() override;
virtual void finalizeReturnProcess(const RankTwoTensor & rotation_increment) override;

private:
///Yield function tolerance (user parameter)
const Real _f_tol;
///dimensionless parameter involving biaxial & uniaxial compressive strengths (user parameter)
const Real _alfa;
///dilatancy factor (user parameter)
const Real _alfa_p;
///stiffness recovery factor (user parameter)
const Real _s0;

///ft_ep_slope_factor_at_zero_ep (user parameter)
const Real _Chi;
///tensile damage at half tensile strength (user parameter)
const Real _Dt;
///yield strength in tension (user parameter)
const Real _ft;
///fracture_energy in tension (user parameter)
const Real _FEt;

///yield strength in compression (user parameter)
const Real _fyc;
///compressive damage at maximum compressive strength (user parameter)
const Real _Dc;
///maximum strength in compression (user parameter)
const Real _fc;
///fracture energy in compression (user parameter)
const Real _FEc;

///@{
/** The following variables are intermediate and are calculated based on the user parameters given
* above */
const Real _at;
const Real _ac;
const Real _zt;
Expand All @@ -59,6 +69,9 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate
const Real _dc_bc;
const Real _ft0;
const Real _fc0;
///@}

/// Intermediate variable calcualted using user parameter tip_smoother
const Real _small_smoother2;

const Real _sqrt3;
Expand All @@ -68,37 +81,90 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate

/// Eigenvectors of the trial stress as a RankTwoTensor, in order to rotate the returned stress back to stress space
RankTwoTensor _eigvecs;


///damage state in tension
MaterialProperty<Real> & _intnl0;
///damage state in compression
MaterialProperty<Real> & _intnl1;
///element length
MaterialProperty<Real> & _ele_len;
///fracture energy in tension
MaterialProperty<Real> & _gt;
///fracture energy in compression
MaterialProperty<Real> & _gc;

///tensile damage
MaterialProperty<Real> & _tD;
///compression damage
MaterialProperty<Real> & _cD;
///damage variable
MaterialProperty<Real> & _D;
///minimum plastic strain
MaterialProperty<Real> & _min_ep;
///mid value of plastic strain
MaterialProperty<Real> & _mid_ep;
///maximum plastic strain
MaterialProperty<Real> & _max_ep;
///damaged minimum principal stress
MaterialProperty<Real> & _sigma0;
///damaged mid value of principal stress
MaterialProperty<Real> & _sigma1;
///damaged maximum principal stress
MaterialProperty<Real> & _sigma2;

Real ft(const std::vector<Real> & intnl) const; /// tensile strength
Real dft(const std::vector<Real> & intnl) const; /// d(tensile strength)/d(intnl)
Real fc(const std::vector<Real> & intnl) const; /// compressive strength
Real dfc(const std::vector<Real> & intnl) const; /// d(compressive strength)/d(intnl)
Real ft(const std::vector<Real> & intnl) const;
/**
* Obtain the tensile strength
* @param intnl (Array containing damage states in tension and compression, respectively)
*/
Real dft(const std::vector<Real> & intnl) const;
/**
* Obtain the partial derivative of the tensile strength to the damage state
* @param intnl (Array containing damage states in tension and compression, respectively)
*/
Real fc(const std::vector<Real> & intnl) const;
/**
* Obtain the conpressive strength
* @param intnl (Array containing damage states in tension and compression, respectively)
*/
Real dfc(const std::vector<Real> & intnl) const;
/**
* Obtain the partial derivative of the compressive strength to the damage state
* @param intnl (Array containing damage states in tension and compression, respectively)
*/
Real beta(const std::vector<Real> & intnl) const;
/**
* beta is a dimensionless constant, which is a component of the yield function
* It is defined in terms of tensile strength, compressive strength, and another
* dimensionless constant alpha (See Eqn. 37 in Lee (1998))
* @param intnl (Array containing damage states in tension and compression, respectively)
*/
Real dbeta0(const std::vector<Real> & intnl) const;
/**
* dbeta0 is a derivative of beta wrt. tensile strength (ft)
* @param intnl (Array containing damage states in tension and compression, respectively)
*/
Real dbeta1(const std::vector<Real> & intnl) const;
void weighfac(const std::vector<Real> & stress_params, Real & wf) const; /// weight factor
void dweighfac(const std::vector<Real> & stress_params,
Real & wf,
std::vector<Real> & dwf) const; /// d(weight factor)/d(stress)
/**
* dbeta1 is a derivative of beta wrt. compressive strength (fc)
* @param intnl (Array containing damage states in tension and compression, respectively)
*/
void weighfac(const std::vector<Real> & stress_params, Real & wf) const;
/**
* weighfac is the weight factor, defined interms of the ratio of sum of principal stresses
* to the sum of absolute value of the principal stresses
* @param stress_params is the array of principal stresses
*/
void dweighfac(const std::vector<Real> & stress_params, Real & wf, std::vector<Real> & dwf) const;
/**
* dweighfac is the derivative of the weight factor
* @param stress_params is the array of principal stresses
*/
Real damageVar(const std::vector<Real> & stress_params, const std::vector<Real> & intnl) const;

/**
* damageVar is the degradation damage variable as defined in Eqn. 43 in Lee (1998)
* @param stress_params is the array of principal stresses
* @param intnl (Array containing damage states in tension and compression, respectively)
* @param r
*/
void computeStressParams(const RankTwoTensor & stress,
std::vector<Real> & stress_params) const override;

Expand Down Expand Up @@ -133,33 +199,71 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate
virtual void flowPotential(const std::vector<Real> & stress_params,
const std::vector<Real> & intnl,
std::vector<Real> & r) const;
/**
* This function calculates the flow potential
* @param stress_params is the array of principal stresses
* @param intnl (Array containing damage states in tension and compression, respectively)
* @param r is the flowpotential
*/

virtual void dflowPotential_dstress(const std::vector<Real> & stress_params,
const std::vector<Real> & intnl,
std::vector<std::vector<Real>> & dr_dstress) const;

/**
* This function calculates the derivative of the flow potential with the stress
* @param stress_params is the array of principal stresses
* @param intnl (Array containing damage states in tension and compression, respectively)
* @param dr_dstress is the dflowpotential
*/
virtual void dflowPotential_dintnl(const std::vector<Real> & stress_params,
const std::vector<Real> & intnl,
std::vector<std::vector<Real>> & dr_dintnl) const;

/**
* This function calculates the derivative of the flow potential with the damage states
* @param stress_params is the array of principal stresses
* @param intnl (Array containing damage states in tension and compression, respectively)
* @param dr_dintnl is the dflowPotential_dintnl
*/
virtual void hardPotential(const std::vector<Real> & stress_params,
const std::vector<Real> & intnl,
std::vector<Real> & h) const;

/**
* This function calculates the hardening potential
* @param stress_params is the array of principal stresses
* @param intnl (Array containing damage states in tension and compression, respectively)
* @param h is the hardPotential
*/
virtual void dhardPotential_dstress(const std::vector<Real> & stress_params,
const std::vector<Real> & intnl,
std::vector<std::vector<Real>> & dh_dsig) const;

/**
* This function calculates the derivative of the hardening potential with the stress
* @param stress_params is the array of principal stresses
* @param intnl (Array containing damage states in tension and compression, respectively)
* @param dh_dsig is the dhardPotential_dstress
*/
virtual void dhardPotential_dintnl(const std::vector<Real> & stress_params,
const std::vector<Real> & intnl,
std::vector<std::vector<Real>> & dh_dintnl) const;

/**
* This function calculates the derivative of the hardening potential with the damage states
* @param stress_params is the array of principal stresses
* @param intnl (Array containing damage states in tension and compression, respectively)
* @param dh_dintnl is the dhardPotential_dintnl
*/
void initialiseVarsV(const std::vector<Real> & trial_stress_params,
const std::vector<Real> & intnl_old,
std::vector<Real> & stress_params,
Real & gaE,
std::vector<Real> & intnl) const;

/**
* This function updates the damage states
* @param stress_params is the array of principal stresses
* @param trial_stress_params is the trial values of the principal stresses
* @param intnl_old (Array containing previous step's damage states in tension and compression,
* respectively)
* @param dh_dintnl is the dhardPotential_dintnl
*/
void setIntnlValuesV(const std::vector<Real> & trial_stress_params,
const std::vector<Real> & current_stress_params,
const std::vector<Real> & intnl_old,
Expand Down
4 changes: 1 addition & 3 deletions src/materials/ComputeMultipleInelasticDamageStress.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ ComputeMultipleInelasticDamageStress::validParams()
{
InputParameters params = ComputeMultipleInelasticStress::validParams();
params.addClassDescription("This ComputeMultipleInelasticStress is to be used with "
"DamagePlasticityStressUpdate");
"DamagePlasticityStressUpdate");
return params;
}

Expand All @@ -39,7 +39,6 @@ ComputeMultipleInelasticDamageStress::computeQpJacobianMult()
{
ComputeMultipleInelasticStress::computeQpJacobianMult();
_Jacobian_mult[_qp] = (1.0 - _D_old[_qp]) * _Jacobian_mult[_qp];
// _Jacobian_mult[_qp] = (1.0 - _D[_qp]) * _Jacobian_mult[_qp];
}

void
Expand All @@ -65,7 +64,6 @@ ComputeMultipleInelasticDamageStress::computeAdmissibleState(
_rotation_increment[_qp],
_stress[_qp],
_stress_old[_qp] / (1.0 - _D_old[_qp]),
// _stress_old[_qp] / (1.0 - _D[_qp]),
_elasticity_tensor[_qp],
_elastic_strain_old[_qp],
_tangent_operator_type == TangentOperatorEnum::nonlinear,
Expand Down
15 changes: 9 additions & 6 deletions src/materials/DamagePlasticityStressUpdate.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,8 +94,7 @@ DamagePlasticityStressUpdate::DamagePlasticityStressUpdate(const InputParameters
_FEc(getParam<Real>("fracture_energy_in_compression")),

_at(1.5 * std::sqrt(1 - _Chi) - 0.5),
_ac((2. * (_fc / _fyc) - 1. + 2. * std::sqrt(std::pow((_fc / _fyc), 2.) - _fc / _fyc))),

_ac((2. * (_fc / _fyc) - 1. + 2. * std::sqrt(Utility::pow<2>(_fc / _fyc) - _fc / _fyc))),
_zt((1. + _at) / _at),
_zc((1. + _ac) / _ac),
_dPhit(_at * (2. + _at)),
Expand All @@ -108,7 +107,6 @@ DamagePlasticityStressUpdate::DamagePlasticityStressUpdate(const InputParameters
((1. - _Dt) * std::pow((_zt - _sqrtPhit_max / _at), (1. - _dt_bt)) * _sqrtPhit_max)),
_fc0(_fc / ((1. - _Dc) * std::pow((_zc - _sqrtPhic_max / _ac), (1. - _dc_bc)) * _sqrtPhic_max)),
_small_smoother2(Utility::pow<2>(getParam<Real>("tip_smoother"))),


_sqrt3(std::sqrt(3.)),
_perfect_guess(getParam<bool>("perfect_guess")),
Expand All @@ -133,6 +131,11 @@ DamagePlasticityStressUpdate::DamagePlasticityStressUpdate(const InputParameters
void
DamagePlasticityStressUpdate::initQpStatefulProperties()
{
/* There are multiple ways to determine the correct element length (or the characteristic length)
used, the following commented lines show several different options. Some other options are
still being considered. In this code, we define the element length as the cube root of the
element volume */

// if (_current_elem->n_vertices() < 3)
// _ele_len[_qp] = _current_elem->length(0, 1);
// else if (_current_elem->n_vertices() < 5)
Expand Down Expand Up @@ -318,7 +321,7 @@ DamagePlasticityStressUpdate::flowPotential(const std::vector<Real> & stress_par
Real D = damageVar(stress_params, intnl);

for (unsigned int i = 0; i < _num_sp; ++i)
r[i] = (_alfa_p + d_sqrt_2J2(i, i)) * std::pow((1. - D), 1);
r[i] = (_alfa_p + d_sqrt_2J2(i, i)) * (1. - D);
}

void
Expand Down Expand Up @@ -360,7 +363,7 @@ DamagePlasticityStressUpdate::dflowPotential_dstress(
for (unsigned i = 0; i < _num_sp; ++i)
for (unsigned j = 0; j < (i + 1); ++j)
{
dr_dstress[i][i] = J2 < _f_tol ? 0. : dfp(i, i, j, j) * std::pow((1. - D), 2);
dr_dstress[i][i] = J2 < _f_tol ? 0. : dfp(i, i, j, j) * Utility::pow<2>(1. - D);
if (i != j)
dr_dstress[j][i] = dr_dstress[i][j];
}
Expand Down Expand Up @@ -582,7 +585,7 @@ DamagePlasticityStressUpdate::beta(const std::vector<Real> & intnl) const
Real
DamagePlasticityStressUpdate::dbeta0(const std::vector<Real> & intnl) const
{
return -(1. - _alfa) * fc(intnl) * dft(intnl) / std::pow(ft(intnl), 2.);
return -(1. - _alfa) * fc(intnl) * dft(intnl) / Utility::pow<2>(ft(intnl));
}

Real
Expand Down

0 comments on commit 7c51277

Please sign in to comment.