Added a simple linear viscous flow rule

Argonne-National-Laboratory · Jul 1, 2022 · 1418f4e · 1418f4e
1 parent b0193bf
commit 1418f4e
Show file tree

Hide file tree

Showing 6 changed files with 260 additions and 18 deletions.
diff --git a/doc/sphinx/vp_flow.rst b/doc/sphinx/vp_flow.rst
@@ -35,6 +35,7 @@ Implementations
 
    vp_flow/superimposed
    vp_flow/perzyna
+   vp_flow/linear
    vp_flow/chaboche
    vp_flow/yaguchi
 

diff --git a/doc/sphinx/vp_flow/linear.rst b/doc/sphinx/vp_flow/linear.rst
@@ -0,0 +1,34 @@
+Linear vicous flow rule
+=======================
+
+Overview
+--------
+
+This simple model provides a linear viscous response of the type
+
+.. math::
+
+   \dot{\gamma} = \frac{3}{2} \frac{f\left(\bm{\sigma}, \mathbf{0}, T\right)}{\eta\left(T\right)}
+
+   \mathbf{g}_{\gamma} = \frac{\partial f}{\partial \bm{\sigma}} \left( \bm{\sigma}, \mathbf{0}, T  \right)
+
+where :math:`f` is a flow surface
+
+The model does not maintain internal variables.
+
+Parameters
+----------
+
+.. csv-table::
+   :header: "Parameter", "Object type", "Description", "Default"
+   :widths: 12, 30, 50, 8
+
+   ``surface``, :cpp:class:`neml::YieldSurface`, Flow surface interface, No
+   ``eta``, :cpp:class:`neml::Interpolate`, Drag stress, No
+
+Class description
+-----------------
+
+.. doxygenclass:: neml::LinearViscousFlow
+   :members:
+   :undoc-members:
diff --git a/include/visco_flow.h b/include/visco_flow.h
@@ -299,6 +299,65 @@ class NEML_EXPORT PerzynaFlowRule : public ViscoPlasticFlowRule {
 
 static Register<PerzynaFlowRule> regPerzynaFlowRule;
 
+/// Linear viscous perfect plasticity
+class NEML_EXPORT LinearViscousFlow : public ViscoPlasticFlowRule {
+ public:
+  /// Parameters: just a surface and a drag stress
+  LinearViscousFlow(ParameterSet & params);
+
+  /// String type for the object system
+  static std::string type();
+  /// Default parameters
+  static std::unique_ptr<NEMLObject> initialize(ParameterSet & params);
+  /// Initialize from parameters
+  static ParameterSet parameters();
+
+  /// Number of history variables
+  virtual size_t nhist() const;
+  /// Initialize history at time zero
+  virtual void init_hist(double * const h) const;
+
+  /// Scalar strain rate
+  virtual void y(const double* const s, const double* const alpha, double T,
+                double & yv) const;
+  /// Derivative of y wrt stress
+  virtual void dy_ds(const double* const s, const double* const alpha, double T,
+                double * const dyv) const;
+  /// Derivative of y wrt history
+  virtual void dy_da(const double* const s, const double* const alpha, double T,
+                double * const dyv) const;
+
+  /// Flow rule proportional to the scalar strain rate
+  virtual void g(const double * const s, const double * const alpha, double T,
+                double * const gv) const;
+  /// Derivative of g wrt stress
+  virtual void dg_ds(const double * const s, const double * const alpha, double T,
+                double * const dgv) const;
+  /// Derivative of g wrt history
+  virtual void dg_da(const double * const s, const double * const alpha, double T,
+               double * const dgv) const;
+
+  /// Hardening rule proportional to the scalar strain rate
+  virtual void h(const double * const s, const double * const alpha, double T,
+                double * const hv) const;
+  /// Derivative of h wrt stress
+  virtual void dh_ds(const double * const s, const double * const alpha, double T,
+                double * const dhv) const;
+  /// Derivative of h wrt history
+  virtual void dh_da(const double * const s, const double * const alpha, double T,
+                double * const dhv) const;
+
+ protected:
+  std::vector<double> fake_hist_() const;
+
+ private:
+  std::shared_ptr<YieldSurface> surface_;
+  std::shared_ptr<Interpolate> eta_;
+};
+
+static Register<LinearViscousFlow> regLinearViscousFLow;
+
+
 /// Various Chaboche type fluidity models.
 //
 //  These depend only on the equivalent plastic strain

diff --git a/src/visco_flow.cxx b/src/visco_flow.cxx
@@ -227,13 +227,15 @@ void SuperimposedViscoPlasticFlowRule::dg_ds(const double * const s,
     for (size_t j = 0; j < 36; j++) dgv[j] += yi * dgvi[j];
     outer_update(gi, 6, dyvi, 6, dgv);
   }
-  for (size_t i = 0; i < 36; i++) dgv[i] /= yv;
+  if (yv > 0)
+    for (size_t i = 0; i < 36; i++) dgv[i] /= yv;
 
   double dy[6];
   dy_ds(s, alpha, T, dy);
   double gv[6];
   g(s, alpha, T, gv);
-  for (size_t i = 0; i < 6; i++) gv[i] /= yv;
+  if (yv > 0)
+    for (size_t i = 0; i < 6; i++) gv[i] /= yv;
 
   outer_update_minus(gv, 6, dy, 6, dgv);
 }
@@ -270,15 +272,17 @@ void SuperimposedViscoPlasticFlowRule::dg_da(const double * const s,
     delete [] dgvi;
     delete [] dyi;
   }
-
-  for (size_t i = 0; i < 6*nhist(); i++) dgv[i] /= yv;
+
+  if (yv > 0)
+    for (size_t i = 0; i < 6*nhist(); i++) dgv[i] /= yv;
 
   double * dy = new double [nhist()];
   double gv[6];
   dy_da(s, alpha, T, dy);
   g(s, alpha, T, gv);
-
-  for (size_t i = 0; i < 6; i++) gv[i] /= yv;
+
+  if (yv > 0)
+    for (size_t i = 0; i < 6; i++) gv[i] /= yv;
 
   outer_update_minus(gv, 6, dy, nhist(), dgv);
 
@@ -327,13 +331,16 @@ void SuperimposedViscoPlasticFlowRule::dg_ds_time(const double * const s,
     for (size_t j = 0; j < 36; j++) dgv[j] += yi * dgvi[j];
     outer_update(gi, 6, dyvi, 6, dgv);
   }
-  for (size_t i = 0; i < 36; i++) dgv[i] /= yv;
+  if (yv > 0)
+    for (size_t i = 0; i < 36; i++) dgv[i] /= yv;
 
   double dy[6];
   dy_ds(s, alpha, T, dy);
   double gv[6];
   g_time(s, alpha, T, gv);
-  for (size_t i = 0; i < 6; i++) gv[i] /= yv;
+
+  if (yv > 0)
+    for (size_t i = 0; i < 6; i++) gv[i] /= yv;
 
   outer_update_minus(gv, 6, dy, 6, dgv);
 }
@@ -370,15 +377,17 @@ void SuperimposedViscoPlasticFlowRule::dg_da_time(const double * const s,
     delete [] dgvi;
     delete [] dyi;
   }
-
-  for (size_t i = 0; i < 6*nhist(); i++) dgv[i] /= yv;
+
+  if (yv > 0)
+    for (size_t i = 0; i < 6*nhist(); i++) dgv[i] /= yv;
 
   double * dy = new double [nhist()];
   double gv[6];
   dy_da(s, alpha, T, dy);
   g_time(s, alpha, T, gv);
-
-  for (size_t i = 0; i < 6; i++) gv[i] /= yv;
+
+  if (yv > 0)
+    for (size_t i = 0; i < 6; i++) gv[i] /= yv;
 
   outer_update_minus(gv, 6, dy, nhist(), dgv);
 
@@ -427,13 +436,15 @@ void SuperimposedViscoPlasticFlowRule::dg_ds_temp(const double * const s,
     for (size_t j = 0; j < 36; j++) dgv[j] += yi * dgvi[j];
     outer_update(gi, 6, dyvi, 6, dgv);
   }
-  for (size_t i = 0; i < 36; i++) dgv[i] /= yv;
+  if (yv > 0)
+    for (size_t i = 0; i < 36; i++) dgv[i] /= yv;
 
   double dy[6];
   dy_ds(s, alpha, T, dy);
   double gv[6];
   g_temp(s, alpha, T, gv);
-  for (size_t i = 0; i < 6; i++) gv[i] /= yv;
+  if (yv > 0)
+    for (size_t i = 0; i < 6; i++) gv[i] /= yv;
 
   outer_update_minus(gv, 6, dy, 6, dgv);
 }
@@ -470,15 +481,17 @@ void SuperimposedViscoPlasticFlowRule::dg_da_temp(const double * const s,
     delete [] dgvi;
     delete [] dyi;
   }
-
-  for (size_t i = 0; i < 6*nhist(); i++) dgv[i] /= yv;
+
+  if (yv > 0)
+    for (size_t i = 0; i < 6*nhist(); i++) dgv[i] /= yv;
 
   double * dy = new double [nhist()];
   double gv[6];
   dy_da(s, alpha, T, dy);
   g_temp(s, alpha, T, gv);
-
-  for (size_t i = 0; i < 6; i++) gv[i] /= yv;
+
+  if (yv > 0)
+    for (size_t i = 0; i < 6; i++) gv[i] /= yv;
 
   outer_update_minus(gv, 6, dy, nhist(), dgv);
 
@@ -853,6 +866,120 @@ void PerzynaFlowRule::dh_da(const double * const s, const double * const alpha,
   mat_mat(nhist(), nhist(), nhist(), dd, jac, dhv);
 }
 
+LinearViscousFlow::LinearViscousFlow(ParameterSet & params) :
+    ViscoPlasticFlowRule(params),
+    surface_(params.get_object_parameter<YieldSurface>("surface")),
+    eta_(params.get_object_parameter<Interpolate>("eta"))
+{
+
+}
+
+std::string LinearViscousFlow::type()
+{
+  return "LinearViscousFlow";
+}
+
+ParameterSet LinearViscousFlow::parameters()
+{
+  ParameterSet pset(LinearViscousFlow::type());
+
+  pset.add_parameter<NEMLObject>("surface");
+  pset.add_parameter<NEMLObject>("eta");
+
+  return pset;
+}
+
+std::unique_ptr<NEMLObject> LinearViscousFlow::initialize(ParameterSet & params)
+{
+  return neml::make_unique<LinearViscousFlow>(params); 
+}
+
+size_t LinearViscousFlow::nhist() const
+{
+  return 0;
+}
+
+void LinearViscousFlow::init_hist(double * const h) const
+{
+  return;
+}
+
+// Zero history of the correct size
+std::vector<double> LinearViscousFlow::fake_hist_() const
+{
+  std::vector<double> h(surface_->nhist());
+  std::fill(h.begin(), h.end(), 0.0);
+  return h;
+}
+
+// Rate rule
+void LinearViscousFlow::y(const double* const s, const double* const alpha, double T,
+              double & yv) const
+{
+  auto h = fake_hist_();
+
+  double fv;
+  surface_->f(s, &h[0], T, fv);
+
+  yv = 3.0 / 2.0 * fv / eta_->value(T);
+}
+
+void LinearViscousFlow::dy_ds(const double* const s, const double* const alpha, double T,
+              double * const dyv) const
+{
+  auto h = fake_hist_();
+
+  surface_->df_ds(s, &h[0], T, dyv);
+  for (size_t i = 0; i < 6; i++)
+    dyv[i] = 3.0 / 2.0 * dyv[i] / eta_->value(T);
+}
+
+void LinearViscousFlow::dy_da(const double* const s, const double* const alpha, double T,
+              double * const dyv) const
+{
+  return;
+}
+
+// Flow rule
+void LinearViscousFlow::g(const double * const s, const double * const alpha, double T,
+              double * const gv) const
+{
+  auto h = fake_hist_();
+  surface_->df_ds(s, &h[0], T, gv);
+}
+
+void LinearViscousFlow::dg_ds(const double * const s, const double * const alpha, double T,
+              double * const dgv) const
+{
+  auto h = fake_hist_();
+  surface_->df_dsds(s, &h[0], T, dgv);
+}
+
+void LinearViscousFlow::dg_da(const double * const s, const double * const alpha, double T,
+             double * const dgv) const
+{
+  return;
+}
+
+// Hardening rule
+void LinearViscousFlow::h(const double * const s, const double * const alpha, double T,
+              double * const hv) const
+{
+  return;
+}
+
+void LinearViscousFlow::dh_ds(const double * const s, const double * const alpha, double T,
+              double * const dhv) const
+{
+  return;
+}
+
+void LinearViscousFlow::dh_da(const double * const s, const double * const alpha, double T,
+              double * const dhv) const
+{
+  return;
+}
+
 FluidityModel::FluidityModel(ParameterSet & params) :
     NEMLObject(params)
 {

diff --git a/src/visco_flow_wrap.cxx b/src/visco_flow_wrap.cxx
@@ -187,6 +187,14 @@ PYBIND11_MODULE(visco_flow, m) {
       }))
       ;
 
+  py::class_<LinearViscousFlow, ViscoPlasticFlowRule, std::shared_ptr<LinearViscousFlow>>(m, "LinearViscousFlow")
+    .def(py::init([](py::args args, py::kwargs kwargs)
+      {
+        return create_object_python<LinearViscousFlow>(args, kwargs, {"surface",
+                                                     "eta"});
+      }))
+      ;
+
   py::class_<GFlow, NEMLObject, std::shared_ptr<GFlow>>(m, "GFlow")
       .def("g", &GFlow::g, "g function in Perzyna model")
       .def("dg", &GFlow::dg, "Derivative of g wrt f")

diff --git a/test/test_visco_flow.py b/test/test_visco_flow.py
@@ -267,6 +267,19 @@ def test_superposition(self):
   def gen_hist(self):
     return np.array([0.01, 0.02])
 
+class TestLinearViscous(unittest.TestCase, CommonFlowRule):
+  def setUp(self):
+    self.eta = 20.0
+
+    self.surface = surfaces.IsoJ2()
+    self.model = visco_flow.LinearViscousFlow(self.surface, self.eta)
+
+    self.hist0 = np.zeros((0,))
+    self.T = 300.0
+
+  def gen_hist(self):
+    return np.zeros((0,))
+
 class TestPerzynaIsoJ2Voce(unittest.TestCase, CommonFlowRule):
   def setUp(self):
     self.s0 = 180.0