1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/aeromodal.cc,v 1.64 2017/01/12 14:45:58 masarati Exp $ */
2 /*
3 * MBDyn (C) is a multibody analysis code.
4 * http://www.mbdyn.org
5 *
6 * Copyright (C) 1996-2017
7 *
8 * Pierangelo Masarati <masarati@aero.polimi.it>
9 * Paolo Mantegazza <mantegazza@aero.polimi.it>
10 *
11 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12 * via La Masa, 34 - 20156 Milano, Italy
13 * http://www.aero.polimi.it
14 *
15 * Changing this copyright notice is forbidden.
16 *
17 * This program is free software; you can redistribute it and/or modify
18 * it under the terms of the GNU General Public License as published by
19 * the Free Software Foundation (version 2 of the License).
20 *
21 *
22 * This program is distributed in the hope that it will be useful,
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 * GNU General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License
28 * along with this program; if not, write to the Free Software
29 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30 */
31 /* Elemento aerodinamico modale */
32
33 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
34
35 #include <cerrno>
36
37 #include "aeromodal.h"
38 #include "dataman.h"
39 #include "Rot.hh"
40
41 /* AerodynamicModal - begin */
42
AerodynamicModal(unsigned int uLabel,const StructNode * pN,const Modal * pMJ,const Mat3x3 & RaTmp,const DofOwner * pDO,doublereal Cd,const int NModal,const int NAero,RigidF_t rgF,const int Gust,const doublereal Vff,SpMapMatrixHandler * pAMat,FullMatrixHandler * pBMat,FullMatrixHandler * pCMat,FullMatrixHandler * pD0Mat,FullMatrixHandler * pD1Mat,FullMatrixHandler * pD2Mat,flag fout)43 AerodynamicModal::AerodynamicModal(unsigned int uLabel,
44 const StructNode* pN,
45 const Modal* pMJ,
46 const Mat3x3& RaTmp,
47 const DofOwner* pDO,
48 doublereal Cd,
49 const int NModal,
50 const int NAero,
51 RigidF_t rgF,
52 const int Gust,
53 const doublereal Vff,
54 SpMapMatrixHandler* pAMat,
55 FullMatrixHandler* pBMat,
56 FullMatrixHandler* pCMat,
57 FullMatrixHandler* pD0Mat,
58 FullMatrixHandler* pD1Mat,
59 FullMatrixHandler* pD2Mat,
60 flag fout)
61 : Elem(uLabel, fout),
62 AerodynamicElem(uLabel, pDO, fout),
63 InitialAssemblyElem(uLabel, fout),
64 pModalNode(pN), pModalJoint(pMJ),
65 Ra(RaTmp),
66 Chord(Cd),
67 NStModes(NModal),
68 NAeroStates(NAero),
69 NGust(Gust),
70 pA(pAMat), pB(pBMat), pC(pCMat),
71 pD0(pD0Mat), pD1(pD1Mat), pD2(pD2Mat),
72 pq(0), pqPrime(0), pqSec(0),
73 pxa(0), pxaPrime(0),
74 pgs(0), pgsPrime(0),
75 gustVff(Vff), gustXi(0.707),
76 RigidF(rgF)
77 {
78 DEBUGCOUTFNAME("AerodynamicModal::AerodynamicModal");
79
80 R0 = pModalNode->GetRCurr()*Ra;
81 P0 = R0.MulTV(pModalNode->GetXCurr());
82
83 SAFENEWWITHCONSTRUCTOR(pq, MyVectorHandler,
84 MyVectorHandler(NStModes+RigidF));
85 SAFENEWWITHCONSTRUCTOR(pqPrime, MyVectorHandler,
86 MyVectorHandler(NStModes+RigidF));
87 SAFENEWWITHCONSTRUCTOR(pqSec, MyVectorHandler,
88 MyVectorHandler(NStModes+RigidF));
89 SAFENEWWITHCONSTRUCTOR(pxa, MyVectorHandler,
90 MyVectorHandler(NAeroStates));
91 SAFENEWWITHCONSTRUCTOR(pxaPrime, MyVectorHandler,
92 MyVectorHandler(NAeroStates));
93 SAFENEWWITHCONSTRUCTOR(pgs, MyVectorHandler,
94 MyVectorHandler(2*NGust));
95 SAFENEWWITHCONSTRUCTOR(pgsPrime, MyVectorHandler,
96 MyVectorHandler(2*NGust));
97
98 ASSERT(pModalNode != 0);
99 ASSERT(pModalNode->GetNodeType() == Node::STRUCTURAL);
100 }
101
~AerodynamicModal(void)102 AerodynamicModal::~AerodynamicModal(void)
103 {
104 DEBUGCOUTFNAME("AerodynamicModal::~AerodynamicModal");
105 if (pq != 0) {
106 SAFEDELETE(pq);
107 }
108 if (pqPrime != 0) {
109 SAFEDELETE(pqPrime);
110 }
111 if (pqSec != 0) {
112 SAFEDELETE(pqSec);
113 }
114 if (pxa != 0) {
115 SAFEDELETE(pxa);
116 }
117 if (pxaPrime != 0) {
118 SAFEDELETE(pxaPrime);
119 }
120 if (pA != 0) {
121 SAFEDELETE(pA);
122 }
123 if (pB != 0) {
124 SAFEDELETE(pB);
125 }
126 if (pC != 0) {
127 SAFEDELETE(pC);
128 }
129 if (pD0 != 0) {
130 SAFEDELETE(pD0);
131 }
132 if (pD1 != 0) {
133 SAFEDELETE(pD1);
134 }
135 if (pD2 != 0) {
136 SAFEDELETE(pD2);
137 }
138 }
139
140 /* Scrive il contributo dell'elemento al file di restart */
141 std::ostream&
Restart(std::ostream & out) const142 AerodynamicModal::Restart(std::ostream& out) const
143 {
144 return out << " /* aerodynamic modal: not implemented yet */" << std::endl;
145 }
146
147 /* assemblaggio jacobiano */
148 VariableSubMatrixHandler&
AssJac(VariableSubMatrixHandler & WorkMat,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)149 AerodynamicModal::AssJac(VariableSubMatrixHandler& WorkMat,
150 doublereal dCoef,
151 const VectorHandler& XCurr,
152 const VectorHandler& XPrimeCurr )
153 {
154 DEBUGCOUT("Entering AerodynamicModal::AssJac()" << std::endl);
155
156 FullSubMatrixHandler& WM = WorkMat.SetFull();
157
158 integer iNumRows = 0;
159 integer iNumCols = 0;
160 WorkSpaceDim(&iNumRows, &iNumCols);
161 WM.ResizeReset(iNumRows, iNumCols);
162
163 integer iModalIndex = pModalJoint->iGetModalIndex();
164 for (unsigned int iCnt = 1; iCnt <= NStModes; iCnt++) {
165 WM.PutRowIndex(iCnt, iModalIndex+NStModes+iCnt);
166 WM.PutColIndex(iCnt, iModalIndex+iCnt);
167 WM.PutColIndex(NStModes+iCnt, iModalIndex+NStModes+iCnt);
168 }
169
170 integer iAeroIndex = iGetFirstIndex();
171 for (unsigned int iCnt = 1; iCnt <= NAeroStates+2*NGust; iCnt++) {
172 WM.PutRowIndex(NStModes+iCnt, iAeroIndex+iCnt);
173 WM.PutColIndex(2*NStModes+iCnt, iAeroIndex+iCnt);
174 }
175
176 /* Dati del nodo rigido */
177 const Vec3& X0(pModalNode->GetXCurr());
178 const Vec3& V0(pModalNode->GetVCurr());
179 const Mat3x3& Rn(pModalNode->GetRCurr());
180 Mat3x3 RR(Rn*Ra);
181 Vec3 Vr(V0);
182
183 doublereal rho, vs, p, T;
184 GetAirProps(X0, rho, vs, p, T); /* p, T are not used yet */
185
186 Vec3 VTmp(Zero3);
187 if (fGetAirVelocity(VTmp, X0)) {
188 Vr -= VTmp;
189 }
190
191 /* velocit� nel riferimento nodale aerodinamico */
192 VTmp = RR.MulTV(Vr);
193 doublereal nV = std::abs(VTmp(1));
194 /* doublereal CV=Chord/(2*nV); */
195 doublereal qd = 0.5*rho*nV*nV;
196 doublereal qd1 = 0.25*rho*nV*Chord; /* qd*CV */
197 doublereal qd2 = 0.125*rho*Chord*Chord; /* qd*CV*CV */
198
199 /* parte deformabile :
200 *
201 * | ||aP|
202 * | || |
203 * | cKae Mae+cCae -cqC ||bP|
204 * | || |
205 * |-c(1/CV)B 0 I-c(1/CV)A ||xa|
206 *
207 * con
208 * Mae=-qd CV^2 D2
209 * Cae=-qd CV D1
210 * Kae=-qd D0
211 */
212 for (unsigned int i = 1; i <= NStModes; i++) {
213 for (unsigned int j = 1; j <= NStModes; j++) {
214 WM.IncCoef(i, j,
215 -dCoef*qd*pD0->operator()(RigidF + i, RigidF + j));
216 WM.IncCoef(i, j + NStModes,
217 -qd2*pD2->operator()(RigidF + i, RigidF + j)
218 -dCoef*qd1*pD1->operator()(RigidF + i, RigidF + j));
219 }
220 }
221
222 for (unsigned int j = 1; j <= NAeroStates; j++) {
223 for (unsigned int i = 1; i <= NStModes; i++) {
224 WM.IncCoef(i, j + 2*NStModes,
225 -dCoef*qd*pC->operator()(RigidF + i, j));
226 WM.IncCoef(j + NStModes, i,
227 -dCoef*(2*nV/Chord)*pB->operator()(j, RigidF + i));
228 }
229 }
230
231 for (unsigned int j = 1; j <= NAeroStates; j++) {
232 for (unsigned int i = 1; i <= NAeroStates; i++) {
233 WM.IncCoef(i + NStModes, j + 2*NStModes,
234 1.*(i == j) - dCoef*(2*nV/Chord)*pA->operator()(i, j));
235 }
236 }
237
238 if (NGust) {
239 for (unsigned int i = 1; i <= NStModes; i++) {
240 WM.IncCoef(i, 2*NStModes + NAeroStates + 1,
241 -dCoef*qd*pD0->operator()(RigidF + i, RigidF + NStModes + 1));
242 WM.IncCoef(i, 2*NStModes + NAeroStates + 3,
243 -dCoef*qd*pD0->operator()(RigidF + i, RigidF + NStModes + 2));
244 WM.IncCoef(i, 2*NStModes + NAeroStates + 2,
245 -qd2*pD2->operator()(RigidF + i, RigidF + NStModes + 1)
246 -dCoef*qd1*pD1->operator()(RigidF + i, RigidF + NStModes + 1));
247 WM.IncCoef(i, 2*NStModes + NAeroStates + 4,
248 -qd2*pD2->operator()(RigidF + i, RigidF + NStModes + 2)
249 -dCoef*qd1*pD1->operator()(RigidF + i, RigidF + NStModes + 2));
250 }
251
252 for (unsigned int i = 1; i <= NAeroStates; i++) {
253 WM.IncCoef(i + NStModes, 2*NStModes + NAeroStates + 1,
254 -dCoef*(2*nV/Chord)*pB->operator()(i, RigidF + NStModes + 1));
255 WM.IncCoef(i + NStModes, 2*NStModes + NAeroStates + 3,
256 -dCoef*(2*nV/Chord)*pB->operator()(i, RigidF + NStModes + 2));
257 }
258
259 for (unsigned int i = 0; i < NGust; i++) {
260 WM.IncCoef(i*2 + 1 + NStModes+NAeroStates,
261 i*2 + 1 + 2*NStModes + NAeroStates, 1.);
262 WM.IncCoef(i*2 + 1 + NStModes+NAeroStates,
263 i*2 + 2 + 2*NStModes + NAeroStates, -dCoef);
264 WM.IncCoef(i*2 + 2 + NStModes+NAeroStates,
265 i*2 + 1 + 2*NStModes+NAeroStates,
266 dCoef*gustVff*gustVff);
267 WM.IncCoef(i*2 + 2 + NStModes + NAeroStates,
268 i*2 + 2 + 2*NStModes + NAeroStates,
269 1. + dCoef*2*gustXi*gustVff);
270 }
271 }
272
273 return WorkMat;
274 }
275
276
277 SubVectorHandler&
AssRes(SubVectorHandler & WorkVec,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)278 AerodynamicModal::AssRes(SubVectorHandler& WorkVec,
279 doublereal dCoef,
280 const VectorHandler& XCurr,
281 const VectorHandler& XPrimeCurr)
282 {
283 DEBUGCOUTFNAME("AerodynamicModal::AssRes");
284 WorkVec.ResizeReset(RigidF + NStModes + NAeroStates + 2*NGust);
285
286 if (RigidF) {
287 integer iFirstIndex = pModalNode->iGetFirstMomentumIndex();
288 for (int iCnt = 1; iCnt <= RigidF; iCnt++) {
289 WorkVec.PutRowIndex(iCnt, iFirstIndex + iCnt);
290 }
291
292 const Vec3& X0(pModalNode->GetXCurr());
293 const Mat3x3& Rn(pModalNode->GetRCurr());
294
295 Mat3x3 RR(Rn*Ra);
296
297 // q
298 pq->Put(1, RR.MulTV(X0) - P0);
299 pq->Put(4, RotManip::VecRot(RR.MulMT(R0))); // nota: wrappa; se serve si pu� eliminare
300
301 // dot{q}
302 const Vec3& V0(pModalNode->GetVCurr());
303 const Vec3& W0(pModalNode->GetWCurr());
304 pqPrime->Put(1, RR.MulTV(V0 + X0.Cross(W0)));
305 pqPrime->Put(4, RR.MulTV(W0));
306
307 // ddot{q}
308 const Vec3& XPP0(pModalNode->GetXPPCurr());
309 const Vec3& WP0(pModalNode->GetWPCurr());
310 pqSec->Put(1, RR.MulTV(XPP0)); // verificare: non mancano i termini di trasporto ecc?
311 pqSec->Put(4, RR.MulTV(WP0));
312 }
313
314 integer iModalIndex = pModalJoint->iGetModalIndex();
315
316 for (unsigned int iCnt = 1; iCnt <= NStModes; iCnt++) {
317 WorkVec.PutRowIndex(RigidF + iCnt, iModalIndex + NStModes + iCnt);
318 }
319
320 integer iAeroIndex = iGetFirstIndex();
321 for (unsigned int iCnt = 1; iCnt <= NAeroStates + 2*NGust; iCnt++) {
322 WorkVec.PutRowIndex(RigidF + NStModes+iCnt, iAeroIndex + iCnt);
323 }
324
325 /* Recupera i vettori {a} e {aP} e {aS} (deformate modali) */
326 // FIXME: get this info from modal joint?
327 for (unsigned int iCnt = 1; iCnt <= NStModes; iCnt++) {
328 pq->PutCoef(iCnt + RigidF, XCurr(iModalIndex + iCnt));
329 pqPrime->PutCoef(iCnt + RigidF, XPrimeCurr(iModalIndex + iCnt));
330 pqSec->PutCoef(iCnt + RigidF, XPrimeCurr(iModalIndex + NStModes + iCnt));
331 }
332
333 /* Recupera i vettori {xa} e {xaP} */
334 for (unsigned int iCnt = 1; iCnt <= NAeroStates; iCnt++) {
335 pxa->PutCoef(iCnt, XCurr(iAeroIndex + iCnt));
336 pxaPrime->PutCoef(iCnt, XPrimeCurr(iAeroIndex + iCnt));
337 }
338
339 for (unsigned int iCnt = 1; iCnt <= 2*NGust; iCnt++) {
340 pgs->PutCoef(iCnt, XCurr(iAeroIndex + NAeroStates + iCnt));
341 pgsPrime->PutCoef(iCnt, XPrimeCurr(iAeroIndex + NAeroStates + iCnt));
342 }
343
344 AssVec(WorkVec);
345
346 return WorkVec;
347 }
348
349 SubVectorHandler&
InitialAssRes(SubVectorHandler & WorkVec,const VectorHandler & XCurr)350 AerodynamicModal::InitialAssRes(SubVectorHandler& WorkVec,
351 const VectorHandler& XCurr)
352 {
353 DEBUGCOUTFNAME("AerodynamicModal::AssRes");
354 const Vec3& X0(pModalNode->GetXCurr());
355 const Mat3x3& Rn(pModalNode->GetRCurr());
356 R0 = Rn*Ra;
357 P0 = R0.MulTV(X0);
358
359 WorkVec.ResizeReset(RigidF + NStModes + NAeroStates + 2*NGust);
360
361 if (RigidF) {
362 integer iFirstIndex = pModalNode->iGetFirstIndex();
363 for (int iCnt = 1; iCnt <= RigidF; iCnt++) {
364 WorkVec.PutRowIndex(iCnt, iFirstIndex + RigidF + iCnt);
365 }
366
367 const Vec3& X0(pModalNode->GetXCurr());
368 const Mat3x3& Rn(pModalNode->GetRCurr());
369 Mat3x3 RR(Rn*Ra);
370
371 // q
372 pq->Put(1, RR.MulTV(X0) - P0);
373
374 Vec3 g(CGR_Rot::Param, RR.MulMT(R0));
375 doublereal d(g.Norm());
376 if (d > std::numeric_limits<doublereal>::epsilon()) {
377 pq->Put(4, g*(2./d*atan(d/2.)));
378
379 } else {
380 pq->Put(4, Zero3);
381 }
382
383 // dot{q}
384 const Vec3& V0(pModalNode->GetVCurr());
385 const Vec3& W0(pModalNode->GetWCurr());
386 pqPrime->Put(1, RR.MulTV(V0));
387 pqPrime->Put(4, RR.MulTV(W0));
388
389 // ddot{q}
390 pqPrime->Put(1, Zero3);
391 pqPrime->Put(4, Zero3);
392 }
393
394 integer iModalIndex = pModalJoint->iGetModalIndex();
395
396 for (unsigned int iCnt = 1; iCnt <= NStModes; iCnt++) {
397 WorkVec.PutRowIndex(RigidF + iCnt, iModalIndex + NStModes + iCnt);
398 }
399
400 integer iAeroIndex = iGetFirstIndex();
401 for (unsigned int iCnt = 1; iCnt <= NAeroStates + 2*NGust; iCnt++) {
402 WorkVec.PutRowIndex(RigidF + NStModes + iCnt, iAeroIndex + iCnt);
403 }
404
405 /* Recupera i vettori {a} e {aP} e {aS} (deformate modali) */
406 for (unsigned int iCnt = 1; iCnt <= NStModes; iCnt++) {
407 pq->PutCoef(iCnt + RigidF, XCurr(iModalIndex + iCnt));
408 pqPrime->PutCoef(iCnt + RigidF, 0);
409 pqSec->PutCoef(iCnt + RigidF, 0);
410 }
411
412 /* Recupera i vettori {xa} e {xaP} */
413 for (unsigned int iCnt = 1; iCnt <= NAeroStates; iCnt++) {
414 pxa->PutCoef(iCnt, XCurr(iAeroIndex + iCnt));
415 pxaPrime->PutCoef(iCnt, 0);
416 }
417
418 for (unsigned int iCnt = 1; iCnt <= 2*NGust; iCnt++) {
419 pgs->PutCoef(iCnt, XCurr(iAeroIndex + NAeroStates + iCnt));
420 pgsPrime->PutCoef(iCnt, 0.);
421 }
422
423 AssVec(WorkVec);
424
425 return WorkVec;
426 }
427
428 /* assemblaggio residuo */
429 void
AssVec(SubVectorHandler & WorkVec)430 AerodynamicModal::AssVec(SubVectorHandler& WorkVec)
431 {
432 DEBUGCOUTFNAME("AerodynamicModal::AssVec");
433
434 /* Dati del nodo rigido */
435 const Vec3& X0(pModalNode->GetXCurr());
436 Vec3 V0(pModalNode->GetVCurr());
437 const Mat3x3& Rn(pModalNode->GetRCurr());
438 Mat3x3 RR(Rn*Ra);
439 Vec3 Vr(V0);
440
441 Vr(1) = 0.0; // perturbazione di velocita', verificare
442
443 doublereal rho, vs, p, T;
444 GetAirProps(X0, rho, vs, p, T); /* p, T are not used yet */
445
446 Vec3 Vair(Zero3);
447 if (fGetAirVelocity(Vair, X0)) {
448 Vr -= Vair;
449 }
450 /* velocita' nel riferimento nodale aerodinamico */
451 V0 = RR.MulTV(Vr);
452 doublereal nV = std::abs(V0(1));
453 MyVectorHandler TmpA(NAeroStates);
454 TmpA.Reset();
455 MyVectorHandler Tmp(RigidF + NStModes);
456 Tmp.Reset();
457 MyVectorHandler TmpP(RigidF + NStModes);
458 TmpP.Reset();
459 MyVectorHandler TmpS(RigidF + NStModes);
460 TmpS.Reset();
461
462 // FIXME: use matrix/vector product?
463 for (unsigned int i = 1; i <= NAeroStates; i++) {
464 for (unsigned int j = 1; j <= NStModes + RigidF; j++) {
465 TmpA.IncCoef(i, pB->operator()(i, j)*(pq->operator()(j)));
466 }
467 }
468
469 pA->MatVecIncMul(TmpA, *pxa);
470
471 /* doublereal CV = Chord/(2*nV); */
472 doublereal qd = 0.5*rho*nV*nV;
473 doublereal qd1 = 0.25*rho*nV*Chord; /* qd*CV */
474 doublereal qd2 = 0.125*rho*Chord*Chord; /* qd*CV*CV */
475
476 for (unsigned int i = 1; i <= NAeroStates; i++) {
477 WorkVec.IncCoef(RigidF + NStModes + i, -pxaPrime->operator()(i) + (2*nV/Chord)*(TmpA(i)));
478 }
479
480 for (unsigned int i = 1; i <= NStModes + RigidF; i++) {
481 for (unsigned int j = 1; j <= NStModes + RigidF; j++) {
482 Tmp.IncCoef(i, pD0->operator()(i, j)*(pq->operator()(j)));
483 TmpP.IncCoef(i, pD1->operator()(i, j)*(pqPrime->operator()(j)));
484 TmpS.IncCoef(i, pD2->operator()(i, j)*(pqSec->operator()(j)));
485 }
486 }
487
488 for (unsigned int i = 1; i <= NStModes + RigidF; i++) {
489 for (unsigned int j = 1; j <= NAeroStates; j++) {
490 Tmp.IncCoef(i, pC->operator()(i, j)*(pxa->operator()(j)));
491 }
492 }
493
494 if (RigidF) {
495 Vec3 F(Zero3);
496 Vec3 M(Zero3);
497
498 for (unsigned int i = 1; i <= 3; i++) {
499 F.Put(i, qd*Tmp(i) + qd1*TmpP(i) + qd2*TmpS(i));
500 M.Put(i, qd*Tmp(i + 3) + qd1*TmpP(i + 3) + qd2*TmpS(i + 3));
501 }
502
503 // std::cout << F << std::endl;
504 F = RR*F;
505 M = RR*M;
506
507 #if 0
508 // tolti perche' altrimenti inizia un Dutch Roll
509 F(1) = 0.;
510 F(2) = 0.;
511 M(1) = 0.;
512 M(3) = 0.;
513 #endif
514
515 WorkVec.Add(1, F);
516 WorkVec.Add(4, M);
517 }
518
519 for (unsigned int i = 1 + RigidF; i <= NStModes + RigidF; i++) {
520 WorkVec.IncCoef(i, qd*Tmp(i) + qd1*TmpP(i) + qd2*TmpS(i));
521 }
522
523 if (NGust) {
524 #if 0 /* unused */
525 doublereal Vyg = Vair(2)/nV;
526 doublereal Vzg = Vair(3)/nV;
527 #endif
528 for (unsigned int i = 1; i <= NAeroStates; i++) {
529 WorkVec.IncCoef(RigidF + NStModes + i,
530 (2*nV/Chord)*pB->operator()(i, RigidF + NStModes + 1)*pgs->operator()(1)
531 + (2*nV/Chord)*pB->operator()(i, RigidF + NStModes + 2)*pgs->operator()(3));
532 }
533
534 if (RigidF) {
535 Vec3 F(Zero3);
536 Vec3 M(Zero3);
537
538 for (unsigned int i = 1; i <= 3; i++) {
539 F.Put(i, qd*pD0->operator()(i, RigidF + NStModes + 1)*pgs->operator()(1)
540 + qd*pD0->operator()(i, RigidF + NStModes + 2)*pgs->operator()(3)
541 + qd1*pD1->operator()(i, RigidF + NStModes + 1)*pgs->operator()(2)
542 + qd1*pD1->operator()(i, RigidF + NStModes + 2)*pgs->operator()(4)
543 + qd2*pD2->operator()(i, RigidF + NStModes + 1)*pgsPrime->operator()(2)
544 + qd2*pD2->operator()(i, RigidF + NStModes + 2)*pgsPrime->operator()(4));
545 M.Put(i, qd*pD0->operator()(i + 3, RigidF + NStModes + 1)*pgs->operator()(1)
546 + qd*pD0->operator()(i + 3, RigidF + NStModes + 2)*pgs->operator()(3)
547 + qd1*pD1->operator()(i + 3, RigidF + NStModes + 1)*pgs->operator()(2)
548 + qd1*pD1->operator()(i + 3, RigidF + NStModes + 2)*pgs->operator()(4)
549 + qd2*pD2->operator()(i + 3, RigidF + NStModes + 1)*pgsPrime->operator()(2)
550 + qd2*pD2->operator()(i + 3, RigidF + NStModes + 2)*pgsPrime->operator()(4));
551 }
552
553 // std::cout << F << std::endl;
554 F = RR*F;
555 M = RR*M;
556 WorkVec.Add(1, F);
557 WorkVec.Add(4, M);
558 }
559
560 for (unsigned int i = 1 + RigidF; i <= NStModes + RigidF; i++) {
561 WorkVec.IncCoef(i, qd*pD0->operator()(i, RigidF + NStModes + 1)*pgs->operator()(1)
562 + qd*pD0->operator()(i, RigidF + NStModes + 2)*pgs->operator()(3)
563 + qd1*pD1->operator()(i, RigidF + NStModes + 1)*pgs->operator()(2)
564 + qd1*pD1->operator()(i, RigidF + NStModes + 2)*pgs->operator()(4)
565 + qd2*pD2->operator()(i, RigidF + NStModes + 1)*pgsPrime->operator()(2)
566 + qd2*pD2->operator()(i, RigidF + NStModes + 2)*pgsPrime->operator()(4));
567 }
568
569 for (unsigned int i = 0; i < NGust; i++) {
570 WorkVec.IncCoef(RigidF + i*2 + 1 + NStModes + NAeroStates,
571 -pgsPrime->operator()(1 + i*2) + pgs->operator()(2 + i*2));
572 if (nV != 0) {
573 WorkVec.IncCoef(RigidF + i*2 + 2 + NStModes + NAeroStates,
574 -pgsPrime->operator()(2 + i*2)
575 - gustVff*gustVff*pgs->operator()(1 + i*2)
576 - 2*gustXi*gustVff*pgs->operator()(2 + i*2)
577 + gustVff*gustVff*Vair(2 + i)/nV);
578
579 } else {
580 WorkVec.IncCoef(RigidF + i*2 + 2 + NStModes + NAeroStates,
581 -pgsPrime->operator()(2 + i*2)
582 - gustVff*gustVff*pgs->operator()(1 + i*2)
583 - 2*gustXi*gustVff*pgs->operator()(2 + i*2));
584 }
585 }
586
587 // std::cout << X0 << std::endl;
588 // std::cout << WorkVec(3) << std::endl;
589 }
590 }
591
592 /* output; si assume che ogni tipo di elemento sappia, attraverso
593 * l'OutputHandler, dove scrivere il proprio output */
594 void
Output(OutputHandler & OH) const595 AerodynamicModal::Output(OutputHandler& OH) const
596 {
597 if (bToBeOutput()) {
598 OH.AeroModals() << std::setw(8) << GetLabel() << " ";
599 for (unsigned int iCnt = 1; iCnt <= NAeroStates; iCnt++) {
600 OH.AeroModals() << " " << pxa->operator()(iCnt);
601 }
602 for (unsigned int iCnt = 1; iCnt <= NAeroStates + 2*NGust; iCnt++) {
603 OH.AeroModals() << " " << pxaPrime->operator()(iCnt);
604 }
605 OH.AeroModals() << std::endl;
606 }
607 }
608
609 /* AerodynamicModal - end */
610
611 Elem *
ReadAerodynamicModal(DataManager * pDM,MBDynParser & HP,const DofOwner * pDO,unsigned int uLabel)612 ReadAerodynamicModal(DataManager* pDM,
613 MBDynParser& HP,
614 const DofOwner* pDO,
615 unsigned int uLabel)
616 {
617 DEBUGCOUTFNAME("ReadAerodynamicModal(" << uLabel << ")" << std::endl);
618
619 /* formato dell'input:
620 *
621 * label,
622 * modal node,
623 * reference modal joint,
624 * rotation of the aerodynamic reference in respect of nodal reference,
625 * reference cord,
626 * number of aerodynamic states,
627 * {rigid, gust, gust filter cut-off frequency}
628 * file name containing state space model matrices;
629 */
630
631 /* giunto modale collegato */
632 const Modal *pModalJoint = pDM->ReadElem<const Modal, const Joint, Elem::JOINT>(HP);
633 ReferenceFrame RF;
634 const StructNode *pModalNode = pModalJoint->pGetModalNode();
635 if (pModalNode) {
636 RF = ReferenceFrame(pModalNode);
637
638 } else {
639 RF = ::AbsRefFrame;
640 }
641
642 Mat3x3 Ra(HP.GetRotRel(RF));
643
644 doublereal Chord = HP.GetReal();
645
646 unsigned int AeroN = HP.GetInt();
647 /* numero modi e FEM */
648 unsigned int NModes = pModalJoint->uGetNModes();
649
650 AerodynamicModal::RigidF_t rigidF = AerodynamicModal::NO_RIGID;
651 if (HP.IsKeyWord("rigid")) {
652 rigidF = AerodynamicModal::RIGID;
653 NModes += rigidF;
654 }
655
656 /* Eventuale raffica */
657 unsigned int GustN = 0;
658 doublereal Vff = 0.;
659 if (HP.IsKeyWord("gust")) {
660 GustN = 2;
661 Vff = HP.GetReal();
662 }
663
664 /* apre il file contenente le matrici A B C D0 D1 D2 */
665 const char *sFileData = HP.GetFileName();
666 std::ifstream fdat(sFileData);
667 DEBUGCOUT("Reading Aerodynamic State Space Matrices from file '"
668 << sFileData << '\'' << std::endl);
669 if (!fdat) {
670 int save_errno = errno;
671 silent_cerr(std::endl << "Unable to open file \"" << sFileData << "\" "
672 "(" << save_errno << ": " << strerror(save_errno) << ")" << std::endl);
673 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
674 }
675
676 SpMapMatrixHandler* pAMat = 0;
677 FullMatrixHandler* pBMat = 0;
678 FullMatrixHandler* pCMat = 0;
679 FullMatrixHandler* pD0Mat = 0;
680 FullMatrixHandler* pD1Mat = 0;
681 FullMatrixHandler* pD2Mat = 0;
682 SAFENEWWITHCONSTRUCTOR(pAMat, SpMapMatrixHandler, SpMapMatrixHandler(AeroN, AeroN));
683 pAMat->Reset();
684 SAFENEWWITHCONSTRUCTOR(pBMat, FullMatrixHandler, FullMatrixHandler(AeroN, NModes + GustN));
685 SAFENEWWITHCONSTRUCTOR(pCMat, FullMatrixHandler, FullMatrixHandler(NModes, AeroN));
686 SAFENEWWITHCONSTRUCTOR(pD0Mat, FullMatrixHandler, FullMatrixHandler(NModes, NModes + GustN));
687 SAFENEWWITHCONSTRUCTOR(pD1Mat, FullMatrixHandler, FullMatrixHandler(NModes, NModes + GustN));
688 SAFENEWWITHCONSTRUCTOR(pD2Mat, FullMatrixHandler, FullMatrixHandler(NModes, NModes + GustN));
689
690 doublereal d;
691 char str[BUFSIZ];
692
693 /* parsing del file */
694 while (!fdat.eof()) {
695 fdat.getline(str, sizeof(str));
696
697 /* legge il primo blocco (HEADER) */
698 if (!strncmp("*** MATRIX A", str, STRLENOF("*** MATRIX A"))) {
699 for (unsigned int iCnt = 1; iCnt <= AeroN; iCnt++) {
700 for (unsigned int jCnt = 1; jCnt <= AeroN; jCnt++) {
701 fdat >> d;
702 pAMat->PutCoef(iCnt, jCnt, d);
703 }
704 }
705
706 /* legge il primo blocco (HEADER) */
707 } else if (!strncmp("*** MATRIX B", str, STRLENOF("*** MATRIX B"))) {
708 for (unsigned int iCnt = 1; iCnt <= AeroN; iCnt++) {
709 for (unsigned int jCnt = 1; jCnt <= NModes + GustN; jCnt++) {
710 fdat >> d;
711 pBMat->PutCoef(iCnt, jCnt, d);
712 }
713 }
714
715 /* legge il primo blocco (HEADER) */
716 } else if (!strncmp("*** MATRIX C", str, STRLENOF("*** MATRIX C"))) {
717 for (unsigned int iCnt = 1; iCnt <= NModes; iCnt++) {
718 for (unsigned int jCnt = 1; jCnt <= AeroN; jCnt++) {
719 fdat >> d;
720 pCMat->PutCoef(iCnt, jCnt, d);
721 }
722 }
723
724 /* legge il primo blocco (HEADER) */
725 } else if (!strncmp("*** MATRIX D0", str, STRLENOF("*** MATRIX D0"))) {
726 for (unsigned int iCnt = 1; iCnt <= NModes; iCnt++) {
727 for (unsigned int jCnt = 1; jCnt <= NModes + GustN; jCnt++) {
728 fdat >> d;
729 pD0Mat->PutCoef(iCnt, jCnt, d);
730 }
731 }
732
733 /* legge il primo blocco (HEADER) */
734 } else if (!strncmp("*** MATRIX D1", str, STRLENOF("*** MATRIX D1"))) {
735 for (unsigned int iCnt = 1; iCnt <= NModes; iCnt++) {
736 for (unsigned int jCnt = 1; jCnt <= NModes + GustN; jCnt++) {
737 fdat >> d;
738 pD1Mat->PutCoef(iCnt, jCnt, d);
739 }
740 }
741
742 /* legge il primo blocco (HEADER) */
743 } else if (!strncmp("*** MATRIX D2", str, STRLENOF("*** MATRIX D2"))) {
744 for (unsigned int iCnt = 1; iCnt <= NModes; iCnt++) {
745 for (unsigned int jCnt = 1; jCnt <= NModes + GustN; jCnt++) {
746 fdat >> d;
747 pD2Mat->PutCoef(iCnt, jCnt, d);
748 }
749 }
750 }
751 }
752 fdat.close();
753
754 flag fOut = pDM->fReadOutput(HP, Elem::AEROMODAL);
755
756 Elem* pEl = 0;
757 SAFENEWWITHCONSTRUCTOR(pEl,
758 AerodynamicModal,
759 AerodynamicModal(uLabel, pModalNode, pModalJoint,
760 Ra, pDO, Chord, NModes - rigidF,
761 AeroN, rigidF, GustN, Vff,
762 pAMat, pBMat, pCMat, pD0Mat, pD1Mat, pD2Mat, fOut));
763
764 /* Se non c'e' il punto e virgola finale */
765 if (HP.IsArg()) {
766 silent_cerr("semicolon expected at line "
767 << HP.GetLineData() << std::endl);
768 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
769 }
770
771 return pEl;
772 } /* End of DataManager::ReadAerodynamicModal() */
773
774