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 1070 : StressDivergenceRZTensors::validParams()
19 : {
20 1070 : InputParameters params = StressDivergenceTensors::validParams();
21 1070 : params.addClassDescription(
22 : "Calculate stress divergence for an axisymmetric problem in cylindrical coordinates.");
23 2140 : 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 1070 : params.set<bool>("use_displaced_mesh") = true;
29 1070 : return params;
30 0 : }
31 :
32 535 : StressDivergenceRZTensors::StressDivergenceRZTensors(const InputParameters & parameters)
33 535 : : StressDivergenceTensors(parameters)
34 : {
35 535 : }
36 :
37 : void
38 493 : StressDivergenceRZTensors::initialSetup()
39 : {
40 : // check if any of the eigenstrains provide derivatives wrt variables that are not coupled
41 1110 : for (auto eigenstrain_name : getParam<std::vector<MaterialPropertyName>>("eigenstrain_names"))
42 372 : validateNonlinearCoupling<RankTwoTensor>(eigenstrain_name);
43 :
44 493 : 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 493 : if (getBlockCoordSystem() == Moose::COORD_RZ && _fe_problem.getAxisymmetricRadialCoord() != 0)
49 2 : mooseError("rz_coord_axis=Y is the only supported option for StressDivergenceRZTensors");
50 491 : }
51 :
52 : Real
53 17322228 : StressDivergenceRZTensors::computeQpResidual()
54 : {
55 : Real div = 0.0;
56 17322228 : if (_component == 0)
57 : {
58 8672850 : div = _grad_test[_i][_qp](0) * _stress[_qp](0, 0) +
59 8672850 : +(_test[_i][_qp] / _q_point[_qp](0)) * _stress[_qp](2, 2) +
60 8672850 : +_grad_test[_i][_qp](1) * _stress[_qp](0, 1); // stress_{rz}
61 :
62 : // volumetric locking correction
63 8672850 : if (_volumetric_locking_correction)
64 1614688 : div += (_avg_grad_test[_i][0] - _grad_test[_i][_qp](0) - _test[_i][_qp] / _q_point[_qp](0)) *
65 807344 : (_stress[_qp].trace()) / 3.0;
66 : }
67 8649378 : else if (_component == 1)
68 : {
69 8649378 : div = _grad_test[_i][_qp](1) * _stress[_qp](1, 1) +
70 8649378 : +_grad_test[_i][_qp](0) * _stress[_qp](1, 0); // stress_{zr}
71 :
72 : // volumetric locking correction
73 8649378 : if (_volumetric_locking_correction)
74 807344 : 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 17322228 : return div;
80 : }
81 :
82 : Real
83 18070148 : StressDivergenceRZTensors::computeQpJacobian()
84 : {
85 18070148 : return calculateJacobian(_component, _component);
86 : }
87 :
88 : Real
89 4680672 : StressDivergenceRZTensors::computeQpOffDiagJacobian(unsigned int jvar)
90 : {
91 7166928 : for (unsigned int i = 0; i < _ndisp; ++i)
92 7069648 : if (jvar == _disp_var[i])
93 4583392 : return calculateJacobian(_component, i);
94 :
95 : // bail out if jvar is not coupled
96 97280 : 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 59264 : RankTwoTensor total_deigenstrain;
102 118528 : for (const auto deigenstrain_darg : _deigenstrain_dargs[cvar])
103 59264 : total_deigenstrain += (*deigenstrain_darg)[_qp];
104 :
105 : Real jac = 0.0;
106 59264 : if (_component == 0)
107 : {
108 118528 : for (unsigned k = 0; k < LIBMESH_DIM; ++k)
109 355584 : for (unsigned l = 0; l < LIBMESH_DIM; ++l)
110 266688 : jac -= (_grad_test[_i][_qp](0) * _Jacobian_mult[_qp](0, 0, k, l) +
111 266688 : _test[_i][_qp] / _q_point[_qp](0) * _Jacobian_mult[_qp](2, 2, k, l) +
112 266688 : _grad_test[_i][_qp](1) * _Jacobian_mult[_qp](0, 1, k, l)) *
113 266688 : total_deigenstrain(k, l);
114 29632 : return jac * _phi[_j][_qp];
115 : }
116 29632 : else if (_component == 1)
117 : {
118 118528 : for (unsigned k = 0; k < LIBMESH_DIM; ++k)
119 355584 : for (unsigned l = 0; l < LIBMESH_DIM; ++l)
120 266688 : jac -= (_grad_test[_i][_qp](1) * _Jacobian_mult[_qp](1, 1, k, l) +
121 266688 : _grad_test[_i][_qp](0) * _Jacobian_mult[_qp](1, 0, k, l)) *
122 266688 : total_deigenstrain(k, l);
123 29632 : return jac * _phi[_j][_qp];
124 : }
125 :
126 : return 0.0;
127 : }
128 :
129 : Real
130 22653540 : 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 22653540 : if (ivar == 0) // Case grad_test for x, requires contributions from stress_xx, stress_xy, and
136 : // stress_zz
137 : {
138 11333330 : test(0) = _grad_test[_i][_qp](0);
139 11333330 : test(1) = _grad_test[_i][_qp](1);
140 11333330 : test_z(2) = _test[_i][_qp] / _q_point[_qp](0);
141 : }
142 : else // Case grad_test for y
143 : {
144 11320210 : test(0) = _grad_test[_i][_qp](0);
145 11320210 : test(1) = _grad_test[_i][_qp](1);
146 : }
147 :
148 22653540 : if (jvar == 0)
149 : {
150 11333330 : phi(0) = _grad_phi[_j][_qp](0);
151 11333330 : phi(1) = _grad_phi[_j][_qp](1);
152 11333330 : phi_z(2) = _phi[_j][_qp] / _q_point[_qp](0);
153 : }
154 : else
155 : {
156 11320210 : phi(0) = _grad_phi[_j][_qp](0);
157 11320210 : phi(1) = _grad_phi[_j][_qp](1);
158 : }
159 :
160 22653540 : if (ivar == 0 &&
161 : jvar == 0) // Case when both phi and test are functions of x and z; requires four terms
162 : {
163 9041634 : const Real first_sum = ElasticityTensorTools::elasticJacobian(
164 9041634 : _Jacobian_mult[_qp], ivar, jvar, test, phi); // test_x and phi_x
165 9041634 : const Real second_sum = ElasticityTensorTools::elasticJacobian(
166 9041634 : _Jacobian_mult[_qp], 2, 2, test_z, phi_z); // test_z and phi_z
167 9041634 : const Real mixed_sum1 = ElasticityTensorTools::elasticJacobian(
168 9041634 : _Jacobian_mult[_qp], ivar, 2, test, phi_z); // test_x and phi_z
169 9041634 : const Real mixed_sum2 = ElasticityTensorTools::elasticJacobian(
170 9041634 : _Jacobian_mult[_qp], 2, jvar, test_z, phi); // test_z and phi_x
171 :
172 9041634 : first_term = first_sum + second_sum + mixed_sum1 + mixed_sum2;
173 : }
174 13611906 : else if (ivar == 0 && jvar == 1)
175 : {
176 2291696 : const Real first_sum = ElasticityTensorTools::elasticJacobian(
177 2291696 : _Jacobian_mult[_qp], ivar, jvar, test, phi); // test_x and phi_y
178 2291696 : const Real mixed_sum2 = ElasticityTensorTools::elasticJacobian(
179 2291696 : _Jacobian_mult[_qp], 2, jvar, test_z, phi); // test_z and phi_y
180 :
181 2291696 : first_term = first_sum + mixed_sum2;
182 : }
183 11320210 : else if (ivar == 1 && jvar == 0)
184 : {
185 2291696 : const Real second_sum = ElasticityTensorTools::elasticJacobian(
186 2291696 : _Jacobian_mult[_qp], ivar, jvar, test, phi); // test_y and phi_x
187 2291696 : const Real mixed_sum1 = ElasticityTensorTools::elasticJacobian(
188 2291696 : _Jacobian_mult[_qp], ivar, 2, test, phi_z); // test_y and phi_z
189 :
190 2291696 : first_term = second_sum + mixed_sum1;
191 : }
192 9028514 : else if (ivar == 1 && jvar == 1)
193 9028514 : first_term = ElasticityTensorTools::elasticJacobian(
194 9028514 : _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 22653540 : if (_volumetric_locking_correction)
203 : {
204 : RealGradient new_test(2, 0.0);
205 : RealGradient new_phi(2, 0.0);
206 :
207 7547008 : new_test(0) = _grad_test[_i][_qp](0) + _test[_i][_qp] / _q_point[_qp](0);
208 7547008 : new_test(1) = _grad_test[_i][_qp](1);
209 7547008 : new_phi(0) = _grad_phi[_j][_qp](0) + _phi[_j][_qp] / _q_point[_qp](0);
210 7547008 : new_phi(1) = _grad_phi[_j][_qp](1);
211 :
212 : // Bvol^T_i * C * Bvol_j
213 7547008 : val += _Jacobian_mult[_qp].sum3x3() * (_avg_grad_test[_i][ivar] - new_test(ivar)) *
214 7547008 : (_avg_grad_phi[_j][jvar] - new_phi(jvar)) / 3.0;
215 :
216 : // B^T_i * C * Bvol_j
217 7547008 : RealGradient sum_3x1 = _Jacobian_mult[_qp].sum3x1();
218 7547008 : if (ivar == 0 && jvar == 0)
219 2010176 : val += (sum_3x1(0) * test(0) + sum_3x1(2) * test_z(2)) * (_avg_grad_phi[_j][0] - new_phi(0));
220 5536832 : else if (ivar == 0 && jvar == 1)
221 1763328 : val += (sum_3x1(0) * test(0) + sum_3x1(2) * test_z(2)) * (_avg_grad_phi[_j][1] - new_phi(1));
222 3773504 : else if (ivar == 1 && jvar == 0)
223 1763328 : val += sum_3x1(1) * test(1) * (_avg_grad_phi[_j][0] - new_phi(0));
224 : else
225 2010176 : 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 7547008 : if (jvar == 0)
230 15094016 : for (unsigned int i = 0; i < 3; ++i)
231 11320512 : val +=
232 11320512 : (_Jacobian_mult[_qp](i, i, 0, 0) * phi(0) + _Jacobian_mult[_qp](i, i, 0, 1) * phi(1) +
233 11320512 : _Jacobian_mult[_qp](i, i, 2, 2) * phi_z(2)) *
234 11320512 : (_avg_grad_test[_i][ivar] - new_test(ivar));
235 3773504 : else if (jvar == 1)
236 15094016 : for (unsigned int i = 0; i < 3; ++i)
237 11320512 : val +=
238 11320512 : (_Jacobian_mult[_qp](i, i, 0, 1) * phi(0) + _Jacobian_mult[_qp](i, i, 1, 1) * phi(1)) *
239 11320512 : (_avg_grad_test[_i][ivar] - new_test(ivar));
240 : }
241 :
242 22653540 : return val / 3.0 + first_term;
243 : }
244 :
245 : void
246 273944 : StressDivergenceRZTensors::computeAverageGradientTest()
247 : {
248 : // calculate volume averaged value of shape function derivative
249 273944 : _avg_grad_test.resize(_test.size());
250 1369720 : for (_i = 0; _i < _test.size(); ++_i)
251 : {
252 1095776 : _avg_grad_test[_i].resize(2);
253 1095776 : _avg_grad_test[_i][_component] = 0.0;
254 5478880 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
255 : {
256 4383104 : if (_component == 0)
257 2191552 : _avg_grad_test[_i][_component] +=
258 2191552 : (_grad_test[_i][_qp](_component) + _test[_i][_qp] / _q_point[_qp](0)) * _JxW[_qp] *
259 2191552 : _coord[_qp];
260 : else
261 2191552 : _avg_grad_test[_i][_component] += _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
262 : }
263 1095776 : _avg_grad_test[_i][_component] /= _current_elem_volume;
264 : }
265 273944 : }
266 :
267 : void
268 173026 : StressDivergenceRZTensors::computeAverageGradientPhi()
269 : {
270 173026 : _avg_grad_phi.resize(_phi.size());
271 865130 : for (_i = 0; _i < _phi.size(); ++_i)
272 : {
273 692104 : _avg_grad_phi[_i].resize(2);
274 2076312 : for (unsigned int component = 0; component < 2; ++component)
275 : {
276 1384208 : _avg_grad_phi[_i][component] = 0.0;
277 6921040 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
278 : {
279 5536832 : if (component == 0)
280 2768416 : _avg_grad_phi[_i][component] +=
281 5536832 : (_grad_phi[_i][_qp](component) + _phi[_i][_qp] / _q_point[_qp](0)) * _JxW[_qp] *
282 2768416 : _coord[_qp];
283 : else
284 2768416 : _avg_grad_phi[_i][component] += _grad_phi[_i][_qp](component) * _JxW[_qp] * _coord[_qp];
285 : }
286 1384208 : _avg_grad_phi[_i][component] /= _current_elem_volume;
287 : }
288 : }
289 173026 : }
|
