Line data Source code
1 : // MASTODON includes
2 : #include "LinearSpring.h"
3 :
4 : // MOOSE includes
5 : #include "MooseMesh.h"
6 : #include "Assembly.h"
7 : #include "NonlinearSystem.h"
8 : #include "MooseVariable.h"
9 : #include "RankTwoTensor.h"
10 :
11 : // libmesh includes
12 : #include "libmesh/quadrature.h"
13 : #include "libmesh/utility.h"
14 :
15 : registerMooseObject("MastodonApp", LinearSpring);
16 :
17 : InputParameters
18 24 : LinearSpring::validParams()
19 : {
20 24 : InputParameters params = Material::validParams();
21 24 : params.addClassDescription(
22 : "Compute the deformations, forces and stiffness matrix of a two-noded spring element.");
23 48 : params.addRequiredCoupledVar(
24 : "rotations",
25 : "The rotation variables appropriate for the simulation geometry and coordinate system.");
26 48 : params.addRequiredCoupledVar(
27 : "displacements",
28 : "The displacement variables appropriate for the simulation geometry and coordinate system.");
29 48 : params.addRequiredParam<RealGradient>("y_orientation",
30 : "Orientation of the y direction along "
31 : "which, Ky is provided. This should be "
32 : "perpendicular to the axis of the spring.");
33 48 : params.addRequiredCoupledVar("kx", "Axial stiffness of the spring");
34 48 : params.addRequiredCoupledVar("ky", "Shear stiffness in the y direction of the spring.");
35 48 : params.addRequiredCoupledVar("kz", "Shear stiffness in the z direction of the spring.");
36 48 : params.addRequiredCoupledVar("krx", "Torsional stiffness of the spring.");
37 48 : params.addRequiredCoupledVar("kry", "Rotational stiffness in the y direction of the spring.");
38 48 : params.addRequiredCoupledVar("krz", "Rotational stiffness in the z direction of the spring.");
39 48 : params.set<MooseEnum>("constant_on") = "ELEMENT"; // set _qp to 0
40 24 : return params;
41 0 : }
42 :
43 18 : LinearSpring::LinearSpring(const InputParameters & parameters)
44 : : Material(parameters),
45 18 : _nrot(coupledComponents("rotations")),
46 18 : _ndisp(coupledComponents("displacements")),
47 18 : _rot_num(3),
48 18 : _disp_num(3),
49 18 : _deformations(declareProperty<ColumnMajorMatrix>("deformations")),
50 18 : _rotations(declareProperty<ColumnMajorMatrix>("rotations")),
51 18 : _kx(coupledValue("kx")),
52 18 : _ky(coupledValue("ky")),
53 18 : _kz(coupledValue("kz")),
54 18 : _krx(coupledValue("krx")),
55 18 : _kry(coupledValue("kry")),
56 18 : _krz(coupledValue("krz")),
57 18 : _spring_forces_global(declareProperty<ColumnMajorMatrix>("global_forces")),
58 18 : _spring_moments_global(declareProperty<ColumnMajorMatrix>("global_moments")),
59 18 : _kdd(declareProperty<ColumnMajorMatrix>("displacement_stiffness_matrix")),
60 18 : _krr(declareProperty<ColumnMajorMatrix>("rotation_stiffness_matrix")),
61 18 : _total_global_to_local_rotation(
62 54 : declareProperty<ColumnMajorMatrix>("total_global_to_local_rotation"))
63 : {
64 : // Checking for consistency between length of the provided displacements and rotations vector
65 18 : if (_ndisp != _nrot)
66 0 : mooseError("LinearSpring: The number of variables supplied in 'displacements' "
67 : "and 'rotations' input parameters must be equal.");
68 :
69 : // Fetch coupled variables and gradients (as stateful properties if necessary)
70 72 : for (unsigned int i = 0; i < _ndisp; ++i)
71 : {
72 108 : MooseVariable * disp_variable = getVar("displacements", i);
73 54 : _disp_num[i] = disp_variable->number(); // Displacement variable numbers in MOOSE
74 :
75 108 : MooseVariable * rot_variable = getVar("rotations", i);
76 54 : _rot_num[i] = rot_variable->number(); // Rotation variable numbers in MOOSE
77 : }
78 18 : }
79 :
80 : void
81 0 : LinearSpring::initQpStatefulProperties()
82 : {
83 0 : _deformations[_qp].reshape(3, 1);
84 0 : _rotations[_qp].reshape(3, 1);
85 0 : _spring_forces_global[_qp].reshape(3, 1);
86 0 : _spring_forces_global[_qp].zero();
87 0 : _spring_moments_global[_qp].reshape(3, 1);
88 0 : _kdd[_qp].reshape(3, 3);
89 0 : _krr[_qp].reshape(3, 3);
90 0 : _total_global_to_local_rotation[_qp].reshape(3, 3);
91 0 : _original_global_to_local_rotation.reshape(3, 3);
92 0 : _global_disp0.reshape(3, 1);
93 0 : _global_disp1.reshape(3, 1);
94 0 : _global_rot0.reshape(3, 1);
95 0 : _global_rot1.reshape(3, 1);
96 0 : }
97 :
98 : void
99 19200 : LinearSpring::computeQpProperties()
100 : {
101 : // Compute initial orientation and length of the spring in global coordinate system
102 : // Fetch the two nodes of the link element
103 : std::vector<const Node *> node;
104 57600 : for (unsigned int i = 0; i < 2; ++i)
105 38400 : node.push_back(_current_elem->node_ptr(i));
106 : RealGradient x_orientation;
107 76800 : for (unsigned int i = 0; i < _ndisp; ++i)
108 57600 : x_orientation(i) = (*node[1])(i) - (*node[0])(i);
109 19200 : x_orientation /= x_orientation.norm(); // Normalizing with length to get orientation
110 :
111 : // Get y orientation of the spring in global coordinate system
112 38400 : RealGradient y_orientation = getParam<RealGradient>("y_orientation");
113 : Real dot = x_orientation * y_orientation;
114 :
115 : // Check if x and y orientations are perpendicular
116 19200 : if (abs(dot) > 1e-4)
117 0 : mooseError("Error in LinearSpring: y_orientation should be perpendicular to "
118 : "the axis of the beam.");
119 :
120 : // Calculate z orientation in the global coordinate system as a cross product of the x and y
121 : // orientations
122 19200 : RealGradient z_orientation = x_orientation.cross(y_orientation);
123 :
124 : // Calculate the rotation matrix from global to spring local configuration at t = 0
125 19200 : _original_global_to_local_rotation(0, 0) = x_orientation(0);
126 19200 : _original_global_to_local_rotation(0, 1) = x_orientation(1);
127 19200 : _original_global_to_local_rotation(1, 0) = y_orientation(0);
128 19200 : _original_global_to_local_rotation(1, 1) = y_orientation(1);
129 19200 : _original_global_to_local_rotation(0, 2) = x_orientation(2);
130 19200 : _original_global_to_local_rotation(1, 2) = y_orientation(2);
131 19200 : _original_global_to_local_rotation(2, 0) = z_orientation(0);
132 19200 : _original_global_to_local_rotation(2, 1) = z_orientation(1);
133 19200 : _original_global_to_local_rotation(2, 2) = z_orientation(2);
134 :
135 : // _qp = 0
136 19200 : computeDeformations();
137 :
138 : // computing spring forces
139 19200 : computeForces();
140 :
141 : // computing spring stiffness matrix
142 19200 : computeStiffnessMatrix();
143 19200 : }
144 :
145 : void
146 19200 : LinearSpring::computeTotalRotation()
147 : {
148 19200 : _qp = 0;
149 :
150 : // Currently this forumation is limited to small deformations in the spring,
151 : // namely, it is assumed that there are no rigid body rotations in the spring,
152 : // and that the total rotation matrix from global to local coordinates
153 : // (calculated below) remains the same as the one at t = 0 throughout the
154 : // duration of the simulation.
155 19200 : _total_global_to_local_rotation[_qp] = _original_global_to_local_rotation;
156 19200 : }
157 :
158 : void
159 19200 : LinearSpring::computeDeformations()
160 : {
161 : // fetch the two end nodes for _current_elem
162 : std::vector<const Node *> node;
163 57600 : for (unsigned int i = 0; i < 2; ++i)
164 38400 : node.push_back(_current_elem->node_ptr(i));
165 :
166 : // Fetch the solution for the two end nodes at time t
167 19200 : NonlinearSystemBase & nonlinear_sys = _fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0);
168 19200 : const NumericVector<Number> & sol = *nonlinear_sys.currentSolution();
169 :
170 : // Calculating global displacements and rotations at the end nodes
171 76800 : for (unsigned int i = 0; i < _ndisp; ++i)
172 : {
173 57600 : _global_disp0(i) = sol(node[0]->dof_number(nonlinear_sys.number(), _disp_num[i], 0));
174 57600 : _global_disp1(i) = sol(node[1]->dof_number(nonlinear_sys.number(), _disp_num[i], 0));
175 57600 : _global_rot0(i) = sol(node[0]->dof_number(nonlinear_sys.number(), _rot_num[i], 0));
176 57600 : _global_rot1(i) = sol(node[1]->dof_number(nonlinear_sys.number(), _rot_num[i], 0));
177 : }
178 :
179 : // Convert spring nodal displacements and rotations from global coordinate system to local
180 : // coordinate system First, compute total rotation
181 19200 : computeTotalRotation();
182 19200 : _local_disp0 = _total_global_to_local_rotation[_qp] * _global_disp0;
183 19200 : _local_disp1 = _total_global_to_local_rotation[_qp] * _global_disp1;
184 19200 : _local_rot0 = _total_global_to_local_rotation[_qp] * _global_rot0;
185 19200 : _local_rot1 = _total_global_to_local_rotation[_qp] * _global_rot1;
186 :
187 : // Calculating spring deformations and rotations in the local
188 : // coordinate system. Deformations and rotations are assumed to be constant
189 : // through the length of the spring.
190 19200 : _deformations[_qp] = _local_disp1 - _local_disp0;
191 19200 : _rotations[_qp] = _local_rot1 - _local_rot0;
192 19200 : }
193 :
194 : void
195 19200 : LinearSpring::computeForces()
196 : // spring forces = Kdd * deformations
197 : // spring moments = Krr * rotations
198 : {
199 : // forces
200 19200 : _spring_forces_local(0) = _kx[_qp] * _deformations[_qp](0);
201 19200 : _spring_forces_local(1) = _ky[_qp] * _deformations[_qp](1);
202 19200 : _spring_forces_local(2) = _kz[_qp] * _deformations[_qp](2);
203 : // convert local forces to global
204 19200 : _spring_forces_global[_qp] =
205 19200 : _total_global_to_local_rotation[_qp].transpose() * _spring_forces_local;
206 :
207 : // moments
208 19200 : _spring_moments_local(0) = _krx[_qp] * _rotations[_qp](0);
209 19200 : _spring_moments_local(1) = _kry[_qp] * _rotations[_qp](1);
210 19200 : _spring_moments_local(2) = _krz[_qp] * _rotations[_qp](2);
211 : // convert local moments to global
212 19200 : _spring_moments_global[_qp] =
213 19200 : _total_global_to_local_rotation[_qp].transpose() * _spring_moments_local;
214 19200 : }
215 :
216 : void
217 19200 : LinearSpring::computeStiffnessMatrix()
218 : {
219 : // The stiffness matrix is of the form
220 : // | kdd kdr |
221 : // | krd krr |
222 : // where kdd, krr, krd and kdr are all RankTwoTensors (3x3 matrix) and
223 : // matrix symmetry is assumed, namely, krd = kdr.transpose()
224 : // This implementation of the spring element has only diagonal stiffness
225 : // terms. Therefore, Kdr and Krd are zero.
226 :
227 : // calculating deformation stiffnesses
228 19200 : _kdd_local(0, 0) = _kx[_qp];
229 19200 : _kdd_local(1, 1) = _ky[_qp];
230 19200 : _kdd_local(2, 2) = _kz[_qp];
231 : // convert stiffness matrix from local to global
232 19200 : _kdd[_qp] = _total_global_to_local_rotation[_qp].transpose() * _kdd_local *
233 19200 : _total_global_to_local_rotation[_qp];
234 :
235 : // calculating rotational stiffness
236 19200 : _krr_local(0, 0) = _krx[_qp];
237 19200 : _krr_local(1, 1) = _kry[_qp];
238 19200 : _krr_local(2, 2) = _krz[_qp];
239 : // convert stiffness matrix from local to global
240 19200 : _krr[_qp] = _total_global_to_local_rotation[_qp].transpose() * _krr_local *
241 19200 : _total_global_to_local_rotation[_qp];
242 19200 : }
|
