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