From 3eecc9b26174a986087b26f8411efff4d6d5b0db Mon Sep 17 00:00:00 2001 From: Prithivirajan Veerappan Date: Sun, 15 Oct 2023 20:41:35 -0600 Subject: [PATCH 1/7] Major Changes : Constants (ft0,fc0,at), Functions(f,fbar,df,dfbar), WIP --- .../materials/DamagePlasticityStressUpdate.h | 65 ++++++++-- src/materials/DamagePlasticityStressUpdate.C | 121 +++++++++++++++--- 2 files changed, 159 insertions(+), 27 deletions(-) diff --git a/include/materials/DamagePlasticityStressUpdate.h b/include/materials/DamagePlasticityStressUpdate.h index e8b47c5a..f454a6e1 100644 --- a/include/materials/DamagePlasticityStressUpdate.h +++ b/include/materials/DamagePlasticityStressUpdate.h @@ -38,13 +38,15 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate ///stiffness recovery factor (user parameter) const Real _s0; ///ft_ep_slope_factor_at_zero_ep (user parameter) - const Real _Chi; + const Real _dstress_dep; ///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; + ///Area under the stress and plastic strain in uniaxial tension (user parameter) + const Real _Area_S_ep; ///yield strength in compression (user parameter) const Real _fyc; ///compressive damage at maximum compressive strength (user parameter) @@ -53,6 +55,10 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate const Real _fc; ///fracture energy in compression (user parameter) const Real _FEc; + /// Maximum stress in tension without damage + const Real _ft0; + /// Maximum stress in compression without damage + const Real _fc0; ///@{ /** The following variables are intermediate and are calculated based on the user parameters given @@ -67,8 +73,7 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate const Real _sqrtPhic_max; const Real _dt_bt; const Real _dc_bc; - const Real _ft0; - const Real _fc0; + ///@} /// Intermediate variable calculated using user parameter tip_smoother @@ -111,33 +116,77 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate ///damaged maximum principal stress MaterialProperty & _sigma2; /** - * Obtain the tensile strength + * Obtain the undamaged tensile strength + * @param intnl (Array containing damage states in tension and compression, respectively) + * @return value of ft (tensile strength) + */ + Real ftbar(const std::vector & intnl) const; + /** + * Obtain the partial derivative of the undamaged tensile strength to the damage state + * @param intnl (Array containing damage states in tension and compression, respectively) + * @return value of dft (partial derivative of the tensile strength to the damage state) + */ + Real dfbar_dkappa(const std::vector & f0, + const std::vector & a, + const std::vector & exponent, + const std::vector & kappa) const; + + /** + * Obtain the partial derivative of the undamaged tensile strength to the damage state + * @param intnl (Array containing damage states in tension and compression, respectively) + * @return value of dft (partial derivative of the tensile strength to the damage state) + */ + Real dftbar(const std::vector & intnl) const; + + /** + * Obtain the undamaged conpressive strength + * @param intnl (Array containing damage states in tension and compression, respectively) + * @return value of fc (compressive strength) + */ + Real fcbar(const std::vector & intnl) const; + + /** + * Obtain the partial derivative of the undamaged compressive strength to the damage state + * @param intnl (Array containing damage states in tension and compression, respectively) + * @return value of dfc + */ + Real dfcbar(const std::vector & intnl) const; + /** + * Obtain the damaged tensile strength * @param intnl (Array containing damage states in tension and compression, respectively) * @return value of ft (tensile strength) */ Real ft(const std::vector & intnl) const; /** - * Obtain the partial derivative of the tensile strength to the damage state + * Obtain the partial derivative of the damaged tensile strength to the damage state * @param intnl (Array containing damage states in tension and compression, respectively) * @return value of dft (partial derivative of the tensile strength to the damage state) */ Real dft(const std::vector & intnl) const; /** - * Obtain the conpressive strength + * Obtain the damaged compressive strength * @param intnl (Array containing damage states in tension and compression, respectively) - * @return value of fc (conpressive strength) + * @return value of fc (compressive strength) */ Real fc(const std::vector & intnl) const; /** - * Obtain the partial derivative of the compressive strength to the damage state + * Obtain the partial derivative of the damaged compressive strength to the damage state * @param intnl (Array containing damage states in tension and compression, respectively) * @return value of dfc */ Real dfc(const std::vector & intnl) const; + /** + * Obtain the partial derivative of the undamaged strength to the damage state + * @param intnl (Array containing damage states in tension and compression, respectively) + * @return value of dft (partial derivative of the tensile strength to the damage state) + */ + Real df_dkappa(const std::vector & f0, + const std::vector & exponent, + const std::vector & kappa) 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 diff --git a/src/materials/DamagePlasticityStressUpdate.C b/src/materials/DamagePlasticityStressUpdate.C index d9006449..28756a90 100644 --- a/src/materials/DamagePlasticityStressUpdate.C +++ b/src/materials/DamagePlasticityStressUpdate.C @@ -48,6 +48,10 @@ DamagePlasticityStressUpdate::validParams() params.addRangeCheckedParam("fracture_energy_in_tension", "fracture_energy_in_tension >= 0", "Fracture energy of concrete in uniaxial tension"); + params.addRangeCheckedParam( + "Area_under_stress_ep", + "Area_under_stress_ep >= 0", + "Area under the stress and plastic strain in uniaxial tension"); params.addRangeCheckedParam("yield_strength_in_compression", "yield_strength_in_compression >= 0", @@ -83,17 +87,20 @@ DamagePlasticityStressUpdate::DamagePlasticityStressUpdate(const InputParameters _alfa_p(getParam("dilatancy_factor")), _s0(getParam("stiff_recovery_factor")), - _Chi(getParam("ft_ep_slope_factor_at_zero_ep")), + _dstress_dep(getParam("ft_ep_slope_factor_at_zero_ep")), _Dt(getParam("tensile_damage_at_half_tensile_strength")), _ft(getParam("yield_strength_in_tension")), _FEt(getParam("fracture_energy_in_tension")), + _Area_S_ep(getParam("Area_under_stress_ep")), _fyc(getParam("yield_strength_in_compression")), _Dc(getParam("compressive_damage_at_max_compressive_strength")), _fc(getParam("maximum_strength_in_compression")), _FEc(getParam("fracture_energy_in_compression")), - _at(1.5 * std::sqrt(1 - _Chi) - 0.5), + _ft0(_ft), + _fc0(_fyc), + _at(std::sqrt(9 / 4 + (2 * _Area_S_ep * _dstress_dep) / (_ft0 * _ft0)) - 0.5), _ac((2. * (_fc / _fyc) - 1. + 2. * std::sqrt(Utility::pow<2>(_fc / _fyc) - _fc / _fyc))), _zt((1. + _at) / _at), _zc((1. + _ac) / _ac), @@ -103,9 +110,7 @@ DamagePlasticityStressUpdate::DamagePlasticityStressUpdate(const InputParameters _sqrtPhic_max((1. + _ac) / 2.), _dt_bt(log(1. - _Dt) / log((1. + _at - std::sqrt(1. + _at * _at)) / (2. * _at))), _dc_bc(log(1. - _Dc) / log((1. + _ac) / (2. * _ac))), - _ft0(0.5 * _ft / - ((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.)), @@ -318,10 +323,8 @@ DamagePlasticityStressUpdate::flowPotential(const std::vector & stress_par RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0, 0, 0) .dsecondInvariant(); - Real D = damageVar(stress_params, intnl); - for (unsigned int i = 0; i < _num_sp; ++i) - r[i] = (_alfa_p + d_sqrt_2J2(i, i)) * (1. - D); + r[i] = (_alfa_p + d_sqrt_2J2(i, i)); } void @@ -358,12 +361,10 @@ DamagePlasticityStressUpdate::dflowPotential_dstress( for (unsigned l = 0; l < 3; ++l) dfp(i, j, k, l) += pre * dII(i, j) * dII(k, l); - Real D = damageVar(stress_params, intnl); - 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) * Utility::pow<2>(1. - D); + dr_dstress[i][i] = J2 < _f_tol ? 0. : dfp(i, i, j, j); if (i != j) dr_dstress[j][i] = dr_dstress[i][j]; } @@ -532,17 +533,31 @@ DamagePlasticityStressUpdate::setIntnlDerivativesV(const std::vector & tri } Real -DamagePlasticityStressUpdate::ft(const std::vector & intnl) const +DamagePlasticityStressUpdate::ftbar(const std::vector & intnl) const { Real sqrtPhi_t = std::sqrt(1. + _at * (2. + _at) * intnl[0]); if (_zt > sqrtPhi_t / _at) return _ft0 * std::pow(_zt - sqrtPhi_t / _at, (1. - _dt_bt)) * sqrtPhi_t; else - return _ft0 * 1.E-6; // A very small number (instead of zero) is used for end of softening + return 1.E-6; // A very small number (instead of zero) is used for end of softening } Real -DamagePlasticityStressUpdate::dft(const std::vector & intnl) const +DamagePlasticityStressUpdate::fbar(const std::vector & f0, + const std::vector & a, + const std::vector & exponent, + const std::vector & kappa) const +{ + Real phi = 1. + a * (2. + a) * kappa; + Real sqrt_phi = std::sqrt(Phi); + Real v = sqrt_phi; + Real u = (1 + a) / a - sqrt_phi / a; + Real fbar = f0 * std::pow(u, exponent) * v; + return (u > 0) ? fbar : 1.E-6; +} + +Real +DamagePlasticityStressUpdate::dftbar(const std::vector & intnl) const { Real sqrtPhi_t = std::sqrt(1. + _at * (2. + _at) * intnl[0]); if (_zt > sqrtPhi_t / _at) @@ -553,7 +568,7 @@ DamagePlasticityStressUpdate::dft(const std::vector & intnl) const } Real -DamagePlasticityStressUpdate::fc(const std::vector & intnl) const +DamagePlasticityStressUpdate::fcbar(const std::vector & intnl) const { Real sqrtPhi_c; if (intnl[1] < 1.0) @@ -564,7 +579,23 @@ DamagePlasticityStressUpdate::fc(const std::vector & intnl) const } Real -DamagePlasticityStressUpdate::dfc(const std::vector & intnl) const +DamagePlasticityStressUpdate::dfbar_dkappa(const std::vector & f0, + const std::vector & a, + const std::vector & exponent, + const std::vector & kappa) const +{ + Real phi = 1. + a * (2. + a) * kappa; + Real sqrt_phi = std::sqrt(Phi); + Real v = sqrt_phi; + Real u = (1 + a) / a - sqrt_phi / a; + Real dv_dphi = 1. / (2 * v); + Real du_dphi = -(1 / (2 * a)) * exponent * std::pow(u, exponent - 1) * (1 / v); + Real dfbar_dkappa = f0 * (u * dv_dphi + v * du_dphi) * dphi_dkappa; + return (u > 0) ? dfbar_dkappa : 0.; +} + +Real +DamagePlasticityStressUpdate::dfcbar(const std::vector & intnl) const { if (intnl[1] < 1.0) { @@ -576,22 +607,74 @@ DamagePlasticityStressUpdate::dfc(const std::vector & intnl) const return 0.; } +Real +DamagePlasticityStressUpdate::ft(const std::vector & intnl) const +{ + Real Phi_t = 1. + _at * (2. + _at) * intnl[0]; + Real sqrtPhi_t = std::sqrt(Phi_t); + if (_zt * sqrtPhi_t > Phi_t / _at) + return _ft0 * (_zt * sqrtPhi_t - Phi_t / _at); + else + return 1.E-6; // A very small number (instead of zero) is used for end of softening +} + +Real +DamagePlasticityStressUpdate::fc(const std::vector & intnl) const +{ + Real Phi_c, sqrtPhi_c; + + if (intnl[1] < 1.0) + Phi_c = 1. + _ac * (2. + _ac) * intnl[1]; + sqrtPhi_c = std::sqrt(Phi_c); + else Phi_c = 1. + _ac * (2. + _ac) * 0.999999; + sqrtPhi_c = std::sqrt(Phi_c); + return _fc0 * (_zc * sqrtPhi_c - Phi_c / _ac); +} + +Real +DamagePlasticityStressUpdate::f(const std::vector & f0, + const std::vector & a, + const std::vector & kappa) const +{ + Real phi = 1. + a * (2. + a) * kappa; + Real sqrt_phi = std::sqrt(Phi); + Real v = phi; + Real u = (1 + a) * sqrt_phi; + Real f = (f0 / a) * (u - v); + return (u > v) ? f : 0.; +} + +Real +DamagePlasticityStressUpdate::df_dkappa(const std::vector & f0, + const std::vector & a, + const std::vector & kappa) const +{ + Real phi = 1. + a * (2. + a) * kappa; + Real sqrt_phi = std::sqrt(Phi); + Real v = phi; + Real u = (1 + a) * sqrt_phi; + Real dv_dphi = 1.; + Real du_dphi = (1 + a) / (2 * sqrt_phi); + Real df_dkappa = (f0 / a) * (du_dphi - dv_dphi) * dphi_dkappa; + return (u > v) ? df_dkappa : 0.; +} + Real DamagePlasticityStressUpdate::beta(const std::vector & intnl) const { - return (1. - _alfa) * fc(intnl) / ft(intnl) - (1. + _alfa); + return (1. - _alfa) * fcbar(intnl) / ftbar(intnl) - (1. + _alfa); } Real DamagePlasticityStressUpdate::dbeta0(const std::vector & intnl) const { - return -(1. - _alfa) * fc(intnl) * dft(intnl) / Utility::pow<2>(ft(intnl)); + return -(1. - _alfa) * fcbar(intnl) * dftbar(intnl) / Utility::pow<2>(ftbar(intnl)); } Real DamagePlasticityStressUpdate::dbeta1(const std::vector & intnl) const { - return dfc(intnl) / ft(intnl) * (1. - _alfa); + return dfcbar(intnl) / ftbar(intnl) * (1. - _alfa); } void From 1013d4aca491ebd2a8dfc089ccd6a0821650ade1 Mon Sep 17 00:00:00 2001 From: Prithivirajan Veerappan Date: Tue, 14 Nov 2023 11:07:02 -0700 Subject: [PATCH 2/7] changes as per f and fbar, and derivatives (beta) --- .../materials/DamagePlasticityStressUpdate.h | 123 +++++++------ src/materials/DamagePlasticityStressUpdate.C | 172 ++++++++++-------- 2 files changed, 162 insertions(+), 133 deletions(-) diff --git a/include/materials/DamagePlasticityStressUpdate.h b/include/materials/DamagePlasticityStressUpdate.h index f454a6e1..598f16aa 100644 --- a/include/materials/DamagePlasticityStressUpdate.h +++ b/include/materials/DamagePlasticityStressUpdate.h @@ -116,68 +116,74 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate ///damaged maximum principal stress MaterialProperty & _sigma2; /** - * Obtain the undamaged tensile strength + * Obtain the undamaged strength * @param intnl (Array containing damage states in tension and compression, respectively) * @return value of ft (tensile strength) */ - Real ftbar(const std::vector & intnl) const; + Real fbar(const std::vector & f0, + const std::vector & a, + const std::vector & exponent, + const std::vector & kappa) const; + // Real ftbar(const std::vector & intnl) const; /** * Obtain the partial derivative of the undamaged tensile strength to the damage state * @param intnl (Array containing damage states in tension and compression, respectively) * @return value of dft (partial derivative of the tensile strength to the damage state) */ Real dfbar_dkappa(const std::vector & f0, - const std::vector & a, - const std::vector & exponent, - const std::vector & kappa) const; - - /** - * Obtain the partial derivative of the undamaged tensile strength to the damage state - * @param intnl (Array containing damage states in tension and compression, respectively) - * @return value of dft (partial derivative of the tensile strength to the damage state) - */ - Real dftbar(const std::vector & intnl) const; - - /** - * Obtain the undamaged conpressive strength - * @param intnl (Array containing damage states in tension and compression, respectively) - * @return value of fc (compressive strength) - */ - Real fcbar(const std::vector & intnl) const; - - /** - * Obtain the partial derivative of the undamaged compressive strength to the damage state - * @param intnl (Array containing damage states in tension and compression, respectively) - * @return value of dfc - */ - Real dfcbar(const std::vector & intnl) const; - /** - * Obtain the damaged tensile strength - * @param intnl (Array containing damage states in tension and compression, respectively) - * @return value of ft (tensile strength) - */ - Real ft(const std::vector & intnl) const; - - /** - * Obtain the partial derivative of the damaged tensile strength to the damage state - * @param intnl (Array containing damage states in tension and compression, respectively) - * @return value of dft (partial derivative of the tensile strength to the damage state) - */ - Real dft(const std::vector & intnl) const; - - /** - * Obtain the damaged compressive strength - * @param intnl (Array containing damage states in tension and compression, respectively) - * @return value of fc (compressive strength) - */ - Real fc(const std::vector & intnl) const; - - /** - * Obtain the partial derivative of the damaged compressive strength to the damage state - * @param intnl (Array containing damage states in tension and compression, respectively) - * @return value of dfc - */ - Real dfc(const std::vector & intnl) const; + const std::vector & a, + const std::vector & exponent, + const std::vector & kappa) const; + + // /** + // * Obtain the partial derivative of the undamaged tensile strength to the damage state + // * @param intnl (Array containing damage states in tension and compression, respectively) + // * @return value of dft (partial derivative of the tensile strength to the damage state) + // */ + // Real dftbar(const std::vector & intnl) const; + + // /** + // * Obtain the undamaged conpressive strength + // * @param intnl (Array containing damage states in tension and compression, respectively) + // * @return value of fc (compressive strength) + // */ + // Real fcbar(const std::vector & intnl) const; + + // /** + // * Obtain the partial derivative of the undamaged compressive strength to the damage state + // * @param intnl (Array containing damage states in tension and compression, respectively) + // * @return value of dfc + // */ + // Real dfcbar(const std::vector & intnl) const; + // /** + // * Obtain the damaged tensile strength + // * @param intnl (Array containing damage states in tension and compression, respectively) + // * @return value of ft (tensile strength) + // */ + Real f(const std::vector & f0, + const std::vector & a, + const std::vector & kappa) const; + + // /** + // * Obtain the partial derivative of the damaged tensile strength to the damage state + // * @param intnl (Array containing damage states in tension and compression, respectively) + // * @return value of dft (partial derivative of the tensile strength to the damage state) + // */ + // Real dft(const std::vector & intnl) const; + + // /** + // * Obtain the damaged compressive strength + // * @param intnl (Array containing damage states in tension and compression, respectively) + // * @return value of fc (compressive strength) + // */ + // Real fc(const std::vector & intnl) const; + + // /** + // * Obtain the partial derivative of the damaged compressive strength to the damage state + // * @param intnl (Array containing damage states in tension and compression, respectively) + // * @return value of dfc + // */ + // Real dfc(const std::vector & intnl) const; /** * Obtain the partial derivative of the undamaged strength to the damage state @@ -185,8 +191,8 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate * @return value of dft (partial derivative of the tensile strength to the damage state) */ Real df_dkappa(const std::vector & f0, - const std::vector & exponent, - const std::vector & kappa) const; + const std::vector & a, + const std::vector & kappa) 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 @@ -196,6 +202,13 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate */ Real beta(const std::vector & intnl) const; + /** + * dbeta_dkappa is a derivative of beta wrt. kappa (plastic damage variable) + * @param intnl (Array containing damage states in tension and compression, respectively) + * @return value of dbeta_dkappa + */ + Real dbeta_dkappa(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) diff --git a/src/materials/DamagePlasticityStressUpdate.C b/src/materials/DamagePlasticityStressUpdate.C index 28756a90..ea12715e 100644 --- a/src/materials/DamagePlasticityStressUpdate.C +++ b/src/materials/DamagePlasticityStressUpdate.C @@ -265,7 +265,7 @@ DamagePlasticityStressUpdate::yieldFunctionValuesV(const std::vector & str yf[0] = 1. / (1. - _alfa) * (_alfa * I1 + _sqrt3 * std::sqrt(J2) + beta(intnl) * (stress_params[2] < 0. ? 0. : stress_params[2])) - - fc(intnl); + fbar(_fc0, _ac, 1. - _dc_bc, _intnl1); } void @@ -287,7 +287,7 @@ DamagePlasticityStressUpdate::computeAllQV(const std::vector & stress_para all_q[0].f = 1. / (1. - _alfa) * (_alfa * I1 + _sqrt3 * std::sqrt(J2) + beta(intnl) * (stress_params[2] < 0. ? 0. : stress_params[2])) - - fc(intnl); + fbar(_fc0, _ac, 1. - _dc_bc, _intnl1); for (unsigned i = 0; i < _num_sp; ++i) if (J2 == 0) @@ -390,8 +390,8 @@ DamagePlasticityStressUpdate::hardPotential(const std::vector & stress_par weighfac(stress_params, wf); std::vector r(3); flowPotential(stress_params, intnl, r); - h[0] = wf * ft(intnl) / _gt[_qp] * r[2]; - h[1] = -(1. - wf) * fc(intnl) / _gc[_qp] * r[0]; + h[0] = wf * f(_ft0, _at, _intnl0) / _gt[_qp] * r[2]; + h[1] = -(1. - wf) * f(_fc0, _ac, _intnl1) / _gc[_qp] * r[0]; } void @@ -421,15 +421,19 @@ DamagePlasticityStressUpdate::dhardPotential_dintnl( const std::vector & intnl, std::vector> & dh_dintnl) const { - Real wf; + Real wf, dft, dfc; weighfac(stress_params, wf); + + dft = df_dkappa(_ft0, _at, _intnl0); + dfc = df_dkappa(_fc0, _ac, _intnl1); + std::vector r(3); flowPotential(stress_params, intnl, r); - dh_dintnl[0][0] = wf * dft(intnl) / _gt[_qp] * r[2]; + dh_dintnl[0][0] = wf * dft / _gt[_qp] * r[2]; dh_dintnl[0][1] = 0.; dh_dintnl[1][0] = 0.; - dh_dintnl[1][1] = -(1 - wf) * dfc(intnl) / _gc[_qp] * r[0]; + dh_dintnl[1][1] = -(1 - wf) * dfc / _gc[_qp] * r[0]; } void @@ -532,15 +536,15 @@ DamagePlasticityStressUpdate::setIntnlDerivativesV(const std::vector & tri dintnl[i][j] = dlam_dsig[j] * h[i] + lam * dh_dsig[i][j]; } -Real -DamagePlasticityStressUpdate::ftbar(const std::vector & intnl) const -{ - Real sqrtPhi_t = std::sqrt(1. + _at * (2. + _at) * intnl[0]); - if (_zt > sqrtPhi_t / _at) - return _ft0 * std::pow(_zt - sqrtPhi_t / _at, (1. - _dt_bt)) * sqrtPhi_t; - else - return 1.E-6; // A very small number (instead of zero) is used for end of softening -} +// Real +// DamagePlasticityStressUpdate::ftbar(const std::vector & intnl) const +// { +// Real sqrtPhi_t = std::sqrt(1. + _at * (2. + _at) * intnl[0]); +// if (_zt > sqrtPhi_t / _at) +// return _ft0 * std::pow(_zt - sqrtPhi_t / _at, (1. - _dt_bt)) * sqrtPhi_t; +// else +// return 1.E-6; // A very small number (instead of zero) is used for end of softening +// } Real DamagePlasticityStressUpdate::fbar(const std::vector & f0, @@ -556,28 +560,6 @@ DamagePlasticityStressUpdate::fbar(const std::vector & f0, return (u > 0) ? fbar : 1.E-6; } -Real -DamagePlasticityStressUpdate::dftbar(const std::vector & intnl) const -{ - Real sqrtPhi_t = std::sqrt(1. + _at * (2. + _at) * intnl[0]); - if (_zt > sqrtPhi_t / _at) - return _ft0 * _dPhit / (2 * sqrtPhi_t) * std::pow(_zt - sqrtPhi_t / _at, -_dt_bt) * - (_zt - (2. - _dt_bt) * sqrtPhi_t / _at); - else - return 0.; -} - -Real -DamagePlasticityStressUpdate::fcbar(const std::vector & intnl) const -{ - Real sqrtPhi_c; - if (intnl[1] < 1.0) - sqrtPhi_c = std::sqrt(1. + _ac * (2. + _ac) * intnl[1]); - else - sqrtPhi_c = std::sqrt(1. + _ac * (2. + _ac) * 0.99); - return _fc0 * std::pow((_zc - sqrtPhi_c / _ac), (1. - _dc_bc)) * sqrtPhi_c; -} - Real DamagePlasticityStressUpdate::dfbar_dkappa(const std::vector & f0, const std::vector & a, @@ -594,43 +576,6 @@ DamagePlasticityStressUpdate::dfbar_dkappa(const std::vector & f0, return (u > 0) ? dfbar_dkappa : 0.; } -Real -DamagePlasticityStressUpdate::dfcbar(const std::vector & intnl) const -{ - if (intnl[1] < 1.0) - { - Real sqrtPhi_c = std::sqrt(1. + _ac * (2. + _ac) * intnl[1]); - return _fc0 * _dPhic / (2. * sqrtPhi_c) * std::pow(_zc - sqrtPhi_c / _ac, -_dc_bc) * - (_zc - (2. - _dc_bc) * sqrtPhi_c / _ac); - } - else - return 0.; -} - -Real -DamagePlasticityStressUpdate::ft(const std::vector & intnl) const -{ - Real Phi_t = 1. + _at * (2. + _at) * intnl[0]; - Real sqrtPhi_t = std::sqrt(Phi_t); - if (_zt * sqrtPhi_t > Phi_t / _at) - return _ft0 * (_zt * sqrtPhi_t - Phi_t / _at); - else - return 1.E-6; // A very small number (instead of zero) is used for end of softening -} - -Real -DamagePlasticityStressUpdate::fc(const std::vector & intnl) const -{ - Real Phi_c, sqrtPhi_c; - - if (intnl[1] < 1.0) - Phi_c = 1. + _ac * (2. + _ac) * intnl[1]; - sqrtPhi_c = std::sqrt(Phi_c); - else Phi_c = 1. + _ac * (2. + _ac) * 0.999999; - sqrtPhi_c = std::sqrt(Phi_c); - return _fc0 * (_zc * sqrtPhi_c - Phi_c / _ac); -} - Real DamagePlasticityStressUpdate::f(const std::vector & f0, const std::vector & a, @@ -659,22 +604,93 @@ DamagePlasticityStressUpdate::df_dkappa(const std::vector & f0, return (u > v) ? df_dkappa : 0.; } +// Real +// DamagePlasticityStressUpdate::dftbar(const std::vector & intnl) const +// { +// Real sqrtPhi_t = std::sqrt(1. + _at * (2. + _at) * intnl[0]); +// if (_zt > sqrtPhi_t / _at) +// return _ft0 * _dPhit / (2 * sqrtPhi_t) * std::pow(_zt - sqrtPhi_t / _at, -_dt_bt) * +// (_zt - (2. - _dt_bt) * sqrtPhi_t / _at); +// else +// return 0.; +// } + +// Real +// DamagePlasticityStressUpdate::fcbar(const std::vector & intnl) const +// { +// Real sqrtPhi_c; +// if (intnl[1] < 1.0) +// sqrtPhi_c = std::sqrt(1. + _ac * (2. + _ac) * intnl[1]); +// else +// sqrtPhi_c = std::sqrt(1. + _ac * (2. + _ac) * 0.99); +// return _fc0 * std::pow((_zc - sqrtPhi_c / _ac), (1. - _dc_bc)) * sqrtPhi_c; +// } + +// Real +// DamagePlasticityStressUpdate::dfcbar(const std::vector & intnl) const +// { +// if (intnl[1] < 1.0) +// { +// Real sqrtPhi_c = std::sqrt(1. + _ac * (2. + _ac) * intnl[1]); +// return _fc0 * _dPhic / (2. * sqrtPhi_c) * std::pow(_zc - sqrtPhi_c / _ac, -_dc_bc) * +// (_zc - (2. - _dc_bc) * sqrtPhi_c / _ac); +// } +// else +// return 0.; +// } + +// Real +// DamagePlasticityStressUpdate::ft(const std::vector & intnl) const +// { +// Real Phi_t = 1. + _at * (2. + _at) * intnl[0]; +// Real sqrtPhi_t = std::sqrt(Phi_t); +// if (_zt * sqrtPhi_t > Phi_t / _at) +// return _ft0 * (_zt * sqrtPhi_t - Phi_t / _at); +// else +// return 1.E-6; // A very small number (instead of zero) is used for end of softening +// } + +// Real +// DamagePlasticityStressUpdate::fc(const std::vector & intnl) const +// { +// Real Phi_c, sqrtPhi_c; + +// if (intnl[1] < 1.0) +// Phi_c = 1. + _ac * (2. + _ac) * intnl[1]; +// sqrtPhi_c = std::sqrt(Phi_c); +// else Phi_c = 1. + _ac * (2. + _ac) * 0.999999; +// sqrtPhi_c = std::sqrt(Phi_c); +// return _fc0 * (_zc * sqrtPhi_c - Phi_c / _ac); +// } + Real DamagePlasticityStressUpdate::beta(const std::vector & intnl) const { - return (1. - _alfa) * fcbar(intnl) / ftbar(intnl) - (1. + _alfa); + Real fcbar, ftbar; + fcbar = fbar(_fc0, _ac, 1. - _dc_bc, _intnl1); + ftbar = fbar(_ft0, _at, 1. - _dt_bt, _intnl0); + + return (1. - _alfa) * (fcbar / ftbar) - (1. + _alfa); } Real DamagePlasticityStressUpdate::dbeta0(const std::vector & intnl) const { - return -(1. - _alfa) * fcbar(intnl) * dftbar(intnl) / Utility::pow<2>(ftbar(intnl)); + Real fcbar, ftbar, dftbar ; + fcbar = fbar(_fc0, _ac, 1. - _dc_bc, _intnl1); + ftbar = fbar(_ft0, _at, 1. - _dt_bt, _intnl0); + dftbar = dfbar_dkappa(_ft0, _at, 1. - _dt_bt, _intnl0); + return -(1. - _alfa) * fcbar * dftbar / Utility::pow<2>(ftbar); } Real DamagePlasticityStressUpdate::dbeta1(const std::vector & intnl) const { - return dfcbar(intnl) / ftbar(intnl) * (1. - _alfa); + Real fcbar, ftbar, dfcbar; + fcbar = fbar(_fc0, _ac, 1. - _dc_bc, _intnl1); + ftbar = fbar(_ft0, _at, 1. - _dt_bt, _intnl0); + dfcbar = dfbar_dkappa(_fc0, _ac, 1. - _dc_bc, _intnl1); + return dfcbar / ftbar * (1. - _alfa); } void From e508432f991a1b070c80c66f4e39754ee1be6cb9 Mon Sep 17 00:00:00 2001 From: Prithivirajan Veerappan Date: Mon, 8 Jan 2024 15:42:16 -0700 Subject: [PATCH 3/7] Refined declarations and definitions --- .../materials/DamagePlasticityStressUpdate.h | 30 +++--- src/materials/DamagePlasticityStressUpdate.C | 93 ++++++++++--------- 2 files changed, 65 insertions(+), 58 deletions(-) diff --git a/include/materials/DamagePlasticityStressUpdate.h b/include/materials/DamagePlasticityStressUpdate.h index 598f16aa..27a1e6d2 100644 --- a/include/materials/DamagePlasticityStressUpdate.h +++ b/include/materials/DamagePlasticityStressUpdate.h @@ -46,7 +46,7 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate ///fracture_energy in tension (user parameter) const Real _FEt; ///Area under the stress and plastic strain in uniaxial tension (user parameter) - const Real _Area_S_ep; + const Real _area_S_ep; ///yield strength in compression (user parameter) const Real _fyc; ///compressive damage at maximum compressive strength (user parameter) @@ -120,20 +120,20 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate * @param intnl (Array containing damage states in tension and compression, respectively) * @return value of ft (tensile strength) */ - Real fbar(const std::vector & f0, - const std::vector & a, - const std::vector & exponent, - const std::vector & kappa) const; + Real fbar(const Real & f0, + const Real & a, + const Real & exponent, + const Real & kappa) const; // Real ftbar(const std::vector & intnl) const; /** * Obtain the partial derivative of the undamaged tensile strength to the damage state * @param intnl (Array containing damage states in tension and compression, respectively) * @return value of dft (partial derivative of the tensile strength to the damage state) */ - Real dfbar_dkappa(const std::vector & f0, - const std::vector & a, - const std::vector & exponent, - const std::vector & kappa) const; + Real dfbar_dkappa(const Real & f0, + const Real & a, + const Real & exponent, + const Real & kappa) const; // /** // * Obtain the partial derivative of the undamaged tensile strength to the damage state @@ -160,9 +160,9 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate // * @param intnl (Array containing damage states in tension and compression, respectively) // * @return value of ft (tensile strength) // */ - Real f(const std::vector & f0, - const std::vector & a, - const std::vector & kappa) const; + Real f(const Real & f0, + const Real & a, + const Real & kappa) const; // /** // * Obtain the partial derivative of the damaged tensile strength to the damage state @@ -190,9 +190,9 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate * @param intnl (Array containing damage states in tension and compression, respectively) * @return value of dft (partial derivative of the tensile strength to the damage state) */ - Real df_dkappa(const std::vector & f0, - const std::vector & a, - const std::vector & kappa) const; + Real df_dkappa(const Real & f0, + const Real & a, + const Real & kappa) 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 diff --git a/src/materials/DamagePlasticityStressUpdate.C b/src/materials/DamagePlasticityStressUpdate.C index ea12715e..4a233d8c 100644 --- a/src/materials/DamagePlasticityStressUpdate.C +++ b/src/materials/DamagePlasticityStressUpdate.C @@ -37,7 +37,7 @@ DamagePlasticityStressUpdate::validParams() "stiffness recovery factor"); params.addRangeCheckedParam( "ft_ep_slope_factor_at_zero_ep", - "ft_ep_slope_factor_at_zero_ep <= 1 & ft_ep_slope_factor_at_zero_ep >= 0", + "ft_ep_slope_factor_at_zero_ep <0", "slope of ft vs plastic strain curve at zero plastic strain"); params.addRequiredParam( "tensile_damage_at_half_tensile_strength", @@ -49,8 +49,8 @@ DamagePlasticityStressUpdate::validParams() "fracture_energy_in_tension >= 0", "Fracture energy of concrete in uniaxial tension"); params.addRangeCheckedParam( - "Area_under_stress_ep", - "Area_under_stress_ep >= 0", + "area_under_stress_ep", + "area_under_stress_ep >= 0", "Area under the stress and plastic strain in uniaxial tension"); params.addRangeCheckedParam("yield_strength_in_compression", @@ -91,7 +91,7 @@ DamagePlasticityStressUpdate::DamagePlasticityStressUpdate(const InputParameters _Dt(getParam("tensile_damage_at_half_tensile_strength")), _ft(getParam("yield_strength_in_tension")), _FEt(getParam("fracture_energy_in_tension")), - _Area_S_ep(getParam("Area_under_stress_ep")), + _area_S_ep(getParam("area_under_stress_ep")), _fyc(getParam("yield_strength_in_compression")), _Dc(getParam("compressive_damage_at_max_compressive_strength")), @@ -100,7 +100,7 @@ DamagePlasticityStressUpdate::DamagePlasticityStressUpdate(const InputParameters _ft0(_ft), _fc0(_fyc), - _at(std::sqrt(9 / 4 + (2 * _Area_S_ep * _dstress_dep) / (_ft0 * _ft0)) - 0.5), + _at(std::sqrt(9. / 4. + (2 * _area_S_ep * _dstress_dep) / (_ft0 * _ft0)) - 0.5), _ac((2. * (_fc / _fyc) - 1. + 2. * std::sqrt(Utility::pow<2>(_fc / _fyc) - _fc / _fyc))), _zt((1. + _at) / _at), _zc((1. + _ac) / _ac), @@ -265,7 +265,7 @@ DamagePlasticityStressUpdate::yieldFunctionValuesV(const std::vector & str yf[0] = 1. / (1. - _alfa) * (_alfa * I1 + _sqrt3 * std::sqrt(J2) + beta(intnl) * (stress_params[2] < 0. ? 0. : stress_params[2])) - - fbar(_fc0, _ac, 1. - _dc_bc, _intnl1); + fbar(_fc0, _ac, 1. - _dc_bc, intnl[1]); } void @@ -277,6 +277,11 @@ DamagePlasticityStressUpdate::computeAllQV(const std::vector & stress_para Real J2 = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0, 0, 0) .secondInvariant(); RankTwoTensor d_sqrt_3J2; + + Real fcbar, dfcbar; + fcbar = fbar(_fc0, _ac, 1. - _dc_bc, intnl[1]); + dfcbar = dfbar_dkappa(_fc0, _ac, 1. - _dc_bc, intnl[1]); + if (J2 == 0) d_sqrt_3J2 = RankTwoTensor(); else @@ -287,7 +292,7 @@ DamagePlasticityStressUpdate::computeAllQV(const std::vector & stress_para all_q[0].f = 1. / (1. - _alfa) * (_alfa * I1 + _sqrt3 * std::sqrt(J2) + beta(intnl) * (stress_params[2] < 0. ? 0. : stress_params[2])) - - fbar(_fc0, _ac, 1. - _dc_bc, _intnl1); + fcbar; for (unsigned i = 0; i < _num_sp; ++i) if (J2 == 0) @@ -301,7 +306,7 @@ DamagePlasticityStressUpdate::computeAllQV(const std::vector & stress_para 1. / (1. - _alfa) * (dbeta0(intnl) * (stress_params[2] < 0. ? 0. : stress_params[2])); all_q[0].df_di[1] = 1. / (1. - _alfa) * (dbeta1(intnl) * (stress_params[2] < 0. ? 0. : stress_params[2])) - - dfc(intnl); + dfcbar; flowPotential(stress_params, intnl, all_q[0].dg); dflowPotential_dstress(stress_params, intnl, all_q[0].d2g); @@ -386,12 +391,14 @@ DamagePlasticityStressUpdate::hardPotential(const std::vector & stress_par const std::vector & intnl, std::vector & h) const { - Real wf; + Real wf, ft, fc; weighfac(stress_params, wf); std::vector r(3); + ft = f(_ft0, _at, intnl[0]); + fc = f(_fc0, _ac, intnl[1]); flowPotential(stress_params, intnl, r); - h[0] = wf * f(_ft0, _at, _intnl0) / _gt[_qp] * r[2]; - h[1] = -(1. - wf) * f(_fc0, _ac, _intnl1) / _gc[_qp] * r[0]; + h[0] = wf * ft / _gt[_qp] * r[2]; + h[1] = -(1. - wf) * fc / _gc[_qp] * r[0]; } void @@ -399,9 +406,11 @@ DamagePlasticityStressUpdate::dhardPotential_dstress(const std::vector & s const std::vector & intnl, std::vector> & dh_dsig) const { - Real wf; + Real wf, ft, fc; std::vector dwf(3); dweighfac(stress_params, wf, dwf); + ft = f(_ft0, _at, intnl[0]); + fc = f(_fc0, _ac, intnl[1]); std::vector r(3); flowPotential(stress_params, intnl, r); @@ -410,8 +419,8 @@ DamagePlasticityStressUpdate::dhardPotential_dstress(const std::vector & s for (unsigned i = 0; i < _num_sp; ++i) { - dh_dsig[0][i] = (wf * dr_dsig[2][i] + dwf[i] * r[2]) * ft(intnl) / _gt[_qp]; - dh_dsig[1][i] = -((1. - wf) * dr_dsig[0][i] - dwf[i] * r[0]) * fc(intnl) / _gc[_qp]; + dh_dsig[0][i] = (wf * dr_dsig[2][i] + dwf[i] * r[2]) * ft / _gt[_qp]; + dh_dsig[1][i] = -((1. - wf) * dr_dsig[0][i] - dwf[i] * r[0]) * fc / _gc[_qp]; } } @@ -424,8 +433,8 @@ DamagePlasticityStressUpdate::dhardPotential_dintnl( Real wf, dft, dfc; weighfac(stress_params, wf); - dft = df_dkappa(_ft0, _at, _intnl0); - dfc = df_dkappa(_fc0, _ac, _intnl1); + dft = df_dkappa(_ft0, _at, intnl[0]); + dfc = df_dkappa(_fc0, _ac, intnl[1]); std::vector r(3); flowPotential(stress_params, intnl, r); @@ -547,13 +556,13 @@ DamagePlasticityStressUpdate::setIntnlDerivativesV(const std::vector & tri // } Real -DamagePlasticityStressUpdate::fbar(const std::vector & f0, - const std::vector & a, - const std::vector & exponent, - const std::vector & kappa) const +DamagePlasticityStressUpdate::fbar(const Real & f0, + const Real & a, + const Real & exponent, + const Real & kappa) const { Real phi = 1. + a * (2. + a) * kappa; - Real sqrt_phi = std::sqrt(Phi); + Real sqrt_phi = std::sqrt(phi); Real v = sqrt_phi; Real u = (1 + a) / a - sqrt_phi / a; Real fbar = f0 * std::pow(u, exponent) * v; @@ -561,13 +570,14 @@ DamagePlasticityStressUpdate::fbar(const std::vector & f0, } Real -DamagePlasticityStressUpdate::dfbar_dkappa(const std::vector & f0, - const std::vector & a, - const std::vector & exponent, - const std::vector & kappa) const +DamagePlasticityStressUpdate::dfbar_dkappa(const Real & f0, + const Real & a, + const Real & exponent, + const Real & kappa) const { Real phi = 1. + a * (2. + a) * kappa; - Real sqrt_phi = std::sqrt(Phi); + Real dphi_dkappa = a * (2. + a); + Real sqrt_phi = std::sqrt(phi); Real v = sqrt_phi; Real u = (1 + a) / a - sqrt_phi / a; Real dv_dphi = 1. / (2 * v); @@ -577,12 +587,10 @@ DamagePlasticityStressUpdate::dfbar_dkappa(const std::vector & f0, } Real -DamagePlasticityStressUpdate::f(const std::vector & f0, - const std::vector & a, - const std::vector & kappa) const +DamagePlasticityStressUpdate::f(const Real & f0, const Real & a, const Real & kappa) const { Real phi = 1. + a * (2. + a) * kappa; - Real sqrt_phi = std::sqrt(Phi); + Real sqrt_phi = std::sqrt(phi); Real v = phi; Real u = (1 + a) * sqrt_phi; Real f = (f0 / a) * (u - v); @@ -590,12 +598,11 @@ DamagePlasticityStressUpdate::f(const std::vector & f0, } Real -DamagePlasticityStressUpdate::df_dkappa(const std::vector & f0, - const std::vector & a, - const std::vector & kappa) const +DamagePlasticityStressUpdate::df_dkappa(const Real & f0, const Real & a, const Real & kappa) const { Real phi = 1. + a * (2. + a) * kappa; - Real sqrt_phi = std::sqrt(Phi); + Real dphi_dkappa = a * (2. + a); + Real sqrt_phi = std::sqrt(phi); Real v = phi; Real u = (1 + a) * sqrt_phi; Real dv_dphi = 1.; @@ -667,8 +674,8 @@ Real DamagePlasticityStressUpdate::beta(const std::vector & intnl) const { Real fcbar, ftbar; - fcbar = fbar(_fc0, _ac, 1. - _dc_bc, _intnl1); - ftbar = fbar(_ft0, _at, 1. - _dt_bt, _intnl0); + fcbar = fbar(_fc0, _ac, 1. - _dc_bc, intnl[1]); + ftbar = fbar(_ft0, _at, 1. - _dt_bt, intnl[0]); return (1. - _alfa) * (fcbar / ftbar) - (1. + _alfa); } @@ -676,10 +683,10 @@ DamagePlasticityStressUpdate::beta(const std::vector & intnl) const Real DamagePlasticityStressUpdate::dbeta0(const std::vector & intnl) const { - Real fcbar, ftbar, dftbar ; - fcbar = fbar(_fc0, _ac, 1. - _dc_bc, _intnl1); - ftbar = fbar(_ft0, _at, 1. - _dt_bt, _intnl0); - dftbar = dfbar_dkappa(_ft0, _at, 1. - _dt_bt, _intnl0); + Real fcbar, ftbar, dftbar; + fcbar = fbar(_fc0, _ac, 1. - _dc_bc, intnl[1]); + ftbar = fbar(_ft0, _at, 1. - _dt_bt, intnl[0]); + dftbar = dfbar_dkappa(_ft0, _at, 1. - _dt_bt, intnl[0]); return -(1. - _alfa) * fcbar * dftbar / Utility::pow<2>(ftbar); } @@ -687,9 +694,9 @@ Real DamagePlasticityStressUpdate::dbeta1(const std::vector & intnl) const { Real fcbar, ftbar, dfcbar; - fcbar = fbar(_fc0, _ac, 1. - _dc_bc, _intnl1); - ftbar = fbar(_ft0, _at, 1. - _dt_bt, _intnl0); - dfcbar = dfbar_dkappa(_fc0, _ac, 1. - _dc_bc, _intnl1); + fcbar = fbar(_fc0, _ac, 1. - _dc_bc, intnl[1]); + ftbar = fbar(_ft0, _at, 1. - _dt_bt, intnl[0]); + dfcbar = dfbar_dkappa(_fc0, _ac, 1. - _dc_bc, intnl[1]); return dfcbar / ftbar * (1. - _alfa); } From 752d3a39f8212c0b71907f8b06389ad25c90f5cd Mon Sep 17 00:00:00 2001 From: Prithivirajan Veerappan Date: Thu, 16 May 2024 10:03:49 -0600 Subject: [PATCH 4/7] Changes in at calc, fbar derivative, gamma calc and derivative, dpotential_dstess --- .../materials/DamagePlasticityStressUpdate.h | 4 +- src/materials/DamagePlasticityStressUpdate.C | 205 +++++------------- 2 files changed, 58 insertions(+), 151 deletions(-) diff --git a/include/materials/DamagePlasticityStressUpdate.h b/include/materials/DamagePlasticityStressUpdate.h index 27a1e6d2..ebb76b6e 100644 --- a/include/materials/DamagePlasticityStressUpdate.h +++ b/include/materials/DamagePlasticityStressUpdate.h @@ -38,15 +38,13 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate ///stiffness recovery factor (user parameter) const Real _s0; ///ft_ep_slope_factor_at_zero_ep (user parameter) - const Real _dstress_dep; + 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; - ///Area under the stress and plastic strain in uniaxial tension (user parameter) - const Real _area_S_ep; ///yield strength in compression (user parameter) const Real _fyc; ///compressive damage at maximum compressive strength (user parameter) diff --git a/src/materials/DamagePlasticityStressUpdate.C b/src/materials/DamagePlasticityStressUpdate.C index 4a233d8c..6a2f4d85 100644 --- a/src/materials/DamagePlasticityStressUpdate.C +++ b/src/materials/DamagePlasticityStressUpdate.C @@ -37,8 +37,9 @@ DamagePlasticityStressUpdate::validParams() "stiffness recovery factor"); params.addRangeCheckedParam( "ft_ep_slope_factor_at_zero_ep", - "ft_ep_slope_factor_at_zero_ep <0", + "ft_ep_slope_factor_at_zero_ep <= 1 & ft_ep_slope_factor_at_zero_ep >= 0", "slope of ft vs plastic strain curve at zero plastic strain"); + params.addRequiredParam( "tensile_damage_at_half_tensile_strength", "fraction of the elastic recovery slope in tension at 0.5*ft0 after yielding"); @@ -48,11 +49,6 @@ DamagePlasticityStressUpdate::validParams() params.addRangeCheckedParam("fracture_energy_in_tension", "fracture_energy_in_tension >= 0", "Fracture energy of concrete in uniaxial tension"); - params.addRangeCheckedParam( - "area_under_stress_ep", - "area_under_stress_ep >= 0", - "Area under the stress and plastic strain in uniaxial tension"); - params.addRangeCheckedParam("yield_strength_in_compression", "yield_strength_in_compression >= 0", "Absolute yield compressice strength"); @@ -87,12 +83,10 @@ DamagePlasticityStressUpdate::DamagePlasticityStressUpdate(const InputParameters _alfa_p(getParam("dilatancy_factor")), _s0(getParam("stiff_recovery_factor")), - _dstress_dep(getParam("ft_ep_slope_factor_at_zero_ep")), + _Chi(getParam("ft_ep_slope_factor_at_zero_ep")), _Dt(getParam("tensile_damage_at_half_tensile_strength")), _ft(getParam("yield_strength_in_tension")), _FEt(getParam("fracture_energy_in_tension")), - _area_S_ep(getParam("area_under_stress_ep")), - _fyc(getParam("yield_strength_in_compression")), _Dc(getParam("compressive_damage_at_max_compressive_strength")), _fc(getParam("maximum_strength_in_compression")), @@ -100,14 +94,10 @@ DamagePlasticityStressUpdate::DamagePlasticityStressUpdate(const InputParameters _ft0(_ft), _fc0(_fyc), - _at(std::sqrt(9. / 4. + (2 * _area_S_ep * _dstress_dep) / (_ft0 * _ft0)) - 0.5), + _at(1.5 * std::sqrt(1 - _Chi) - 0.5), _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)), - _dPhic(_ac * (2. + _ac)), - _sqrtPhit_max((1. + _at + std::sqrt(1. + _at * _at)) / 2.), - _sqrtPhic_max((1. + _ac) / 2.), _dt_bt(log(1. - _Dt) / log((1. + _at - std::sqrt(1. + _at * _at)) / (2. * _at))), _dc_bc(log(1. - _Dc) / log((1. + _ac) / (2. * _ac))), @@ -315,7 +305,7 @@ DamagePlasticityStressUpdate::computeAllQV(const std::vector & stress_para void DamagePlasticityStressUpdate::flowPotential(const std::vector & stress_params, - const std::vector & intnl, + const std::vector & /* intnl */, std::vector & r) const { Real J2 = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0, 0, 0) @@ -335,7 +325,7 @@ DamagePlasticityStressUpdate::flowPotential(const std::vector & stress_par void DamagePlasticityStressUpdate::dflowPotential_dstress( const std::vector & stress_params, - const std::vector & intnl, + const std::vector & /* intnl */, std::vector> & dr_dstress) const { Real J2 = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0, 0, 0) @@ -344,35 +334,27 @@ DamagePlasticityStressUpdate::dflowPotential_dstress( .dsecondInvariant(); RankTwoTensor d_sqrt_2J2; RankFourTensor dfp; - Real pre; + RankTwoTensor d2J2_dsigi_dsigj = + RankTwoTensor(2. / 3., 2. / 3., 2. / 3., -1. / 3., -1. / 3., -1. / 3.); + std::vector dJ2_dsigi(3); + for (unsigned i = 0; i < 3; ++i) + dJ2_dsigi[i] = + (2 * stress_params[i] - stress_params[(i + 1) % 3] - stress_params[(i + 2) % 3]) / 3; + if (J2 == 0) { - d_sqrt_2J2 = RankTwoTensor(); - dfp = RankFourTensor(); - pre = 0; + for (unsigned i = 0; i < 3; ++i) + for (unsigned j = 0; j < 3; ++j) + dr_dstress[i][j] = 0.0; + } else { - d_sqrt_2J2 = 0.5 * std::sqrt(2.0 / J2) * dII; - dfp = 0.5 * std::sqrt(2.0 / J2) * - RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0, 0, 0) - .d2secondInvariant(); - pre = -0.25 * std::sqrt(2.0) * std::pow(J2, -1.5); + for (unsigned i = 0; i < 3; ++i) + for (unsigned j = 0; j < 3; ++j) + dr_dstress[i][j] = 0.5 * (std::sqrt(2.0 / J2) * d2J2_dsigi_dsigj(i, j) - + (1 / std::sqrt(2)) * std::pow(J2, -1.5) * dJ2_dsigi[i]*dJ2_dsigi[j]); } - - for (unsigned i = 0; i < 3; ++i) - for (unsigned j = 0; j < 3; ++j) - for (unsigned k = 0; k < 3; ++k) - for (unsigned l = 0; l < 3; ++l) - dfp(i, j, k, l) += pre * dII(i, j) * dII(k, l); - - 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); - if (i != j) - dr_dstress[j][i] = dr_dstress[i][j]; - } } void @@ -461,32 +443,31 @@ DamagePlasticityStressUpdate::setIntnlValuesV(const std::vector & trial_st const std::vector & intnl_old, std::vector & intnl) const { + RankTwoTensor sigma_trial = RankTwoTensor( + trial_stress_params[0], trial_stress_params[1], trial_stress_params[2], 0, 0, 0); Real I1_trial = trial_stress_params[0] + trial_stress_params[1] + trial_stress_params[2]; - Real J2_trial = - RankTwoTensor(trial_stress_params[0], trial_stress_params[1], trial_stress_params[2], 0, 0, 0) - .secondInvariant(); - Real invsqrt2J2_trial = 1. / std::sqrt(2. * J2_trial); + RankTwoTensor Identity_tensor = RankTwoTensor(1, 1, 1, 0, 0, 0); + RankTwoTensor sigmadev_trial = sigma_trial - (I1_trial / 3.) * Identity_tensor; + Real norm_sigmadev_trial = sigmadev_trial.L2norm(); + Real G = 0.5 * (_Eij[0][0] - _Eij[0][1]); // Lame's mu Real K = _Eij[0][1] + 2. * G / 3.; // Bulk modulus - Real C1 = (2. * G * invsqrt2J2_trial); - Real C2 = -(I1_trial / 3. * G * invsqrt2J2_trial - 3. * K * _alfa_p); Real C3 = 3. * K * _alfa_p; + RankTwoTensor denominator_tensor = + (2 * G / norm_sigmadev_trial) * sigmadev_trial + C3 * Identity_tensor; + + RankTwoTensor dsig = RankTwoTensor(trial_stress_params[0] - current_stress_params[0], trial_stress_params[1] - current_stress_params[1], trial_stress_params[2] - current_stress_params[2], 0., 0., 0.); - RankTwoTensor fac = J2_trial < _f_tol ? C3 * RankTwoTensor(1., 1., 1., 0., 0., 0.) - : RankTwoTensor(C1 * trial_stress_params[0] - C2, - C1 * trial_stress_params[1] - C2, - C1 * trial_stress_params[2] - C2, - 0., - 0., - 0.); - - Real lam = dsig.L2norm() / fac.L2norm(); + +// Implementing Eqn. (21) of Lee & Fenves, 2001, with backsubstituting eqn. (22) + Real lam = dsig.L2norm() / denominator_tensor.L2norm(); + std::vector h(2); hardPotential(current_stress_params, intnl_old, h); @@ -500,38 +481,33 @@ DamagePlasticityStressUpdate::setIntnlDerivativesV(const std::vector & tri const std::vector & intnl, std::vector> & dintnl) const { + RankTwoTensor sigma_trial = RankTwoTensor( + trial_stress_params[0], trial_stress_params[1], trial_stress_params[2], 0, 0, 0); Real I1_trial = trial_stress_params[0] + trial_stress_params[1] + trial_stress_params[2]; - Real J2_trial = - RankTwoTensor(trial_stress_params[0], trial_stress_params[1], trial_stress_params[2], 0, 0, 0) - .secondInvariant(); - Real invsqrt2J2_trial = 1. / std::sqrt(2. * J2_trial); + RankTwoTensor Identity_tensor = RankTwoTensor(1, 1, 1, 0, 0, 0); + RankTwoTensor sigmadev_trial = sigma_trial - (I1_trial / 3.) * Identity_tensor; + Real norm_sigmadev_trial = sigmadev_trial.L2norm(); + Real G = 0.5 * (_Eij[0][0] - _Eij[0][1]); // Lame's mu Real K = _Eij[0][1] + 2. * G / 3.; // Bulk modulus - Real C1 = (2. * G * invsqrt2J2_trial); - Real C2 = -(I1_trial / 3. * G * invsqrt2J2_trial - 3. * K * _alfa_p); Real C3 = 3. * K * _alfa_p; + RankTwoTensor denominator_tensor = + 2 * G * (sigmadev_trial / norm_sigmadev_trial) + C3 * Identity_tensor; RankTwoTensor dsig = RankTwoTensor(trial_stress_params[0] - current_stress_params[0], trial_stress_params[1] - current_stress_params[1], trial_stress_params[2] - current_stress_params[2], 0., 0., 0.); - RankTwoTensor fac = J2_trial < _f_tol ? C3 * RankTwoTensor(1., 1., 1., 0., 0., 0.) - : RankTwoTensor(C1 * trial_stress_params[0] - C2, - C1 * trial_stress_params[1] - C2, - C1 * trial_stress_params[2] - C2, - 0., - 0., - 0.); - Real lam = dsig.L2norm() / fac.L2norm(); + Real lam = dsig.L2norm() / denominator_tensor.L2norm(); std::vector dlam_dsig(3); for (unsigned i = 0; i < _num_sp; ++i) dlam_dsig[i] = dsig.L2norm() == 0. ? 0. : -(trial_stress_params[i] - current_stress_params[i]) / - (dsig.L2norm() * fac.L2norm()); + (dsig.L2norm() * denominator_tensor.L2norm()); std::vector h(2); hardPotential(current_stress_params, intnl, h); @@ -545,16 +521,6 @@ DamagePlasticityStressUpdate::setIntnlDerivativesV(const std::vector & tri dintnl[i][j] = dlam_dsig[j] * h[i] + lam * dh_dsig[i][j]; } -// Real -// DamagePlasticityStressUpdate::ftbar(const std::vector & intnl) const -// { -// Real sqrtPhi_t = std::sqrt(1. + _at * (2. + _at) * intnl[0]); -// if (_zt > sqrtPhi_t / _at) -// return _ft0 * std::pow(_zt - sqrtPhi_t / _at, (1. - _dt_bt)) * sqrtPhi_t; -// else -// return 1.E-6; // A very small number (instead of zero) is used for end of softening -// } - Real DamagePlasticityStressUpdate::fbar(const Real & f0, const Real & a, @@ -565,8 +531,8 @@ DamagePlasticityStressUpdate::fbar(const Real & f0, Real sqrt_phi = std::sqrt(phi); Real v = sqrt_phi; Real u = (1 + a) / a - sqrt_phi / a; - Real fbar = f0 * std::pow(u, exponent) * v; - return (u > 0) ? fbar : 1.E-6; + Real cal_fbar = f0 * std::pow(u, exponent) * v; + return (u > 0) ? cal_fbar : 1.E-6; } Real @@ -581,9 +547,11 @@ DamagePlasticityStressUpdate::dfbar_dkappa(const Real & f0, Real v = sqrt_phi; Real u = (1 + a) / a - sqrt_phi / a; Real dv_dphi = 1. / (2 * v); - Real du_dphi = -(1 / (2 * a)) * exponent * std::pow(u, exponent - 1) * (1 / v); - Real dfbar_dkappa = f0 * (u * dv_dphi + v * du_dphi) * dphi_dkappa; - return (u > 0) ? dfbar_dkappa : 0.; + Real du_dphi = -(1 / (2 * a)) * (1 / sqrt_phi); + Real cal_dfbar_dkappa = + f0 * (std::pow(u, exponent) * dv_dphi + exponent * std::pow(u, exponent - 1) * v * du_dphi) * + dphi_dkappa; + return (u > 0) ? cal_dfbar_dkappa : 0.; } Real @@ -593,8 +561,8 @@ DamagePlasticityStressUpdate::f(const Real & f0, const Real & a, const Real & ka Real sqrt_phi = std::sqrt(phi); Real v = phi; Real u = (1 + a) * sqrt_phi; - Real f = (f0 / a) * (u - v); - return (u > v) ? f : 0.; + Real cal_f = (f0 / a) * (u - v); + return (u > v) ? cal_f : 1.E-6; } Real @@ -607,68 +575,9 @@ DamagePlasticityStressUpdate::df_dkappa(const Real & f0, const Real & a, const R Real u = (1 + a) * sqrt_phi; Real dv_dphi = 1.; Real du_dphi = (1 + a) / (2 * sqrt_phi); - Real df_dkappa = (f0 / a) * (du_dphi - dv_dphi) * dphi_dkappa; - return (u > v) ? df_dkappa : 0.; -} - -// Real -// DamagePlasticityStressUpdate::dftbar(const std::vector & intnl) const -// { -// Real sqrtPhi_t = std::sqrt(1. + _at * (2. + _at) * intnl[0]); -// if (_zt > sqrtPhi_t / _at) -// return _ft0 * _dPhit / (2 * sqrtPhi_t) * std::pow(_zt - sqrtPhi_t / _at, -_dt_bt) * -// (_zt - (2. - _dt_bt) * sqrtPhi_t / _at); -// else -// return 0.; -// } - -// Real -// DamagePlasticityStressUpdate::fcbar(const std::vector & intnl) const -// { -// Real sqrtPhi_c; -// if (intnl[1] < 1.0) -// sqrtPhi_c = std::sqrt(1. + _ac * (2. + _ac) * intnl[1]); -// else -// sqrtPhi_c = std::sqrt(1. + _ac * (2. + _ac) * 0.99); -// return _fc0 * std::pow((_zc - sqrtPhi_c / _ac), (1. - _dc_bc)) * sqrtPhi_c; -// } - -// Real -// DamagePlasticityStressUpdate::dfcbar(const std::vector & intnl) const -// { -// if (intnl[1] < 1.0) -// { -// Real sqrtPhi_c = std::sqrt(1. + _ac * (2. + _ac) * intnl[1]); -// return _fc0 * _dPhic / (2. * sqrtPhi_c) * std::pow(_zc - sqrtPhi_c / _ac, -_dc_bc) * -// (_zc - (2. - _dc_bc) * sqrtPhi_c / _ac); -// } -// else -// return 0.; -// } - -// Real -// DamagePlasticityStressUpdate::ft(const std::vector & intnl) const -// { -// Real Phi_t = 1. + _at * (2. + _at) * intnl[0]; -// Real sqrtPhi_t = std::sqrt(Phi_t); -// if (_zt * sqrtPhi_t > Phi_t / _at) -// return _ft0 * (_zt * sqrtPhi_t - Phi_t / _at); -// else -// return 1.E-6; // A very small number (instead of zero) is used for end of softening -// } - -// Real -// DamagePlasticityStressUpdate::fc(const std::vector & intnl) const -// { -// Real Phi_c, sqrtPhi_c; - -// if (intnl[1] < 1.0) -// Phi_c = 1. + _ac * (2. + _ac) * intnl[1]; -// sqrtPhi_c = std::sqrt(Phi_c); -// else Phi_c = 1. + _ac * (2. + _ac) * 0.999999; -// sqrtPhi_c = std::sqrt(Phi_c); -// return _fc0 * (_zc * sqrtPhi_c - Phi_c / _ac); -// } + Real cal_df_dkappa = (f0 / a) * (du_dphi - dv_dphi) * dphi_dkappa; + return (u > v) ? cal_df_dkappa : 0.; +} Real DamagePlasticityStressUpdate::beta(const std::vector & intnl) const From 6e8b1a211aa4d121dc4b13c71855fbd7d4606048 Mon Sep 17 00:00:00 2001 From: Prithivirajan Veerappan Date: Tue, 2 Jul 2024 21:21:18 -0600 Subject: [PATCH 5/7] Final cleanup and gold files update ref #184 --- .../materials/DamagePlasticityStressUpdate.md | 125 +++++++++--------- .../materials/DamagePlasticityStressUpdate.h | 72 +--------- src/materials/DamagePlasticityStressUpdate.C | 22 +-- .../gold/dilatancy_out.csv | 90 ++++++------- .../gold/shear_test_out.csv | 56 ++++---- .../gold/uniaxial_compression_out.csv | 90 ++++++------- .../gold/uniaxial_tension_2elem_out.csv | 25 ---- .../gold/uniaxial_tension_out.csv | 98 +++++++------- 8 files changed, 239 insertions(+), 339 deletions(-) delete mode 100644 test/tests/damage_plasticity_model/gold/uniaxial_tension_2elem_out.csv diff --git a/doc/content/source/materials/DamagePlasticityStressUpdate.md b/doc/content/source/materials/DamagePlasticityStressUpdate.md index 689cd492..714a149a 100644 --- a/doc/content/source/materials/DamagePlasticityStressUpdate.md +++ b/doc/content/source/materials/DamagePlasticityStressUpdate.md @@ -2,22 +2,24 @@ The [!cite](lee1996theory) model accounts for the independent damage in tension and compression. It also accounts for degradation of the elastic modulus of the concrete as the loading goes beyond yielding in either tension or compression. The model uses the incremental theory of plasticity and decomposes the total strain, $\boldsymbol{\varepsilon}$, into elastic strain, $\boldsymbol{\varepsilon}^{e}$, and plastic strain, $\boldsymbol{\varepsilon}^{p}$, as follows \begin{equation} - \boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}^{e} + \boldsymbol{\varepsilon}^{p} \label{eps_def} + \label{straindecomposition} + \boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}^{e} + \boldsymbol{\varepsilon}^{p} \end{equation} where bold symbol represents a vectoral or tensorial quantity. The relation between elastic strain and the stress, $\boldsymbol{\sigma}$, is given by \begin{equation} - \boldsymbol{\varepsilon}^{e} = \boldsymbol{\mathfrak{E}}^{-1}:\boldsymbol{\sigma} \label{eps_e_def} + \label{elasticstrain} + \boldsymbol{\varepsilon}^{e} = \boldsymbol{\mathfrak{E}}^{-1}:\boldsymbol{\sigma} \end{equation} -where $\boldsymbol{\mathfrak{E}}$ is the elasticity tensor. Using Eqs. \eqref{eps_def}-\eqref{eps_e_def}, the relation between $\boldsymbol{\sigma}$ and $\boldsymbol{\varepsilon}^{p}$ is expressed as +where $\boldsymbol{\mathfrak{E}}$ is the elasticity tensor. Using [straindecomposition] and [elasticstrain], the relation between $\boldsymbol{\sigma}$ and $\boldsymbol{\varepsilon}^{p}$ is expressed as \begin{equation} - \boldsymbol{\sigma} = \boldsymbol{\mathfrak{E}}:\left(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{e}\right) + \boldsymbol{\sigma} = \boldsymbol{\mathfrak{E}}:\left(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{p}\right) \end{equation} Since the model considers the effect of damage in elastic stiffness, an effective stress, $\boldsymbol{\sigma}^{e}$, is defined, where the stress for a given strain always corresponds to the undamaged elastic stiffness of the material, $\boldsymbol{\mathfrak{E}}_{0}$ The relation between $\boldsymbol{\sigma}^{e}$, $\boldsymbol{\varepsilon}$, and $\boldsymbol{\varepsilon}^{p}$ is given by \begin{equation} - \boldsymbol{\sigma}^e = \boldsymbol{\mathfrak{E}}_0:\left(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{e}\right) + \boldsymbol{\sigma}^e = \boldsymbol{\mathfrak{E}}_0:\left(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{p}\right) \end{equation} To consider the degradation of reinforced-concrete structures, an isotropic damage was considered in concrete material. Hence, the relation between $\boldsymbol{\sigma}^e$ and $\boldsymbol{\sigma}$ can be established by @@ -26,7 +28,8 @@ the isotropic scalar degradation damage variable, D, as follows \boldsymbol{\sigma} = \left(1-D\right)\boldsymbol{\sigma}^e \label{sigma_def} \end{equation} \begin{equation} - \boldsymbol{\sigma} = \left(1-D\right)\boldsymbol{\mathfrak{E}}_0:\left(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{e}\right)\label{sigma_def2} + \label{sigma_def2} + \boldsymbol{\sigma} = \left(1-D\right)\boldsymbol{\mathfrak{E}}_0:\left(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{e}\right) \end{equation} The Damage Plasticity Model has various attributes to define the mechanical behavior of concrete in tension and compression such as the yield function, plastic potential, strength of material @@ -38,9 +41,10 @@ sections. ## Yield Function The yield function, $\mathfrak{F}$ is a function of $\boldsymbol{\sigma}$, the strength of the material in uniaxial tension, $f_t$, and the strength of the material in uniaxial compression, $f_c$. It was used to describe the admissible stress space. For this implementation, the yield function in stress space is defined as follows -\begin{equation} \label{yf} +\begin{equation} +\label{yf} \begin{gathered} - \mathfrak{F}\left(\boldsymbol{\sigma},f_t,f_c\right) = \frac{1}{1-\alpha} \\ + \mathfrak{F}\left(\boldsymbol{\sigma},f_t,f_c\right) = \frac{1}{1-\alpha} \left(\alpha I_1 + \sqrt{3J_2} + \beta\left(\boldsymbol{\kappa}\right)<{\hat{\boldsymbol{\sigma}}_{max}}>\right) - f_c\left(\boldsymbol{\kappa}\right) \end{gathered} \end{equation} @@ -53,26 +57,30 @@ relates tensile, $f_t\left(\boldsymbol{\kappa}\right)$, and compressive, $f_c\le function of a vector of damage variable, $\boldsymbol{\kappa} = \{\kappa_t, \kappa_c\}$ and $\kappa_t$ and $\kappa_c$ are the damage variables in tension and compression, respectively. -The implementation first solves the given problem in the effective stress space and then transform the effective stress to stress space using Eq. \eqref{sigma_def2}. Thus, the yield strength of the concrete under uniaxial loading is expressed as effective yield strength as follows +The implementation first solves the given problem in the effective stress space and then transform the effective stress to stress space using [sigma_def2]. Thus, the yield strength of the concrete under uniaxial loading is expressed as effective yield strength as follows \begin{equation} - f_t\left(\boldsymbol{\kappa}\right) = \left(1-D_t \left(\kappa_t\right)\right)f_{t}^{e}\left(\kappa_t\right) \label{ft} + \label{ft} + f_t\left(\boldsymbol{\kappa}\right) = \left(1-D_t \left(\kappa_t\right)\right)f_{t}^{e}\left(\kappa_t\right) \end{equation} \begin{equation} - f_c\left(\boldsymbol{\kappa}\right) = \left(1-D_c \left(\kappa_c\right)\right)f_{c}^{e}\left(\kappa_c\right) \label{fc} + \label{fc} + f_c\left(\boldsymbol{\kappa}\right) = \left(1-D_c \left(\kappa_c\right)\right)f_{c}^{e}\left(\kappa_c\right) \end{equation} where $f_{t}^{e}$ and $f_{c}^{e}$ are the yield strength of the concrete in tension and compression, respectively and $D_t$ and $D_c$ are the degradation damage variables in -tension and compression, respectively such that $0\leq D_t$\textless 1 and $0\leq D_c$\textless 1. +tension and compression, respectively such that $0\leq D_t\leq 1$ and $0\leq D_c\leq 1$. The scalar degradation damage variable is expressed in terms of $D_t$ and $D_c$ as follows \begin{equation} D\left(\boldsymbol{\kappa}\right) = 1-\left(1-D_t\left(\kappa_t\right)\right)\left(1-D_c\left(\kappa_c\right)\right) \label{D} \end{equation} -Hence, for uniaxial tension, $D=D_t$, while for uniaxial compression, $D=D_c$.The yield strength for multi-axial loading, i.e., Eqs. \eqref{ft}-\eqref{fc}, can be rewritten as +Hence, for uniaxial tension, $D=D_t$, while for uniaxial compression, $D=D_c$. The yield strength for multi-axial loading, i.e., [ft] and [fc], can be rewritten as \begin{equation} - f_t\left(\boldsymbol{\kappa}\right) = \left(1-D\left(\boldsymbol{\kappa}\right)\right)f_{t}^{e}\left(\kappa_t\right) \label{ft_new} + \label{ft_new} + f_t\left(\boldsymbol{\kappa}\right) = \left(1-D\left(\boldsymbol{\kappa}\right)\right)f_{t}^{e}\left(\kappa_t\right) \end{equation} \begin{equation} - f_c\left(\boldsymbol{\kappa}\right) = \left(1-D\left(\boldsymbol{\kappa}\right)\right)f_{c}^{e}\left(\kappa_c\right) \label{fc_new} + \label{fc_new} + f_c\left(\boldsymbol{\kappa}\right) = \left(1-D\left(\boldsymbol{\kappa}\right)\right)f_{c}^{e}\left(\kappa_c\right) \end{equation} Similarly, the first invariant of $\boldsymbol{\sigma}^e$, $I_1^e$, and second invariant of the deviatoric component of $\boldsymbol{\sigma}^e$, $J_2^e$, can be rewritten in terms of $I_1$ and $J_2$ as follows \begin{equation} @@ -81,13 +89,15 @@ Similarly, the first invariant of $\boldsymbol{\sigma}^e$, $I_1^e$, and second i \begin{equation} J_2^e = \left(1-D\left(\boldsymbol{\kappa}\right)\right)^2J_2 \label{J2e} \end{equation} -Since $D$ \textless 1, the maximum principal effective stress ${\hat{\boldsymbol{\sigma}}_{max}}^e$ is expressed in the terms of ${\hat{\boldsymbol{\sigma}}_{max}}$ as follows +The maximum principal effective stress ${\hat{\boldsymbol{\sigma}}_{max}}^e$ is expressed in the terms of ${\hat{\boldsymbol{\sigma}}_{max}}$ as follows \begin{equation} - {\hat{\boldsymbol{\sigma}}_{max}}^e = \left(1-D\left(\boldsymbol{\kappa}\right)\right){\hat{\boldsymbol{\sigma}}_{max}} \label{sig_max_e} + \label{sig_max_e} + {\hat{\boldsymbol{\sigma}}_{max}}^e = \left(1-D\left(\boldsymbol{\kappa}\right)\right){\hat{\boldsymbol{\sigma}}_{max}} \end{equation} Consequently, yield function $\left(\mathfrak{F}\left(\boldsymbol{\sigma},f_t,f_c\right)\right)$ is a homogenous -function, i.e., $x \mathfrak{F}\left(\boldsymbol{\sigma},f_t,f_c\right) = \mathfrak{F}\left(x \boldsymbol{\sigma},x f_t,x f_c\right)$ Hence, using Eqs. \eqref{ft_new}-\eqref{sig_max_e}, the yield function in the effective stress space was obtained by multiplying by a factor $\left(1-D\right)$ of both sides of Eq. \eqref{yf}, as follows -\begin{equation}\label{yf_e} +function, i.e., $x \mathfrak{F}\left(\boldsymbol{\sigma},f_t,f_c\right) = \mathfrak{F}\left(x \boldsymbol{\sigma},x f_t,x f_c\right)$ Hence, using [ft_new] and [sig_max_e], the yield function in the effective stress space was obtained by multiplying by a factor $\left(1-D\right)$ of both sides of [yf], as follows +\begin{equation} +\label{yf_e} \begin{gathered} \mathfrak{F}\left(\boldsymbol{\sigma}^e,f_t^e,f_c^e\right) = \frac{1}{1-\alpha} \\ \left(\alpha I_1^e + \sqrt{3J_2^e} + \beta\left(\boldsymbol{\kappa}\right)<{\hat{\boldsymbol{\sigma}}_{max}}^e>\right) - f_c^e\left(\boldsymbol{\kappa}\right) @@ -111,7 +121,8 @@ dilatancy of concrete, and $\dot{\gamma}$ is the plastic consistency parameter. Since the concrete shows strain-softening in tension and strain hardening and softening in compression, the concrete strength is expressed as a combination of two exponential functions as follows \begin{equation} - f_N = f_{N0} \left(\left(1+a_N\right) e^{-b_N \varepsilon^p}- a_N e^{-2b_N \varepsilon^p}\right) \label{fN} +\label{fN} + f_N = f_{N0} \left(\left(1+a_N\right) e^{-b_N \varepsilon^p}- a_N e^{-2b_N \varepsilon^p}\right) \end{equation} where $f_{N0}$ is the initial yield stress of the material, $N = t$, for the uniaxial tension, $N = c$, for uniaxial compression, $a_N$ and $b_N$, are the material constants @@ -119,12 +130,14 @@ that describe the softening and hardening behavior of the concrete. Similarly, t degradation of the elastic modulus is also expressed as another exponential function as follows \begin{equation} - D_N = 1 - e^{-d_N \varepsilon^p} \label{DN} +\label{DN} + D_N = 1 - e^{-d_N \varepsilon^p} \end{equation} where $d_N$ is a constant that determine the rate of degradation of $\boldsymbol{\mathfrak{E}}$ with the increase in plastic strain. The strength of the material in the effective stress space was -obtained using Eqs. \eqref{ft_new}-\eqref{fc_new}, and \eqref{fN}-\eqref{DN}, as follows -\begin{equation}\label{fNe} +obtained using [ft_new], [fc_new], [fN], and [DN], as follows +\begin{equation} + \label{fNe} f_N^e = f_{N0} \left(\left(1+a_N\right) \left(e^{-b_N \varepsilon^p}\right)^{1-\frac{d_N}{b_N}}- a_N \left(e^{-b_N \varepsilon^p}\right)^{2-\frac{d_N}{b_N}}\right) \end{equation} @@ -144,7 +157,7 @@ Thus, the plastic strain can be presented in terms of damage variable as follows \begin{equation} \varepsilon^p = \frac{1}{b_N} \log{\frac{\sqrt{\Phi_N}}{a_N}} \label{eps_p} \end{equation} -where $\Phi_N = 1 + a_N \left(2+a_N \right)\kappa_N$. Using Eqs. \eqref{fN} and \eqref{eps_p}, the +where $\Phi_N = 1 + a_N \left(2+a_N \right)\kappa_N$. Using [fN] and [eps_p], the strength of the concrete can be expressed in terms of the damage variable as follows \begin{equation} f_N = f_{N0} \frac{1+a_N-\sqrt{\Phi_N\left(\kappa_N\right)}}{a_N}\sqrt{\Phi_N\left(\kappa_N\right)} \label{fN_new} @@ -160,7 +173,7 @@ Thus, the strength of the material and degradation damage variable in the effect where $a_N$, $b_N$,and $d_N$ are the modeling parameters, which are evaluated from the material properties. Since the maximum compressive strength of concrete, $f_{cm}$, was used as a material property, $f_{cm}$ was obtained in terms of $a_c$ by finding maximum value of -compressive strength in Eq. \eqref{fNe} as follows +compressive strength in [fNe] as follows \begin{equation} f_{cm} = \frac{f_{c0}\left(1+a_c\right)^2}{4a_c} \label{fcm} \end{equation} @@ -182,21 +195,25 @@ as follows \begin{equation} \left(\frac{d\sigma}{d\varepsilon^p}\right)_{\varepsilon^p=0} = f_{t0}b_t\left(a_t-1\right) \label{slope} \end{equation} -Thus, $a_t$ was obtained using Eqs. \eqref{bt}-\eqref{slope} as follows +Thus, $a_t$ was obtained using [bt]-[slope] as follows \begin{equation} - a_t = \sqrt{\frac{9}{4}+\frac{2\frac{G_t}{l_t} \left(\frac{d\sigma}{d\varepsilon^p}\right)_{\varepsilon^p=0}}{f_{t0}^2}}\label{at} + \label{a_t} + a_t = \sqrt{\frac{9}{4}+\frac{2\frac{G_t}{l_t} \left(\frac{d\sigma}{d\varepsilon^p}\right)_{\varepsilon^p=0}}{f_{t0}^2}} \end{equation} -The minimum slope of the $\sigma$ versus $\varepsilon^p$ curve is +To obtain a real value of $a_t$, the quantity inside the square root must be $\geq$ 0. Therefore, the minimum possible slope of the $\sigma$ versus $\varepsilon^p$ curve is $\left(\left(\frac{d\sigma}{d\varepsilon^p}\right)_{\varepsilon^p=0}\right)_{min}= -\frac{9}{8}\frac{f_{t0}^2}{\frac{G_t}{l_t}}$, which is a function of the characteristic length in tension. Therefore, a mesh independent slope parameter $\omega\in\left(0,1\right)$, is defined such that \begin{equation} \left(\frac{d\boldsymbol{\sigma}}{d\varepsilon^p}\right)_{\varepsilon^p=0} = \omega \left(\left(\frac{d\sigma}{d\varepsilon^p}\right)_{\varepsilon^p=0}\right)_{min} \label{slope_new} \end{equation} -Using Eqs. \eqref{at}-\eqref{slope_new}, $a_t$ is rewritten as follows +Using [a_t] and [slope_new], $a_t$ is rewritten as follows \begin{equation} a_t = \frac{3}{2}\sqrt{1-\omega}-\frac{1}{2}\label{at_new} \end{equation} + +Note that $\omega$ is a fitting parameter that must be provided by the user. + The ratio of $\frac{d_c}{b_c}$ was obtained by specifying degradation values for uniaxial compression case from experiments. If the degradation in the elastic modulus is known, denoted as $\widetilde{D}_c$, when the concrete is unloaded from $\sigma =f_{cm}$, then $\frac{d_c}{b_c}$ will be obtained using the following relation @@ -214,9 +231,7 @@ Similarly, if degradation in the elastic modulus is known, denoted as $\widetild \frac{d_t}{b_t} = \frac{\log\left(1-\widetilde{D}_t\right)}{\log\left(\frac{1+a_t-\sqrt{1+a_t^2}}{2a_t}\right)} \label{Dt_ft0} \end{equation} Thus, material modeling parameters $a_N$,$b_N$, and $d_N$ were obtained, which were used in -defining the strength of concrete in both tension and compression as given in Eq. -\eqref{fNe_new}. These parameters are also used to define the degradation damage variable in -both tension and compression as indicated in Eq. \eqref{DN_new}. +defining the strength of concrete in both tension and compression as given in [fNe_new]. These parameters are also used to define the degradation damage variable in both tension and compression as indicated in [DN_new]. ## Hardening Potential @@ -226,7 +241,7 @@ The vector of two damage variables, $\boldsymbol{\kappa}=\{\kappa_t, \kappa_c\}$ \end{equation} The evolution of the damage variable is expressed in terms of the evolution of $\boldsymbol{\varepsilon}^p$ as follows \begin{equation} - \dot{\boldsymbol{\kappa}} = \frac{1}{g_N}f_N^e\left(\kappa_N\right)\dot{\boldsymbol{\varepsilon}^p} \label{kappa_ep} + \dot{\boldsymbol{\kappa}} = \frac{1}{g_N}f_N\left(\kappa_N\right)\dot{\boldsymbol{\varepsilon}^p} \label{kappa_ep} \end{equation} where $g_N$ is dissipated energy density during the process of cracking. The scalar $\dot{\boldsymbol{\varepsilon}^p}$, is extended to multi-dimensional case as follows \begin{equation} @@ -241,7 +256,7 @@ where $\delta_{ij}$ is the Dirac delta function and $\hat{\boldsymbol{\sigma}^e} \end{cases} \end{equation} $\dot{\varepsilon}^{p}_{max}$ and $\dot{\varepsilon}^{p}_{min}$ -are the maximum and minimum principal plastic strain, respectively. From Eqs. \eqref{kappa_ep} - \eqref{r_sige}, the evolution of $\boldsymbol{\kappa}$ was obtained as +are the maximum and minimum principal plastic strain, respectively. From [kappa_ep] - [r_sige], the evolution of $\boldsymbol{\kappa}$ was obtained as \begin{equation} \dot{\boldsymbol{\kappa}} = \boldsymbol{h}\left(\hat{\boldsymbol{\sigma}^e}\right):\dot{\boldsymbol{\varepsilon}}^{\hat{p}} \label{kappa_h_ep} \end{equation} @@ -249,12 +264,12 @@ where \begin{equation}\label{h} \boldsymbol{h}\left(\hat{\boldsymbol{\sigma}^e}\right)= \begin{bmatrix} - \frac{r\left(\hat{\boldsymbol{\sigma}^e}\right)}{g_t}f_t^e\left(\kappa_t\right)&0&0\\ + \frac{r\left(\hat{\boldsymbol{\sigma}^e}\right)}{g_t}f_t\left(\kappa_t\right)&0&0\\ 0&1&0\\ - 0&0&\frac{1-r\left(\hat{\boldsymbol{\sigma}^e}\right)}{g_c}f_c^e\left(\kappa_c\right)\\ + 0&0&\frac{1-r\left(\hat{\boldsymbol{\sigma}^e}\right)}{g_c}f_c\left(\kappa_c\right)\\ \end{bmatrix} \end{equation} -and ‘:’ represents products of two matrices. Hence, $H\left(\boldsymbol{\sigma}^e,\boldsymbol{\kappa}\right)$ in Eq. \eqref{kappa} was obtained as follows +and ‘:’ represents products of two matrices. Hence, $H\left(\boldsymbol{\sigma}^e,\boldsymbol{\kappa}\right)$ in [kappa] was obtained as follows \begin{equation} H\left(\boldsymbol{\sigma}^e, \boldsymbol{\kappa}\right) = \boldsymbol{h}\cdot \nabla_{\hat{\boldsymbol{\sigma}^e}}\Phi\left(\hat{\boldsymbol{\sigma}^e}\right) \label{H_def} \end{equation} @@ -292,11 +307,11 @@ During the plastic corrector step, the returned effective stress should satisfy \mathfrak{F}\left(\boldsymbol{\sigma}^e,f_t^e,f_c^e\right) = 0 \end{split} \end{equation} -As per flow rule in Eq. \eqref{flowRule}, the plastic corrector step, i.e., Eq. \eqref{plasticCorrector} can be rewritten as +As per flow rule in [flowRule], the plastic corrector step, i.e., [plasticCorrector] can be rewritten as \begin{equation} \boldsymbol{\sigma^e}_{n+1} = \boldsymbol{\sigma}_{n+1}^{e^{tr}}-\dot{\gamma}\left(2G\frac{\boldsymbol{s}_{n+1}^e}{\|\boldsymbol{s}_{n+1}^e\|} + 3K\alpha_p\boldsymbol{I}\right) \label{returnMap1} \end{equation} -where $G$ is shear modulus and $K$ is bulk modulus. After separating the volumetric and deviatoric components from Eq. \eqref{returnMap1} following relations can be obtained +where $G$ is shear modulus and $K$ is bulk modulus. After separating the volumetric and deviatoric components from [returnMap1] following relations can be obtained \begin{equation} I_{1|n+1} = I_{1|n+1}^{e^{tr}} - 9K\alpha \alpha_p \dot{\gamma} \label{stressRelation1} \end{equation} @@ -306,36 +321,18 @@ where $G$ is shear modulus and $K$ is bulk modulus. After separating the volumet {\|\boldsymbol{s}^{e}_{n+1}\|} = {\|\boldsymbol{s}_{n+1}^{e^{tr}}\|} - 2G\dot{\gamma} \end{gathered} \end{equation} -Using Eqs. \eqref{stressRelation1} and \eqref{stressRelation2}, Eq. \eqref{returnMap1} can be written as +Using [stressRelation1] and [stressRelation2], [returnMap1] can be written as \begin{equation} \boldsymbol{\sigma}_{n+1}^e = \boldsymbol{\sigma}_{n+1}^{e^{tr}}-\dot{\gamma}\left(2G\frac{\boldsymbol{s}^{e^{tr}}_{n+1}}{\|\boldsymbol{s}_{{n+1}}^{e^{tr}} \|}+ 3K\alpha_p\boldsymbol{I}\right) \label{returnMap2} \end{equation} -In case of plastic deformation, the returned state of stress should lie on the yield surface as per Kuhn-Tucker conditions (Eq. \eqref{khunTuckerConditions}, therefore $\mathfrak{F}\left(\boldsymbol{\sigma}_{n+1}^e,f_t^e,f_c^e\right) = 0$, i.e., -\begin{equation} \label{yfnext} - \begin{gathered} - \alpha I_{1|n+1}^e + \sqrt{3J_{2|n+1}^e} + \beta\left(\boldsymbol{\kappa}\right)<\hat{\boldsymbol{\sigma}}^e_{n+1|max}> \\ - \left(1-\alpha\right)f_c^e\left(\boldsymbol{\kappa}\right) = 0 - \end{gathered} -\end{equation} -Using Eq. \eqref{stressRelation1}, \eqref{stressRelation2}, and \eqref{returnMap2}, Eq. \eqref{yfnext} can be written as -\begin{equation} \label{yfzero} - \begin{gathered} - \alpha\left(I_{1|n+1}^{e^{tr}} - 9K\alpha \alpha_p \dot{\gamma}\right) + - \left(\sqrt{\frac{3}{2}}\|\boldsymbol{s}_{{n+1}}^{e^{tr}}\| - \sqrt{6}G\dot{\gamma}\right)\\+ - \beta\left(\boldsymbol{\kappa}\right)<\hat{\boldsymbol{\sigma}}^e_{n+1|max}> - \left(1-\alpha\right)f_c^e\left(\boldsymbol{\kappa}\right) = 0 - \end{gathered} -\end{equation} -Thus, the plastic multiplier can be by solving Eq. \eqref{yfzero} as -\begin{equation}\label{gammaDef} - \dot{\gamma} = - \begin{cases} - \frac{\alpha I_{1|n+1}^{e^{tr}}+\sqrt{\frac{3}{2}}\|\boldsymbol{s}_{{n+1}}^{e^{tr}}\|-\left(1-\alpha\right)f_c^e\left(\boldsymbol{\kappa}\right)} - {9K \alpha_p + \sqrt{6}G}, & \text{if $\sigma_{m|n+1}^e < 0$}\\ - \frac{\alpha I_{1|n+1}^{e^{tr}}+\sqrt{\frac{3}{2}}\|\boldsymbol{s}_{{n+1}}^{e^{tr}}\|+\beta\left(\boldsymbol{\kappa}\right) \sigma_{m|n+1}^{e^{tr}}-\left(1-\alpha\right)f_c^e\left(\boldsymbol{\kappa}\right)} - {9K \alpha_p + \sqrt{6}G + \beta\left(\boldsymbol{\kappa}\right)\left(2G\frac{s^{e^{tr}}_{m|n+1}}{\|\boldsymbol{s}_{{n+1}}^{e^{tr}} \|}+ 3K\alpha_p\right)}, & \text{otherwise}. - \end{cases} + +$\dot{\gamma}$ is calculated as: + +\begin{equation} \label{plasticparameter} + \dot{\gamma}=\frac{\|\boldsymbol{\sigma}_{n+1}^e - \boldsymbol{\sigma}_{n+1}^{e^{tr}}\|}{\|2G\frac{\boldsymbol{s}^{e^{tr}}_{n+1}}{\|\boldsymbol{s}_{{n+1}}^{e^{tr}} \|}+ 3K\alpha_p\boldsymbol{I}\|} \end{equation} -where $\sigma_{m|n+1}^e$, $\sigma_{m|n+1}^{e^{tr}}$, and $s^{e^{tr}}_{m|n+1}$ are the $m^{th}$ component of the $\hat{\boldsymbol{\sigma}}_{n+1}^e$, $\boldsymbol{\sigma}_{n+1}^{e^{tr}}$, and $\boldsymbol{s}^{e^{tr}}_{n+1}$, respectively, which corresponds to maximum principal effective stress in $\left(n+1\right)^{th}$ step. Eq. \eqref{gammaDef} is solved iteratively. +The plastic parameter, [plasticparameter], is evaluated during each iteration of the return mapping algorithm as the current stress is being updated. !syntax parameters /Materials/DamagePlasticityStressUpdate diff --git a/include/materials/DamagePlasticityStressUpdate.h b/include/materials/DamagePlasticityStressUpdate.h index ebb76b6e..c7cf6b6c 100644 --- a/include/materials/DamagePlasticityStressUpdate.h +++ b/include/materials/DamagePlasticityStressUpdate.h @@ -65,10 +65,6 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate const Real _ac; const Real _zt; const Real _zc; - const Real _dPhit; - const Real _dPhic; - const Real _sqrtPhit_max; - const Real _sqrtPhic_max; const Real _dt_bt; const Real _dc_bc; @@ -118,79 +114,28 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate * @param intnl (Array containing damage states in tension and compression, respectively) * @return value of ft (tensile strength) */ - Real fbar(const Real & f0, - const Real & a, - const Real & exponent, - const Real & kappa) const; + Real fbar(const Real & f0, const Real & a, const Real & exponent, const Real & kappa) const; // Real ftbar(const std::vector & intnl) const; /** * Obtain the partial derivative of the undamaged tensile strength to the damage state * @param intnl (Array containing damage states in tension and compression, respectively) * @return value of dft (partial derivative of the tensile strength to the damage state) */ - Real dfbar_dkappa(const Real & f0, - const Real & a, - const Real & exponent, - const Real & kappa) const; + Real + dfbar_dkappa(const Real & f0, const Real & a, const Real & exponent, const Real & kappa) const; - // /** - // * Obtain the partial derivative of the undamaged tensile strength to the damage state - // * @param intnl (Array containing damage states in tension and compression, respectively) - // * @return value of dft (partial derivative of the tensile strength to the damage state) - // */ - // Real dftbar(const std::vector & intnl) const; - - // /** - // * Obtain the undamaged conpressive strength - // * @param intnl (Array containing damage states in tension and compression, respectively) - // * @return value of fc (compressive strength) - // */ - // Real fcbar(const std::vector & intnl) const; - - // /** - // * Obtain the partial derivative of the undamaged compressive strength to the damage state - // * @param intnl (Array containing damage states in tension and compression, respectively) - // * @return value of dfc - // */ - // Real dfcbar(const std::vector & intnl) const; - // /** // * Obtain the damaged tensile strength // * @param intnl (Array containing damage states in tension and compression, respectively) // * @return value of ft (tensile strength) // */ - Real f(const Real & f0, - const Real & a, - const Real & kappa) const; - - // /** - // * Obtain the partial derivative of the damaged tensile strength to the damage state - // * @param intnl (Array containing damage states in tension and compression, respectively) - // * @return value of dft (partial derivative of the tensile strength to the damage state) - // */ - // Real dft(const std::vector & intnl) const; - - // /** - // * Obtain the damaged compressive strength - // * @param intnl (Array containing damage states in tension and compression, respectively) - // * @return value of fc (compressive strength) - // */ - // Real fc(const std::vector & intnl) const; - - // /** - // * Obtain the partial derivative of the damaged compressive strength to the damage state - // * @param intnl (Array containing damage states in tension and compression, respectively) - // * @return value of dfc - // */ - // Real dfc(const std::vector & intnl) const; + Real f(const Real & f0, const Real & a, const Real & kappa) const; /** * Obtain the partial derivative of the undamaged strength to the damage state * @param intnl (Array containing damage states in tension and compression, respectively) * @return value of dft (partial derivative of the tensile strength to the damage state) */ - Real df_dkappa(const Real & f0, - const Real & a, - const Real & kappa) const; + Real df_dkappa(const Real & f0, const Real & a, const Real & kappa) 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 @@ -200,13 +145,6 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate */ Real beta(const std::vector & intnl) const; - /** - * dbeta_dkappa is a derivative of beta wrt. kappa (plastic damage variable) - * @param intnl (Array containing damage states in tension and compression, respectively) - * @return value of dbeta_dkappa - */ - Real dbeta_dkappa(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) diff --git a/src/materials/DamagePlasticityStressUpdate.C b/src/materials/DamagePlasticityStressUpdate.C index 6a2f4d85..0cfc8883 100644 --- a/src/materials/DamagePlasticityStressUpdate.C +++ b/src/materials/DamagePlasticityStressUpdate.C @@ -130,15 +130,6 @@ DamagePlasticityStressUpdate::initQpStatefulProperties() 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) - // _ele_len[_qp] = (_current_elem->length(0, 1) + _current_elem->length(1, 2)) / 2.; - // else - // _ele_len[_qp] = - // (_current_elem->length(0, 1) + _current_elem->length(1, 2) + _current_elem->length(0, 4)) - // / 3.; _ele_len[_qp] = std::cbrt(_current_elem->volume()); _gt[_qp] = _FEt / _ele_len[_qp]; @@ -346,14 +337,14 @@ DamagePlasticityStressUpdate::dflowPotential_dstress( for (unsigned i = 0; i < 3; ++i) for (unsigned j = 0; j < 3; ++j) dr_dstress[i][j] = 0.0; - } else { for (unsigned i = 0; i < 3; ++i) for (unsigned j = 0; j < 3; ++j) - dr_dstress[i][j] = 0.5 * (std::sqrt(2.0 / J2) * d2J2_dsigi_dsigj(i, j) - - (1 / std::sqrt(2)) * std::pow(J2, -1.5) * dJ2_dsigi[i]*dJ2_dsigi[j]); + dr_dstress[i][j] = + 0.5 * (std::sqrt(2.0 / J2) * d2J2_dsigi_dsigj(i, j) - + (1 / std::sqrt(2)) * std::pow(J2, -1.5) * dJ2_dsigi[i] * dJ2_dsigi[j]); } } @@ -457,7 +448,6 @@ DamagePlasticityStressUpdate::setIntnlValuesV(const std::vector & trial_st RankTwoTensor denominator_tensor = (2 * G / norm_sigmadev_trial) * sigmadev_trial + C3 * Identity_tensor; - RankTwoTensor dsig = RankTwoTensor(trial_stress_params[0] - current_stress_params[0], trial_stress_params[1] - current_stress_params[1], trial_stress_params[2] - current_stress_params[2], @@ -465,8 +455,8 @@ DamagePlasticityStressUpdate::setIntnlValuesV(const std::vector & trial_st 0., 0.); -// Implementing Eqn. (21) of Lee & Fenves, 2001, with backsubstituting eqn. (22) - Real lam = dsig.L2norm() / denominator_tensor.L2norm(); + // Implementing Eqn. (21) of Lee & Fenves, 2001, with backsubstituting eqn. (22) + Real lam = dsig.L2norm() / denominator_tensor.L2norm(); std::vector h(2); hardPotential(current_stress_params, intnl_old, h); @@ -493,7 +483,7 @@ DamagePlasticityStressUpdate::setIntnlDerivativesV(const std::vector & tri Real C3 = 3. * K * _alfa_p; RankTwoTensor denominator_tensor = - 2 * G * (sigmadev_trial / norm_sigmadev_trial) + C3 * Identity_tensor; + 2 * G * (sigmadev_trial / norm_sigmadev_trial) + C3 * Identity_tensor; RankTwoTensor dsig = RankTwoTensor(trial_stress_params[0] - current_stress_params[0], trial_stress_params[1] - current_stress_params[1], trial_stress_params[2] - current_stress_params[2], diff --git a/test/tests/damage_plasticity_model/gold/dilatancy_out.csv b/test/tests/damage_plasticity_model/gold/dilatancy_out.csv index 77c7ca92..83b8a04c 100644 --- a/test/tests/damage_plasticity_model/gold/dilatancy_out.csv +++ b/test/tests/damage_plasticity_model/gold/dilatancy_out.csv @@ -5,48 +5,48 @@ time,displacement_x,e_xx,ep_xx,react_x,s_xx,volumetric_strain 30,-0.0003,-0.0003,0,9.5100000000004,-9.5100000000004,-0.00019198156917977 40,-0.0004,-0.0004,0,12.68,-12.68,-0.00025596723479604 50,-0.0005,-0.0005,0,15.85,-15.85,-0.00031994880546093 -60,-0.0006,-0.0006,-1.1665738232695e-05,18.650196118397,-18.650196118397,-0.00035614282315743 -70,-0.0007,-0.0007,-6.1828854256181e-05,20.158118921705,-20.158118921705,-0.00030064680681963 -80,-0.0008,-0.0008,-0.00011331162648892,21.334432816519,-21.334432816519,-0.00024200461336688 -90,-0.0009,-0.0009,-0.00016614913708567,22.416123211693,-22.416123211693,-0.00018013186395605 -100,-0.001,-0.001,-0.00022037903585689,23.403594041213,-23.403594041213,-0.00011493827993692 -110,-0.0011,-0.0011,-0.00027604130414689,24.29730204034,-24.29730204034,-4.6327966773618e-05 -120,-0.0012,-0.0012,-0.00033317799381014,25.097776879489,-25.097776879489,2.5799966407725e-05 -130,-0.0013,-0.0013,-0.00039183297149043,25.80564063209,-25.80564063209,0.00010155080142793 -140,-0.0014,-0.0014,-0.00045205165978624,26.421626612503,-26.421626612503,0.00018103359893895 -150,-0.0015,-0.0015,-0.00051388076630347,26.946597815556,-26.946597815556,0.00026436055778145 -160,-0.0016,-0.0016,-0.00057736797596275,27.381566352239,-27.381566352239,0.0003516465822615 -170,-0.0017,-0.0017,-0.0006425616575271,27.727706334389,-27.727706334389,0.00044300735151637 -180,-0.0018,-0.0018,-0.00070951052885508,27.986379747075,-27.986379747075,0.00053856047597312 -190,-0.0019,-0.0019,-0.00077826318801117,28.159148034979,-28.159148034979,0.00063842301136541 -200,-0.002,-0.002,-0.00084886769514124,28.247791272165,-28.247791272165,0.00074271083270294 -210,-0.0021,-0.0021,-0.00092137106591391,28.254324170102,-28.254324170102,0.00085153743048116 -220,-0.0022,-0.0022,-0.00099581870620686,28.181011331419,-28.181011331419,0.00096501256609294 -230,-0.0023,-0.0023,-0.0010722537858888,28.030380980245,-28.030380980245,0.0010832407810357 -240,-0.0024,-0.0024,-0.0011507165516958,27.805236645313,-27.805236645313,0.0012063197596193 -250,-0.0025,-0.0025,-0.0012312435818667,27.508666191805,-27.508666191805,0.0013343385512292 -260,-0.0026,-0.0026,-0.0013138669884325,27.144047528988,-27.144047528988,0.0014673756658661 -270,-0.0027,-0.0027,-0.0013986135768529,26.715050269669,-26.715050269669,0.0016054970657726 -280,-0.0028,-0.0028,-0.0014855039769868,26.225632593792,-26.225632593792,0.0017487540861813 -290,-0.0029,-0.0029,-0.0015745517640636,25.680032580799,-25.680032580799,0.0018971813294333 -300,-0.003,-0.003,-0.0016657625931362,25.082753334829,-25.082753334829,0.0020507945882438 -310,-0.0031,-0.0031,-0.0017591333751467,24.438541341851,-24.438541341851,0.0022095888650772 -320,-0.0032,-0.0032,-0.0018546515267932,23.752357675995,-23.752357675995,0.0023735365643793 -330,-0.0033,-0.0033,-0.0019522943910391,23.029341963938,-23.029341963938,0.0025425862783206 -340,-0.0034,-0.0034,-0.0020520284412335,22.2747677891,-22.2747677891,0.0027166599366719 -350,-0.0035,-0.0035,-0.0021538095487603,21.493998405543,-21.493998405543,0.0028956551718606 -360,-0.0036,-0.0036,-0.0022575823253257,20.692422407394,-20.692422407394,0.0030794421478335 -370,-0.0037,-0.0037,-0.0023632805093535,19.875400849711,-19.875400849711,0.003267865088308 -380,-0.0038,-0.0038,-0.0024708273331579,19.048202325047,-19.048202325047,0.0034607430605169 -390,-0.0039,-0.0039,-0.0025801361706607,18.215940108345,-18.215940108345,0.0036578714780258 -400,-0.004,-0.004,-0.0026911114471323,17.383510747769,-17.383510747769,0.0038590242449542 -410,-0.0041,-0.0041,-0.0028036497814008,16.555536432998,-16.555536432998,0.0040639564579886 -420,-0.0042,-0.0042,-0.0029176413247561,15.736313397885,-15.736313397885,0.0042724075817051 -430,-0.0043,-0.0043,-0.0030329712499113,14.929768288804,-14.929768288804,0.0044841049865935 -440,-0.0044,-0.0044,-0.00314952133541,14.139424020279,-14.139424020279,0.0046987677201216 -450,-0.0045,-0.0045,-0.00326717163969,13.368375742274,-13.368375742274,0.0049161108236475 -460,-0.0046,-0.0046,-0.0033858019253512,12.61927840845,-12.61927840845,0.0051358474282963 -470,-0.0047,-0.0047,-0.0035052933649444,11.894346691443,-11.894346691443,0.0053576947783152 -480,-0.0048,-0.0048,-0.0036255298401004,11.195360464619,-11.195360464619,0.0055813760528074 -490,-0.0049,-0.0049,-0.0037463991756173,10.523683243842,-10.523683243842,0.0058066236436187 -500,-0.005,-0.005,-0.0038677940117025,9.8802875310769,-9.8802875310769,0.0060331809640377 +60,-0.0006,-0.0006,-9.1669120668388e-06,18.729408887481,-18.729408887481,-0.00036209417652089 +70,-0.0007,-0.0007,-5.1100405080019e-05,20.480388180155,-20.480388180155,-0.00032619960446634 +80,-0.0008,-0.0008,-9.4379488537372e-05,21.842176200533,-21.842176200533,-0.00028709878756816 +90,-0.0009,-0.0009,-0.00013911416050504,23.082819869485,-23.082819869485,-0.00024452934243235 +100,-0.001,-0.001,-0.00018541745675529,24.199780448795,-24.199780448795,-0.00019822157075589 +110,-0.0011,-0.0011,-0.00023340663909691,25.190720917952,-25.190720917952,-0.00014789561727513 +120,-0.0012,-0.0012,-0.00028320471414772,26.053504111012,-26.053504111012,-9.3257838422622e-05 +130,-0.0013,-0.0013,-0.00033494147858477,26.786206261363,-26.786206261363,-3.3998296411264e-05 +140,-0.0014,-0.0014,-0.00038875445895819,27.387143977919,-27.387143977919,3.021149822291e-05 +150,-0.0015,-0.0015,-0.00044478885987997,27.854934856356,-27.854934856356,9.9718484357636e-05 +160,-0.0016,-0.0016,-0.00050319898547507,28.188536109209,-28.188536109209,0.00017489146115257 +170,-0.0017,-0.0017,-0.00056414738965333,28.387329985863,-28.387329985863,0.00025611909080281 +180,-0.0018,-0.0018,-0.00062780478900611,28.451219898281,-28.451219898281,0.00034380972109638 +190,-0.0019,-0.0019,-0.00069434920301473,28.380736556975,-28.380736556975,0.00043838936953322 +200,-0.002,-0.002,-0.00076396445061699,28.177165326922,-28.177165326922,0.00054029817543122 +210,-0.0021,-0.0021,-0.00083683783532906,27.842690489329,-27.842690489329,0.0006499849189312 +220,-0.0022,-0.0022,-0.00091315682537807,27.380553973888,-27.380553973888,0.00076789914334285 +230,-0.0023,-0.0023,-0.00099310453316238,26.795223576689,-26.795223576689,0.00089448041134865 +240,-0.0024,-0.0024,-0.0010768538250872,26.092562056488,-26.092562056488,0.0010301442880194 +250,-0.0025,-0.0025,-0.0011645599632525,25.279983939445,-25.279983939445,0.0011752648101582 +260,-0.0026,-0.0026,-0.0012563518152646,24.366581592785,-24.366581592785,0.0013301535216756 +270,-0.0027,-0.0027,-0.0013523217001089,23.363199984937,-23.363199984937,0.001495035229127 +280,-0.0028,-0.0028,-0.001452515656565,22.282412851501,-22.282412851501,0.0016700247317505 +290,-0.0029,-0.0029,-0.0015569213245906,21.138421960984,-21.138421960984,0.0018550978156373 +300,-0.003,-0.003,-0.0016654606744881,19.946792804303,-19.946792804303,0.0020500737708264 +310,-0.0031,-0.0031,-0.0017779831789587,18.724049172024,-18.724049172024,0.0022545989241167 +320,-0.0032,-0.0032,-0.0018942641465703,17.487143393173,-17.487143393173,0.0024681424560244 +330,-0.0033,-0.0033,-0.0020140082947401,16.252815970977,-16.252815970977,0.0026900047066953 +340,-0.0034,-0.0034,-0.0021368589938802,15.036913337511,-15.036913337511,0.0029193390247211 +350,-0.0035,-0.0035,-0.0022624126097281,13.853738724886,-13.853738724886,0.0031551858172942 +360,-0.0036,-0.0036,-0.0023902364343348,12.715514313872,-12.715514313872,0.003396515221098 +370,-0.0037,-0.0037,-0.0025198880640689,11.632017261926,-11.632017261926,0.0036422732956947 +380,-0.0038,-0.0038,-0.0026509339358262,10.610422180391,-10.610422180391,0.0038914262781964 +390,-0.0039,-0.0039,-0.002782965089357,9.6553470172073,-9.6553470172072,0.0041429982781065 +400,-0.004,-0.004,-0.0029156089333603,8.7690683071801,-8.7690683071801,0.00439609947789 +410,-0.0041,-0.0041,-0.0030485366039898,7.9518528448446,-7.9518528448446,0.0046499438372609 +420,-0.0042,-0.0042,-0.0031814662031486,7.2023481521898,-7.2023481521898,0.0049038569688302 +430,-0.0043,-0.0043,-0.0033141626550219,6.5179810814301,-6.5179810814301,0.0051572759334462 +440,-0.0044,-0.0044,-0.003446435109213,5.8953274140784,-5.8953274140784,0.0054097431621782 +450,-0.0045,-0.0045,-0.0035781328009164,5.3304303250977,-5.3304303250977,0.0056608966738658 +460,-0.0046,-0.0046,-0.0037091401361174,4.8190585889557,-4.8190585889557,0.0059104584208118 +470,-0.0047,-0.0047,-0.003839371578319,4.356904820951,-4.356904820951,0.0061582221404615 +480,-0.0048,-0.0048,-0.0039687667248337,3.9397297528545,-3.9397297528545,0.0064040416422462 +490,-0.0049,-0.0049,-0.0040972858027194,3.5634612215178,-3.5634612215178,0.0066478200821194 +500,-0.005,-0.005,-0.0042249056959746,3.2242571376929,-3.2242571376929,0.0068895004943903 diff --git a/test/tests/damage_plasticity_model/gold/shear_test_out.csv b/test/tests/damage_plasticity_model/gold/shear_test_out.csv index 671ae84a..e73d2d32 100644 --- a/test/tests/damage_plasticity_model/gold/shear_test_out.csv +++ b/test/tests/damage_plasticity_model/gold/shear_test_out.csv @@ -2,31 +2,31 @@ time,e_xy,ep_xy,s_xy 0,0,0,0 5,7.8182642983415e-06,0,0.66498115657902 10,1.5636528591598e-05,0,1.3299623130426 -15,2.3452074951896e-05,4.8983821438774e-07,1.8597699615448 -20,3.1171096317802e-05,2.3402390354366e-06,1.9957411207585 -25,3.8953715744496e-05,4.3144052368041e-06,2.1193699480148 -30,4.6800867741694e-05,6.3955307618849e-06,2.2428277546278 -35,5.5655695121865e-05,1.2095885208344e-05,2.2884800606222 -40,6.4887705610769e-05,1.8474922611195e-05,2.3236101673253 -45,7.440502468223e-05,2.5198091613366e-05,2.356566968123 -50,8.4187054613461e-05,3.2241992754993e-05,2.3882556131892 -55,9.401159776781e-05,3.9971695561668e-05,2.4101616596486 -60,0.00010338033680737,4.9085263419077e-05,2.4042945560395 -65,0.00011281439236663,5.8362558366631e-05,2.3950062148185 -70,0.00012232009850832,6.7761648832007e-05,2.38432596789 -75,0.00013189086765114,7.7262519498049e-05,2.3725090336921 -80,0.00014152159167585,8.6850076139939e-05,2.3597322824429 -85,0.00015120842100567,9.6513149578552e-05,2.3461300399041 -90,0.00016092874359371,0.00010623041065843,2.3330556963112 -95,0.00017068567253073,0.00011599606146946,2.3198767221114 -100,0.00018047911045853,0.00012580544978692,2.3064133062834 -105,0.00019030664536916,0.00013565398059374,2.2927233230542 -110,0.00020016527362133,0.00014553747877069,2.2789204272634 -115,0.00021003576134283,0.00015544217447879,2.2662878314405 -120,0.00021992865864995,0.00016536932057557,2.2537113938821 -125,0.00022984240675046,0.00017531700685573,2.2411960375094 -130,0.00023977530424017,0.00018528382299711,2.2287385663326 -135,0.00024972616547397,0.00019526828146045,2.2163479234135 -140,0.0002596939415043,0.00020526911796477,2.2040309148733 -145,0.00026967769966424,0.00021528525206936,2.1917922967591 -150,0.00027967660824062,0.00022531575412505,2.1796352296801 +15,2.3452942579244e-05,4.909775947702e-07,1.85972220034 +20,3.1170444537696e-05,2.3446630533557e-06,1.9957658984007 +25,3.8937175168089e-05,4.3177448402497e-06,2.1201180140516 +30,4.6761633253159e-05,6.3934989443464e-06,2.243936989742 +35,5.5582665760452e-05,1.2072106128636e-05,2.2898334885283 +40,6.4772321822859e-05,1.8421679998463e-05,2.324907142761 +45,7.4238096905006e-05,2.5103788886935e-05,2.3575535955192 +50,8.3958852172671e-05,3.2094434251767e-05,2.3887363856108 +55,9.3705845077172e-05,3.9767139045216e-05,2.4097318565037 +60,0.00010298810192046,4.8798317682579e-05,2.4030047870696 +65,0.00011232658518956,5.7977735123915e-05,2.3929423386959 +70,0.00012172662902616,6.726665258888e-05,2.3814088104059 +75,0.00013118235998486,7.6645852258031e-05,2.3686419411157 +80,0.00014068949359589,8.6101114941365e-05,2.3548051286361 +85,0.00015024516716526,9.5622313507776e-05,2.3400146078255 +90,0.00015983110415471,0.00010519139385668,2.3254071517662 +95,0.00016944979956828,0.00011480314553222,2.3105040499412 +100,0.00017910192938049,0.00012445399327563,2.2951463105256 +105,0.00018878612455347,0.00013414054703901,2.2793882927109 +110,0.00019850137767844,0.00014386047434706,2.2632728477185 +115,0.00020823526544992,0.00015360354219724,2.2478069011565 +120,0.00021799171392252,0.0001633701666081,2.2323898751356 +125,0.00022777173485429,0.00017316021549842,2.2168522134588 +130,0.00023757484926882,0.00018297296864154,2.2012137034622 +135,0.00024740068262291,0.00019280793054675,2.1854926435378 +140,0.00025724894758144,0.00020266477573483,2.16970505669 +145,0.00026711942680453,0.00021254330678078,2.1538650970437 +150,0.00027701195853931,0.00022244342128784,2.1379853693726 diff --git a/test/tests/damage_plasticity_model/gold/uniaxial_compression_out.csv b/test/tests/damage_plasticity_model/gold/uniaxial_compression_out.csv index 93bb413f..77b90094 100644 --- a/test/tests/damage_plasticity_model/gold/uniaxial_compression_out.csv +++ b/test/tests/damage_plasticity_model/gold/uniaxial_compression_out.csv @@ -5,48 +5,48 @@ time,displacement_x,e_xx,ep_xx,react_x,s_xx,volumetric_strain 30,-0.0003,-0.0003,0,9.5100000000004,-9.5100000000004,-0.00019198156917977 40,-0.0004,-0.0004,0,12.68,-12.68,-0.00025596723479604 50,-0.0005,-0.0005,0,15.85,-15.85,-0.00031994880546093 -60,-0.0006,-0.0006,-1.1804722007126e-05,18.64579033546,-18.64579033546,-0.00036249125782661 -70,-0.0007,-0.0007,-6.2975427612293e-05,20.122808012171,-20.122808012171,-0.00033355073055652 -80,-0.0008,-0.0008,-0.0001154167741537,21.277335768097,-21.277335768097,-0.0003023022102745 -90,-0.0009,-0.0009,-0.00016916425553668,22.340882112175,-22.340882112175,-0.00026868076341113 -100,-0.001,-0.001,-0.0002242550492167,23.313815876581,-23.313815876581,-0.00023261864528201 -110,-0.0011,-0.0011,-0.00028072790806763,24.196578701063,-24.196578701063,-0.00019404521771726 -120,-0.0012,-0.0012,-0.0003386229880583,24.989699557207,-24.989699557207,-0.0001528872622274 -130,-0.0013,-0.0013,-0.00039798165755841,25.693809552899,-25.693809552899,-0.00010906932520027 -140,-0.0014,-0.0014,-0.000458846287862,26.309656701158,-26.309656701158,-6.2514096568211e-05 -150,-0.0015,-0.0015,-0.00052126002186276,26.838120648805,-26.838120648805,-1.3142827843815e-05 -160,-0.0016,-0.0016,-0.00058526651672606,27.280227348863,-27.280227348863,3.9124202758867e-05 -170,-0.0017,-0.0017,-0.00065090964072649,27.637164790429,-27.637164790429,9.4367394857242e-05 -180,-0.0018,-0.0018,-0.000718233178804,27.910292510559,-27.910292510559,0.00015266623725263 -190,-0.0019,-0.0019,-0.00078728048789792,28.101163262761,-28.101163262761,0.00021410034415714 -200,-0.002,-0.002,-0.00085809403408486,28.211530809999,-28.211530809999,0.00027874745080192 -210,-0.0021,-0.0021,-0.00093071497365328,28.243364506337,-28.243364506337,0.00034668296338003 -220,-0.0022,-0.0022,-0.0010051826581814,28.198860492669,-28.198860492669,0.00041797906124197 -230,-0.0023,-0.0023,-0.0010815340934391,28.080451641813,-28.080451641813,0.00049270371458576 -240,-0.0024,-0.0024,-0.0011598033541029,27.89081552488,-27.89081552488,0.00057091962142053 -250,-0.0025,-0.0025,-0.0012400209582391,27.632879910271,-27.632879910271,0.00065268307068678 -260,-0.0026,-0.0026,-0.0013222132081332,27.309825267881,-27.309825267881,0.000738042743399 -270,-0.0027,-0.0027,-0.0014064015068765,26.925083730841,-26.925083730841,0.0008270384687743 -280,-0.0028,-0.0028,-0.0014926016632254,26.48233397031,-26.48233397031,0.00091969995796282 -290,-0.0029,-0.0029,-0.0015808232004642,25.985491471064,-25.985491471064,0.0010160455438672 -300,-0.003,-0.003,-0.001671068688159,25.438693763315,-25.438693763315,0.0011160809612831 -310,-0.0031,-0.0031,-0.0017633331185362,24.846280273874,-24.846280273874,0.00121979820682 -320,-0.0032,-0.0032,-0.0018576034195297,24.212766856934,-24.212766856934,0.0013271749262602 -330,-0.0033,-0.0033,-0.0019538576499078,23.542811758095,-23.542811758095,0.0014381715405911 -340,-0.0034,-0.0034,-0.0020520654025411,22.841188210973,-22.841188210973,0.0015527349162456 -350,-0.0035,-0.0035,-0.0021521866610541,22.112727957278,-22.112727957278,0.0016707927475175 -360,-0.0036,-0.0036,-0.0022541726183581,21.362290944166,-21.362290944166,0.0017922578233942 -370,-0.0037,-0.0037,-0.0023579652160688,20.594705632046,-20.594705632046,0.0019170256393675 -380,-0.0038,-0.0038,-0.0024634977828491,19.81472524027,-19.81472524027,0.0020449760713206 -390,-0.0039,-0.0039,-0.0025706954950948,19.02697458799,-19.02697458799,0.0021759741236631 -400,-0.004,-0.004,-0.0026794760919484,18.235901478961,-18.235901478961,0.0023098712174732 -410,-0.0041,-0.0041,-0.0027897506588635,17.445730487522,-17.445730487522,0.0024465069783212 -420,-0.0042,-0.0042,-0.0029014250713089,16.660420227761,-16.660420227761,0.0025857104013669 -430,-0.0043,-0.0043,-0.0030144004958485,15.883629423217,-15.883629423217,0.0027273028563726 -440,-0.0044,-0.0044,-0.0031285750309957,15.118687704021,-15.118687704021,0.002871099716808 -450,-0.0045,-0.0045,-0.0032438448492614,14.368572974311,-14.368572974311,0.0030169127559689 -460,-0.0046,-0.0046,-0.0033601054375264,13.635899019672,-13.635899019672,0.00316455240073 -470,-0.0047,-0.0047,-0.0034772527912414,12.922910524265,-12.922910524265,0.0033138299034297 -480,-0.0048,-0.0048,-0.0035951845279389,12.231485835291,-12.231485835291,0.0034645593689513 -490,-0.0049,-0.0049,-0.0037138008917802,11.56314678824,-11.56314678824,0.0036165595854004 -500,-0.005,-0.005,-0.0038330056280617,10.91907465858,-10.91907465858,0.0037696556197346 +60,-0.0006,-0.0006,-9.1667955882632e-06,18.729412579852,-18.729412579852,-0.00036728125827523 +70,-0.0007,-0.0007,-5.1100412364366e-05,20.480389077331,-20.480389077331,-0.00035511418545198 +80,-0.0008,-0.0008,-9.4379540204427e-05,21.842174603733,-21.842174603733,-0.0003405037092421 +90,-0.0009,-0.0009,-0.00013911418877067,23.082818551871,-23.082818551871,-0.00032324995706001 +100,-0.001,-0.001,-0.0001854174728184,24.199779744781,-24.199779744781,-0.00030314738395343 +110,-0.0011,-0.0011,-0.0002334066488026,25.190720539876,-25.190720539876,-0.00027998288402609 +120,-0.0012,-0.0012,-0.00028320472057042,26.05350390551,-26.05350390551,-0.0002535328674077 +130,-0.0013,-0.0013,-0.00033494148333174,26.78620614686,-26.78620614686,-0.00022356135601476 +140,-0.0014,-0.0014,-0.0003887543255599,27.387147629906,-27.387147629906,-0.00018981851751732 +150,-0.0015,-0.0015,-0.00044478888576989,27.854936026822,-27.854936026822,-0.00015203946130693 +160,-0.0016,-0.0016,-0.00050319899537279,28.188535608744,-28.188535608744,-0.00010994432940659 +170,-0.0017,-0.0017,-0.0005641473931745,28.387329865581,-28.387329865581,-6.3236984168658e-05 +180,-0.0018,-0.0018,-0.00062780479144562,28.451219893474,-28.451219893474,-1.1606688089194e-05 +190,-0.0019,-0.0019,-0.00069434920690641,28.380736534535,-28.380736534535,4.5270699125188e-05 +200,-0.002,-0.002,-0.00076396445647247,28.177165242495,-28.177165242495,0.00010773000767905 +210,-0.0021,-0.0021,-0.00083683784250223,27.842690353245,-27.842690353245,0.00017611257340167 +220,-0.0022,-0.0022,-0.00091315683270075,27.380553822188,-27.380553822188,0.00025076022693904 +230,-0.0023,-0.0023,-0.0009931045394922,26.795223448964,-26.795223448964,0.00033200710030856 +240,-0.0024,-0.0024,-0.001076853829627,26.092561981741,-26.092561981741,0.00042016894110053 +250,-0.0025,-0.0025,-0.0011645599658147,25.279983927708,-25.279983927708,0.00051552975297753 +260,-0.0026,-0.0026,-0.0012563518162686,24.366581634074,-24.366581634074,0.00061832582395116 +270,-0.0027,-0.0027,-0.0013523216443961,23.363201026152,-23.363201026152,0.00072872716037731 +280,-0.0028,-0.0028,-0.0014525156531681,22.282413615625,-22.282413615625,0.00084682027487815 +290,-0.0029,-0.0029,-0.0015569213239817,21.138422083796,-21.138422083796,0.00097258512568854 +300,-0.003,-0.003,-0.0016654606729723,19.94679290568,-19.94679290568,0.0011058827886525 +310,-0.0031,-0.0031,-0.0017779831754749,18.724049304057,-18.724049304057,0.0012464425033285 +320,-0.0032,-0.0032,-0.0018942641394318,17.4871435838,-17.4871435838,0.0013938586064373 +330,-0.0033,-0.0033,-0.0020140082816395,16.252816255658,-16.252816255658,0.0015475969149708 +340,-0.0034,-0.0034,-0.0021368589721236,15.036913753372,-15.036913753372,0.0017070114176256 +350,-0.0035,-0.0035,-0.0022624125765996,13.853739302921,-13.853739302921,0.0018713702522481 +360,-0.0036,-0.0036,-0.0023902363875508,12.715515071623,-12.715515071623,0.0020398882370614 +370,-0.0037,-0.0037,-0.0025198880019489,11.632018201222,-11.632018201222,0.0022117620687592 +380,-0.0038,-0.0038,-0.0026509338576133,10.610423286885,-10.610423286885,0.0023862040288534 +390,-0.0039,-0.0039,-0.0027829649954276,9.655348261142,-9.655348261142,0.0025624706791514 +400,-0.004,-0.004,-0.0029156088249252,8.7690696496587,-8.7690696496587,0.0027398843147504 +410,-0.0041,-0.0041,-0.0030485364829342,7.951854243986,-7.951854243986,0.0029178464178032 +420,-0.0042,-0.0042,-0.0031814660717804,7.2023495676928,-7.2023495676928,0.0030958436230102 +430,-0.0043,-0.0043,-0.0033141629696839,6.5179794772875,-6.5179794772875,0.0032734483574832 +440,-0.0044,-0.0044,-0.0034464354416034,5.8953243913278,-5.8953243913278,0.0034503109174754 +450,-0.0045,-0.0045,-0.0035781331462184,5.3304274028279,-5.3304274028279,0.003626156817484 +460,-0.0046,-0.0046,-0.0037091404899291,4.8190558047763,-4.8190558047763,0.0038007747408817 +470,-0.0047,-0.0047,-0.003839371936794,4.3569021994199,-4.3569021994199,0.0039740080444517 +480,-0.0048,-0.0048,-0.0039687670846935,3.939727307777,-3.939727307777,0.0041457458269218 +490,-0.0049,-0.0049,-0.00409728616129,3.5634589580962,-3.5634589580962,0.0043159147532896 +500,-0.005,-0.005,-0.0042249060510952,3.2242550548464,-3.2242550548464,0.004484471838974 diff --git a/test/tests/damage_plasticity_model/gold/uniaxial_tension_2elem_out.csv b/test/tests/damage_plasticity_model/gold/uniaxial_tension_2elem_out.csv deleted file mode 100644 index 8d33b01c..00000000 --- a/test/tests/damage_plasticity_model/gold/uniaxial_tension_2elem_out.csv +++ /dev/null @@ -1,25 +0,0 @@ -time,displacement_x,e_xx,ep_xx,react_x,s_xx,volumetric_strain -0,0,0,0,0,0,0 -1,5e-05,5.0000000001209e-05,0,-1.5850000000321,1.5850000000321,3.2000512005526e-05 -2,0.0001,0.00010000000000003,0,-3.1700000000012,3.1700000000012,6.4002048043743e-05 -3,0.00015,0.00014999999727429,4.0024734223106e-05,-3.4862156892916,3.4862156892916,9.6779051148882e-05 -4,0.0002,0.00017862734293029,6.8093756626783e-05,-3.477471609198,3.477471609198,0.00011216165953351 -5,0.00025,0.00024835041190431,0.000137933469697,-3.4652098234512,3.4652098234512,0.0001616852321851 -6,0.0003,0.00026993896006456,0.00015888007847208,-3.4291886168407,3.4291886168407,0.00017133095578839 -7,0.00035,0.00034485935625532,0.00023401341905145,-3.3921352348948,3.3921352348948,0.00022512601060176 -8,0.0004,0.00035654059152294,0.0002450169275363,-3.35653185387,3.35653185387,0.0002265865480974 -9,0.00045,0.00043833198349069,0.00032705678256184,-3.3201883999425,3.3201883999425,0.00028619595227733 -10,0.0005,0.00044651293525115,0.00033495631513795,-3.2792710308472,3.2792710308472,0.00028450762708399 -11,0.00055,0.00052715305030193,0.00041541055592181,-3.2502722981684,3.2502722981684,0.00034349724648966 -12,0.0006,0.0005332924201731,0.00042146720973885,-3.2069341423991,3.2069341423991,0.00034001250270793 -13,0.00065,0.00060642256236638,0.00049410815746748,-3.1829649589623,3.1829649589623,0.00039312366765298 -14,0.0007,0.00061084634713792,0.00049870146528674,-3.1374831784249,3.1374831784249,0.0003883105411033 -15,0.00075,0.0006645560288364,0.00055144486396568,-3.118863533145,3.118863533145,0.00042589491311074 -16,0.0008,0.00066769219744388,0.00055502022061305,-3.0715413732282,3.0715413732282,0.00042044367651783 -17,0.00085,0.00068831671615083,0.00057436933338461,-3.0558252843467,3.0558252843467,0.0004314913621756 -18,0.0009,0.00069078433673531,0.00057748426642578,-3.0082898057619,3.0082898057619,0.0004264043942068 -19,0.00095,0.00070155032206198,0.00058809635710821,-2.9782137952719,2.9782137952719,0.00042946948977682 -20,0.001,0.00070359512299156,0.0005907971189515,-2.9350739662572,2.9350739662572,0.00042504800025697 -21,0.00105,0.00070982616304493,0.0005969673208269,-2.903402488752,2.903402488752,0.00042476536078784 -22,0.0011,0.00071277747533545,0.00060028276453205,-2.8668857778067,2.8668857778067,0.0004218749400221 -23,0.00115,0.00071611407176282,0.00060388799583745,-2.8403058869124,2.8403058869124,0.00041987152776998 diff --git a/test/tests/damage_plasticity_model/gold/uniaxial_tension_out.csv b/test/tests/damage_plasticity_model/gold/uniaxial_tension_out.csv index 647140fd..16422c59 100644 --- a/test/tests/damage_plasticity_model/gold/uniaxial_tension_out.csv +++ b/test/tests/damage_plasticity_model/gold/uniaxial_tension_out.csv @@ -1,52 +1,52 @@ time,displacement_x,e_xx,ep_xx,react_x,s_xx,volumetric_strain 0,0,0,0,0,0,0 2,0.0001,0.0001,0,-3.1699999999999,3.1699999999999,6.4002048043577e-05 -4,0.0002,0.0002,8.9681412657498e-05,-3.4970986694694,3.4970986694694,0.00012974332005933 -6,0.0003,0.0003,0.00018914297676655,-3.4141368184928,3.4141368184928,0.00019567795022701 -8,0.0004,0.0004,0.00028866489214183,-3.3206475816851,3.3206475816851,0.00026161801110303 -10,0.0005,0.0005,0.00038824436782359,-3.2279013119298,3.2279013119298,0.00032756354008456 -12,0.0006,0.0006,0.00048787872917602,-3.1360782486459,3.1360782486459,0.00039351448492519 -14,0.0007,0.0007,0.00058756540621693,-3.0453370463787,3.0453370463787,0.00045947079552255 -16,0.0008,0.0008,0.00068730185017397,-2.9558165711376,2.9558165711376,0.00052543256741 -18,0.0009,0.0009,0.0007870858579702,-2.8676375604281,2.8676375604281,0.0005913994353004 -20,0.001,0.001,0.00088691507170409,-2.7809040805176,2.7809040805176,0.0006573715348468 -22,0.0011,0.0011,0.00098678730899121,-2.6957052204399,2.6957052204399,0.00072334882337133 -24,0.0012,0.0012,0.0010867003547575,-2.6121161796742,2.6121161796742,0.0007893313770424 -26,0.0013,0.0013,0.0011866524940094,-2.5301995856713,2.5301995856713,0.00085531886808354 -28,0.0014,0.0014,0.0012866419364226,-2.4500065326682,2.4500065326682,0.00092131148098096 -30,0.0015,0.0015,0.0013866658893711,-2.3715780055438,2.3715780055438,0.00098730904650979 -32,0.0016,0.0016,0.0014867235005769,-2.2949456732058,2.2949456732058,0.0010533116629232 -34,0.0017,0.0017,0.0015868128119294,-2.2201320197247,2.2201320197247,0.0011193192527386 -36,0.0018,0.0018,0.0016869322049859,-2.1471525918744,2.1471525918744,0.0011853317842476 -38,0.0019,0.0019,0.0017870801305661,-2.0760158597098,2.0760158597098,0.0012513492271478 -40,0.002,0.002,0.0018872551061383,-2.0067241340424,2.0067241340424,0.0013173715524852 -42,0.0021,0.0021,0.0019874557132744,-1.9392742189882,1.9392742189882,0.0013833987325986 -44,0.0022,0.0022,0.0020876805951768,-1.8736580114884,1.8736580114884,0.001449430741068 -46,0.0023,0.0023,0.0021879284542768,-1.8098630515141,1.8098630515141,0.0015154675526621 -48,0.0024,0.0024,0.0022881980499046,-1.7478730264498,1.7478730264498,0.0015815091432911 -50,0.0025,0.0025,0.0023884881960325,-1.6876682329443,1.6876682329443,0.0016475554899593 -52,0.0026,0.0026,0.002488797763518,-1.6292259997676,1.6292259997676,0.0017136068218411 -54,0.0027,0.0027,0.0025891256604954,-1.5725210801648,1.5725210801648,0.0017796626206013 -56,0.0028,0.0028,0.0026894708562166,-1.5175259836082,1.5175259836082,0.0018457231109581 -58,0.0029,0.0029,0.0027898323620786,-1.4642113133651,1.4642113133651,0.0019117882740176 -60,0.003,0.003,0.0028902092338994,-1.4125460569405,1.4125460569405,0.001977858091754 -62,0.0031,0.0031,0.0029906005700804,-1.3624978488801,1.3624978488801,0.0020439325469717 -64,0.0032,0.0032,0.0030910055039585,-1.3140332106806,1.3140332106806,0.0021100113633035 -66,0.0033,0.0033,0.0031914232255192,-1.2671177626704,1.2671177626704,0.002176095047314 -68,0.0034,0.0034,0.0032918529447778,-1.2217164388913,1.2217164388913,0.0022421833226036 -70,0.0035,0.0035,0.0033922939133948,-1.1777936498391,1.1777936498391,0.002308276174863 -72,0.0036,0.0036,0.0034927454173893,-1.1353134423902,1.1353134423902,0.0023743735904589 -74,0.0037,0.0037,0.0035932067756743,-1.0942396456863,1.0942396456863,0.0024404755564074 -76,0.0038,0.0038,0.0036936773386507,-1.0545359995697,1.0545359995697,0.0025065820603414 -78,0.0039,0.0039,0.0037941564868545,-1.0161662692946,1.0161662692946,0.0025726930904868 -80,0.004,0.004,0.0038946436296585,-0.97909434766253,0.97909434766253,0.0026388086356353 -82,0.0041,0.0041,0.0039951382040248,-0.9432843456424,0.9432843456424,0.0027049286851193 -84,0.0042,0.0042,0.0040956396733079,-0.90870067245746,0.90870067245746,0.0027710532287879 -86,0.0043,0.0043,0.0041961475261056,-0.87530810604835,0.87530810604835,0.0028371822569846 -88,0.0044,0.0044,0.0042966612751568,-0.8430718547509,0.8430718547509,0.0029033157605249 -90,0.0045,0.0045,0.0043971804562844,-0.8119576109642,0.8119576109642,0.0029694537306755 -92,0.0046,0.0046,0.0044977046273811,-0.78193159752422,0.78193159752422,0.0030355961591337 -94,0.0047,0.0047,0.0045982333674382,-0.75296060744272,0.75296060744272,0.0031017430380103 -96,0.0048,0.0048,0.0046987662756136,-0.7250120376196,0.7250120376196,0.0031678943598068 -98,0.0049,0.0049,0.0047993029703398,-0.69805391708855,0.69805391708855,0.0032340501174029 -100,0.005,0.005,0.0048998430884682,-0.67205493031169,0.67205493031169,0.0033002103040352 +4,0.0002,0.0002,8.9758064365535e-05,-3.4946693596017,3.4946693596017,0.00012974458960002 +6,0.0003,0.0003,0.00018928852816498,-3.4242928840277,3.4242928840277,0.00019568051563401 +8,0.0004,0.0004,0.00028886340701436,-3.3448433484916,3.3448433484916,0.00026162164798715 +10,0.0005,0.0005,0.00038848092772624,-3.2658812469577,3.2658812469577,0.00032756795280586 +12,0.0006,0.0006,0.00048813938908643,-3.187523628263,3.187523628263,0.00039351939762144 +14,0.0007,0.0007,0.00058783714795083,-3.1098756378058,3.1098756378058,0.00045947595107942 +16,0.0008,0.0008,0.00068757261768049,-3.0330313773949,3.0330313773949,0.00052543758291046 +18,0.0009,0.0009,0.00078734426658164,-2.9570747026478,2.9570747026478,0.00059140426390103 +20,0.001,0.001,0.00088715061635529,-2.8820799721032,2.8820799721032,0.00065737596586279 +22,0.0011,0.0011,0.00098699024056038,-2.8081127504315,2.8081127504315,0.00072335266160439 +24,0.0012,0.0012,0.0010868617630938,-2.7352304680633,2.7352304680633,0.00078933432490125 +26,0.0013,0.0013,0.0011867638566903,-2.663483039488,2.663483039488,0.00085532093046692 +28,0.0014,0.0014,0.0012866952414446,-2.5929134423978,2.5929134423978,0.00092131245392624 +30,0.0015,0.0015,0.0013866546833579,-2.5235582597798,2.5235582597798,0.00098730887178511 +32,0.0016,0.0016,0.0014866409929105,-2.4554481869726,2.4554481869726,0.0010533101614052 +34,0.0017,0.0017,0.0015866530236584,-2.3886085057012,2.3886085057012,0.0011193163009762 +36,0.0018,0.0018,0.0016866896708536,-2.3230595269688,2.3230595269688,0.0011853272694891 +38,0.0019,0.0019,0.0017867498701298,-2.2588170037684,2.2588170037684,0.0012513430467116 +40,0.002,0.002,0.0018868325961767,-2.195892516883,2.195892516883,0.0013173636131623 +42,0.0021,0.0021,0.0019869368614583,-2.1342938343109,2.1342938343109,0.0013833889500865 +44,0.0022,0.0022,0.0020870617149631,-2.0740252460421,2.0740252460421,0.0014494190394321 +46,0.0023,0.0023,0.0021872062409864,-2.0150878757039,2.0150878757039,0.001515453863826 +48,0.0024,0.0024,0.0022873695579442,-1.9574799704671,1.9574799704671,0.0015814934065503 +50,0.0025,0.0025,0.002387550817218,-1.9011971705387,1.9011971705387,0.0016475376515239 +52,0.0026,0.0026,0.002487749202032,-1.8462327594915,1.8462327594915,0.001713586583276 +54,0.0027,0.0027,0.0025879639263608,-1.7925778966196,1.7925778966196,0.001779640186929 +56,0.0028,0.0028,0.0026881942338685,-1.7402218324434,1.7402218324434,0.0018456984481752 +58,0.0029,0.0029,0.0027884393968784,-1.6891521084257,1.6891521084257,0.00191176135326 +60,0.003,0.003,0.0028886987153729,-1.6393547419025,1.6393547419025,0.0019778288889603 +62,0.0031,0.0031,0.0029889715160234,-1.5908143971777,1.5908143971777,0.0020439010425668 +64,0.0032,0.0032,0.0030892571512496,-1.5435145436733,1.5435145436733,0.0021099778018665 +66,0.0033,0.0033,0.0031895549983077,-1.497437601979,1.497437601979,0.0021760591551234 +68,0.0034,0.0034,0.0032898644584066,-1.4525650785958,1.4525650785958,0.0022421450910641 +70,0.0035,0.0035,0.0033901849558527,-1.4088776901201,1.4088776901201,0.0023082355988588 +72,0.0036,0.0036,0.0034905159372209,-1.3663554775712,1.3663554775712,0.002374330668107 +74,0.0037,0.0037,0.0035908568705531,-1.3249779115262,1.3249779115262,0.0024404302888217 +76,0.0038,0.0038,0.0036912072445822,-1.2847239886809,1.2847239886809,0.0025065344514146 +78,0.0039,0.0039,0.0037915665679817,-1.2455723204238,1.2455723204238,0.0025726431466804 +80,0.004,0.004,0.0038919343686399,-1.2075012139692,1.2075012139692,0.0026387563657861 +82,0.0041,0.0041,0.0039923101929586,-1.1704887465654,1.1704887465654,0.0027048741002531 +84,0.0042,0.0042,0.0040926936051749,-1.1345128332615,1.1345128332615,0.002770996341948 +86,0.0043,0.0043,0.0041930841867059,-1.0995512886826,1.0995512886826,0.0028371230830675 +88,0.0044,0.0044,0.0042934815355163,-1.0655818832414,1.0655818832414,0.0029032543161278 +90,0.0045,0.0045,0.0043938852655068,-1.0325823941797,1.0325823941797,0.002969390033952 +92,0.0046,0.0046,0.004494295005924,-1.0005306518149,1.0005306518149,0.003035530229659 +94,0.0047,0.0047,0.0045947104007906,-0.96940458133726,0.96940458133726,0.0031016748966524 +96,0.0048,0.0048,0.0046951311083556,-0.93918224048461,0.93918224048461,0.0031678240286102 +98,0.0049,0.0049,0.0047955568005633,-0.90984185339939,0.90984185339939,0.0032339776194745 +100,0.005,0.005,0.0048959871625415,-0.88136184095056,0.88136184095056,0.0033001356634412 From fadd3b667965e195a752912b3b4a7756080784e4 Mon Sep 17 00:00:00 2001 From: Prithivirajan Veerappan Date: Tue, 30 Jul 2024 16:31:15 -0600 Subject: [PATCH 6/7] Addressed review comments - Minor changes ref #184 --- .../materials/DamagePlasticityStressUpdate.h | 1 - src/materials/DamagePlasticityStressUpdate.C | 26 ++++++++++++++----- 2 files changed, 19 insertions(+), 8 deletions(-) diff --git a/include/materials/DamagePlasticityStressUpdate.h b/include/materials/DamagePlasticityStressUpdate.h index c7cf6b6c..ee4b6681 100644 --- a/include/materials/DamagePlasticityStressUpdate.h +++ b/include/materials/DamagePlasticityStressUpdate.h @@ -231,7 +231,6 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate * @return dr_dstress (dflowpotential_dstress) */ 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 damage states diff --git a/src/materials/DamagePlasticityStressUpdate.C b/src/materials/DamagePlasticityStressUpdate.C index 0cfc8883..2e6efcb9 100644 --- a/src/materials/DamagePlasticityStressUpdate.C +++ b/src/materials/DamagePlasticityStressUpdate.C @@ -130,6 +130,16 @@ DamagePlasticityStressUpdate::initQpStatefulProperties() 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) + // _ele_len[_qp] = (_current_elem->length(0, 1) + _current_elem->length(1, 2)) / 2.; + // else + // _ele_len[_qp] = + // (_current_elem->length(0, 1) + _current_elem->length(1, 2) + _current_elem->length(0, 4)) + // / 3.; + _ele_len[_qp] = std::cbrt(_current_elem->volume()); _gt[_qp] = _FEt / _ele_len[_qp]; @@ -290,7 +300,7 @@ DamagePlasticityStressUpdate::computeAllQV(const std::vector & stress_para dfcbar; flowPotential(stress_params, intnl, all_q[0].dg); - dflowPotential_dstress(stress_params, intnl, all_q[0].d2g); + dflowPotential_dstress(stress_params, all_q[0].d2g); dflowPotential_dintnl(stress_params, intnl, all_q[0].d2g_di); } @@ -315,9 +325,7 @@ DamagePlasticityStressUpdate::flowPotential(const std::vector & stress_par void DamagePlasticityStressUpdate::dflowPotential_dstress( - const std::vector & stress_params, - const std::vector & /* intnl */, - std::vector> & dr_dstress) const + const std::vector & stress_params, std::vector> & dr_dstress) const { Real J2 = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0, 0, 0) .secondInvariant(); @@ -388,7 +396,7 @@ DamagePlasticityStressUpdate::dhardPotential_dstress(const std::vector & s std::vector r(3); flowPotential(stress_params, intnl, r); std::vector> dr_dsig(3, std::vector(3)); - dflowPotential_dstress(stress_params, intnl, dr_dsig); + dflowPotential_dstress(stress_params, dr_dsig); for (unsigned i = 0; i < _num_sp; ++i) { @@ -522,7 +530,9 @@ DamagePlasticityStressUpdate::fbar(const Real & f0, Real v = sqrt_phi; Real u = (1 + a) / a - sqrt_phi / a; Real cal_fbar = f0 * std::pow(u, exponent) * v; - return (u > 0) ? cal_fbar : 1.E-6; + return (u > 0) ? cal_fbar : 1.E-6; // The minimum value for the strength parameter is 1.E-6, as a + // zero value can cause numerical instabilities due to + // derivatives and use of logarithmic functions. } Real @@ -552,7 +562,9 @@ DamagePlasticityStressUpdate::f(const Real & f0, const Real & a, const Real & ka Real v = phi; Real u = (1 + a) * sqrt_phi; Real cal_f = (f0 / a) * (u - v); - return (u > v) ? cal_f : 1.E-6; + return (u > v) ? cal_f : 1.E-6; // The minimum value for the strength parameter is 1.E-6, as a + // zero value can cause numerical instabilities due to derivatives + // and use of logarithmic functions. } Real From 6e457242ca2bdb3fc34048cbb987da71a38e1594 Mon Sep 17 00:00:00 2001 From: Ben Spencer Date: Wed, 31 Jul 2024 17:53:40 -0600 Subject: [PATCH 7/7] Minor updates to address remaining review comments ref #184 --- .../source/materials/DamagePlasticityStressUpdate.md | 6 +++--- src/materials/DamagePlasticityStressUpdate.C | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/content/source/materials/DamagePlasticityStressUpdate.md b/doc/content/source/materials/DamagePlasticityStressUpdate.md index 714a149a..e8acbb04 100644 --- a/doc/content/source/materials/DamagePlasticityStressUpdate.md +++ b/doc/content/source/materials/DamagePlasticityStressUpdate.md @@ -57,7 +57,7 @@ relates tensile, $f_t\left(\boldsymbol{\kappa}\right)$, and compressive, $f_c\le function of a vector of damage variable, $\boldsymbol{\kappa} = \{\kappa_t, \kappa_c\}$ and $\kappa_t$ and $\kappa_c$ are the damage variables in tension and compression, respectively. -The implementation first solves the given problem in the effective stress space and then transform the effective stress to stress space using [sigma_def2]. Thus, the yield strength of the concrete under uniaxial loading is expressed as effective yield strength as follows +The implementation first solves the given problem in the effective stress space and then transforms the effective stress to stress space using [sigma_def2]. Thus, the yield strength of the concrete under uniaxial loading is expressed as effective yield strength as follows \begin{equation} \label{ft} f_t\left(\boldsymbol{\kappa}\right) = \left(1-D_t \left(\kappa_t\right)\right)f_{t}^{e}\left(\kappa_t\right) @@ -307,11 +307,11 @@ During the plastic corrector step, the returned effective stress should satisfy \mathfrak{F}\left(\boldsymbol{\sigma}^e,f_t^e,f_c^e\right) = 0 \end{split} \end{equation} -As per flow rule in [flowRule], the plastic corrector step, i.e., [plasticCorrector] can be rewritten as +Per the flow rule in [flowRule], the plastic corrector step, i.e., [plasticCorrector] can be rewritten as \begin{equation} \boldsymbol{\sigma^e}_{n+1} = \boldsymbol{\sigma}_{n+1}^{e^{tr}}-\dot{\gamma}\left(2G\frac{\boldsymbol{s}_{n+1}^e}{\|\boldsymbol{s}_{n+1}^e\|} + 3K\alpha_p\boldsymbol{I}\right) \label{returnMap1} \end{equation} -where $G$ is shear modulus and $K$ is bulk modulus. After separating the volumetric and deviatoric components from [returnMap1] following relations can be obtained +where $G$ is the shear modulus and $K$ is the bulk modulus. After separating the volumetric and deviatoric components from [returnMap1] the following relations can be obtained \begin{equation} I_{1|n+1} = I_{1|n+1}^{e^{tr}} - 9K\alpha \alpha_p \dot{\gamma} \label{stressRelation1} \end{equation} diff --git a/src/materials/DamagePlasticityStressUpdate.C b/src/materials/DamagePlasticityStressUpdate.C index 2e6efcb9..ece5c04f 100644 --- a/src/materials/DamagePlasticityStressUpdate.C +++ b/src/materials/DamagePlasticityStressUpdate.C @@ -350,9 +350,9 @@ DamagePlasticityStressUpdate::dflowPotential_dstress( { for (unsigned i = 0; i < 3; ++i) for (unsigned j = 0; j < 3; ++j) - dr_dstress[i][j] = - 0.5 * (std::sqrt(2.0 / J2) * d2J2_dsigi_dsigj(i, j) - - (1 / std::sqrt(2)) * std::pow(J2, -1.5) * dJ2_dsigi[i] * dJ2_dsigi[j]); + dr_dstress[i][j] = 0.5 * (std::sqrt(2.0 / J2) * d2J2_dsigi_dsigj(i, j) - + (1 / std::sqrt(2)) * 1.0 / std::sqrt(Utility::pow<3>(J2)) * + dJ2_dsigi[i] * dJ2_dsigi[j]); } }