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 "NodalRotationalInertia.h"
11 : #include "MooseVariable.h"
12 : #include "AuxiliarySystem.h"
13 : #include "MooseMesh.h"
14 : #include "TimeIntegrator.h"
15 : #include "FEProblemBase.h"
16 :
17 : registerMooseObject("SolidMechanicsApp", NodalRotationalInertia);
18 :
19 : InputParameters
20 126 : NodalRotationalInertia::validParams()
21 : {
22 126 : InputParameters params = TimeNodalKernel::validParams();
23 126 : params.addClassDescription("Calculates the inertial torques and inertia proportional damping "
24 : "corresponding to the nodal rotational inertia.");
25 252 : params.addCoupledVar("rotations", "rotational displacement variables");
26 252 : params.addCoupledVar("rotational_velocities", "rotational velocity variables");
27 252 : params.addCoupledVar("rotational_accelerations", "rotational acceleration variables");
28 252 : params.addRangeCheckedParam<Real>(
29 : "beta", "beta>0.0", "beta parameter for Newmark Time integration");
30 252 : params.addRangeCheckedParam<Real>(
31 : "gamma", "gamma>0.0", "gamma parameter for Newmark Time integration");
32 378 : params.addRangeCheckedParam<Real>("eta",
33 252 : 0.0,
34 : "eta>=0.0",
35 : "Constant real number defining the eta parameter for"
36 : "Rayleigh damping.");
37 378 : params.addRangeCheckedParam<Real>("alpha",
38 252 : 0.0,
39 : "alpha >= -0.3333 & alpha <= 0.0",
40 : "Alpha parameter for mass dependent numerical damping induced "
41 : "by HHT time integration scheme");
42 252 : params.addRequiredRangeCheckedParam<Real>(
43 : "Ixx", "Ixx>0.0", "Moment of inertia in the x direction.");
44 252 : params.addRequiredRangeCheckedParam<Real>(
45 : "Iyy", "Iyy>0.0", "Moment of inertia in the y direction.");
46 252 : params.addRequiredRangeCheckedParam<Real>(
47 : "Izz", "Izz>0.0", "Moment of inertia in the z direction.");
48 252 : params.addParam<Real>("Ixy", 0.0, "Moment of inertia in the xy direction.");
49 252 : params.addParam<Real>("Ixz", 0.0, "Moment of inertia in the xz direction.");
50 252 : params.addParam<Real>("Iyz", 0.0, "Moment of inertia in the yz direction.");
51 252 : params.addParam<RealGradient>(
52 : "x_orientation", "Unit vector along the x direction if different from global x direction.");
53 252 : params.addParam<RealGradient>(
54 : "y_orientation", "Unit vector along the y direction if different from global y direction.");
55 252 : params.addRequiredRangeCheckedParam<unsigned int>(
56 : "component",
57 : "component<3",
58 : "An integer corresponding to the direction "
59 : "the variable this nodalkernel acts in. (0 for rot_x, "
60 : "1 for rot_y, and 2 for rot_z).");
61 126 : return params;
62 0 : }
63 :
64 64 : NodalRotationalInertia::NodalRotationalInertia(const InputParameters & parameters)
65 : : TimeNodalKernel(parameters),
66 64 : _has_beta(isParamValid("beta")),
67 128 : _has_gamma(isParamValid("gamma")),
68 128 : _has_rot_velocities(isParamValid("rotational_velocities")),
69 128 : _has_rot_accelerations(isParamValid("rotational_accelerations")),
70 128 : _has_x_orientation(isParamValid("x_orientation")),
71 128 : _has_y_orientation(isParamValid("y_orientation")),
72 64 : _nrot(coupledComponents("rotations")),
73 64 : _rot(_nrot),
74 64 : _rot_old(_nrot),
75 64 : _rot_vel_num(_nrot),
76 64 : _rot_accel_num(_nrot),
77 64 : _rot_variables(_nrot),
78 64 : _rot_accel(_nrot),
79 64 : _rot_vel(_nrot),
80 64 : _rot_vel_old(_nrot),
81 110 : _beta(_has_beta ? getParam<Real>("beta") : 0.1),
82 110 : _gamma(_has_gamma ? getParam<Real>("gamma") : 0.1),
83 128 : _eta(getParam<Real>("eta")),
84 128 : _alpha(getParam<Real>("alpha")),
85 128 : _component(getParam<unsigned int>("component")),
86 64 : _rot_dot_residual(_nrot),
87 64 : _rot_vel_old_value(_nrot),
88 64 : _rot_dotdot_residual(_nrot),
89 192 : _time_integrator(_sys.getTimeIntegrator(_var.number()))
90 : {
91 64 : if (_has_beta && _has_gamma && _has_rot_velocities && _has_rot_accelerations)
92 : {
93 46 : _aux_sys = &(_fe_problem.getAuxiliarySystem());
94 92 : if (coupledComponents("rotational_velocities") != _nrot ||
95 90 : coupledComponents("rotational_accelerations") != _nrot)
96 2 : mooseError(
97 : "NodalRotationalInertia: rotational_velocities and rotational_accelerations should "
98 : "be same size "
99 : "as rotations.");
100 :
101 176 : for (unsigned int i = 0; i < _nrot; ++i)
102 : {
103 132 : MooseVariable * rot_var = getVar("rotations", i);
104 132 : MooseVariable * rot_vel_var = getVar("rotational_velocities", i);
105 132 : MooseVariable * rot_accel_var = getVar("rotational_accelerations", i);
106 :
107 132 : _rot[i] = &rot_var->dofValues();
108 132 : _rot_old[i] = &rot_var->dofValuesOld();
109 :
110 132 : _rot_vel_num[i] = rot_vel_var->number();
111 132 : _rot_accel_num[i] = rot_accel_var->number();
112 :
113 132 : _rot_variables[i] = coupled("rotations", i);
114 : }
115 : }
116 18 : else if (!_has_beta && !_has_gamma && !_has_rot_velocities && !_has_rot_accelerations)
117 : {
118 72 : for (unsigned int i = 0; i < _nrot; ++i)
119 : {
120 54 : MooseVariable * rot_var = getVar("rotations", i);
121 54 : _rot_vel_old_value[i] = &rot_var->dofValuesDotOld();
122 54 : _rot_dot_residual[i] =
123 54 : &coupledVectorTagDofValue("rotations", _time_integrator.uDotFactorTag(), i);
124 54 : _rot_dotdot_residual[i] =
125 54 : &coupledVectorTagDofValue("rotations", _time_integrator.uDotDotFactorTag(), i);
126 :
127 54 : if (i == 0)
128 : {
129 18 : _du_dot_du = &rot_var->dofValuesDuDotDu();
130 18 : _du_dotdot_du = &rot_var->dofValuesDuDotDotDu();
131 : }
132 : }
133 : }
134 : else
135 0 : mooseError("NodalRotationalInertia: Either all or none of `beta`, `gamma`, "
136 : "`rotational_velocities` and `rotational_accelerations` should be provided as "
137 : "input.");
138 :
139 : // Store inertia values in inertia tensor
140 : _inertia.zero();
141 124 : _inertia(0, 0) = getParam<Real>("Ixx");
142 124 : _inertia(0, 1) = getParam<Real>("Ixy");
143 124 : _inertia(0, 2) = getParam<Real>("Ixz");
144 62 : _inertia(1, 0) = _inertia(0, 1);
145 124 : _inertia(1, 1) = getParam<Real>("Iyy");
146 124 : _inertia(1, 2) = getParam<Real>("Iyz");
147 62 : _inertia(2, 0) = _inertia(0, 2);
148 62 : _inertia(2, 1) = _inertia(1, 2);
149 124 : _inertia(2, 2) = getParam<Real>("Izz");
150 :
151 : // Check for positive definiteness of matrix using Sylvester's criterion.
152 62 : const Real det1 = _inertia(0, 0);
153 62 : const Real det2 = _inertia(0, 0) * _inertia(1, 1) - _inertia(0, 1) * _inertia(1, 0);
154 62 : const Real det3 = _inertia.det();
155 :
156 62 : if (det1 < 1e-6 || det2 < 1e-6 || det3 < 1e-6)
157 2 : mooseError("NodalRotationalInertia: The moment of inertia tensor should be positive definite.");
158 60 : if (_has_x_orientation && _has_y_orientation)
159 : {
160 8 : const RealGradient x_orientation = getParam<RealGradient>("x_orientation");
161 8 : const RealGradient y_orientation = getParam<RealGradient>("y_orientation");
162 :
163 4 : if (std::abs(x_orientation.norm() - 1.0) > 1e-6 || std::abs(y_orientation.norm() - 1.0) > 1e-6)
164 2 : mooseError("NodalRotationalInertia: x_orientation and y_orientation must be unit vectors.");
165 :
166 2 : Real sum = x_orientation(0) * y_orientation(0) + x_orientation(1) * y_orientation(1) +
167 2 : x_orientation(2) * y_orientation(2);
168 :
169 2 : if (std::abs(sum) > 1e-4)
170 2 : mooseError("NodalRotationalInertia: x_orientation and y_orientation should be perpendicular "
171 : "to each other.");
172 :
173 : // Calculate z orientation as a cross product of the x and y orientations
174 : RealGradient z_orientation;
175 0 : z_orientation(0) = (x_orientation(1) * y_orientation(2) - x_orientation(2) * y_orientation(1));
176 0 : z_orientation(1) = (x_orientation(2) * y_orientation(0) - x_orientation(0) * y_orientation(2));
177 0 : z_orientation(2) = (x_orientation(0) * y_orientation(1) - x_orientation(1) * y_orientation(0));
178 :
179 : // Rotation matrix from global to original beam local configuration
180 0 : RankTwoTensor orientation;
181 0 : orientation(0, 0) = x_orientation(0);
182 0 : orientation(0, 1) = x_orientation(1);
183 0 : orientation(0, 2) = x_orientation(2);
184 0 : orientation(1, 0) = y_orientation(0);
185 0 : orientation(1, 1) = y_orientation(1);
186 0 : orientation(1, 2) = y_orientation(2);
187 0 : orientation(2, 0) = z_orientation(0);
188 0 : orientation(2, 1) = z_orientation(1);
189 0 : orientation(2, 2) = z_orientation(2);
190 :
191 0 : RankTwoTensor global_inertia = orientation.transpose() * _inertia * orientation;
192 0 : _inertia = global_inertia;
193 0 : }
194 56 : else if ((_has_x_orientation && !_has_y_orientation) ||
195 54 : (!_has_x_orientation && _has_y_orientation))
196 2 : mooseError("NodalRotationalInertia: Both x_orientation and y_orientation should be provided if "
197 : "x_orientation or "
198 : "y_orientation is different from global x or y direction, respectively.");
199 :
200 : // Check for Explicit and alpha parameter
201 54 : if (_alpha != 0 && _time_integrator.isExplicit())
202 0 : mooseError("NodalRotationalInertia: HHT time integration parameter can only be used with "
203 : "Newmark-Beta time integration.");
204 :
205 : // Check for Explicit and beta parameter
206 54 : if (_has_beta != 0 && _time_integrator.isExplicit())
207 0 : mooseError("NodalRotationalInertia: beta time integration parameter can only be used with "
208 : "Newmark-Beta time integrator.");
209 54 : }
210 :
211 : Real
212 3606 : NodalRotationalInertia::computeQpResidual()
213 : {
214 3606 : if (_dt == 0)
215 : return 0;
216 : else
217 : {
218 3606 : if (_has_beta)
219 : {
220 2400 : const NumericVector<Number> & aux_sol_old = _aux_sys->solutionOld();
221 :
222 9600 : for (unsigned int i = 0; i < _nrot; ++i)
223 : {
224 7200 : _rot_vel_old[i] =
225 7200 : aux_sol_old(_current_node->dof_number(_aux_sys->number(), _rot_vel_num[i], 0));
226 : const Real rot_accel_old =
227 7200 : aux_sol_old(_current_node->dof_number(_aux_sys->number(), _rot_accel_num[i], 0));
228 :
229 7200 : _rot_accel[i] = 1.0 / _beta *
230 7200 : ((((*_rot[i])[_qp] - (*_rot_old[i])[_qp]) / (_dt * _dt)) -
231 7200 : _rot_vel_old[i] / _dt - rot_accel_old * (0.5 - _beta));
232 7200 : _rot_vel[i] =
233 7200 : _rot_vel_old[i] + (_dt * (1.0 - _gamma)) * rot_accel_old + _gamma * _dt * _rot_accel[i];
234 : }
235 :
236 : Real res = 0.0;
237 9600 : for (unsigned int i = 0; i < _nrot; ++i)
238 7200 : res += _inertia(_component, i) * (_rot_accel[i] + _rot_vel[i] * _eta * (1.0 + _alpha) -
239 7200 : _alpha * _eta * _rot_vel_old[i]);
240 :
241 2400 : return res;
242 : }
243 : else
244 : {
245 : // All cases (Explicit, implicit and implicit with HHT)
246 : // Note that _alpha is ensured to be zero with explicit integration
247 : Real res = 0.0;
248 4824 : for (unsigned int i = 0; i < _nrot; ++i)
249 3618 : res += _inertia(_component, i) * ((*_rot_dotdot_residual[i])[_qp] +
250 3618 : (*_rot_dot_residual[i])[_qp] * _eta * (1.0 + _alpha) -
251 3618 : _alpha * _eta * (*_rot_vel_old_value[i])[_qp]);
252 1206 : return res;
253 : }
254 : }
255 : }
256 :
257 : Real
258 1806 : NodalRotationalInertia::computeQpJacobian()
259 : {
260 1806 : if (_dt == 0)
261 : return 0.0;
262 : else
263 : {
264 1806 : if (_has_beta)
265 1200 : return _inertia(_component, _component) / (_beta * _dt * _dt) +
266 1200 : _eta * (1.0 + _alpha) * _inertia(_component, _component) * _gamma / _beta / _dt;
267 606 : else if (_time_integrator.isExplicit())
268 : // for explicit central difference integration, _eta does not appear in the
269 : // Jacobian (mass matrix), and alpha is zero
270 0 : return _inertia(_component, _component) * (*_du_dotdot_du)[_qp];
271 : else
272 : // for NewmarkBeta time integrator
273 606 : return _inertia(_component, _component) * (*_du_dotdot_du)[_qp] +
274 606 : _eta * (1.0 + _alpha) * _inertia(_component, _component) * (*_du_dot_du)[_qp];
275 : }
276 : }
277 :
278 : Real
279 3612 : NodalRotationalInertia::computeQpOffDiagJacobian(unsigned int jvar)
280 : {
281 : unsigned int coupled_component = 0;
282 : bool rot_coupled = false;
283 :
284 9648 : for (unsigned int i = 0; i < _nrot; ++i)
285 : {
286 8436 : if (jvar == _rot_variables[i])
287 : {
288 : coupled_component = i;
289 : rot_coupled = true;
290 : break;
291 : }
292 : }
293 :
294 3612 : if (_dt == 0)
295 : return 0.0;
296 3612 : else if (rot_coupled)
297 : {
298 2400 : if (_has_beta)
299 2400 : return _inertia(_component, coupled_component) / (_beta * _dt * _dt) +
300 2400 : _eta * (1.0 + _alpha) * _inertia(_component, coupled_component) * _gamma / _beta / _dt;
301 : else
302 0 : return _inertia(_component, coupled_component) * (*_du_dotdot_du)[_qp] +
303 0 : _eta * (1.0 + _alpha) * _inertia(_component, coupled_component) * (*_du_dot_du)[_qp];
304 : }
305 : else
306 : return 0.0;
307 : }
|
