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 566 : NodalTranslationalInertia::validParams()
22 : {
23 566 : InputParameters params = TimeNodalKernel::validParams();
24 566 : params.addClassDescription("Computes the inertial forces and mass proportional damping terms "
25 : "corresponding to nodal mass.");
26 1132 : params.addCoupledVar("velocity", "velocity variable");
27 1132 : params.addCoupledVar("acceleration", "acceleration variable");
28 1132 : params.addRangeCheckedParam<Real>(
29 : "beta", "beta>0.0", "beta parameter for Newmark Time integration");
30 1132 : params.addRangeCheckedParam<Real>(
31 : "gamma", "gamma>0.0", "gamma parameter for Newmark Time integration");
32 1698 : params.addRangeCheckedParam<Real>("eta",
33 1132 : 0.0,
34 : "eta>=0.0",
35 : "Constant real number defining the eta parameter for "
36 : "Rayleigh damping.");
37 1698 : params.addRangeCheckedParam<Real>("alpha",
38 1132 : 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 1132 : params.addParam<Real>("mass", "Mass associated with the node");
43 1132 : params.addParam<FileName>(
44 : "nodal_mass_file",
45 : "The file containing the nodal positions and the corresponding nodal masses.");
46 566 : return params;
47 0 : }
48 :
49 274 : NodalTranslationalInertia::NodalTranslationalInertia(const InputParameters & parameters)
50 : : TimeNodalKernel(parameters),
51 274 : _has_mass(isParamValid("mass")),
52 548 : _has_beta(isParamValid("beta")),
53 548 : _has_gamma(isParamValid("gamma")),
54 548 : _has_velocity(isParamValid("velocity")),
55 548 : _has_acceleration(isParamValid("acceleration")),
56 548 : _has_nodal_mass_file(isParamValid("nodal_mass_file")),
57 452 : _mass(_has_mass ? getParam<Real>("mass") : 0.0),
58 458 : _beta(_has_beta ? getParam<Real>("beta") : 0.1),
59 458 : _gamma(_has_gamma ? getParam<Real>("gamma") : 0.1),
60 548 : _eta(getParam<Real>("eta")),
61 548 : _alpha(getParam<Real>("alpha")),
62 548 : _time_integrator(_sys.getTimeIntegrator(_var.number()))
63 : {
64 274 : if (_has_beta && _has_gamma && _has_velocity && _has_acceleration)
65 : {
66 184 : _u_old = &(_var.dofValuesOld());
67 184 : _aux_sys = &(_fe_problem.getAuxiliarySystem());
68 :
69 368 : MooseVariable * vel_variable = getVar("velocity", 0);
70 184 : _vel_num = vel_variable->number();
71 368 : MooseVariable * accel_variable = getVar("acceleration", 0);
72 184 : _accel_num = accel_variable->number();
73 184 : }
74 90 : else if (!_has_beta && !_has_gamma && !_has_velocity && !_has_acceleration)
75 : {
76 90 : _du_dot_du = &(_var.duDotDu());
77 90 : _du_dotdot_du = &(_var.duDotDotDu());
78 90 : _u_dot_old = &(_var.dofValuesDotOld());
79 :
80 90 : addFEVariableCoupleableVectorTag(_time_integrator.uDotFactorTag());
81 90 : addFEVariableCoupleableVectorTag(_time_integrator.uDotDotFactorTag());
82 :
83 90 : _u_dot_factor = &_var.vectorTagDofValue(_time_integrator.uDotFactorTag());
84 90 : _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 274 : if (!_has_nodal_mass_file && !_has_mass)
91 2 : mooseError(
92 : "NodalTranslationalInertia: Please provide either mass or nodal_mass_file as input.");
93 272 : 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 268 : if (_has_nodal_mass_file)
98 : {
99 282 : MooseUtils::DelimitedFileReader nodal_mass_file(getParam<FileName>("nodal_mass_file"));
100 : nodal_mass_file.setHeaderFlag(MooseUtils::DelimitedFileReader::HeaderFlag::OFF);
101 94 : nodal_mass_file.read();
102 94 : std::vector<std::vector<Real>> data = nodal_mass_file.getData();
103 94 : if (data.size() != 4)
104 4 : mooseError("NodalTranslationalInertia: The number of columns in ",
105 : getParam<FileName>("nodal_mass_file"),
106 : " should be 4.");
107 92 : unsigned int node_found = 0;
108 92 : const std::set<BoundaryID> bnd_ids = BoundaryRestrictable::boundaryIDs();
109 202 : for (auto & bnd_id : bnd_ids)
110 : {
111 110 : const std::vector<dof_id_type> & bnd_node_set = _mesh.getNodeList(bnd_id);
112 796 : for (auto & bnd_node : bnd_node_set)
113 : {
114 686 : const Node & node = _mesh.nodeRef(bnd_node);
115 686 : _node_id_to_mass[bnd_node] = 0.0;
116 :
117 3622 : for (unsigned int i = 0; i < data[0].size(); ++i)
118 : {
119 1910 : if (MooseUtils::absoluteFuzzyEqual(data[0][i], node(0), 1e-6) &&
120 4740 : MooseUtils::absoluteFuzzyEqual(data[1][i], node(1), 1e-6) &&
121 : MooseUtils::absoluteFuzzyEqual(data[2][i], node(2), 1e-6))
122 : {
123 686 : _node_id_to_mass[bnd_node] = data[3][i];
124 686 : node_found += 1;
125 686 : break;
126 : }
127 : }
128 : }
129 : }
130 92 : 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 90 : }
139 :
140 : // Check for Explicit and alpha parameter
141 264 : 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 264 : 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 264 : }
151 :
152 : Real
153 33366 : NodalTranslationalInertia::computeQpResidual()
154 : {
155 33366 : if (_dt == 0)
156 : return 0;
157 : else
158 : {
159 : Real mass = 0.0;
160 33366 : if (_has_mass)
161 15606 : mass = _mass;
162 : else
163 : {
164 17760 : if (_node_id_to_mass.find(_current_node->id()) != _node_id_to_mass.end())
165 17760 : 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 33366 : if (_has_beta)
172 : {
173 : mooseAssert(_beta > 0.0, "NodalTranslationalInertia: Beta parameter should be positive.");
174 16800 : const NumericVector<Number> & aux_sol_old = _aux_sys->solutionOld();
175 :
176 16800 : const Real vel_old = aux_sol_old(_current_node->dof_number(_aux_sys->number(), _vel_num, 0));
177 : const Real accel_old =
178 16800 : aux_sol_old(_current_node->dof_number(_aux_sys->number(), _accel_num, 0));
179 :
180 : const Real accel =
181 16800 : 1. / _beta *
182 16800 : (((_u[_qp] - (*_u_old)[_qp]) / (_dt * _dt)) - vel_old / _dt - accel_old * (0.5 - _beta));
183 16800 : const Real vel = vel_old + (_dt * (1 - _gamma)) * accel_old + _gamma * _dt * accel;
184 16800 : 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 16566 : return mass * ((*_u_dotdot_factor)[_qp] + (*_u_dot_factor)[_qp] * _eta * (1.0 + _alpha) -
191 16566 : _alpha * _eta * (*_u_dot_old)[_qp]);
192 : }
193 : }
194 :
195 : Real
196 18402 : NodalTranslationalInertia::computeQpJacobian()
197 : {
198 : mooseAssert(_beta > 0.0, "NodalTranslationalInertia: Beta parameter should be positive.");
199 :
200 18402 : if (_dt == 0)
201 : return 0;
202 : else
203 : {
204 : Real mass = 0.0;
205 18402 : if (_has_mass)
206 7806 : mass = _mass;
207 : else
208 : {
209 10596 : if (_node_id_to_mass.find(_current_node->id()) != _node_id_to_mass.end())
210 10596 : 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 18402 : if (_has_beta)
216 8400 : return mass / (_beta * _dt * _dt) + _eta * (1 + _alpha) * mass * _gamma / _beta / _dt;
217 : else
218 10002 : return mass * (*_du_dotdot_du)[_qp] + _eta * (1.0 + _alpha) * mass * (*_du_dot_du)[_qp];
219 : }
220 : }
|
