diff --git a/include/utils/LynxRheologyStructures.h b/include/utils/LynxRheologyStructures.h index 7b8bf28619822fd2f6b78fe5922bd06f846225e5..024b80063ba23592d737d6967befece6700b4fc3 100644 --- a/include/utils/LynxRheologyStructures.h +++ b/include/utils/LynxRheologyStructures.h @@ -227,7 +227,7 @@ struct damage_plasticity Real _alpha2; // square of alpha parameter for capped yield Real _dxi_cr_dp; // derivative wrt pressure of the critical strain ratio Real _dxi_cr_dq; // derivative wrt eqv_stress of the critical strain ratio - Real _dmu2_dxi_cr; // derivative wrt the critical strain ratio of the square of alpha parameter + Real _dalpha2_dxi_cr; // derivative wrt the critical strain ratio of the square of alpha parameter damage_plasticity() : _xi0(-std::sqrt(3.0)), _gamma(0.0), @@ -248,7 +248,7 @@ struct damage_plasticity _alpha2(0.0), _dxi_cr_dp(0.0), _dxi_cr_dq(0.0), - _dmu2_dxi_cr(0.0) + _dalpha2_dxi_cr(0.0) { } void fill(const Real xi0, diff --git a/src/materials/LynxDamageDeformation.C b/src/materials/LynxDamageDeformation.C index 1a5d6de4cb0590f0cee0697e2da70b19b99f16d7..ec9c49dcc423ed03c32a5938d59cfa901cd2f861 100644 --- a/src/materials/LynxDamageDeformation.C +++ b/src/materials/LynxDamageDeformation.C @@ -327,8 +327,8 @@ LynxDamageDeformation::convexPlasticIncrement(Real & vol_plastic_incr, Real & eq { _stress_corr_p = (_damage_plasticity->_q_tr != 0.0) - ? _damage_plasticity->_q_tr - - 3.0 * _G[_qp] * eqv_plastic_incr / _damage_plasticity->_q_tr + ? (_damage_plasticity->_q_tr - + 3.0 * _G[_qp] * eqv_plastic_incr) / _damage_plasticity->_q_tr : 1.0; Real dF2_dp0 = dConvexPlasticYield2_dp(p_0, q_0) / dConvexPlasticYield2(rho_0); Real dF2_dq0 = dConvexPlasticYield2_dq(p_0, q_0) / dConvexPlasticYield2(rho_0); @@ -386,10 +386,10 @@ LynxDamageDeformation::computeDamageProperties(const Real & pressure, const Real Utility::pow<2>(eqv_stress - _damage_plasticity->_q_r)); // Projection direction - _damage_plasticity->_rp = - (_damage_plasticity->_p_tr - _damage_plasticity->_p_r) / _damage_plasticity->_rho_tr; - _damage_plasticity->_rq = - (_damage_plasticity->_q_tr - _damage_plasticity->_q_r) / _damage_plasticity->_rho_tr; + _damage_plasticity->_rp = (_damage_plasticity->_rho_tr != 0.0) ? + (_damage_plasticity->_p_tr - _damage_plasticity->_p_r) / _damage_plasticity->_rho_tr : 0.0; + _damage_plasticity->_rq = (_damage_plasticity->_rho_tr != 0.0) ? + (_damage_plasticity->_q_tr - _damage_plasticity->_q_r) / _damage_plasticity->_rho_tr : 0.0; } } @@ -401,6 +401,7 @@ LynxDamageDeformation::updateDamageParameters() Real porous_coeff = averageProperty(_porous_coeff) + averageProperty(_porous_coeff_linear) * _porosity[_qp]; Real p_cr = (porous_coeff != 0.0) ? _K[_qp] * (std::sqrt(3.0) + xi0) / (3.0 * porous_coeff) : 0.0; + Real k0 = std::sqrt(3.0) * averageProperty(_cohesion) * std::cos(averageProperty(_friction_angle) * libMesh::pi / 180.0); @@ -441,7 +442,7 @@ LynxDamageDeformation::updateDamageConvexParameters(const Real & pressure, const eqv_stress * (1.0 - pq_term) : 0.0; - _damage_plasticity->_dmu2_dxi_cr = + _damage_plasticity->_dalpha2_dxi_cr = (_damage_plasticity->_xi_cr != 0.0) ? -12.0 * Utility::pow<2>(_G[_qp] / _K[_qp]) / Utility::pow<3>(_damage_plasticity->_xi_cr) : 0.0; @@ -450,16 +451,23 @@ LynxDamageDeformation::updateDamageConvexParameters(const Real & pressure, const Real LynxDamageDeformation::convexReferencePressure() { - Real p_proj = - (1.0 - Utility::pow<2>(_damage_plasticity->_alpha0) * _K[_qp] / - (3.0 * _G[_qp] + Utility::pow<2>(_damage_plasticity->_alpha0) * _K[_qp])) * - _damage_plasticity->_p_tr + - _damage_plasticity->_alpha0 * _K[_qp] * _damage_plasticity->_q_tr / - (3.0 * _G[_qp] + Utility::pow<2>(_damage_plasticity->_alpha0) * _K[_qp]); - - return ((p_proj > 0.0) && (p_proj < _damage_plasticity->_p_cr)) - ? p_proj - : 2.0 * _damage_plasticity->_p_cr / 3.0; + // Real p_proj = + // (1.0 - Utility::pow<2>(_damage_plasticity->_alpha0) * _K[_qp] / + // (3.0 * _G[_qp] + Utility::pow<2>(_damage_plasticity->_alpha0) * _K[_qp])) * + // _damage_plasticity->_p_tr + + // _damage_plasticity->_alpha0 * _K[_qp] * _damage_plasticity->_q_tr / + // (3.0 * _G[_qp] + Utility::pow<2>(_damage_plasticity->_alpha0) * _K[_qp]); + + Real p_proj = _damage_plasticity->_q_tr / _damage_plasticity->_alpha0 - _damage_plasticity->_p_k; + + if (_damage_plasticity->_q_tr == 0.0) + return _damage_plasticity->_p_tr; + else if (p_proj < -_damage_plasticity->_p_k) + return -_damage_plasticity->_p_k; + else if (p_proj > _damage_plasticity->_p_cr) + return 2.0 * _damage_plasticity->_p_cr / 3.0; + else + return p_proj; } Real @@ -480,10 +488,10 @@ LynxDamageDeformation::dConvexPlasticYield2_dp(const Real & pressure, const Real updateDamageConvexParameters(pressure, eqv_stress); - return -(_damage_plasticity->_dmu2_dxi_cr * _damage_plasticity->_dxi_cr_dp * + return -(_damage_plasticity->_dalpha2_dxi_cr * _damage_plasticity->_dxi_cr_dp * (pressure + _damage_plasticity->_p_k) + 2.0 * _damage_plasticity->_alpha2) * - pressure * + (pressure + _damage_plasticity->_p_k) * ((0.0 < pressure + _damage_plasticity->_p_k) - (pressure + _damage_plasticity->_p_k < 0.0)); } @@ -497,7 +505,7 @@ LynxDamageDeformation::dConvexPlasticYield2_dq(const Real & pressure, const Real updateDamageConvexParameters(pressure, eqv_stress); return 2.0 * eqv_stress - - _damage_plasticity->_dmu2_dxi_cr * _damage_plasticity->_dxi_cr_dq * + _damage_plasticity->_dalpha2_dxi_cr * _damage_plasticity->_dxi_cr_dq * Utility::pow<2>( pressure + _damage_plasticity->_p_k); // * diff --git a/test/tests/plastic/gold/convex_compression_out.e b/test/tests/plastic/gold/convex_compression_out.e index d1c6acc51b20bd79f6683a7f09469e897cd22ea6..407e5cecb0a265e7bf3eb7c19df5fd99cb56b91f 100644 Binary files a/test/tests/plastic/gold/convex_compression_out.e and b/test/tests/plastic/gold/convex_compression_out.e differ diff --git a/test/tests/plastic/gold/convex_oedometer_out.e b/test/tests/plastic/gold/convex_oedometer_out.e index 34bed275e312d3de8642efad7e8084b6fc81ef52..ec8c14c044416ff2432709c57490d57b6511fe0d 100644 Binary files a/test/tests/plastic/gold/convex_oedometer_out.e and b/test/tests/plastic/gold/convex_oedometer_out.e differ diff --git a/test/tests/plastic/gold/convex_oedometer_vp_out.e b/test/tests/plastic/gold/convex_oedometer_vp_out.e index 2306111eeb6f55a46501e258171b5ee3a31c8084..053c911031a51f48ca25780abf251d20c2be96c4 100644 Binary files a/test/tests/plastic/gold/convex_oedometer_vp_out.e and b/test/tests/plastic/gold/convex_oedometer_vp_out.e differ