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