Line data Source code
1 : /*************************************************/
2 : /* DO NOT MODIFY THIS HEADER */
3 : /* */
4 : /* MASTODON */
5 : /* */
6 : /* (c) 2015 Battelle Energy Alliance, LLC */
7 : /* ALL RIGHTS RESERVED */
8 : /* */
9 : /* Prepared by Battelle Energy Alliance, LLC */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* See COPYRIGHT for full restrictions */
13 : /*************************************************/
14 :
15 : // MASTODON includes
16 : #include "ComputeIsolatorDeformation.h"
17 :
18 : // MOOSE includes
19 : #include "MooseMesh.h"
20 : #include "Assembly.h"
21 : #include "NonlinearSystem.h"
22 : #include "MooseVariable.h"
23 : #include "ColumnMajorMatrix.h"
24 :
25 : // libmesh includes
26 : #include "libmesh/quadrature.h"
27 : #include "libmesh/utility.h"
28 :
29 : registerMooseObject("MastodonApp", ComputeIsolatorDeformation);
30 :
31 : InputParameters
32 165 : ComputeIsolatorDeformation::validParams()
33 : {
34 165 : InputParameters params = Material::validParams();
35 165 : params.addClassDescription(
36 : "Compute the deformations and rotations in a two-noded isolator element.");
37 330 : params.addRequiredCoupledVar(
38 : "rotations",
39 : "The rotation variables appropriate for the simulation geometry and coordinate system.");
40 330 : params.addRequiredCoupledVar(
41 : "displacements",
42 : "The displacement variables appropriate for the simulation geometry and coordinate system.");
43 330 : params.addRequiredCoupledVar("velocities", "Translational velocity variables.");
44 330 : params.addRequiredParam<RealGradient>("y_orientation",
45 : "Orientation of the local Y direction along "
46 : "which, Ky is provided. This should be "
47 : "perpendicular to the axis of the isolator.");
48 330 : params.addParam<MaterialPropertyName>("sd_ratio", 0.5, "Shear distance ratio.");
49 330 : params.set<MooseEnum>("constant_on") = "ELEMENT"; // sets _qp to 0. Material properties are
50 : // assumed to be constant throughout the length
51 : // of the element.
52 165 : return params;
53 0 : }
54 :
55 124 : ComputeIsolatorDeformation::ComputeIsolatorDeformation(const InputParameters & parameters)
56 : : Material(parameters),
57 124 : _nrot(coupledComponents("rotations")),
58 124 : _ndisp(coupledComponents("displacements")),
59 124 : _nvel(coupledComponents("velocities")),
60 124 : _rot_num(3),
61 124 : _disp_num(3),
62 124 : _vel_num(3),
63 248 : _sD(getMaterialProperty<Real>("sd_ratio")),
64 124 : _local_disp(declareProperty<ColumnMajorMatrix>("local_deformations")),
65 124 : _basic_disp(declareProperty<ColumnMajorMatrix>("deformations")),
66 124 : _basic_disp_old(declareProperty<ColumnMajorMatrix>("old_deformations")),
67 124 : _basic_vel(declareProperty<ColumnMajorMatrix>("deformation_rates")),
68 124 : _basic_vel_old(declareProperty<ColumnMajorMatrix>("old_deformation_rates")),
69 124 : _original_gl(declareProperty<ColumnMajorMatrix>("original_global_to_local_transformation")),
70 124 : _total_gl(declareProperty<ColumnMajorMatrix>("total_global_to_local_transformation")),
71 124 : _total_lb(declareProperty<ColumnMajorMatrix>("total_local_to_basic_transformation")),
72 248 : _length(declareProperty<Real>("initial_isolator_length"))
73 : {
74 :
75 : // Checking for consistency between length of the provided displacements and rotations vector
76 124 : if (_ndisp != _nrot)
77 1 : mooseError("Error in ComputeIsolatorDeformation block, ",
78 : name(),
79 : ". The number of variables supplied in 'displacements' "
80 : "and 'rotations' input parameters must be equal.");
81 :
82 : // Fetch coupled variables and gradients (as stateful properties if necessary)
83 492 : for (unsigned int i = 0; i < _ndisp; ++i)
84 : {
85 738 : MooseVariable * disp_variable = getVar("displacements", i);
86 369 : _disp_num[i] = disp_variable->number(); // Displacement variable numbers in MOOSE
87 :
88 738 : MooseVariable * vel_variable = getVar("velocities", i);
89 369 : _vel_num[i] = vel_variable->number();
90 :
91 738 : MooseVariable * rot_variable = getVar("rotations", i);
92 369 : _rot_num[i] = rot_variable->number(); // Rotation variable numbers in MOOSE
93 : }
94 123 : }
95 :
96 : void
97 44425 : ComputeIsolatorDeformation::computeQpProperties()
98 : {
99 : // Compute initial orientation and length of the isolator in global coordinate system
100 : // Fetch the two nodes of the link element
101 : std::vector<const Node *> node;
102 133275 : for (unsigned int i = 0; i < 2; ++i)
103 88850 : node.push_back(_current_elem->node_ptr(i));
104 : RealGradient x_orientation;
105 177700 : for (unsigned int i = 0; i < _ndisp; ++i)
106 133275 : x_orientation(i) = (*node[1])(i) - (*node[0])(i);
107 44425 : _length[_qp] = x_orientation.norm();
108 44425 : if (_length[_qp] == 0.0)
109 0 : mooseError("Error in ComputeIsolatorDeformation block, ",
110 : name(),
111 : ". Isolator element cannot be of zero length.");
112 : x_orientation /= _length[_qp]; // Normalizing with length to get orientation
113 :
114 : // Get y orientation of the isolator in global coordinate system
115 88850 : RealGradient y_orientation = getParam<RealGradient>("y_orientation");
116 : Real dot = x_orientation * y_orientation;
117 :
118 : // Check if x and y orientations are perpendicular
119 44425 : if (abs(dot) > 1e-4)
120 1 : mooseError("Error in ComputeIsolatorDeformation block, ",
121 : name(),
122 : ". 'y_orientation' should be perpendicular to ",
123 : "the axis of the isolator.");
124 :
125 : // Calculate z orientation in the global coordinate system as a cross product of the x and y
126 : // orientations
127 44424 : RealGradient z_orientation = x_orientation.cross(y_orientation);
128 :
129 : // Create rotation matrix from global to local coordinate system
130 44424 : _original_gl[_qp].reshape(12, 12);
131 44424 : _original_gl[_qp].zero();
132 44424 : _original_gl[_qp](0, 0) = _original_gl[_qp](3, 3) = _original_gl[_qp](6, 6) =
133 44424 : _original_gl[_qp](9, 9) = x_orientation(0);
134 44424 : _original_gl[_qp](0, 1) = _original_gl[_qp](3, 4) = _original_gl[_qp](6, 7) =
135 44424 : _original_gl[_qp](9, 10) = x_orientation(1);
136 44424 : _original_gl[_qp](0, 2) = _original_gl[_qp](3, 5) = _original_gl[_qp](6, 8) =
137 44424 : _original_gl[_qp](9, 11) = x_orientation(2);
138 44424 : _original_gl[_qp](1, 0) = _original_gl[_qp](4, 3) = _original_gl[_qp](7, 6) =
139 44424 : _original_gl[_qp](10, 9) = y_orientation(0);
140 44424 : _original_gl[_qp](1, 1) = _original_gl[_qp](4, 4) = _original_gl[_qp](7, 7) =
141 44424 : _original_gl[_qp](10, 10) = y_orientation(1);
142 44424 : _original_gl[_qp](1, 2) = _original_gl[_qp](4, 5) = _original_gl[_qp](7, 8) =
143 44424 : _original_gl[_qp](10, 11) = y_orientation(2);
144 44424 : _original_gl[_qp](2, 0) = _original_gl[_qp](5, 3) = _original_gl[_qp](8, 6) =
145 44424 : _original_gl[_qp](11, 9) = z_orientation(0);
146 44424 : _original_gl[_qp](2, 1) = _original_gl[_qp](5, 4) = _original_gl[_qp](8, 7) =
147 44424 : _original_gl[_qp](11, 10) = z_orientation(1);
148 44424 : _original_gl[_qp](2, 2) = _original_gl[_qp](5, 5) = _original_gl[_qp](8, 8) =
149 44424 : _original_gl[_qp](11, 11) = z_orientation(2);
150 :
151 : // Create rotation matrix from local to basic system (linear)
152 44424 : _total_lb[_qp].reshape(6, 12);
153 44424 : _total_lb[_qp].zero();
154 44424 : _total_lb[_qp](0, 0) = _total_lb[_qp](1, 1) = _total_lb[_qp](2, 2) = _total_lb[_qp](3, 3) =
155 44424 : _total_lb[_qp](4, 4) = _total_lb[_qp](5, 5) = -1.0;
156 44424 : _total_lb[_qp](0, 6) = _total_lb[_qp](1, 7) = _total_lb[_qp](2, 8) = _total_lb[_qp](3, 9) =
157 44424 : _total_lb[_qp](4, 10) = _total_lb[_qp](5, 11) = 1.0;
158 44424 : _total_lb[_qp](1, 5) = -_sD[0] * _length[_qp];
159 44424 : _total_lb[_qp](1, 11) = -(1.0 - _sD[0]) * _length[_qp];
160 44424 : _total_lb[_qp](2, 4) = -_total_lb[_qp](1, 5);
161 44424 : _total_lb[_qp](2, 10) = -_total_lb[_qp](1, 11);
162 :
163 44424 : computeTotalRotation();
164 44424 : computeDeformation();
165 44424 : }
166 :
167 : void
168 44424 : ComputeIsolatorDeformation::computeTotalRotation()
169 : {
170 : std::size_t _qp = 0;
171 : // Currently, this formulation is limited to small deformations in the isolator,
172 : // namely, it is assumed that there are no rigid body rotations in the isolator,
173 : // and that the total rotation matrix from global to local coordinates
174 : // (calculated below) remains the same as the one at t = 0 throughout the
175 : // duration of the simulation.
176 44424 : _total_gl[_qp] = _original_gl[_qp];
177 44424 : }
178 :
179 : void
180 44424 : ComputeIsolatorDeformation::computeDeformation()
181 : {
182 : // fetch the two end nodes for _current_elem
183 : std::vector<const Node *> node;
184 133272 : for (unsigned int i = 0; i < 2; ++i)
185 88848 : node.push_back(_current_elem->node_ptr(i));
186 :
187 : // Fetch the solution for the two end nodes at time t
188 44424 : NonlinearSystemBase & nonlinear_sys = _fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0);
189 44424 : const NumericVector<Number> & sol = *nonlinear_sys.currentSolution();
190 44424 : const NumericVector<Number> & sol_old = nonlinear_sys.solutionOld();
191 :
192 44424 : AuxiliarySystem & aux = _fe_problem.getAuxiliarySystem();
193 44424 : const NumericVector<Number> & aux_sol = *aux.currentSolution();
194 44424 : const NumericVector<Number> & aux_sol_old = aux.solutionOld();
195 :
196 : // Calculating global displacements (including rotations) and velocities (rotational velocity
197 : // terms are zero) 12 x 1 matrix with first six rows corresponding to node 0 dofs and next six to
198 : // node 1 dofs
199 44424 : ColumnMajorMatrix global_disp(12, 1);
200 44424 : ColumnMajorMatrix global_vel(12, 1);
201 44424 : ColumnMajorMatrix global_disp_old(12, 1);
202 44424 : ColumnMajorMatrix global_vel_old(12, 1);
203 177696 : for (unsigned int i = 0; i < _ndisp; ++i)
204 : {
205 133272 : global_disp(i) =
206 133272 : sol(node[0]->dof_number(nonlinear_sys.number(), _disp_num[i], 0)); // node 0 displacements
207 133272 : global_disp(i + 3) =
208 133272 : sol(node[0]->dof_number(nonlinear_sys.number(), _rot_num[i], 0)); // node 0 rotations
209 133272 : global_disp(i + 6) =
210 133272 : sol(node[1]->dof_number(nonlinear_sys.number(), _disp_num[i], 0)); // node 1 displacements
211 133272 : global_disp(i + 9) =
212 133272 : sol(node[1]->dof_number(nonlinear_sys.number(), _rot_num[i], 0)); // node 1 rotations
213 133272 : global_disp_old(i) = sol_old(
214 : node[0]->dof_number(nonlinear_sys.number(), _disp_num[i], 0)); // node 0 displacements
215 133272 : global_disp_old(i + 3) =
216 133272 : sol_old(node[0]->dof_number(nonlinear_sys.number(), _rot_num[i], 0)); // node 0 rotations
217 133272 : global_disp_old(i + 6) = sol_old(
218 : node[1]->dof_number(nonlinear_sys.number(), _disp_num[i], 0)); // node 1 displacements
219 133272 : global_disp_old(i + 9) =
220 133272 : sol_old(node[1]->dof_number(nonlinear_sys.number(), _rot_num[i], 0)); // node 1 rotations
221 133272 : global_vel(i) = aux_sol(node[0]->dof_number(
222 : aux.number(), _vel_num[i], 0)); // node 0 velocities; rotational velocity terms remain zero
223 133272 : global_vel(i + 6) = aux_sol_old(node[1]->dof_number(
224 : aux.number(), _vel_num[i], 0)); // node 1 velocities; rotational velocity terms remain zero
225 133272 : global_vel_old(i) = aux_sol_old(node[0]->dof_number(
226 : aux.number(), _vel_num[i], 0)); // node 0 velocities; rotational velocity terms remain zero
227 133272 : global_vel_old(i + 6) = aux_sol_old(node[1]->dof_number(
228 : aux.number(), _vel_num[i], 0)); // node 1 velocities; rotational velocity terms remain zero
229 : }
230 :
231 : // Converting global deformations and deformation rates to the isolator
232 : // basic system.
233 44424 : _local_disp[_qp].reshape(12, 1);
234 44424 : _basic_disp[_qp].reshape(6, 1);
235 44424 : _basic_vel[_qp].reshape(6, 1);
236 44424 : _basic_disp_old[_qp].reshape(6, 1);
237 44424 : _basic_vel_old[_qp].reshape(6, 1);
238 :
239 44424 : _local_disp[_qp] = _total_gl[_qp] * global_disp;
240 44424 : _basic_disp[_qp] = _total_lb[_qp] * _total_gl[_qp] * global_disp;
241 44424 : _basic_vel[_qp] = _total_lb[_qp] * _total_gl[_qp] * global_vel;
242 44424 : _basic_disp_old[_qp] = _total_lb[_qp] * _total_gl[_qp] * global_disp_old;
243 88848 : _basic_vel_old[_qp] = _total_lb[_qp] * _total_gl[_qp] * global_vel_old;
244 44424 : }
|
