Line data Source code
1 : // MASTODON includes
2 : #include "ComputeFVDamperElasticity.h"
3 :
4 : // MOOSE includes
5 : #include "MooseMesh.h"
6 : #include "Assembly.h"
7 : #include "NonlinearSystem.h"
8 : #include "MooseVariable.h"
9 : #include "ColumnMajorMatrix.h"
10 : #include "MathUtils.h"
11 :
12 : // libmesh includes
13 : #include "libmesh/quadrature.h"
14 : #include "libmesh/utility.h"
15 :
16 : registerMooseObject("MastodonApp", ComputeFVDamperElasticity);
17 :
18 : InputParameters
19 60 : ComputeFVDamperElasticity::validParams()
20 : {
21 60 : InputParameters params = Material::validParams();
22 60 : params.addClassDescription("Compute the deformations, forces and the stiffness "
23 : "matrix corresponding to a two-node nonlinear fluid "
24 : "viscous damper element.");
25 120 : params.addRequiredCoupledVar("displacements",
26 : "The displacement variables appropriate for the simulation of "
27 : "geometry and coordinate system.");
28 120 : params.addRequiredParam<RealGradient>("y_orientation",
29 : "Orientation of the y direction along "
30 : "with Ky is provided. This should be "
31 : "perpendicular to the axis of the Damper.");
32 :
33 : // Input parameters
34 120 : params.addRequiredRangeCheckedParam<Real>("cd", "cd > 0.0", "Damping co-efficient.");
35 120 : params.addRequiredRangeCheckedParam<Real>(
36 : "alpha", "alpha > 0.0", "Velocity exponent of the damper.");
37 120 : params.addRequiredRangeCheckedParam<Real>(
38 : "k", "k > 0.0", "Axial stiffness of the damper assembly.");
39 120 : params.addRequiredRangeCheckedParam<Real>(
40 : "gamma", "gamma > 0.0", "Gamma parameter of Newmark algorithm.");
41 120 : params.addRequiredRangeCheckedParam<Real>(
42 : "beta", "beta > 0.0", "Beta parameter of Newmark algorithm.");
43 120 : params.addParam<Real>("rel_tol", 1e-6, "Relative tolerance for error in adaptive algorithm");
44 120 : params.addParam<Real>("abs_tol", 1e-6, "Absolute tolerance for error in adaptive algorithm");
45 120 : params.set<MooseEnum>("constant_on") = "ELEMENT";
46 60 : return params;
47 0 : }
48 :
49 45 : ComputeFVDamperElasticity::ComputeFVDamperElasticity(const InputParameters & parameters)
50 : : Material(parameters),
51 45 : _ndisp(coupledComponents("displacements")),
52 45 : _disp_num(3),
53 45 : _length(declareProperty<Real>("initial_Damper_length")),
54 90 : _cd(getParam<Real>("cd")),
55 90 : _alpha(getParam<Real>("alpha")),
56 90 : _k(getParam<Real>("k")),
57 90 : _gamma(getParam<Real>("gamma")),
58 90 : _beta(getParam<Real>("beta")),
59 90 : _rel_tol(getParam<Real>("rel_tol")),
60 90 : _abs_tol(getParam<Real>("abs_tol")),
61 45 : _basic_def(declareProperty<Real>("basic_deformation")),
62 90 : _basic_def_old(getMaterialPropertyOld<Real>("basic_deformation")),
63 45 : _vel(declareProperty<Real>("velocity")),
64 90 : _vel_old(getMaterialPropertyOld<Real>("velocity")),
65 45 : _Fb(declareProperty<Real>("basic_force")),
66 90 : _Fb_old(getMaterialPropertyOld<Real>("basic_force")),
67 45 : _Fl(declareProperty<ColumnMajorMatrix>("local_forces")),
68 45 : _Fg(declareProperty<ColumnMajorMatrix>("global_forces")),
69 45 : _Kb(declareProperty<Real>("basic_stiffness")),
70 45 : _Kl(declareProperty<ColumnMajorMatrix>("local_stiffness_matrix")),
71 45 : _Kg(declareProperty<ColumnMajorMatrix>("global_stiffness_matrix")),
72 90 : _total_gl(declareProperty<ColumnMajorMatrix>("total_global_to_local_transformation"))
73 :
74 : {
75 : // Fetch coupled variables (as stateful properties if necessary)
76 180 : for (unsigned int i = 0; i < _ndisp; ++i)
77 : {
78 270 : MooseVariable * disp_variable = getVar("displacements", i);
79 135 : _disp_num[i] = disp_variable->number(); // Displacement variable numbers in MOOSE
80 : }
81 45 : }
82 :
83 : void
84 40 : ComputeFVDamperElasticity::initQpStatefulProperties()
85 : {
86 40 : _vel[_qp] = 0;
87 40 : }
88 :
89 : void
90 18552 : ComputeFVDamperElasticity::computeQpProperties()
91 : {
92 : // Compute transformation matrices
93 18552 : computeTransformationMatrix();
94 :
95 : // Compute the global displacement values
96 18552 : computeDeformation();
97 :
98 : // Initialize the stiffness matrices and force vectors
99 18552 : initializeFVdamper();
100 :
101 : // Compute axial forces and stiffness terms
102 18552 : computeAxial();
103 :
104 : // Finalize forces and stiffness matrix and convert them into global co-ordinate system
105 18552 : finalize();
106 18552 : }
107 :
108 : void
109 18552 : ComputeFVDamperElasticity::computeTransformationMatrix()
110 : {
111 : // Fetch the two nodes of the link element
112 : std::vector<const Node *> node;
113 55656 : for (unsigned int i = 0; i < 2; ++i)
114 37104 : node.push_back(_current_elem->node_ptr(i));
115 :
116 : // Defining orientation of damper (direction cosines)
117 : RealGradient x_orientation;
118 74208 : for (unsigned int i = 0; i < _ndisp; ++i)
119 55656 : x_orientation(i) = (*node[1])(i) - (*node[0])(i);
120 :
121 : // Length of the damper element
122 18552 : _length[_qp] = x_orientation.norm();
123 :
124 18552 : if (_length[_qp] == 0.0)
125 0 : mooseError("Error in ComputeFVDamperElasticity material block, ",
126 : name(),
127 : ". Damper element cannot be of zero length.");
128 :
129 : // Normalizing the orientation vector
130 : x_orientation /= _length[_qp];
131 :
132 : // Get y orientation of the Damper in global coordinate system
133 37104 : RealGradient y_orientation = getParam<RealGradient>("y_orientation");
134 : Real dot = x_orientation * y_orientation;
135 :
136 : // Check if x and y orientations are perpendicular
137 18552 : if (abs(dot) > 1e-4)
138 0 : mooseError("Error in ComputeFVDamperElasticity material block, ",
139 : name(),
140 : ". y_orientation should be perpendicular to the axis of the Damper.");
141 :
142 : // Transformation matrix from global to local coordinate system
143 18552 : _total_gl[_qp].reshape(2, 6);
144 18552 : _total_gl[_qp].zero();
145 18552 : _total_gl[_qp](0, 0) = _total_gl[_qp](1, 3) = x_orientation(0); // direction cosine in x
146 18552 : _total_gl[_qp](0, 1) = _total_gl[_qp](1, 4) = x_orientation(1); // direction cosine in y
147 18552 : _total_gl[_qp](0, 2) = _total_gl[_qp](1, 5) = x_orientation(2); // direction cosine in z
148 18552 : }
149 :
150 : void
151 18552 : ComputeFVDamperElasticity::computeDeformation()
152 : {
153 : // Fetch the two end nodes for _current_elem
154 : std::vector<const Node *> node;
155 55656 : for (unsigned int i = 0; i < 2; ++i)
156 37104 : node.push_back(_current_elem->node_ptr(i));
157 :
158 : // Fetch the solution for the two end nodes at current time
159 18552 : NonlinearSystemBase & nonlinear_sys = _fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0);
160 18552 : const NumericVector<Number> & sol = *nonlinear_sys.currentSolution();
161 18552 : const NumericVector<Number> & sol_old = nonlinear_sys.solutionOld();
162 :
163 : // Calculating global displacements 6 x 1 matrix with
164 : // first three rows corresponding to node 0 dofs and next three to node 1 dofs
165 18552 : ColumnMajorMatrix global_disp(6, 1);
166 18552 : ColumnMajorMatrix global_disp_old(6, 1);
167 :
168 74208 : for (unsigned int i = 0; i < _ndisp; ++i)
169 : {
170 55656 : global_disp(i) =
171 55656 : sol(node[0]->dof_number(nonlinear_sys.number(), _disp_num[i], 0)); // node 0 displacements
172 55656 : global_disp(i + 3) =
173 55656 : sol(node[1]->dof_number(nonlinear_sys.number(), _disp_num[i], 0)); // node 1 displacements
174 55656 : global_disp_old(i) = sol_old(
175 : node[0]->dof_number(nonlinear_sys.number(), _disp_num[i], 0)); // node 0 displacements
176 55656 : global_disp_old(i + 3) = sol_old(
177 : node[1]->dof_number(nonlinear_sys.number(), _disp_num[i], 0)); // node 1 displacements
178 : }
179 :
180 18552 : _basic_def[_qp] = (_total_gl[_qp] * global_disp)(1, 0) - (_total_gl[_qp] * global_disp)(0, 0);
181 18552 : }
182 :
183 : void
184 18552 : ComputeFVDamperElasticity::initializeFVdamper()
185 : {
186 : // initialize basic force and stiffness terms
187 18552 : _Kb[_qp] = 0;
188 18552 : _Fb[_qp] = 0;
189 :
190 : // Initialize local force and stiffness terms
191 18552 : _Kl[_qp].reshape(2, 2);
192 18552 : _Fl[_qp].reshape(2, 1);
193 18552 : }
194 :
195 : void
196 18552 : ComputeFVDamperElasticity::computeAxial()
197 : {
198 18552 : _vel[_qp] = _vel_old[_qp] + (_gamma / _beta) * (((_basic_def[_qp] - _basic_def_old[_qp]) / _dt) -
199 : (_vel_old[_qp]));
200 :
201 : // adaptive iterative algorithm starts
202 :
203 : Real t_vel = _vel[_qp];
204 18552 : Real accel = (_vel[_qp] - _vel_old[_qp]) / _dt;
205 18552 : Real t_fd = _Fb_old[_qp];
206 : Real s = 1;
207 : Real nit = 0;
208 : Real s_tot = 0;
209 : Real yt, eps, error;
210 :
211 : do
212 : {
213 225772 : Real h = s * _dt;
214 225772 : Real t_vel_new = t_vel + accel * h;
215 :
216 225772 : dormandPrince(t_vel, t_vel_new, t_fd, h, yt, eps, error);
217 :
218 225772 : if (eps <= _rel_tol || fabs(eps) <= _abs_tol)
219 : {
220 : t_vel = t_vel_new;
221 196612 : t_fd = yt;
222 196612 : s_tot = s_tot + s;
223 : }
224 : else
225 : {
226 29160 : nit = nit + 1;
227 29160 : s = pow(2, -nit);
228 : }
229 :
230 225772 : } while (s_tot < 1.0);
231 :
232 : // set axial force
233 18552 : _Fb[_qp] = t_fd;
234 :
235 : // set axial stiffness
236 18552 : if (_vel[_qp] != 0)
237 : {
238 18430 : Real k_damp = _cd * MathUtils::sign(_vel[_qp]) * _alpha * pow(fabs(_vel[_qp]), _alpha - 1.0) *
239 18430 : (_gamma / _beta) / _dt;
240 18430 : Real k_spring = _k;
241 :
242 : // springs in series connection
243 18430 : _Kb[_qp] = k_damp * k_spring / (k_damp + k_spring);
244 : }
245 18552 : }
246 :
247 : void
248 225772 : ComputeFVDamperElasticity::dormandPrince(
249 : Real t_vel, Real t_vel_new, Real t_fd, Real h, Real & yt, Real & eps, Real & error)
250 : {
251 225772 : Real k1 = fdot(t_vel, t_fd) * h;
252 225772 : Real k2 = fdot((t_vel_new - t_vel) * (1.0 / 5.0) + t_vel, t_fd + (1.0 / 5.0) * k1) * h;
253 225772 : Real k3 = fdot((t_vel_new - t_vel) * (3.0 / 10.0) + t_vel,
254 225772 : t_fd + (3.0 / 40.0) * k1 + (9.0 / 40.0) * k2) *
255 225772 : h;
256 225772 : Real k4 = fdot((t_vel_new - t_vel) * (4.0 / 5.0) + t_vel,
257 225772 : t_fd + (44.0 / 45.0) * k1 + (-56.0 / 15.0) * k2 + (32.0 / 9.0) * k3) *
258 225772 : h;
259 225772 : Real k5 = fdot((t_vel_new - t_vel) * (8.0 / 9.0) + t_vel,
260 225772 : t_fd + (19372.0 / 6561.0) * k1 + (-25360.0 / 2187.0) * k2 +
261 225772 : (64448.0 / 6561.0) * k3 + (-212.0 / 729.0) * k4) *
262 225772 : h;
263 225772 : Real k6 = fdot((t_vel_new - t_vel) * (1.0) + t_vel,
264 225772 : t_fd + (9017.0 / 3168.0) * k1 + (-355.0 / 33.0) * k2 + (46732.0 / 5247.0) * k3 +
265 225772 : (49.0 / 176.0) * k4 + (-5103.0 / 18656.0) * k5) *
266 225772 : h;
267 :
268 225772 : yt = t_fd + (35.0 / 384.0) * k1 + (500.0 / 1113.0) * k3 + (125.0 / 192.0) * k4 +
269 225772 : (-2187.0 / 6784.0) * k5 + (11.0 / 84.0) * k6;
270 :
271 225772 : Real k7 = fdot((t_vel_new - t_vel) * (1.0) + t_vel, yt) * h;
272 :
273 225772 : error = (71.0 / 57600.0) * k1 + (-71.0 / 16695.0) * k3 + (71.0 / 1920.0) * k4 +
274 225772 : (-17253.0 / 339200.) * k5 + (22.0 / 525.0) * k6 + (-1.0 / 40.0) * k7;
275 :
276 225772 : eps = fabs(error);
277 225772 : }
278 :
279 : Real
280 1580404 : ComputeFVDamperElasticity::fdot(Real v, Real fd)
281 : {
282 1580404 : return (v - (MathUtils::sign(fd) * pow(fabs(fd) / _cd, 1.0 / _alpha))) * _k;
283 : }
284 :
285 : void
286 18552 : ComputeFVDamperElasticity::finalize()
287 : {
288 : // calculate local forces
289 18552 : _Fl[_qp](0, 0) = -_Fb[_qp]; // force at node 1
290 18552 : _Fl[_qp](1, 0) = _Fb[_qp]; // force at node 2
291 :
292 : // calculate global forces from local system
293 18552 : _Fg[_qp] = _total_gl[_qp].transpose() * _Fl[_qp];
294 :
295 : // populate the local stiffness matrix
296 55656 : for (unsigned int i = 0; i < 2; ++i)
297 111312 : for (unsigned int j = 0; j < 2; j++)
298 111312 : _Kl[_qp](i, j) = (i == j ? 1 : -1) * _Kb[_qp];
299 :
300 : // convert local stiffness matrix to global coordinate system
301 18552 : _Kg[_qp] = _total_gl[_qp].transpose() * _Kl[_qp] * _total_gl[_qp];
302 18552 : }
|
