From ebf9875489f35da368ceb59810661458cd68a4c2 Mon Sep 17 00:00:00 2001 From: "Mark C. Messner" Date: Fri, 1 Jul 2022 15:23:37 -0500 Subject: [PATCH] Added an option to log transform the critical work relation --- doc/sphinx/damage/work.rst | 8 ++++- include/damage.h | 7 ++++ src/damage.cxx | 35 +++++++++++++----- test/test_damage.py | 74 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 115 insertions(+), 9 deletions(-) diff --git a/doc/sphinx/damage/work.rst b/doc/sphinx/damage/work.rst index c43b89a6..c71365ca 100644 --- a/doc/sphinx/damage/work.rst +++ b/doc/sphinx/damage/work.rst @@ -22,6 +22,11 @@ In principle, the critical work might also depend on temperatures. However, at least for one material (Alloy 617) the temperature dependence is relatively negligible. +If requested the model will first take the log of the work rate +before passing it to the :math:`W_{crit}` function and uses :math:`10^f` of the +returned value. This means the user is providing the function on a log-log +scale. + Parameters ---------- @@ -32,7 +37,8 @@ Parameters ``elastic``, :cpp:class:`neml::LinearElasticModel`, Elasticity model, No ``Wcrit``, :cpp:class:`neml::Interpolate`, Critical work, No ``n``, :code:`double`, Damage exponent, No - ``eps``, :code:`double`, Numerical offset``, ``1.0e-30`` + ``eps``, :code:`double`, Numerical offset, ``1.0e-30`` + ``log``, :code:`bool`, Log transform the work to failure relation, ``false`` Class description ----------------- diff --git a/include/damage.h b/include/damage.h index 9669c22f..5153d737 100644 --- a/include/damage.h +++ b/include/damage.h @@ -653,11 +653,18 @@ class NEML_EXPORT WorkDamage: public ScalarDamage { const double * const stress_np1, const double * const stress_n, double T_np1, double T_n, double t_np1, double t_n, double d_np1, double d_n) const; + + /// Get the critical work, accounting for log scaling if requested + double Wcrit(double Wdot) const; + + /// Get the derivative of the critical work, accounting for log scaling + double dWcrit(double Wdot) const; protected: std::shared_ptr Wcrit_; double n_; double eps_; + bool log_; }; static Register regWorkDamage; diff --git a/src/damage.cxx b/src/damage.cxx index 2a706247..68c0821c 100644 --- a/src/damage.cxx +++ b/src/damage.cxx @@ -1274,7 +1274,8 @@ WorkDamage::WorkDamage(ParameterSet & params) : ScalarDamage(params), Wcrit_(params.get_object_parameter("Wcrit")), n_(params.get_parameter("n")), - eps_(params.get_parameter("eps")) + eps_(params.get_parameter("eps")), + log_(params.get_parameter("log")) { } @@ -1291,6 +1292,7 @@ ParameterSet WorkDamage::parameters() pset.add_parameter("Wcrit"); pset.add_parameter("n"); pset.add_optional_parameter("eps", 1e-30); + pset.add_optional_parameter("log", false); return pset; } @@ -1321,7 +1323,7 @@ void WorkDamage::damage(double d_np1, double d_n, double dt = t_np1 - t_n; - double val = Wcrit_->value(wrate); + double val = Wcrit(wrate); *dd = d_n + n_ * std::pow(std::fabs(d_np1), (n_-1.0)/n_) * wrate * dt / val; @@ -1346,8 +1348,8 @@ void WorkDamage::ddamage_dd(double d_np1, double d_n, return; } - double val = Wcrit_->value(wrate); - double deriv = Wcrit_->derivative(wrate); + double val = Wcrit(wrate); + double deriv = dWcrit(wrate); double dt = t_np1 - t_n; double S[36]; @@ -1381,8 +1383,8 @@ void WorkDamage::ddamage_de(double d_np1, double d_n, return; } - double val = Wcrit_->value(wrate); - double dval = Wcrit_->derivative(wrate); + double val = Wcrit(wrate); + double dval = dWcrit(wrate); double fact = n_ * std::pow(std::fabs(d_np1), (n_-1.0)/n_) / val * (1.0 - wrate / val * dval) * (1.0 - std::fabs(d_np1)); @@ -1407,8 +1409,8 @@ void WorkDamage::ddamage_ds(double d_np1, double d_n, return; } - double val = Wcrit_->value(wrate); - double dval = Wcrit_->derivative(wrate); + double val = Wcrit(wrate); + double dval = dWcrit(wrate); double S[36]; elastic_->S(T_np1, S); @@ -1466,6 +1468,23 @@ double WorkDamage::workrate( return fabs(dot_vec(stress_np1, dp, 6) / dt * (1.0 - d_np1)); } +double WorkDamage::Wcrit(double Wdot) const +{ + if (log_) + return std::pow(10.0, Wcrit_->value(std::log10(Wdot))); + + return Wcrit_->value(Wdot); +} + +double WorkDamage::dWcrit(double Wdot) const +{ + if (log_) + return std::pow(10.0, Wcrit_->value(std::log10(Wdot))) * + Wcrit_->derivative(std::log10(Wdot)) / Wdot; + + return Wcrit_->derivative(Wdot); +} + PowerLawDamage::PowerLawDamage(ParameterSet & params) : StandardScalarDamage(params), A_(params.get_object_parameter("A")), diff --git a/test/test_damage.py b/test/test_damage.py index 9128c795..731e6725 100644 --- a/test/test_damage.py +++ b/test/test_damage.py @@ -275,6 +275,80 @@ def test_definition(self): self.assertAlmostEqual(damage, should) +class TestWorkDamageLog(unittest.TestCase, CommonScalarDamageModel, + CommonDamagedModel): + def setUp(self): + self.E = 92000.0 + self.nu = 0.3 + + self.s0 = 180.0 + self.Kp = 1000.0 + self.H = 1000.0 + + self.elastic = elasticity.IsotropicLinearElasticModel(self.E, "youngs", + self.nu, "poissons") + + surface = surfaces.IsoKinJ2() + iso = hardening.LinearIsotropicHardeningRule(self.s0, self.Kp) + kin = hardening.LinearKinematicHardeningRule(self.H) + hrule = hardening.CombinedHardeningRule(iso, kin) + + flow = ri_flow.RateIndependentAssociativeFlow(surface, hrule) + + self.bmodel = models.SmallStrainRateIndependentPlasticity(self.elastic, + flow) + + self.fn = interpolate.PolynomialInterpolate([0.1]) + self.n = 2.1 + + self.dmodel = damage.WorkDamage(self.elastic, self.fn, self.n, log = True) + + self.model = damage.NEMLScalarDamagedModel_sd(self.elastic, self.bmodel, + self.dmodel) + + self.stress = np.array([100,-50.0,300.0,-99,50.0,125.0]) * 0.75 + self.T = 100.0 + + self.s_np1 = self.stress + self.s_n = np.array([-25,150,250,-25,-100,25]) * 0.0 + + self.d_np1 = 0.2 + self.d_n = 0.1 + + self.e_np1 = np.array([0.1,-0.01,0.15,-0.05,-0.1,0.15]) * 0.75 + self.e_n = np.array([-0.05,0.025,-0.1,0.2,0.11,0.13]) * 0.0 + + self.T_np1 = self.T + self.T_n = 90.0 + + self.t_np1 = 5000.0 + self.t_n = 0.9 + + self.u_n = 0.0 + self.p_n = 0.0 + + # This is a rather boring baseline history state to probe, but I can't + # think of a better way to get a "generic" history from a generic model + self.hist_n = np.array([self.d_n] + list(self.bmodel.init_store())) + self.x_trial = np.array([50,-25,150,-150,190,100.0] + [0.41]) + + self.nsteps = 10 + self.etarget = np.array([0.1,-0.025,0.02,0.015,-0.02,-0.05]) + self.ttarget = 2.0 + + self.ee = np.dot(self.elastic.S(self.T_np1), self.s_np1*(1.0-self.d_np1) - self.s_n*(1.0-self.d_n)) + self.de = self.e_np1 - self.e_n + self.dp = self.de - self.ee + self.dt = self.t_np1 - self.t_n + self.Wdot = np.dot(self.s_np1*(1.0-self.d_np1), self.dp) / self.dt + + def test_definition(self): + damage = self.dmodel.damage(self.d_np1, self.d_n, self.e_np1, self.e_n, + self.s_np1, self.s_n, self.T_np1, self.T_n, self.t_np1, self.t_n) + should = self.d_n + self.n * self.d_np1**((self.n-1)/self.n) * self.Wdot * self.dt / np.exp(self.fn.value(np.log(self.Wdot))) + + self.assertAlmostEqual(damage, should) + class TestClassicalDamage(unittest.TestCase, CommonScalarDamageModel, CommonDamagedModel, CommonScalarDamageRate): def setUp(self):