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 "StressDivergenceRZTensors.h"
11 : #include "Assembly.h"
12 : #include "ElasticityTensorTools.h"
13 : #include "libmesh/quadrature.h"
14 :
15 : registerMooseObject("SolidMechanicsApp", StressDivergenceRZTensors);
16 :
17 : InputParameters
18 594 : StressDivergenceRZTensors::validParams()
19 : {
20 594 : InputParameters params = StressDivergenceTensors::validParams();
21 594 : params.addClassDescription(
22 : "Calculate stress divergence for an axisymmetric problem in cylindrical coordinates.");
23 1188 : params.addRequiredRangeCheckedParam<unsigned int>(
24 : "component",
25 : "component < 2",
26 : "An integer corresponding to the direction the variable this kernel acts in. (0 "
27 : "refers to the radial and 1 to the axial displacement.)");
28 594 : params.set<bool>("use_displaced_mesh") = true;
29 594 : return params;
30 0 : }
31 :
32 297 : StressDivergenceRZTensors::StressDivergenceRZTensors(const InputParameters & parameters)
33 297 : : StressDivergenceTensors(parameters)
34 : {
35 297 : }
36 :
37 : void
38 273 : StressDivergenceRZTensors::initialSetup()
39 : {
40 : // check if any of the eigenstrains provide derivatives wrt variables that are not coupled
41 616 : for (auto eigenstrain_name : getParam<std::vector<MaterialPropertyName>>("eigenstrain_names"))
42 210 : validateNonlinearCoupling<RankTwoTensor>(eigenstrain_name);
43 :
44 273 : if (getBlockCoordSystem() != Moose::COORD_RZ)
45 0 : mooseError("The coordinate system in the Problem block must be set to RZ for axisymmetric "
46 : "geometries.");
47 :
48 273 : if (getBlockCoordSystem() == Moose::COORD_RZ && _fe_problem.getAxisymmetricRadialCoord() != 0)
49 1 : mooseError("rz_coord_axis=Y is the only supported option for StressDivergenceRZTensors");
50 272 : }
51 :
52 : Real
53 11192556 : StressDivergenceRZTensors::computeQpResidual()
54 : {
55 : Real div = 0.0;
56 11192556 : if (_component == 0)
57 : {
58 5603718 : div = _grad_test[_i][_qp](0) * _stress[_qp](0, 0) +
59 5603718 : +(_test[_i][_qp] / _q_point[_qp](0)) * _stress[_qp](2, 2) +
60 5603718 : +_grad_test[_i][_qp](1) * _stress[_qp](0, 1); // stress_{rz}
61 :
62 : // volumetric locking correction
63 5603718 : if (_volumetric_locking_correction)
64 1058784 : div += (_avg_grad_test[_i][0] - _grad_test[_i][_qp](0) - _test[_i][_qp] / _q_point[_qp](0)) *
65 529392 : (_stress[_qp].trace()) / 3.0;
66 : }
67 5588838 : else if (_component == 1)
68 : {
69 5588838 : div = _grad_test[_i][_qp](1) * _stress[_qp](1, 1) +
70 5588838 : +_grad_test[_i][_qp](0) * _stress[_qp](1, 0); // stress_{zr}
71 :
72 : // volumetric locking correction
73 5588838 : if (_volumetric_locking_correction)
74 529392 : div += (_avg_grad_test[_i][1] - _grad_test[_i][_qp](1)) * (_stress[_qp].trace()) / 3.0;
75 : }
76 : else
77 0 : mooseError("Invalid component for this AxisymmetricRZ problem.");
78 :
79 11192556 : return div;
80 : }
81 :
82 : Real
83 11637080 : StressDivergenceRZTensors::computeQpJacobian()
84 : {
85 11637080 : return calculateJacobian(_component, _component);
86 : }
87 :
88 : Real
89 3021376 : StressDivergenceRZTensors::computeQpOffDiagJacobian(unsigned int jvar)
90 : {
91 4624032 : for (unsigned int i = 0; i < _ndisp; ++i)
92 4562720 : if (jvar == _disp_var[i])
93 2960064 : return calculateJacobian(_component, i);
94 :
95 : // bail out if jvar is not coupled
96 61312 : if (getJvarMap()[jvar] < 0)
97 : return 0.0;
98 :
99 : // off-diagonal Jacobian with respect to any other coupled variable
100 : const unsigned int cvar = mapJvarToCvar(jvar);
101 35968 : RankTwoTensor total_deigenstrain;
102 71936 : for (const auto deigenstrain_darg : _deigenstrain_dargs[cvar])
103 35968 : total_deigenstrain += (*deigenstrain_darg)[_qp];
104 :
105 : Real jac = 0.0;
106 35968 : if (_component == 0)
107 : {
108 71936 : for (unsigned k = 0; k < LIBMESH_DIM; ++k)
109 215808 : for (unsigned l = 0; l < LIBMESH_DIM; ++l)
110 161856 : jac -= (_grad_test[_i][_qp](0) * _Jacobian_mult[_qp](0, 0, k, l) +
111 161856 : _test[_i][_qp] / _q_point[_qp](0) * _Jacobian_mult[_qp](2, 2, k, l) +
112 161856 : _grad_test[_i][_qp](1) * _Jacobian_mult[_qp](0, 1, k, l)) *
113 161856 : total_deigenstrain(k, l);
114 17984 : return jac * _phi[_j][_qp];
115 : }
116 17984 : else if (_component == 1)
117 : {
118 71936 : for (unsigned k = 0; k < LIBMESH_DIM; ++k)
119 215808 : for (unsigned l = 0; l < LIBMESH_DIM; ++l)
120 161856 : jac -= (_grad_test[_i][_qp](1) * _Jacobian_mult[_qp](1, 1, k, l) +
121 161856 : _grad_test[_i][_qp](0) * _Jacobian_mult[_qp](1, 0, k, l)) *
122 161856 : total_deigenstrain(k, l);
123 17984 : return jac * _phi[_j][_qp];
124 : }
125 :
126 : return 0.0;
127 : }
128 :
129 : Real
130 14597144 : StressDivergenceRZTensors::calculateJacobian(unsigned int ivar, unsigned int jvar)
131 : {
132 : // B^T_i * C * B_j
133 : RealGradient test, test_z, phi, phi_z;
134 : Real first_term = 0.0;
135 14597144 : if (ivar == 0) // Case grad_test for x, requires contributions from stress_xx, stress_xy, and
136 : // stress_zz
137 : {
138 7302860 : test(0) = _grad_test[_i][_qp](0);
139 7302860 : test(1) = _grad_test[_i][_qp](1);
140 7302860 : test_z(2) = _test[_i][_qp] / _q_point[_qp](0);
141 : }
142 : else // Case grad_test for y
143 : {
144 7294284 : test(0) = _grad_test[_i][_qp](0);
145 7294284 : test(1) = _grad_test[_i][_qp](1);
146 : }
147 :
148 14597144 : if (jvar == 0)
149 : {
150 7302860 : phi(0) = _grad_phi[_j][_qp](0);
151 7302860 : phi(1) = _grad_phi[_j][_qp](1);
152 7302860 : phi_z(2) = _phi[_j][_qp] / _q_point[_qp](0);
153 : }
154 : else
155 : {
156 7294284 : phi(0) = _grad_phi[_j][_qp](0);
157 7294284 : phi(1) = _grad_phi[_j][_qp](1);
158 : }
159 :
160 14597144 : if (ivar == 0 &&
161 : jvar == 0) // Case when both phi and test are functions of x and z; requires four terms
162 : {
163 5822828 : const Real first_sum = ElasticityTensorTools::elasticJacobian(
164 5822828 : _Jacobian_mult[_qp], ivar, jvar, test, phi); // test_x and phi_x
165 5822828 : const Real second_sum = ElasticityTensorTools::elasticJacobian(
166 5822828 : _Jacobian_mult[_qp], 2, 2, test_z, phi_z); // test_z and phi_z
167 5822828 : const Real mixed_sum1 = ElasticityTensorTools::elasticJacobian(
168 5822828 : _Jacobian_mult[_qp], ivar, 2, test, phi_z); // test_x and phi_z
169 5822828 : const Real mixed_sum2 = ElasticityTensorTools::elasticJacobian(
170 5822828 : _Jacobian_mult[_qp], 2, jvar, test_z, phi); // test_z and phi_x
171 :
172 5822828 : first_term = first_sum + second_sum + mixed_sum1 + mixed_sum2;
173 : }
174 8774316 : else if (ivar == 0 && jvar == 1)
175 : {
176 1480032 : const Real first_sum = ElasticityTensorTools::elasticJacobian(
177 1480032 : _Jacobian_mult[_qp], ivar, jvar, test, phi); // test_x and phi_y
178 1480032 : const Real mixed_sum2 = ElasticityTensorTools::elasticJacobian(
179 1480032 : _Jacobian_mult[_qp], 2, jvar, test_z, phi); // test_z and phi_y
180 :
181 1480032 : first_term = first_sum + mixed_sum2;
182 : }
183 7294284 : else if (ivar == 1 && jvar == 0)
184 : {
185 1480032 : const Real second_sum = ElasticityTensorTools::elasticJacobian(
186 1480032 : _Jacobian_mult[_qp], ivar, jvar, test, phi); // test_y and phi_x
187 1480032 : const Real mixed_sum1 = ElasticityTensorTools::elasticJacobian(
188 1480032 : _Jacobian_mult[_qp], ivar, 2, test, phi_z); // test_y and phi_z
189 :
190 1480032 : first_term = second_sum + mixed_sum1;
191 : }
192 5814252 : else if (ivar == 1 && jvar == 1)
193 5814252 : first_term = ElasticityTensorTools::elasticJacobian(
194 5814252 : _Jacobian_mult[_qp], ivar, jvar, test, phi); // test_y and phi_y
195 : else
196 0 : mooseError("Invalid component in Jacobian Calculation");
197 :
198 : Real val = 0.0;
199 : // volumetric locking correction
200 : // K = Bbar^T_i * C * Bbar^T_j where Bbar = B + Bvol
201 : // K = B^T_i * C * B_j + Bvol^T_i * C * Bvol_j + B^T_i * C * Bvol_j + Bvol^T_i * C * B_j
202 14597144 : if (_volumetric_locking_correction)
203 : {
204 : RealGradient new_test(2, 0.0);
205 : RealGradient new_phi(2, 0.0);
206 :
207 5001856 : new_test(0) = _grad_test[_i][_qp](0) + _test[_i][_qp] / _q_point[_qp](0);
208 5001856 : new_test(1) = _grad_test[_i][_qp](1);
209 5001856 : new_phi(0) = _grad_phi[_j][_qp](0) + _phi[_j][_qp] / _q_point[_qp](0);
210 5001856 : new_phi(1) = _grad_phi[_j][_qp](1);
211 :
212 : // Bvol^T_i * C * Bvol_j
213 5001856 : val += _Jacobian_mult[_qp].sum3x3() * (_avg_grad_test[_i][ivar] - new_test(ivar)) *
214 5001856 : (_avg_grad_phi[_j][jvar] - new_phi(jvar)) / 3.0;
215 :
216 : // B^T_i * C * Bvol_j
217 5001856 : RealGradient sum_3x1 = _Jacobian_mult[_qp].sum3x1();
218 5001856 : if (ivar == 0 && jvar == 0)
219 1325504 : val += (sum_3x1(0) * test(0) + sum_3x1(2) * test_z(2)) * (_avg_grad_phi[_j][0] - new_phi(0));
220 3676352 : else if (ivar == 0 && jvar == 1)
221 1175424 : val += (sum_3x1(0) * test(0) + sum_3x1(2) * test_z(2)) * (_avg_grad_phi[_j][1] - new_phi(1));
222 2500928 : else if (ivar == 1 && jvar == 0)
223 1175424 : val += sum_3x1(1) * test(1) * (_avg_grad_phi[_j][0] - new_phi(0));
224 : else
225 1325504 : val += sum_3x1(1) * test(1) * (_avg_grad_phi[_j][1] - new_phi(1));
226 :
227 : // Bvol^T_i * C * B_j
228 : // val = trace (C * B_j) *(avg_grad_test[_i][ivar] - new_test(ivar))
229 5001856 : if (jvar == 0)
230 10003712 : for (unsigned int i = 0; i < 3; ++i)
231 7502784 : val +=
232 7502784 : (_Jacobian_mult[_qp](i, i, 0, 0) * phi(0) + _Jacobian_mult[_qp](i, i, 0, 1) * phi(1) +
233 7502784 : _Jacobian_mult[_qp](i, i, 2, 2) * phi_z(2)) *
234 7502784 : (_avg_grad_test[_i][ivar] - new_test(ivar));
235 2500928 : else if (jvar == 1)
236 10003712 : for (unsigned int i = 0; i < 3; ++i)
237 7502784 : val +=
238 7502784 : (_Jacobian_mult[_qp](i, i, 0, 1) * phi(0) + _Jacobian_mult[_qp](i, i, 1, 1) * phi(1)) *
239 7502784 : (_avg_grad_test[_i][ivar] - new_test(ivar));
240 : }
241 :
242 14597144 : return val / 3.0 + first_term;
243 : }
244 :
245 : void
246 181060 : StressDivergenceRZTensors::computeAverageGradientTest()
247 : {
248 : // calculate volume averaged value of shape function derivative
249 181060 : _avg_grad_test.resize(_test.size());
250 905300 : for (_i = 0; _i < _test.size(); ++_i)
251 : {
252 724240 : _avg_grad_test[_i].resize(2);
253 724240 : _avg_grad_test[_i][_component] = 0.0;
254 3621200 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
255 : {
256 2896960 : if (_component == 0)
257 1448480 : _avg_grad_test[_i][_component] +=
258 1448480 : (_grad_test[_i][_qp](_component) + _test[_i][_qp] / _q_point[_qp](0)) * _JxW[_qp] *
259 1448480 : _coord[_qp];
260 : else
261 1448480 : _avg_grad_test[_i][_component] += _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
262 : }
263 724240 : _avg_grad_test[_i][_component] /= _current_elem_volume;
264 : }
265 181060 : }
266 :
267 : void
268 114886 : StressDivergenceRZTensors::computeAverageGradientPhi()
269 : {
270 114886 : _avg_grad_phi.resize(_phi.size());
271 574430 : for (_i = 0; _i < _phi.size(); ++_i)
272 : {
273 459544 : _avg_grad_phi[_i].resize(2);
274 1378632 : for (unsigned int component = 0; component < 2; ++component)
275 : {
276 919088 : _avg_grad_phi[_i][component] = 0.0;
277 4595440 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
278 : {
279 3676352 : if (component == 0)
280 1838176 : _avg_grad_phi[_i][component] +=
281 3676352 : (_grad_phi[_i][_qp](component) + _phi[_i][_qp] / _q_point[_qp](0)) * _JxW[_qp] *
282 1838176 : _coord[_qp];
283 : else
284 1838176 : _avg_grad_phi[_i][component] += _grad_phi[_i][_qp](component) * _JxW[_qp] * _coord[_qp];
285 : }
286 919088 : _avg_grad_phi[_i][component] /= _current_elem_volume;
287 : }
288 : }
289 114886 : }
|
