Line data Source code
1 : // Copyright 2024, UChicago Argonne, LLC
2 : // All Rights Reserved
3 : // Software Name: NEML2 -- the New Engineering material Model Library, version 2
4 : // By: Argonne National Laboratory
5 : // OPEN SOURCE LICENSE (MIT)
6 : //
7 : // Permission is hereby granted, free of charge, to any person obtaining a copy
8 : // of this software and associated documentation files (the "Software"), to deal
9 : // in the Software without restriction, including without limitation the rights
10 : // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
11 : // copies of the Software, and to permit persons to whom the Software is
12 : // furnished to do so, subject to the following conditions:
13 : //
14 : // The above copyright notice and this permission notice shall be included in
15 : // all copies or substantial portions of the Software.
16 : //
17 : // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
18 : // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
19 : // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
20 : // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
21 : // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
22 : // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
23 : // THE SOFTWARE.
24 :
25 : #include "neml2/models/solid_mechanics/KocksMeckingFlowViscosity.h"
26 : #include "neml2/tensors/Scalar.h"
27 : #include "neml2/tensors/functions/pow.h"
28 : #include "neml2/tensors/functions/exp.h"
29 :
30 : namespace neml2
31 : {
32 : register_NEML2_object(KocksMeckingFlowViscosity);
33 :
34 : OptionSet
35 2 : KocksMeckingFlowViscosity::expected_options()
36 : {
37 2 : OptionSet options = Model::expected_options();
38 2 : options.doc() =
39 : "Calculates the temperature-dependent flow viscosity for a Perzyna-type model using the "
40 : "Kocks-Mecking model. The value is \\f$ \\eta = \\exp{B} \\mu "
41 : "\\dot{\\varepsilon}_0^\\frac{-k T A}{\\mu b^3} \\f$ with \\f$ \\mu "
42 : "\\f$ the shear modulus, \\f$ \\dot{\\varepsilon}_0 \\f$ a reference strain rate, \\f$ b "
43 : "\\f$ the Burgers vector, "
44 : "\\f$ k\\f$ the Boltzmann constant, "
45 : "\\f$ T \\f$ absolute temperature, \\f$ A \\f$ the Kocks-Mecking slope parameter, and \\f$ B "
46 2 : "\\f$ the Kocks-Mecking intercept parameter.";
47 :
48 4 : options.set<bool>("define_second_derivatives") = true;
49 :
50 4 : options.set_parameter<TensorName<Scalar>>("A");
51 4 : options.set("A").doc() = "The Kocks-Mecking slope parameter";
52 4 : options.set_parameter<TensorName<Scalar>>("B");
53 4 : options.set("B").doc() = "The Kocks-Mecking intercept parameter";
54 4 : options.set_parameter<TensorName<Scalar>>("shear_modulus");
55 4 : options.set("shear_modulus").doc() = "The shear modulus";
56 :
57 4 : options.set<double>("eps0");
58 4 : options.set("eps0").doc() = "The reference strain rate";
59 :
60 4 : options.set<double>("k");
61 4 : options.set("k").doc() = "Boltzmann constant";
62 4 : options.set<double>("b");
63 2 : options.set("b").doc() = "The Burgers vector";
64 :
65 6 : options.set_input("temperature") = VariableName(FORCES, "T");
66 2 : options.set("temperature").doc() = "Absolute temperature";
67 :
68 2 : return options;
69 0 : }
70 :
71 2 : KocksMeckingFlowViscosity::KocksMeckingFlowViscosity(const OptionSet & options)
72 : : Model(options),
73 6 : _A(declare_parameter<Scalar>("A", "A", /*allow_nonlinear=*/true)),
74 8 : _B(declare_parameter<Scalar>("B", "B", /*allow_nonlinear=*/true)),
75 10 : _mu(declare_parameter<Scalar>("mu", "shear_modulus", /*allow_nonlinear=*/true)),
76 4 : _eps0(options.get<double>("eps0")),
77 4 : _k(options.get<double>("k")),
78 10 : _b3(options.get<double>("b") * options.get<double>("b") * options.get<double>("b")),
79 2 : _T(declare_input_variable<Scalar>("temperature")),
80 4 : _eta(declare_output_variable<Scalar>(VariableName(PARAMETERS, name())))
81 : {
82 2 : }
83 :
84 : void
85 6 : KocksMeckingFlowViscosity::set_value(bool out, bool dout_din, bool d2out_din2)
86 : {
87 6 : auto post = pow(_eps0, _k * _T * _A / (_mu * _b3));
88 :
89 6 : if (out)
90 4 : _eta = exp(_B) * _mu * post;
91 :
92 6 : if (dout_din)
93 : {
94 3 : if (_T.is_dependent())
95 3 : _eta.d(_T) = _A * exp(_B) * _k * std::log(_eps0) * post / _b3;
96 :
97 9 : if (const auto * const A = nl_param("A"))
98 2 : _eta.d(*A) = exp(_B) * _k * _T * std::log(_eps0) * post / _b3;
99 :
100 9 : if (const auto * const B = nl_param("B"))
101 2 : _eta.d(*B) = exp(_B) * _mu * post;
102 :
103 9 : if (const auto * const mu = nl_param("mu"))
104 2 : _eta.d(*mu) = exp(_B) * post * (1.0 - _A * _k * _T * std::log(_eps0) / (_b3 * _mu));
105 : }
106 :
107 6 : if (d2out_din2)
108 : {
109 : // T
110 2 : if (_T.is_dependent())
111 : {
112 2 : _eta.d(_T, _T) = pow(_A * _k * std::log(_eps0) / _b3, 2.0) * exp(_B) * post / _mu;
113 :
114 6 : if (const auto * const A = nl_param("A"))
115 4 : _eta.d(_T, *A) = exp(_B) * _k * std::log(_eps0) * post / _b3 *
116 5 : (std::log(_eps0) * _A * _k * _T / (_b3 * _mu) + 1.0);
117 :
118 6 : if (const auto * const B = nl_param("B"))
119 1 : _eta.d(_T, *B) = _A * exp(_B) * _k * std::log(_eps0) * post / _b3;
120 :
121 6 : if (const auto * const mu = nl_param("mu"))
122 1 : _eta.d(_T, *mu) = -pow(_A * _k * std::log(_eps0) / (_b3 * _mu), 2.0) * exp(_B) * _T * post;
123 : }
124 :
125 : // A
126 6 : if (const auto * const A = nl_param("A"))
127 : {
128 1 : if (_T.is_dependent())
129 4 : _eta.d(*A, _T) = exp(_B) * _k * std::log(_eps0) * post / _b3 *
130 5 : (std::log(_eps0) * _A * _k * _T / (_b3 * _mu) + 1.0);
131 :
132 1 : _eta.d(*A, *A) = exp(_B) * pow(_k * _T * std::log(_eps0) / _b3, 2.0) * post / _mu;
133 :
134 3 : if (const auto * const B = nl_param("B"))
135 1 : _eta.d(*A, *B) = exp(_B) * _k * _T * std::log(_eps0) * post / _b3;
136 :
137 3 : if (const auto * const mu = nl_param("mu"))
138 1 : _eta.d(*A, *mu) = -_A * exp(_B) * pow(_k * _T * std::log(_eps0) / (_b3 * _mu), 2.0) * post;
139 : }
140 :
141 : // B
142 6 : if (const auto * const B = nl_param("B"))
143 : {
144 1 : if (_T.is_dependent())
145 1 : _eta.d(*B, _T) = _A * exp(_B) * _k * std::log(_eps0) * post / _b3;
146 :
147 3 : if (const auto * const A = nl_param("A"))
148 1 : _eta.d(*B, *A) = exp(_B) * _k * _T * std::log(_eps0) * post / _b3;
149 :
150 1 : _eta.d(*B, *B) = exp(_B) * _mu * post;
151 :
152 3 : if (const auto * const mu = nl_param("mu"))
153 2 : _eta.d(*B, *mu) =
154 3 : exp(_B) * post * (_b3 * _mu - _A * _k * _T * std::log(_eps0)) / (_b3 * _mu);
155 : }
156 :
157 : // mu
158 6 : if (const auto * const mu = nl_param("mu"))
159 : {
160 1 : if (_T.is_dependent())
161 1 : _eta.d(*mu, _T) = -exp(_B) * pow(_A * _k * std::log(_eps0) / (_b3 * _mu), 2.0) * _T * post;
162 :
163 3 : if (const auto * const A = nl_param("A"))
164 1 : _eta.d(*mu, *A) = -_A * exp(_B) * pow(_k * _T * std::log(_eps0) / (_b3 * _mu), 2.0) * post;
165 :
166 3 : if (const auto * const B = nl_param("B"))
167 2 : _eta.d(*mu, *B) =
168 3 : exp(_B) * post * (_b3 * _mu - _A * _k * _T * std::log(_eps0)) / (_b3 * _mu);
169 :
170 2 : _eta.d(*mu, *mu) =
171 3 : -pow(_A * _k * _T * std::log(_eps0) / (_b3 * _mu), 2.0) * exp(_B) * post / _mu;
172 : }
173 : }
174 6 : }
175 : } // namespace neml2
|