From 1418f4e6371b02c33933426a518b523390b2ca9b Mon Sep 17 00:00:00 2001 From: "Mark C. Messner" Date: Fri, 1 Jul 2022 14:58:35 -0500 Subject: [PATCH] Added a simple linear viscous flow rule --- doc/sphinx/vp_flow.rst | 1 + doc/sphinx/vp_flow/linear.rst | 34 +++++++ include/visco_flow.h | 59 ++++++++++++ src/visco_flow.cxx | 163 ++++++++++++++++++++++++++++++---- src/visco_flow_wrap.cxx | 8 ++ test/test_visco_flow.py | 13 +++ 6 files changed, 260 insertions(+), 18 deletions(-) create mode 100644 doc/sphinx/vp_flow/linear.rst diff --git a/doc/sphinx/vp_flow.rst b/doc/sphinx/vp_flow.rst index bf534c57..7c3bd732 100644 --- 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 new file mode 100644 index 00000000..7e17f273 --- /dev/null +++ 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 index 45316614..0c5683eb 100644 --- a/include/visco_flow.h +++ b/include/visco_flow.h @@ -299,6 +299,65 @@ class NEML_EXPORT PerzynaFlowRule : public ViscoPlasticFlowRule { static Register 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 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 fake_hist_() const; + + private: + std::shared_ptr surface_; + std::shared_ptr eta_; +}; + +static Register 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 index e34f1b5e..2933151d 100644 --- 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("surface")), + eta_(params.get_object_parameter("eta")) +{ + +} + +std::string LinearViscousFlow::type() +{ + return "LinearViscousFlow"; +} + +ParameterSet LinearViscousFlow::parameters() +{ + ParameterSet pset(LinearViscousFlow::type()); + + pset.add_parameter("surface"); + pset.add_parameter("eta"); + + return pset; +} + +std::unique_ptr LinearViscousFlow::initialize(ParameterSet & params) +{ + return neml::make_unique(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 LinearViscousFlow::fake_hist_() const +{ + std::vector 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 index 2cf28396..2f172e47 100644 --- a/src/visco_flow_wrap.cxx +++ b/src/visco_flow_wrap.cxx @@ -187,6 +187,14 @@ PYBIND11_MODULE(visco_flow, m) { })) ; + py::class_>(m, "LinearViscousFlow") + .def(py::init([](py::args args, py::kwargs kwargs) + { + return create_object_python(args, kwargs, {"surface", + "eta"}); + })) + ; + py::class_>(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 index deeb5b87..e6cb3416 100644 --- 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