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 354 : InertialForceBeam::validParams()
23 : {
24 354 : InputParameters params = TimeKernel::validParams();
25 354 : 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 354 : params.set<bool>("use_displaced_mesh") = true;
29 708 : params.addRequiredCoupledVar(
30 : "rotations",
31 : "The rotational variables appropriate for the simulation geometry and coordinate system");
32 708 : params.addRequiredCoupledVar(
33 : "displacements",
34 : "The displacement variables appropriate for the simulation geometry and coordinate system");
35 708 : params.addCoupledVar("velocities", "Translational velocity variables");
36 708 : params.addCoupledVar("accelerations", "Translational acceleration variables");
37 708 : params.addCoupledVar("rotational_velocities", "Rotational velocity variables");
38 708 : params.addCoupledVar("rotational_accelerations", "Rotational acceleration variables");
39 708 : params.addRangeCheckedParam<Real>(
40 : "beta", "beta>0.0", "beta parameter for Newmark Time integration");
41 708 : params.addRangeCheckedParam<Real>(
42 : "gamma", "gamma>0.0", "gamma parameter for Newmark Time integration");
43 708 : params.addParam<MaterialPropertyName>("eta",
44 708 : 0.0,
45 : "Name of material property or a constant real "
46 : "number defining the eta parameter for the "
47 : "Rayleigh damping.");
48 1062 : params.addRangeCheckedParam<Real>("alpha",
49 708 : 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 708 : params.addParam<MaterialPropertyName>(
54 : "density",
55 : "density",
56 : "Name of Material Property or a constant real number defining the density of the beam.");
57 708 : params.addRequiredCoupledVar("area", "Variable containing cross-section area");
58 708 : params.addCoupledVar("Ay", 0.0, "Variable containing first moment of area about y axis");
59 708 : params.addCoupledVar("Az", 0.0, "Variable containing first moment of area about z axis");
60 708 : params.addCoupledVar("Ix",
61 : "Variable containing second moment of area about x axis. Defaults to Iy+Iz");
62 708 : params.addRequiredCoupledVar("Iy", "Variable containing second moment of area about y axis");
63 708 : params.addRequiredCoupledVar("Iz", "Variable containing second moment of area about z axis");
64 708 : 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 354 : return params;
71 0 : }
72 :
73 184 : InertialForceBeam::InertialForceBeam(const InputParameters & parameters)
74 : : TimeKernel(parameters),
75 184 : _has_beta(isParamValid("beta")),
76 368 : _has_gamma(isParamValid("gamma")),
77 368 : _has_velocities(isParamValid("velocities")),
78 368 : _has_rot_velocities(isParamValid("rotational_velocities")),
79 368 : _has_accelerations(isParamValid("accelerations")),
80 368 : _has_rot_accelerations(isParamValid("rotational_accelerations")),
81 368 : _has_Ix(isParamValid("Ix")),
82 368 : _density(getMaterialProperty<Real>("density")),
83 184 : _nrot(coupledComponents("rotations")),
84 184 : _ndisp(coupledComponents("displacements")),
85 184 : _rot_num(_nrot),
86 184 : _disp_num(_ndisp),
87 184 : _vel_num(_ndisp),
88 184 : _accel_num(_ndisp),
89 184 : _rot_vel_num(_nrot),
90 184 : _rot_accel_num(_nrot),
91 184 : _area(coupledValue("area")),
92 184 : _Ay(coupledValue("Ay")),
93 184 : _Az(coupledValue("Az")),
94 184 : _Ix(_has_Ix ? coupledValue("Ix") : _zero),
95 184 : _Iy(coupledValue("Iy")),
96 184 : _Iz(coupledValue("Iz")),
97 332 : _beta(_has_beta ? getParam<Real>("beta") : 0.1),
98 332 : _gamma(_has_gamma ? getParam<Real>("gamma") : 0.1),
99 368 : _eta(getMaterialProperty<Real>("eta")),
100 368 : _alpha(getParam<Real>("alpha")),
101 184 : _original_local_config(getMaterialPropertyByName<RankTwoTensor>("initial_rotation")),
102 184 : _original_length(getMaterialPropertyByName<Real>("original_length")),
103 368 : _component(getParam<unsigned int>("component")),
104 184 : _local_force(2),
105 368 : _local_moment(2)
106 : {
107 : // Checking for consistency between the length of the provided rotations and displacements vector
108 184 : if (_ndisp != _nrot)
109 2 : mooseError("InertialForceBeam: The number of variables supplied in 'displacements' and "
110 : "'rotations' must match.");
111 :
112 182 : if (_has_beta && _has_gamma && _has_velocities && _has_accelerations && _has_rot_velocities &&
113 146 : _has_rot_accelerations)
114 : {
115 292 : if ((coupledComponents("velocities") != _ndisp) ||
116 292 : (coupledComponents("accelerations") != _ndisp) ||
117 436 : (coupledComponents("rotational_velocities") != _ndisp) ||
118 290 : (coupledComponents("rotational_accelerations") != _ndisp))
119 2 : 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 576 : for (unsigned int i = 0; i < _ndisp; ++i)
126 : {
127 864 : MooseVariable * vel_variable = getVar("velocities", i);
128 432 : _vel_num[i] = vel_variable->number();
129 :
130 864 : MooseVariable * accel_variable = getVar("accelerations", i);
131 432 : _accel_num[i] = accel_variable->number();
132 :
133 864 : MooseVariable * rot_vel_variable = getVar("rotational_velocities", i);
134 432 : _rot_vel_num[i] = rot_vel_variable->number();
135 :
136 864 : MooseVariable * rot_accel_variable = getVar("rotational_accelerations", i);
137 432 : _rot_accel_num[i] = rot_accel_variable->number();
138 : }
139 : }
140 36 : else if (!_has_beta && !_has_gamma && !_has_velocities && !_has_accelerations &&
141 36 : !_has_rot_velocities && !_has_rot_accelerations)
142 : {
143 36 : _du_dot_du = &coupledDotDu("displacements", 0);
144 36 : _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 720 : for (unsigned int i = 0; i < _ndisp; ++i)
153 : {
154 1080 : MooseVariable * disp_variable = getVar("displacements", i);
155 540 : _disp_num[i] = disp_variable->number();
156 :
157 1080 : MooseVariable * rot_variable = getVar("rotations", i);
158 540 : _rot_num[i] = rot_variable->number();
159 : }
160 180 : }
161 :
162 : void
163 576000 : InertialForceBeam::computeResidual()
164 : {
165 576000 : prepareVectorTag(_assembly, _var.number());
166 :
167 576000 : if (_dt != 0.0)
168 : {
169 : // fetch the two end nodes for _current_elem
170 : std::vector<const Node *> node;
171 1728000 : for (unsigned int i = 0; i < 2; ++i)
172 1152000 : node.push_back(_current_elem->node_ptr(i));
173 :
174 : // Fetch the solution for the two end nodes at time t
175 576000 : NonlinearSystemBase & nonlinear_sys = _fe_problem.getNonlinearSystemBase(_sys.number());
176 :
177 576000 : if (_has_beta)
178 : {
179 564000 : const NumericVector<Number> & sol = *nonlinear_sys.currentSolution();
180 564000 : const NumericVector<Number> & sol_old = nonlinear_sys.solutionOld();
181 :
182 564000 : AuxiliarySystem & aux = _fe_problem.getAuxiliarySystem();
183 564000 : 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 2256000 : for (unsigned int i = 0; i < _ndisp; ++i)
189 : {
190 : // obtain delta displacement
191 1692000 : unsigned int dof_index_0 = node[0]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
192 1692000 : unsigned int dof_index_1 = node[1]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
193 1692000 : const Real disp_0 = sol(dof_index_0) - sol_old(dof_index_0);
194 1692000 : const Real disp_1 = sol(dof_index_1) - sol_old(dof_index_1);
195 :
196 1692000 : dof_index_0 = node[0]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
197 1692000 : dof_index_1 = node[1]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
198 1692000 : const Real rot_0 = sol(dof_index_0) - sol_old(dof_index_0);
199 1692000 : 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 1692000 : _vel_old_0(i) = aux_sol_old(node[0]->dof_number(aux.number(), _vel_num[i], 0));
204 1692000 : _vel_old_1(i) = aux_sol_old(node[1]->dof_number(aux.number(), _vel_num[i], 0));
205 1692000 : const Real accel_old_0 = aux_sol_old(node[0]->dof_number(aux.number(), _accel_num[i], 0));
206 1692000 : const Real accel_old_1 = aux_sol_old(node[1]->dof_number(aux.number(), _accel_num[i], 0));
207 :
208 1692000 : _rot_vel_old_0(i) = aux_sol_old(node[0]->dof_number(aux.number(), _rot_vel_num[i], 0));
209 1692000 : _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 1692000 : aux_sol_old(node[0]->dof_number(aux.number(), _rot_accel_num[i], 0));
212 : const Real rot_accel_old_1 =
213 1692000 : aux_sol_old(node[1]->dof_number(aux.number(), _rot_accel_num[i], 0));
214 :
215 1692000 : _accel_0(i) =
216 1692000 : 1. / _beta * (disp_0 / (_dt * _dt) - _vel_old_0(i) / _dt - accel_old_0 * (0.5 - _beta));
217 1692000 : _accel_1(i) =
218 1692000 : 1. / _beta * (disp_1 / (_dt * _dt) - _vel_old_1(i) / _dt - accel_old_1 * (0.5 - _beta));
219 1692000 : _rot_accel_0(i) =
220 1692000 : 1. / _beta *
221 1692000 : (rot_0 / (_dt * _dt) - _rot_vel_old_0(i) / _dt - rot_accel_old_0 * (0.5 - _beta));
222 1692000 : _rot_accel_1(i) =
223 1692000 : 1. / _beta *
224 1692000 : (rot_1 / (_dt * _dt) - _rot_vel_old_1(i) / _dt - rot_accel_old_1 * (0.5 - _beta));
225 :
226 1692000 : _vel_0(i) = _vel_old_0(i) + (_dt * (1 - _gamma)) * accel_old_0 + _gamma * _dt * _accel_0(i);
227 1692000 : _vel_1(i) = _vel_old_1(i) + (_dt * (1 - _gamma)) * accel_old_1 + _gamma * _dt * _accel_1(i);
228 1692000 : _rot_vel_0(i) = _rot_vel_old_0(i) + (_dt * (1 - _gamma)) * rot_accel_old_0 +
229 1692000 : _gamma * _dt * _rot_accel_0(i);
230 1692000 : _rot_vel_1(i) = _rot_vel_old_1(i) + (_dt * (1 - _gamma)) * rot_accel_old_1 +
231 1692000 : _gamma * _dt * _rot_accel_1(i);
232 : }
233 : }
234 : else
235 : {
236 12000 : 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 12000 : 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 12000 : 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 12000 : const NumericVector<Number> & vel = *nonlinear_sys.solutionUDot();
251 12000 : const NumericVector<Number> & vel_old = *nonlinear_sys.solutionUDotOld();
252 12000 : const NumericVector<Number> & accel = *nonlinear_sys.solutionUDotDot();
253 :
254 48000 : for (unsigned int i = 0; i < _ndisp; ++i)
255 : {
256 : // translational velocities and accelerations
257 36000 : unsigned int dof_index_0 = node[0]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
258 36000 : unsigned int dof_index_1 = node[1]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
259 36000 : _vel_0(i) = vel(dof_index_0);
260 36000 : _vel_1(i) = vel(dof_index_1);
261 36000 : _vel_old_0(i) = vel_old(dof_index_0);
262 36000 : _vel_old_1(i) = vel_old(dof_index_1);
263 36000 : _accel_0(i) = accel(dof_index_0);
264 36000 : _accel_1(i) = accel(dof_index_1);
265 :
266 : // rotational velocities and accelerations
267 36000 : dof_index_0 = node[0]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
268 36000 : dof_index_1 = node[1]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
269 36000 : _rot_vel_0(i) = vel(dof_index_0);
270 36000 : _rot_vel_1(i) = vel(dof_index_1);
271 36000 : _rot_vel_old_0(i) = vel_old(dof_index_0);
272 36000 : _rot_vel_old_1(i) = vel_old(dof_index_1);
273 36000 : _rot_accel_0(i) = accel(dof_index_0);
274 36000 : _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 576000 : _local_vel_old_0 = _original_local_config[0] * _vel_old_0;
281 576000 : _local_vel_old_1 = _original_local_config[0] * _vel_old_1;
282 576000 : _local_vel_0 = _original_local_config[0] * _vel_0;
283 576000 : _local_vel_1 = _original_local_config[0] * _vel_1;
284 576000 : _local_accel_0 = _original_local_config[0] * _accel_0;
285 576000 : _local_accel_1 = _original_local_config[0] * _accel_1;
286 :
287 576000 : _local_rot_vel_old_0 = _original_local_config[0] * _rot_vel_old_0;
288 576000 : _local_rot_vel_old_1 = _original_local_config[0] * _rot_vel_old_1;
289 576000 : _local_rot_vel_0 = _original_local_config[0] * _rot_vel_0;
290 576000 : _local_rot_vel_1 = _original_local_config[0] * _rot_vel_1;
291 576000 : _local_rot_accel_0 = _original_local_config[0] * _rot_accel_0;
292 576000 : _local_rot_accel_1 = _original_local_config[0] * _rot_accel_1;
293 :
294 : // local residual
295 2304000 : for (unsigned int i = 0; i < _ndisp; ++i)
296 : {
297 1728000 : if (_component < 3)
298 : {
299 864000 : _local_force[0](i) = _density[0] * _area[0] * _original_length[0] / 3.0 *
300 864000 : (_local_accel_0(i) + _local_accel_1(i) / 2.0 +
301 864000 : _eta[0] * (1.0 + _alpha) * (_local_vel_0(i) + _local_vel_1(i) / 2.0) -
302 864000 : _alpha * _eta[0] * (_local_vel_old_0(i) + _local_vel_old_1(i) / 2.0));
303 864000 : _local_force[1](i) = _density[0] * _area[0] * _original_length[0] / 3.0 *
304 864000 : (_local_accel_1(i) + _local_accel_0(i) / 2.0 +
305 864000 : _eta[0] * (1.0 + _alpha) * (_local_vel_1(i) + _local_vel_0(i) / 2.0) -
306 864000 : _alpha * _eta[0] * (_local_vel_old_1(i) + _local_vel_old_0(i) / 2.0));
307 : }
308 :
309 1728000 : if (_component > 2)
310 : {
311 864000 : Real I = _Iy[0] + _Iz[0];
312 864000 : if (_has_Ix && (i == 0))
313 0 : I = _Ix[0];
314 864000 : if (i == 1)
315 : I = _Iz[0];
316 576000 : else if (i == 2)
317 : I = _Iy[0];
318 :
319 864000 : _local_moment[0](i) =
320 864000 : _density[0] * I * _original_length[0] / 3.0 *
321 864000 : (_local_rot_accel_0(i) + _local_rot_accel_1(i) / 2.0 +
322 864000 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_0(i) + _local_rot_vel_1(i) / 2.0) -
323 864000 : _alpha * _eta[0] * (_local_rot_vel_old_0(i) + _local_rot_vel_old_1(i) / 2.0));
324 864000 : _local_moment[1](i) =
325 864000 : _density[0] * I * _original_length[0] / 3.0 *
326 864000 : (_local_rot_accel_1(i) + _local_rot_accel_0(i) / 2.0 +
327 864000 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_1(i) + _local_rot_vel_0(i) / 2.0) -
328 864000 : _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 576000 : if (_component < 3)
335 : {
336 288000 : _local_force[0](0) +=
337 288000 : _density[0] * _original_length[0] / 3.0 *
338 288000 : (_Az[0] * (_local_rot_accel_0(1) + _local_rot_accel_1(1) / 2.0 +
339 288000 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_0(1) + _local_rot_vel_1(1) / 2.0) -
340 288000 : _alpha * _eta[0] * (_local_rot_vel_old_0(1) + _local_rot_vel_old_1(1) / 2.0)) -
341 288000 : _Ay[0] * (_local_rot_accel_0(2) + _local_rot_accel_1(2) / 2.0 +
342 288000 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_0(2) + _local_rot_vel_1(2) / 2.0) -
343 288000 : _alpha * _eta[0] * (_local_rot_vel_old_0(2) + _local_rot_vel_old_1(2) / 2.0)));
344 288000 : _local_force[1](0) +=
345 288000 : _density[0] * _original_length[0] / 3.0 *
346 288000 : (_Az[0] * (_local_rot_accel_1(1) + _local_rot_accel_0(1) / 2.0 +
347 288000 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_1(1) + _local_rot_vel_0(1) / 2.0) -
348 288000 : _alpha * _eta[0] * (_local_rot_vel_old_1(1) + _local_rot_vel_old_0(1) / 2.0)) -
349 288000 : _Ay[0] * (_local_rot_accel_1(2) + _local_rot_accel_0(2) / 2.0 +
350 288000 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_1(2) + _local_rot_vel_0(2) / 2.0) -
351 288000 : _alpha * _eta[0] * (_local_rot_vel_old_1(2) + _local_rot_vel_old_0(2) / 2.0)));
352 :
353 288000 : _local_force[0](1) +=
354 288000 : -_density[0] * _original_length[0] / 3.0 * _Az[0] *
355 288000 : (_local_rot_accel_0(0) + _local_rot_accel_1(0) / 2.0 +
356 288000 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_0(0) + _local_rot_vel_1(0) / 2.0) -
357 288000 : _alpha * _eta[0] * (_local_rot_vel_old_0(0) + _local_rot_vel_old_1(0) / 2.0));
358 288000 : _local_force[1](1) +=
359 288000 : -_density[0] * _original_length[0] / 3.0 * _Az[0] *
360 288000 : (_local_rot_accel_1(0) + _local_rot_accel_0(0) / 2.0 +
361 288000 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_1(0) + _local_rot_vel_0(0) / 2.0) -
362 288000 : _alpha * _eta[0] * (_local_rot_vel_old_1(0) + _local_rot_vel_old_0(0) / 2.0));
363 :
364 288000 : _local_force[0](2) +=
365 288000 : _density[0] * _original_length[0] / 3.0 * _Ay[0] *
366 288000 : (_local_rot_accel_0(0) + _local_rot_accel_1(0) / 2.0 +
367 288000 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_0(0) + _local_rot_vel_1(0) / 2.0) -
368 288000 : _alpha * _eta[0] * (_local_rot_vel_old_0(0) + _local_rot_vel_old_1(0) / 2.0));
369 288000 : _local_force[1](2) +=
370 288000 : _density[0] * _original_length[0] / 3.0 * _Ay[0] *
371 288000 : (_local_rot_accel_1(0) + _local_rot_accel_0(0) / 2.0 +
372 288000 : _eta[0] * (1.0 + _alpha) * (_local_rot_vel_1(0) + _local_rot_vel_0(0) / 2.0) -
373 288000 : _alpha * _eta[0] * (_local_rot_vel_old_1(0) + _local_rot_vel_old_0(0) / 2.0));
374 : }
375 : else
376 : {
377 288000 : _local_moment[0](0) += _density[0] * _original_length[0] / 3.0 *
378 288000 : (-_Az[0] * (_local_accel_0(1) + _local_accel_1(1) / 2.0) +
379 288000 : _Ay[0] * (_local_accel_0(1) + _local_accel_1(1) / 2.0));
380 288000 : _local_moment[1](0) += _density[0] * _original_length[0] / 3.0 *
381 288000 : (-_Az[0] * (_local_accel_1(1) + _local_accel_0(1) / 2.0) +
382 288000 : _Ay[0] * (_local_accel_1(1) + _local_accel_0(1) / 2.0));
383 :
384 288000 : _local_moment[0](1) += _density[0] * _original_length[0] / 3.0 * _Az[0] *
385 288000 : (_local_accel_0(0) + _local_accel_1(0) / 2.0);
386 288000 : _local_moment[1](1) += _density[0] * _original_length[0] / 3.0 * _Az[0] *
387 288000 : (_local_accel_1(0) + _local_accel_0(0) / 2.0);
388 :
389 288000 : _local_moment[0](2) += -_density[0] * _original_length[0] / 3.0 * _Ay[0] *
390 288000 : (_local_accel_0(0) + _local_accel_1(0) / 2.0);
391 288000 : _local_moment[1](2) += -_density[0] * _original_length[0] / 3.0 * _Ay[0] *
392 288000 : (_local_accel_1(0) + _local_accel_0(0) / 2.0);
393 : }
394 :
395 : // Global force and moments
396 576000 : if (_component < 3)
397 : {
398 288000 : _global_force_0 = _original_local_config[0] * _local_force[0];
399 288000 : _global_force_1 = _original_local_config[0] * _local_force[1];
400 288000 : _local_re(0) = _global_force_0(_component);
401 288000 : _local_re(1) = _global_force_1(_component);
402 : }
403 : else
404 : {
405 288000 : _global_moment_0 = _original_local_config[0] * _local_moment[0];
406 288000 : _global_moment_1 = _original_local_config[0] * _local_moment[1];
407 288000 : _local_re(0) = _global_moment_0(_component - 3);
408 288000 : _local_re(1) = _global_moment_1(_component - 3);
409 : }
410 : }
411 :
412 576000 : accumulateTaggedLocalResidual();
413 :
414 576000 : 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 576000 : }
421 :
422 : void
423 294000 : InertialForceBeam::computeJacobian()
424 : {
425 294000 : 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 294000 : if (_has_beta)
431 288000 : factor = 1.0 / (_beta * _dt * _dt) + _eta[0] * (1.0 + _alpha) * _gamma / _beta / _dt;
432 : else
433 6000 : factor = (*_du_dotdot_du)[0] + _eta[0] * (1.0 + _alpha) * (*_du_dot_du)[0];
434 :
435 882000 : for (unsigned int i = 0; i < _test.size(); ++i)
436 : {
437 1764000 : for (unsigned int j = 0; j < _phi.size(); ++j)
438 : {
439 1176000 : if (_component < 3)
440 882000 : _local_ke(i, j) = (i == j ? 1.0 / 3.0 : 1.0 / 6.0) * _density[0] * _area[0] *
441 588000 : _original_length[0] * factor;
442 : else if (_component > 2)
443 : {
444 588000 : RankTwoTensor I;
445 588000 : if (_has_Ix)
446 0 : I(0, 0) = _Ix[0];
447 : else
448 588000 : I(0, 0) = _Iy[0] + _Iz[0];
449 588000 : I(1, 1) = _Iz[0];
450 588000 : I(2, 2) = _Iy[0];
451 :
452 : // conversion from local config to global coordinate system
453 588000 : RankTwoTensor Ig = _original_local_config[0].transpose() * I * _original_local_config[0];
454 :
455 882000 : _local_ke(i, j) = (i == j ? 1.0 / 3.0 : 1.0 / 6.0) * _density[0] *
456 588000 : Ig(_component - 3, _component - 3) * _original_length[0] * factor;
457 : }
458 : }
459 : }
460 :
461 294000 : accumulateTaggedLocalMatrix();
462 :
463 294000 : 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 294000 : }
475 :
476 : void
477 1764000 : 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 1764000 : if (_has_beta)
483 1728000 : factor = 1.0 / (_beta * _dt * _dt) + _eta[0] * (1.0 + _alpha) * _gamma / _beta / _dt;
484 : else
485 36000 : factor = (*_du_dotdot_du)[0] + _eta[0] * (1.0 + _alpha) * (*_du_dot_du)[0];
486 :
487 1764000 : if (jvar_num == _var.number())
488 294000 : computeJacobian();
489 : else
490 : {
491 : unsigned int coupled_component = 0;
492 : bool disp_coupled = false;
493 : bool rot_coupled = false;
494 :
495 4998000 : for (unsigned int i = 0; i < _ndisp; ++i)
496 : {
497 3969000 : if (jvar_num == _disp_num[i] && _component > 2)
498 : {
499 : coupled_component = i;
500 : disp_coupled = true;
501 : break;
502 : }
503 : }
504 :
505 4410000 : for (unsigned int i = 0; i < _nrot; ++i)
506 : {
507 3675000 : if (jvar_num == _rot_num[i])
508 : {
509 735000 : coupled_component = i + 3;
510 : rot_coupled = true;
511 735000 : break;
512 : }
513 : }
514 :
515 1470000 : prepareMatrixTag(_assembly, _var.number(), jvar_num);
516 :
517 1470000 : if (disp_coupled || rot_coupled)
518 : {
519 3528000 : for (unsigned int i = 0; i < _test.size(); ++i)
520 : {
521 7056000 : for (unsigned int j = 0; j < _phi.size(); ++j)
522 : {
523 4704000 : if (_component < 3 && coupled_component > 2)
524 : {
525 1764000 : RankTwoTensor A;
526 1764000 : A(0, 1) = _Az[0];
527 1764000 : A(0, 2) = -_Ay[0];
528 1764000 : A(1, 0) = -_Az[0];
529 1764000 : A(2, 0) = _Ay[0];
530 :
531 : // conversion from local config to global coordinate system
532 : const RankTwoTensor Ag =
533 1764000 : _original_local_config[0].transpose() * A * _original_local_config[0];
534 :
535 2646000 : _local_ke(i, j) += (i == j ? 1.0 / 3.0 : 1.0 / 6.0) * _density[0] *
536 1764000 : Ag(_component, coupled_component - 3) * _original_length[0] * factor;
537 1764000 : }
538 2940000 : else if (_component > 2 && coupled_component < 3)
539 : {
540 1764000 : RankTwoTensor A;
541 1764000 : A(0, 1) = -_Az[0];
542 1764000 : A(0, 2) = _Ay[0];
543 1764000 : A(1, 0) = _Az[0];
544 1764000 : A(2, 0) = -_Ay[0];
545 :
546 : // conversion from local config to global coordinate system
547 : const RankTwoTensor Ag =
548 1764000 : _original_local_config[0].transpose() * A * _original_local_config[0];
549 :
550 2646000 : _local_ke(i, j) += (i == j ? 1.0 / 3.0 : 1.0 / 6.0) * _density[0] *
551 1764000 : Ag(_component - 3, coupled_component) * _original_length[0] * factor;
552 1764000 : }
553 1176000 : else if (_component > 2 && coupled_component > 2)
554 : {
555 1176000 : RankTwoTensor I;
556 1176000 : if (_has_Ix)
557 0 : I(0, 0) = _Ix[0];
558 : else
559 1176000 : I(0, 0) = _Iy[0] + _Iz[0];
560 1176000 : I(1, 1) = _Iz[0];
561 1176000 : I(2, 2) = _Iy[0];
562 :
563 : // conversion from local config to global coordinate system
564 : const RankTwoTensor Ig =
565 1176000 : _original_local_config[0].transpose() * I * _original_local_config[0];
566 :
567 1764000 : _local_ke(i, j) += (i == j ? 1.0 / 3.0 : 1.0 / 6.0) * _density[0] *
568 1176000 : Ig(_component - 3, coupled_component - 3) * _original_length[0] *
569 : factor;
570 : }
571 : }
572 : }
573 : }
574 :
575 1470000 : accumulateTaggedLocalMatrix();
576 : }
577 1764000 : }
|
