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 "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("SolidMechanicsApp", StressDivergenceTensors);
22 :
23 : InputParameters
24 20908 : StressDivergenceTensors::validParams()
25 : {
26 20908 : InputParameters params = JvarMapKernelInterface<ALEKernel>::validParams();
27 20908 : params.addClassDescription("Stress divergence kernel for the Cartesian coordinate system");
28 41816 : 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 41816 : 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 41816 : 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 41816 : 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 41816 : params.addCoupledVar("out_of_plane_strain",
48 : "The name of the out_of_plane_strain variable used in the "
49 : "WeakPlaneStress kernel.");
50 41816 : MooseEnum out_of_plane_direction("x y z", "z");
51 41816 : 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 41816 : params.addParam<std::string>("base_name", "Material property base name");
56 20908 : params.set<bool>("use_displaced_mesh") = false;
57 41816 : params.addParam<bool>(
58 41816 : "use_finite_deform_jacobian", false, "Jacobian for corotational finite strain");
59 41816 : params.addParam<bool>("volumetric_locking_correction",
60 41816 : false,
61 : "Set to false to turn off volumetric locking correction");
62 20908 : return params;
63 20908 : }
64 :
65 10424 : StressDivergenceTensors::StressDivergenceTensors(const InputParameters & parameters)
66 : : JvarMapKernelInterface<ALEKernel>(parameters),
67 11324 : _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
68 20848 : _use_finite_deform_jacobian(getParam<bool>("use_finite_deform_jacobian")),
69 20848 : _stress(getMaterialPropertyByName<RankTwoTensor>(_base_name + "stress")),
70 31272 : _Jacobian_mult(getMaterialPropertyByName<RankFourTensor>(_base_name + "Jacobian_mult")),
71 20848 : _component(getParam<unsigned int>("component")),
72 10424 : _ndisp(coupledComponents("displacements")),
73 10424 : _disp_var(_ndisp),
74 10424 : _out_of_plane_strain_coupled(isCoupled("out_of_plane_strain")),
75 10424 : _out_of_plane_strain(_out_of_plane_strain_coupled ? &coupledValue("out_of_plane_strain")
76 : : nullptr),
77 10424 : _out_of_plane_strain_var(_out_of_plane_strain_coupled ? coupled("out_of_plane_strain") : 0),
78 20848 : _out_of_plane_direction(getParam<MooseEnum>("out_of_plane_direction")),
79 20848 : _use_displaced_mesh(getParam<bool>("use_displaced_mesh")),
80 10424 : _avg_grad_test(_test.size(), std::vector<Real>(3, 0.0)),
81 10424 : _avg_grad_phi(_phi.size(), std::vector<Real>(3, 0.0)),
82 31272 : _volumetric_locking_correction(getParam<bool>("volumetric_locking_correction"))
83 : {
84 : // get coupled displacements
85 39936 : for (unsigned int i = 0; i < _ndisp; ++i)
86 29512 : _disp_var[i] = coupled("displacements", i);
87 :
88 : // fetch eigenstrain derivatives
89 : const auto nvar = _coupled_moose_vars.size();
90 10424 : _deigenstrain_dargs.resize(nvar);
91 41619 : for (std::size_t i = 0; i < nvar; ++i)
92 64259 : for (auto eigenstrain_name : getParam<std::vector<MaterialPropertyName>>("eigenstrain_names"))
93 3738 : _deigenstrain_dargs[i].push_back(&getMaterialPropertyDerivative<RankTwoTensor>(
94 1869 : eigenstrain_name, _coupled_moose_vars[i]->name()));
95 :
96 : // Checking for consistency between mesh size and length of the provided displacements vector
97 10424 : 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 10424 : 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 10424 : if (_use_finite_deform_jacobian)
104 : {
105 176 : _deformation_gradient =
106 176 : &getMaterialProperty<RankTwoTensor>(_base_name + "deformation_gradient");
107 176 : _deformation_gradient_old =
108 176 : &getMaterialPropertyOld<RankTwoTensor>(_base_name + "deformation_gradient");
109 352 : _rotation_increment = &getMaterialProperty<RankTwoTensor>(_base_name + "rotation_increment");
110 : }
111 :
112 : // Error if volumetric locking correction is turned on for 1D problems
113 10424 : 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 10424 : 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 10424 : }
122 :
123 : void
124 9894 : StressDivergenceTensors::initialSetup()
125 : {
126 : // check if any of the eigenstrains provide derivatives wrt variables that are not coupled
127 20374 : for (auto eigenstrain_name : getParam<std::vector<MaterialPropertyName>>("eigenstrain_names"))
128 1758 : validateNonlinearCoupling<RankTwoTensor>(eigenstrain_name);
129 :
130 : // make sure the coordinate system is cartesioan
131 9894 : 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 9894 : }
135 :
136 : void
137 13303022 : StressDivergenceTensors::computeResidual()
138 : {
139 13303022 : prepareVectorTag(_assembly, _var.number());
140 :
141 13303022 : if (_volumetric_locking_correction)
142 7175221 : computeAverageGradientTest();
143 :
144 13303022 : precalculateResidual();
145 103014606 : for (_i = 0; _i < _test.size(); ++_i)
146 807755632 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
147 718044048 : _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
148 :
149 13303022 : accumulateTaggedLocalResidual();
150 :
151 13303022 : if (_has_save_in)
152 : {
153 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
154 398504 : for (const auto & var : _save_in)
155 199252 : var->sys().solution().add_vector(_local_re, var->dofIndices());
156 : }
157 13303022 : }
158 :
159 : Real
160 689027108 : StressDivergenceTensors::computeQpResidual()
161 : {
162 689027108 : Real residual = _stress[_qp].row(_component) * _grad_test[_i][_qp];
163 : // volumetric locking correction
164 689027108 : if (_volumetric_locking_correction)
165 301837408 : residual += _stress[_qp].trace() / 3.0 *
166 301837408 : (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component));
167 :
168 689027108 : if (_ndisp != 3 && _out_of_plane_strain_coupled && _use_displaced_mesh)
169 : {
170 68576 : const Real out_of_plane_thickness = std::exp((*_out_of_plane_strain)[_qp]);
171 68576 : residual *= out_of_plane_thickness;
172 : }
173 :
174 689027108 : return residual;
175 : }
176 :
177 : void
178 1681168 : StressDivergenceTensors::computeJacobian()
179 : {
180 1681168 : if (_volumetric_locking_correction)
181 : {
182 709504 : computeAverageGradientTest();
183 709504 : computeAverageGradientPhi();
184 : }
185 :
186 1681168 : if (_use_finite_deform_jacobian)
187 : {
188 21966 : _finite_deform_Jacobian_mult.resize(_qrule->n_points());
189 :
190 191102 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
191 169136 : computeFiniteDeformJacobian();
192 :
193 21966 : ALEKernel::computeJacobian();
194 : }
195 : else
196 1659202 : Kernel::computeJacobian();
197 1681168 : }
198 :
199 : void
200 2904079 : StressDivergenceTensors::computeOffDiagJacobian(const unsigned int jvar)
201 : {
202 2904079 : if (_volumetric_locking_correction)
203 : {
204 931048 : computeAverageGradientPhi();
205 931048 : computeAverageGradientTest();
206 : }
207 :
208 2904079 : if (_use_finite_deform_jacobian)
209 : {
210 44676 : _finite_deform_Jacobian_mult.resize(_qrule->n_points());
211 :
212 388900 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
213 344224 : computeFiniteDeformJacobian();
214 :
215 44676 : ALEKernel::computeOffDiagJacobian(jvar);
216 : }
217 : else
218 2859403 : Kernel::computeOffDiagJacobian(jvar);
219 2904079 : }
220 :
221 : Real
222 826105026 : StressDivergenceTensors::computeQpJacobian()
223 : {
224 826105026 : if (_use_finite_deform_jacobian)
225 10756928 : return ElasticityTensorTools::elasticJacobian(_finite_deform_Jacobian_mult[_qp],
226 : _component,
227 10756928 : _component,
228 10756928 : _grad_test[_i][_qp],
229 10756928 : _grad_phi_undisplaced[_j][_qp]);
230 :
231 815348098 : Real sum_C3x3 = _Jacobian_mult[_qp].sum3x3();
232 815348098 : RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum3x1();
233 :
234 : Real jacobian = 0.0;
235 : // B^T_i * C * B_j
236 815348098 : jacobian += ElasticityTensorTools::elasticJacobian(
237 815348098 : _Jacobian_mult[_qp], _component, _component, _grad_test[_i][_qp], _grad_phi[_j][_qp]);
238 :
239 815348098 : 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 223337472 : jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
247 223337472 : (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 9.0;
248 :
249 : // B^T_i * C * Bvol_j
250 223337472 : jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
251 223337472 : (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 3.0;
252 :
253 : // Bvol^T_i * C * B_j
254 223337472 : RankTwoTensor phi;
255 223337472 : switch (_component)
256 : {
257 77204992 : case 0:
258 77204992 : phi(0, 0) = _grad_phi[_j][_qp](0);
259 77204992 : phi(0, 1) = phi(1, 0) = _grad_phi[_j][_qp](1);
260 77204992 : phi(0, 2) = phi(2, 0) = _grad_phi[_j][_qp](2);
261 77204992 : break;
262 :
263 77204992 : case 1:
264 77204992 : phi(1, 1) = _grad_phi[_j][_qp](1);
265 77204992 : phi(0, 1) = phi(1, 0) = _grad_phi[_j][_qp](0);
266 77204992 : phi(1, 2) = phi(2, 1) = _grad_phi[_j][_qp](2);
267 77204992 : break;
268 :
269 68927488 : case 2:
270 68927488 : phi(2, 2) = _grad_phi[_j][_qp](2);
271 68927488 : phi(0, 2) = phi(2, 0) = _grad_phi[_j][_qp](0);
272 68927488 : phi(1, 2) = phi(2, 1) = _grad_phi[_j][_qp](1);
273 68927488 : break;
274 : }
275 :
276 223337472 : jacobian += (_Jacobian_mult[_qp] * phi).trace() *
277 223337472 : (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
278 : }
279 :
280 815348098 : if (_ndisp != 3 && _out_of_plane_strain_coupled && _use_displaced_mesh)
281 : {
282 54528 : const Real out_of_plane_thickness = std::exp((*_out_of_plane_strain)[_qp]);
283 54528 : jacobian *= out_of_plane_thickness;
284 : }
285 :
286 : return jacobian;
287 : }
288 :
289 : Real
290 937346352 : StressDivergenceTensors::computeQpOffDiagJacobian(unsigned int jvar)
291 : {
292 : // off-diagonal Jacobian with respect to a coupled displacement component
293 1887398304 : for (unsigned int coupled_component = 0; coupled_component < _ndisp; ++coupled_component)
294 1869451808 : if (jvar == _disp_var[coupled_component])
295 : {
296 919399856 : if (_out_of_plane_direction != 2)
297 : {
298 4096 : if (coupled_component == _out_of_plane_direction)
299 0 : continue;
300 : }
301 :
302 919399856 : if (_use_finite_deform_jacobian)
303 13973312 : return ElasticityTensorTools::elasticJacobian(_finite_deform_Jacobian_mult[_qp],
304 13973312 : _component,
305 : coupled_component,
306 13973312 : _grad_test[_i][_qp],
307 919399856 : _grad_phi_undisplaced[_j][_qp]);
308 :
309 905426544 : const Real sum_C3x3 = _Jacobian_mult[_qp].sum3x3();
310 905426544 : const RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum3x1();
311 : Real jacobian = 0.0;
312 :
313 : // B^T_i * C * B_j
314 1810853088 : jacobian += ElasticityTensorTools::elasticJacobian(_Jacobian_mult[_qp],
315 905426544 : _component,
316 : coupled_component,
317 905426544 : _grad_test[_i][_qp],
318 905426544 : _grad_phi[_j][_qp]);
319 :
320 905426544 : 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 227556608 : jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
328 227556608 : (_avg_grad_phi[_j][coupled_component] - _grad_phi[_j][_qp](coupled_component)) /
329 : 9.0;
330 :
331 : // B^T_i * C * Bvol_j
332 227556608 : jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
333 227556608 : (_avg_grad_phi[_j][coupled_component] - _grad_phi[_j][_qp](coupled_component)) /
334 : 3.0;
335 :
336 : // Bvol^T_i * C * B_i
337 227556608 : RankTwoTensor phi;
338 910226432 : for (unsigned int i = 0; i < 3; ++i)
339 682669824 : phi(coupled_component, i) = _grad_phi[_j][_qp](i);
340 :
341 227556608 : jacobian += (_Jacobian_mult[_qp] * phi).trace() *
342 227556608 : (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
343 : }
344 :
345 905426544 : return jacobian;
346 : }
347 :
348 : // off-diagonal Jacobian with respect to a coupled out_of_plane_strain variable
349 17946496 : 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 17945472 : 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 1438720 : RankTwoTensor total_deigenstrain;
361 1441792 : for (const auto deigenstrain_darg : _deigenstrain_dargs[cvar])
362 3072 : total_deigenstrain += (*deigenstrain_darg)[_qp];
363 :
364 1438720 : return -((_Jacobian_mult[_qp] * total_deigenstrain) *
365 1438720 : _grad_test[_i][_qp])(_component)*_phi[_j][_qp];
366 : }
367 :
368 : void
369 513360 : StressDivergenceTensors::computeFiniteDeformJacobian()
370 : {
371 : usingTensorIndices(i_, j_, k_, l_);
372 : const auto I = RankTwoTensor::Identity();
373 513360 : const RankFourTensor I2 = I.times<i_, k_, j_, l_>(I);
374 :
375 : // Bring back to unrotated config
376 : const RankTwoTensor unrotated_stress =
377 513360 : (*_rotation_increment)[_qp].transpose() * _stress[_qp] * (*_rotation_increment)[_qp];
378 :
379 : // Incremental deformation gradient Fhat
380 : const RankTwoTensor Fhat =
381 513360 : (*_deformation_gradient)[_qp] * (*_deformation_gradient_old)[_qp].inverse();
382 513360 : const RankTwoTensor Fhatinv = Fhat.inverse();
383 :
384 513360 : const RankTwoTensor rot_times_stress = (*_rotation_increment)[_qp] * unrotated_stress;
385 : const RankFourTensor dstress_drot =
386 513360 : 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 513360 : (*_rotation_increment)[_qp].times<i_, k_, j_, l_>((*_rotation_increment)[_qp]);
389 513360 : const RankFourTensor drot_dUhatinv = Fhat.times<i_, k_, j_, l_>(I);
390 :
391 513360 : const RankTwoTensor A = I - Fhatinv;
392 :
393 : // Ctilde = Chat^-1 - I
394 513360 : const RankTwoTensor Ctilde = A * A.transpose() - A - A.transpose();
395 : const RankFourTensor dCtilde_dFhatinv =
396 513360 : -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 1026720 : 0.5 * I2 - 1.0 / 8.0 * (I.times<i_, k_, j_, l_>(Ctilde) + Ctilde.times<i_, k_, j_, l_>(I));
403 513360 : RankFourTensor drot_dFhatinv = drot_dUhatinv * dUhatinv_dCtilde * dCtilde_dFhatinv;
404 :
405 513360 : drot_dFhatinv -= Fhat.times<i_, k_, j_, l_>((*_rotation_increment)[_qp].transpose());
406 513360 : _finite_deform_Jacobian_mult[_qp] = dstress_drot * drot_dFhatinv;
407 :
408 : const RankFourTensor dstrain_increment_dCtilde =
409 1026720 : -0.5 * I2 + 0.25 * (I.times<i_, k_, j_, l_>(Ctilde) + Ctilde.times<i_, k_, j_, l_>(I));
410 513360 : _finite_deform_Jacobian_mult[_qp] +=
411 513360 : rot_rank_four * _Jacobian_mult[_qp] * dstrain_increment_dCtilde * dCtilde_dFhatinv;
412 513360 : _finite_deform_Jacobian_mult[_qp] += Fhat.times<j_, k_, i_, l_>(_stress[_qp]);
413 :
414 513360 : const RankFourTensor dFhat_dFhatinv = -Fhat.times<i_, k_, j_, l_>(Fhat.transpose());
415 513360 : const RankTwoTensor dJ_dFhatinv = dFhat_dFhatinv.innerProductTranspose(Fhat.ddet());
416 :
417 : // Component from Jacobian derivative
418 513360 : _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 513360 : const RankTwoTensor Finv = (*_deformation_gradient)[_qp].inverse();
422 513360 : const RankFourTensor dFhatinv_dGradu = -Fhatinv.times<i_, k_, j_, l_>(Finv.transpose());
423 513360 : _finite_deform_Jacobian_mult[_qp] = _finite_deform_Jacobian_mult[_qp] * dFhatinv_dGradu;
424 513360 : }
425 :
426 : void
427 8634713 : StressDivergenceTensors::computeAverageGradientTest()
428 : {
429 : // Calculate volume averaged value of shape function derivative
430 8634713 : _avg_grad_test.resize(_test.size());
431 63073273 : for (_i = 0; _i < _test.size(); ++_i)
432 : {
433 54438560 : _avg_grad_test[_i].resize(3);
434 54438560 : _avg_grad_test[_i][_component] = 0.0;
435 431390464 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
436 376951904 : _avg_grad_test[_i][_component] += _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
437 :
438 54438560 : _avg_grad_test[_i][_component] /= _current_elem_volume;
439 : }
440 8634713 : }
441 :
442 : void
443 1525666 : StressDivergenceTensors::computeAverageGradientPhi()
444 : {
445 : // Calculate volume average derivatives for phi
446 1525666 : _avg_grad_phi.resize(_phi.size());
447 11853650 : for (_i = 0; _i < _phi.size(); ++_i)
448 : {
449 10327984 : _avg_grad_phi[_i].resize(3);
450 39434592 : for (unsigned int component = 0; component < _mesh.dimension(); ++component)
451 : {
452 29106608 : _avg_grad_phi[_i][component] = 0.0;
453 246940720 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
454 217834112 : _avg_grad_phi[_i][component] += _grad_phi[_i][_qp](component) * _JxW[_qp] * _coord[_qp];
455 :
456 29106608 : _avg_grad_phi[_i][component] /= _current_elem_volume;
457 : }
458 : }
459 1525666 : }
|
