Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
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("TensorMechanicsApp", StressDivergenceRZTensors);
16 :
17 : InputParameters
18 476 : StressDivergenceRZTensors::validParams()
19 : {
20 476 : InputParameters params = StressDivergenceTensors::validParams();
21 476 : params.addClassDescription(
22 : "Calculate stress divergence for an axisymmetric problem in cylindrical coordinates.");
23 952 : 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 476 : params.set<bool>("use_displaced_mesh") = true;
29 476 : return params;
30 0 : }
31 :
32 238 : StressDivergenceRZTensors::StressDivergenceRZTensors(const InputParameters & parameters)
33 238 : : StressDivergenceTensors(parameters)
34 : {
35 238 : }
36 :
37 : void
38 220 : StressDivergenceRZTensors::initialSetup()
39 : {
40 : // check if any of the eigenstrains provide derivatives wrt variables that are not coupled
41 494 : for (auto eigenstrain_name : getParam<std::vector<MaterialPropertyName>>("eigenstrain_names"))
42 162 : validateNonlinearCoupling<RankTwoTensor>(eigenstrain_name);
43 :
44 220 : 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 220 : if (getBlockCoordSystem() == Moose::COORD_RZ && _fe_problem.getAxisymmetricRadialCoord() != 0)
49 1 : mooseError("rz_coord_axis=Y is the only supported option for StressDivergenceRZTensors");
50 219 : }
51 :
52 : Real
53 6574524 : StressDivergenceRZTensors::computeQpResidual()
54 : {
55 : Real div = 0.0;
56 6574524 : if (_component == 0)
57 : {
58 3292446 : div = _grad_test[_i][_qp](0) * _stress[_qp](0, 0) +
59 3292446 : +(_test[_i][_qp] / _q_point[_qp](0)) * _stress[_qp](2, 2) +
60 3292446 : +_grad_test[_i][_qp](1) * _stress[_qp](0, 1); // stress_{rz}
61 :
62 : // volumetric locking correction
63 3292446 : if (_volumetric_locking_correction)
64 719168 : div += (_avg_grad_test[_i][0] - _grad_test[_i][_qp](0) - _test[_i][_qp] / _q_point[_qp](0)) *
65 359584 : (_stress[_qp].trace()) / 3.0;
66 : }
67 3282078 : else if (_component == 1)
68 : {
69 3282078 : div = _grad_test[_i][_qp](1) * _stress[_qp](1, 1) +
70 3282078 : +_grad_test[_i][_qp](0) * _stress[_qp](1, 0); // stress_{zr}
71 :
72 : // volumetric locking correction
73 3282078 : if (_volumetric_locking_correction)
74 359584 : 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 6574524 : return div;
80 : }
81 :
82 : Real
83 6433068 : StressDivergenceRZTensors::computeQpJacobian()
84 : {
85 6433068 : return calculateJacobian(_component, _component);
86 : }
87 :
88 : Real
89 1659296 : StressDivergenceRZTensors::computeQpOffDiagJacobian(unsigned int jvar)
90 : {
91 2542896 : for (unsigned int i = 0; i < _ndisp; ++i)
92 2506928 : if (jvar == _disp_var[i])
93 1623328 : return calculateJacobian(_component, i);
94 :
95 : // bail out if jvar is not coupled
96 35968 : 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 23296 : RankTwoTensor total_deigenstrain;
102 46592 : for (const auto deigenstrain_darg : _deigenstrain_dargs[cvar])
103 23296 : total_deigenstrain += (*deigenstrain_darg)[_qp];
104 :
105 : Real jac = 0.0;
106 23296 : if (_component == 0)
107 : {
108 46592 : for (unsigned k = 0; k < LIBMESH_DIM; ++k)
109 139776 : for (unsigned l = 0; l < LIBMESH_DIM; ++l)
110 104832 : jac -= (_grad_test[_i][_qp](0) * _Jacobian_mult[_qp](0, 0, k, l) +
111 104832 : _test[_i][_qp] / _q_point[_qp](0) * _Jacobian_mult[_qp](2, 2, k, l) +
112 104832 : _grad_test[_i][_qp](1) * _Jacobian_mult[_qp](0, 1, k, l)) *
113 104832 : total_deigenstrain(k, l);
114 11648 : return jac * _phi[_j][_qp];
115 : }
116 11648 : else if (_component == 1)
117 : {
118 46592 : for (unsigned k = 0; k < LIBMESH_DIM; ++k)
119 139776 : for (unsigned l = 0; l < LIBMESH_DIM; ++l)
120 104832 : jac -= (_grad_test[_i][_qp](1) * _Jacobian_mult[_qp](1, 1, k, l) +
121 104832 : _grad_test[_i][_qp](0) * _Jacobian_mult[_qp](1, 0, k, l)) *
122 104832 : total_deigenstrain(k, l);
123 11648 : return jac * _phi[_j][_qp];
124 : }
125 :
126 : return 0.0;
127 : }
128 :
129 : Real
130 8056396 : 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 8056396 : if (ivar == 0) // Case grad_test for x, requires contributions from stress_xx, stress_xy, and
136 : // stress_zz
137 : {
138 4030470 : test(0) = _grad_test[_i][_qp](0);
139 4030470 : test(1) = _grad_test[_i][_qp](1);
140 4030470 : test_z(2) = _test[_i][_qp] / _q_point[_qp](0);
141 : }
142 : else // Case grad_test for y
143 : {
144 4025926 : test(0) = _grad_test[_i][_qp](0);
145 4025926 : test(1) = _grad_test[_i][_qp](1);
146 : }
147 :
148 8056396 : if (jvar == 0)
149 : {
150 4030470 : phi(0) = _grad_phi[_j][_qp](0);
151 4030470 : phi(1) = _grad_phi[_j][_qp](1);
152 4030470 : phi_z(2) = _phi[_j][_qp] / _q_point[_qp](0);
153 : }
154 : else
155 : {
156 4025926 : phi(0) = _grad_phi[_j][_qp](0);
157 4025926 : phi(1) = _grad_phi[_j][_qp](1);
158 : }
159 :
160 8056396 : if (ivar == 0 &&
161 : jvar == 0) // Case when both phi and test are functions of x and z; requires four terms
162 : {
163 3218806 : const Real first_sum = ElasticityTensorTools::elasticJacobian(
164 3218806 : _Jacobian_mult[_qp], ivar, jvar, test, phi); // test_x and phi_x
165 3218806 : const Real second_sum = ElasticityTensorTools::elasticJacobian(
166 3218806 : _Jacobian_mult[_qp], 2, 2, test_z, phi_z); // test_z and phi_z
167 3218806 : const Real mixed_sum1 = ElasticityTensorTools::elasticJacobian(
168 3218806 : _Jacobian_mult[_qp], ivar, 2, test, phi_z); // test_x and phi_z
169 3218806 : const Real mixed_sum2 = ElasticityTensorTools::elasticJacobian(
170 3218806 : _Jacobian_mult[_qp], 2, jvar, test_z, phi); // test_z and phi_x
171 :
172 3218806 : first_term = first_sum + second_sum + mixed_sum1 + mixed_sum2;
173 : }
174 4837590 : else if (ivar == 0 && jvar == 1)
175 : {
176 811664 : const Real first_sum = ElasticityTensorTools::elasticJacobian(
177 811664 : _Jacobian_mult[_qp], ivar, jvar, test, phi); // test_x and phi_y
178 811664 : const Real mixed_sum2 = ElasticityTensorTools::elasticJacobian(
179 811664 : _Jacobian_mult[_qp], 2, jvar, test_z, phi); // test_z and phi_y
180 :
181 811664 : first_term = first_sum + mixed_sum2;
182 : }
183 4025926 : else if (ivar == 1 && jvar == 0)
184 : {
185 811664 : const Real second_sum = ElasticityTensorTools::elasticJacobian(
186 811664 : _Jacobian_mult[_qp], ivar, jvar, test, phi); // test_y and phi_x
187 811664 : const Real mixed_sum1 = ElasticityTensorTools::elasticJacobian(
188 811664 : _Jacobian_mult[_qp], ivar, 2, test, phi_z); // test_y and phi_z
189 :
190 811664 : first_term = second_sum + mixed_sum1;
191 : }
192 3214262 : else if (ivar == 1 && jvar == 1)
193 3214262 : first_term = ElasticityTensorTools::elasticJacobian(
194 3214262 : _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 8056396 : if (_volumetric_locking_correction)
203 : {
204 : RealGradient new_test(2, 0.0);
205 : RealGradient new_phi(2, 0.0);
206 :
207 2545152 : new_test(0) = _grad_test[_i][_qp](0) + _test[_i][_qp] / _q_point[_qp](0);
208 2545152 : new_test(1) = _grad_test[_i][_qp](1);
209 2545152 : new_phi(0) = _grad_phi[_j][_qp](0) + _phi[_j][_qp] / _q_point[_qp](0);
210 2545152 : new_phi(1) = _grad_phi[_j][_qp](1);
211 :
212 : // Bvol^T_i * C * Bvol_j
213 2545152 : val += _Jacobian_mult[_qp].sum3x3() * (_avg_grad_test[_i][ivar] - new_test(ivar)) *
214 2545152 : (_avg_grad_phi[_j][jvar] - new_phi(jvar)) / 3.0;
215 :
216 : // B^T_i * C * Bvol_j
217 2545152 : RealGradient sum_3x1 = _Jacobian_mult[_qp].sum3x1();
218 2545152 : if (ivar == 0 && jvar == 0)
219 684672 : val += (sum_3x1(0) * test(0) + sum_3x1(2) * test_z(2)) * (_avg_grad_phi[_j][0] - new_phi(0));
220 1860480 : else if (ivar == 0 && jvar == 1)
221 587904 : val += (sum_3x1(0) * test(0) + sum_3x1(2) * test_z(2)) * (_avg_grad_phi[_j][1] - new_phi(1));
222 1272576 : else if (ivar == 1 && jvar == 0)
223 587904 : val += sum_3x1(1) * test(1) * (_avg_grad_phi[_j][0] - new_phi(0));
224 : else
225 684672 : 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 2545152 : if (jvar == 0)
230 5090304 : for (unsigned int i = 0; i < 3; ++i)
231 3817728 : val +=
232 3817728 : (_Jacobian_mult[_qp](i, i, 0, 0) * phi(0) + _Jacobian_mult[_qp](i, i, 0, 1) * phi(1) +
233 3817728 : _Jacobian_mult[_qp](i, i, 2, 2) * phi_z(2)) *
234 3817728 : (_avg_grad_test[_i][ivar] - new_test(ivar));
235 1272576 : else if (jvar == 1)
236 5090304 : for (unsigned int i = 0; i < 3; ++i)
237 3817728 : val +=
238 3817728 : (_Jacobian_mult[_qp](i, i, 0, 1) * phi(0) + _Jacobian_mult[_qp](i, i, 1, 1) * phi(1)) *
239 3817728 : (_avg_grad_test[_i][ivar] - new_test(ivar));
240 : }
241 :
242 8056396 : return val / 3.0 + first_term;
243 : }
244 :
245 : void
246 103088 : StressDivergenceRZTensors::computeAverageGradientTest()
247 : {
248 : // calculate volume averaged value of shape function derivative
249 103088 : _avg_grad_test.resize(_test.size());
250 515440 : for (_i = 0; _i < _test.size(); ++_i)
251 : {
252 412352 : _avg_grad_test[_i].resize(2);
253 412352 : _avg_grad_test[_i][_component] = 0.0;
254 2061760 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
255 : {
256 1649408 : if (_component == 0)
257 824704 : _avg_grad_test[_i][_component] +=
258 824704 : (_grad_test[_i][_qp](_component) + _test[_i][_qp] / _q_point[_qp](0)) * _JxW[_qp] *
259 824704 : _coord[_qp];
260 : else
261 824704 : _avg_grad_test[_i][_component] += _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
262 : }
263 412352 : _avg_grad_test[_i][_component] /= _current_elem_volume;
264 : }
265 103088 : }
266 :
267 : void
268 58140 : StressDivergenceRZTensors::computeAverageGradientPhi()
269 : {
270 58140 : _avg_grad_phi.resize(_phi.size());
271 290700 : for (_i = 0; _i < _phi.size(); ++_i)
272 : {
273 232560 : _avg_grad_phi[_i].resize(2);
274 697680 : for (unsigned int component = 0; component < 2; ++component)
275 : {
276 465120 : _avg_grad_phi[_i][component] = 0.0;
277 2325600 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
278 : {
279 1860480 : if (component == 0)
280 930240 : _avg_grad_phi[_i][component] +=
281 930240 : (_grad_phi[_i][_qp](component) + _phi[_i][_qp] / _q_point[_qp](0)) * _JxW[_qp] *
282 930240 : _coord[_qp];
283 : else
284 930240 : _avg_grad_phi[_i][component] += _grad_phi[_i][_qp](component) * _JxW[_qp] * _coord[_qp];
285 : }
286 465120 : _avg_grad_phi[_i][component] /= _current_elem_volume;
287 : }
288 : }
289 58140 : }
|
