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