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 1312 : AbaqusUMATStress::validParams()
26 : {
27 1312 : InputParameters params = ComputeGeneralStressBase::validParams();
28 1312 : params.addClassDescription("Coupling material to use Abaqus UMAT models in MOOSE");
29 2624 : params.addRequiredParam<FileName>(
30 : "plugin", "The path to the compiled dynamic library for the plugin you want to use");
31 2624 : 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 2624 : params.addRequiredParam<std::vector<Real>>(
37 : "constant_properties", "Constant mechanical and thermal material properties (PROPS)");
38 2624 : params.addRequiredParam<unsigned int>("num_state_vars",
39 : "The number of state variables this UMAT is going to use");
40 2624 : params.addCoupledVar("temperature", 0.0, "Coupled temperature");
41 2624 : params.addCoupledVar("external_fields",
42 : "The external fields that can be used in the UMAT subroutine");
43 2624 : params.addParam<std::vector<MaterialPropertyName>>("external_properties", {}, "");
44 2624 : params.addParam<MooseEnum>("decomposition_method",
45 2624 : ComputeFiniteStrain::decompositionType(),
46 : "Method to calculate the strain kinematics.");
47 2624 : params.addParam<bool>(
48 : "use_displaced_mesh",
49 2624 : false,
50 : "Whether or not this object should use the "
51 : "displaced mesh for computing displacements and quantities based on the deformed state.");
52 2624 : params.addParam<UserObjectName>(
53 : "analysis_step_user_object",
54 : "The AnalysisStepUserObject that provides times from simulation loading steps.");
55 2624 : params.addParam<RealVectorValue>(
56 : "orientation",
57 : "Euler angles that describe the orientation of the local material coordinate system.");
58 1312 : 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 984 : AbaqusUMATStress::AbaqusUMATStress(const InputParameters & parameters)
66 : : ComputeGeneralStressBase(parameters),
67 1968 : _plugin(getParam<FileName>("plugin")),
68 1968 : _library(_plugin + std::string("-") + QUOTE(METHOD) + ".plugin"),
69 984 : _umat(_library.getFunction<umat_t>("umat_")),
70 1968 : _aqNSTATV(getParam<unsigned int>("num_state_vars")),
71 984 : _aqSTATEV(_aqNSTATV),
72 1968 : _aqPROPS(getParam<std::vector<Real>>("constant_properties")),
73 984 : _aqNPROPS(_aqPROPS.size()),
74 1968 : _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
75 1968 : _total_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "total_strain")),
76 1968 : _strain_increment(getOptionalMaterialProperty<RankTwoTensor>(_base_name + "strain_increment")),
77 984 : _jacobian_mult(declareProperty<RankFourTensor>(_base_name + "Jacobian_mult")),
78 1968 : _Fbar(getOptionalMaterialProperty<RankTwoTensor>(_base_name + "deformation_gradient")),
79 1968 : _Fbar_old(getOptionalMaterialPropertyOld<RankTwoTensor>(_base_name + "deformation_gradient")),
80 984 : _state_var(declareProperty<std::vector<Real>>(_base_name + "state_var")),
81 1968 : _state_var_old(getMaterialPropertyOld<std::vector<Real>>(_base_name + "state_var")),
82 984 : _elastic_strain_energy(declareProperty<Real>(_base_name + "elastic_strain_energy")),
83 984 : _plastic_dissipation(declareProperty<Real>(_base_name + "plastic_dissipation")),
84 984 : _creep_dissipation(declareProperty<Real>(_base_name + "creep_dissipation")),
85 984 : _material_timestep(declareProperty<Real>(_base_name + "material_timestep_limit")),
86 984 : _rotation_increment(
87 984 : getOptionalMaterialProperty<RankTwoTensor>(_base_name + "rotation_increment")),
88 984 : _rotation_increment_old(
89 984 : getOptionalMaterialPropertyOld<RankTwoTensor>(_base_name + "rotation_increment")),
90 984 : _temperature(coupledValue("temperature")),
91 984 : _temperature_old(coupledValueOld("temperature")),
92 1968 : _external_fields(isCoupled("external_fields") ? coupledValues("external_fields")
93 : : std::vector<const VariableValue *>{}),
94 1968 : _external_fields_old(isCoupled("external_fields") ? coupledValuesOld("external_fields")
95 : : std::vector<const VariableValue *>{}),
96 984 : _number_external_fields(_external_fields.size()),
97 2952 : _external_property_names(getParam<std::vector<MaterialPropertyName>>("external_properties")),
98 984 : _number_external_properties(_external_property_names.size()),
99 984 : _external_properties(_number_external_properties),
100 984 : _external_properties_old(_number_external_properties),
101 1968 : _use_one_based_indexing(getParam<bool>("use_one_based_indexing")),
102 1968 : _use_orientation(isParamValid("orientation")),
103 1968 : _R(_use_orientation ? getParam<RealVectorValue>("orientation") : RealVectorValue(0.0)),
104 984 : _total_rotation(declareProperty<RankTwoTensor>("total_rotation")),
105 1968 : _total_rotation_old(getMaterialPropertyOld<RankTwoTensor>("total_rotation")),
106 984 : _decomposition_method(
107 4920 : getParam<MooseEnum>("decomposition_method").getEnum<ComputeFiniteStrain::DecompMethod>())
108 : {
109 984 : if (!_use_one_based_indexing)
110 36 : 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 1020 : for (std::size_t i = 0; i < _number_external_properties; ++i)
118 : {
119 36 : _external_properties[i] = &getMaterialProperty<Real>(_external_property_names[i]);
120 36 : _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 984 : _aqNTENS = 6; // Size of the stress or strain component array (NDI+NSHR)
125 984 : _aqNSHR = 3; // Number of engineering shear stress components
126 984 : _aqNDI = 3; // Number of direct stress components (always 3)
127 :
128 984 : _aqDDSDDT.resize(_aqNTENS);
129 984 : _aqDRPLDE.resize(_aqNTENS);
130 984 : _aqSTRAN.resize(_aqNTENS);
131 984 : _aqDFGRD0.resize(9);
132 984 : _aqDFGRD1.resize(9);
133 984 : _aqDROT.resize(9);
134 984 : _aqSTRESS.resize(_aqNTENS);
135 984 : _aqDDSDDE.resize(_aqNTENS * _aqNTENS);
136 984 : _aqDSTRAN.resize(_aqNTENS);
137 984 : _aqPREDEF.resize(_number_external_fields + _number_external_properties);
138 984 : _aqDPRED.resize(_number_external_fields + _number_external_properties);
139 984 : }
140 :
141 : void
142 980 : 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 980 : if (!_strain_increment || !_Fbar || !_Fbar_old || !_rotation_increment)
148 2 : mooseError("AbaqusUMATStress '",
149 2 : 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 1956 : if (!isParamSetByUser("analysis_step_user_object"))
157 960 : getAnalysisStepUserObject(_fe_problem, _step_user_object, name());
158 : else
159 18 : _step_user_object = &getUserObject<AnalysisStepUserObject>("analysis_step_user_object");
160 978 : }
161 :
162 : void
163 48864 : AbaqusUMATStress::initQpStatefulProperties()
164 : {
165 48864 : ComputeGeneralStressBase::initQpStatefulProperties();
166 :
167 : // Initialize state variable vector
168 48864 : _state_var[_qp].resize(_aqNSTATV);
169 228608 : for (const auto i : make_range(_aqNSTATV))
170 179744 : _state_var[_qp][i] = 0.0;
171 :
172 : // Initialize total rotation tensor
173 48864 : _total_rotation[_qp] = _R;
174 48864 : }
175 :
176 : void
177 458328 : AbaqusUMATStress::computeProperties()
178 : {
179 : // current element "number"
180 458328 : _aqNOEL = _current_elem->id() + (_use_one_based_indexing ? 1 : 0);
181 :
182 : // characteristic element length
183 458328 : _aqCELENT = std::pow(_current_elem->volume(), 1.0 / _current_elem->dim());
184 :
185 458328 : if (!_step_user_object)
186 : {
187 : // Value of total time at the beginning of the current increment
188 456768 : _aqTIME[0] = _t - _dt;
189 : }
190 : else
191 : {
192 1560 : const unsigned int start_time_step = _step_user_object->getStep(_t - _dt);
193 1560 : const Real step_time = _step_user_object->getStartTime(start_time_step);
194 : // Value of step time at the beginning of the current increment
195 1560 : _aqTIME[0] = step_time;
196 : }
197 : // Value of total time at the beginning of the current increment
198 458328 : _aqTIME[1] = _t - _dt;
199 :
200 : // Time increment
201 458328 : _aqDTIME = _dt;
202 :
203 : // Fill unused characters with spaces (Fortran)
204 458328 : std::fill(_aqCMNAME, _aqCMNAME + 80, ' ');
205 458328 : std::memcpy(_aqCMNAME, name().c_str(), name().size());
206 :
207 458328 : ComputeGeneralStressBase::computeProperties();
208 458328 : }
209 :
210 : void
211 2112116 : 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 2112116 : RankTwoTensor FBar_old_fortran = _Fbar_old[_qp].transpose();
216 2112116 : 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 2112116 : RankTwoTensor DROT_fortran;
221 2112116 : if (_use_orientation)
222 : {
223 32400 : DROT_fortran = RankTwoTensor::Identity();
224 : }
225 : else
226 : {
227 2079716 : if (_decomposition_method == ComputeFiniteStrain::DecompMethod::HughesWinget)
228 14368 : DROT_fortran = _rotation_increment[_qp].transpose();
229 : else
230 2065348 : 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 2112116 : RankTwoTensor stress_old = _stress_old[_qp];
239 2112116 : RankTwoTensor total_strain_old = _total_strain_old[_qp];
240 2112116 : RankTwoTensor strain_increment = _strain_increment[_qp];
241 :
242 : // check if we need to rotate to intermediate configuration
243 2112116 : if (_use_orientation)
244 : {
245 : // keep track of total rotation
246 32400 : _total_rotation[_qp] = _rotation_increment[_qp] * _total_rotation_old[_qp];
247 : // rotate stress/strain/increment from reference configuration to intermediate configuration
248 32400 : stress_old.rotate(_total_rotation_old[_qp].transpose());
249 32400 : total_strain_old.rotate(_total_rotation_old[_qp].transpose());
250 :
251 32400 : if (_decomposition_method == ComputeFiniteStrain::DecompMethod::HughesWinget)
252 25920 : strain_increment.rotate(_total_rotation[_qp].transpose());
253 : else
254 6480 : strain_increment.rotate(_total_rotation_old[_qp].transpose());
255 : }
256 2079716 : else if (_decomposition_method == ComputeFiniteStrain::DecompMethod::HughesWinget)
257 : // rotate old stress to reference configuration
258 14368 : stress_old.rotate(_rotation_increment[_qp]);
259 :
260 : // copy because UMAT does not guarantee constness
261 21121160 : for (const auto i : make_range(9))
262 : {
263 19009044 : _aqDFGRD0[i] = myDFGRD0[i];
264 19009044 : _aqDFGRD1[i] = myDFGRD1[i];
265 19009044 : _aqDROT[i] = myDROT[i];
266 : }
267 :
268 : // Recover "old" state variables
269 8062452 : for (const auto i : make_range(_aqNSTATV))
270 5950336 : _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 14784812 : for (const auto i : make_range(_aqNTENS))
280 : {
281 12672696 : const auto a = component[i].first;
282 12672696 : const auto b = component[i].second;
283 12672696 : _aqSTRESS[i] = stress_old(a, b);
284 12672696 : _aqSTRAN[i] = total_strain_old(a, b) * strain_factor[i];
285 12672696 : _aqDSTRAN[i] = strain_increment(a, b) * strain_factor[i];
286 : }
287 :
288 : // current coordinates
289 8448464 : for (const auto i : make_range(Moose::dim))
290 6336348 : _aqCOORDS[i] = _q_point[_qp](i);
291 :
292 : // zero out Jacobian contribution
293 78148292 : for (const auto i : make_range(_aqNTENS * _aqNTENS))
294 76036176 : _aqDDSDDE[i] = 0.0;
295 :
296 : // Set PNEWDT initially to a large value
297 2112116 : _aqPNEWDT = std::numeric_limits<Real>::max();
298 :
299 : // Temperature
300 2112116 : _aqTEMP = _temperature_old[_qp];
301 :
302 : // Temperature increment
303 2112116 : _aqDTEMP = _temperature[_qp] - _temperature_old[_qp];
304 :
305 2800580 : for (const auto i : make_range(_number_external_fields))
306 : {
307 : // External field at this step
308 688464 : _aqPREDEF[i] = (*_external_fields_old[i])[_qp];
309 :
310 : // External field increments
311 688464 : _aqDPRED[i] = (*_external_fields[i])[_qp] - (*_external_fields_old[i])[_qp];
312 : }
313 :
314 2242740 : for (const auto i : make_range(_number_external_properties))
315 : {
316 : // External property at this step
317 130624 : _aqPREDEF[i + _number_external_fields] = (*_external_properties_old[i])[_qp];
318 :
319 : // External property increments
320 130624 : _aqDPRED[i + _number_external_fields] =
321 130624 : (*_external_properties[i])[_qp] - (*_external_properties_old[i])[_qp];
322 : }
323 :
324 : // Layer number (not supported)
325 2112116 : _aqLAYER = -1;
326 :
327 : // Section point number within the layer (not supported)
328 2112116 : _aqKSPT = -1;
329 :
330 : // Increment number
331 2112116 : _aqKINC = _t_step;
332 2112116 : _aqKSTEP = 1;
333 :
334 : // integration point number
335 2112116 : _aqNPT = _qp + (_use_one_based_indexing ? 1 : 0);
336 :
337 : // Connection to extern statement
338 2112116 : _umat(_aqSTRESS.data(),
339 : _aqSTATEV.data(),
340 : _aqDDSDDE.data(),
341 2112116 : &_elastic_strain_energy[_qp],
342 2112116 : &_plastic_dissipation[_qp],
343 2112116 : &_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 2112116 : _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 8062452 : for (int i = 0; i < _aqNSTATV; ++i)
378 5950336 : _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 2112116 : _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 2112116 : _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 8448464 : for (const auto i : make_range(N))
397 25345392 : for (const auto j : make_range(N))
398 76036176 : for (const auto k : make_range(N))
399 228108528 : for (const auto l : make_range(N))
400 : {
401 171081396 : if (i == j)
402 57027132 : _jacobian_mult[_qp](i, j, k, l) =
403 57027132 : k == l ? _aqDDSDDE[i * ntens + k] : _aqDDSDDE[i * ntens + k + nskip + l];
404 : else
405 : // i!=j
406 114054264 : _jacobian_mult[_qp](i, j, k, l) =
407 114054264 : k == l ? _aqDDSDDE[(nskip + i + j) * ntens + k]
408 76036176 : : _aqDDSDDE[(nskip + i + j) * ntens + k + nskip + l];
409 : }
410 :
411 : // check if we need to rotate from intermediate reference frame
412 2112116 : if (_use_orientation)
413 : {
414 : // rotate to current configuration
415 32400 : _stress[_qp].rotate(_total_rotation[_qp]);
416 32400 : _jacobian_mult[_qp].rotate(_total_rotation[_qp]);
417 : }
418 2079716 : else if (_decomposition_method != ComputeFiniteStrain::DecompMethod::HughesWinget)
419 2065348 : _stress[_qp].rotate(_rotation_increment[_qp]);
420 2112116 : }
|
