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