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