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