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 "StressDivergenceTensors.h"
11 :
12 : // MOOSE includes
13 : #include "ElasticityTensorTools.h"
14 : #include "Material.h"
15 : #include "MooseMesh.h"
16 : #include "MooseVariable.h"
17 : #include "SystemBase.h"
18 :
19 : #include "libmesh/quadrature.h"
20 :
21 : registerMooseObject("TensorMechanicsApp", StressDivergenceTensors);
22 :
23 : InputParameters
24 15607 : StressDivergenceTensors::validParams()
25 : {
26 15607 : InputParameters params = JvarMapKernelInterface<ALEKernel>::validParams();
27 15607 : params.addClassDescription("Stress divergence kernel for the Cartesian coordinate system");
28 31214 : params.addRequiredRangeCheckedParam<unsigned int>("component",
29 : "component < 3",
30 : "An integer corresponding to the direction "
31 : "the variable this kernel acts in. (0 for x, "
32 : "1 for y, 2 for z)");
33 31214 : params.addRequiredCoupledVar("displacements",
34 : "The string of displacements suitable for the problem statement");
35 :
36 : // maybe this should be deprecated in favor of args
37 31214 : params.addCoupledVar("temperature",
38 : "The name of the temperature variable used in the "
39 : "ComputeThermalExpansionEigenstrain. (Not required for "
40 : "simulations without temperature coupling.)");
41 :
42 31214 : params.addParam<std::vector<MaterialPropertyName>>(
43 : "eigenstrain_names",
44 : {},
45 : "List of eigenstrains used in the strain calculation. Used for computing their derivatives "
46 : "for off-diagonal Jacobian terms.");
47 31214 : params.addCoupledVar("out_of_plane_strain",
48 : "The name of the out_of_plane_strain variable used in the "
49 : "WeakPlaneStress kernel.");
50 31214 : MooseEnum out_of_plane_direction("x y z", "z");
51 31214 : params.addParam<MooseEnum>(
52 : "out_of_plane_direction",
53 : out_of_plane_direction,
54 : "The direction of the out_of_plane_strain variable used in the WeakPlaneStress kernel.");
55 31214 : params.addParam<std::string>("base_name", "Material property base name");
56 15607 : params.set<bool>("use_displaced_mesh") = false;
57 31214 : params.addParam<bool>(
58 31214 : "use_finite_deform_jacobian", false, "Jacobian for corotational finite strain");
59 31214 : params.addParam<bool>("volumetric_locking_correction",
60 31214 : false,
61 : "Set to false to turn off volumetric locking correction");
62 15607 : return params;
63 15607 : }
64 :
65 7637 : StressDivergenceTensors::StressDivergenceTensors(const InputParameters & parameters)
66 : : JvarMapKernelInterface<ALEKernel>(parameters),
67 8357 : _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
68 15274 : _use_finite_deform_jacobian(getParam<bool>("use_finite_deform_jacobian")),
69 15274 : _stress(getMaterialPropertyByName<RankTwoTensor>(_base_name + "stress")),
70 22911 : _Jacobian_mult(getMaterialPropertyByName<RankFourTensor>(_base_name + "Jacobian_mult")),
71 15274 : _component(getParam<unsigned int>("component")),
72 7637 : _ndisp(coupledComponents("displacements")),
73 7637 : _disp_var(_ndisp),
74 7637 : _out_of_plane_strain_coupled(isCoupled("out_of_plane_strain")),
75 7637 : _out_of_plane_strain(_out_of_plane_strain_coupled ? &coupledValue("out_of_plane_strain")
76 : : nullptr),
77 7637 : _out_of_plane_strain_var(_out_of_plane_strain_coupled ? coupled("out_of_plane_strain") : 0),
78 15274 : _out_of_plane_direction(getParam<MooseEnum>("out_of_plane_direction")),
79 15274 : _use_displaced_mesh(getParam<bool>("use_displaced_mesh")),
80 7637 : _avg_grad_test(_test.size(), std::vector<Real>(3, 0.0)),
81 7637 : _avg_grad_phi(_phi.size(), std::vector<Real>(3, 0.0)),
82 30548 : _volumetric_locking_correction(getParam<bool>("volumetric_locking_correction"))
83 : {
84 : // get coupled displacements
85 29264 : for (unsigned int i = 0; i < _ndisp; ++i)
86 21627 : _disp_var[i] = coupled("displacements", i);
87 :
88 : // fetch eigenstrain derivatives
89 : const auto nvar = _coupled_moose_vars.size();
90 7637 : _deigenstrain_dargs.resize(nvar);
91 30660 : for (std::size_t i = 0; i < nvar; ++i)
92 47306 : for (auto eigenstrain_name : getParam<std::vector<MaterialPropertyName>>("eigenstrain_names"))
93 3780 : _deigenstrain_dargs[i].push_back(&getMaterialPropertyDerivative<RankTwoTensor>(
94 1260 : eigenstrain_name, _coupled_moose_vars[i]->name()));
95 :
96 : // Checking for consistency between mesh size and length of the provided displacements vector
97 7637 : if (_out_of_plane_direction != 2 && _ndisp != 3)
98 0 : mooseError("For 2D simulations where the out-of-plane direction is x or y coordinate "
99 : "directions the number of supplied displacements must be three.");
100 7637 : else if (_out_of_plane_direction == 2 && _ndisp != _mesh.dimension())
101 0 : mooseError("The number of displacement variables supplied must match the mesh dimension");
102 :
103 7637 : if (_use_finite_deform_jacobian)
104 : {
105 140 : _deformation_gradient =
106 140 : &getMaterialProperty<RankTwoTensor>(_base_name + "deformation_gradient");
107 140 : _deformation_gradient_old =
108 140 : &getMaterialPropertyOld<RankTwoTensor>(_base_name + "deformation_gradient");
109 280 : _rotation_increment = &getMaterialProperty<RankTwoTensor>(_base_name + "rotation_increment");
110 : }
111 :
112 : // Error if volumetric locking correction is turned on for 1D problems
113 7637 : if (_ndisp == 1 && _volumetric_locking_correction)
114 0 : mooseError("Volumetric locking correction should be set to false for 1-D problems.");
115 :
116 : // Generate warning when volumetric locking correction is used with second order elements
117 7637 : if (_mesh.hasSecondOrderElements() && _volumetric_locking_correction)
118 0 : mooseWarning("Volumteric locking correction is not required for second order elements. Using "
119 : "volumetric locking with second order elements could cause zigzag patterns in "
120 : "stresses and strains.");
121 7637 : }
122 :
123 : void
124 7216 : StressDivergenceTensors::initialSetup()
125 : {
126 : // check if any of the eigenstrains provide derivatives wrt variables that are not coupled
127 14831 : for (auto eigenstrain_name : getParam<std::vector<MaterialPropertyName>>("eigenstrain_names"))
128 1197 : validateNonlinearCoupling<RankTwoTensor>(eigenstrain_name);
129 :
130 : // make sure the coordinate system is cartesioan
131 7216 : if (getBlockCoordSystem() != Moose::COORD_XYZ)
132 0 : mooseError("The coordinate system in the Problem block must be set to XYZ for cartesian "
133 : "geometries.");
134 7216 : }
135 :
136 : void
137 8857279 : StressDivergenceTensors::computeResidual()
138 : {
139 8857279 : prepareVectorTag(_assembly, _var.number());
140 :
141 8857279 : if (_volumetric_locking_correction)
142 4645524 : computeAverageGradientTest();
143 :
144 8857279 : precalculateResidual();
145 68896495 : for (_i = 0; _i < _test.size(); ++_i)
146 533251164 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
147 473211948 : _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
148 :
149 8857279 : accumulateTaggedLocalResidual();
150 :
151 8857279 : if (_has_save_in)
152 : {
153 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
154 63904 : for (const auto & var : _save_in)
155 31952 : var->sys().solution().add_vector(_local_re, var->dofIndices());
156 : }
157 8857279 : }
158 :
159 : Real
160 453179328 : StressDivergenceTensors::computeQpResidual()
161 : {
162 453179328 : Real residual = _stress[_qp].row(_component) * _grad_test[_i][_qp];
163 : // volumetric locking correction
164 453179328 : if (_volumetric_locking_correction)
165 206716480 : residual += _stress[_qp].trace() / 3.0 *
166 206716480 : (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component));
167 :
168 453179328 : if (_ndisp != 3 && _out_of_plane_strain_coupled && _use_displaced_mesh)
169 : {
170 48064 : const Real out_of_plane_thickness = std::exp((*_out_of_plane_strain)[_qp]);
171 48064 : residual *= out_of_plane_thickness;
172 : }
173 :
174 453179328 : return residual;
175 : }
176 :
177 : void
178 911271 : StressDivergenceTensors::computeJacobian()
179 : {
180 911271 : if (_volumetric_locking_correction)
181 : {
182 369604 : computeAverageGradientTest();
183 369604 : computeAverageGradientPhi();
184 : }
185 :
186 911271 : if (_use_finite_deform_jacobian)
187 : {
188 13290 : _finite_deform_Jacobian_mult.resize(_qrule->n_points());
189 :
190 116090 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
191 102800 : computeFiniteDeformJacobian();
192 :
193 13290 : ALEKernel::computeJacobian();
194 : }
195 : else
196 897981 : Kernel::computeJacobian();
197 911271 : }
198 :
199 : void
200 1535048 : StressDivergenceTensors::computeOffDiagJacobian(const unsigned int jvar)
201 : {
202 1535048 : if (_volumetric_locking_correction)
203 : {
204 422054 : computeAverageGradientPhi();
205 422054 : computeAverageGradientTest();
206 : }
207 :
208 1535048 : if (_use_finite_deform_jacobian)
209 : {
210 26028 : _finite_deform_Jacobian_mult.resize(_qrule->n_points());
211 :
212 227212 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
213 201184 : computeFiniteDeformJacobian();
214 :
215 26028 : ALEKernel::computeOffDiagJacobian(jvar);
216 : }
217 : else
218 1509020 : Kernel::computeOffDiagJacobian(jvar);
219 1535048 : }
220 :
221 : Real
222 441230586 : StressDivergenceTensors::computeQpJacobian()
223 : {
224 441230586 : if (_use_finite_deform_jacobian)
225 6578240 : return ElasticityTensorTools::elasticJacobian(_finite_deform_Jacobian_mult[_qp],
226 : _component,
227 6578240 : _component,
228 6578240 : _grad_test[_i][_qp],
229 6578240 : _grad_phi_undisplaced[_j][_qp]);
230 :
231 434652346 : Real sum_C3x3 = _Jacobian_mult[_qp].sum3x3();
232 434652346 : RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum3x1();
233 :
234 : Real jacobian = 0.0;
235 : // B^T_i * C * B_j
236 434652346 : jacobian += ElasticityTensorTools::elasticJacobian(
237 434652346 : _Jacobian_mult[_qp], _component, _component, _grad_test[_i][_qp], _grad_phi[_j][_qp]);
238 :
239 434652346 : if (_volumetric_locking_correction)
240 : {
241 : // jacobian = Bbar^T_i * C * Bbar_j where Bbar = B + Bvol
242 : // jacobian = B^T_i * C * B_j + Bvol^T_i * C * Bvol_j + Bvol^T_i * C * B_j + B^T_i * C *
243 : // Bvol_j
244 :
245 : // Bvol^T_i * C * Bvol_j
246 114005248 : jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
247 114005248 : (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 9.0;
248 :
249 : // B^T_i * C * Bvol_j
250 114005248 : jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
251 114005248 : (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 3.0;
252 :
253 : // Bvol^T_i * C * B_j
254 114005248 : RankTwoTensor phi;
255 114005248 : switch (_component)
256 : {
257 39487616 : case 0:
258 39487616 : phi(0, 0) = _grad_phi[_j][_qp](0);
259 39487616 : phi(0, 1) = phi(1, 0) = _grad_phi[_j][_qp](1);
260 39487616 : phi(0, 2) = phi(2, 0) = _grad_phi[_j][_qp](2);
261 39487616 : break;
262 :
263 39487616 : case 1:
264 39487616 : phi(1, 1) = _grad_phi[_j][_qp](1);
265 39487616 : phi(0, 1) = phi(1, 0) = _grad_phi[_j][_qp](0);
266 39487616 : phi(1, 2) = phi(2, 1) = _grad_phi[_j][_qp](2);
267 39487616 : break;
268 :
269 35030016 : case 2:
270 35030016 : phi(2, 2) = _grad_phi[_j][_qp](2);
271 35030016 : phi(0, 2) = phi(2, 0) = _grad_phi[_j][_qp](0);
272 35030016 : phi(1, 2) = phi(2, 1) = _grad_phi[_j][_qp](1);
273 35030016 : break;
274 : }
275 :
276 114005248 : jacobian += (_Jacobian_mult[_qp] * phi).trace() *
277 114005248 : (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
278 : }
279 :
280 434652346 : if (_ndisp != 3 && _out_of_plane_strain_coupled && _use_displaced_mesh)
281 : {
282 24064 : const Real out_of_plane_thickness = std::exp((*_out_of_plane_strain)[_qp]);
283 24064 : jacobian *= out_of_plane_thickness;
284 : }
285 :
286 : return jacobian;
287 : }
288 :
289 : Real
290 507580960 : StressDivergenceTensors::computeQpOffDiagJacobian(unsigned int jvar)
291 : {
292 : // off-diagonal Jacobian with respect to a coupled displacement component
293 1027008856 : for (unsigned int coupled_component = 0; coupled_component < _ndisp; ++coupled_component)
294 1015145560 : if (jvar == _disp_var[coupled_component])
295 : {
296 495717664 : if (_out_of_plane_direction != 2)
297 : {
298 4096 : if (coupled_component == _out_of_plane_direction)
299 0 : continue;
300 : }
301 :
302 495717664 : if (_use_finite_deform_jacobian)
303 8166464 : return ElasticityTensorTools::elasticJacobian(_finite_deform_Jacobian_mult[_qp],
304 8166464 : _component,
305 : coupled_component,
306 8166464 : _grad_test[_i][_qp],
307 495717664 : _grad_phi_undisplaced[_j][_qp]);
308 :
309 487551200 : const Real sum_C3x3 = _Jacobian_mult[_qp].sum3x3();
310 487551200 : const RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum3x1();
311 : Real jacobian = 0.0;
312 :
313 : // B^T_i * C * B_j
314 975102400 : jacobian += ElasticityTensorTools::elasticJacobian(_Jacobian_mult[_qp],
315 487551200 : _component,
316 : coupled_component,
317 487551200 : _grad_test[_i][_qp],
318 487551200 : _grad_phi[_j][_qp]);
319 :
320 487551200 : if (_volumetric_locking_correction)
321 : {
322 : // jacobian = Bbar^T_i * C * Bbar_j where Bbar = B + Bvol
323 : // jacobian = B^T_i * C * B_j + Bvol^T_i * C * Bvol_j + Bvol^T_i * C * B_j + B^T_i * C *
324 : // Bvol_j
325 :
326 : // Bvol^T_i * C * Bvol_j
327 98606080 : jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
328 98606080 : (_avg_grad_phi[_j][coupled_component] - _grad_phi[_j][_qp](coupled_component)) /
329 : 9.0;
330 :
331 : // B^T_i * C * Bvol_j
332 98606080 : jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
333 98606080 : (_avg_grad_phi[_j][coupled_component] - _grad_phi[_j][_qp](coupled_component)) /
334 : 3.0;
335 :
336 : // Bvol^T_i * C * B_i
337 98606080 : RankTwoTensor phi;
338 394424320 : for (unsigned int i = 0; i < 3; ++i)
339 295818240 : phi(coupled_component, i) = _grad_phi[_j][_qp](i);
340 :
341 98606080 : jacobian += (_Jacobian_mult[_qp] * phi).trace() *
342 98606080 : (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
343 : }
344 :
345 487551200 : return jacobian;
346 : }
347 :
348 : // off-diagonal Jacobian with respect to a coupled out_of_plane_strain variable
349 11863296 : if (_out_of_plane_strain_coupled && jvar == _out_of_plane_strain_var)
350 1024 : return _Jacobian_mult[_qp](
351 1024 : _component, _component, _out_of_plane_direction, _out_of_plane_direction) *
352 1024 : _grad_test[_i][_qp](_component) * _phi[_j][_qp];
353 :
354 : // bail out if jvar is not coupled
355 11862272 : if (getJvarMap()[jvar] < 0)
356 : return 0.0;
357 :
358 : // off-diagonal Jacobian with respect to any other coupled variable
359 : const unsigned int cvar = mapJvarToCvar(jvar);
360 1936384 : RankTwoTensor total_deigenstrain;
361 1939456 : for (const auto deigenstrain_darg : _deigenstrain_dargs[cvar])
362 3072 : total_deigenstrain += (*deigenstrain_darg)[_qp];
363 :
364 1936384 : return -((_Jacobian_mult[_qp] * total_deigenstrain) *
365 1936384 : _grad_test[_i][_qp])(_component)*_phi[_j][_qp];
366 : }
367 :
368 : void
369 303984 : StressDivergenceTensors::computeFiniteDeformJacobian()
370 : {
371 : usingTensorIndices(i_, j_, k_, l_);
372 : const auto I = RankTwoTensor::Identity();
373 303984 : const RankFourTensor I2 = I.times<i_, k_, j_, l_>(I);
374 :
375 : // Bring back to unrotated config
376 : const RankTwoTensor unrotated_stress =
377 303984 : (*_rotation_increment)[_qp].transpose() * _stress[_qp] * (*_rotation_increment)[_qp];
378 :
379 : // Incremental deformation gradient Fhat
380 : const RankTwoTensor Fhat =
381 303984 : (*_deformation_gradient)[_qp] * (*_deformation_gradient_old)[_qp].inverse();
382 303984 : const RankTwoTensor Fhatinv = Fhat.inverse();
383 :
384 303984 : const RankTwoTensor rot_times_stress = (*_rotation_increment)[_qp] * unrotated_stress;
385 : const RankFourTensor dstress_drot =
386 303984 : I.times<i_, k_, j_, l_>(rot_times_stress) + I.times<j_, k_, i_, l_>(rot_times_stress);
387 : const RankFourTensor rot_rank_four =
388 303984 : (*_rotation_increment)[_qp].times<i_, k_, j_, l_>((*_rotation_increment)[_qp]);
389 303984 : const RankFourTensor drot_dUhatinv = Fhat.times<i_, k_, j_, l_>(I);
390 :
391 303984 : const RankTwoTensor A = I - Fhatinv;
392 :
393 : // Ctilde = Chat^-1 - I
394 303984 : const RankTwoTensor Ctilde = A * A.transpose() - A - A.transpose();
395 : const RankFourTensor dCtilde_dFhatinv =
396 303984 : -I.times<i_, k_, j_, l_>(A) - I.times<j_, k_, i_, l_>(A) + I2 + I.times<j_, k_, i_, l_>(I);
397 :
398 : // Second order approximation of Uhat - consistent with strain increment definition
399 : // const RankTwoTensor Uhat = I - 0.5 * Ctilde - 3.0/8.0 * Ctilde * Ctilde;
400 :
401 : RankFourTensor dUhatinv_dCtilde =
402 607968 : 0.5 * I2 - 1.0 / 8.0 * (I.times<i_, k_, j_, l_>(Ctilde) + Ctilde.times<i_, k_, j_, l_>(I));
403 303984 : RankFourTensor drot_dFhatinv = drot_dUhatinv * dUhatinv_dCtilde * dCtilde_dFhatinv;
404 :
405 303984 : drot_dFhatinv -= Fhat.times<i_, k_, j_, l_>((*_rotation_increment)[_qp].transpose());
406 303984 : _finite_deform_Jacobian_mult[_qp] = dstress_drot * drot_dFhatinv;
407 :
408 : const RankFourTensor dstrain_increment_dCtilde =
409 607968 : -0.5 * I2 + 0.25 * (I.times<i_, k_, j_, l_>(Ctilde) + Ctilde.times<i_, k_, j_, l_>(I));
410 303984 : _finite_deform_Jacobian_mult[_qp] +=
411 303984 : rot_rank_four * _Jacobian_mult[_qp] * dstrain_increment_dCtilde * dCtilde_dFhatinv;
412 303984 : _finite_deform_Jacobian_mult[_qp] += Fhat.times<j_, k_, i_, l_>(_stress[_qp]);
413 :
414 303984 : const RankFourTensor dFhat_dFhatinv = -Fhat.times<i_, k_, j_, l_>(Fhat.transpose());
415 303984 : const RankTwoTensor dJ_dFhatinv = dFhat_dFhatinv.innerProductTranspose(Fhat.ddet());
416 :
417 : // Component from Jacobian derivative
418 303984 : _finite_deform_Jacobian_mult[_qp] += _stress[_qp].times<i_, j_, k_, l_>(dJ_dFhatinv);
419 :
420 : // Derivative of Fhatinv w.r.t. undisplaced coordinates
421 303984 : const RankTwoTensor Finv = (*_deformation_gradient)[_qp].inverse();
422 303984 : const RankFourTensor dFhatinv_dGradu = -Fhatinv.times<i_, k_, j_, l_>(Finv.transpose());
423 303984 : _finite_deform_Jacobian_mult[_qp] = _finite_deform_Jacobian_mult[_qp] * dFhatinv_dGradu;
424 303984 : }
425 :
426 : void
427 5334094 : StressDivergenceTensors::computeAverageGradientTest()
428 : {
429 : // Calculate volume averaged value of shape function derivative
430 5334094 : _avg_grad_test.resize(_test.size());
431 39714014 : for (_i = 0; _i < _test.size(); ++_i)
432 : {
433 34379920 : _avg_grad_test[_i].resize(3);
434 34379920 : _avg_grad_test[_i][_component] = 0.0;
435 276247952 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
436 241868032 : _avg_grad_test[_i][_component] += _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
437 :
438 34379920 : _avg_grad_test[_i][_component] /= _current_elem_volume;
439 : }
440 5334094 : }
441 :
442 : void
443 733518 : StressDivergenceTensors::computeAverageGradientPhi()
444 : {
445 : // Calculate volume average derivatives for phi
446 733518 : _avg_grad_phi.resize(_phi.size());
447 5618862 : for (_i = 0; _i < _phi.size(); ++_i)
448 : {
449 4885344 : _avg_grad_phi[_i].resize(3);
450 18558576 : for (unsigned int component = 0; component < _mesh.dimension(); ++component)
451 : {
452 13673232 : _avg_grad_phi[_i][component] = 0.0;
453 115196688 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
454 101523456 : _avg_grad_phi[_i][component] += _grad_phi[_i][_qp](component) * _JxW[_qp] * _coord[_qp];
455 :
456 13673232 : _avg_grad_phi[_i][component] /= _current_elem_volume;
457 : }
458 : }
459 733518 : }
|
