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 "InertialForceBeam.h"
11 : #include "SubProblem.h"
12 : #include "libmesh/utility.h"
13 : #include "MooseVariable.h"
14 : #include "Assembly.h"
15 : #include "NonlinearSystem.h"
16 : #include "AuxiliarySystem.h"
17 : #include "MooseMesh.h"
18 :
19 : registerMooseObject("SolidMechanicsApp", InertialForceBeam);
20 :
21 : InputParameters
22 232 : InertialForceBeam::validParams()
23 : {
24 232 : InputParameters params = TimeKernel::validParams();
25 232 : params.addClassDescription("Calculates the residual for the inertial force/moment and the "
26 : "contribution of mass dependent Rayleigh damping and HHT time "
27 : "integration scheme.");
28 232 : params.set<bool>("use_displaced_mesh") = true;
29 464 : params.addRequiredCoupledVar(
30 : "rotations",
31 : "The rotational variables appropriate for the simulation geometry and coordinate system");
32 464 : params.addRequiredCoupledVar(
33 : "displacements",
34 : "The displacement variables appropriate for the simulation geometry and coordinate system");
35 464 : params.addCoupledVar("velocities", "Translational velocity variables");
36 464 : params.addCoupledVar("accelerations", "Translational acceleration variables");
37 464 : params.addCoupledVar("rotational_velocities", "Rotational velocity variables");
38 464 : params.addCoupledVar("rotational_accelerations", "Rotational acceleration variables");
39 464 : params.addRangeCheckedParam<Real>(
40 : "beta", "beta>0.0", "beta parameter for Newmark Time integration");
41 464 : params.addRangeCheckedParam<Real>(
42 : "gamma", "gamma>0.0", "gamma parameter for Newmark Time integration");
43 464 : params.addParam<MaterialPropertyName>("eta",
44 464 : 0.0,
45 : "Name of material property or a constant real "
46 : "number defining the eta parameter for the "
47 : "Rayleigh damping.");
48 696 : params.addRangeCheckedParam<Real>("alpha",
49 464 : 0.0,
50 : "alpha >= -0.3333 & alpha <= 0.0",
51 : "Alpha parameter for mass dependent numerical damping induced "
52 : "by HHT time integration scheme");
53 464 : params.addParam<MaterialPropertyName>(
54 : "density",
55 : "density",
56 : "Name of Material Property or a constant real number defining the density of the beam.");
57 464 : params.addRequiredCoupledVar("area", "Variable containing cross-section area");
58 464 : params.addCoupledVar("Ay", 0.0, "Variable containing first moment of area about y axis");
59 464 : params.addCoupledVar("Az", 0.0, "Variable containing first moment of area about z axis");
60 464 : params.addCoupledVar("Ix",
61 : "Variable containing second moment of area about x axis. Defaults to Iy+Iz");
62 464 : params.addRequiredCoupledVar("Iy", "Variable containing second moment of area about y axis");
63 464 : params.addRequiredCoupledVar("Iz", "Variable containing second moment of area about z axis");
64 464 : params.addRequiredRangeCheckedParam<unsigned int>(
65 : "component",
66 : "component<6",
67 : "An integer corresponding to the direction "
68 : "the variable this kernel acts in. (0 for disp_x, "
69 : "1 for disp_y, 2 for disp_z, 3 for rot_x, 4 for rot_y and 5 for rot_z)");
70 232 : return params;
71 0 : }
72 :
73 122 : InertialForceBeam::InertialForceBeam(const InputParameters & parameters)
74 : : TimeKernel(parameters),
75 122 : _has_beta(isParamValid("beta")),
76 244 : _has_gamma(isParamValid("gamma")),
77 244 : _has_velocities(isParamValid("velocities")),
78 244 : _has_rot_velocities(isParamValid("rotational_velocities")),
79 244 : _has_accelerations(isParamValid("accelerations")),
80 244 : _has_rot_accelerations(isParamValid("rotational_accelerations")),
81 244 : _has_Ix(isParamValid("Ix")),
82 244 : _density(getMaterialProperty<Real>("density")),
83 122 : _nrot(coupledComponents("rotations")),
84 122 : _ndisp(coupledComponents("displacements")),
85 122 : _rot_num(_nrot),
86 122 : _disp_num(_ndisp),
87 122 : _vel_num(_ndisp),
88 122 : _accel_num(_ndisp),
89 122 : _rot_vel_num(_nrot),
90 122 : _rot_accel_num(_nrot),
91 122 : _area(coupledValue("area")),
92 122 : _Ay(coupledValue("Ay")),
93 122 : _Az(coupledValue("Az")),
94 122 : _Ix(_has_Ix ? coupledValue("Ix") : _zero),
95 122 : _Iy(coupledValue("Iy")),
96 122 : _Iz(coupledValue("Iz")),
97 220 : _beta(_has_beta ? getParam<Real>("beta") : 0.1),
98 220 : _gamma(_has_gamma ? getParam<Real>("gamma") : 0.1),
99 244 : _eta(getMaterialProperty<Real>("eta")),
100 244 : _alpha(getParam<Real>("alpha")),
101 122 : _original_local_config(getMaterialPropertyByName<RankTwoTensor>("initial_rotation")),
102 122 : _original_length(getMaterialPropertyByName<Real>("original_length")),
103 244 : _component(getParam<unsigned int>("component")),
104 122 : _local_force(2),
105 244 : _local_moment(2)
106 : {
107 : // Checking for consistency between the length of the provided rotations and displacements vector
108 122 : if (_ndisp != _nrot)
109 1 : mooseError("InertialForceBeam: The number of variables supplied in 'displacements' and "
110 : "'rotations' must match.");
111 :
112 121 : if (_has_beta && _has_gamma && _has_velocities && _has_accelerations && _has_rot_velocities &&
113 97 : _has_rot_accelerations)
114 : {
115 194 : if ((coupledComponents("velocities") != _ndisp) ||
116 194 : (coupledComponents("accelerations") != _ndisp) ||
117 290 : (coupledComponents("rotational_velocities") != _ndisp) ||
118 193 : (coupledComponents("rotational_accelerations") != _ndisp))
119 1 : mooseError(
120 : "InertialForceBeam: The number of variables supplied in 'velocities', "
121 : "'accelerations', 'rotational_velocities' and 'rotational_accelerations' must match "
122 : "the number of displacement variables.");
123 :
124 : // fetch coupled velocities and accelerations
125 384 : for (unsigned int i = 0; i < _ndisp; ++i)
126 : {
127 576 : MooseVariable * vel_variable = getVar("velocities", i);
128 288 : _vel_num[i] = vel_variable->number();
129 :
130 576 : MooseVariable * accel_variable = getVar("accelerations", i);
131 288 : _accel_num[i] = accel_variable->number();
132 :
133 576 : MooseVariable * rot_vel_variable = getVar("rotational_velocities", i);
134 288 : _rot_vel_num[i] = rot_vel_variable->number();
135 :
136 576 : MooseVariable * rot_accel_variable = getVar("rotational_accelerations", i);
137 288 : _rot_accel_num[i] = rot_accel_variable->number();
138 : }
139 : }
140 24 : else if (!_has_beta && !_has_gamma && !_has_velocities && !_has_accelerations &&
141 24 : !_has_rot_velocities && !_has_rot_accelerations)
142 : {
143 24 : _du_dot_du = &coupledDotDu("displacements", 0);
144 24 : _du_dotdot_du = &coupledDotDotDu("displacements", 0);
145 : }
146 : else
147 0 : mooseError("InertialForceBeam: Either all or none of `beta`, `gamma`, `velocities`, "
148 : "`accelerations`, `rotational_velocities` and `rotational_accelerations` should be "
149 : "provided as input.");
150 :
151 : // fetch coupled displacements and rotations
152 480 : for (unsigned int i = 0; i < _ndisp; ++i)
153 : {
154 720 : MooseVariable * disp_variable = getVar("displacements", i);
155 360 : _disp_num[i] = disp_variable->number();
156 :
157 720 : MooseVariable * rot_variable = getVar("rotations", i);
158 360 : _rot_num[i] = rot_variable->number();
159 : }
160 120 : }
161 :
162 : void
163 430680 : InertialForceBeam::computeResidual()
164 : {
165 430680 : prepareVectorTag(_assembly, _var.number());
166 :
167 430680 : if (_dt != 0.0)
168 : {
169 : // fetch the two end nodes for _current_elem
170 : std::vector<const Node *> node;
171 1292040 : for (unsigned int i = 0; i < 2; ++i)
172 861360 : node.push_back(_current_elem->node_ptr(i));
173 :
174 : // Fetch the solution for the two end nodes at time t
175 430680 : NonlinearSystemBase & nonlinear_sys = _fe_problem.getNonlinearSystemBase(_sys.number());
176 :
177 430680 : if (_has_beta)
178 : {
179 421560 : const NumericVector<Number> & sol = *nonlinear_sys.currentSolution();
180 421560 : const NumericVector<Number> & sol_old = nonlinear_sys.solutionOld();
181 :
182 421560 : AuxiliarySystem & aux = _fe_problem.getAuxiliarySystem();
183 421560 : const NumericVector<Number> & aux_sol_old = aux.solutionOld();
184 :
185 : mooseAssert(_beta > 0.0, "InertialForceBeam: Beta parameter should be positive.");
186 : mooseAssert(_eta[0] >= 0.0, "InertialForceBeam: Eta parameter should be non-negative.");
187 :
188 1686240 : for (unsigned int i = 0; i < _ndisp; ++i)
189 : {
190 : // obtain delta displacement
191 1264680 : unsigned int dof_index_0 = node[0]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
192 1264680 : unsigned int dof_index_1 = node[1]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
193 1264680 : const Real disp_0 = sol(dof_index_0) - sol_old(dof_index_0);
194 1264680 : const Real disp_1 = sol(dof_index_1) - sol_old(dof_index_1);
195 :
196 1264680 : dof_index_0 = node[0]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
197 1264680 : dof_index_1 = node[1]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
198 1264680 : const Real rot_0 = sol(dof_index_0) - sol_old(dof_index_0);
199 1264680 : const Real rot_1 = sol(dof_index_1) - sol_old(dof_index_1);
200 :
201 : // obtain new translational and rotational velocities and accelerations using newmark-beta
202 : // time integration
203 1264680 : _vel_old_0(i) = aux_sol_old(node[0]->dof_number(aux.number(), _vel_num[i], 0));
204 1264680 : _vel_old_1(i) = aux_sol_old(node[1]->dof_number(aux.number(), _vel_num[i], 0));
205 1264680 : const Real accel_old_0 = aux_sol_old(node[0]->dof_number(aux.number(), _accel_num[i], 0));
206 1264680 : const Real accel_old_1 = aux_sol_old(node[1]->dof_number(aux.number(), _accel_num[i], 0));
207 :
208 1264680 : _rot_vel_old_0(i) = aux_sol_old(node[0]->dof_number(aux.number(), _rot_vel_num[i], 0));
209 1264680 : _rot_vel_old_1(i) = aux_sol_old(node[1]->dof_number(aux.number(), _rot_vel_num[i], 0));
210 : const Real rot_accel_old_0 =
211 1264680 : aux_sol_old(node[0]->dof_number(aux.number(), _rot_accel_num[i], 0));
212 : const Real rot_accel_old_1 =
213 1264680 : aux_sol_old(node[1]->dof_number(aux.number(), _rot_accel_num[i], 0));
214 :
215 1264680 : _accel_0(i) =
216 1264680 : 1. / _beta * (disp_0 / (_dt * _dt) - _vel_old_0(i) / _dt - accel_old_0 * (0.5 - _beta));
217 1264680 : _accel_1(i) =
218 1264680 : 1. / _beta * (disp_1 / (_dt * _dt) - _vel_old_1(i) / _dt - accel_old_1 * (0.5 - _beta));
219 1264680 : _rot_accel_0(i) =
220 1264680 : 1. / _beta *
221 1264680 : (rot_0 / (_dt * _dt) - _rot_vel_old_0(i) / _dt - rot_accel_old_0 * (0.5 - _beta));
222 1264680 : _rot_accel_1(i) =
223 1264680 : 1. / _beta *
224 1264680 : (rot_1 / (_dt * _dt) - _rot_vel_old_1(i) / _dt - rot_accel_old_1 * (0.5 - _beta));
225 :
226 1264680 : _vel_0(i) = _vel_old_0(i) + (_dt * (1 - _gamma)) * accel_old_0 + _gamma * _dt * _accel_0(i);
227 1264680 : _vel_1(i) = _vel_old_1(i) + (_dt * (1 - _gamma)) * accel_old_1 + _gamma * _dt * _accel_1(i);
228 1264680 : _rot_vel_0(i) = _rot_vel_old_0(i) + (_dt * (1 - _gamma)) * rot_accel_old_0 +
229 1264680 : _gamma * _dt * _rot_accel_0(i);
230 1264680 : _rot_vel_1(i) = _rot_vel_old_1(i) + (_dt * (1 - _gamma)) * rot_accel_old_1 +
231 1264680 : _gamma * _dt * _rot_accel_1(i);
232 : }
233 : }
234 : else
235 : {
236 9120 : if (!nonlinear_sys.solutionUDot())
237 0 : mooseError("InertialForceBeam: Time derivative of solution (`u_dot`) is not stored. Please "
238 : "set uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
239 :
240 9120 : if (!nonlinear_sys.solutionUDotOld())
241 0 : mooseError("InertialForceBeam: Old time derivative of solution (`u_dot_old`) is not "
242 : "stored. Please set uDotOldRequested() to true in FEProblemBase before "
243 : "requesting `u_dot_old`.");
244 :
245 9120 : if (!nonlinear_sys.solutionUDotDot())
246 0 : mooseError("InertialForceBeam: Second time derivative of solution (`u_dotdot`) is not "
247 : "stored. Please set uDotDotRequested() to true in FEProblemBase before "
248 : "requesting `u_dotdot`.");
249 :
250 9120 : const NumericVector<Number> & vel = *nonlinear_sys.solutionUDot();
251 9120 : const NumericVector<Number> & vel_old = *nonlinear_sys.solutionUDotOld();
252 9120 : const NumericVector<Number> & accel = *nonlinear_sys.solutionUDotDot();
253 :
254 36480 : for (unsigned int i = 0; i < _ndisp; ++i)
255 : {
256 : // translational velocities and accelerations
257 27360 : unsigned int dof_index_0 = node[0]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
258 27360 : unsigned int dof_index_1 = node[1]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
259 27360 : _vel_0(i) = vel(dof_index_0);
260 27360 : _vel_1(i) = vel(dof_index_1);
261 27360 : _vel_old_0(i) = vel_old(dof_index_0);
262 27360 : _vel_old_1(i) = vel_old(dof_index_1);
263 27360 : _accel_0(i) = accel(dof_index_0);
264 27360 : _accel_1(i) = accel(dof_index_1);
265 :
266 : // rotational velocities and accelerations
267 27360 : dof_index_0 = node[0]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
268 27360 : dof_index_1 = node[1]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
269 27360 : _rot_vel_0(i) = vel(dof_index_0);
270 27360 : _rot_vel_1(i) = vel(dof_index_1);
271 27360 : _rot_vel_old_0(i) = vel_old(dof_index_0);
272 27360 : _rot_vel_old_1(i) = vel_old(dof_index_1);
273 27360 : _rot_accel_0(i) = accel(dof_index_0);
274 27360 : _rot_accel_1(i) = accel(dof_index_1);
275 : }
276 : }
277 :
278 : // transform translational and rotational velocities and accelerations to the initial local
279 : // configuration of the beam
280 430680 : _local_vel_old_0 = _original_local_config[0] * _vel_old_0;
281 430680 : _local_vel_old_1 = _original_local_config[0] * _vel_old_1;
282 430680 : _local_vel_0 = _original_local_config[0] * _vel_0;
283 430680 : _local_vel_1 = _original_local_config[0] * _vel_1;
284 430680 : _local_accel_0 = _original_local_config[0] * _accel_0;
285 430680 : _local_accel_1 = _original_local_config[0] * _accel_1;
286 :
287 430680 : _local_rot_vel_old_0 = _original_local_config[0] * _rot_vel_old_0;
288 430680 : _local_rot_vel_old_1 = _original_local_config[0] * _rot_vel_old_1;
289 430680 : _local_rot_vel_0 = _original_local_config[0] * _rot_vel_0;
290 430680 : _local_rot_vel_1 = _original_local_config[0] * _rot_vel_1;
291 430680 : _local_rot_accel_0 = _original_local_config[0] * _rot_accel_0;
292 430680 : _local_rot_accel_1 = _original_local_config[0] * _rot_accel_1;
293 :
294 : // local residual
295 1722720 : for (unsigned int i = 0; i < _ndisp; ++i)
296 : {
297 1292040 : if (_component < 3)
298 : {
299 646020 : _local_force[0](i) = _density[0] * _area[0] * _original_length[0] / 3.0 *
300 646020 : (_local_accel_0(i) + _local_accel_1(i) / 2.0 +
301 646020 : _eta[0] * (1.0 + _alpha) * (_local_vel_0(i) + _local_vel_1(i) / 2.0) -
302 646020 : _alpha * _eta[0] * (_local_vel_old_0(i) + _local_vel_old_1(i) / 2.0));
303 646020 : _local_force[1](i) = _density[0] * _area[0] * _original_length[0] / 3.0 *
304 646020 : (_local_accel_1(i) + _local_accel_0(i) / 2.0 +
305 646020 : _eta[0] * (1.0 + _alpha) * (_local_vel_1(i) + _local_vel_0(i) / 2.0) -
306 646020 : _alpha * _eta[0] * (_local_vel_old_1(i) + _local_vel_old_0(i) / 2.0));
307 : }
308 :
309 1292040 : if (_component > 2)
310 : {
311 646020 : Real I = _Iy[0] + _Iz[0];
312 646020 : if (_has_Ix && (i == 0))
313 0 : I = _Ix[0];
314 646020 : if (i == 1)
315 : I = _Iz[0];
316 430680 : else if (i == 2)
317 : I = _Iy[0];
318 :
319 646020 : _local_moment[0](i) =
320 646020 : _density[0] * I * _original_length[0] / 3.0 *
321 646020 : (_local_rot_accel_0(i) + _local_rot_accel_1(i) / 2.0 +
322 646020 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_0(i) + _local_rot_vel_1(i) / 2.0) -
323 646020 : _alpha * _eta[0] * (_local_rot_vel_old_0(i) + _local_rot_vel_old_1(i) / 2.0));
324 646020 : _local_moment[1](i) =
325 646020 : _density[0] * I * _original_length[0] / 3.0 *
326 646020 : (_local_rot_accel_1(i) + _local_rot_accel_0(i) / 2.0 +
327 646020 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_1(i) + _local_rot_vel_0(i) / 2.0) -
328 646020 : _alpha * _eta[0] * (_local_rot_vel_old_1(i) + _local_rot_vel_old_0(i) / 2.0));
329 : }
330 : }
331 :
332 : // If Ay or Az are non-zero, contribution of rotational accelerations to translational forces
333 : // and vice versa have to be added
334 430680 : if (_component < 3)
335 : {
336 215340 : _local_force[0](0) +=
337 215340 : _density[0] * _original_length[0] / 3.0 *
338 215340 : (_Az[0] * (_local_rot_accel_0(1) + _local_rot_accel_1(1) / 2.0 +
339 215340 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_0(1) + _local_rot_vel_1(1) / 2.0) -
340 215340 : _alpha * _eta[0] * (_local_rot_vel_old_0(1) + _local_rot_vel_old_1(1) / 2.0)) -
341 215340 : _Ay[0] * (_local_rot_accel_0(2) + _local_rot_accel_1(2) / 2.0 +
342 215340 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_0(2) + _local_rot_vel_1(2) / 2.0) -
343 215340 : _alpha * _eta[0] * (_local_rot_vel_old_0(2) + _local_rot_vel_old_1(2) / 2.0)));
344 215340 : _local_force[1](0) +=
345 215340 : _density[0] * _original_length[0] / 3.0 *
346 215340 : (_Az[0] * (_local_rot_accel_1(1) + _local_rot_accel_0(1) / 2.0 +
347 215340 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_1(1) + _local_rot_vel_0(1) / 2.0) -
348 215340 : _alpha * _eta[0] * (_local_rot_vel_old_1(1) + _local_rot_vel_old_0(1) / 2.0)) -
349 215340 : _Ay[0] * (_local_rot_accel_1(2) + _local_rot_accel_0(2) / 2.0 +
350 215340 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_1(2) + _local_rot_vel_0(2) / 2.0) -
351 215340 : _alpha * _eta[0] * (_local_rot_vel_old_1(2) + _local_rot_vel_old_0(2) / 2.0)));
352 :
353 215340 : _local_force[0](1) +=
354 215340 : -_density[0] * _original_length[0] / 3.0 * _Az[0] *
355 215340 : (_local_rot_accel_0(0) + _local_rot_accel_1(0) / 2.0 +
356 215340 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_0(0) + _local_rot_vel_1(0) / 2.0) -
357 215340 : _alpha * _eta[0] * (_local_rot_vel_old_0(0) + _local_rot_vel_old_1(0) / 2.0));
358 215340 : _local_force[1](1) +=
359 215340 : -_density[0] * _original_length[0] / 3.0 * _Az[0] *
360 215340 : (_local_rot_accel_1(0) + _local_rot_accel_0(0) / 2.0 +
361 215340 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_1(0) + _local_rot_vel_0(0) / 2.0) -
362 215340 : _alpha * _eta[0] * (_local_rot_vel_old_1(0) + _local_rot_vel_old_0(0) / 2.0));
363 :
364 215340 : _local_force[0](2) +=
365 215340 : _density[0] * _original_length[0] / 3.0 * _Ay[0] *
366 215340 : (_local_rot_accel_0(0) + _local_rot_accel_1(0) / 2.0 +
367 215340 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_0(0) + _local_rot_vel_1(0) / 2.0) -
368 215340 : _alpha * _eta[0] * (_local_rot_vel_old_0(0) + _local_rot_vel_old_1(0) / 2.0));
369 215340 : _local_force[1](2) +=
370 215340 : _density[0] * _original_length[0] / 3.0 * _Ay[0] *
371 215340 : (_local_rot_accel_1(0) + _local_rot_accel_0(0) / 2.0 +
372 215340 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_1(0) + _local_rot_vel_0(0) / 2.0) -
373 215340 : _alpha * _eta[0] * (_local_rot_vel_old_1(0) + _local_rot_vel_old_0(0) / 2.0));
374 : }
375 : else
376 : {
377 215340 : _local_moment[0](0) += _density[0] * _original_length[0] / 3.0 *
378 215340 : (-_Az[0] * (_local_accel_0(1) + _local_accel_1(1) / 2.0) +
379 215340 : _Ay[0] * (_local_accel_0(1) + _local_accel_1(1) / 2.0));
380 215340 : _local_moment[1](0) += _density[0] * _original_length[0] / 3.0 *
381 215340 : (-_Az[0] * (_local_accel_1(1) + _local_accel_0(1) / 2.0) +
382 215340 : _Ay[0] * (_local_accel_1(1) + _local_accel_0(1) / 2.0));
383 :
384 215340 : _local_moment[0](1) += _density[0] * _original_length[0] / 3.0 * _Az[0] *
385 215340 : (_local_accel_0(0) + _local_accel_1(0) / 2.0);
386 215340 : _local_moment[1](1) += _density[0] * _original_length[0] / 3.0 * _Az[0] *
387 215340 : (_local_accel_1(0) + _local_accel_0(0) / 2.0);
388 :
389 215340 : _local_moment[0](2) += -_density[0] * _original_length[0] / 3.0 * _Ay[0] *
390 215340 : (_local_accel_0(0) + _local_accel_1(0) / 2.0);
391 215340 : _local_moment[1](2) += -_density[0] * _original_length[0] / 3.0 * _Ay[0] *
392 215340 : (_local_accel_1(0) + _local_accel_0(0) / 2.0);
393 : }
394 :
395 : // Global force and moments
396 430680 : if (_component < 3)
397 : {
398 215340 : _global_force_0 = _original_local_config[0] * _local_force[0];
399 215340 : _global_force_1 = _original_local_config[0] * _local_force[1];
400 215340 : _local_re(0) = _global_force_0(_component);
401 215340 : _local_re(1) = _global_force_1(_component);
402 : }
403 : else
404 : {
405 215340 : _global_moment_0 = _original_local_config[0] * _local_moment[0];
406 215340 : _global_moment_1 = _original_local_config[0] * _local_moment[1];
407 215340 : _local_re(0) = _global_moment_0(_component - 3);
408 215340 : _local_re(1) = _global_moment_1(_component - 3);
409 : }
410 430680 : }
411 :
412 430680 : accumulateTaggedLocalResidual();
413 :
414 430680 : if (_has_save_in)
415 : {
416 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
417 0 : for (unsigned int i = 0; i < _save_in.size(); ++i)
418 0 : _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
419 : }
420 430680 : }
421 :
422 : void
423 218340 : InertialForceBeam::computeJacobian()
424 : {
425 218340 : prepareMatrixTag(_assembly, _var.number(), _var.number());
426 :
427 : mooseAssert(_beta > 0.0, "InertialForceBeam: Beta parameter should be positive.");
428 :
429 : Real factor = 0.0;
430 218340 : if (_has_beta)
431 213780 : factor = 1.0 / (_beta * _dt * _dt) + _eta[0] * (1.0 + _alpha) * _gamma / _beta / _dt;
432 : else
433 4560 : factor = (*_du_dotdot_du)[0] + _eta[0] * (1.0 + _alpha) * (*_du_dot_du)[0];
434 :
435 655020 : for (unsigned int i = 0; i < _test.size(); ++i)
436 : {
437 1310040 : for (unsigned int j = 0; j < _phi.size(); ++j)
438 : {
439 873360 : if (_component < 3)
440 655020 : _local_ke(i, j) = (i == j ? 1.0 / 3.0 : 1.0 / 6.0) * _density[0] * _area[0] *
441 436680 : _original_length[0] * factor;
442 : else if (_component > 2)
443 : {
444 436680 : RankTwoTensor I;
445 436680 : if (_has_Ix)
446 0 : I(0, 0) = _Ix[0];
447 : else
448 436680 : I(0, 0) = _Iy[0] + _Iz[0];
449 436680 : I(1, 1) = _Iz[0];
450 436680 : I(2, 2) = _Iy[0];
451 :
452 : // conversion from local config to global coordinate system
453 436680 : RankTwoTensor Ig = _original_local_config[0].transpose() * I * _original_local_config[0];
454 :
455 655020 : _local_ke(i, j) = (i == j ? 1.0 / 3.0 : 1.0 / 6.0) * _density[0] *
456 436680 : Ig(_component - 3, _component - 3) * _original_length[0] * factor;
457 : }
458 : }
459 : }
460 :
461 218340 : accumulateTaggedLocalMatrix();
462 :
463 218340 : if (_has_diag_save_in)
464 : {
465 : unsigned int rows = _local_ke.m();
466 0 : DenseVector<Number> diag(rows);
467 0 : for (unsigned int i = 0; i < rows; ++i)
468 0 : diag(i) = _local_ke(i, i);
469 :
470 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
471 0 : for (unsigned int i = 0; i < _diag_save_in.size(); ++i)
472 0 : _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
473 : }
474 218340 : }
475 :
476 : void
477 1310040 : InertialForceBeam::computeOffDiagJacobian(const unsigned int jvar_num)
478 : {
479 : mooseAssert(_beta > 0.0, "InertialForceBeam: Beta parameter should be positive.");
480 :
481 : Real factor = 0.0;
482 1310040 : if (_has_beta)
483 1282680 : factor = 1.0 / (_beta * _dt * _dt) + _eta[0] * (1.0 + _alpha) * _gamma / _beta / _dt;
484 : else
485 27360 : factor = (*_du_dotdot_du)[0] + _eta[0] * (1.0 + _alpha) * (*_du_dot_du)[0];
486 :
487 1310040 : if (jvar_num == _var.number())
488 218340 : computeJacobian();
489 : else
490 : {
491 : unsigned int coupled_component = 0;
492 : bool disp_coupled = false;
493 : bool rot_coupled = false;
494 :
495 3711780 : for (unsigned int i = 0; i < _ndisp; ++i)
496 : {
497 2947590 : if (jvar_num == _disp_num[i] && _component > 2)
498 : {
499 : coupled_component = i;
500 : disp_coupled = true;
501 : break;
502 : }
503 : }
504 :
505 3275100 : for (unsigned int i = 0; i < _nrot; ++i)
506 : {
507 2729250 : if (jvar_num == _rot_num[i])
508 : {
509 545850 : coupled_component = i + 3;
510 : rot_coupled = true;
511 545850 : break;
512 : }
513 : }
514 :
515 1091700 : prepareMatrixTag(_assembly, _var.number(), jvar_num);
516 :
517 1091700 : if (disp_coupled || rot_coupled)
518 : {
519 2620080 : for (unsigned int i = 0; i < _test.size(); ++i)
520 : {
521 5240160 : for (unsigned int j = 0; j < _phi.size(); ++j)
522 : {
523 3493440 : if (_component < 3 && coupled_component > 2)
524 : {
525 1310040 : RankTwoTensor A;
526 1310040 : A(0, 1) = _Az[0];
527 1310040 : A(0, 2) = -_Ay[0];
528 1310040 : A(1, 0) = -_Az[0];
529 1310040 : A(2, 0) = _Ay[0];
530 :
531 : // conversion from local config to global coordinate system
532 : const RankTwoTensor Ag =
533 1310040 : _original_local_config[0].transpose() * A * _original_local_config[0];
534 :
535 1965060 : _local_ke(i, j) += (i == j ? 1.0 / 3.0 : 1.0 / 6.0) * _density[0] *
536 1310040 : Ag(_component, coupled_component - 3) * _original_length[0] * factor;
537 1310040 : }
538 2183400 : else if (_component > 2 && coupled_component < 3)
539 : {
540 1310040 : RankTwoTensor A;
541 1310040 : A(0, 1) = -_Az[0];
542 1310040 : A(0, 2) = _Ay[0];
543 1310040 : A(1, 0) = _Az[0];
544 1310040 : A(2, 0) = -_Ay[0];
545 :
546 : // conversion from local config to global coordinate system
547 : const RankTwoTensor Ag =
548 1310040 : _original_local_config[0].transpose() * A * _original_local_config[0];
549 :
550 1965060 : _local_ke(i, j) += (i == j ? 1.0 / 3.0 : 1.0 / 6.0) * _density[0] *
551 1310040 : Ag(_component - 3, coupled_component) * _original_length[0] * factor;
552 1310040 : }
553 873360 : else if (_component > 2 && coupled_component > 2)
554 : {
555 873360 : RankTwoTensor I;
556 873360 : if (_has_Ix)
557 0 : I(0, 0) = _Ix[0];
558 : else
559 873360 : I(0, 0) = _Iy[0] + _Iz[0];
560 873360 : I(1, 1) = _Iz[0];
561 873360 : I(2, 2) = _Iy[0];
562 :
563 : // conversion from local config to global coordinate system
564 : const RankTwoTensor Ig =
565 873360 : _original_local_config[0].transpose() * I * _original_local_config[0];
566 :
567 1310040 : _local_ke(i, j) += (i == j ? 1.0 / 3.0 : 1.0 / 6.0) * _density[0] *
568 873360 : Ig(_component - 3, coupled_component - 3) * _original_length[0] *
569 : factor;
570 : }
571 : }
572 : }
573 : }
574 :
575 1091700 : accumulateTaggedLocalMatrix();
576 : }
577 1310040 : }
|
