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 "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("SolidMechanicsApp", StressDivergenceBeam);
24 :
25 : InputParameters
26 3757 : StressDivergenceBeam::validParams()
27 : {
28 3757 : InputParameters params = Kernel::validParams();
29 3757 : params.addClassDescription("Quasi-static and dynamic stress divergence kernel for Beam element");
30 7514 : 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 7514 : params.addRequiredCoupledVar(
36 : "displacements",
37 : "The displacements appropriate for the simulation geometry and coordinate system");
38 7514 : params.addRequiredCoupledVar(
39 : "rotations", "The rotations appropriate for the simulation geometry and coordinate system");
40 7514 : params.addParam<MaterialPropertyName>(
41 : "zeta",
42 7514 : 0.0,
43 : "Name of material property or a constant real number defining the zeta parameter for the "
44 : "Rayleigh damping.");
45 11271 : params.addRangeCheckedParam<Real>(
46 7514 : "alpha", 0.0, "alpha >= -0.3333 & alpha <= 0.0", "alpha parameter for HHT time integration");
47 :
48 3757 : params.set<bool>("use_displaced_mesh") = true;
49 3757 : return params;
50 0 : }
51 :
52 2264 : StressDivergenceBeam::StressDivergenceBeam(const InputParameters & parameters)
53 : : Kernel(parameters),
54 2264 : _component(getParam<unsigned int>("component")),
55 2264 : _ndisp(coupledComponents("displacements")),
56 2264 : _disp_var(_ndisp),
57 2264 : _nrot(coupledComponents("rotations")),
58 2264 : _rot_var(_nrot),
59 4528 : _force(getMaterialPropertyByName<RealVectorValue>("forces")),
60 2264 : _moment(getMaterialPropertyByName<RealVectorValue>("moments")),
61 2264 : _K11(getMaterialPropertyByName<RankTwoTensor>("Jacobian_11")),
62 2264 : _K22(getMaterialPropertyByName<RankTwoTensor>("Jacobian_22")),
63 2264 : _K22_cross(getMaterialPropertyByName<RankTwoTensor>("Jacobian_22_cross")),
64 2264 : _K21_cross(getMaterialPropertyByName<RankTwoTensor>("Jacobian_12")),
65 2264 : _K21(getMaterialPropertyByName<RankTwoTensor>("Jacobian_21")),
66 2264 : _original_length(getMaterialPropertyByName<Real>("original_length")),
67 2264 : _total_rotation(getMaterialPropertyByName<RankTwoTensor>("total_rotation")),
68 4528 : _zeta(getMaterialProperty<Real>("zeta")),
69 4528 : _alpha(getParam<Real>("alpha")),
70 6792 : _isDamped(getParam<MaterialPropertyName>("zeta") != "0.0" || std::abs(_alpha) > 0.0),
71 4528 : _force_old(_isDamped ? &getMaterialPropertyOld<RealVectorValue>("forces") : nullptr),
72 4528 : _moment_old(_isDamped ? &getMaterialPropertyOld<RealVectorValue>("moments") : nullptr),
73 4528 : _total_rotation_old(_isDamped ? &getMaterialPropertyOld<RankTwoTensor>("total_rotation")
74 : : nullptr),
75 2390 : _force_older(std::abs(_alpha) > 0.0 ? &getMaterialPropertyOlder<RealVectorValue>("forces")
76 : : nullptr),
77 2390 : _moment_older(std::abs(_alpha) > 0.0 ? &getMaterialPropertyOlder<RealVectorValue>("moments")
78 : : nullptr),
79 2264 : _total_rotation_older(std::abs(_alpha) > 0.0
80 2390 : ? &getMaterialPropertyOlder<RankTwoTensor>("total_rotation")
81 : : nullptr),
82 2264 : _global_force_res(0),
83 2264 : _global_moment_res(0),
84 2264 : _force_local_t(0),
85 2264 : _moment_local_t(0),
86 2264 : _local_force_res(0),
87 4528 : _local_moment_res(0)
88 : {
89 2264 : if (_ndisp != _nrot)
90 2 : mooseError(
91 : "StressDivergenceBeam: The number of displacement and rotation variables should be same.");
92 :
93 9048 : for (unsigned int i = 0; i < _ndisp; ++i)
94 6786 : _disp_var[i] = coupled("displacements", i);
95 :
96 9048 : for (unsigned int i = 0; i < _nrot; ++i)
97 6786 : _rot_var[i] = coupled("rotations", i);
98 2262 : }
99 :
100 : void
101 1758846 : StressDivergenceBeam::computeResidual()
102 : {
103 1758846 : prepareVectorTag(_assembly, _var.number());
104 :
105 : mooseAssert(_local_re.size() == 2,
106 : "StressDivergenceBeam: Beam element must have two nodes only.");
107 :
108 1758846 : _global_force_res.resize(_test.size());
109 1758846 : _global_moment_res.resize(_test.size());
110 :
111 1758846 : 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 1758846 : if (_isDamped && _dt > 0.0)
116 1758846 : computeDynamicTerms(_global_force_res, _global_moment_res);
117 :
118 5276538 : for (_i = 0; _i < _test.size(); ++_i)
119 : {
120 3517692 : if (_component < 3)
121 1758846 : _local_re(_i) = _global_force_res[_i](_component);
122 : else
123 1758846 : _local_re(_i) = _global_moment_res[_i](_component - 3);
124 : }
125 :
126 1758846 : accumulateTaggedLocalResidual();
127 :
128 1758846 : 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 1758846 : }
135 :
136 : void
137 724920 : StressDivergenceBeam::computeJacobian()
138 : {
139 724920 : prepareMatrixTag(_assembly, _var.number(), _var.number());
140 :
141 2174760 : for (unsigned int i = 0; i < _test.size(); ++i)
142 : {
143 4349520 : for (unsigned int j = 0; j < _phi.size(); ++j)
144 : {
145 2899680 : if (_component < 3)
146 2174760 : _local_ke(i, j) = (i == j ? 1 : -1) * _K11[0](_component, _component);
147 : else
148 : {
149 1449840 : if (i == j)
150 724920 : _local_ke(i, j) = _K22[0](_component - 3, _component - 3);
151 : else
152 724920 : _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 724920 : if (_isDamped && _dt > 0.0)
159 724920 : _local_ke *= (1.0 + _alpha + (1.0 + _alpha) * _zeta[0] / _dt);
160 :
161 724920 : accumulateTaggedLocalMatrix();
162 :
163 724920 : 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 724920 : }
175 :
176 : void
177 4349520 : StressDivergenceBeam::computeOffDiagJacobian(const unsigned int jvar_num)
178 : {
179 4349520 : if (jvar_num == _var.number())
180 724920 : computeJacobian();
181 : else
182 : {
183 : unsigned int coupled_component = 0;
184 : bool disp_coupled = false;
185 : bool rot_coupled = false;
186 :
187 10873800 : for (unsigned int i = 0; i < _ndisp; ++i)
188 : {
189 9061500 : if (jvar_num == _disp_var[i])
190 : {
191 : coupled_component = i;
192 : disp_coupled = true;
193 : break;
194 : }
195 : }
196 :
197 10873800 : for (unsigned int i = 0; i < _nrot; ++i)
198 : {
199 9061500 : if (jvar_num == _rot_var[i])
200 : {
201 1812300 : coupled_component = i + 3;
202 : rot_coupled = true;
203 1812300 : break;
204 : }
205 : }
206 :
207 3624600 : prepareMatrixTag(_assembly, _var.number(), jvar_num);
208 :
209 3624600 : if (disp_coupled || rot_coupled)
210 : {
211 10873800 : for (unsigned int i = 0; i < _test.size(); ++i)
212 : {
213 21747600 : for (unsigned int j = 0; j < _phi.size(); ++j)
214 : {
215 14498400 : if (_component < 3 && coupled_component < 3)
216 4349520 : _local_ke(i, j) += (i == j ? 1 : -1) * _K11[0](_component, coupled_component);
217 11598720 : else if (_component < 3 && coupled_component > 2)
218 : {
219 4349520 : if (i == 0)
220 2174760 : _local_ke(i, j) += _K21[0](coupled_component - 3, _component);
221 : else
222 2174760 : _local_ke(i, j) += _K21_cross[0](coupled_component - 3, _component);
223 : }
224 7249200 : else if (_component > 2 && coupled_component < 3)
225 : {
226 4349520 : if (j == 0)
227 2174760 : _local_ke(i, j) += _K21[0](_component - 3, coupled_component);
228 : else
229 2174760 : _local_ke(i, j) += _K21_cross[0](_component - 3, coupled_component);
230 : }
231 : else
232 : {
233 2899680 : if (i == j)
234 1449840 : _local_ke(i, j) += _K22[0](_component - 3, coupled_component - 3);
235 : else
236 1449840 : _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 3624600 : if (_isDamped && _dt > 0.0)
244 3624600 : _local_ke *= (1.0 + _alpha + (1.0 + _alpha) * _zeta[0] / _dt);
245 :
246 3624600 : accumulateTaggedLocalMatrix();
247 : }
248 4349520 : }
249 :
250 : void
251 1758846 : 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 1758846 : std::vector<RealVectorValue> global_force_res_old(_test.size());
256 1758846 : std::vector<RealVectorValue> global_moment_res_old(_test.size());
257 1758846 : 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 3517692 : std::vector<RealVectorValue> global_force_res_older(_test.size());
262 1758846 : std::vector<RealVectorValue> global_moment_res_older(_test.size());
263 :
264 1758846 : if (std::abs(_alpha) > 0.0)
265 45360 : 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 5276538 : for (_i = 0; _i < _test.size(); ++_i)
274 : {
275 3517692 : global_force_res[_i] =
276 3517692 : global_force_res[_i] * (1.0 + _alpha + (1.0 + _alpha) * _zeta[0] / _dt) -
277 3517692 : global_force_res_old[_i] * (_alpha + (1.0 + 2.0 * _alpha) * _zeta[0] / _dt) +
278 3517692 : global_force_res_older[_i] * (_alpha * _zeta[0] / _dt);
279 3517692 : global_moment_res[_i] =
280 3517692 : global_moment_res[_i] * (1.0 + _alpha + (1.0 + _alpha) * _zeta[0] / _dt) -
281 3517692 : global_moment_res_old[_i] * (_alpha + (1.0 + 2.0 * _alpha) * _zeta[0] / _dt) +
282 3517692 : global_moment_res_older[_i] * (_alpha * _zeta[0] / _dt);
283 : }
284 1758846 : }
285 :
286 : void
287 3563052 : 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 3563052 : _force_local_t.resize(_qrule->n_points());
295 3563052 : _moment_local_t.resize(_qrule->n_points());
296 3563052 : _local_force_res.resize(_test.size());
297 3563052 : _local_moment_res.resize(_test.size());
298 :
299 : // convert forces/moments from global coordinate system to current beam local configuration
300 10689156 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
301 : {
302 7126104 : _force_local_t[_qp] = (*total_rotation)[0] * (*force)[_qp];
303 7126104 : _moment_local_t[_qp] = (*total_rotation)[0] * (*moment)[_qp];
304 : }
305 :
306 : // residual for displacement variables
307 10689156 : for (_i = 0; _i < _test.size(); ++_i)
308 : {
309 7126104 : _local_force_res[_i] = a;
310 28504416 : for (unsigned int component = 0; component < 3; ++component)
311 : {
312 64134936 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
313 42756624 : _local_force_res[_i](component) +=
314 64134936 : (_i == 0 ? -1 : 1) * _force_local_t[_qp](component) * 0.5;
315 : }
316 : }
317 :
318 : // residual for rotation variables
319 10689156 : for (_i = 0; _i < _test.size(); ++_i)
320 : {
321 7126104 : _local_moment_res[_i] = a;
322 28504416 : for (unsigned int component = 3; component < 6; ++component)
323 : {
324 64134936 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
325 : {
326 42756624 : if (component == 3)
327 14252208 : _local_moment_res[_i](component - 3) +=
328 21378312 : (_i == 0 ? -1 : 1) * _moment_local_t[_qp](0) * 0.5;
329 28504416 : else if (component == 4)
330 14252208 : _local_moment_res[_i](component - 3) +=
331 21378312 : (_i == 0 ? -1 : 1) * _moment_local_t[_qp](1) * 0.5 +
332 14252208 : _force_local_t[_qp](2) * 0.25 * _original_length[0];
333 : else
334 14252208 : _local_moment_res[_i](component - 3) +=
335 21378312 : (_i == 0 ? -1 : 1) * _moment_local_t[_qp](2) * 0.5 -
336 14252208 : _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 10689156 : for (_i = 0; _i < _test.size(); ++_i)
344 : {
345 7126104 : global_force_res[_i] = (*total_rotation)[0].transpose() * _local_force_res[_i];
346 7126104 : global_moment_res[_i] = (*total_rotation)[0].transpose() * _local_moment_res[_i];
347 : }
348 3563052 : }
|
