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 20764 : StressDivergenceTensors::validParams()
25 : {
26 20764 : InputParameters params = JvarMapKernelInterface<ALEKernel>::validParams();
27 20764 : params.addClassDescription("Stress divergence kernel for the Cartesian coordinate system");
28 41528 : 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 41528 : 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 41528 : 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 41528 : 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 41528 : params.addCoupledVar("out_of_plane_strain",
48 : "The name of the out_of_plane_strain variable used in the "
49 : "WeakPlaneStress kernel.");
50 41528 : MooseEnum out_of_plane_direction("x y z", "z");
51 41528 : 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 41528 : params.addParam<std::string>("base_name", "Material property base name");
56 20764 : params.set<bool>("use_displaced_mesh") = false;
57 41528 : params.addParam<bool>(
58 41528 : "use_finite_deform_jacobian", false, "Jacobian for corotational finite strain");
59 41528 : params.addParam<bool>("volumetric_locking_correction",
60 41528 : false,
61 : "Set to false to turn off volumetric locking correction");
62 20764 : return params;
63 20764 : }
64 :
65 10352 : StressDivergenceTensors::StressDivergenceTensors(const InputParameters & parameters)
66 : : JvarMapKernelInterface<ALEKernel>(parameters),
67 11252 : _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
68 20704 : _use_finite_deform_jacobian(getParam<bool>("use_finite_deform_jacobian")),
69 20704 : _stress(getMaterialPropertyByName<RankTwoTensor>(_base_name + "stress")),
70 31056 : _Jacobian_mult(getMaterialPropertyByName<RankFourTensor>(_base_name + "Jacobian_mult")),
71 20704 : _component(getParam<unsigned int>("component")),
72 10352 : _ndisp(coupledComponents("displacements")),
73 10352 : _disp_var(_ndisp),
74 10352 : _out_of_plane_strain_coupled(isCoupled("out_of_plane_strain")),
75 10352 : _out_of_plane_strain(_out_of_plane_strain_coupled ? &coupledValue("out_of_plane_strain")
76 : : nullptr),
77 10352 : _out_of_plane_strain_var(_out_of_plane_strain_coupled ? coupled("out_of_plane_strain") : 0),
78 20704 : _out_of_plane_direction(getParam<MooseEnum>("out_of_plane_direction")),
79 20704 : _use_displaced_mesh(getParam<bool>("use_displaced_mesh")),
80 10352 : _avg_grad_test(_test.size(), std::vector<Real>(3, 0.0)),
81 10352 : _avg_grad_phi(_phi.size(), std::vector<Real>(3, 0.0)),
82 31056 : _volumetric_locking_correction(getParam<bool>("volumetric_locking_correction"))
83 : {
84 : // get coupled displacements
85 39684 : for (unsigned int i = 0; i < _ndisp; ++i)
86 29332 : _disp_var[i] = coupled("displacements", i);
87 :
88 : // fetch eigenstrain derivatives
89 : const auto nvar = _coupled_moose_vars.size();
90 10352 : _deigenstrain_dargs.resize(nvar);
91 41317 : for (std::size_t i = 0; i < nvar; ++i)
92 63655 : for (auto eigenstrain_name : getParam<std::vector<MaterialPropertyName>>("eigenstrain_names"))
93 3450 : _deigenstrain_dargs[i].push_back(&getMaterialPropertyDerivative<RankTwoTensor>(
94 1725 : eigenstrain_name, _coupled_moose_vars[i]->name()));
95 :
96 : // Checking for consistency between mesh size and length of the provided displacements vector
97 10352 : 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 10352 : 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 10352 : 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 10352 : 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 10352 : 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 10352 : }
122 :
123 : void
124 9828 : StressDivergenceTensors::initialSetup()
125 : {
126 : // check if any of the eigenstrains provide derivatives wrt variables that are not coupled
127 20206 : for (auto eigenstrain_name : getParam<std::vector<MaterialPropertyName>>("eigenstrain_names"))
128 1650 : validateNonlinearCoupling<RankTwoTensor>(eigenstrain_name);
129 :
130 : // make sure the coordinate system is cartesioan
131 9828 : 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 9828 : }
135 :
136 : void
137 13205483 : StressDivergenceTensors::computeResidual()
138 : {
139 13205483 : prepareVectorTag(_assembly, _var.number());
140 :
141 13205483 : if (_volumetric_locking_correction)
142 7082830 : computeAverageGradientTest();
143 :
144 13205483 : precalculateResidual();
145 102040467 : for (_i = 0; _i < _test.size(); ++_i)
146 799481080 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
147 710646096 : _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
148 :
149 13205483 : accumulateTaggedLocalResidual();
150 :
151 13205483 : if (_has_save_in)
152 : {
153 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
154 96936 : for (const auto & var : _save_in)
155 48468 : var->sys().solution().add_vector(_local_re, var->dofIndices());
156 : }
157 13205483 : }
158 :
159 : Real
160 681628868 : StressDivergenceTensors::computeQpResidual()
161 : {
162 681628868 : Real residual = _stress[_qp].row(_component) * _grad_test[_i][_qp];
163 : // volumetric locking correction
164 681628868 : if (_volumetric_locking_correction)
165 294554656 : residual += _stress[_qp].trace() / 3.0 *
166 294554656 : (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component));
167 :
168 681628868 : if (_ndisp != 3 && _out_of_plane_strain_coupled && _use_displaced_mesh)
169 : {
170 68288 : const Real out_of_plane_thickness = std::exp((*_out_of_plane_strain)[_qp]);
171 68288 : residual *= out_of_plane_thickness;
172 : }
173 :
174 681628868 : return residual;
175 : }
176 :
177 : void
178 1685403 : StressDivergenceTensors::computeJacobian()
179 : {
180 1685403 : if (_volumetric_locking_correction)
181 : {
182 710308 : computeAverageGradientTest();
183 710308 : computeAverageGradientPhi();
184 : }
185 :
186 1685403 : if (_use_finite_deform_jacobian)
187 : {
188 22008 : _finite_deform_Jacobian_mult.resize(_qrule->n_points());
189 :
190 191480 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
191 169472 : computeFiniteDeformJacobian();
192 :
193 22008 : ALEKernel::computeJacobian();
194 : }
195 : else
196 1663395 : Kernel::computeJacobian();
197 1685403 : }
198 :
199 : void
200 2938448 : StressDivergenceTensors::computeOffDiagJacobian(const unsigned int jvar)
201 : {
202 2938448 : if (_volumetric_locking_correction)
203 : {
204 938824 : computeAverageGradientPhi();
205 938824 : computeAverageGradientTest();
206 : }
207 :
208 2938448 : 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 2893772 : Kernel::computeOffDiagJacobian(jvar);
219 2938448 : }
220 :
221 : Real
222 823452866 : StressDivergenceTensors::computeQpJacobian()
223 : {
224 823452866 : if (_use_finite_deform_jacobian)
225 10778432 : return ElasticityTensorTools::elasticJacobian(_finite_deform_Jacobian_mult[_qp],
226 : _component,
227 10778432 : _component,
228 10778432 : _grad_test[_i][_qp],
229 10778432 : _grad_phi_undisplaced[_j][_qp]);
230 :
231 812674434 : Real sum_C3x3 = _Jacobian_mult[_qp].sum3x3();
232 812674434 : RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum3x1();
233 :
234 : Real jacobian = 0.0;
235 : // B^T_i * C * B_j
236 812674434 : jacobian += ElasticityTensorTools::elasticJacobian(
237 812674434 : _Jacobian_mult[_qp], _component, _component, _grad_test[_i][_qp], _grad_phi[_j][_qp]);
238 :
239 812674434 : 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 220534272 : jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
247 220534272 : (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 9.0;
248 :
249 : // B^T_i * C * Bvol_j
250 220534272 : jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
251 220534272 : (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 3.0;
252 :
253 : // Bvol^T_i * C * B_j
254 220534272 : RankTwoTensor phi;
255 220534272 : switch (_component)
256 : {
257 76346624 : case 0:
258 76346624 : phi(0, 0) = _grad_phi[_j][_qp](0);
259 76346624 : phi(0, 1) = phi(1, 0) = _grad_phi[_j][_qp](1);
260 76346624 : phi(0, 2) = phi(2, 0) = _grad_phi[_j][_qp](2);
261 76346624 : break;
262 :
263 76346624 : case 1:
264 76346624 : phi(1, 1) = _grad_phi[_j][_qp](1);
265 76346624 : phi(0, 1) = phi(1, 0) = _grad_phi[_j][_qp](0);
266 76346624 : phi(1, 2) = phi(2, 1) = _grad_phi[_j][_qp](2);
267 76346624 : break;
268 :
269 67841024 : case 2:
270 67841024 : phi(2, 2) = _grad_phi[_j][_qp](2);
271 67841024 : phi(0, 2) = phi(2, 0) = _grad_phi[_j][_qp](0);
272 67841024 : phi(1, 2) = phi(2, 1) = _grad_phi[_j][_qp](1);
273 67841024 : break;
274 : }
275 :
276 220534272 : jacobian += (_Jacobian_mult[_qp] * phi).trace() *
277 220534272 : (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
278 : }
279 :
280 812674434 : 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 941310512 : StressDivergenceTensors::computeQpOffDiagJacobian(unsigned int jvar)
291 : {
292 : // off-diagonal Jacobian with respect to a coupled displacement component
293 1896368928 : for (unsigned int coupled_component = 0; coupled_component < _ndisp; ++coupled_component)
294 1877242784 : if (jvar == _disp_var[coupled_component])
295 : {
296 922184368 : if (_out_of_plane_direction != 2)
297 : {
298 4096 : if (coupled_component == _out_of_plane_direction)
299 0 : continue;
300 : }
301 :
302 922184368 : 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 922184368 : _grad_phi_undisplaced[_j][_qp]);
308 :
309 908211056 : const Real sum_C3x3 = _Jacobian_mult[_qp].sum3x3();
310 908211056 : const RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum3x1();
311 : Real jacobian = 0.0;
312 :
313 : // B^T_i * C * B_j
314 1816422112 : jacobian += ElasticityTensorTools::elasticJacobian(_Jacobian_mult[_qp],
315 908211056 : _component,
316 : coupled_component,
317 908211056 : _grad_test[_i][_qp],
318 908211056 : _grad_phi[_j][_qp]);
319 :
320 908211056 : 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 230210816 : jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
328 230210816 : (_avg_grad_phi[_j][coupled_component] - _grad_phi[_j][_qp](coupled_component)) /
329 : 9.0;
330 :
331 : // B^T_i * C * Bvol_j
332 230210816 : jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
333 230210816 : (_avg_grad_phi[_j][coupled_component] - _grad_phi[_j][_qp](coupled_component)) /
334 : 3.0;
335 :
336 : // Bvol^T_i * C * B_i
337 230210816 : RankTwoTensor phi;
338 920843264 : for (unsigned int i = 0; i < 3; ++i)
339 690632448 : phi(coupled_component, i) = _grad_phi[_j][_qp](i);
340 :
341 230210816 : jacobian += (_Jacobian_mult[_qp] * phi).trace() *
342 230210816 : (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
343 : }
344 :
345 908211056 : return jacobian;
346 : }
347 :
348 : // off-diagonal Jacobian with respect to a coupled out_of_plane_strain variable
349 19126144 : 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 19125120 : 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 513696 : StressDivergenceTensors::computeFiniteDeformJacobian()
370 : {
371 : usingTensorIndices(i_, j_, k_, l_);
372 : const auto I = RankTwoTensor::Identity();
373 513696 : const RankFourTensor I2 = I.times<i_, k_, j_, l_>(I);
374 :
375 : // Bring back to unrotated config
376 : const RankTwoTensor unrotated_stress =
377 513696 : (*_rotation_increment)[_qp].transpose() * _stress[_qp] * (*_rotation_increment)[_qp];
378 :
379 : // Incremental deformation gradient Fhat
380 : const RankTwoTensor Fhat =
381 513696 : (*_deformation_gradient)[_qp] * (*_deformation_gradient_old)[_qp].inverse();
382 513696 : const RankTwoTensor Fhatinv = Fhat.inverse();
383 :
384 513696 : const RankTwoTensor rot_times_stress = (*_rotation_increment)[_qp] * unrotated_stress;
385 : const RankFourTensor dstress_drot =
386 513696 : 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 513696 : (*_rotation_increment)[_qp].times<i_, k_, j_, l_>((*_rotation_increment)[_qp]);
389 513696 : const RankFourTensor drot_dUhatinv = Fhat.times<i_, k_, j_, l_>(I);
390 :
391 513696 : const RankTwoTensor A = I - Fhatinv;
392 :
393 : // Ctilde = Chat^-1 - I
394 513696 : const RankTwoTensor Ctilde = A * A.transpose() - A - A.transpose();
395 : const RankFourTensor dCtilde_dFhatinv =
396 513696 : -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 1027392 : 0.5 * I2 - 1.0 / 8.0 * (I.times<i_, k_, j_, l_>(Ctilde) + Ctilde.times<i_, k_, j_, l_>(I));
403 513696 : RankFourTensor drot_dFhatinv = drot_dUhatinv * dUhatinv_dCtilde * dCtilde_dFhatinv;
404 :
405 513696 : drot_dFhatinv -= Fhat.times<i_, k_, j_, l_>((*_rotation_increment)[_qp].transpose());
406 513696 : _finite_deform_Jacobian_mult[_qp] = dstress_drot * drot_dFhatinv;
407 :
408 : const RankFourTensor dstrain_increment_dCtilde =
409 1027392 : -0.5 * I2 + 0.25 * (I.times<i_, k_, j_, l_>(Ctilde) + Ctilde.times<i_, k_, j_, l_>(I));
410 513696 : _finite_deform_Jacobian_mult[_qp] +=
411 513696 : rot_rank_four * _Jacobian_mult[_qp] * dstrain_increment_dCtilde * dCtilde_dFhatinv;
412 513696 : _finite_deform_Jacobian_mult[_qp] += Fhat.times<j_, k_, i_, l_>(_stress[_qp]);
413 :
414 513696 : const RankFourTensor dFhat_dFhatinv = -Fhat.times<i_, k_, j_, l_>(Fhat.transpose());
415 513696 : const RankTwoTensor dJ_dFhatinv = dFhat_dFhatinv.innerProductTranspose(Fhat.ddet());
416 :
417 : // Component from Jacobian derivative
418 513696 : _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 513696 : const RankTwoTensor Finv = (*_deformation_gradient)[_qp].inverse();
422 513696 : const RankFourTensor dFhatinv_dGradu = -Fhatinv.times<i_, k_, j_, l_>(Finv.transpose());
423 513696 : _finite_deform_Jacobian_mult[_qp] = _finite_deform_Jacobian_mult[_qp] * dFhatinv_dGradu;
424 513696 : }
425 :
426 : void
427 8550902 : StressDivergenceTensors::computeAverageGradientTest()
428 : {
429 : // Calculate volume averaged value of shape function derivative
430 8550902 : _avg_grad_test.resize(_test.size());
431 62176318 : for (_i = 0; _i < _test.size(); ++_i)
432 : {
433 53625416 : _avg_grad_test[_i].resize(3);
434 53625416 : _avg_grad_test[_i][_component] = 0.0;
435 423501544 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
436 369876128 : _avg_grad_test[_i][_component] += _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
437 :
438 53625416 : _avg_grad_test[_i][_component] /= _current_elem_volume;
439 : }
440 8550902 : }
441 :
442 : void
443 1534246 : StressDivergenceTensors::computeAverageGradientPhi()
444 : {
445 : // Calculate volume average derivatives for phi
446 1534246 : _avg_grad_phi.resize(_phi.size());
447 11902358 : for (_i = 0; _i < _phi.size(); ++_i)
448 : {
449 10368112 : _avg_grad_phi[_i].resize(3);
450 39566592 : for (unsigned int component = 0; component < _mesh.dimension(); ++component)
451 : {
452 29198480 : _avg_grad_phi[_i][component] = 0.0;
453 247539472 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
454 218340992 : _avg_grad_phi[_i][component] += _grad_phi[_i][_qp](component) * _JxW[_qp] * _coord[_qp];
455 :
456 29198480 : _avg_grad_phi[_i][component] /= _current_elem_volume;
457 : }
458 : }
459 1534246 : }
|
