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/GTNYieldFunction.h"
26 : #include "neml2/tensors/Scalar.h"
27 : #include "neml2/tensors/functions/pow.h"
28 : #include "neml2/tensors/functions/cosh.h"
29 : #include "neml2/tensors/functions/sinh.h"
30 :
31 : namespace neml2
32 : {
33 : register_NEML2_object(GTNYieldFunction);
34 :
35 : OptionSet
36 2 : GTNYieldFunction::expected_options()
37 : {
38 2 : OptionSet options = Model::expected_options();
39 2 : options.doc() =
40 : "Gurson-Tvergaard-Needleman yield function for poroplasticity. The yield function is defined "
41 : "as \\f$ f = \\left( \\frac{\\bar{\\sigma}}{\\sigma_y + k} \\right)^2 + 2 q_1 \\phi \\cosh "
42 : "\\left( \\frac{1}{2} q_2 \\frac{3\\sigma_h-\\sigma_s}{\\sigma_y + k} \\right) - \\left( q_3 "
43 : "\\phi^2 + 1 \\right) \\f$, where \\f$ \\bar{\\sigma} \\f$ is the von Mises stress, \\f$ "
44 : "\\sigma_y \\f$ is the yield stress, \\f$ k \\f$ is isotropic hardening, \\f$ \\phi \\f$ is "
45 : "the porosity, \\f$ \\sigma_h \\f$ is the hydrostatic stress, and \\f$ \\sigma_s \\f$ is the "
46 : "void growth back stress (sintering stress). \\f$ q_1 \\f$, \\f$ q_2 \\f$, and \\f$ q_3 \\f$ "
47 2 : "are parameters controlling the yield mechanisms.";
48 :
49 4 : options.set<bool>("define_second_derivatives") = true;
50 :
51 4 : options.set_parameter<TensorName<Scalar>>("yield_stress");
52 4 : options.set("yield_stress").doc() = "Yield stress";
53 :
54 4 : options.set_parameter<TensorName<Scalar>>("q1");
55 6 : options.set("q1").doc() =
56 2 : "Parameter controlling the balance/competition between plastic flow and void evolution.";
57 :
58 4 : options.set_parameter<TensorName<Scalar>>("q2");
59 4 : options.set("q2").doc() = "Void evolution rate";
60 :
61 4 : options.set_parameter<TensorName<Scalar>>("q3");
62 2 : options.set("q3").doc() = "Pore pressure";
63 :
64 6 : options.set_input("flow_invariant") = VariableName(STATE, "internal", "se");
65 2 : options.set("flow_invariant").doc() = "Effective stress driving plastic flow";
66 :
67 6 : options.set_input("poro_invariant") = VariableName(STATE, "internal", "sp");
68 4 : options.set("poro_invariant").doc() = "Effective stress driving porous flow";
69 :
70 4 : options.set_input("isotropic_hardening");
71 2 : options.set("isotropic_hardening").doc() = "Isotropic hardening";
72 :
73 6 : options.set_input("void_fraction") = VariableName(STATE, "internal", "f");
74 2 : options.set("void_fraction").doc() = "Void fraction (porosity)";
75 :
76 6 : options.set_output("yield_function") = VariableName(STATE, "internal", "fp");
77 2 : options.set("yield_function").doc() = "Yield function";
78 :
79 2 : return options;
80 0 : }
81 :
82 3 : GTNYieldFunction::GTNYieldFunction(const OptionSet & options)
83 : : Model(options),
84 3 : _f(declare_output_variable<Scalar>("yield_function")),
85 3 : _se(declare_input_variable<Scalar>("flow_invariant")),
86 3 : _sp(declare_input_variable<Scalar>("poro_invariant")),
87 3 : _phi(declare_input_variable<Scalar>("void_fraction")),
88 6 : _h(options.get<VariableName>("isotropic_hardening").empty()
89 3 : ? nullptr
90 3 : : &declare_input_variable<Scalar>("isotropic_hardening")),
91 12 : _s0(declare_parameter<Scalar>("sy", "yield_stress", /*allow_nonlinear=*/true)),
92 12 : _q1(declare_parameter<Scalar>("q1", "q1", /*allow_nonlinear=*/true)),
93 12 : _q2(declare_parameter<Scalar>("q2", "q2", /*allow_nonlinear=*/true)),
94 15 : _q3(declare_parameter<Scalar>("q3", "q3", /*allow_nonlinear=*/true))
95 : {
96 3 : }
97 :
98 : void
99 9 : GTNYieldFunction::set_value(bool out, bool dout_din, bool d2out_din2)
100 : {
101 : // Flow stress (depending on whether isotropic hardening is provided)
102 9 : const auto sf = _h ? _s0 + (*_h) : _s0;
103 :
104 9 : if (out)
105 14 : _f = pow(_se / sf, 2.0) + 2 * _q1 * _phi * cosh(_q2 / 2.0 * _sp / sf) -
106 21 : (1.0 + _q3 * pow(Scalar(_phi), 2.0));
107 :
108 9 : if (dout_din)
109 : {
110 5 : if (_se.is_dependent())
111 5 : _f.d(_se) = 2.0 * _se / pow(sf, 2.0);
112 :
113 5 : if (_sp.is_dependent())
114 5 : _f.d(_sp) = _q1 * _phi * _q2 / sf * sinh(_q2 / 2.0 * _sp / sf);
115 :
116 5 : if (_phi.is_dependent())
117 5 : _f.d(_phi) = 2.0 * _q1 * cosh(_q2 / 2.0 * _sp / sf) - 2.0 * _q3 * _phi;
118 :
119 5 : if (_h)
120 20 : _f.d(*_h) = -2 * pow(Scalar(_se), 2.0) / pow(sf, 3.0) -
121 25 : _q1 * _phi * _q2 * _sp / pow(sf, 2.0) * sinh(_q2 / 2.0 * _sp / sf);
122 :
123 : // Handle the case of variable coupling
124 15 : if (const auto * const sy = nl_param("sy"))
125 8 : _f.d(*sy) = -2 * pow(Scalar(_se), 2.0) / pow(sf, 3.0) -
126 10 : _q1 * _phi * _q2 * _sp / pow(sf, 2.0) * sinh(_q2 / 2.0 * _sp / sf);
127 :
128 15 : if (const auto * const q1 = nl_param("q1"))
129 2 : _f.d(*q1) = 2.0 * _phi * cosh(_q2 / 2.0 * _sp / sf);
130 :
131 15 : if (const auto * const q2 = nl_param("q2"))
132 2 : _f.d(*q2) = _q1 * _phi * _sp / sf * sinh(_q2 / 2.0 * _sp / sf);
133 :
134 15 : if (const auto * const q3 = nl_param("q3"))
135 2 : _f.d(*q3) = -pow(Scalar(_phi), 2.0);
136 : }
137 :
138 9 : if (d2out_din2)
139 : {
140 6 : const auto * const sy = nl_param("sy");
141 6 : const auto * const q1 = nl_param("q1");
142 6 : const auto * const q2 = nl_param("q2");
143 6 : const auto * const q3 = nl_param("q3");
144 :
145 : ////////////////////////////////////////////////////////////////////////////////////////////////
146 : //
147 : // The GTN yield function can be expressed as
148 : //
149 : // f(se, sp, phi, h; sy, q1, q2, q3)
150 : //
151 : // - Arguments before the semicolon are variables
152 : // - Arguments after the semicolon are (nonlinear) parameters
153 : // - Derivatives w.r.t. the first three arguments se, sp, and phi are mandatory
154 : // - Derivatives w.r.t. the rest of the arguments are optional
155 : //
156 : ////////////////////////////////////////////////////////////////////////////////////////////////
157 :
158 : ////////////////////////////////////////////////////////////////////////////////////////////////
159 : //
160 : // The second derivative is nothing but a big matrix. We will fill out the matrix row by row, in
161 : // the order of the arguments listed above.
162 : //
163 : // Rows will separated by big fences like this.
164 : //
165 : ////////////////////////////////////////////////////////////////////////////////////////////////
166 :
167 : ////////////////////////////////////////////////////////////////////////////////////////////////
168 : //
169 : // f(se, sp, phi, h; sy, q1, q2, q3)
170 : //
171 : // se: Flow invariant
172 : //
173 : ////////////////////////////////////////////////////////////////////////////////////////////////
174 3 : if (_se.is_dependent())
175 : {
176 3 : _f.d(_se, _se) = 2.0 / pow(sf, 2.0);
177 :
178 3 : if (_h)
179 3 : _f.d(_se, *_h) = -4.0 * _se / pow(sf, 3.0);
180 :
181 3 : if (sy)
182 1 : _f.d(_se, *sy) = -4.0 * _se / pow(sf, 3.0);
183 : }
184 :
185 : ////////////////////////////////////////////////////////////////////////////////////////////////
186 : //
187 : // f(se, sp, phi, h; sy, q1, q2, q3)
188 : //
189 : // sp: Poro invariant
190 : //
191 : ////////////////////////////////////////////////////////////////////////////////////////////////
192 3 : if (_sp.is_dependent())
193 : {
194 6 : _f.d(_sp, _sp) =
195 9 : _phi * _q1 * pow(_q2, 2.0) / (2.0 * pow(sf, 2.0)) * cosh(_q2 / 2.0 * _sp / sf);
196 :
197 3 : if (_phi.is_dependent())
198 3 : _f.d(_sp, _phi) = _q1 * _q2 * sinh(_q2 / 2.0 * _sp / sf) / sf;
199 :
200 3 : if (_h)
201 6 : _f.d(_sp, *_h) =
202 6 : -_phi * _q1 * _q2 *
203 12 : (_q2 * _sp * cosh(_q2 / 2.0 * _sp / sf) + 2 * sf * sinh(_q2 / 2.0 * _sp / sf)) /
204 15 : (2 * pow(sf, 3.0));
205 3 : if (sy)
206 2 : _f.d(_sp, *sy) =
207 2 : -_phi * _q1 * _q2 *
208 4 : (_q2 * _sp * cosh(_q2 / 2.0 * _sp / sf) + 2 * sf * sinh(_q2 / 2.0 * _sp / sf)) /
209 5 : (2 * pow(sf, 3.0));
210 :
211 3 : if (q1)
212 1 : _f.d(_sp, *q1) = _phi * _q2 * sinh(_q2 / 2.0 * _sp / sf) / sf;
213 :
214 3 : if (q2)
215 2 : _f.d(_sp, *q2) =
216 2 : _phi * _q1 *
217 4 : (_q2 * _sp * cosh(_q2 / 2.0 * _sp / sf) + 2.0 * sf * sinh(_q2 / 2.0 * _sp / sf)) /
218 5 : (2.0 * pow(sf, 2.0));
219 : }
220 :
221 : ////////////////////////////////////////////////////////////////////////////////////////////////
222 : //
223 : // f(se, sp, phi, h; sy, q1, q2, q3)
224 : //
225 : // phi: Void fraction
226 : //
227 : ////////////////////////////////////////////////////////////////////////////////////////////////
228 3 : if (_phi.is_dependent())
229 : {
230 3 : if (_sp.is_dependent())
231 3 : _f.d(_phi, _sp) = _q1 * _q2 * sinh(_q2 / 2.0 * _sp / sf) / sf;
232 :
233 3 : _f.d(_phi, _phi) = -2.0 * _q3;
234 :
235 3 : if (_h)
236 3 : _f.d(_phi, *_h) = -_q1 * _q2 * _sp * sinh(_q2 / 2.0 * _sp / sf) / pow(sf, 2.0);
237 :
238 3 : if (sy)
239 1 : _f.d(_phi, *sy) = -_q1 * _q2 * _sp * sinh(_q2 / 2.0 * _sp / sf) / pow(sf, 2.0);
240 :
241 3 : if (q1)
242 1 : _f.d(_phi, *q1) = 2 * cosh(_q2 / 2.0 * _sp / sf);
243 :
244 3 : if (q2)
245 1 : _f.d(_phi, *q2) = _q1 * _sp * sinh(_q2 / 2.0 * _sp / sf) / sf;
246 :
247 3 : if (q3)
248 1 : _f.d(_phi, *q3) = -2.0 * _phi;
249 : }
250 :
251 : ////////////////////////////////////////////////////////////////////////////////////////////////
252 : //
253 : // f(se, sp, phi, h; sy, q1, q2, q3)
254 : //
255 : // h: (Optional) isotropic hardening
256 : //
257 : ////////////////////////////////////////////////////////////////////////////////////////////////
258 3 : if (_h)
259 : {
260 3 : if (_se.is_dependent())
261 3 : _f.d(*_h, _se) = -4.0 * _se / pow(sf, 3.0);
262 :
263 3 : if (_sp.is_dependent())
264 6 : _f.d(*_h, _sp) =
265 6 : -_phi * _q1 * _q2 *
266 12 : (_q2 * _sp * cosh(_q2 / 2.0 * _sp / sf) + 2.0 * sf * sinh(_q2 / 2.0 * _sp / sf)) /
267 15 : (2.0 * pow(sf, 3.0));
268 :
269 3 : if (_phi.is_dependent())
270 3 : _f.d(*_h, _phi) = -_q1 * _q2 * _sp * sinh(_q2 / 2.0 * _sp / sf) / pow(sf, 2.0);
271 :
272 12 : _f.d(*_h, *_h) = (12 * pow(Scalar(_se), 2.0) + _phi * _q1 * _q2 * _sp *
273 6 : (_q2 * _sp * cosh(_q2 / 2.0 * _sp / sf) +
274 12 : 4.0 * sf * sinh(_q2 / 2.0 * _sp / sf))) /
275 15 : (2 * pow(sf, 4.0));
276 :
277 3 : if (sy)
278 2 : _f.d(*_h, *sy) =
279 2 : (12 * pow(Scalar(_se), 2.0) +
280 2 : _phi * _q1 * _q2 * _sp *
281 4 : (_q2 * _sp * cosh(_q2 / 2.0 * _sp / sf) + 4.0 * sf * sinh(_q2 / 2.0 * _sp / sf))) /
282 5 : (2 * pow(sf, 4.0));
283 :
284 3 : if (q1)
285 1 : _f.d(*_h, *q1) = -_phi * _q2 * _sp * sinh(_q2 / 2.0 * _sp / sf) / pow(sf, 2.0);
286 :
287 3 : if (q2)
288 2 : _f.d(*_h, *q2) =
289 2 : -_phi * _q1 * _sp *
290 4 : (_q2 * _sp * cosh(_q2 / 2.0 * _sp / sf) + 2.0 * sf * sinh(_q2 / 2.0 * _sp / sf)) /
291 5 : (2 * pow(sf, 3.0));
292 : }
293 :
294 : ////////////////////////////////////////////////////////////////////////////////////////////////
295 : //
296 : // f(se, sp, phi, h; sy, q1, q2, q3)
297 : //
298 : // sy: (Optionally nonlinear) yield stress
299 : //
300 : ////////////////////////////////////////////////////////////////////////////////////////////////
301 3 : if (sy)
302 : {
303 1 : if (_se.is_dependent())
304 1 : _f.d(*sy, _se) = -4.0 * _se / pow(sf, 3.0);
305 :
306 1 : if (_phi.is_dependent())
307 1 : _f.d(*sy, _phi) = -_q1 * _q2 * _sp * sinh(_q2 / 2.0 * _sp / sf) / pow(sf, 2.0);
308 :
309 1 : if (_h)
310 2 : _f.d(*sy, *_h) =
311 2 : (12 * pow(Scalar(_se), 2.0) +
312 2 : _phi * _q1 * _q2 * _sp *
313 4 : (_q2 * _sp * cosh(_q2 / 2.0 * _sp / sf) + 4.0 * sf * sinh(_q2 / 2.0 * _sp / sf))) /
314 5 : (2 * pow(sf, 4.0));
315 :
316 4 : _f.d(*sy, *sy) = (12 * pow(Scalar(_se), 2.0) + _phi * _q1 * _q2 * _sp *
317 2 : (_q2 * _sp * cosh(_q2 / 2.0 * _sp / sf) +
318 4 : 4.0 * sf * sinh(_q2 / 2.0 * _sp / sf))) /
319 5 : (2 * pow(sf, 4.0));
320 :
321 1 : if (q1)
322 1 : _f.d(*sy, *q1) = -_phi * _q2 * _sp * sinh(_q2 / 2.0 * _sp / sf) / pow(sf, 2.0);
323 :
324 1 : if (q2)
325 2 : _f.d(*sy, *q2) =
326 2 : -_phi * _q1 * _sp *
327 4 : (_q2 * _sp * cosh(_q2 / 2.0 * _sp / sf) + 2.0 * sf * sinh(_q2 / 2.0 * _sp / sf)) /
328 5 : (2 * pow(sf, 3.0));
329 : }
330 :
331 : ////////////////////////////////////////////////////////////////////////////////////////////////
332 : //
333 : // f(se, sp, phi, h; sy, q1, q2, q3)
334 : //
335 : // q1: (Optionally nonlinear) GTN parameter q1
336 : //
337 : ////////////////////////////////////////////////////////////////////////////////////////////////
338 3 : if (q1)
339 : {
340 1 : if (_sp.is_dependent())
341 1 : _f.d(*q1, _sp) = _phi * _q2 * sinh(_q2 / 2.0 * _sp / sf) / sf;
342 :
343 1 : if (_phi.is_dependent())
344 1 : _f.d(*q1, _phi) = 2.0 * cosh(_q2 / 2.0 * _sp / sf);
345 :
346 1 : if (_h)
347 1 : _f.d(*q1, *_h) = -_phi * _q2 * _sp * sinh(_q2 / 2.0 * _sp / sf) / pow(sf, 2.0);
348 :
349 1 : if (sy)
350 1 : _f.d(*q1, *sy) = -_phi * _q2 * _sp * sinh(_q2 / 2.0 * _sp / sf) / pow(sf, 2.0);
351 :
352 1 : if (q2)
353 1 : _f.d(*q1, *q2) = _phi * _sp * sinh(_q2 / 2.0 * _sp / sf) / sf;
354 : }
355 :
356 : ////////////////////////////////////////////////////////////////////////////////////////////////
357 : //
358 : // f(se, sp, phi, h; sy, q1, q2, q3)
359 : //
360 : // q2: (Optionally nonlinear) GTN parameter q2
361 : //
362 : ////////////////////////////////////////////////////////////////////////////////////////////////
363 3 : if (q2)
364 : {
365 1 : if (_sp.is_dependent())
366 2 : _f.d(*q2, _sp) =
367 2 : _phi * _q1 *
368 4 : (_q2 * _sp * cosh(_q2 / 2.0 * _sp / sf) + 2 * sf * sinh(_q2 / 2.0 * _sp / sf)) /
369 5 : (2 * pow(sf, 2.0));
370 :
371 1 : if (_phi.is_dependent())
372 1 : _f.d(*q2, _phi) = _q1 * _sp * sinh(_q2 / 2.0 * _sp / sf) / sf;
373 :
374 1 : if (_h)
375 2 : _f.d(*q2, *_h) =
376 2 : -_phi * _q1 * _sp *
377 4 : (_q2 * _sp * cosh(_q2 / 2.0 * _sp / sf) + 2 * sf * sinh(_q2 / 2.0 * _sp / sf)) /
378 5 : (2 * pow(sf, 3.0));
379 :
380 1 : if (sy)
381 2 : _f.d(*q2, *sy) =
382 2 : -_phi * _q1 * _sp *
383 4 : (_q2 * _sp * cosh(_q2 / 2.0 * _sp / sf) + 2 * sf * sinh(_q2 / 2.0 * _sp / sf)) /
384 5 : (2 * pow(sf, 3.0));
385 :
386 1 : if (q1)
387 1 : _f.d(*q2, *q1) = _phi * _sp * sinh(_q2 / 2.0 * _sp / sf) / sf;
388 :
389 2 : _f.d(*q2, *q2) =
390 3 : _phi * _q1 * pow(Scalar(_sp), 2.0) * cosh(_q2 / 2.0 * _sp / sf) / (2 * pow(sf, 2.0));
391 : }
392 :
393 : ////////////////////////////////////////////////////////////////////////////////////////////////
394 : //
395 : // f(se, sp, phi, h; sy, q1, q2, q3)
396 : //
397 : // q3: (Optionally nonlinear) GTN parameter q3
398 : //
399 : ////////////////////////////////////////////////////////////////////////////////////////////////
400 3 : if (q3)
401 : {
402 1 : if (_phi.is_dependent())
403 1 : _f.d(*q3, _phi) = -2.0 * _phi;
404 : }
405 : }
406 9 : }
407 : } // namespace neml2
|
