Repository: Argonne-National-Laboratory/neml (language: python)

Issue:
Add prefactor to Chaboche's flow rule
A scaling prefactor is useful when converting material parameters from macroscopic to microscopic


Discussion / hints:

Fix patch:
diff --git a/doc/sphinx/vp_flow/chaboche.rst b/doc/sphinx/vp_flow/chaboche.rst
index 8f08c18a..d734f179 100644
--- a/doc/sphinx/vp_flow/chaboche.rst
+++ b/doc/sphinx/vp_flow/chaboche.rst
@@ -12,15 +12,16 @@ The model is defined by:
 
 .. math::
 
-   \dot{\gamma} = \sqrt{\frac{3}{2}} \left\langle \frac{f\left(\bm{\sigma}, \mathbf{q}\left(\bm{\alpha}\right), T\right)}{\sqrt{2/3}\eta\left(\bar{\varepsilon}_{vp}, T\right)}\right\rangle^n
+   \dot{\gamma} = \dot{\gamma}_0 \sqrt{\frac{3}{2}} \left\langle \frac{f\left(\bm{\sigma}, \mathbf{q}\left(\bm{\alpha}\right), T\right)}{\sqrt{2/3}\eta\left(\bar{\varepsilon}_{vp}, T\right)}\right\rangle^n
 
-   \mathbf{g}_{\gamma} = \frac{\partial f}{\partial \bm{\sigma}} 
+   \mathbf{g}_{\gamma} = \frac{\partial f}{\partial \bm{\sigma}}
       \left( \bm{\sigma}, \mathbf{q}\left(\bm{\alpha}\right), T  \right)
 
 The time and temperature rate contributions to the flow function are zero.
-The rate sensitivity exponent is generally temperature dependent; the 
+The rate sensitivity exponent is generally temperature dependent; the prefactor
+:math:`\dot_{\gamma}_0` can be temperature dependent; the
 :ref:`fluidity-model` :math:`\eta` can depend on both temperature and inelastic strain.
-The hardening model is defined by a NonassociativeHardening model. 
+The hardening model is defined by a NonassociativeHardening model.
 The time and temperature rate contributions can be non-zero.
 
 The model maintains the history variables defined by the hardening model.
@@ -30,12 +31,13 @@ Parameters
 
 .. csv-table::
    :header: "Parameter", "Object type", "Description", "Default"
-   :widths: 12, 30, 50, 8
+   :widths: 15, 30, 50, 8
 
    ``surface``, :cpp:class:`neml::YieldSurface`, Flow surface interface, No
    ``hardening``, :cpp:class:`neml::NonAssociativeHardeningRule`, Hardening rule interface, No
    ``fluidity``, :cpp:class:`neml::FluidityModel`, Fluidity definition, No
    ``n``, :cpp:class:`neml::Interpolate`, Rate sensitivity exponent, No
+   ``prefactor``, :cpp:class:`neml::Interpolate`, Prefactor, 1.
 
 Class description
 -----------------
diff --git a/src/visco_flow.cxx b/src/visco_flow.cxx
index bbd7a1c9..06b1940a 100644
--- a/src/visco_flow.cxx
+++ b/src/visco_flow.cxx
@@ -476,41 +476,39 @@ double SaturatingFluidity::deta(double a, double T) const
   return A * b * exp(-b * a);
 }
 
-ChabocheFlowRule::ChabocheFlowRule(std::shared_ptr<YieldSurface> surface,
-                                   std::shared_ptr<NonAssociativeHardening> hardening,
-                                   std::shared_ptr<FluidityModel> fluidity,
-                                   std::shared_ptr<Interpolate> n) :
-    surface_(surface), hardening_(hardening), fluidity_(fluidity), n_(n),
-    recovery_(false)
-{
-  
-}
+ChabocheFlowRule::ChabocheFlowRule(
+    std::shared_ptr<YieldSurface> surface,
+    std::shared_ptr<NonAssociativeHardening> hardening,
+    std::shared_ptr<FluidityModel> fluidity, std::shared_ptr<Interpolate> n,
+    std::shared_ptr<Interpolate> prefactor)
+    : surface_(surface), hardening_(hardening), fluidity_(fluidity), n_(n),
+      prefactor_(prefactor), recovery_(false) {}
 
 std::string ChabocheFlowRule::type()
 {
   return "ChabocheFlowRule";
 }
 
-ParameterSet ChabocheFlowRule::parameters()
-{
+ParameterSet ChabocheFlowRule::parameters() {
   ParameterSet pset(ChabocheFlowRule::type());
 
   pset.add_parameter<NEMLObject>("surface");
   pset.add_parameter<NEMLObject>("hardening");
   pset.add_parameter<NEMLObject>("fluidity");
   pset.add_parameter<NEMLObject>("n");
+  pset.add_optional_parameter<NEMLObject>(
+      "prefactor", std::make_shared<ConstantInterpolate>(1.0));
 
   return pset;
 }
 
-std::unique_ptr<NEMLObject> ChabocheFlowRule::initialize(ParameterSet & params)
-{
+std::unique_ptr<NEMLObject> ChabocheFlowRule::initialize(ParameterSet &params) {
   return neml::make_unique<ChabocheFlowRule>(
       params.get_object_parameter<YieldSurface>("surface"),
       params.get_object_parameter<NonAssociativeHardening>("hardening"),
       params.get_object_parameter<FluidityModel>("fluidity"),
-      params.get_object_parameter<Interpolate>("n")
-      ); 
+      params.get_object_parameter<Interpolate>("n"),
+      params.get_object_parameter<Interpolate>("prefactor"));
 }
 
 size_t ChabocheFlowRule::nhist() const
@@ -540,10 +538,9 @@ int ChabocheFlowRule::y(const double* const s, const double* const alpha, double
   if (ier != SUCCESS) return ier;
 
   if (fv > 0.0) {
-    double eta = sqrt(2.0/3.0) * fluidity_->eta(alpha[0], T);
-    yv = sqrt(3.0/2.0) * pow(fv/eta, n_->value(T));
-  }
-  else {
+    double eta = sqrt(2.0 / 3.0) * fluidity_->eta(alpha[0], T);
+    yv = sqrt(3.0 / 2.0) * pow(fv / eta, n_->value(T)) * prefactor_->value(T);
+  } else {
     yv = 0.0;
   }
 
@@ -566,10 +563,13 @@ int ChabocheFlowRule::dy_ds(const double* const s, const double* const alpha, do
 
   if (fv > 0.0) {
     ier = surface_->df_ds(s, q, T, dyv);
-    if (ier != SUCCESS) return ier;
-    double eta = sqrt(2.0/3.0) * fluidity_->eta(alpha[0], T);
-    double mv = sqrt(3.0/2.0) * pow(fv/eta, n_->value(T) - 1.0) * n_->value(T) / eta;
-    for (int i=0; i<6; i++) dyv[i] *= mv;
+    if (ier != SUCCESS)
+      return ier;
+    double eta = sqrt(2.0 / 3.0) * fluidity_->eta(alpha[0], T);
+    double mv = sqrt(3.0 / 2.0) * pow(fv / eta, n_->value(T) - 1.0) *
+                n_->value(T) / eta * prefactor_->value(T);
+    for (int i = 0; i < 6; i++)
+      dyv[i] *= mv;
   }
 
   return 0;
@@ -604,10 +604,14 @@ int ChabocheFlowRule::dy_da(const double* const s, const double* const alpha, do
     if (ier != SUCCESS) return ier;
 
     double eta = sqrt(2.0/3.0) * fluidity_->eta(alpha[0], T);
-    double mv = sqrt(3.0/2.0) * pow(fv/eta, n_->value(T) - 1.0) * n_->value(T) / eta;
-    for (size_t i=0; i<nhist(); i++) dyv[i] *= mv;
+    double mv = sqrt(3.0 / 2.0) * pow(fv / eta, n_->value(T) - 1.0) *
+                n_->value(T) / eta * prefactor_->value(T);
+    for (size_t i = 0; i < nhist(); i++)
+      dyv[i] *= mv;
+
+    double mv2 = -sqrt(3.0 / 2.0) * fv * pow(fv / eta, n_->value(T) - 1.0) *
+                 n_->value(T) / (eta * eta) * prefactor_->value(T);
 
-    double mv2 = -sqrt(3.0/2.0) * fv * pow(fv/eta, n_->value(T) - 1.0) * n_->value(T) / (eta * eta);
     double deta = sqrt(2.0/3.0) * fluidity_->deta(alpha[0], T);
     dyv[0] += deta * mv2;
   }
diff --git a/src/visco_flow.h b/src/visco_flow.h
index e7964dc7..fab12c47 100644
--- a/src/visco_flow.h
+++ b/src/visco_flow.h
@@ -270,7 +270,8 @@ class NEML_EXPORT ChabocheFlowRule: public ViscoPlasticFlowRule {
   ChabocheFlowRule(std::shared_ptr<YieldSurface> surface,
                    std::shared_ptr<NonAssociativeHardening> hardening,
                    std::shared_ptr<FluidityModel> fluidity,
-                   std::shared_ptr<Interpolate> n);
+                   std::shared_ptr<Interpolate> n,
+                   std::shared_ptr<Interpolate> prefactor);
 
   /// String type for the object system
   static std::string type();
@@ -339,6 +340,7 @@ class NEML_EXPORT ChabocheFlowRule: public ViscoPlasticFlowRule {
   std::shared_ptr<NonAssociativeHardening> hardening_;
   std::shared_ptr<FluidityModel> fluidity_;
   const std::shared_ptr<const Interpolate> n_;
+  const std::shared_ptr<const Interpolate> prefactor_;
   const bool recovery_;
 };
 


Test patch:
diff --git a/test/test_visco_flow.py b/test/test_visco_flow.py
index 3dfe0609..11789155 100644
--- a/test/test_visco_flow.py
+++ b/test/test_visco_flow.py
@@ -282,7 +282,7 @@ def gen_hist(self):
     h[13:19] = make_dev(h[13:19])
     return h
 
-class TestChebocheFlow(unittest.TestCase):
+class TestChabocheFlow(unittest.TestCase):
   def setUp(self):
     self.n = 20.0
     self.eta = 108.0
@@ -428,6 +428,199 @@ def test_Xdot(self):
     self.assertTrue(np.allclose(dX, dX_calc))
 
 
+class TestChabocheFlowWithPrefactor(unittest.TestCase):
+  def setUp(self):
+    self.n = 20.0
+    self.eta = 108.0
+    self.sY = 89.0
+    self.prefactor = np.asarray([2.])
+    self.prefactor = 2.
+
+    self.Q = 165.0
+    self.b = 12.0
+
+    self.m = 3
+
+    C1 = 80.0e3
+    C2 = 14.02e3
+    C3 = 3.333e3
+
+    y1 = 0.9e3
+    y2 = 1.5e3
+    y3 = 1.0
+
+    surface = surfaces.IsoKinJ2()
+    self.iso = hardening.VoceIsotropicHardeningRule(self.sY, self.Q, self.b)
+    self.cs = [C1, C2, C3]
+    self.gs = [y1, y2, y3]
+    self.As = [0.0, 0.0, 0.0]
+    self.ns = [1.0, 1.0, 1.0]
+    self.gmodels = [hardening.ConstantGamma(g) for g in self.gs]
+    self.hmodel = hardening.Chaboche(self.iso, self.cs, self.gmodels,
+        self.As, self.ns)
+
+    self.fluidity = visco_flow.ConstantFluidity(self.eta)
+
+    self.hist0 = np.zeros((19,))
+    self.T = 300.0
+
+    self.model = visco_flow.ChabocheFlowRule(
+        surface, self.hmodel, self.fluidity, self.n, prefactor = self.prefactor)
+
+  def gen_stress(self):
+    return np.array([300.0,-200.0,100.0,50.0,150.0,300.0])
+
+  def gen_hist(self):
+    h = np.array(range(1,self.m*6+2)) / (self.m*7)
+    h[1:self.m*6+1] = (2.0 * h[1:self.m*6+1] - 1.0) * 200
+    for i in range(self.m):
+      h[1+i*6:1+(i+1)*6] = make_dev(h[1+i*6:1+(i+1)*6])
+    return h
+
+  def test_y(self):
+    s = self.gen_stress()
+    h = self.gen_hist()
+    y = self.model.y(s, h, self.T)
+
+    sdev = make_dev(s)
+    X = sum(h[1+i*6:1+(i+1)*6] for i in range(self.m))
+    ih = self.iso.q(h[0:1], self.T)[0]
+
+    tv = np.sqrt(3.0/2.0) * la.norm(sdev + X) + ih
+
+    if tv > 0.0:
+      pd = self.prefactor*(tv / self.eta)**self.n
+    else:
+      pd = 0.0
+
+    gd = pd / np.sqrt(2.0/3.0)
+
+    self.assertTrue(np.isclose(gd, y))
+
+  def test_ep(self):
+    s = self.gen_stress()
+    h = self.gen_hist()
+    y = self.model.y(s, h, self.T)
+
+    ep_model = self.model.g(s, h, self.T) * y
+
+    sdev = make_dev(s)
+    X = sum(h[1+i*6:1+(i+1)*6] for i in range(self.m))
+    ih = self.iso.q(h[0:1], self.T)[0]
+
+    tv = np.sqrt(3.0/2.0) * la.norm(sdev + X) + ih
+
+    if tv > 0.0:
+      pd = self.prefactor*(tv / self.eta)**self.n
+    else:
+      pd = 0.0
+
+    ep_calc = 3.0/2.0 * pd * (sdev + X) / (np.sqrt(3.0/2.0) * la.norm(sdev+X))
+
+    self.assertTrue(np.allclose(ep_calc, ep_model))
+
+  def test_rdot(self):
+    s = self.gen_stress()
+    h = self.gen_hist()
+    y = self.model.y(s, h, self.T)
+    adot = self.model.h(s,h,self.T)[0] * y
+    dr = -self.b * np.exp(-self.b * h[0]) * self.Q * adot
+
+    sdev = make_dev(s)
+    X = sum(h[1+i*6:1+(i+1)*6] for i in range(self.m))
+    ih = self.iso.q(h[0:1], self.T)[0]
+
+    tv = np.sqrt(3.0/2.0) * la.norm(sdev + X) + ih
+
+    if tv > 0.0:
+      pd = self.prefactor*(tv / self.eta)**self.n
+    else:
+      pd = 0.0
+
+    q = self.hmodel.q(h, self.T)
+
+    Rtot = -self.sY - self.Q * (1.0 - np.exp(-self.b*h[0]))
+
+    self.assertTrue(np.isclose(Rtot, q[0]))
+
+    R_act = -(Rtot + self.sY)
+
+    dr_direct = -self.b*(self.Q - R_act) * pd
+
+    self.assertTrue(np.isclose(dr_direct, dr))
+
+  def test_Xdot(self):
+    s = self.gen_stress()
+    h = self.gen_hist()
+    y = self.model.y(s, h, self.T)
+    dX = self.model.h(s, h, self.T)[1:] * y
+
+    sdev = make_dev(s)
+    X = sum(h[1+i*6:1+(i+1)*6] for i in range(self.m))
+    ih = self.iso.q(h[0:1], self.T)[0]
+
+    tv = np.sqrt(3.0/2.0) * la.norm(sdev + X) + ih
+
+    if tv > 0.0:
+      pd = self.prefactor*(tv / self.eta)**self.n
+    else:
+      pd = 0.0
+
+    ep_calc = 3.0/2.0 * pd * (sdev + X) / (np.sqrt(3.0/2.0) * la.norm(sdev+X))
+
+    dX_calc = np.zeros((6*self.m,))
+
+    for i in range(self.m):
+      dX_calc[i*6:(i+1)*6] = -2.0/3.0 * self.cs[i] * ep_calc - self.gs[i] * h[1+i*6:1+(i+1)*6] * pd
+
+    self.assertTrue(np.allclose(dX, dX_calc))
+
+class TestChabocheModelWithFactor(unittest.TestCase, CommonFlowRule):
+  def setUp(self):
+    n = 20.0
+    eta = 108.0
+    sY = 89.0
+
+    Q = 165.0
+    b = 12.0
+    
+    self.m = 3
+
+    C1 = 80.0e3
+    C2 = 14.02e3
+    C3 = 3.333e3
+
+    y1 = 0.9e3
+    y2 = 1.5e3
+    y3 = 1.0
+
+    surface = surfaces.IsoKinJ2()
+    iso = hardening.VoceIsotropicHardeningRule(sY, Q, b)
+    cs = [C1, C2, C3]
+    gs = [y1, y2, y3]
+    gmodels = [hardening.ConstantGamma(g) for g in gs]
+    A = [0.0, 0.0, 0.0]
+    ae = [1.0, 1.0, 1.0]
+
+    self.prefactor = 2.0
+
+    hmodel = hardening.Chaboche(iso, cs, gmodels, A, ae)
+
+    fluidity = visco_flow.ConstantFluidity(eta)
+
+    self.hist0 = np.zeros((19,))
+    self.T = 300.0
+
+    self.model = visco_flow.ChabocheFlowRule(
+        surface, hmodel, fluidity, n, prefactor = self.prefactor)
+
+  def gen_hist(self):
+    h = np.array([0.05,20,-30,40.0,5.0,2.0,40.0,-10,-20,10,50,30,10,-50,60,70,15,-15,10])
+    h[1:7] = make_dev(h[1:7])
+    h[7:13] = make_dev(h[7:13])
+    h[13:19] = make_dev(h[13:19])
+    return h
+
 class CommonFluidity(object):
   """
     Common fluidity model tests