1 // =============================================================================
2 // PROJECT CHRONO - http://projectchrono.org
3 //
4 // Copyright (c) 2014 projectchrono.org
5 // All rights reserved.
6 //
7 // Use of this source code is governed by a BSD-style license that can be found
8 // in the LICENSE file at the top level of the distribution and at
9 // http://projectchrono.org/license-chrono.txt.
10 //
11 // =============================================================================
12 // Authors: Alessandro Tasora, Radu Serban
13 // =============================================================================
14
15 //#define BEAM_VERBOSE
16
17 #include "chrono/fea/ChElementBeamTaperedTimoshenko.h"
18
19 namespace chrono {
20 namespace fea {
21
ChElementBeamTaperedTimoshenko()22 ChElementBeamTaperedTimoshenko::ChElementBeamTaperedTimoshenko()
23 : q_refrotA(QUNIT),
24 q_refrotB(QUNIT),
25 q_element_abs_rot(QUNIT),
26 q_element_ref_rot(QUNIT),
27 force_symmetric_stiffness(false),
28 disable_corotate(false),
29 use_geometric_stiffness(true),
30 use_Rc(true),
31 use_Rs(true) {
32 nodes.resize(2);
33
34 Km.setZero(this->GetNdofs(), this->GetNdofs());
35 Kg.setZero(this->GetNdofs(), this->GetNdofs());
36 M.setZero(this->GetNdofs(), this->GetNdofs());
37 Rm.setZero(this->GetNdofs(), this->GetNdofs());
38 Ri.setZero(this->GetNdofs(), this->GetNdofs());
39 Ki.setZero(this->GetNdofs(), this->GetNdofs());
40
41 T.setZero(this->GetNdofs(), this->GetNdofs());
42 Rs.setIdentity(6, 6);
43 Rc.setIdentity(6, 6);
44 }
45
SetNodes(std::shared_ptr<ChNodeFEAxyzrot> nodeA,std::shared_ptr<ChNodeFEAxyzrot> nodeB)46 void ChElementBeamTaperedTimoshenko::SetNodes(std::shared_ptr<ChNodeFEAxyzrot> nodeA,
47 std::shared_ptr<ChNodeFEAxyzrot> nodeB) {
48 assert(nodeA);
49 assert(nodeB);
50
51 nodes[0] = nodeA;
52 nodes[1] = nodeB;
53 std::vector<ChVariables*> mvars;
54 mvars.push_back(&nodes[0]->Variables());
55 mvars.push_back(&nodes[1]->Variables());
56 Kmatr.SetVariables(mvars);
57 }
58
ShapeFunctionsTimoshenko(ShapeFunctionGroup & NN,double eta)59 void ChElementBeamTaperedTimoshenko::ShapeFunctionsTimoshenko(ShapeFunctionGroup& NN, double eta) {
60 // eta = 2 * x/L;
61 // x = (-L/2, L/2), hence eta = (-1, 1)
62
63 double L = this->length;
64 double LL = L * L;
65 double LLL = LL * L;
66 double phiy = this->tapered_section->GetAverageSectionParameters()->phiy;
67 double phiz = this->tapered_section->GetAverageSectionParameters()->phiz;
68
69 double eta2 = eta * eta;
70 double eta3 = eta2 * eta;
71
72 double ay = 1. / (1. + phiy);
73 double by = phiy / (1. + phiy);
74 double az = 1. / (1. + phiz);
75 double bz = phiz / (1. + phiz);
76
77 // axial, torsion
78 double Nx1 = 1. / 2. * (1 - eta);
79 double Nx2 = 1. / 2. * (1 + eta);
80 // shear
81 double Nsu1 = 1. / 2. * (1 - eta);
82 double Nsu2 = 1. / 2. * (1 + eta);
83 double Nsr = -1. / 4. * (1 + eta) * L;
84 // bending
85 double Nbu1 = 1. / 4. * (eta3 - 3. * eta + 2.);
86 double Nbu2 = 1. / 4. * (-eta3 + 3. * eta + 2.);
87 double Nbr1 = 1. / 8. * (eta3 - eta2 - eta + 1.) * L;
88 double Nbr2 = 1. / 8. * (-eta2 + 2. * eta + 3.) * L;
89 double Nbr3 = 1. / 8. * (eta3 + eta2 - eta - 1.) * L;
90 double Nbr4 = 1. / 8. * (eta2 + 2. * eta + 1.) * L;
91
92 // shape functions, reduced from three-dimensional shape functions.
93 // hence, only valid for the points at centerline
94 ShapeFunctionN N;
95 N.setZero();
96 // x
97 N(0, 0) = Nx1;
98 N(0, 6) = Nx2;
99 // y
100 N(1, 1) = Nbu1 * ay + Nsu1 * by;
101 N(1, 5) = Nbr1 * ay + Nbr2 * by + Nsr * by;
102 N(1, 7) = Nbu2 * ay + Nsu2 * by;
103 N(1, 11) = Nbr3 * ay + Nbr4 * by + Nsr * by;
104 // z
105 N(2, 2) = Nbu1 * az + Nsu1 * bz;
106 N(2, 4) = -Nbr1 * az - Nbr2 * bz - Nsr * bz;
107 N(2, 8) = Nbu2 * az + Nsu2 * bz;
108 N(2, 10) = -Nbr3 * az - Nbr4 * bz - Nsr * bz;
109 // rotx
110 N(3, 3) = Nx1;
111 N(3, 9) = Nx2;
112
113 // Variables named with initial 'k' means block of shape function matrix
114 // Variables named with initial 'd' means derivatives respect to x
115 // 'd' --> first derivative
116 // 'dd' --> second derivative
117 // 'ddd' --> third derivative
118
119 // some important blocks from shape function matrix
120 // bending y
121 SFBlock kNby;
122 kNby.setZero();
123 kNby << Nbu1 * ay, Nbr1 * ay + Nbr2 * by, Nbu2 * ay, Nbr3 * ay + Nbr4 * by;
124 // shear y
125 SFBlock kNsy;
126 kNsy.setZero();
127 kNsy << Nsu1 * by, Nsr * by, Nsu2 * by, Nsr * by;
128 // bending z
129 SFBlock kNbz;
130 kNbz.setZero();
131 kNbz << Nbu1 * az, -Nbr1 * az - Nbr2 * bz, Nbu2 * az, -Nbr3 * az - Nbr4 * bz;
132 // shear z
133 SFBlock kNsz;
134 kNsz.setZero();
135 kNsz << Nsu1 * bz, -Nsr * bz, Nsu2 * bz, -Nsr * bz;
136 // axial, torsion
137 ChMatrixNM<double, 1, 2> kNx;
138 kNx.setZero();
139 kNx << Nx1, Nx2;
140
141 // first derivatives
142 double dNx1_dx = -1. / L;
143 double dNx2_dx = 1. / L;
144 double dNsu1_dx = -1. / L;
145 double dNsu2_dx = 1. / L;
146 double dNsr_dx = -1. / 2.;
147
148 double dNbu1_dx = 1. / (2. * L) * (3. * eta2 - 3.);
149 double dNbu2_dx = 1. / (2. * L) * (-3. * eta2 + 3.);
150
151 double dNbr1_dx = 1. / 4. * (3. * eta2 - 2. * eta - 1.);
152 double dNbr2_dx = 1. / 2. * (1. - eta);
153 double dNbr3_dx = 1. / 4. * (3. * eta2 + 2. * eta - 1.);
154 double dNbr4_dx = 1. / 2. * (1. + eta);
155
156 SFBlock dkNby_dx;
157 dkNby_dx.setZero();
158 dkNby_dx << dNbu1_dx * ay, dNbr1_dx * ay + dNbr2_dx * by, dNbu2_dx * ay, dNbr3_dx * ay + dNbr4_dx * by;
159
160 SFBlock dkNsy_dx;
161 dkNsy_dx.setZero();
162 dkNsy_dx << dNsu1_dx * by, dNsr_dx * by, dNsu2_dx * by, dNsr_dx * by;
163
164 SFBlock dkNbz_dx;
165 dkNbz_dx.setZero();
166 dkNbz_dx << dNbu1_dx * az, -dNbr1_dx * az - dNbr2_dx * bz, dNbu2_dx * az, -dNbr3_dx * az - dNbr4_dx * bz;
167
168 SFBlock dkNsz_dx;
169 dkNsz_dx.setZero();
170 dkNsz_dx << dNsu1_dx * bz, -dNsr_dx * bz, dNsu2_dx * bz, -dNsr_dx * bz;
171
172 ChMatrixNM<double, 1, 2> dkNx_dx;
173 dkNx_dx.setZero();
174 dkNx_dx << dNx1_dx, dNx2_dx;
175
176 // Continue to fill in the shape functions to finish it.
177 // roty
178 N(4, 2) = -dNbu1_dx * az; // - dNsu1_dx * bz;
179 N(4, 4) = dNbr1_dx * az + dNbr2_dx * bz; // + dNsr_dx * bz;
180 N(4, 8) = -dNbu2_dx * az; // - dNsu2_dx * bz;
181 N(4, 10) = dNbr3_dx * az + dNbr4_dx * bz; // + dNsr_dx * bz;
182 // rotz
183 N(5, 1) = dNbu1_dx * ay; // + dNsu1_dx * by;
184 N(5, 5) = dNbr1_dx * ay + dNbr2_dx * by; // + dNsr_dx * by;
185 N(5, 7) = dNbu2_dx * ay; // + dNsu2_dx * by;
186 N(5, 11) = dNbr3_dx * ay + dNbr4_dx * by; // + dNsr_dx * by;
187
188 // Second derivatives
189 double ddNbu1_dx = 6. * eta / LL;
190 double ddNbu2_dx = -6. * eta / LL;
191
192 double ddNbr1_dx = 1. / L * (3. * eta - 1.);
193 double ddNbr2_dx = -1. / L;
194 double ddNbr3_dx = 1. / L * (3. * eta + 1.);
195 double ddNbr4_dx = 1. / L;
196
197 SFBlock ddkNby_dx;
198 ddkNby_dx.setZero();
199 ddkNby_dx << ddNbu1_dx * ay, ddNbr1_dx * ay + ddNbr2_dx * by, ddNbu2_dx * ay, ddNbr3_dx * ay + ddNbr4_dx * by;
200
201 SFBlock ddkNbz_dx;
202 ddkNbz_dx.setZero();
203 ddkNbz_dx << ddNbu1_dx * az, -ddNbr1_dx * az - ddNbr2_dx * bz, ddNbu2_dx * az, -ddNbr3_dx * az - ddNbr4_dx * bz;
204
205 // Third derivatives
206 double dddNbu1_dx = 12. / LLL;
207 double dddNbu2_dx = -12. / LLL;
208
209 double dddNbr1_dx = 6. / LL;
210 double dddNbr2_dx = 0.;
211 double dddNbr3_dx = 6. / LL;
212 double dddNbr4_dx = 0.;
213
214 SFBlock dddkNby_dx;
215 dddkNby_dx.setZero();
216 dddkNby_dx << dddNbu1_dx * ay, dddNbr1_dx * ay + dddNbr2_dx * by, dddNbu2_dx * ay,
217 dddNbr3_dx * ay + dddNbr4_dx * by;
218
219 SFBlock dddkNbz_dx;
220 dddkNbz_dx.setZero();
221 dddkNbz_dx << dddNbu1_dx * az, -dddNbr1_dx * az - dddNbr2_dx * bz, dddNbu2_dx * az,
222 -dddNbr3_dx * az - dddNbr4_dx * bz;
223
224 ShapeFunction5Blocks SFblk = std::make_tuple(kNby, kNsy, kNbz, kNsz, kNx);
225 ShapeFunction5Blocks SFblk1d = std::make_tuple(dkNby_dx, dkNsy_dx, dkNbz_dx, dkNsz_dx, dkNx_dx);
226 ShapeFunction2Blocks SFblk2d = std::make_tuple(ddkNby_dx, ddkNbz_dx);
227 ShapeFunction2Blocks SFblk3d = std::make_tuple(dddkNby_dx, dddkNbz_dx);
228
229 std::get<0>(NN) = N;
230 std::get<1>(NN) = SFblk;
231 std::get<2>(NN) = SFblk1d;
232 std::get<3>(NN) = SFblk2d;
233 std::get<4>(NN) = SFblk3d;
234 }
235
Update()236 void ChElementBeamTaperedTimoshenko::Update() {
237 // parent class update:
238 ChElementGeneric::Update();
239
240 // always keep updated the rotation matrix A:
241 this->UpdateRotation();
242 }
243
UpdateRotation()244 void ChElementBeamTaperedTimoshenko::UpdateRotation() {
245 ChMatrix33<> A0(this->q_element_ref_rot);
246
247 ChMatrix33<> Aabs;
248 if (this->disable_corotate) {
249 Aabs = A0;
250 q_element_abs_rot = q_element_ref_rot;
251 } else {
252 ChVector<> mXele_w = nodes[1]->Frame().GetPos() - nodes[0]->Frame().GetPos();
253 // propose Y_w as absolute dir of the Y axis of A node, removing the effect of Aref-to-A rotation if any:
254 // Y_w = [R Aref->w ]*[R Aref->A ]'*{0,1,0}
255 ChVector<> myele_wA = nodes[0]->Frame().GetRot().Rotate(q_refrotA.RotateBack(ChVector<>(0, 1, 0)));
256 // propose Y_w as absolute dir of the Y axis of B node, removing the effect of Bref-to-B rotation if any:
257 // Y_w = [R Bref->w ]*[R Bref->B ]'*{0,1,0}
258 ChVector<> myele_wB = nodes[1]->Frame().GetRot().Rotate(q_refrotB.RotateBack(ChVector<>(0, 1, 0)));
259 // Average the two Y directions to have midpoint torsion (ex -30?torsion A and +30?torsion B= 0?
260 ChVector<> myele_w = (myele_wA + myele_wB).GetNormalized();
261 Aabs.Set_A_Xdir(mXele_w, myele_w);
262 q_element_abs_rot = Aabs.Get_A_quaternion();
263 }
264
265 this->A = A0.transpose() * Aabs;
266 }
267
GetStateBlock(ChVectorDynamic<> & mD)268 void ChElementBeamTaperedTimoshenko::GetStateBlock(ChVectorDynamic<>& mD) {
269 mD.resize(12);
270
271 ChVector<> delta_rot_dir;
272 double delta_rot_angle;
273
274 // Node 0, displacement (in local element frame, corotated back)
275 // d = [Atw]' Xt - [A0w]'X0
276 ChVector<> displ = this->q_element_abs_rot.RotateBack(nodes[0]->Frame().GetPos()) -
277 this->q_element_ref_rot.RotateBack(nodes[0]->GetX0().GetPos());
278 mD.segment(0, 3) = displ.eigen();
279
280 // Node 0, x,y,z small rotations (in local element frame)
281 ChQuaternion<> q_delta0 = q_element_abs_rot.GetConjugate() % nodes[0]->Frame().GetRot() % q_refrotA.GetConjugate();
282 // note, for small incremental rotations this is opposite of ChNodeFEAxyzrot::VariablesQbIncrementPosition
283 q_delta0.Q_to_AngAxis(delta_rot_angle, delta_rot_dir);
284
285 if (delta_rot_angle > CH_C_PI)
286 delta_rot_angle -= CH_C_2PI; // no 0..360 range, use -180..+180
287
288 mD.segment(3, 3) = delta_rot_angle * delta_rot_dir.eigen();
289
290 // Node 1, displacement (in local element frame, corotated back)
291 // d = [Atw]' Xt - [A0w]'X0
292 displ = this->q_element_abs_rot.RotateBack(nodes[1]->Frame().GetPos()) -
293 this->q_element_ref_rot.RotateBack(nodes[1]->GetX0().GetPos());
294 mD.segment(6, 3) = displ.eigen();
295
296 // Node 1, x,y,z small rotations (in local element frame)
297 ChQuaternion<> q_delta1 = q_element_abs_rot.GetConjugate() % nodes[1]->Frame().GetRot() % q_refrotB.GetConjugate();
298 // note, for small incremental rotations this is opposite of ChNodeFEAxyzrot::VariablesQbIncrementPosition
299 q_delta1.Q_to_AngAxis(delta_rot_angle, delta_rot_dir);
300
301 if (delta_rot_angle > CH_C_PI)
302 delta_rot_angle -= CH_C_2PI; // no 0..360 range, use -180..+180
303
304 mD.segment(9, 3) = delta_rot_angle * delta_rot_dir.eigen();
305 }
306
GetField_dt(ChVectorDynamic<> & mD_dt)307 void ChElementBeamTaperedTimoshenko::GetField_dt(ChVectorDynamic<>& mD_dt) {
308 mD_dt.resize(12);
309
310 // Node 0, velocity (in local element frame, corotated back by A' )
311 mD_dt.segment(0, 3) = q_element_abs_rot.RotateBack(nodes[0]->Frame().GetPos_dt()).eigen();
312
313 // Node 0, x,y,z ang.velocity (in local element frame, corotated back by A' )
314 mD_dt.segment(3, 3) = q_element_abs_rot.RotateBack(nodes[0]->Frame().GetWvel_par()).eigen();
315
316 // Node 1, velocity (in local element frame, corotated back by A' )
317 mD_dt.segment(6, 3) = q_element_abs_rot.RotateBack(nodes[1]->Frame().GetPos_dt()).eigen();
318
319 // Node 1, x,y,z ang.velocity (in local element frame, corotated back by A' )
320 mD_dt.segment(9, 3) = q_element_abs_rot.RotateBack(nodes[1]->Frame().GetWvel_par()).eigen();
321 }
322
GetField_dtdt(ChVectorDynamic<> & mD_dtdt)323 void ChElementBeamTaperedTimoshenko::GetField_dtdt(ChVectorDynamic<>& mD_dtdt) {
324 mD_dtdt.resize(12);
325
326 // Node 0, accelaration (in local element frame, corotated back by A' )
327 mD_dtdt.segment(0, 3) = q_element_abs_rot.RotateBack(nodes[0]->Frame().GetPos_dtdt()).eigen();
328
329 // Node 0, x,y,z ang.accelaration (in local element frame, corotated back by A' )
330 mD_dtdt.segment(3, 3) = q_element_abs_rot.RotateBack(nodes[0]->Frame().GetWacc_par()).eigen();
331
332 // Node 1, accelaration (in local element frame, corotated back by A' )
333 mD_dtdt.segment(6, 3) = q_element_abs_rot.RotateBack(nodes[1]->Frame().GetPos_dtdt()).eigen();
334
335 // Node 1, x,y,z ang.accelaration (in local element frame, corotated back by A' )
336 mD_dtdt.segment(9, 3) = q_element_abs_rot.RotateBack(nodes[1]->Frame().GetWacc_par()).eigen();
337 }
338
ComputeTransformMatrix()339 void ChElementBeamTaperedTimoshenko::ComputeTransformMatrix() {
340 double alpha1 = this->tapered_section->GetSectionA()->GetSectionRotation();
341 double Cy1 = this->tapered_section->GetSectionA()->GetCentroidY();
342 double Cz1 = this->tapered_section->GetSectionA()->GetCentroidZ();
343 double Sy1 = this->tapered_section->GetSectionA()->GetShearCenterY();
344 double Sz1 = this->tapered_section->GetSectionA()->GetShearCenterZ();
345
346 double alpha2 = this->tapered_section->GetSectionB()->GetSectionRotation();
347 double Cy2 = this->tapered_section->GetSectionB()->GetCentroidY();
348 double Cz2 = this->tapered_section->GetSectionB()->GetCentroidZ();
349 double Sy2 = this->tapered_section->GetSectionB()->GetShearCenterY();
350 double Sz2 = this->tapered_section->GetSectionB()->GetShearCenterZ();
351
352 double L = this->length;
353
354 // In case the section is rotated:
355 ChMatrix33<> RotsectA;
356 RotsectA.Set_A_Rxyz(ChVector<>(alpha1, 0, 0));
357 ChMatrix33<> RotsectB;
358 RotsectB.Set_A_Rxyz(ChVector<>(alpha2, 0, 0));
359 ChMatrixNM<double, 12, 12> Rotsect;
360 Rotsect.setZero();
361 Rotsect.block<3, 3>(0, 0) = RotsectA;
362 Rotsect.block<3, 3>(3, 3) = RotsectA;
363 Rotsect.block<3, 3>(6, 6) = RotsectB;
364 Rotsect.block<3, 3>(9, 9) = RotsectB;
365
366 // In case the section has a centroid displacement:
367 // double dCx = 0;
368 double dCy = Cy2 - Cy1;
369 double dCz = Cz2 - Cz1;
370 double LN = L;
371 double LG = pow(LN * LN - dCy * dCy - dCz * dCz, 0.5);
372 double LA = LG;
373 double LB = pow(LG * LG + dCy * dCy, 0.5);
374 ChMatrix33<> rc;
375 rc.setIdentity();
376 rc(0, 0) = LA / LN;
377 rc(0, 1) = dCy / LB;
378 rc(0, 2) = LA * dCz / (LN * LB);
379 rc(1, 0) = -dCy / LN;
380 rc(1, 1) = LA / LB;
381 rc(1, 2) = -dCy * dCz / (LN * LB);
382 rc(2, 0) = -dCz / LN;
383 rc(2, 1) = 0.0;
384 rc(2, 2) = LB / LN;
385 // ChMatrixNM<double, 6, 6> Rc;
386 // Rc.setIdentity();
387 if (use_Rc) {
388 this->Rc.block<3, 3>(0, 0) = rc;
389 this->Rc.block<3, 3>(3, 3) = rc;
390 }
391
392 ChMatrixNM<double, 6, 6> Tc1;
393 Tc1.setIdentity();
394 Tc1(0, 4) = Cz1;
395 Tc1(0, 5) = -Cy1;
396 Tc1(1, 3) = -Cz1;
397 Tc1(2, 3) = Cy1;
398
399 ChMatrixNM<double, 6, 6> Tc2;
400 Tc2.setIdentity();
401 Tc2(0, 4) = Cz2;
402 Tc2(0, 5) = -Cy2;
403 Tc2(1, 3) = -Cz2;
404 Tc2(2, 3) = Cy2;
405
406 ChMatrixNM<double, 12, 12> Tc;
407 Tc.setIdentity();
408 Tc.block<6, 6>(0, 0) = Tc1 * Rc;
409 Tc.block<6, 6>(6, 6) = Tc2 * Rc;
410
411 // In case the section has a shear center displacement:
412 // The shear center offset is respect to the centerline of beam element.
413 /*Sy1 = Sy1 - Cy1; // Unnecessary to do this substraction
414 Sz1 = Sz1 - Cz1;
415 Sy2 = Sy2 - Cy2;
416 Sz2 = Sz2 - Cz2;*/
417 double dSy = Sy1 - Sy2;
418 double dSz = Sz1 - Sz2;
419
420 double Lsy = pow(LG * LG + dSy * dSy, 0.5);
421 double Lsyz = pow(LG * LG + dSy * dSy + dSz * dSz, 0.5);
422 double C1 = Lsyz / LG;
423 double C2 = dSy * Lsyz / (LG * Lsy);
424 double C3 = dSz / Lsy;
425
426 ChMatrixNM<double, 6, 6> Ts1;
427 Ts1.setIdentity();
428 Ts1(1, 3) = -Sz1;
429 Ts1(2, 3) = Sy1;
430
431 ChMatrixNM<double, 6, 6> Ts2;
432 Ts2.setIdentity();
433 Ts2(1, 3) = -Sz2;
434 Ts2(2, 3) = Sy2;
435
436 // ChMatrixNM<double, 6, 6> Rs;
437 // Rs.setIdentity();
438 if (use_Rs) {
439 this->Rs(3, 3) = C1;
440 this->Rs(4, 3) = C2;
441 this->Rs(5, 3) = C3;
442 }
443
444 ChMatrixNM<double, 12, 12> Ts;
445 Ts.setIdentity();
446 Ts.block<6, 6>(0, 0) = Rs * Ts1;
447 Ts.block<6, 6>(6, 6) = Rs * Ts2;
448
449 this->T = Rotsect * Ts * Tc;
450 }
451
ComputeTransformMatrixAtPoint(ChMatrixDynamic<> & mT,const double eta)452 void ChElementBeamTaperedTimoshenko::ComputeTransformMatrixAtPoint(ChMatrixDynamic<>& mT, const double eta) {
453 mT.resize(6, 6);
454
455 double alpha1 = this->tapered_section->GetSectionA()->GetSectionRotation();
456 double Cy1 = this->tapered_section->GetSectionA()->GetCentroidY();
457 double Cz1 = this->tapered_section->GetSectionA()->GetCentroidZ();
458 double Sy1 = this->tapered_section->GetSectionA()->GetShearCenterY();
459 double Sz1 = this->tapered_section->GetSectionA()->GetShearCenterZ();
460
461 double alpha2 = this->tapered_section->GetSectionB()->GetSectionRotation();
462 double Cy2 = this->tapered_section->GetSectionB()->GetCentroidY();
463 double Cz2 = this->tapered_section->GetSectionB()->GetCentroidZ();
464 double Sy2 = this->tapered_section->GetSectionB()->GetShearCenterY();
465 double Sz2 = this->tapered_section->GetSectionB()->GetShearCenterZ();
466
467 // double L = this->length;
468
469 // The shear center offset is respect to the centerline of beam element.
470 /*Sy1 = Sy1 - Cy1; // Unnecessary to do this substraction
471 Sz1 = Sz1 - Cz1;
472 Sy2 = Sy2 - Cy2;
473 Sz2 = Sz2 - Cz2;*/
474
475 // calculate the orientation angle and centers by linear interpolation
476 double Nx1 = (1. / 2.) * (1 - eta);
477 double Nx2 = (1. / 2.) * (1 + eta);
478 double alpha = Nx1 * alpha1 + Nx2 * alpha2;
479 double Cy = Nx1 * Cy1 + Nx2 * Cy2;
480 double Cz = Nx1 * Cz1 + Nx2 * Cz2;
481 double Sy = Nx1 * Sy1 + Nx2 * Sy2;
482 double Sz = Nx1 * Sz1 + Nx2 * Sz2;
483
484 // In case the section is rotated:
485 ChMatrix33<> RotsectA;
486 RotsectA.Set_A_Rxyz(ChVector<>(alpha, 0, 0));
487 ChMatrixNM<double, 6, 6> Rotsect;
488 Rotsect.setZero();
489 Rotsect.block<3, 3>(0, 0) = RotsectA;
490 Rotsect.block<3, 3>(3, 3) = RotsectA;
491
492 // In case the section has a centroid displacement:
493 ChMatrixNM<double, 6, 6> Tc;
494 Tc.setIdentity();
495 Tc(0, 4) = Cz;
496 Tc(0, 5) = -Cy;
497 Tc(1, 3) = -Cz;
498 Tc(2, 3) = Cy;
499
500 Tc = Tc * this->Rc;
501
502 // In case the section has a shear center displacement:
503 ChMatrixNM<double, 6, 6> Ts;
504 Ts.setIdentity();
505 Ts(1, 3) = -Sz;
506 Ts(2, 3) = Sy;
507
508 Ts = this->Rs * Ts;
509
510 // the transformation matrix at point eta:
511 mT = Rotsect * Ts * Tc;
512 }
513
ComputeStiffnessMatrix()514 void ChElementBeamTaperedTimoshenko::ComputeStiffnessMatrix() {
515 assert(tapered_section);
516
517 double L = this->length;
518 double LL = L * L;
519 double LLL = LL * L;
520
521 double EA = this->tapered_section->GetAverageSectionParameters()->EA;
522 double GJ = this->tapered_section->GetAverageSectionParameters()->GJ;
523 // double GAyy = this->tapered_section->GetAverageSectionParameters()->GAyy;
524 // double GAzz = this->tapered_section->GetAverageSectionParameters()->GAzz;
525 double EIyy = this->tapered_section->GetAverageSectionParameters()->EIyy;
526 double EIzz = this->tapered_section->GetAverageSectionParameters()->EIzz;
527
528 double phiy = this->tapered_section->GetAverageSectionParameters()->phiy;
529 double phiz = this->tapered_section->GetAverageSectionParameters()->phiz;
530
531 double az = 12.0 * EIzz / (LLL * (1 + phiy));
532 double ay = 12.0 * EIyy / (LLL * (1 + phiz));
533 double cz = 6.0 * EIzz / (LL * (1 + phiy));
534 double cy = 6.0 * EIyy / (LL * (1 + phiz));
535 double ez = (4 + phiy) * EIzz / (L * (1 + phiy));
536 double ey = (4 + phiz) * EIyy / (L * (1 + phiz));
537 double fz = (2 - phiy) * EIzz / (L * (1 + phiy));
538 double fy = (2 - phiz) * EIyy / (L * (1 + phiz));
539
540 Km(0, 0) = EA / L;
541 Km(1, 1) = az;
542 Km(2, 2) = ay;
543 Km(3, 3) = GJ / L;
544 Km(4, 4) = ey;
545 Km(5, 5) = ez;
546 Km(6, 6) = EA / L;
547 Km(7, 7) = az;
548 Km(8, 8) = ay;
549 Km(9, 9) = GJ / L;
550 Km(10, 10) = ey;
551 Km(11, 11) = ez;
552
553 Km(0, 6) = -EA / L;
554 Km(1, 7) = -az;
555 Km(2, 8) = -ay;
556 Km(3, 9) = -GJ / L;
557 Km(4, 10) = fy;
558 Km(5, 11) = fz;
559
560 Km(4, 8) = cy;
561 Km(5, 7) = -cz;
562 Km(1, 11) = cz;
563 Km(2, 10) = -cy;
564
565 Km(1, 5) = cz;
566 Km(2, 4) = -cy;
567 Km(7, 11) = -cz;
568 Km(8, 10) = cy;
569
570 // symmetric part;
571 for (int r = 0; r < 12; r++)
572 for (int c = r + 1; c < 12; c++)
573 Km(c, r) = Km(r, c);
574
575 Km = this->T.transpose() * Km * this->T;
576 }
577
ComputeDampingMatrix()578 void ChElementBeamTaperedTimoshenko::ComputeDampingMatrix() {
579 assert(tapered_section);
580
581 double L = this->length;
582 double LL = L * L;
583 double LLL = LL * L;
584
585 double EA = this->tapered_section->GetAverageSectionParameters()->EA;
586 double GJ = this->tapered_section->GetAverageSectionParameters()->GJ;
587 // double GAyy = this->tapered_section->GetAverageSectionParameters()->GAyy; // not used here
588 // double GAzz = this->tapered_section->GetAverageSectionParameters()->GAzz; // not used here
589 double EIyy = this->tapered_section->GetAverageSectionParameters()->EIyy;
590 double EIzz = this->tapered_section->GetAverageSectionParameters()->EIzz;
591
592 double phiy = this->tapered_section->GetAverageSectionParameters()->phiy;
593 double phiz = this->tapered_section->GetAverageSectionParameters()->phiz;
594
595 double bx2 = pow(this->tapered_section->GetAverageSectionParameters()->rdamping_coeff.bx, 2.0);
596 double by2 = pow(this->tapered_section->GetAverageSectionParameters()->rdamping_coeff.by, 2.0);
597 double bz2 = pow(this->tapered_section->GetAverageSectionParameters()->rdamping_coeff.bz, 2.0);
598 double bt2 = pow(this->tapered_section->GetAverageSectionParameters()->rdamping_coeff.bt, 2.0);
599 double rdamping_alpha = this->tapered_section->GetAverageSectionParameters()->rdamping_coeff.alpha;
600
601 // Correction for the structural damping in the shear deformation, to improve the numerical stability in long-time
602 // simulation. It might be helpful, also possible to be useless at all.
603 double cc2 = pow(this->tapered_section->GetAverageSectionParameters()->artificial_factor_for_shear_damping, 2.0);
604 double ccphiy = (1 + phiy * cc2) / (1. + phiy);
605 double ccphiz = (1 + phiz * cc2) / (1. + phiz);
606 double ccbz = (EIyy / L * phiz * bz2 * (cc2 - 1)) / (1. + phiz);
607 double ccby = (EIzz / L * phiy * by2 * (cc2 - 1)) / (1. + phiy);
608
609 double az = 12.0 * EIzz / (LLL * (1. + phiy));
610 double ay = 12.0 * EIyy / (LLL * (1. + phiz));
611 double cz = 6.0 * EIzz / (LL * (1. + phiy));
612 double cy = 6.0 * EIyy / (LL * (1. + phiz));
613 double ez = (4.0 + phiy) * EIzz / (L * (1. + phiy));
614 double ey = (4.0 + phiz) * EIyy / (L * (1. + phiz));
615 double fz = (2.0 - phiy) * EIzz / (L * (1. + phiy));
616 double fy = (2.0 - phiz) * EIyy / (L * (1. + phiz));
617
618 Rm(0, 0) = EA / L * bx2;
619 Rm(1, 1) = az * by2 * ccphiy;
620 Rm(2, 2) = ay * bz2 * ccphiz;
621 Rm(3, 3) = GJ / L * bt2;
622 Rm(4, 4) = ey * bz2 * ccphiz;
623 Rm(5, 5) = ez * by2 * ccphiy;
624 Rm(6, 6) = EA / L * bx2;
625 Rm(7, 7) = az * by2 * ccphiy;
626 Rm(8, 8) = ay * bz2 * ccphiz;
627 Rm(9, 9) = GJ / L * bt2;
628 Rm(10, 10) = ey * bz2 * ccphiz;
629 Rm(11, 11) = ez * by2 * ccphiy;
630
631 Rm(0, 6) = -EA / L * bx2;
632 Rm(1, 7) = -az * by2 * ccphiy;
633 Rm(2, 8) = -ay * bz2 * ccphiz;
634 Rm(3, 9) = -GJ / L * bt2;
635 Rm(4, 10) = fy * bz2 * ccphiz;
636 Rm(5, 11) = fz * by2 * ccphiy;
637
638 Rm(4, 8) = cy * bz2 * ccphiz;
639 Rm(5, 7) = -cz * by2 * ccphiy;
640 Rm(1, 11) = cz * by2 * ccphiy;
641 Rm(2, 10) = -cy * bz2 * ccphiz;
642
643 Rm(1, 5) = cz * by2 * ccphiy;
644 Rm(2, 4) = -cy * bz2 * ccphiz;
645 Rm(7, 11) = -cz * by2 * ccphiy;
646 Rm(8, 10) = cy * bz2 * ccphiz;
647
648 Rm(4, 4) += -ccbz;
649 Rm(5, 5) += -ccby;
650 Rm(10, 4) += ccbz;
651 Rm(4, 10) += ccbz;
652 Rm(11, 5) += ccby;
653 Rm(5, 11) += ccby;
654 Rm(10, 10) += -ccbz;
655 Rm(11, 11) += -ccby;
656
657 // symmetric part;
658 for (int r = 0; r < 12; r++)
659 for (int c = r + 1; c < 12; c++)
660 Rm(c, r) = Rm(r, c);
661
662 // The stiffness-proportional term.
663 // The geometric stiffness Kg is NOT considered in the damping matrix Rm.
664 Rm = this->T.transpose() * Rm * this->T;
665
666 // The mass-proportional term
667 if (this->tapered_section->GetLumpedMassMatrixType()) {
668 double node_multiplier_fact = 0.5 * this->length;
669 Rm += rdamping_alpha * this->M * node_multiplier_fact;
670 } else {
671 Rm += rdamping_alpha * this->M;
672 }
673 }
674
ComputeGeometricStiffnessMatrix()675 void ChElementBeamTaperedTimoshenko::ComputeGeometricStiffnessMatrix() {
676 assert(tapered_section);
677
678 // Compute the local geometric stiffness matrix Kg without the P multiplication term, that is Kg*(1/P),
679 // so that it is a constant matrix for performance reasons.
680 // We used the analytical values from
681 // [1] "Nonlinear finite element analysis of elastic frames", Kuo Mo, Hsiao Fang, Yu Hou
682 // Computers & Structures Volume 26, Issue 4, 1987, Pages 693-701
683 // For the Reddy or timoshenko more detailed case with higher order terms, look also to:
684 // [2] "A Unified Approach to the Timoshenko Geometric Stiffness Matrix Considering Higher-Order Terms in the
685 // Strain Tensor"
686 // https://www.scielo.br/pdf/lajss/v16n4/1679-7825-lajss-16-04-e185.pdf
687 // Look also at: https://enercalc.com/3d_help/toc161394033.html or in Plesha, Malkus, Cook �Concepts and
688 // Applications of Finite Element Analysis? or in W. McGuire & R.H. Gallagher & R.D. Ziemian, �Matrix Structural
689 // Analysis?
690
691 // double L = this->tapered_section->GetLength();
692 double L = this->length;
693
694 // double P_L = 1. / L;
695 double P6_5L_y = 6. / (5. * L); // optional [2]: ...+ 12*IzA /(L*L*L);
696 double P6_5L_z = 6. / (5. * L); // optional [2]: ...+ 12*IyA /(L*L*L);
697 double P_10_y = 1. / (10.); // optional [2]: ...+ 6*IzA /(L*L);
698 double P_10_z = 1. / (10.); // optional [2]: ...+ 6*IyA /(L*L);
699 double PL2_15_y = 2. * L / (15.); // optional [2]: ...+ 4*IzA /(L);
700 double PL2_15_z = 2. * L / (15.); // optional [2]: ...+ 4*IyA /(L);
701 double PL_30_y = L / (30.); // optional [2]: ...+ 2*IyA /(L);
702 double PL_30_z = L / (30.); // optional [2]: ...+ 2*IyA /(L);
703 /* DONOT use the axial terms:
704 this->Kg(0, 0) = P_L;
705 this->Kg(6, 6) = P_L;
706 this->Kg(0, 6) = -P_L;
707 */
708
709 this->Kg(1, 1) = P6_5L_y;
710 this->Kg(1, 5) = P_10_y;
711 this->Kg(1, 7) = -P6_5L_y;
712 this->Kg(1, 11) = P_10_y;
713
714 this->Kg(2, 2) = P6_5L_z;
715 this->Kg(2, 4) = -P_10_z;
716 this->Kg(2, 8) = -P6_5L_z;
717 this->Kg(2, 10) = -P_10_z;
718
719 this->Kg(4, 4) = PL2_15_y;
720 this->Kg(4, 8) = P_10_y;
721 this->Kg(4, 10) = -PL_30_y;
722
723 this->Kg(5, 5) = PL2_15_z;
724 this->Kg(5, 7) = -P_10_z;
725 this->Kg(5, 11) = -PL_30_z;
726
727 this->Kg(7, 7) = P6_5L_y;
728 this->Kg(7, 11) = -P_10_y;
729
730 this->Kg(8, 8) = P6_5L_z;
731 this->Kg(8, 10) = P_10_y;
732
733 this->Kg(10, 10) = PL2_15_y;
734
735 this->Kg(11, 11) = PL2_15_z;
736
737 // symmetric part;
738 for (int r = 0; r < 12; r++)
739 for (int c = r + 1; c < 12; c++)
740 Kg(c, r) = Kg(r, c);
741
742 Kg = this->T.transpose() * Kg * this->T;
743 }
744
ComputeKiRimatricesLocal(bool inertial_damping,bool inertial_stiffness)745 void ChElementBeamTaperedTimoshenko::ComputeKiRimatricesLocal(bool inertial_damping, bool inertial_stiffness) {
746 assert(tapered_section);
747
748 if (inertial_damping) {
749 ChVector<> mW_A = this->GetNodeA()->GetWvel_loc();
750 ChVector<> mW_B = this->GetNodeB()->GetWvel_loc();
751 ChMatrixNM<double, 12, 12> mH;
752 this->tapered_section->ComputeInertiaDampingMatrix(mH, mW_A, mW_B);
753 this->Ri = 0.5 * this->length * mH;
754 } else {
755 this->Ri.setZero(this->GetNdofs(), this->GetNdofs());
756 }
757
758 if (inertial_stiffness) {
759 // current angular velocity of section of node A, in material frame
760 ChVector<> mWvel_A = this->GetNodeA()->GetWvel_loc();
761 // current angular acceleration of section of node A, in material frame
762 ChVector<> mWacc_A = this->GetNodeA()->GetWacc_loc();
763 // current acceleration of section of node A, in material frame)
764 ChVector<> mXacc_A = this->GetNodeA()->TransformDirectionParentToLocal(this->GetNodeA()->GetPos_dtdt());
765 // current angular velocity of section of node B, in material frame
766 ChVector<> mWvel_B = this->GetNodeB()->GetWvel_loc();
767 // current angular acceleration of section of node B, in material frame
768 ChVector<> mWacc_B = this->GetNodeB()->GetWacc_loc();
769 // current acceleration of section of node B, in material frame
770 ChVector<> mXacc_B = this->GetNodeB()->TransformDirectionParentToLocal(this->GetNodeB()->GetPos_dtdt());
771
772 ChMatrixNM<double, 12, 12> mH;
773 this->tapered_section->ComputeInertiaStiffnessMatrix(mH, mWvel_A, mWacc_A, mXacc_A, mWvel_B, mWacc_B, mXacc_B);
774 this->Ki = 0.5 * this->length * mH;
775 } else {
776 this->Ki.setZero(this->GetNdofs(), this->GetNdofs());
777 }
778 }
779
SetupInitial(ChSystem * system)780 void ChElementBeamTaperedTimoshenko::SetupInitial(ChSystem* system) {
781 assert(tapered_section);
782
783 // Compute rest length, mass:
784 this->length = (nodes[1]->GetX0().GetPos() - nodes[0]->GetX0().GetPos()).Length();
785 this->mass = 0.5 * this->length * this->tapered_section->GetSectionA()->GetMassPerUnitLength() +
786 0.5 * this->length * this->tapered_section->GetSectionB()->GetMassPerUnitLength();
787
788 // Compute initial rotation
789 ChMatrix33<> A0;
790 ChVector<> mXele = nodes[1]->GetX0().GetPos() - nodes[0]->GetX0().GetPos();
791 ChVector<> myele = nodes[0]->GetX0().GetA().Get_A_Yaxis();
792 A0.Set_A_Xdir(mXele, myele);
793 q_element_ref_rot = A0.Get_A_quaternion();
794
795 // Compute local mass matrix
796 // It could be lumped or consistent mass matrix, depends on SetLumpedMassMatrix(true/false)
797 // If it is lumped mass matrix, you need to multiple 0.5 * length to obtain the final mass matrix
798 // For consistent mass matrix, don't need to multiple anything.
799 this->tapered_section->ComputeInertiaMatrix(this->M);
800
801 // Compute transformation matrix
802 ComputeTransformMatrix();
803
804 // Compute local stiffness matrix:
805 ComputeStiffnessMatrix();
806
807 // Compute local geometric stiffness matrix normalized by pull force P: Kg/P
808 ComputeGeometricStiffnessMatrix();
809
810 // Compute local damping matrix:
811 ComputeDampingMatrix();
812 }
813
ComputeKRMmatricesGlobal(ChMatrixRef H,double Kfactor,double Rfactor,double Mfactor)814 void ChElementBeamTaperedTimoshenko::ComputeKRMmatricesGlobal(ChMatrixRef H,
815 double Kfactor,
816 double Rfactor,
817 double Mfactor) {
818 assert((H.rows() == 12) && (H.cols() == 12));
819 assert(tapered_section);
820
821 //
822 // The K stiffness matrix and R damping matrix of this element:
823 //
824
825 if (Kfactor || Rfactor) {
826 // Corotational K stiffness:
827 ChMatrixDynamic<> CK(12, 12);
828 ChMatrixDynamic<> CKCt(12, 12); // the global, corotated, K matrix
829
830 //
831 // Corotate local stiffness matrix
832 //
833
834 ChMatrix33<> Atoabs(this->q_element_abs_rot);
835 ChMatrix33<> AtolocwelA(this->GetNodeA()->Frame().GetRot().GetConjugate() % this->q_element_abs_rot);
836 ChMatrix33<> AtolocwelB(this->GetNodeB()->Frame().GetRot().GetConjugate() % this->q_element_abs_rot);
837 std::vector<ChMatrix33<>*> R;
838 R.push_back(&Atoabs);
839 R.push_back(&AtolocwelA);
840 R.push_back(&Atoabs);
841 R.push_back(&AtolocwelB);
842
843 if (this->use_geometric_stiffness) {
844 // compute Px tension of the beam along centerline, using temporary but fast data structures:
845 ChVectorDynamic<> displ(this->GetNdofs());
846 this->GetStateBlock(displ);
847 double Px = -this->Km.row(0) * displ;
848 // ChVector<> mFo, mTo;
849 // this->EvaluateSectionForceTorque(0, mFo, mTo); // for double checking the Px value
850 // GetLog() << " Px = " << Px << " Px_eval = " << mFo.x() << " \n";
851
852 // corotate Km + Kg (where Kg = this->Kg * Px)
853 ChMatrixCorotation::ComputeCK(this->Km + Px * this->Kg, R, 4, CK);
854 } else {
855 ChMatrixCorotation::ComputeCK(this->Km, R, 4, CK);
856 }
857 ChMatrixCorotation::ComputeKCt(CK, R, 4, CKCt);
858
859 // For strict symmetry, copy L=U because the computations above might
860 // lead to small errors because of numerical roundoff even with force_symmetric_stiffness
861 if (force_symmetric_stiffness) {
862 for (int row = 0; row < CKCt.rows() - 1; ++row)
863 for (int col = row + 1; col < CKCt.cols(); ++col)
864 CKCt(row, col) = CKCt(col, row);
865 }
866
867 //// RADU
868 //// Check if the above can be done with the one-liner:
869 ////CKCt.triangularView<Eigen::Upper>() = CKCt.transpose();
870
871 // The code below may be quite slow!!
872 ChMatrixDynamic<> CR(12, 12);
873 ChMatrixDynamic<> CRCt(12, 12); // the global, corotated, R matrix
874 ChMatrixCorotation::ComputeCK(this->Rm, R, 4, CR);
875 ChMatrixCorotation::ComputeKCt(CR, R, 4, CRCt);
876 if (force_symmetric_stiffness) {
877 for (int row = 0; row < CRCt.rows() - 1; ++row)
878 for (int col = row + 1; col < CRCt.cols(); ++col)
879 CRCt(row, col) = CRCt(col, row);
880 }
881
882 //// For K stiffness matrix and R matrix: scale by factors
883 // CKCt *= Kfactor + Rfactor * this->tapered_section->GetBeamRaleyghDamping();
884
885 H.block(0, 0, 12, 12) = CKCt * Kfactor + CRCt * Rfactor;
886
887 // Add inertial stiffness matrix and inertial damping matrix (gyroscopic damping),
888 // if enabled in this beam element.
889 // These matrices are not symmetric. Also note.
890 if (this->tapered_section->compute_inertia_damping_matrix ||
891 this->tapered_section->compute_inertia_stiffness_matrix) {
892 ChMatrixNM<double, 12, 12> matr_loc;
893 ChMatrixNM<double, 12, 12> KRi_loc;
894 KRi_loc.setZero();
895 // A lumped version of the inertial damping/stiffness matrix computation is used here, on a per-node basis:
896 double node_multiplier_fact_R = 0.5 * length * Rfactor;
897 double node_multiplier_fact_K = 0.5 * length * Kfactor;
898
899 ///< current angular velocity of section of node A, in material frame
900 ChVector<> mWvel_A = this->GetNodeA()->GetWvel_loc();
901 ///< current angular acceleration of section of node A, in material frame
902 ChVector<> mWacc_A = this->GetNodeA()->GetWacc_loc();
903 ///< current acceleration of section of node A, in material frame)
904 ChVector<> mXacc_A = this->GetNodeA()->TransformDirectionParentToLocal(this->GetNodeA()->GetPos_dtdt());
905 ///< current angular velocity of section of node B, in material frame
906 ChVector<> mWvel_B = this->GetNodeB()->GetWvel_loc();
907 ///< current angular acceleration of section of node B, in material frame
908 ChVector<> mWacc_B = this->GetNodeB()->GetWacc_loc();
909 ///< current acceleration of section of node B, in material frame
910 ChVector<> mXacc_B = this->GetNodeB()->TransformDirectionParentToLocal(this->GetNodeB()->GetPos_dtdt());
911
912 if (this->tapered_section->compute_inertia_damping_matrix) {
913 this->tapered_section->ComputeInertiaDampingMatrix(matr_loc, mWvel_A, mWacc_A);
914 KRi_loc += matr_loc * node_multiplier_fact_R;
915 // check the inertial damping matrix, it should be skew symmetric
916 // if ((matr_loc + matr_loc.transpose()).norm() < 1E-9) {
917 // std::cout << "Gyroscopic damping matrix is skew symmetric, that's correct!" << std::endl;
918 //} else {
919 // std::cout << "Gyroscopic damping matrix is NOT skew symmetric, something is wrong!" << std::endl;
920 //}
921 }
922
923 if (this->tapered_section->compute_inertia_stiffness_matrix) {
924 this->tapered_section->ComputeInertiaStiffnessMatrix(matr_loc, mWvel_A, mWacc_A, mXacc_A, mWvel_B,
925 mWacc_B, mXacc_B);
926 KRi_loc += matr_loc * node_multiplier_fact_K;
927 // check the inertial stiffness matrix, is it must be skew symmetric as the stiffness matrix?
928 // if ((matr_loc + matr_loc.transpose()).norm() < 1E-9) {
929 // std::cout << "Inertial stiffness matrix is skew symmetric" << std::endl;
930 //} else {
931 // std::cout << "Inertial stiffness matrix is NOT skew symmetric" << std::endl;
932 //}
933 }
934
935 for (int i = 0; i < nodes.size(); ++i) {
936 int stride = i * 6;
937 // corotate the local damping and stiffness matrices (at once, already scaled) into absolute one
938 // H.block<3, 3>(stride, stride ) += nodes[i]->GetA() * KRi_loc.block<3, 3>(0,0) *
939 // (nodes[i]->GetA().transpose()); // NOTE: not needed as KRi_loc.block<3, 3>(0,0) is null by
940 // construction
941 H.block<3, 3>(stride + 3, stride + 3) += KRi_loc.block<3, 3>(3, 3);
942 H.block<3, 3>(stride, stride + 3) += nodes[i]->GetA() * KRi_loc.block<3, 3>(0, 3);
943 // H.block<3, 3>(stride+3, stride) += KRi_loc.block<3, 3>(3,0) *
944 // (nodes[i]->GetA().transpose()); // NOTE: not needed as KRi_loc.block<3, 3>(3,0) is null by
945 // construction
946 }
947 }
948 } else
949 H.setZero();
950
951 //
952 // The M mass matrix of this element:
953 //
954
955 if (Mfactor) {
956 ChMatrixDynamic<> Mloc(12, 12);
957 Mloc.setZero();
958 ChMatrix33<> Mxw;
959
960 // ChMatrixNM<double, 12, 12> sectional_mass = this->M; // it could be consistent mass matrix, depends on
961 // SetLumpedMassMatrix(true/false)
962
963 if (this->tapered_section->GetLumpedMassMatrixType()) {
964 double node_multiplier_fact = 0.5 * this->length * Mfactor;
965 for (int i = 0; i < nodes.size(); ++i) {
966 int stride = i * 6;
967 // if there is no mass center offset, the upper right and lower left blocks need not be rotated,
968 // hence it can be the simple (constant) expression
969 // Mloc.block<6, 6>(stride, stride) += sectional_mass * node_multiplier_fact;
970 // but the more general case needs the rotations, hence:
971 Mloc.block<3, 3>(stride, stride) += this->M.block<3, 3>(stride, stride) * node_multiplier_fact;
972 Mloc.block<3, 3>(stride + 3, stride + 3) +=
973 this->M.block<3, 3>(stride + 3, stride + 3) * node_multiplier_fact;
974 Mxw = nodes[i]->GetA() * this->M.block<3, 3>(stride, stride + 3) * node_multiplier_fact;
975 Mloc.block<3, 3>(stride, stride + 3) += Mxw;
976 Mloc.block<3, 3>(stride + 3, stride) += Mxw.transpose();
977 }
978 // ..rather do this because lumped mass matrix does not need rotation transf.
979 H.block(0, 0, 12, 12) += Mloc;
980
981 } else {
982 Mloc = this->M * Mfactor;
983
984 // The following would be needed if consistent mass matrix is used, but...
985 ChMatrix33<> Atoabs(this->q_element_abs_rot);
986 ChMatrix33<> AtolocwelA(this->GetNodeA()->Frame().GetRot().GetConjugate() % this->q_element_abs_rot);
987 ChMatrix33<> AtolocwelB(this->GetNodeB()->Frame().GetRot().GetConjugate() % this->q_element_abs_rot);
988 std::vector<ChMatrix33<>*> R;
989 R.push_back(&Atoabs);
990 R.push_back(&AtolocwelA);
991 R.push_back(&Atoabs);
992 R.push_back(&AtolocwelB);
993
994 ChMatrixNM<double, 12, 12> CK;
995 ChMatrixNM<double, 12, 12> CKCt;
996 ChMatrixCorotation::ComputeCK(Mloc, R, 4, CK);
997 ChMatrixCorotation::ComputeKCt(CK, R, 4, CKCt);
998
999 H.block(0, 0, 12, 12) += CKCt;
1000 }
1001
1002 //// TODO better per-node lumping, or 4x4 consistent mass matrices, maybe with integration if not uniform
1003 // materials.
1004 }
1005 }
1006
ComputeInternalForces(ChVectorDynamic<> & Fi)1007 void ChElementBeamTaperedTimoshenko::ComputeInternalForces(ChVectorDynamic<>& Fi) {
1008 assert(Fi.size() == 12);
1009 assert(tapered_section);
1010
1011 // set up vector of nodal displacements and small rotations (in local element system)
1012 ChVectorDynamic<> displ(12);
1013 this->GetStateBlock(displ);
1014
1015 // [local Internal Forces] = [Klocal] * displ + [Rlocal] * displ_dt
1016 ChVectorDynamic<> Fi_local = this->Km * displ;
1017
1018 // set up vector of nodal velocities (in local element system)
1019 ChVectorDynamic<> displ_dt(12);
1020 this->GetField_dt(displ_dt);
1021
1022 // ChMatrixDynamic<> FiR_local = this->tapered_section->GetBeamRaleyghDamping() * this->Km * displ_dt;
1023 ChMatrixDynamic<> FiR_local = this->Rm * displ_dt;
1024
1025 Fi_local += FiR_local;
1026 Fi_local *= -1.0;
1027
1028 //
1029 // Corotate local internal loads
1030 //
1031
1032 // Fi = C * Fi_local with C block-diagonal rotations A , for nodal forces in abs. frame
1033 ChMatrix33<> Atoabs(this->q_element_abs_rot);
1034 ChMatrix33<> AtolocwelA(this->GetNodeA()->Frame().GetRot().GetConjugate() % this->q_element_abs_rot);
1035 ChMatrix33<> AtolocwelB(this->GetNodeB()->Frame().GetRot().GetConjugate() % this->q_element_abs_rot);
1036 std::vector<ChMatrix33<>*> R;
1037 R.push_back(&Atoabs);
1038 R.push_back(&AtolocwelA);
1039 R.push_back(&Atoabs);
1040 R.push_back(&AtolocwelB);
1041 ChMatrixCorotation::ComputeCK(Fi_local, R, 4, Fi);
1042
1043 // Add also inertial quadratic terms: gyroscopic and centrifugal
1044
1045 // CASE OF LUMPED MASS - fast
1046 double node_multiplier = 0.5 * length;
1047 ChVector<> mFcent_i;
1048 ChVector<> mTgyro_i;
1049 for (int i = 0; i < nodes.size(); ++i) {
1050 // int stride = i * 6;
1051 if (i == 0) {
1052 this->tapered_section->GetSectionA()->ComputeQuadraticTerms(mFcent_i, mTgyro_i, nodes[i]->GetWvel_loc());
1053 } else { // i==1
1054 this->tapered_section->GetSectionB()->ComputeQuadraticTerms(mFcent_i, mTgyro_i, nodes[i]->GetWvel_loc());
1055 }
1056 ChQuaternion<> q_i(nodes[i]->GetRot());
1057 Fi.segment(i * 6, 3) -= node_multiplier * (nodes[i]->GetA() * mFcent_i).eigen();
1058 Fi.segment(3 + i * 6, 3) -= node_multiplier * mTgyro_i.eigen();
1059 }
1060
1061 #ifdef BEAM_VERBOSE
1062 GetLog() << "\nInternal forces (local): \n";
1063 for (int c = 0; c < 6; c++)
1064 GetLog() << FiK_local(c) << " ";
1065 GetLog() << "\n";
1066 for (int c = 6; c < 12; c++)
1067 GetLog() << FiK_local(c) << " ";
1068 GetLog() << "\n\nInternal forces (ABS) : \n";
1069 for (int c = 0; c < 6; c++)
1070 GetLog() << Fi(c) << " ";
1071 GetLog() << "\n";
1072 for (int c = 6; c < 12; c++)
1073 GetLog() << Fi(c) << " ";
1074 GetLog() << "\n";
1075 #endif
1076 }
1077
ComputeInternalForces(ChVectorDynamic<> & Fi,bool Mfactor,bool Kfactor,bool Rfactor,bool Gfactor)1078 void ChElementBeamTaperedTimoshenko::ComputeInternalForces(ChVectorDynamic<>& Fi,
1079 bool Mfactor,
1080 bool Kfactor,
1081 bool Rfactor,
1082 bool Gfactor) {
1083 assert(Fi.size() == 12);
1084 assert(tapered_section);
1085
1086 ChVectorDynamic<> FiMKR_local(12);
1087 FiMKR_local.setZero();
1088
1089 if (Mfactor && true) {
1090 // set up vector of nodal accelerations (in local element system)
1091 ChVectorDynamic<> displ_dtdt(12);
1092 this->GetField_dtdt(displ_dtdt);
1093
1094 ChVectorDynamic<> FiM_local = this->M * displ_dtdt; // this->M is the local mass matrix of element
1095 FiMKR_local -= FiM_local;
1096 }
1097
1098 if (Kfactor) {
1099 // set up vector of nodal displacements and small rotations (in local element system)
1100 ChVectorDynamic<> displ(12);
1101 this->GetStateBlock(displ);
1102
1103 // [local Internal Forces] = [Klocal] * displ + [Rlocal] * displ_dt
1104 ChVectorDynamic<> FiK_local = this->Km * displ;
1105 FiMKR_local -= FiK_local;
1106 }
1107 if (Rfactor) {
1108 // set up vector of nodal velocities (in local element system)
1109 ChVectorDynamic<> displ_dt(12);
1110 this->GetField_dt(displ_dt);
1111
1112 // ChMatrixDynamic<> FiR_local = this->tapered_section->GetBeamRaleyghDamping() * this->Km * displ_dt;
1113 ChMatrixDynamic<> FiR_local = this->Rm * displ_dt;
1114 FiMKR_local -= FiR_local;
1115 }
1116
1117 //
1118 // Corotate local internal loads
1119 //
1120
1121 // Fi = C * Fi_local with C block-diagonal rotations A , for nodal forces in abs. frame
1122 ChMatrix33<> Atoabs(this->q_element_abs_rot);
1123 ChMatrix33<> AtolocwelA(this->GetNodeA()->Frame().GetRot().GetConjugate() % this->q_element_abs_rot);
1124 ChMatrix33<> AtolocwelB(this->GetNodeB()->Frame().GetRot().GetConjugate() % this->q_element_abs_rot);
1125 std::vector<ChMatrix33<>*> R;
1126 R.push_back(&Atoabs);
1127 R.push_back(&AtolocwelA);
1128 R.push_back(&Atoabs);
1129 R.push_back(&AtolocwelB);
1130 ChMatrixCorotation::ComputeCK(FiMKR_local, R, 4, Fi);
1131
1132 // DO NOT USE THIS PIECE OF CODE.
1133 // if (Mfactor && false) { // Is this correct? Or the code in lines 1029~1036
1134 // // set up vector of nodal accelarations (in absolute element system)
1135 // ChVectorDynamic<> displ_dtdt(12);
1136 // displ_dtdt.segment(0, 3) = nodes[0]->Frame().GetPos_dtdt().eigen();
1137 // // displ_dtdt.segment(3, 3) = nodes[0]->Frame().GetWacc_par().eigen(); //但是之后dynamics求解发散
1138 // displ_dtdt.segment(3, 3) = nodes[0]->Frame().GetWacc_loc().eigen(); //但是之后dynamics求解发散
1139
1140 // displ_dtdt.segment(6, 3) = nodes[1]->Frame().GetPos_dtdt().eigen();
1141 // // displ_dtdt.segment(9, 3) = nodes[1]->Frame().GetWacc_par().eigen();
1142 // displ_dtdt.segment(9, 3) = nodes[1]->Frame().GetWacc_loc().eigen();
1143 // // this->GetField_dtdt(displ_dtdt);
1144
1145 // ChMatrixDynamic<> Mabs(12, 12);
1146 // Mabs.setZero();
1147 // this->ComputeKRMmatricesGlobal(Mabs, false, false,
1148 // true); // Mabs is the mass matrix in absolute coordinate for this element
1149
1150 // ChVectorDynamic<> FiM_abs = Mabs * displ_dtdt;
1151 // Fi -= FiM_abs;
1152 //}
1153
1154 if (Gfactor) {
1155 // Add also inertial quadratic terms: gyroscopic and centrifugal
1156
1157 // CASE OF LUMPED MASS - fast
1158 double node_multiplier = 0.5 * length;
1159 ChVector<> mFcent_i;
1160 ChVector<> mTgyro_i;
1161 for (int i = 0; i < nodes.size(); ++i) {
1162 // int stride = i * 6;
1163 if (i == 0) {
1164 this->tapered_section->GetSectionA()->ComputeQuadraticTerms(mFcent_i, mTgyro_i,
1165 nodes[i]->GetWvel_loc());
1166 } else { // i==1
1167 this->tapered_section->GetSectionB()->ComputeQuadraticTerms(mFcent_i, mTgyro_i,
1168 nodes[i]->GetWvel_loc());
1169 }
1170 ChQuaternion<> q_i(nodes[i]->GetRot());
1171 Fi.segment(i * 6, 3) -= node_multiplier * (nodes[i]->GetA() * mFcent_i).eigen();
1172 Fi.segment(3 + i * 6, 3) -= node_multiplier * mTgyro_i.eigen();
1173 }
1174 }
1175
1176 #ifdef BEAM_VERBOSE
1177 GetLog() << "\nInternal forces (local): \n";
1178 for (int c = 0; c < 6; c++)
1179 GetLog() << FiK_local(c) << " ";
1180 GetLog() << "\n";
1181 for (int c = 6; c < 12; c++)
1182 GetLog() << FiK_local(c) << " ";
1183 GetLog() << "\n\nInternal forces (ABS) : \n";
1184 for (int c = 0; c < 6; c++)
1185 GetLog() << Fi(c) << " ";
1186 GetLog() << "\n";
1187 for (int c = 6; c < 12; c++)
1188 GetLog() << Fi(c) << " ";
1189 GetLog() << "\n";
1190 #endif
1191 }
1192
ComputeGravityForces(ChVectorDynamic<> & Fg,const ChVector<> & G_acc)1193 void ChElementBeamTaperedTimoshenko::ComputeGravityForces(ChVectorDynamic<>& Fg, const ChVector<>& G_acc) {
1194 // no so efficient... a temporary mass matrix here:
1195 ChMatrixDynamic<> mM(12, 12);
1196 this->ComputeMmatrixGlobal(mM);
1197
1198 // a vector of G accelerations for the two nodes (for translation degrees of freedom)
1199 ChVectorDynamic<> mG(12);
1200 mG.setZero();
1201 mG.segment(0, 3) = G_acc.eigen();
1202 mG.segment(6, 3) = G_acc.eigen();
1203
1204 // Gravity forces as M*g, always works, regardless of the way M
1205 // is computed (lumped or consistent, with offset center of mass or centered, etc.)
1206 // [Maybe one can replace this function with a faster ad-hoc implementation in case of lumped masses.]
1207 Fg = mM * mG;
1208
1209 //// TODO for the lumped mass matrix case, the mM * mG product can be unrolled into few multiplications as mM mostly
1210 /// zero, and same for mG
1211 }
1212
EvaluateSectionDisplacement(const double eta,ChVector<> & u_displ,ChVector<> & u_rotaz)1213 void ChElementBeamTaperedTimoshenko::EvaluateSectionDisplacement(const double eta,
1214 ChVector<>& u_displ,
1215 ChVector<>& u_rotaz) {
1216 ChVectorDynamic<> displ(this->GetNdofs());
1217 this->GetStateBlock(displ);
1218 // No transformation for the displacement of two nodes,
1219 // so the section displacement is evaluated at the centerline of beam
1220
1221 ShapeFunctionGroup NN;
1222 ShapeFunctionsTimoshenko(NN, eta);
1223 ShapeFunction5Blocks sfblk = std::get<1>(NN);
1224 ShapeFunction5Blocks sfblk1d = std::get<2>(NN);
1225 auto kNby = std::get<0>(sfblk);
1226 auto kNsy = std::get<1>(sfblk);
1227 auto kNbz = std::get<2>(sfblk);
1228 auto kNsz = std::get<3>(sfblk);
1229 auto kNx = std::get<4>(sfblk);
1230 auto dkNby = std::get<0>(sfblk1d);
1231 auto dkNbz = std::get<2>(sfblk1d);
1232
1233 ChVectorN<double, 4> qey;
1234 qey << displ(1), displ(5), displ(7), displ(11);
1235 ChVectorN<double, 4> qez;
1236 qez << displ(2), displ(4), displ(8), displ(10);
1237 ChVectorN<double, 2> qeux;
1238 qeux << displ(0), displ(6);
1239 ChVectorN<double, 2> qerx;
1240 qerx << displ(3), displ(9);
1241
1242 u_displ.x() = kNx * qeux;
1243 u_displ.y() = (kNby + kNsy) * qey;
1244 u_displ.z() = (kNbz + kNsz) * qez;
1245 u_rotaz.x() = kNx * qerx;
1246 u_rotaz.y() = -dkNbz * qez;
1247 u_rotaz.z() = dkNby * qey;
1248 }
1249
EvaluateSectionFrame(const double eta,ChVector<> & point,ChQuaternion<> & rot)1250 void ChElementBeamTaperedTimoshenko::EvaluateSectionFrame(const double eta, ChVector<>& point, ChQuaternion<>& rot) {
1251 ChVector<> u_displ;
1252 ChVector<> u_rotaz;
1253 double Nx1 = (1. / 2.) * (1 - eta);
1254 double Nx2 = (1. / 2.) * (1 + eta);
1255
1256 this->EvaluateSectionDisplacement(eta, u_displ, u_rotaz);
1257
1258 // Since d = [Atw]' Xt - [A0w]'X0 , so we have
1259 // Xt = [Atw] (d + [A0w]'X0)
1260
1261 point = this->q_element_abs_rot.Rotate(u_displ +
1262 this->q_element_ref_rot.RotateBack(Nx1 * this->nodes[0]->GetX0().GetPos() +
1263 Nx2 * this->nodes[1]->GetX0().GetPos()));
1264
1265 ChQuaternion<> msectionrot;
1266 msectionrot.Q_from_AngAxis(u_rotaz.Length(), u_rotaz.GetNormalized());
1267 rot = this->q_element_abs_rot % msectionrot;
1268 }
1269
EvaluateSectionForceTorque(const double eta,ChVector<> & Fforce,ChVector<> & Mtorque)1270 void ChElementBeamTaperedTimoshenko::EvaluateSectionForceTorque(const double eta,
1271 ChVector<>& Fforce,
1272 ChVector<>& Mtorque) {
1273 assert(tapered_section);
1274
1275 ChVectorDynamic<> displ(this->GetNdofs());
1276 this->GetStateBlock(displ);
1277
1278 ChVectorDynamic<> displ_ec = this->T * displ; // transform the displacement of two nodes to elastic axis
1279
1280 ChVectorN<double, 4> qey;
1281 qey << displ_ec(1), displ_ec(5), displ_ec(7), displ_ec(11);
1282 ChVectorN<double, 4> qez;
1283 qez << displ_ec(2), displ_ec(4), displ_ec(8), displ_ec(10);
1284 ChVectorN<double, 2> qeux;
1285 qeux << displ_ec(0), displ_ec(6);
1286 ChVectorN<double, 2> qerx;
1287 qerx << displ_ec(3), displ_ec(9);
1288
1289 ShapeFunctionGroup NN;
1290 ShapeFunctionsTimoshenko(NN, eta);
1291 ShapeFunction5Blocks sfblk1d = std::get<2>(NN);
1292 ShapeFunction2Blocks sfblk2d = std::get<3>(NN);
1293 ShapeFunction2Blocks sfblk3d = std::get<4>(NN);
1294 auto dkNx1 = std::get<4>(sfblk1d);
1295 auto dkNsy = std::get<1>(sfblk1d);
1296 auto dkNsz = std::get<3>(sfblk1d);
1297 auto ddkNby = std::get<0>(sfblk2d);
1298 auto ddkNbz = std::get<1>(sfblk2d);
1299 auto dddkNby = std::get<0>(sfblk3d);
1300 auto dddkNbz = std::get<1>(sfblk3d);
1301
1302 double EA = this->tapered_section->GetAverageSectionParameters()->EA;
1303 double GJ = this->tapered_section->GetAverageSectionParameters()->GJ;
1304 double GAyy = this->tapered_section->GetAverageSectionParameters()->GAyy;
1305 double GAzz = this->tapered_section->GetAverageSectionParameters()->GAzz;
1306 double EIyy = this->tapered_section->GetAverageSectionParameters()->EIyy;
1307 double EIzz = this->tapered_section->GetAverageSectionParameters()->EIzz;
1308
1309 double eps = 1.0e-3;
1310 bool use_shear_stain = true; // As default, use shear strain to evaluate the shear forces
1311 if (abs(GAyy) < eps ||
1312 abs(GAzz) < eps) { // Sometimes, the user will not input GAyy, GAzz, so GAyy, GAzz may be zero.
1313 use_shear_stain = false;
1314 }
1315
1316 // generalized strains/curvatures;
1317 ChVectorN<double, 6> sect_ek;
1318 sect_ek(0) = dkNx1 * qeux; // ux
1319 sect_ek(3) = dkNx1 * qerx; // rotx
1320 sect_ek(4) = -ddkNbz * qez; // roty
1321 sect_ek(5) = ddkNby * qey; // rotz
1322
1323 if (use_shear_stain) {
1324 // Strictly speaking, dkNsy * qey,dkNsz * qez are the real shear strains.
1325 sect_ek(1) = dkNsy * qey; // gamma_y = dkNsy * qey; Fy = GAyy * gamma_y;
1326 sect_ek(2) = dkNsz * qez; // gamma_z = dkNsz * qez; Fz = GAzz * gamma_z;
1327 } else {
1328 // Calculate the shear strain through third-differentian of bending displacement
1329 sect_ek(1) = -dddkNby * qey; // Fy == -EIzz * dddkNby * qey == GAyy * dkNsy * qey;
1330 sect_ek(2) = -dddkNbz * qez; // Fz == -EIyy * dddkNbz * qez == GAzz * dkNsz * qez;
1331 }
1332
1333 if (false) // section->alpha ==0 && section->Cy ==0 && section->Cz==0 && section->Sy==0 && section->Sz==0)
1334 {
1335 // Fast computation:
1336 Fforce.x() = EA * sect_ek(0);
1337 if (use_shear_stain) {
1338 Fforce.y() = GAyy * sect_ek(1);
1339 Fforce.z() = GAzz * sect_ek(2);
1340 } else {
1341 Fforce.y() = EIzz * sect_ek(1);
1342 Fforce.z() = EIyy * sect_ek(2);
1343 }
1344
1345 Mtorque.x() = GJ * sect_ek(3);
1346 Mtorque.y() = EIyy * sect_ek(4);
1347 Mtorque.z() = EIzz * sect_ek(5);
1348 } else {
1349 // Constitutive matrix of the beam:
1350 ChMatrixNM<double, 6, 6> Klaw_d;
1351 Klaw_d.setZero();
1352 Klaw_d(0, 0) = EA;
1353 Klaw_d(3, 3) = GJ;
1354 Klaw_d(4, 4) = EIyy;
1355 Klaw_d(5, 5) = EIzz;
1356
1357 if (use_shear_stain) {
1358 Klaw_d(1, 1) = GAyy;
1359 Klaw_d(2, 2) = GAzz;
1360 } else {
1361 Klaw_d(1, 1) = EIzz;
1362 Klaw_d(2, 2) = EIyy;
1363 }
1364
1365 ChMatrixDynamic<> Teta;
1366 ComputeTransformMatrixAtPoint(Teta, eta);
1367
1368 // ..unrolled rotated constitutive matrix..
1369 ChMatrixNM<double, 6, 6> Klaw_r;
1370 Klaw_r.setZero();
1371 Klaw_r = Teta.transpose() * Klaw_d;
1372
1373 // .. compute wrench = Klaw_r * sect_ek
1374 ChVectorN<double, 6> wrench = Klaw_r * sect_ek;
1375 Fforce = wrench.segment(0, 3);
1376 Mtorque = wrench.segment(3, 3);
1377 }
1378 }
1379
EvaluateSectionStrain(const double eta,ChVector<> & StrainV_trans,ChVector<> & StrainV_rot)1380 void ChElementBeamTaperedTimoshenko::EvaluateSectionStrain(const double eta,
1381 ChVector<>& StrainV_trans,
1382 ChVector<>& StrainV_rot) {
1383 assert(tapered_section);
1384
1385 ChVectorDynamic<> displ(this->GetNdofs());
1386 this->GetStateBlock(displ);
1387
1388 ChVectorDynamic<> displ_ec = this->T * displ; // transform the displacement of two nodes to elastic axis
1389
1390 ChVectorN<double, 4> qey;
1391 qey << displ_ec(1), displ_ec(5), displ_ec(7), displ_ec(11);
1392 ChVectorN<double, 4> qez;
1393 qez << displ_ec(2), displ_ec(4), displ_ec(8), displ_ec(10);
1394 ChVectorN<double, 2> qeux;
1395 qeux << displ_ec(0), displ_ec(6);
1396 ChVectorN<double, 2> qerx;
1397 qerx << displ_ec(3), displ_ec(9);
1398
1399 ShapeFunctionGroup NN;
1400 ShapeFunctionsTimoshenko(NN, eta);
1401 ShapeFunction5Blocks sfblk1d = std::get<2>(NN);
1402 ShapeFunction2Blocks sfblk2d = std::get<3>(NN);
1403 ShapeFunction2Blocks sfblk3d = std::get<4>(NN);
1404 auto dkNx1 = std::get<4>(sfblk1d);
1405 auto dkNsy = std::get<1>(sfblk1d);
1406 auto dkNsz = std::get<3>(sfblk1d);
1407 auto ddkNby = std::get<0>(sfblk2d);
1408 auto ddkNbz = std::get<1>(sfblk2d);
1409 auto dddkNby = std::get<0>(sfblk3d);
1410 auto dddkNbz = std::get<1>(sfblk3d);
1411
1412 ////double EA = this->tapered_section->GetAverageSectionParameters()->EA;
1413 ////double GJ = this->tapered_section->GetAverageSectionParameters()->GJ;
1414 double GAyy = this->tapered_section->GetAverageSectionParameters()->GAyy;
1415 double GAzz = this->tapered_section->GetAverageSectionParameters()->GAzz;
1416 ////double EIyy = this->tapered_section->GetAverageSectionParameters()->EIyy;
1417 ////double EIzz = this->tapered_section->GetAverageSectionParameters()->EIzz;
1418
1419 double eps = 1.0e-3;
1420 bool use_shear_stain = true; // As default, use shear strain to evaluate the shear forces
1421 if (abs(GAyy) < eps ||
1422 abs(GAzz) < eps) { // Sometimes, the user will not input GAyy, GAzz, so GAyy, GAzz may be zero.
1423 use_shear_stain = false;
1424 }
1425
1426 // generalized strains/curvatures;
1427 ChVectorN<double, 6> sect_ek;
1428 sect_ek(0) = dkNx1 * qeux; // ux
1429 sect_ek(3) = dkNx1 * qerx; // rotx
1430 sect_ek(4) = -ddkNbz * qez; // roty
1431 sect_ek(5) = ddkNby * qey; // rotz
1432
1433 if (use_shear_stain) {
1434 // Strictly speaking, dkNsy * qey,dkNsz * qez are the real shear strains.
1435 sect_ek(1) = dkNsy * qey; // gamma_y = dkNsy * qey; Fy = GAyy * gamma_y;
1436 sect_ek(2) = dkNsz * qez; // gamma_z = dkNsz * qez; Fz = GAzz * gamma_z;
1437 } else {
1438 // Calculate the shear strain through third-differentian of bending displacement
1439 sect_ek(1) = -dddkNby * qey; // Fy == -EIzz * dddkNby * qey == GAyy * dkNsy * qey;
1440 sect_ek(2) = -dddkNbz * qez; // Fz == -EIyy * dddkNbz * qez == GAzz * dkNsz * qez;
1441 }
1442
1443 StrainV_trans = sect_ek.segment(0, 3);
1444 StrainV_rot = sect_ek.segment(3, 3);
1445 }
1446
EvaluateElementStrainEnergy(ChVector<> & StrainEnergyV_trans,ChVector<> & StrainEnergyV_rot)1447 void ChElementBeamTaperedTimoshenko::EvaluateElementStrainEnergy(ChVector<>& StrainEnergyV_trans,
1448 ChVector<>& StrainEnergyV_rot) {
1449 ChVectorDynamic<> displ(this->GetNdofs());
1450 this->GetStateBlock(displ);
1451
1452 ChVectorN<double, 12> strain_energy = 1.0 / 2.0 * displ.asDiagonal() * this->Km * displ;
1453 ChVectorN<double, 6> strain_energy_v;
1454 // double strain_energy_sum = 0;
1455 for (int i = 0; i < 6; i++) {
1456 strain_energy_v(i) = strain_energy(i) + strain_energy(i + 6);
1457 // strain_energy_sum += strain_energy_v(i);
1458 }
1459
1460 StrainEnergyV_trans = strain_energy_v.segment(0, 3);
1461 StrainEnergyV_rot = strain_energy_v.segment(3, 3);
1462 }
1463
EvaluateElementDampingEnergy(ChVector<> & DampingEnergyV_trans,ChVector<> & DampingEnergyV_rot)1464 void ChElementBeamTaperedTimoshenko::EvaluateElementDampingEnergy(ChVector<>& DampingEnergyV_trans,
1465 ChVector<>& DampingEnergyV_rot) {
1466 ChVectorDynamic<> displ_dt(this->GetNdofs());
1467 this->GetField_dt(displ_dt);
1468
1469 ChVectorN<double, 12> damping_energy = 1.0 / 2.0 * displ_dt.asDiagonal() * this->Rm * displ_dt;
1470 ChVectorN<double, 6> damping_energy_v;
1471 // double damping_energy_sum = 0;
1472 for (int i = 0; i < 6; i++) {
1473 damping_energy_v(i) = damping_energy(i) + damping_energy(i + 6);
1474 // damping_energy_sum += damping_energy_v(i);
1475 }
1476
1477 DampingEnergyV_trans = damping_energy_v.segment(0, 3);
1478 DampingEnergyV_rot = damping_energy_v.segment(3, 3);
1479 }
1480
LoadableGetStateBlock_x(int block_offset,ChState & mD)1481 void ChElementBeamTaperedTimoshenko::LoadableGetStateBlock_x(int block_offset, ChState& mD) {
1482 mD.segment(block_offset + 0, 3) = nodes[0]->GetPos().eigen();
1483 mD.segment(block_offset + 3, 4) = nodes[0]->GetRot().eigen();
1484
1485 mD.segment(block_offset + 7, 3) = nodes[1]->GetPos().eigen();
1486 mD.segment(block_offset + 10, 4) = nodes[1]->GetRot().eigen();
1487 }
1488
LoadableGetStateBlock_w(int block_offset,ChStateDelta & mD)1489 void ChElementBeamTaperedTimoshenko::LoadableGetStateBlock_w(int block_offset, ChStateDelta& mD) {
1490 mD.segment(block_offset + 0, 3) = nodes[0]->GetPos_dt().eigen();
1491 mD.segment(block_offset + 3, 3) = nodes[0]->GetWvel_loc().eigen();
1492
1493 mD.segment(block_offset + 6, 3) = nodes[1]->GetPos_dt().eigen();
1494 mD.segment(block_offset + 9, 3) = nodes[1]->GetWvel_loc().eigen();
1495 }
1496
LoadableStateIncrement(const unsigned int off_x,ChState & x_new,const ChState & x,const unsigned int off_v,const ChStateDelta & Dv)1497 void ChElementBeamTaperedTimoshenko::LoadableStateIncrement(const unsigned int off_x,
1498 ChState& x_new,
1499 const ChState& x,
1500 const unsigned int off_v,
1501 const ChStateDelta& Dv) {
1502 nodes[0]->NodeIntStateIncrement(off_x, x_new, x, off_v, Dv);
1503 nodes[1]->NodeIntStateIncrement(off_x + 7, x_new, x, off_v + 6, Dv);
1504 }
1505
LoadableGetVariables(std::vector<ChVariables * > & mvars)1506 void ChElementBeamTaperedTimoshenko::LoadableGetVariables(std::vector<ChVariables*>& mvars) {
1507 mvars.push_back(&this->nodes[0]->Variables());
1508 mvars.push_back(&this->nodes[1]->Variables());
1509 }
1510
ComputeNF(const double U,ChVectorDynamic<> & Qi,double & detJ,const ChVectorDynamic<> & F,ChVectorDynamic<> * state_x,ChVectorDynamic<> * state_w)1511 void ChElementBeamTaperedTimoshenko::ComputeNF(const double U,
1512 ChVectorDynamic<>& Qi,
1513 double& detJ,
1514 const ChVectorDynamic<>& F,
1515 ChVectorDynamic<>* state_x,
1516 ChVectorDynamic<>* state_w) {
1517 ShapeFunctionGroup NN;
1518 ShapeFunctionsTimoshenko(NN, U);
1519 ShapeFunctionN N = std::get<0>(NN);
1520
1521 // eta = 2*x/L;
1522 // ---> Deta/dx = 2./L;
1523 // ---> detJ = dx/Deta = L/2.;
1524 detJ = this->GetRestLength() / 2.0;
1525
1526 Qi = N.transpose() * F;
1527 }
1528
ComputeNF(const double U,const double V,const double W,ChVectorDynamic<> & Qi,double & detJ,const ChVectorDynamic<> & F,ChVectorDynamic<> * state_x,ChVectorDynamic<> * state_w)1529 void ChElementBeamTaperedTimoshenko::ComputeNF(const double U,
1530 const double V,
1531 const double W,
1532 ChVectorDynamic<>& Qi,
1533 double& detJ,
1534 const ChVectorDynamic<>& F,
1535 ChVectorDynamic<>* state_x,
1536 ChVectorDynamic<>* state_w) {
1537 this->ComputeNF(U, Qi, detJ, F, state_x, state_w);
1538 detJ /= 4.0; // because volume
1539 }
1540
GetDensity()1541 double ChElementBeamTaperedTimoshenko::GetDensity() {
1542 double mu1 = this->tapered_section->GetSectionA()->GetMassPerUnitLength();
1543 double mu2 = this->tapered_section->GetSectionB()->GetMassPerUnitLength();
1544 double mu = (mu1 + mu2) / 2.0;
1545
1546 return mu;
1547 }
1548
1549 } // end namespace fea
1550 } // end namespace chrono
1551