Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "AbaqusUMATStress.h"
11 : #include "StepUserObject.h"
12 : #include "Factory.h"
13 : #include "MooseMesh.h"
14 : #include "RankTwoTensor.h"
15 : #include "RankFourTensor.h"
16 : #include "libmesh/int_range.h"
17 : #include <string.h>
18 : #include <algorithm>
19 :
20 : #define QUOTE(macro) stringifyName(macro)
21 :
22 : registerMooseObject("TensorMechanicsApp", AbaqusUMATStress);
23 :
24 : InputParameters
25 560 : AbaqusUMATStress::validParams()
26 : {
27 560 : InputParameters params = ComputeGeneralStressBase::validParams();
28 560 : params.addClassDescription("Coupling material to use Abaqus UMAT models in MOOSE");
29 1120 : params.addRequiredParam<FileName>(
30 : "plugin", "The path to the compiled dynamic library for the plugin you want to use");
31 1120 : params.addRequiredParam<bool>(
32 : "use_one_based_indexing",
33 : "Parameter to control whether indexing for element and integration points as presented to "
34 : "UMAT models is based on 1 (true) or 0 (false). This does not affect internal MOOSE "
35 : "numbering. The option to use 0-based numbering is deprecated and will be removed soon.");
36 1120 : params.addRequiredParam<std::vector<Real>>(
37 : "constant_properties", "Constant mechanical and thermal material properties (PROPS)");
38 1120 : params.addRequiredParam<unsigned int>("num_state_vars",
39 : "The number of state variables this UMAT is going to use");
40 1120 : params.addCoupledVar("temperature", 0.0, "Coupled temperature");
41 1120 : params.addCoupledVar("external_fields",
42 : "The external fields that can be used in the UMAT subroutine");
43 1120 : params.addParam<std::vector<MaterialPropertyName>>("external_properties", {}, "");
44 1120 : params.addParam<MooseEnum>("decomposition_method",
45 1120 : ComputeFiniteStrain::decompositionType(),
46 : "Method to calculate the strain kinematics.");
47 1120 : params.addParam<bool>(
48 : "use_displaced_mesh",
49 1120 : false,
50 : "Whether or not this object should use the "
51 : "displaced mesh for computing displacements and quantities based on the deformed state.");
52 1120 : params.addParam<UserObjectName>(
53 : "step_user_object", "The StepUserObject that provides times from simulation loading steps.");
54 560 : return params;
55 0 : }
56 :
57 : #ifndef METHOD
58 : #error "The METHOD preprocessor symbol must be supplied by the build system."
59 : #endif
60 :
61 420 : AbaqusUMATStress::AbaqusUMATStress(const InputParameters & parameters)
62 : : ComputeGeneralStressBase(parameters),
63 840 : _plugin(getParam<FileName>("plugin")),
64 840 : _library(_plugin + std::string("-") + QUOTE(METHOD) + ".plugin"),
65 420 : _umat(_library.getFunction<umat_t>("umat_")),
66 840 : _aqNSTATV(getParam<unsigned int>("num_state_vars")),
67 420 : _aqSTATEV(_aqNSTATV),
68 840 : _aqPROPS(getParam<std::vector<Real>>("constant_properties")),
69 420 : _aqNPROPS(_aqPROPS.size()),
70 840 : _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
71 840 : _total_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "total_strain")),
72 840 : _strain_increment(getOptionalMaterialProperty<RankTwoTensor>(_base_name + "strain_increment")),
73 420 : _jacobian_mult(declareProperty<RankFourTensor>(_base_name + "Jacobian_mult")),
74 840 : _Fbar(getOptionalMaterialProperty<RankTwoTensor>(_base_name + "deformation_gradient")),
75 840 : _Fbar_old(getOptionalMaterialPropertyOld<RankTwoTensor>(_base_name + "deformation_gradient")),
76 420 : _state_var(declareProperty<std::vector<Real>>(_base_name + "state_var")),
77 840 : _state_var_old(getMaterialPropertyOld<std::vector<Real>>(_base_name + "state_var")),
78 420 : _elastic_strain_energy(declareProperty<Real>(_base_name + "elastic_strain_energy")),
79 420 : _plastic_dissipation(declareProperty<Real>(_base_name + "plastic_dissipation")),
80 420 : _creep_dissipation(declareProperty<Real>(_base_name + "creep_dissipation")),
81 420 : _material_timestep(declareProperty<Real>(_base_name + "material_timestep_limit")),
82 420 : _rotation_increment(
83 420 : getOptionalMaterialProperty<RankTwoTensor>(_base_name + "rotation_increment")),
84 420 : _temperature(coupledValue("temperature")),
85 420 : _temperature_old(coupledValueOld("temperature")),
86 840 : _external_fields(isCoupled("external_fields") ? coupledValues("external_fields")
87 : : std::vector<const VariableValue *>{}),
88 840 : _external_fields_old(isCoupled("external_fields") ? coupledValuesOld("external_fields")
89 : : std::vector<const VariableValue *>{}),
90 420 : _number_external_fields(_external_fields.size()),
91 1260 : _external_property_names(getParam<std::vector<MaterialPropertyName>>("external_properties")),
92 420 : _number_external_properties(_external_property_names.size()),
93 420 : _external_properties(_number_external_properties),
94 420 : _external_properties_old(_number_external_properties),
95 840 : _use_one_based_indexing(getParam<bool>("use_one_based_indexing")),
96 420 : _decomposition_method(
97 2100 : getParam<MooseEnum>("decomposition_method").getEnum<ComputeFiniteStrain::DecompMethod>())
98 : {
99 420 : if (!_use_one_based_indexing)
100 : mooseDeprecated(
101 : "AbaqusUMATStress has transitioned to 1-based indexing in the element (NOEL) and "
102 : "integration point (NPT) numbers to ensure maximum compatibility with legacy UMAT files. "
103 : "Please ensure that any new UMAT plugins using these quantities are using the correct "
104 : "indexing. 0-based indexing will be deprecated soon.");
105 :
106 : // get material properties
107 438 : for (std::size_t i = 0; i < _number_external_properties; ++i)
108 : {
109 18 : _external_properties[i] = &getMaterialProperty<Real>(_external_property_names[i]);
110 18 : _external_properties_old[i] = &getMaterialPropertyOld<Real>(_external_property_names[i]);
111 : }
112 :
113 : // Read mesh dimension and size UMAT arrays (we always size for full 3D)
114 420 : _aqNTENS = 6; // Size of the stress or strain component array (NDI+NSHR)
115 420 : _aqNSHR = 3; // Number of engineering shear stress components
116 420 : _aqNDI = 3; // Number of direct stress components (always 3)
117 :
118 420 : _aqDDSDDT.resize(_aqNTENS);
119 420 : _aqDRPLDE.resize(_aqNTENS);
120 420 : _aqSTRAN.resize(_aqNTENS);
121 420 : _aqDFGRD0.resize(9);
122 420 : _aqDFGRD1.resize(9);
123 420 : _aqDROT.resize(9);
124 420 : _aqSTRESS.resize(_aqNTENS);
125 420 : _aqDDSDDE.resize(_aqNTENS * _aqNTENS);
126 420 : _aqDSTRAN.resize(_aqNTENS);
127 420 : _aqPREDEF.resize(_number_external_fields + _number_external_properties);
128 420 : _aqDPRED.resize(_number_external_fields + _number_external_properties);
129 420 : }
130 :
131 : void
132 418 : AbaqusUMATStress::initialSetup()
133 : {
134 : // The _Fbar, _Fbar_old, and _rotation_increment optional properties are only available when an
135 : // incremental strain formulation is used. If they are not avaliable we advide the user to
136 : // select an incremental formulation.
137 418 : if (!_strain_increment || !_Fbar || !_Fbar_old || !_rotation_increment)
138 1 : mooseError("AbaqusUMATStress '",
139 1 : name(),
140 : "': Incremental strain quantities are not available. You likely are using a total "
141 : "strain formulation. Specify `incremental = true` in the tensor mechanics action, "
142 : "or use ComputeIncrementalSmallStrain in your input file.");
143 :
144 : // Let's automatically detect uos and identify the one we are interested in.
145 : // If there is more than one, we assume something is off and error out.
146 834 : if (!isParamSetByUser("step_user_object"))
147 408 : getStepUserObject(_fe_problem, _step_user_object, name());
148 : else
149 9 : _step_user_object = &getUserObject<StepUserObject>("step_user_object");
150 417 : }
151 :
152 : void
153 23644 : AbaqusUMATStress::initQpStatefulProperties()
154 : {
155 23644 : ComputeGeneralStressBase::initQpStatefulProperties();
156 :
157 : // Initialize state variable vector
158 23644 : _state_var[_qp].resize(_aqNSTATV);
159 113516 : for (const auto i : make_range(_aqNSTATV))
160 89872 : _state_var[_qp][i] = 0.0;
161 23644 : }
162 :
163 : void
164 268362 : AbaqusUMATStress::computeProperties()
165 : {
166 : // current element "number"
167 268362 : _aqNOEL = _current_elem->id() + (_use_one_based_indexing ? 1 : 0);
168 :
169 : // characteristic element length
170 268362 : _aqCELENT = std::pow(_current_elem->volume(), 1.0 / _current_elem->dim());
171 :
172 268362 : if (!_step_user_object)
173 : {
174 : // Value of total time at the beginning of the current increment
175 267612 : _aqTIME[0] = _t - _dt;
176 : }
177 : else
178 : {
179 750 : const unsigned int start_time_step = _step_user_object->getStep(_t - _dt);
180 750 : const Real step_time = _step_user_object->getStartTime(start_time_step);
181 : // Value of step time at the beginning of the current increment
182 750 : _aqTIME[0] = step_time;
183 : }
184 : // Value of total time at the beginning of the current increment
185 268362 : _aqTIME[1] = _t - _dt;
186 :
187 : // Time increment
188 268362 : _aqDTIME = _dt;
189 :
190 : // Fill unused characters with spaces (Fortran)
191 268362 : std::fill(_aqCMNAME, _aqCMNAME + 80, ' ');
192 268362 : std::memcpy(_aqCMNAME, name().c_str(), name().size());
193 :
194 268362 : ComputeGeneralStressBase::computeProperties();
195 268362 : }
196 :
197 : void
198 1166034 : AbaqusUMATStress::computeQpStress()
199 : {
200 : // C uses row-major, whereas Fortran uses column major
201 : // therefore, all unsymmetric matrices must be transposed before passing them to Fortran
202 1166034 : RankTwoTensor FBar_old_fortran = _Fbar_old[_qp].transpose();
203 1166034 : RankTwoTensor FBar_fortran = _Fbar[_qp].transpose();
204 1166034 : RankTwoTensor DROT_fortran = _rotation_increment[_qp].transpose();
205 : const Real * myDFGRD0 = &(FBar_old_fortran(0, 0));
206 : const Real * myDFGRD1 = &(FBar_fortran(0, 0));
207 : const Real * myDROT = &(DROT_fortran(0, 0));
208 :
209 : // copy because UMAT does not guarantee constness
210 11660340 : for (const auto i : make_range(9))
211 : {
212 10494306 : _aqDFGRD0[i] = myDFGRD0[i];
213 10494306 : _aqDFGRD1[i] = myDFGRD1[i];
214 10494306 : _aqDROT[i] = myDROT[i];
215 : }
216 :
217 : // Recover "old" state variables
218 4678682 : for (const auto i : make_range(_aqNSTATV))
219 3512648 : _aqSTATEV[i] = _state_var_old[_qp][i];
220 :
221 : // Pass through updated stress, total strain, and strain increment arrays
222 : static const std::array<Real, 6> strain_factor{{1, 1, 1, 2, 2, 2}};
223 : // Account for difference in vector order convention: yz, xz, xy (MOOSE) vs xy, xz, yz
224 : // (commercial software)
225 : static const std::array<std::pair<unsigned int, unsigned int>, 6> component{
226 : {{0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}}};
227 :
228 : // rotate old stress if HughesWinget
229 1166034 : RankTwoTensor stress_old = _stress_old[_qp];
230 1166034 : if (_decomposition_method == ComputeFiniteStrain::DecompMethod::HughesWinget)
231 1056 : stress_old.rotate(_rotation_increment[_qp]);
232 :
233 8162238 : for (const auto i : make_range(_aqNTENS))
234 : {
235 6996204 : const auto a = component[i].first;
236 6996204 : const auto b = component[i].second;
237 6996204 : _aqSTRESS[i] = stress_old(a, b);
238 6996204 : _aqSTRAN[i] = _total_strain_old[_qp](a, b) * strain_factor[i];
239 6996204 : _aqDSTRAN[i] = _strain_increment[_qp](a, b) * strain_factor[i];
240 : }
241 :
242 : // current coordinates
243 4664136 : for (const auto i : make_range(Moose::dim))
244 3498102 : _aqCOORDS[i] = _q_point[_qp](i);
245 :
246 : // zero out Jacobian contribution
247 43143258 : for (const auto i : make_range(_aqNTENS * _aqNTENS))
248 41977224 : _aqDDSDDE[i] = 0.0;
249 :
250 : // Set PNEWDT initially to a large value
251 1166034 : _aqPNEWDT = std::numeric_limits<Real>::max();
252 :
253 : // Temperature
254 1166034 : _aqTEMP = _temperature_old[_qp];
255 :
256 : // Temperature increment
257 1166034 : _aqDTEMP = _temperature[_qp] - _temperature_old[_qp];
258 :
259 1575862 : for (const auto i : make_range(_number_external_fields))
260 : {
261 : // External field at this step
262 409828 : _aqPREDEF[i] = (*_external_fields_old[i])[_qp];
263 :
264 : // External field increments
265 409828 : _aqDPRED[i] = (*_external_fields[i])[_qp] - (*_external_fields_old[i])[_qp];
266 : }
267 :
268 1233842 : for (const auto i : make_range(_number_external_properties))
269 : {
270 : // External property at this step
271 67808 : _aqPREDEF[i + _number_external_fields] = (*_external_properties_old[i])[_qp];
272 :
273 : // External property increments
274 67808 : _aqDPRED[i + _number_external_fields] =
275 67808 : (*_external_properties[i])[_qp] - (*_external_properties_old[i])[_qp];
276 : }
277 :
278 : // Layer number (not supported)
279 1166034 : _aqLAYER = -1;
280 :
281 : // Section point number within the layer (not supported)
282 1166034 : _aqKSPT = -1;
283 :
284 : // Increment number
285 1166034 : _aqKINC = _t_step;
286 1166034 : _aqKSTEP = 1;
287 :
288 : // integration point number
289 1166034 : _aqNPT = _qp + (_use_one_based_indexing ? 1 : 0);
290 :
291 : // Connection to extern statement
292 1166034 : _umat(_aqSTRESS.data(),
293 : _aqSTATEV.data(),
294 : _aqDDSDDE.data(),
295 1166034 : &_elastic_strain_energy[_qp],
296 1166034 : &_plastic_dissipation[_qp],
297 1166034 : &_creep_dissipation[_qp],
298 : &_aqRPL,
299 : _aqDDSDDT.data(),
300 : _aqDRPLDE.data(),
301 : &_aqDRPLDT,
302 : _aqSTRAN.data(),
303 : _aqDSTRAN.data(),
304 : _aqTIME.data(),
305 : &_aqDTIME,
306 : &_aqTEMP,
307 : &_aqDTEMP,
308 : _aqPREDEF.data(),
309 : _aqDPRED.data(),
310 1166034 : _aqCMNAME,
311 : &_aqNDI,
312 : &_aqNSHR,
313 : &_aqNTENS,
314 : &_aqNSTATV,
315 : _aqPROPS.data(),
316 : &_aqNPROPS,
317 : _aqCOORDS.data(),
318 : _aqDROT.data(),
319 : &_aqPNEWDT,
320 : &_aqCELENT,
321 : _aqDFGRD0.data(),
322 : _aqDFGRD1.data(),
323 : &_aqNOEL,
324 : &_aqNPT,
325 : &_aqLAYER,
326 : &_aqKSPT,
327 : &_aqKSTEP,
328 : &_aqKINC);
329 :
330 : // Update state variables
331 4678682 : for (int i = 0; i < _aqNSTATV; ++i)
332 3512648 : _state_var[_qp][i] = _aqSTATEV[i];
333 :
334 : // Here, we apply UMAT convention: Always multiply _dt by PNEWDT to determine the material time
335 : // step MOOSE time stepper will choose the most limiting of all material time step increments
336 : // provided
337 1166034 : _material_timestep[_qp] = _aqPNEWDT * _dt;
338 :
339 : // Get new stress tensor - UMAT should update stress
340 : // Account for difference in vector order convention: yz, xz, xy (MOOSE) vs xy, xz, yz
341 : // (commercial software)
342 1166034 : _stress[_qp] = RankTwoTensor(
343 1166034 : _aqSTRESS[0], _aqSTRESS[1], _aqSTRESS[2], _aqSTRESS[5], _aqSTRESS[4], _aqSTRESS[3]);
344 :
345 : // Rotate the stress state to the current configuration, unless HughesWinget
346 1166034 : if (_decomposition_method != ComputeFiniteStrain::DecompMethod::HughesWinget)
347 1164978 : _stress[_qp].rotate(_rotation_increment[_qp]);
348 :
349 : // Build Jacobian matrix from UMAT's Voigt non-standard order to fourth order tensor.
350 : const unsigned int N = Moose::dim;
351 : const unsigned int ntens = N * (N + 1) / 2;
352 : const int nskip = N - 1;
353 :
354 4664136 : for (const auto i : make_range(N))
355 13992408 : for (const auto j : make_range(N))
356 41977224 : for (const auto k : make_range(N))
357 125931672 : for (const auto l : make_range(N))
358 : {
359 94448754 : if (i == j)
360 31482918 : _jacobian_mult[_qp](i, j, k, l) =
361 31482918 : k == l ? _aqDDSDDE[i * ntens + k] : _aqDDSDDE[i * ntens + k + nskip + l];
362 : else
363 : // i!=j
364 62965836 : _jacobian_mult[_qp](i, j, k, l) =
365 62965836 : k == l ? _aqDDSDDE[(nskip + i + j) * ntens + k]
366 41977224 : : _aqDDSDDE[(nskip + i + j) * ntens + k + nskip + l];
367 : }
368 1166034 : }
|
