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 "StressDivergenceBeam.h"
11 :
12 : // MOOSE includes
13 : #include "Assembly.h"
14 : #include "Material.h"
15 : #include "MooseVariable.h"
16 : #include "SystemBase.h"
17 : #include "RankTwoTensor.h"
18 : #include "NonlinearSystem.h"
19 : #include "MooseMesh.h"
20 :
21 : #include "libmesh/quadrature.h"
22 :
23 : registerMooseObject("TensorMechanicsApp", StressDivergenceBeam);
24 :
25 : InputParameters
26 1630 : StressDivergenceBeam::validParams()
27 : {
28 1630 : InputParameters params = Kernel::validParams();
29 1630 : params.addClassDescription("Quasi-static and dynamic stress divergence kernel for Beam element");
30 3260 : params.addRequiredParam<unsigned int>(
31 : "component",
32 : "An integer corresponding to the direction "
33 : "the variable this kernel acts in. (0 for disp_x, "
34 : "1 for disp_y, 2 for disp_z, 3 for rot_x, 4 for rot_y and 5 for rot_z)");
35 3260 : params.addRequiredCoupledVar(
36 : "displacements",
37 : "The displacements appropriate for the simulation geometry and coordinate system");
38 3260 : params.addRequiredCoupledVar(
39 : "rotations", "The rotations appropriate for the simulation geometry and coordinate system");
40 3260 : params.addParam<MaterialPropertyName>(
41 : "zeta",
42 3260 : 0.0,
43 : "Name of material property or a constant real number defining the zeta parameter for the "
44 : "Rayleigh damping.");
45 4890 : params.addRangeCheckedParam<Real>(
46 3260 : "alpha", 0.0, "alpha >= -0.3333 & alpha <= 0.0", "alpha parameter for HHT time integration");
47 :
48 1630 : params.set<bool>("use_displaced_mesh") = true;
49 1630 : return params;
50 0 : }
51 :
52 979 : StressDivergenceBeam::StressDivergenceBeam(const InputParameters & parameters)
53 : : Kernel(parameters),
54 979 : _component(getParam<unsigned int>("component")),
55 979 : _ndisp(coupledComponents("displacements")),
56 979 : _disp_var(_ndisp),
57 979 : _nrot(coupledComponents("rotations")),
58 979 : _rot_var(_nrot),
59 1958 : _force(getMaterialPropertyByName<RealVectorValue>("forces")),
60 979 : _moment(getMaterialPropertyByName<RealVectorValue>("moments")),
61 979 : _K11(getMaterialPropertyByName<RankTwoTensor>("Jacobian_11")),
62 979 : _K22(getMaterialPropertyByName<RankTwoTensor>("Jacobian_22")),
63 979 : _K22_cross(getMaterialPropertyByName<RankTwoTensor>("Jacobian_22_cross")),
64 979 : _K21_cross(getMaterialPropertyByName<RankTwoTensor>("Jacobian_12")),
65 979 : _K21(getMaterialPropertyByName<RankTwoTensor>("Jacobian_21")),
66 979 : _original_length(getMaterialPropertyByName<Real>("original_length")),
67 979 : _total_rotation(getMaterialPropertyByName<RankTwoTensor>("total_rotation")),
68 1958 : _zeta(getMaterialProperty<Real>("zeta")),
69 1958 : _alpha(getParam<Real>("alpha")),
70 2937 : _isDamped(getParam<MaterialPropertyName>("zeta") != "0.0" || std::abs(_alpha) > 0.0),
71 1958 : _force_old(_isDamped ? &getMaterialPropertyOld<RealVectorValue>("forces") : nullptr),
72 1958 : _moment_old(_isDamped ? &getMaterialPropertyOld<RealVectorValue>("moments") : nullptr),
73 1958 : _total_rotation_old(_isDamped ? &getMaterialPropertyOld<RankTwoTensor>("total_rotation")
74 : : nullptr),
75 1033 : _force_older(std::abs(_alpha) > 0.0 ? &getMaterialPropertyOlder<RealVectorValue>("forces")
76 : : nullptr),
77 1033 : _moment_older(std::abs(_alpha) > 0.0 ? &getMaterialPropertyOlder<RealVectorValue>("moments")
78 : : nullptr),
79 979 : _total_rotation_older(std::abs(_alpha) > 0.0
80 1033 : ? &getMaterialPropertyOlder<RankTwoTensor>("total_rotation")
81 : : nullptr),
82 979 : _global_force_res(0),
83 979 : _global_moment_res(0),
84 979 : _force_local_t(0),
85 979 : _moment_local_t(0),
86 979 : _local_force_res(0),
87 1958 : _local_moment_res(0)
88 : {
89 979 : if (_ndisp != _nrot)
90 1 : mooseError(
91 : "StressDivergenceBeam: The number of displacement and rotation variables should be same.");
92 :
93 3912 : for (unsigned int i = 0; i < _ndisp; ++i)
94 2934 : _disp_var[i] = coupled("displacements", i);
95 :
96 3912 : for (unsigned int i = 0; i < _nrot; ++i)
97 2934 : _rot_var[i] = coupled("rotations", i);
98 978 : }
99 :
100 : void
101 933660 : StressDivergenceBeam::computeResidual()
102 : {
103 933660 : prepareVectorTag(_assembly, _var.number());
104 :
105 : mooseAssert(_local_re.size() == 2,
106 : "StressDivergenceBeam: Beam element must have two nodes only.");
107 :
108 933660 : _global_force_res.resize(_test.size());
109 933660 : _global_moment_res.resize(_test.size());
110 :
111 933660 : computeGlobalResidual(&_force, &_moment, &_total_rotation, _global_force_res, _global_moment_res);
112 :
113 : // add contributions from stiffness proportional damping (non-zero _zeta) or HHT time integration
114 : // (non-zero _alpha)
115 933660 : if (_isDamped && _dt > 0.0)
116 933660 : computeDynamicTerms(_global_force_res, _global_moment_res);
117 :
118 2800980 : for (_i = 0; _i < _test.size(); ++_i)
119 : {
120 1867320 : if (_component < 3)
121 933660 : _local_re(_i) = _global_force_res[_i](_component);
122 : else
123 933660 : _local_re(_i) = _global_moment_res[_i](_component - 3);
124 : }
125 :
126 933660 : accumulateTaggedLocalResidual();
127 :
128 933660 : if (_has_save_in)
129 : {
130 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
131 0 : for (_i = 0; _i < _save_in.size(); ++_i)
132 0 : _save_in[_i]->sys().solution().add_vector(_local_re, _save_in[_i]->dofIndices());
133 : }
134 933660 : }
135 :
136 : void
137 289572 : StressDivergenceBeam::computeJacobian()
138 : {
139 289572 : prepareMatrixTag(_assembly, _var.number(), _var.number());
140 :
141 868716 : for (unsigned int i = 0; i < _test.size(); ++i)
142 : {
143 1737432 : for (unsigned int j = 0; j < _phi.size(); ++j)
144 : {
145 1158288 : if (_component < 3)
146 868716 : _local_ke(i, j) = (i == j ? 1 : -1) * _K11[0](_component, _component);
147 : else
148 : {
149 579144 : if (i == j)
150 289572 : _local_ke(i, j) = _K22[0](_component - 3, _component - 3);
151 : else
152 289572 : _local_ke(i, j) = _K22_cross[0](_component - 3, _component - 3);
153 : }
154 : }
155 : }
156 :
157 : // scaling factor for Rayliegh damping and HHT time integration
158 289572 : if (_isDamped && _dt > 0.0)
159 289572 : _local_ke *= (1.0 + _alpha + (1.0 + _alpha) * _zeta[0] / _dt);
160 :
161 289572 : accumulateTaggedLocalMatrix();
162 :
163 289572 : if (_has_diag_save_in)
164 : {
165 : unsigned int rows = _local_ke.m();
166 0 : DenseVector<Number> diag(rows);
167 0 : for (unsigned int i = 0; i < rows; ++i)
168 0 : diag(i) = _local_ke(i, i);
169 :
170 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
171 0 : for (unsigned int i = 0; i < _diag_save_in.size(); ++i)
172 0 : _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
173 : }
174 289572 : }
175 :
176 : void
177 1737432 : StressDivergenceBeam::computeOffDiagJacobian(const unsigned int jvar_num)
178 : {
179 1737432 : if (jvar_num == _var.number())
180 289572 : computeJacobian();
181 : else
182 : {
183 : unsigned int coupled_component = 0;
184 : bool disp_coupled = false;
185 : bool rot_coupled = false;
186 :
187 4343580 : for (unsigned int i = 0; i < _ndisp; ++i)
188 : {
189 3619650 : if (jvar_num == _disp_var[i])
190 : {
191 : coupled_component = i;
192 : disp_coupled = true;
193 : break;
194 : }
195 : }
196 :
197 4343580 : for (unsigned int i = 0; i < _nrot; ++i)
198 : {
199 3619650 : if (jvar_num == _rot_var[i])
200 : {
201 723930 : coupled_component = i + 3;
202 : rot_coupled = true;
203 723930 : break;
204 : }
205 : }
206 :
207 1447860 : prepareMatrixTag(_assembly, _var.number(), jvar_num);
208 :
209 1447860 : if (disp_coupled || rot_coupled)
210 : {
211 4343580 : for (unsigned int i = 0; i < _test.size(); ++i)
212 : {
213 8687160 : for (unsigned int j = 0; j < _phi.size(); ++j)
214 : {
215 5791440 : if (_component < 3 && coupled_component < 3)
216 1737432 : _local_ke(i, j) += (i == j ? 1 : -1) * _K11[0](_component, coupled_component);
217 4633152 : else if (_component < 3 && coupled_component > 2)
218 : {
219 1737432 : if (i == 0)
220 868716 : _local_ke(i, j) += _K21[0](coupled_component - 3, _component);
221 : else
222 868716 : _local_ke(i, j) += _K21_cross[0](coupled_component - 3, _component);
223 : }
224 2895720 : else if (_component > 2 && coupled_component < 3)
225 : {
226 1737432 : if (j == 0)
227 868716 : _local_ke(i, j) += _K21[0](_component - 3, coupled_component);
228 : else
229 868716 : _local_ke(i, j) += _K21_cross[0](_component - 3, coupled_component);
230 : }
231 : else
232 : {
233 1158288 : if (i == j)
234 579144 : _local_ke(i, j) += _K22[0](_component - 3, coupled_component - 3);
235 : else
236 579144 : _local_ke(i, j) += _K22_cross[0](_component - 3, coupled_component - 3);
237 : }
238 : }
239 : }
240 : }
241 :
242 : // scaling factor for Rayleigh damping and HHT time integration
243 1447860 : if (_isDamped && _dt > 0.0)
244 1447860 : _local_ke *= (1.0 + _alpha + (1.0 + _alpha) * _zeta[0] / _dt);
245 :
246 1447860 : accumulateTaggedLocalMatrix();
247 : }
248 1737432 : }
249 :
250 : void
251 933660 : StressDivergenceBeam::computeDynamicTerms(std::vector<RealVectorValue> & global_force_res,
252 : std::vector<RealVectorValue> & global_moment_res)
253 : {
254 : mooseAssert(_zeta[0] >= 0.0, "StressDivergenceBeam: Zeta parameter should be non-negative.");
255 933660 : std::vector<RealVectorValue> global_force_res_old(_test.size());
256 933660 : std::vector<RealVectorValue> global_moment_res_old(_test.size());
257 933660 : computeGlobalResidual(
258 : _force_old, _moment_old, _total_rotation_old, global_force_res_old, global_moment_res_old);
259 :
260 : // For HHT calculation, the global force and moment residual from t_older is required
261 933660 : std::vector<RealVectorValue> global_force_res_older(_test.size());
262 933660 : std::vector<RealVectorValue> global_moment_res_older(_test.size());
263 :
264 933660 : if (std::abs(_alpha) > 0.0)
265 27000 : computeGlobalResidual(_force_older,
266 : _moment_older,
267 : _total_rotation_older,
268 : global_force_res_older,
269 : global_moment_res_older);
270 :
271 : // Update the global_force_res and global_moment_res to include HHT and Rayleigh damping
272 : // contributions
273 2800980 : for (_i = 0; _i < _test.size(); ++_i)
274 : {
275 1867320 : global_force_res[_i] =
276 1867320 : global_force_res[_i] * (1.0 + _alpha + (1.0 + _alpha) * _zeta[0] / _dt) -
277 1867320 : global_force_res_old[_i] * (_alpha + (1.0 + 2.0 * _alpha) * _zeta[0] / _dt) +
278 3734640 : global_force_res_older[_i] * (_alpha * _zeta[0] / _dt);
279 1867320 : global_moment_res[_i] =
280 1867320 : global_moment_res[_i] * (1.0 + _alpha + (1.0 + _alpha) * _zeta[0] / _dt) -
281 1867320 : global_moment_res_old[_i] * (_alpha + (1.0 + 2.0 * _alpha) * _zeta[0] / _dt) +
282 1867320 : global_moment_res_older[_i] * (_alpha * _zeta[0] / _dt);
283 : }
284 933660 : }
285 :
286 : void
287 1894320 : StressDivergenceBeam::computeGlobalResidual(const MaterialProperty<RealVectorValue> * force,
288 : const MaterialProperty<RealVectorValue> * moment,
289 : const MaterialProperty<RankTwoTensor> * total_rotation,
290 : std::vector<RealVectorValue> & global_force_res,
291 : std::vector<RealVectorValue> & global_moment_res)
292 : {
293 : RealVectorValue a;
294 1894320 : _force_local_t.resize(_qrule->n_points());
295 1894320 : _moment_local_t.resize(_qrule->n_points());
296 1894320 : _local_force_res.resize(_test.size());
297 1894320 : _local_moment_res.resize(_test.size());
298 :
299 : // convert forces/moments from global coordinate system to current beam local configuration
300 5682960 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
301 : {
302 3788640 : _force_local_t[_qp] = (*total_rotation)[0] * (*force)[_qp];
303 3788640 : _moment_local_t[_qp] = (*total_rotation)[0] * (*moment)[_qp];
304 : }
305 :
306 : // residual for displacement variables
307 5682960 : for (_i = 0; _i < _test.size(); ++_i)
308 : {
309 3788640 : _local_force_res[_i] = a;
310 15154560 : for (unsigned int component = 0; component < 3; ++component)
311 : {
312 34097760 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
313 22731840 : _local_force_res[_i](component) +=
314 34097760 : (_i == 0 ? -1 : 1) * _force_local_t[_qp](component) * 0.5;
315 : }
316 : }
317 :
318 : // residual for rotation variables
319 5682960 : for (_i = 0; _i < _test.size(); ++_i)
320 : {
321 3788640 : _local_moment_res[_i] = a;
322 15154560 : for (unsigned int component = 3; component < 6; ++component)
323 : {
324 34097760 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
325 : {
326 22731840 : if (component == 3)
327 7577280 : _local_moment_res[_i](component - 3) +=
328 11365920 : (_i == 0 ? -1 : 1) * _moment_local_t[_qp](0) * 0.5;
329 15154560 : else if (component == 4)
330 7577280 : _local_moment_res[_i](component - 3) +=
331 11365920 : (_i == 0 ? -1 : 1) * _moment_local_t[_qp](1) * 0.5 +
332 7577280 : _force_local_t[_qp](2) * 0.25 * _original_length[0];
333 : else
334 7577280 : _local_moment_res[_i](component - 3) +=
335 11365920 : (_i == 0 ? -1 : 1) * _moment_local_t[_qp](2) * 0.5 -
336 7577280 : _force_local_t[_qp](1) * 0.25 * _original_length[0];
337 : }
338 : }
339 : }
340 :
341 : // convert residual for each variable from current beam local configuration to global
342 : // configuration
343 5682960 : for (_i = 0; _i < _test.size(); ++_i)
344 : {
345 3788640 : global_force_res[_i] = (*total_rotation)[0].transpose() * _local_force_res[_i];
346 3788640 : global_moment_res[_i] = (*total_rotation)[0].transpose() * _local_moment_res[_i];
347 : }
348 1894320 : }
|
