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