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