Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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 "AnalysisStepUserObject.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("SolidMechanicsApp", AbaqusUMATStress);
23 :
24 : InputParameters
25 832 : AbaqusUMATStress::validParams()
26 : {
27 832 : InputParameters params = ComputeGeneralStressBase::validParams();
28 832 : params.addClassDescription("Coupling material to use Abaqus UMAT models in MOOSE");
29 1664 : params.addRequiredParam<FileName>(
30 : "plugin", "The path to the compiled dynamic library for the plugin you want to use");
31 1664 : 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 1664 : params.addRequiredParam<std::vector<Real>>(
37 : "constant_properties", "Constant mechanical and thermal material properties (PROPS)");
38 1664 : params.addRequiredParam<unsigned int>("num_state_vars",
39 : "The number of state variables this UMAT is going to use");
40 1664 : params.addCoupledVar("temperature", 0.0, "Coupled temperature");
41 1664 : params.addCoupledVar("external_fields",
42 : "The external fields that can be used in the UMAT subroutine");
43 1664 : params.addParam<std::vector<MaterialPropertyName>>("external_properties", {}, "");
44 1664 : params.addParam<MooseEnum>("decomposition_method",
45 1664 : ComputeFiniteStrain::decompositionType(),
46 : "Method to calculate the strain kinematics.");
47 1664 : params.addParam<bool>(
48 : "use_displaced_mesh",
49 1664 : false,
50 : "Whether or not this object should use the "
51 : "displaced mesh for computing displacements and quantities based on the deformed state.");
52 1664 : params.addParam<UserObjectName>(
53 : "analysis_step_user_object",
54 : "The AnalysisStepUserObject that provides times from simulation loading steps.");
55 1664 : params.addParam<RealVectorValue>(
56 : "orientation",
57 : "Euler angles that describe the orientation of the local material coordinate system.");
58 832 : return params;
59 0 : }
60 :
61 : #ifndef METHOD
62 : #error "The METHOD preprocessor symbol must be supplied by the build system."
63 : #endif
64 :
65 624 : AbaqusUMATStress::AbaqusUMATStress(const InputParameters & parameters)
66 : : ComputeGeneralStressBase(parameters),
67 1248 : _plugin(getParam<FileName>("plugin")),
68 1248 : _library(_plugin + std::string("-") + QUOTE(METHOD) + ".plugin"),
69 624 : _umat(_library.getFunction<umat_t>("umat_")),
70 1248 : _aqNSTATV(getParam<unsigned int>("num_state_vars")),
71 624 : _aqSTATEV(_aqNSTATV),
72 1872 : _aqPROPS(getParam<std::vector<Real>>("constant_properties")),
73 624 : _aqNPROPS(_aqPROPS.size()),
74 1248 : _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
75 1248 : _total_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "total_strain")),
76 1248 : _strain_increment(getOptionalMaterialProperty<RankTwoTensor>(_base_name + "strain_increment")),
77 624 : _jacobian_mult(declareProperty<RankFourTensor>(_base_name + "Jacobian_mult")),
78 1248 : _Fbar(getOptionalMaterialProperty<RankTwoTensor>(_base_name + "deformation_gradient")),
79 1248 : _Fbar_old(getOptionalMaterialPropertyOld<RankTwoTensor>(_base_name + "deformation_gradient")),
80 624 : _state_var(declareProperty<std::vector<Real>>(_base_name + "state_var")),
81 1248 : _state_var_old(getMaterialPropertyOld<std::vector<Real>>(_base_name + "state_var")),
82 624 : _elastic_strain_energy(declareProperty<Real>(_base_name + "elastic_strain_energy")),
83 1248 : _elastic_strain_energy_old(getMaterialPropertyOld<Real>(_base_name + "elastic_strain_energy")),
84 624 : _plastic_dissipation(declareProperty<Real>(_base_name + "plastic_dissipation")),
85 1248 : _plastic_dissipation_old(getMaterialPropertyOld<Real>(_base_name + "plastic_dissipation")),
86 624 : _creep_dissipation(declareProperty<Real>(_base_name + "creep_dissipation")),
87 1248 : _creep_dissipation_old(getMaterialPropertyOld<Real>(_base_name + "creep_dissipation")),
88 624 : _material_timestep(declareProperty<Real>(_base_name + "material_timestep_limit")),
89 624 : _rotation_increment(
90 624 : getOptionalMaterialProperty<RankTwoTensor>(_base_name + "rotation_increment")),
91 624 : _rotation_increment_old(
92 624 : getOptionalMaterialPropertyOld<RankTwoTensor>(_base_name + "rotation_increment")),
93 624 : _temperature(coupledValue("temperature")),
94 624 : _temperature_old(coupledValueOld("temperature")),
95 1248 : _external_fields(isCoupled("external_fields") ? coupledValues("external_fields")
96 : : std::vector<const VariableValue *>{}),
97 1248 : _external_fields_old(isCoupled("external_fields") ? coupledValuesOld("external_fields")
98 : : std::vector<const VariableValue *>{}),
99 624 : _number_external_fields(_external_fields.size()),
100 1872 : _external_property_names(getParam<std::vector<MaterialPropertyName>>("external_properties")),
101 624 : _number_external_properties(_external_property_names.size()),
102 624 : _external_properties(_number_external_properties),
103 624 : _external_properties_old(_number_external_properties),
104 1248 : _use_one_based_indexing(getParam<bool>("use_one_based_indexing")),
105 1248 : _use_orientation(isParamValid("orientation")),
106 1248 : _R(_use_orientation ? getParam<RealVectorValue>("orientation") : RealVectorValue(0.0)),
107 624 : _total_rotation(declareProperty<RankTwoTensor>("total_rotation")),
108 1248 : _total_rotation_old(getMaterialPropertyOld<RankTwoTensor>("total_rotation")),
109 624 : _decomposition_method(
110 1248 : getParam<MooseEnum>("decomposition_method").getEnum<ComputeFiniteStrain::DecompMethod>())
111 : {
112 624 : if (!_use_one_based_indexing)
113 24 : mooseDeprecated(
114 : "AbaqusUMATStress has transitioned to 1-based indexing in the element (NOEL) and "
115 : "integration point (NPT) numbers to ensure maximum compatibility with legacy UMAT files. "
116 : "Please ensure that any new UMAT plugins using these quantities are using the correct "
117 : "indexing. 0-based indexing will be deprecated soon.");
118 :
119 : // get material properties
120 648 : for (std::size_t i = 0; i < _number_external_properties; ++i)
121 : {
122 24 : _external_properties[i] = &getMaterialProperty<Real>(_external_property_names[i]);
123 24 : _external_properties_old[i] = &getMaterialPropertyOld<Real>(_external_property_names[i]);
124 : }
125 :
126 : // Read mesh dimension and size UMAT arrays (we always size for full 3D)
127 624 : _aqNTENS = 6; // Size of the stress or strain component array (NDI+NSHR)
128 624 : _aqNSHR = 3; // Number of engineering shear stress components
129 624 : _aqNDI = 3; // Number of direct stress components (always 3)
130 :
131 624 : _aqDDSDDT.resize(_aqNTENS);
132 624 : _aqDRPLDE.resize(_aqNTENS);
133 624 : _aqSTRAN.resize(_aqNTENS);
134 624 : _aqDFGRD0.resize(9);
135 624 : _aqDFGRD1.resize(9);
136 624 : _aqDROT.resize(9);
137 624 : _aqSTRESS.resize(_aqNTENS);
138 624 : _aqDDSDDE.resize(_aqNTENS * _aqNTENS);
139 624 : _aqDSTRAN.resize(_aqNTENS);
140 624 : _aqPREDEF.resize(_number_external_fields + _number_external_properties);
141 624 : _aqDPRED.resize(_number_external_fields + _number_external_properties);
142 624 : }
143 :
144 : void
145 622 : AbaqusUMATStress::initialSetup()
146 : {
147 : // The _Fbar, _Fbar_old, and _rotation_increment optional properties are only available when an
148 : // incremental strain formulation is used. If they are not avaliable we advide the user to
149 : // select an incremental formulation.
150 622 : if (!_strain_increment || !_Fbar || !_Fbar_old || !_rotation_increment)
151 1 : mooseError("AbaqusUMATStress '",
152 : name(),
153 : "': Incremental strain quantities are not available. You likely are using a total "
154 : "strain formulation. Specify `incremental = true` in the tensor mechanics action, "
155 : "or use ComputeIncrementalStrain in your input file.");
156 :
157 : // Let's automatically detect uos and identify the one we are interested in.
158 : // If there is more than one, we assume something is off and error out.
159 1242 : if (!isParamSetByUser("analysis_step_user_object"))
160 609 : getAnalysisStepUserObject(_fe_problem, _step_user_object, name());
161 : else
162 12 : _step_user_object = &getUserObject<AnalysisStepUserObject>("analysis_step_user_object");
163 621 : }
164 :
165 : void
166 36082 : AbaqusUMATStress::initQpStatefulProperties()
167 : {
168 36082 : ComputeGeneralStressBase::initQpStatefulProperties();
169 :
170 : // Initialize state variable vector
171 36082 : _state_var[_qp].resize(_aqNSTATV);
172 170890 : for (const auto i : make_range(_aqNSTATV))
173 134808 : _state_var[_qp][i] = 0.0;
174 :
175 : // Initialize total rotation tensor
176 36082 : _total_rotation[_qp] = _R;
177 36082 : }
178 :
179 : void
180 340534 : AbaqusUMATStress::computeProperties()
181 : {
182 : // current element "number"
183 340534 : _aqNOEL = _current_elem->id() + (_use_one_based_indexing ? 1 : 0);
184 :
185 : // characteristic element length
186 340534 : _aqCELENT = std::pow(_current_elem->volume(), 1.0 / _current_elem->dim());
187 :
188 340534 : if (!_step_user_object)
189 : {
190 : // Value of total time at the beginning of the current increment
191 339352 : _aqTIME[0] = _t - _dt;
192 : }
193 : else
194 : {
195 1182 : const unsigned int start_time_step = _step_user_object->getStep(_t - _dt);
196 1182 : const Real step_time = _step_user_object->getStartTime(start_time_step);
197 : // Value of step time at the beginning of the current increment
198 1182 : _aqTIME[0] = step_time;
199 : }
200 : // Value of total time at the beginning of the current increment
201 340534 : _aqTIME[1] = _t - _dt;
202 :
203 : // Time increment
204 340534 : _aqDTIME = _dt;
205 :
206 : // Fill unused characters with spaces (Fortran)
207 340534 : std::fill(_aqCMNAME, _aqCMNAME + 80, ' ');
208 340534 : std::memcpy(_aqCMNAME, name().c_str(), name().size());
209 :
210 340534 : ComputeGeneralStressBase::computeProperties();
211 340534 : }
212 :
213 : void
214 1574278 : AbaqusUMATStress::computeQpStress()
215 : {
216 : // C uses row-major, whereas Fortran uses column major
217 : // therefore, all unsymmetric matrices must be transposed before passing them to Fortran
218 1574278 : RankTwoTensor FBar_old_fortran = _Fbar_old[_qp].transpose();
219 1574278 : RankTwoTensor FBar_fortran = _Fbar[_qp].transpose();
220 :
221 : // DROT needed by UMAT will depend on kinematics and whether or not an intermediate configuration
222 : // is used
223 1574278 : RankTwoTensor DROT_fortran;
224 1574278 : if (_use_orientation)
225 : {
226 22032 : DROT_fortran = RankTwoTensor::Identity();
227 : }
228 : else
229 : {
230 1552246 : if (_decomposition_method == ComputeFiniteStrain::DecompMethod::HughesWinget)
231 9512 : DROT_fortran = _rotation_increment[_qp].transpose();
232 : else
233 1542734 : DROT_fortran = _rotation_increment_old[_qp].transpose();
234 : }
235 :
236 : const Real * myDFGRD0 = &(FBar_old_fortran(0, 0));
237 : const Real * myDFGRD1 = &(FBar_fortran(0, 0));
238 : const Real * myDROT = &(DROT_fortran(0, 0));
239 :
240 : // More local copies of materials so we can (optionally) rotate
241 1574278 : RankTwoTensor stress_old = _stress_old[_qp];
242 1574278 : RankTwoTensor total_strain_old = _total_strain_old[_qp];
243 1574278 : RankTwoTensor strain_increment = _strain_increment[_qp];
244 :
245 : // check if we need to rotate to intermediate configuration
246 1574278 : if (_use_orientation)
247 : {
248 : // keep track of total rotation
249 22032 : _total_rotation[_qp] = _rotation_increment[_qp] * _total_rotation_old[_qp];
250 : // rotate stress/strain/increment from reference configuration to intermediate configuration
251 22032 : stress_old.rotate(_total_rotation_old[_qp].transpose());
252 22032 : total_strain_old.rotate(_total_rotation_old[_qp].transpose());
253 :
254 22032 : if (_decomposition_method == ComputeFiniteStrain::DecompMethod::HughesWinget)
255 16848 : strain_increment.rotate(_total_rotation[_qp].transpose());
256 : else
257 5184 : strain_increment.rotate(_total_rotation_old[_qp].transpose());
258 : }
259 1552246 : else if (_decomposition_method == ComputeFiniteStrain::DecompMethod::HughesWinget)
260 : // rotate old stress to reference configuration
261 9512 : stress_old.rotate(_rotation_increment[_qp]);
262 :
263 : // copy because UMAT does not guarantee constness
264 15742780 : for (const auto i : make_range(9))
265 : {
266 14168502 : _aqDFGRD0[i] = myDFGRD0[i];
267 14168502 : _aqDFGRD1[i] = myDFGRD1[i];
268 14168502 : _aqDROT[i] = myDROT[i];
269 : }
270 :
271 : // Recover "old" state variables
272 6176570 : for (const auto i : make_range(_aqNSTATV))
273 4602292 : _aqSTATEV[i] = _state_var_old[_qp][i];
274 :
275 : // Recover "old" energy quantities
276 1574278 : _elastic_strain_energy[_qp] = _elastic_strain_energy_old[_qp];
277 1574278 : _plastic_dissipation[_qp] = _plastic_dissipation_old[_qp];
278 1574278 : _creep_dissipation[_qp] = _creep_dissipation_old[_qp];
279 :
280 : // Pass through updated stress, total strain, and strain increment arrays
281 : static const std::array<Real, 6> strain_factor{{1, 1, 1, 2, 2, 2}};
282 : // Account for difference in vector order convention: yz, xz, xy (MOOSE) vs xy, xz, yz
283 : // (commercial software)
284 : static const std::array<std::pair<unsigned int, unsigned int>, 6> component{
285 : {{0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}}};
286 :
287 11019946 : for (const auto i : make_range(_aqNTENS))
288 : {
289 9445668 : const auto a = component[i].first;
290 9445668 : const auto b = component[i].second;
291 9445668 : _aqSTRESS[i] = stress_old(a, b);
292 9445668 : _aqSTRAN[i] = total_strain_old(a, b) * strain_factor[i];
293 9445668 : _aqDSTRAN[i] = strain_increment(a, b) * strain_factor[i];
294 : }
295 :
296 : // current coordinates
297 6297112 : for (const auto i : make_range(Moose::dim))
298 4722834 : _aqCOORDS[i] = _q_point[_qp](i);
299 :
300 : // zero out Jacobian contribution
301 58248286 : for (const auto i : make_range(_aqNTENS * _aqNTENS))
302 56674008 : _aqDDSDDE[i] = 0.0;
303 :
304 : // Set PNEWDT initially to a large value
305 1574278 : _aqPNEWDT = std::numeric_limits<Real>::max();
306 :
307 : // Temperature
308 1574278 : _aqTEMP = _temperature_old[_qp];
309 :
310 : // Temperature increment
311 1574278 : _aqDTEMP = _temperature[_qp] - _temperature_old[_qp];
312 :
313 2058232 : for (const auto i : make_range(_number_external_fields))
314 : {
315 : // External field at this step
316 483954 : _aqPREDEF[i] = (*_external_fields_old[i])[_qp];
317 :
318 : // External field increments
319 483954 : _aqDPRED[i] = (*_external_fields[i])[_qp] - (*_external_fields_old[i])[_qp];
320 : }
321 :
322 1638502 : for (const auto i : make_range(_number_external_properties))
323 : {
324 : // External property at this step
325 64224 : _aqPREDEF[i + _number_external_fields] = (*_external_properties_old[i])[_qp];
326 :
327 : // External property increments
328 64224 : _aqDPRED[i + _number_external_fields] =
329 64224 : (*_external_properties[i])[_qp] - (*_external_properties_old[i])[_qp];
330 : }
331 :
332 : // Layer number (not supported)
333 1574278 : _aqLAYER = -1;
334 :
335 : // Section point number within the layer (not supported)
336 1574278 : _aqKSPT = -1;
337 :
338 : // Increment number
339 1574278 : _aqKINC = _t_step;
340 1574278 : _aqKSTEP = 1;
341 :
342 : // integration point number
343 1574278 : _aqNPT = _qp + (_use_one_based_indexing ? 1 : 0);
344 :
345 : // Connection to extern statement
346 1574278 : _umat(_aqSTRESS.data(),
347 : _aqSTATEV.data(),
348 : _aqDDSDDE.data(),
349 : &_elastic_strain_energy[_qp],
350 : &_plastic_dissipation[_qp],
351 : &_creep_dissipation[_qp],
352 : &_aqRPL,
353 : _aqDDSDDT.data(),
354 : _aqDRPLDE.data(),
355 : &_aqDRPLDT,
356 : _aqSTRAN.data(),
357 : _aqDSTRAN.data(),
358 : _aqTIME.data(),
359 : &_aqDTIME,
360 : &_aqTEMP,
361 : &_aqDTEMP,
362 : _aqPREDEF.data(),
363 : _aqDPRED.data(),
364 1574278 : _aqCMNAME,
365 : &_aqNDI,
366 : &_aqNSHR,
367 : &_aqNTENS,
368 : &_aqNSTATV,
369 : _aqPROPS.data(),
370 : &_aqNPROPS,
371 : _aqCOORDS.data(),
372 : _aqDROT.data(),
373 : &_aqPNEWDT,
374 : &_aqCELENT,
375 : _aqDFGRD0.data(),
376 : _aqDFGRD1.data(),
377 : &_aqNOEL,
378 : &_aqNPT,
379 : &_aqLAYER,
380 : &_aqKSPT,
381 : &_aqKSTEP,
382 : &_aqKINC);
383 :
384 : // Update state variables
385 6176570 : for (int i = 0; i < _aqNSTATV; ++i)
386 4602292 : _state_var[_qp][i] = _aqSTATEV[i];
387 :
388 : // Here, we apply UMAT convention: Always multiply _dt by PNEWDT to determine the material time
389 : // step MOOSE time stepper will choose the most limiting of all material time step increments
390 : // provided
391 1574278 : _material_timestep[_qp] = _aqPNEWDT * _dt;
392 :
393 : // Get new stress tensor - UMAT should update stress
394 : // Account for difference in vector order convention: yz, xz, xy (MOOSE) vs xy, xz, yz
395 : // (commercial software)
396 1574278 : _stress[_qp] = RankTwoTensor(
397 : _aqSTRESS[0], _aqSTRESS[1], _aqSTRESS[2], _aqSTRESS[5], _aqSTRESS[4], _aqSTRESS[3]);
398 :
399 : // Build Jacobian matrix from UMAT's Voigt non-standard order to fourth order tensor.
400 : const unsigned int N = Moose::dim;
401 : const unsigned int ntens = N * (N + 1) / 2;
402 : const int nskip = N - 1;
403 :
404 6297112 : for (const auto i : make_range(N))
405 18891336 : for (const auto j : make_range(N))
406 56674008 : for (const auto k : make_range(N))
407 170022024 : for (const auto l : make_range(N))
408 : {
409 127516518 : if (i == j)
410 42505506 : _jacobian_mult[_qp](i, j, k, l) =
411 42505506 : k == l ? _aqDDSDDE[i * ntens + k] : _aqDDSDDE[i * ntens + k + nskip + l];
412 : else
413 : // i!=j
414 85011012 : _jacobian_mult[_qp](i, j, k, l) =
415 85011012 : k == l ? _aqDDSDDE[(nskip + i + j) * ntens + k]
416 56674008 : : _aqDDSDDE[(nskip + i + j) * ntens + k + nskip + l];
417 : }
418 :
419 : // check if we need to rotate from intermediate reference frame
420 1574278 : if (_use_orientation)
421 : {
422 : // rotate to current configuration
423 22032 : _stress[_qp].rotate(_total_rotation[_qp]);
424 22032 : _jacobian_mult[_qp].rotate(_total_rotation[_qp]);
425 : }
426 1552246 : else if (_decomposition_method != ComputeFiniteStrain::DecompMethod::HughesWinget)
427 1542734 : _stress[_qp].rotate(_rotation_increment[_qp]);
428 1574278 : }
|
