From 7c51277af5caede10ae9a43cab1fcb415be8ef34 Mon Sep 17 00:00:00 2001 From: Prithivirajan Veerappan Date: Wed, 1 Feb 2023 15:47:52 -0700 Subject: [PATCH] Cleanup ref #184 --- .../ComputeMultipleInelasticDamageStress.md | 4 +- .../ComputeMultipleInelasticDamageStress.h | 1 - .../materials/DamagePlasticityStressUpdate.h | 150 +++++++++++++++--- .../ComputeMultipleInelasticDamageStress.C | 4 +- src/materials/DamagePlasticityStressUpdate.C | 15 +- 5 files changed, 140 insertions(+), 34 deletions(-) diff --git a/doc/content/source/materials/ComputeMultipleInelasticDamageStress.md b/doc/content/source/materials/ComputeMultipleInelasticDamageStress.md index c7166ced..6fef69c6 100644 --- a/doc/content/source/materials/ComputeMultipleInelasticDamageStress.md +++ b/doc/content/source/materials/ComputeMultipleInelasticDamageStress.md @@ -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 diff --git a/include/materials/ComputeMultipleInelasticDamageStress.h b/include/materials/ComputeMultipleInelasticDamageStress.h index 00e2e88e..dae6d1c7 100644 --- a/include/materials/ComputeMultipleInelasticDamageStress.h +++ b/include/materials/ComputeMultipleInelasticDamageStress.h @@ -30,7 +30,6 @@ class ComputeMultipleInelasticDamageStress : public ComputeMultipleInelasticStre const MaterialProperty & _D; const MaterialProperty & _D_old; - virtual void computeQpJacobianMult() override; virtual void computeAdmissibleState(unsigned model_number, diff --git a/include/materials/DamagePlasticityStressUpdate.h b/include/materials/DamagePlasticityStressUpdate.h index e10639be..04d64252 100644 --- a/include/materials/DamagePlasticityStressUpdate.h +++ b/include/materials/DamagePlasticityStressUpdate.h @@ -22,9 +22,6 @@ 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: @@ -32,21 +29,34 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate 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; @@ -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; @@ -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 & _intnl0; + ///damage state in compression MaterialProperty & _intnl1; + ///element length MaterialProperty & _ele_len; + ///fracture energy in tension MaterialProperty & _gt; + ///fracture energy in compression MaterialProperty & _gc; - + ///tensile damage MaterialProperty & _tD; + ///compression damage MaterialProperty & _cD; + ///damage variable MaterialProperty & _D; + ///minimum plastic strain MaterialProperty & _min_ep; + ///mid value of plastic strain MaterialProperty & _mid_ep; + ///maximum plastic strain MaterialProperty & _max_ep; + ///damaged minimum principal stress MaterialProperty & _sigma0; + ///damaged mid value of principal stress MaterialProperty & _sigma1; + ///damaged maximum principal stress MaterialProperty & _sigma2; - Real ft(const std::vector & intnl) const; /// tensile strength - Real dft(const std::vector & intnl) const; /// d(tensile strength)/d(intnl) - Real fc(const std::vector & intnl) const; /// compressive strength - Real dfc(const std::vector & intnl) const; /// d(compressive strength)/d(intnl) + Real ft(const std::vector & intnl) const; + /** + * Obtain the tensile strength + * @param intnl (Array containing damage states in tension and compression, respectively) + */ + Real dft(const std::vector & 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 & intnl) const; + /** + * Obtain the conpressive strength + * @param intnl (Array containing damage states in tension and compression, respectively) + */ + Real dfc(const std::vector & 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 & 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 & 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 & intnl) const; - void weighfac(const std::vector & stress_params, Real & wf) const; /// weight factor - void dweighfac(const std::vector & stress_params, - Real & wf, - std::vector & 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 & 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 & stress_params, Real & wf, std::vector & dwf) const; + /** + * dweighfac is the derivative of the weight factor + * @param stress_params is the array of principal stresses + */ Real damageVar(const std::vector & stress_params, const std::vector & 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 & stress_params) const override; @@ -133,33 +199,71 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate virtual void flowPotential(const std::vector & stress_params, const std::vector & intnl, std::vector & 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 & stress_params, const std::vector & intnl, std::vector> & 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 & stress_params, const std::vector & intnl, std::vector> & 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 & stress_params, const std::vector & intnl, std::vector & 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 & stress_params, const std::vector & intnl, std::vector> & 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 & stress_params, const std::vector & intnl, std::vector> & 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 & trial_stress_params, const std::vector & intnl_old, std::vector & stress_params, Real & gaE, std::vector & 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 & trial_stress_params, const std::vector & current_stress_params, const std::vector & intnl_old, diff --git a/src/materials/ComputeMultipleInelasticDamageStress.C b/src/materials/ComputeMultipleInelasticDamageStress.C index a102bad3..08b68099 100644 --- a/src/materials/ComputeMultipleInelasticDamageStress.C +++ b/src/materials/ComputeMultipleInelasticDamageStress.C @@ -22,7 +22,7 @@ ComputeMultipleInelasticDamageStress::validParams() { InputParameters params = ComputeMultipleInelasticStress::validParams(); params.addClassDescription("This ComputeMultipleInelasticStress is to be used with " - "DamagePlasticityStressUpdate"); + "DamagePlasticityStressUpdate"); return params; } @@ -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 @@ -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, diff --git a/src/materials/DamagePlasticityStressUpdate.C b/src/materials/DamagePlasticityStressUpdate.C index a4797fd0..bd06e05a 100644 --- a/src/materials/DamagePlasticityStressUpdate.C +++ b/src/materials/DamagePlasticityStressUpdate.C @@ -94,8 +94,7 @@ DamagePlasticityStressUpdate::DamagePlasticityStressUpdate(const InputParameters _FEc(getParam("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)), @@ -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("tip_smoother"))), - _sqrt3(std::sqrt(3.)), _perfect_guess(getParam("perfect_guess")), @@ -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) @@ -318,7 +321,7 @@ DamagePlasticityStressUpdate::flowPotential(const std::vector & 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 @@ -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]; } @@ -582,7 +585,7 @@ DamagePlasticityStressUpdate::beta(const std::vector & intnl) const Real DamagePlasticityStressUpdate::dbeta0(const std::vector & 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