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 "NodalTranslationalInertia.h"
11 : #include "MooseVariable.h"
12 : #include "AuxiliarySystem.h"
13 : #include "MooseUtils.h"
14 : #include "DelimitedFileReader.h"
15 : #include "TimeIntegrator.h"
16 : #include "FEProblemBase.h"
17 :
18 : registerMooseObject("SolidMechanicsApp", NodalTranslationalInertia);
19 :
20 : InputParameters
21 642 : NodalTranslationalInertia::validParams()
22 : {
23 642 : InputParameters params = TimeNodalKernel::validParams();
24 642 : params.addClassDescription("Computes the inertial forces and mass proportional damping terms "
25 : "corresponding to nodal mass.");
26 1284 : params.addCoupledVar("velocity", "velocity variable");
27 1284 : params.addCoupledVar("acceleration", "acceleration variable");
28 1284 : params.addRangeCheckedParam<Real>(
29 : "beta", "beta>0.0", "beta parameter for Newmark Time integration");
30 1284 : params.addRangeCheckedParam<Real>(
31 : "gamma", "gamma>0.0", "gamma parameter for Newmark Time integration");
32 1926 : params.addRangeCheckedParam<Real>("eta",
33 1284 : 0.0,
34 : "eta>=0.0",
35 : "Constant real number defining the eta parameter for "
36 : "Rayleigh damping.");
37 1926 : params.addRangeCheckedParam<Real>("alpha",
38 1284 : 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 1284 : params.addParam<Real>("mass", "Mass associated with the node");
43 1284 : params.addParam<FileName>(
44 : "nodal_mass_file",
45 : "The file containing the nodal positions and the corresponding nodal masses.");
46 642 : return params;
47 0 : }
48 :
49 313 : NodalTranslationalInertia::NodalTranslationalInertia(const InputParameters & parameters)
50 : : TimeNodalKernel(parameters),
51 313 : _has_mass(isParamValid("mass")),
52 626 : _has_beta(isParamValid("beta")),
53 626 : _has_gamma(isParamValid("gamma")),
54 626 : _has_velocity(isParamValid("velocity")),
55 626 : _has_acceleration(isParamValid("acceleration")),
56 626 : _has_nodal_mass_file(isParamValid("nodal_mass_file")),
57 515 : _mass(_has_mass ? getParam<Real>("mass") : 0.0),
58 521 : _beta(_has_beta ? getParam<Real>("beta") : 0.1),
59 521 : _gamma(_has_gamma ? getParam<Real>("gamma") : 0.1),
60 626 : _eta(getParam<Real>("eta")),
61 626 : _alpha(getParam<Real>("alpha")),
62 626 : _time_integrator(_sys.getTimeIntegrator(_var.number()))
63 : {
64 313 : if (_has_beta && _has_gamma && _has_velocity && _has_acceleration)
65 : {
66 208 : _u_old = &(_var.dofValuesOld());
67 208 : _aux_sys = &(_fe_problem.getAuxiliarySystem());
68 :
69 416 : MooseVariable * vel_variable = getVar("velocity", 0);
70 208 : _vel_num = vel_variable->number();
71 416 : MooseVariable * accel_variable = getVar("acceleration", 0);
72 208 : _accel_num = accel_variable->number();
73 208 : }
74 105 : else if (!_has_beta && !_has_gamma && !_has_velocity && !_has_acceleration)
75 : {
76 105 : _du_dot_du = &(_var.duDotDu());
77 105 : _du_dotdot_du = &(_var.duDotDotDu());
78 105 : _u_dot_old = &(_var.dofValuesDotOld());
79 :
80 105 : addFEVariableCoupleableVectorTag(_time_integrator.uDotFactorTag());
81 105 : addFEVariableCoupleableVectorTag(_time_integrator.uDotDotFactorTag());
82 :
83 105 : _u_dot_factor = &_var.vectorTagDofValue(_time_integrator.uDotFactorTag());
84 105 : _u_dotdot_factor = &_var.vectorTagDofValue(_time_integrator.uDotDotFactorTag());
85 : }
86 : else
87 0 : mooseError("NodalTranslationalInertia: Either all or none of `beta`, `gamma`, `velocity` and "
88 : "`acceleration` should be provided as input.");
89 :
90 313 : if (!_has_nodal_mass_file && !_has_mass)
91 2 : mooseError(
92 : "NodalTranslationalInertia: Please provide either mass or nodal_mass_file as input.");
93 311 : else if (_has_nodal_mass_file && _has_mass)
94 4 : mooseError("NodalTranslationalInertia: Please provide either mass or nodal_mass_file as input, "
95 : "not both.");
96 :
97 307 : if (_has_nodal_mass_file)
98 : {
99 327 : MooseUtils::DelimitedFileReader nodal_mass_file(getParam<FileName>("nodal_mass_file"));
100 : nodal_mass_file.setHeaderFlag(MooseUtils::DelimitedFileReader::HeaderFlag::OFF);
101 109 : nodal_mass_file.read();
102 109 : std::vector<std::vector<Real>> data = nodal_mass_file.getData();
103 109 : if (data.size() != 4)
104 4 : mooseError("NodalTranslationalInertia: The number of columns in ",
105 : getParam<FileName>("nodal_mass_file"),
106 : " should be 4.");
107 107 : unsigned int node_found = 0;
108 107 : const std::set<BoundaryID> bnd_ids = BoundaryRestrictable::boundaryIDs();
109 235 : for (auto & bnd_id : bnd_ids)
110 : {
111 128 : const std::vector<dof_id_type> & bnd_node_set = _mesh.getNodeList(bnd_id);
112 928 : for (auto & bnd_node : bnd_node_set)
113 : {
114 800 : const Node & node = _mesh.nodeRef(bnd_node);
115 800 : _node_id_to_mass[bnd_node] = 0.0;
116 :
117 4225 : for (unsigned int i = 0; i < data[0].size(); ++i)
118 : {
119 2228 : if (MooseUtils::absoluteFuzzyEqual(data[0][i], node(0), 1e-6) &&
120 5529 : MooseUtils::absoluteFuzzyEqual(data[1][i], node(1), 1e-6) &&
121 : MooseUtils::absoluteFuzzyEqual(data[2][i], node(2), 1e-6))
122 : {
123 800 : _node_id_to_mass[bnd_node] = data[3][i];
124 800 : node_found += 1;
125 800 : break;
126 : }
127 : }
128 : }
129 : }
130 107 : if (node_found != data[0].size())
131 2 : mooseError("NodalTranslationalInertia: Out of ",
132 2 : data[0].size(),
133 : " nodal positions in ",
134 : getParam<FileName>("nodal_mass_file"),
135 : " only ",
136 : node_found,
137 : " nodes were found in the boundary.");
138 105 : }
139 :
140 : // Check for Explicit and alpha parameter
141 303 : if (_alpha != 0 && _time_integrator.isExplicit())
142 0 : mooseError("NodalTranslationalInertia: HHT time integration parameter can only be used with "
143 : "Newmark-Beta time integration.");
144 :
145 : // Check if beta and explicit are being used simultaneously
146 303 : if (_has_beta && _time_integrator.isExplicit())
147 0 : mooseError("NodalTranslationalInertia: Newmark-beta integration parameter, beta, cannot be "
148 : "provided along with an explicit time "
149 : "integrator.");
150 303 : }
151 :
152 : Real
153 41856 : NodalTranslationalInertia::computeQpResidual()
154 : {
155 41856 : if (_dt == 0)
156 : return 0;
157 : else
158 : {
159 : Real mass = 0.0;
160 41856 : if (_has_mass)
161 19554 : mass = _mass;
162 : else
163 : {
164 22302 : if (_node_id_to_mass.find(_current_node->id()) != _node_id_to_mass.end())
165 22302 : mass = _node_id_to_mass[_current_node->id()];
166 : else
167 0 : mooseError("NodalTranslationalInertia: Unable to find an entry for the current node in the "
168 : "_node_id_to_mass map.");
169 : }
170 :
171 41856 : if (_has_beta)
172 : {
173 : mooseAssert(_beta > 0.0, "NodalTranslationalInertia: Beta parameter should be positive.");
174 21054 : const NumericVector<Number> & aux_sol_old = _aux_sys->solutionOld();
175 :
176 21054 : const Real vel_old = aux_sol_old(_current_node->dof_number(_aux_sys->number(), _vel_num, 0));
177 : const Real accel_old =
178 21054 : aux_sol_old(_current_node->dof_number(_aux_sys->number(), _accel_num, 0));
179 :
180 : const Real accel =
181 21054 : 1. / _beta *
182 21054 : (((_u[_qp] - (*_u_old)[_qp]) / (_dt * _dt)) - vel_old / _dt - accel_old * (0.5 - _beta));
183 21054 : const Real vel = vel_old + (_dt * (1 - _gamma)) * accel_old + _gamma * _dt * accel;
184 21054 : return mass * (accel + vel * _eta * (1 + _alpha) - _alpha * _eta * vel_old);
185 : }
186 :
187 : else
188 : // all cases (Explicit, implicit and implicit with HHT)
189 : // Note that _alpha is enforced to be zero for explicit integration
190 20802 : return mass * ((*_u_dotdot_factor)[_qp] + (*_u_dot_factor)[_qp] * _eta * (1.0 + _alpha) -
191 20802 : _alpha * _eta * (*_u_dot_old)[_qp]);
192 : }
193 : }
194 :
195 : Real
196 23118 : NodalTranslationalInertia::computeQpJacobian()
197 : {
198 : mooseAssert(_beta > 0.0, "NodalTranslationalInertia: Beta parameter should be positive.");
199 :
200 23118 : if (_dt == 0)
201 : return 0;
202 : else
203 : {
204 : Real mass = 0.0;
205 23118 : if (_has_mass)
206 9780 : mass = _mass;
207 : else
208 : {
209 13338 : if (_node_id_to_mass.find(_current_node->id()) != _node_id_to_mass.end())
210 13338 : mass = _node_id_to_mass[_current_node->id()];
211 : else
212 0 : mooseError("NodalTranslationalInertia: Unable to find an entry for the current node in the "
213 : "_node_id_to_mass map.");
214 : }
215 23118 : if (_has_beta)
216 10527 : return mass / (_beta * _dt * _dt) + _eta * (1 + _alpha) * mass * _gamma / _beta / _dt;
217 : else
218 12591 : return mass * (*_du_dotdot_du)[_qp] + _eta * (1.0 + _alpha) * mass * (*_du_dot_du)[_qp];
219 : }
220 : }
|
