MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
aerodata_impl.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/aerodata_impl.cc,v 1.36 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 
32 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
33 
34 #include "mynewmem.h"
35 #include "aerodata_impl.h"
36 #include "gauss.h"
37 #include "submat.h"
38 #include "aerod2.h"
39 #include "c81data.h"
40 
41 /* STAHRAeroData - begin */
42 
43 STAHRAeroData::STAHRAeroData(int i_p, int i_dim,
45  DriveCaller *ptime)
46 : AeroData(i_p, i_dim, u, ptime),
47 profile(p)
48 {
49  ASSERT(u != AeroData::STEADY ? (ptime != 0) : 1);
50 }
51 
53 {
54  NO_OP;
55 }
56 
57 std::ostream&
58 STAHRAeroData::Restart(std::ostream& out) const
59 {
60  switch (profile) {
61  case 1:
62  out << "NACA0012";
63  break;
64 
65  case 2:
66  out << "RAE9671";
67  break;
68 
69  default:
71  }
72 
73  return RestartUnsteady(out);
74 }
75 
76 int
78  outa_t& OUTA)
79 {
80 
81  switch (unsteadyflag) {
82  case AeroData::HARRIS:
83  case AeroData::BIELAWA:
84  Predict(i, atan2(-W[VY], W[VX]), OUTA.alf1, OUTA.alf2);
85  break;
86 
87  default:
88  break;
89  }
90 
92  __FC_DECL__(aerod2)(const_cast<doublereal *>(W), reinterpret_cast<doublereal *>(&VAM),
93  TNG, reinterpret_cast<doublereal *>(&OUTA), &u, &Omega, &profile);
94 
95  return 0;
96 }
97 
98 int
100 {
101  return AeroData::GetForcesJacForwardDiff_int(i, W, TNG, J, OUTA);
102 }
103 
104 /* STAHRAeroData - end */
105 
106 /* C81AeroData - begin */
107 
108 C81AeroData::C81AeroData(int i_p, int i_dim,
110  const c81_data* d, DriveCaller *ptime)
111 : AeroData(i_p, i_dim, u, ptime),
112 profile(p), data(d)
113 {
114  ASSERT(data != NULL);
115 }
116 
118 {
119  NO_OP;
120 }
121 
122 std::ostream&
123 C81AeroData::Restart(std::ostream& out) const
124 {
125  out << "C81, " << profile;
126 
127  return RestartUnsteady(out);
128 }
129 
130 int
131 C81AeroData::GetForces(int i, const doublereal* W, doublereal* TNG, outa_t& OUTA)
132 {
133  switch (unsteadyflag) {
134  case AeroData::HARRIS:
135  case AeroData::BIELAWA:
136  Predict(i, atan2(-W[VY], W[VX]), OUTA.alf1, OUTA.alf2);
137  break;
138 
139  default:
140  break;
141  }
142 
143  return c81_aerod2_u(const_cast<doublereal *>(W), &VAM, TNG, &OUTA,
144  const_cast<c81_data *>(data), unsteadyflag);
145 }
146 
147 int
149 {
150  return AeroData::GetForcesJacForwardDiff_int(i, W, TNG, J, OUTA);
151 }
152 
153 /* C81AeroData - end */
154 
155 /* C81MultipleAeroData - begin */
156 
158  int i_p, int i_dim,
160  std::vector<unsigned>& p,
161  std::vector<doublereal>& ub,
162  std::vector<const c81_data *>& d,
163  DriveCaller *ptime)
164 : AeroData(i_p, i_dim, u, ptime),
165 profiles(p), upper_bounds(ub), data(d)
166 {
167  ASSERT(!profiles.empty());
168  ASSERT(!upper_bounds.empty());
169  ASSERT(!data.empty());
170  ASSERT(profiles.size() == upper_bounds.size());
171  ASSERT(profiles.size() == data.size());
172 }
173 
175 {
176  NO_OP;
177 }
178 
179 std::ostream&
180 C81MultipleAeroData::Restart(std::ostream& out) const
181 {
182  out << "C81, multiple";
183  for (unsigned i = 0; i < profiles.size(); i++) {
184  out << ", " << profiles[i] << ", " << upper_bounds[i];
185  }
186 
187  return RestartUnsteady(out);
188 }
189 
190 void
192  const doublereal& abscissa,
193  const doublereal& chord,
194  const doublereal& forcepoint,
195  const doublereal& velocitypoint,
196  const doublereal& twist,
197  const doublereal& omega)
198 {
199  ASSERT(abscissa >= -1. && abscissa <= 1.);
200 
201  AeroData::SetSectionData(abscissa, chord, forcepoint, velocitypoint,
202  twist, omega);
203 
204  for (int i = profiles.size() - 1; i--; ) {
205  if (abscissa > upper_bounds[i]) {
206  curr_data = i + 1;
207  return;
208  }
209  }
210 
211  curr_data = 0;
212 }
213 
214 int
216 {
217  switch (unsteadyflag) {
218  case AeroData::HARRIS:
219  case AeroData::BIELAWA:
220  Predict(i, atan2(-W[VY], W[VX]), OUTA.alf1, OUTA.alf2);
221  break;
222 
223  default:
224  break;
225  }
226 
227  ASSERT(i >= 0);
228  ASSERT(unsigned(i) < data.size());
229 
230  return c81_aerod2_u(const_cast<doublereal *>(W), &VAM, TNG, &OUTA,
231  const_cast<c81_data *>(data[curr_data]), unsteadyflag);
232 }
233 
234 int
236 {
237  return AeroData::GetForcesJacForwardDiff_int(i, W, TNG, J, OUTA);
238 }
239 
240 /* C81MultipleAeroData - end */
241 
242 /* C81InterpolatedAeroData - begin */
243 
245  integer i_p, int i_dim,
247  std::vector<unsigned>& p,
248  std::vector<doublereal>& ub,
249  std::vector<const c81_data *>& d,
250  doublereal dcltol,
251  DriveCaller *ptime
252 )
253 : AeroData(i_p, i_dim, u, ptime),
254 profiles(p), upper_bounds(ub), data(d),
255 i_data(i_p*i_dim)
256 {
257  ASSERT(!profiles.empty());
258  ASSERT(!upper_bounds.empty());
259  ASSERT(!data.empty());
260  ASSERT(profiles.size() == upper_bounds.size());
261  ASSERT(profiles.size() == data.size());
262 
263  ASSERT(i_dim >= 1);
264  ASSERT(i_dim <= 3);
265  ASSERT(!i_data.empty());
266 
267  const doublereal bs[3][4] = {
268  { -1., 1. },
269  { -1., 0., 1. },
270  { -1., -1./std::sqrt(3.), 1./std::sqrt(3.), 1. }
271  };
272 
273  int i_point = 0;
274  for (int b = 0; b < i_dim; b++) {
275  GaussDataIterator GDI(i_p);
276  PntWght PW = GDI.GetFirst();
277  do {
278  doublereal dCsi = PW.dGetPnt();
279 
280  doublereal bl = bs[i_dim - 1][b + 1] - bs[i_dim - 1][b];
281  doublereal bm = bs[i_dim - 1][b + 1] + bs[i_dim - 1][b];
282  dCsi = (bm + dCsi*bl)/2.;
283 
284  if (c81_data_merge(data.size(), &data[0], &upper_bounds[0],
285  dCsi, dcltol, &i_data[i_point]))
286  {
287  // function logs error
289  }
290 
291  i_point++;
292  } while (GDI.fGetNext(PW));
293  }
294 }
295 
297 {
298  for (unsigned i = 0; i < i_data.size(); i++) {
300  }
301 }
302 
303 std::ostream&
304 C81InterpolatedAeroData::Restart(std::ostream& out) const
305 {
306  out << "C81, interpolated";
307  for (unsigned i = 0; i < profiles.size(); i++) {
308  out << ", " << profiles[i] << ", " << upper_bounds[i];
309  }
310 
311  return RestartUnsteady(out);
312 }
313 
314 void
316  const doublereal& abscissa,
317  const doublereal& chord,
318  const doublereal& forcepoint,
319  const doublereal& velocitypoint,
320  const doublereal& twist,
321  const doublereal& omega)
322 {
323  ASSERT(abscissa >= -1. && abscissa <= 1.);
324 
325  AeroData::SetSectionData(abscissa, chord, forcepoint, velocitypoint,
326  twist, omega);
327 }
328 
329 int
331 {
332 
333  switch (unsteadyflag) {
334  case AeroData::HARRIS:
335  case AeroData::BIELAWA: {
336  Predict(i, atan2(-W[VY], W[VX]), OUTA.alf1, OUTA.alf2);
337  break;
338  }
339 
340  default:
341  break;
342  }
343 
344  return c81_aerod2_u(const_cast<doublereal *>(W), &VAM, TNG, &OUTA,
345  const_cast<c81_data *>(&i_data[i]), unsteadyflag);
346 }
347 
348 int
350 {
351  return AeroData::GetForcesJacForwardDiff_int(i, W, TNG, J, OUTA);
352 }
353 
354 /* C81InterpolatedAeroData - end */
355 
356 /* TheodorsenAeroData - begin */
357 
358 static const doublereal TheodorsenParams[2][4] = {
359  { 0.165, 0.335, 0.0455, 0.3 },
360  { 0.165, 0.335, 0.041, 0.32 }
361 };
362 
364  int i_p, int i_dim,
365  AeroData *pa,
366  DriveCaller *ptime)
367 : AeroData(i_p, i_dim, STEADY, ptime),
368 iParam(0),
369 d14(0.), d34(0.),
370 chord(-1.),
371 a(0.),
372 A1(0.), A2(0.), b1(0.), b2(0.),
373 alpha_pivot(0), dot_alpha_pivot(0), dot_alpha(0), ddot_alpha(0),
374 cfx_0(0), cfy_0(0), cmz_0(0), clalpha(0),
375 prev_alpha_pivot(0), prev_dot_alpha(0),
376 prev_time(pTime->dGet()),
377 pAeroData(pa)
378 {
380  ASSERT(pAeroData->iGetNumDof() == 0);
381 }
382 
384 {
385  if (alpha_pivot) {
387  }
388 
389  if (pAeroData) {
391  }
392 }
393 
394 std::ostream&
395 TheodorsenAeroData::Restart(std::ostream& out) const
396 {
397  out << "theodorsen";
398 
399  return pAeroData->Restart(out);
400 }
401 
402 void
404 {
405  AeroData::SetAirData(rho, c);
406  pAeroData->SetAirData(rho, c);
407 }
408 
409 void
411  const doublereal& chord,
412  const doublereal& forcepoint,
413  const doublereal& velocitypoint,
414  const doublereal& twist,
415  const doublereal& omega)
416 {
417  AeroData::SetSectionData(abscissa, chord, forcepoint, velocitypoint, twist, omega);
418  pAeroData->SetSectionData(abscissa, chord, forcepoint, velocitypoint, twist, omega);
419 }
420 
421 // aerodynamic models with internal states
422 unsigned int
424 {
425  return 2;
426 }
427 
429 TheodorsenAeroData::GetDofType(unsigned int i) const
430 {
431  return DofOrder::DIFFERENTIAL;
432 }
433 
434 void
436  doublereal dCoef,
437  const VectorHandler& XCurr,
438  const VectorHandler& XPrimeCurr,
439  integer iFirstIndex, integer iFirstSubIndex,
440  int i, const doublereal* W, doublereal* TNG, outa_t& OUTA)
441 {
442  // FIXME: should do it earlier
443  if (alpha_pivot == 0) {
444  doublereal *d = 0;
446 
447 #ifdef HAVE_MEMSET
448  memset(d, '0', 10*GetNumPoints()*sizeof(doublereal));
449 #else // ! HAVE_MEMSET
450  for (int i = 0; i < 10*GetNumPoints(); i++) {
451  d[i] = 0.;
452  }
453 #endif // ! HAVE_MEMSET
454 
455  alpha_pivot = d;
456  d += GetNumPoints();
457  dot_alpha_pivot = d;
458  d += GetNumPoints();
459  dot_alpha = d;
460  d += GetNumPoints();
461  ddot_alpha = d;
462  d += GetNumPoints();
463  cfx_0 = d;
464  d += GetNumPoints();
465  cfy_0 = d;
466  d += GetNumPoints();
467  cmz_0 = d;
468  d += GetNumPoints();
469  clalpha = d;
470  d += GetNumPoints();
471  prev_alpha_pivot = d;
472  d += GetNumPoints();
473  prev_dot_alpha = d;
474  d += GetNumPoints();
475  }
476 
477  doublereal q1 = XCurr(iFirstIndex + 1);
478  doublereal q2 = XCurr(iFirstIndex + 2);
479 
480  doublereal q1p = XPrimeCurr(iFirstIndex + 1);
481  doublereal q2p = XPrimeCurr(iFirstIndex + 2);
482 
484  d34 = VAM.bc_position;
485  if (std::abs(d34 - d14) < std::numeric_limits<doublereal>::epsilon()) {
486  silent_cerr("TheodorsenAeroData::AssRes() [#" << i << "]: aerodynamic center and boundary condition point are almost coincident" << std::endl);
488  }
489  chord = VAM.chord;
490 
491  a = (d14 + d34)/chord;
492 
493  //doublereal Uinf = sqrt(W[VX]*W[VX] + W[VY]*W[VY]);
494  doublereal Uinf = sqrt(W[VX]*W[VX] + (W[VY]+d34*W[WZ])*(W[VY]+d34*W[WZ]));
495 
496  A1 = TheodorsenParams[iParam][0];
497  A2 = TheodorsenParams[iParam][1];
498  b1 = TheodorsenParams[iParam][2];
499  b2 = TheodorsenParams[iParam][3];
500 
501  doublereal u1 = atan2(- W[VY] - W[WZ]*d14, W[VX]);
502  doublereal u2 = atan2(- W[VY] - W[WZ]*d34, W[VX]); //deve essere l'incidenza a 3/4 corda!!!
503 
504  doublereal d = 2*Uinf/chord;
505 
506  doublereal y1 = (A1 + A2)*b1*b2*d*d*q1 + (A1*b1 + A2*b2)*d*q2
507  + (1 - A1 - A2)*u2;
508 
509  doublereal tan_y1 = std::tan(y1);
510  doublereal Vxp2 = Uinf*Uinf - pow(W[VX]*tan_y1, 2)
511  - std::pow(d34*W[WZ], 2);
512 
513  doublereal WW[6];
514  WW[VX] = copysign(std::sqrt(Vxp2), W[VX]);
515  //WW[VY] = -W[VX]*tan_y1 - d34*W[WZ];
516  WW[VY] = -WW[VX]*tan_y1 - d34*W[WZ];
517  WW[VZ] = W[VZ];
518  WW[WX] = W[WX];
519  WW[WY] = W[WY];
520  WW[WZ] = W[WZ];
521 
522  pAeroData->GetForces(i, WW, TNG, OUTA);
523 
524  doublereal UUinf2 = Uinf*Uinf + W[VZ]*W[VZ];
525  doublereal qD = .5*VAM.density*UUinf2;
526  if (qD > std::numeric_limits<doublereal>::epsilon()) {
527  cfx_0[i] = OUTA.cd;
528  cfy_0[i] = OUTA.cl;
529 #if 0
530  cfz_0 = 0.;
531  cmx_0 = 0.;
532  cmy_0 = 0.;
533 #endif
534  cmz_0[i] = OUTA.cm;
535  } else {
536  cfx_0[i] = 0.;
537  cfy_0[i] = 0.;
538 #if 0
539  cfz_0 = 0.;
540  cmx_0 = 0.;
541  cmy_0 = 0.;
542 #endif
543  cmz_0[i] = 0.;
544  }
545 
546  alpha_pivot[i] = (u1*d34 - u2*d14)/(d34 - d14);
547  dot_alpha[i] = Uinf*((u1 - u2)/(d34 - d14));
548 
549  doublereal Delta_t = pTime->dGet() - prev_time;
550  /* controllo che dovrebbe servire solo al primo istante di
551  * di tempo, in cui non calcolo le derivate per differenze
552  * finite per non avere 0 a denominatore. */
553  if (Delta_t > std::numeric_limits<doublereal>::epsilon()) {
554  dot_alpha_pivot[i] = (alpha_pivot[i] - prev_alpha_pivot[i])/Delta_t;
555  ddot_alpha[i] = (dot_alpha[i] - prev_dot_alpha[i])/Delta_t;
556 
557  } else {
558  dot_alpha_pivot[0] = 0.;
559  ddot_alpha[0] = 0.;
560  }
561 
562  clalpha[i] = OUTA.clalpha;
563  if (clalpha[i] > 0.) {
564  //if (Uinf > std::numeric_limits<doublereal>::epsilon()) {
565  if (Uinf > 1.e-5) {
566  doublereal cl = OUTA.cl + clalpha[i]/2/d*(dot_alpha_pivot[i] - a/d*ddot_alpha[i]);
567  doublereal cd = OUTA.cd;
568  doublereal cm = OUTA.cm + clalpha[i]/2/d*(-dot_alpha_pivot[i]/4 + (a/4/d - 1/d/16)*ddot_alpha[i] - dot_alpha[i]/4);
569  doublereal v[3], vp;
570  v[0] = W[0];
571  v[1] = W[1] + d34*W[5];
572  v[2] = W[2] - d34*W[4];
573  vp = sqrt(v[0]*v[0] + v[1]*v[1]);
574  TNG[0] = -qD*chord*(cl*v[1] + cd*v[0])/vp;
575  TNG[1] = qD*chord*(cl*v[0] - cd*v[1])/vp;
576  TNG[5] = qD*chord*chord*cm + d14*TNG[1];
577  } else {
578  TNG[0] = 0.;
579  TNG[1] = 0.;
580  TNG[2] = 0.;
581  TNG[3] = 0.;
582  TNG[4] = 0.;
583  TNG[5] = 0.;
584  }
585  } else {
586  TNG[0] = 0.;
587  TNG[1] = 0.;
588  TNG[2] = 0.;
589  TNG[3] = 0.;
590  TNG[4] = 0.;
591  TNG[5] = 0.;
592  }
593 
594  WorkVec.PutCoef(iFirstSubIndex + 1, -q1p + q2);
595  WorkVec.PutCoef(iFirstSubIndex + 2, -q2p - b1*b2*d*d*q1 - (b1 + b2)*d*q2 + u2);
596 }
597 
598 void
600  doublereal dCoef,
601  const VectorHandler& XCurr,
602  const VectorHandler& XPrimeCurr,
603  integer iFirstIndex, integer iFirstSubIndex,
604  const Mat3xN& vx, const Mat3xN& wx, Mat3xN& fq, Mat3xN& cq,
605  int i, const doublereal* W, doublereal* TNG, Mat6x6& J, outa_t& OUTA)
606 {
607 #ifdef DEBUG_JACOBIAN_UNSTEADY
608  doublereal q1 = XCurr(iFirstIndex + 1);
609  doublereal q2 = XCurr(iFirstIndex + 2);
610 #endif // DEBUG_JACOBIAN_UNSTEADY
611 
612  doublereal Uinf = sqrt(W[VX]*W[VX] + W[VY]*W[VY]);
613  doublereal d = 2*Uinf/chord;
614 #ifdef DEBUG_JACOBIAN_UNSTEADY
615  doublereal UUinf2 = Uinf*Uinf + W[VZ]*W[VZ];
616 #endif // DEBUG_JACOBIAN_UNSTEADY
617 #if 0
618  doublereal qD = .5*VAM.density*UUinf2;
619 #endif
620 
621  // doublereal u1 = atan2(- W[VY] - W[WZ]*d14, W[VX]);
622 #ifdef DEBUG_JACOBIAN_UNSTEADY
623  doublereal u2 = atan2(- W[VY] - W[WZ]*d34, W[VX]);
624 #endif // DEBUG_JACOBIAN_UNSTEADY
625 
626 #if 0
627  doublereal y1 = (A1 + A2)*b1*b2*d*d*q1 + (A1*b1 + A2*b2)*d*q2
628  + (1 - A1 - A2)*u2;
629 #endif
630 
631  WorkMat.IncCoef(iFirstSubIndex + 1, iFirstSubIndex + 1, 1.);
632  WorkMat.IncCoef(iFirstSubIndex + 2, iFirstSubIndex + 2, 1.);
633 
634  WorkMat.IncCoef(iFirstSubIndex + 1, iFirstSubIndex + 2, -dCoef);
635  WorkMat.IncCoef(iFirstSubIndex + 2, iFirstSubIndex + 1, dCoef*b1*b2*d*d);
636  WorkMat.IncCoef(iFirstSubIndex + 2, iFirstSubIndex + 2, dCoef*(b1 + b2)*d);
637 
638 #if 0
639  //if (Uinf > std::numeric_limits<doublereal>::epsilon()) {
640  if (Uinf > 1.e-2) {
641  /* computing the matrix g_{/\tilde{v}} - dimension: 2x3, where 2 is n_states */
642  doublereal dU_V_11, dU_V_12;
643 
644  doublereal dDen34 = W[VX]*W[VX] + W[VY]*W[VY] + W[WZ]*W[WZ]*d34*d34 + 2*W[VY]*W[WZ]*d34;
645 
646  dU_V_11 = (W[VY] + W[WZ]*d34)/dDen34;
647  dU_V_12 = -W[VX]/dDen34;
648 
649  doublereal dG_V_21, dG_V_22;
650 
651  dG_V_21 = (b1*b2*d*4*q1/chord + (b1+b2)*2*q2/chord)*W[VX]/Uinf - dU_V_11;
652  dG_V_22 = (b1*b2*d*4*q1/chord + (b1+b2)*2*q2/chord)*W[VY]/Uinf - dU_V_12;
653 
654 
655  /* computing the matrix g_{/\tilde{\omega}} - dimension: 2x3, where 2 is n_states */
656  doublereal dU_W_13;
657 
658  dU_W_13 = - W[VX]*d34/dDen34;
659 
660  doublereal dG_W_23;
661 
662  dG_W_23 = -dU_W_13;
663 
664  /* assembling the jacobian */
665  int iIndexColumn;
666  /* faccio i calcoli tenendo conto della sparsit√† della matrice g_{\/\tilde{v}} */
667  #if 1
668  for (iIndexColumn = 1; iIndexColumn <= vx.iGetNumCols(); iIndexColumn++){
669  WorkMat.IncCoef(iFirstSubIndex+2, iIndexColumn, dG_V_21*vx(1,iIndexColumn) + dG_V_22*vx(2,iIndexColumn));
670  WorkMat.IncCoef(iFirstSubIndex+2, iIndexColumn, dG_W_23*wx(3,iIndexColumn));
671  }
672  #endif
673  /* computing the matrix fq (f_a_{/q}) 3x2 */
674 
675  /* computing the derivative of the aerodynamic coefficient in the lookup table
676  * using the finite difference method */
677  /* perturbazione per calcolo delle derivare con le differenze finite */
678  doublereal dDeltaY1 = 1.*M_PI/180.;
679  doublereal y1d = y1 + dDeltaY1;
680  doublereal tan_y1 = std::tan(y1d);
681  doublereal Vxp2 = Uinf*Uinf - pow(W[VX]*tan_y1, 2) - std::pow(d34*W[WZ], 2);
682 
683  doublereal WW[6];
684  WW[VX] = copysign(std::sqrt(Vxp2), W[VX]);
685  WW[VY] = -W[VX]*tan_y1 - d34*W[WZ];
686  WW[VZ] = W[VZ];
687  WW[WX] = W[WX];
688  WW[WY] = W[WY];
689  WW[WZ] = W[WZ];
690 
691  c81_aerod2_u(WW, &VAM, TNG, &OUTA, const_cast<c81_data *>(data), unsteadyflag);
692  doublereal cx_1, cy_1, cmz_1;
693  if (qD > std::numeric_limits<doublereal>::epsilon()) {
694  cx_1 = OUTA.cd;
695  cy_1 = OUTA.cl;
696  cmz_1 = OUTA.cm;
697  }else{
698  cx_1 = 0.;
699  cy_1 = 0.;
700  cmz_1 = 0.;
701  }
702 
703  doublereal dCd_alpha = (cx_1-cfx_0[i])/dDeltaY1;
704  doublereal dCl_alpha = (cy_1-cfy_0[i])/dDeltaY1;
705  doublereal dCm_alpha = (cmz_1-cmz_0[i])/dDeltaY1;
706 
707  doublereal C11 = (A1+A2)*b1*b2*d*d;
708  doublereal C12 = (A1*b1+A2*b2)*d;
709 
710  doublereal qDc = qD*chord;
711  /* ATTENZIONE: le derivate aerodinamiche calcolate (vedi tecman)
712  * sono le dirivate di L D e M34 e vanno quindi routate per avere
713  * le derivate di TNG */
714  doublereal sin_alpha = (W[VY]+d34*W[WZ])/sqrt(W[VX]*W[VX]+W[VY]*W[VY]);
715  doublereal cos_alpha = (W[VX])/sqrt(W[VX]*W[VX]+W[VY]*W[VY]);
716 
717  doublereal D_q1 = qDc*dCd_alpha*C11;
718  doublereal D_q2 = qDc*dCd_alpha*C12;
719  doublereal L_q1 = qDc*dCl_alpha*C11;
720  doublereal L_q2 = qDc*dCl_alpha*C12;
721 
722  doublereal M_q1 = qDc*chord*dCm_alpha*C11;
723  doublereal M_q2 = qDc*chord*dCm_alpha*C12;
724  #if 0
725  fq.Put(1, 1, (-L_q1*sin_alpha -D_q1*cos_alpha));
726  fq.Put(1, 2, (-L_q2*sin_alpha -D_q2*cos_alpha));
727  fq.Put(2, 1, (L_q1*cos_alpha -D_q1*sin_alpha));
728  fq.Put(2, 2, (L_q2*cos_alpha -D_q2*sin_alpha));
729 
730  cq.Put(3, 1, (M_q1 +d14*fq(2,1)));
731  cq.Put(3, 2, (M_q2 +d14*fq(2,2)));
732  #endif
733  /* computing the J matrix */
734  /* f_{/\tilde{v}} */
735  doublereal dY_V_11, dY_V_12;
736 
737  /* (C_{/Uinf}*q + D_{/Uinf}*u)*Uinf_{v} */
738  dY_V_11 = (((A1+A2)*b1*b2*d*(4./chord)*q1 + (A1*b1+A2*b2)*(2./chord)*q2)*W[VX]/Uinf) + ((1-A1-A2)*dU_V_11);
739  dY_V_12 = (((A1+A2)*b1*b2*d*(4./chord)*q1 + (A1*b1+A2*b2)*(2./chord)*q2)*W[VY]/Uinf) + ((1-A1-A2)*dU_V_12);
740 
741  doublereal dCfy_Uinf = 0.5*clalpha[i]*(-dot_alpha_pivot[i] + (chord*a*ddot_alpha[i])/Uinf)/(d*Uinf);
742  doublereal rho = VAM.density;
743 
744  doublereal cy = cfy_0[i] + clalpha[i]/2/d*( dot_alpha_pivot[i] - a/d*ddot_alpha[i]);
745 
746  Mat6x6 Jaero = 0.;
747  Jaero(1,1) = rho*chord*cfx_0[i]*W[VX] + qDc*( dCd_alpha*dY_V_11);
748  Jaero(1,2) = rho*chord*cfx_0[i]*W[VY] + qDc*( dCd_alpha*dY_V_12);
749  Jaero(1,3) = rho*chord*cfx_0[i]*W[VZ];
750  Jaero(2,1) = rho*chord*cy*W[VX] + qDc*( dCl_alpha*dY_V_11 + dCfy_Uinf*W[VX]/Uinf);
751  Jaero(2,2) = rho*chord*cy*W[VY] + qDc*( dCl_alpha*dY_V_12 + dCfy_Uinf*W[VY]/Uinf );
752  Jaero(2,3) = rho*chord*cy*W[VZ];
753  Jaero(3,1) = rho*chord*cfz_0*W[VX];
754  Jaero(3,1) = rho*chord*cfz_0*W[VY];
755  Jaero(3,3) = rho*chord*cfz_0*W[VZ];
756 
757  /* f_{/\tilde{omega}} */
758  Jaero(1,6) = qDc*dCd_alpha*(1.-A1-A2)*dU_W_13;
759  Jaero(2,6) = qDc*dCl_alpha*(1.-A1-A2)*dU_W_13;
760 
761  /* c_{/\tilde{v}} */
762  doublereal dCmz_Uinf = 0.5*clalpha[i]*(dot_alpha_pivot[i] - ( ((chord*a)/(Uinf)) - ((chord)/(4*Uinf)) )*ddot_alpha[i] + dot_alpha[i])/(d*4*Uinf);
763 
764  doublereal cmz = cmz_0 + clalpha[i]/2/d*(-dot_alpha_pivot[i]/4 + (a/4/d - 1/d/16)*ddot_alpha[i] - dot_alpha[i]/4);
765 
766  Jaero(4,1) = rho*chord*chord*cmx_0*W[VX];
767  Jaero(4,2) = rho*chord*chord*cmx_0*W[VY];
768  Jaero(4,3) = rho*chord*chord*cmx_0*W[VZ];
769  Jaero(5,1) = rho*chord*chord*cmy_0*W[VX];
770  Jaero(5,2) = rho*chord*chord*cmy_0*W[VY];
771  Jaero(5,3) = rho*chord*chord*cmy_0*W[VZ];
772  Jaero(6,1) = rho*chord*chord*cmz*W[VX] +qDc*chord*( dCm_alpha*dY_V_11 + dCmz_Uinf*W[VX]/Uinf );
773  Jaero(6,2) = rho*chord*chord*cmz*W[VY] +qDc*chord*( dCm_alpha*dY_V_12 + dCmz_Uinf*W[VY]/Uinf );
774  Jaero(6,3) = rho*chord*chord*cmz*W[VZ];
775 
776  /* c_{/\tilde{omega}} */
777  Jaero(6,6) = qDc*chord*dCm_alpha*(1.-A1-A2)*dU_W_13;
778 
779  /* Jaero √® lo jacobiano delle forze aerodinamiche L, D e M34
780  * occorre ruotarlo per avere lo jacobiano corretto J.
781  * Solo le righe 1 2 e 6 sono coivolte nella rotazione */
782  J = Jaero;
783  /* calcolo le forze aerodinamiche */
784  doublereal Drag = qDc*cfx_0[i];
785  doublereal Lift = qDc*cfy_0[i] + qDc*clalpha[i]/2/d*(dot_alpha_pivot[i] - a/d*ddot_alpha[i]);
786  /* calcolo le derivate del cos e sen dell'angolo di rotazione */
787  doublereal den = (W[VX]*W[VX]+W[VY]*W[VY])*sqrt((W[VX]*W[VX]+W[VY]*W[VY]));
788  doublereal sin_alpha_vx = -W[VX]*(W[VY]+d34*W[WZ])/den;
789  doublereal sin_alpha_vy = ( W[VX]*W[VX] - d34*W[WZ]*W[VY] )/den;
790  doublereal sin_alpha_wz = d34/sqrt(W[VX]*W[VX]+W[VY]*W[VY]);
791  doublereal cos_alpha_vx = ( W[VY]*W[VY] )/den;
792  doublereal cos_alpha_vy = -( W[VX]*W[VY] )/den;
793  #if 1
794  J(1,1) = -Jaero(2,1)*sin_alpha -Jaero(1,1)*cos_alpha - Lift*sin_alpha_vx - Drag*cos_alpha_vx;
795  J(1,2) = -Jaero(2,2)*sin_alpha -Jaero(1,2)*cos_alpha - Lift*sin_alpha_vy - Drag*cos_alpha_vy;
796  J(1,3) = -Jaero(2,3)*sin_alpha -Jaero(1,3)*cos_alpha;
797  J(1,6) = -Jaero(2,6)*sin_alpha -Jaero(1,6)*cos_alpha - Lift*sin_alpha_wz;
798 
799  J(2,1) = Jaero(2,1)*cos_alpha -Jaero(1,1)*sin_alpha + Lift*cos_alpha_vx - Drag*sin_alpha_vx;
800  J(2,2) = Jaero(2,2)*cos_alpha -Jaero(1,2)*sin_alpha + Lift*cos_alpha_vy - Drag*sin_alpha_vy;
801  J(2,3) = Jaero(2,3)*cos_alpha -Jaero(1,3)*sin_alpha;
802  J(2,6) = Jaero(2,6)*cos_alpha -Jaero(1,6)*sin_alpha - Drag*sin_alpha_wz;
803 
804  J(6,1) = Jaero(6,1) + d14*J(2,1);
805  J(6,2) = Jaero(6,2) + d14*J(2,2);
806  J(6,3) = Jaero(6,3) + d14*J(2,3);
807  J(6,4) = d14*J(2,4);
808  J(6,5) = d14*J(2,5);
809  J(6,6) = Jaero(6,6) + d14*J(2,6);
810  #endif
811  }
812 #endif
813 #ifdef DEBUG_JACOBIAN_UNSTEADY
814  printf("G/v matrix\n");
815  printf("%lf %lf %lf\n", 0., 0., 0.);
816  printf("%lf %lf %lf\n", dG_V_21, dG_V_22, 0.);
817 
818  printf("G/w matrix\n");
819  printf("%lf %lf %lf\n", 0., 0., 0.);
820  printf("%lf %lf %lf\n", 0., 0., dG_W_23);
821 
822  printf("f/q matrix\n");
823  printf("%lf %lf \n", fq(1,1), fq(1,2));
824  printf("%lf %lf \n", fq(2,1), fq(2,2));
825  printf("%lf %lf \n", 0., 0.);
826 
827  printf("c/q matrix\n");
828  printf("%lf %lf\n", 0., 0.);
829  printf("%lf %lf\n", 0., 0.);
830  printf("%lf %lf\n", cq(3,1), cq(3,2));
831 
832  printf("J matrix\n");
833  for( int iii=1; iii<=6; iii++){
834  for( int jjj=1; jjj<=6; jjj++){
835  printf("%lf ", J(iii,jjj));
836  }
837  printf("\n");
838  }
839 
840  FILE *fd;
841  fd = fopen("X.mat","w");
842  printf("DATI X\n");
843  for( int iii=1; iii<=2; iii++){
844  printf("%lf ", XCurr(iFirstIndex+iii));
845  fprintf(fd,"%15.7e ", XCurr(iFirstIndex+iii));
846  }
847  fclose(fd);
848  fd = fopen("XP.mat","w");
849  printf("\n");
850  printf("DATI XP\n");
851  for( int iii=1; iii<=2; iii++){
852  printf("%lf ", XPrimeCurr(iFirstIndex+iii));
853  fprintf(fd,"%15.7e ", XPrimeCurr(iFirstIndex+iii));
854  }
855  printf("\n");
856  fclose(fd);
857  fd = fopen("W.mat","w");
858  printf("DATI W\n");
859  printf("%lf %lf %lf %lf %lf %lf", W[VX], W[VY], W[VZ], W[WX], W[WY], W[WZ]);
860  fprintf(fd,"%15.7e %15.7e %15.7e %15.7e %15.7e %15.7e", W[VX], W[VY], W[VZ], W[WX], W[WY], W[WZ]);
861  fclose(fd);
862  printf("\n");
863  printf("PARAMETRI \n");
864  printf("d14 %lf ", d14);
865  printf("\nd34 %lf ", d34);
866  printf("\nA1 %lf ", A1);
867  printf("\nA2 %lf ", A2);
868  printf("\nb1 %lf ", b1);
869  printf("\nb2 %lf ", b2);
870  printf("\nchord %lf ", chord);
871  printf("\na %lf ", a);
872  printf("\nrho %lf ", rho);
873  printf("\ncfx_0 %lf ", cfx_0[i]);
874  fd = fopen("cfx_0.mat","w");
875  fprintf(fd,"%15.7e", cfx_0[i]);
876  fclose(fd);
877  printf("\ncfy_0 %lf ", cfy_0[i]);
878  fd = fopen("cfy_0.mat","w");
879  fprintf(fd,"%15.7e ", cfy_0[i]);
880  fclose(fd);
881  printf("\ncmz_0 %lf ", cmz_0[i]);
882  printf("\ndCl_alpha %lf ", dCl_alpha);
883  fd = fopen("dCl_alpha.mat","w");
884  fprintf(fd,"%15.7e", dCl_alpha);
885  fclose(fd);
886  printf("\ndCd_alpha %lf", dCd_alpha);
887  fd = fopen("dCd_alpha.mat","w");
888  fprintf(fd,"%15.7e", dCd_alpha);
889  fclose(fd);
890  printf("\ndCm_alpha %lf", dCm_alpha);
891  printf("\nclalpha %lf ", clalpha[i]);
892  printf("\n q1 q2 i %lf %lf %d", q1, q2, i);
893  fd = fopen("ddot_alpha.mat","w");
894  fprintf(fd,"%15.7e", ddot_alpha[i]);
895  printf("%lf", ddot_alpha[i]);
896  fclose(fd);
897  fd = fopen("dot_alpha.mat","w");
898  fprintf(fd,"%15.7e", dot_alpha[i]);
899  printf("%lf", dot_alpha[i]);
900  fclose(fd);
901  fd = fopen("dot_alpha_pivot.mat","w");
902  fprintf(fd,"%15.7e", dot_alpha_pivot[i]);
903  printf("%lf", dot_alpha_pivot[i]);
904  fclose(fd);
905  getchar();
906 #endif /* DEBUG_JACOBIAN_UNSTEADY */
907 
908  pAeroData->GetForcesJac(i, W, TNG, J, OUTA);
909 
910  // probably, we need to reset fq, cq
911 }
912 
913 void
915  const VectorHandler& X, const VectorHandler& XP )
916 {
917  /* aggiorno i valori delle variabili per il calcolo
918  * della derivata attraverso le differenze finite */
920  prev_dot_alpha[i] = dot_alpha[i];
921 
922  // same for all
923  prev_time = pTime->dGet();
924 }
925 
926 /* TheodorsenAeroData - end */
std::ostream & Restart(std::ostream &out) const
virtual DofOrder::Order GetDofType(unsigned int i) const
int GetForcesJac(int i, const doublereal *W, doublereal *TNG, Mat6x6 &J, outa_t &OUTA)
virtual void AfterConvergence(int i, const VectorHandler &X, const VectorHandler &XP)
doublereal cd
Definition: aerodc81.h:66
virtual ~TheodorsenAeroData(void)
virtual int GetForcesJac(int i, const doublereal *W, doublereal *TNG, Mat6x6 &J, outa_t &OUTA)
const c81_data * data
Definition: aerodata_impl.h:63
doublereal density
Definition: aerodc81.h:108
virtual void SetSectionData(const doublereal &abscissa, const doublereal &chord, const doublereal &forcepoint, const doublereal &velocitypoint, const doublereal &twist, const doublereal &omega=0.)
Definition: aerodata.cc:237
#define M_PI
Definition: gradienttest.cc:67
std::ostream & RestartUnsteady(std::ostream &out) const
Definition: aerodata.cc:252
GradientExpression< BinaryExpr< FuncPow, LhsExpr, RhsExpr > > pow(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2961
doublereal * alpha_pivot
virtual int GetForcesJac(int i, const doublereal *W, doublereal *TNG, Mat6x6 &J, outa_t &OUTA)
Definition: aerodata.cc:384
void Put(int i, integer j, const doublereal &d)
Definition: matvec3n.h:306
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
std::vector< unsigned > profiles
Definition: aerodata_impl.h:90
virtual void SetSectionData(const doublereal &abscissa, const doublereal &chord, const doublereal &forcepoint, const doublereal &velocitypoint, const doublereal &twist, const doublereal &omega=0.)
vam_t VAM
Definition: aerodata.h:117
virtual void AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr, integer iFirstIndex, integer iFirstSubIndex, int i, const doublereal *W, doublereal *TNG, outa_t &OUTA)
std::vector< doublereal > upper_bounds
std::vector< c81_data > i_data
std::ostream & Restart(std::ostream &out) const
void c81_data_destroy(c81_data *data)
Definition: c81data.cc:276
int GetForcesJac(int i, const doublereal *W, doublereal *TNG, Mat6x6 &J, outa_t &OUTA)
flag fGetNext(doublereal &d, integer i=0) const
Definition: gauss.cc:205
virtual unsigned int iGetNumDof(void) const
Definition: aerodata.cc:362
PntWght GetFirst(void) const
Definition: gauss.cc:198
#define SAFEDELETEARR(pnt)
Definition: mynewmem.h:713
virtual int GetForces(int i, const doublereal *W, doublereal *TNG, outa_t &OUTA)
Definition: aerodata.cc:377
STAHRAeroData(int i_p, int i_dim, AeroData::UnsteadyModel u, integer p, DriveCaller *ptime=0)
int GetForces(int i, const doublereal *W, doublereal *TNG, outa_t &OUTA)
virtual ~C81AeroData(void)
void SetSectionData(const doublereal &abscissa, const doublereal &chord, const doublereal &forcepoint, const doublereal &velocitypoint, const doublereal &twist, const doublereal &omega=0.)
C81AeroData(int i_p, int i_dim, AeroData::UnsteadyModel u, integer p, const c81_data *d, DriveCaller *ptime=0)
doublereal alf1
Definition: aerodc81.h:70
int GetForcesJac(int i, const doublereal *W, doublereal *TNG, Mat6x6 &J, outa_t &OUTA)
doublereal bc_position
Definition: aerodc81.h:112
doublereal * cfy_0
doublereal * prev_dot_alpha
#define NO_OP
Definition: myassert.h:74
int c81_aerod2_u(doublereal *W, const vam_t *VAM, doublereal *TNG, outa_t *OUTA, c81_data *data, long unsteadyflag)
Definition: aerodc81.c:243
virtual std::ostream & Restart(std::ostream &out) const
void IncCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: submat.h:683
virtual void SetAirData(const doublereal &rho, const doublereal &c)
Definition: aerodata.cc:214
std::vector< unsigned > profiles
doublereal * dot_alpha_pivot
doublereal prev_time
doublereal * ddot_alpha
doublereal Omega
Definition: aerodata.h:118
virtual std::ostream & Restart(std::ostream &out) const =0
AeroData::UnsteadyModel Unsteady(void) const
Definition: aerodata.cc:208
doublereal * dot_alpha
void Predict(int i, doublereal alpha, doublereal &alf1, doublereal &alf2)
Definition: aerodata.cc:65
int GetNumPoints(void) const
Definition: aerodata.cc:162
virtual int GetForces(int i, const doublereal *W, doublereal *TNG, outa_t &OUTA)
doublereal alf2
Definition: aerodc81.h:71
C81InterpolatedAeroData(int i_p, int i_dim, AeroData::UnsteadyModel u, std::vector< unsigned > &p, std::vector< doublereal > &ub, std::vector< const c81_data * > &d, doublereal dcltol, DriveCaller *ptime=0)
std::vector< doublereal > upper_bounds
Definition: aerodata_impl.h:91
DriveCaller * pTime
Definition: aerodata.h:67
doublereal * cmz_0
doublereal copysign(doublereal x, doublereal y)
Definition: gradient.h:97
int GetForcesJacForwardDiff_int(int i, const doublereal *W, doublereal *TNG, Mat6x6 &J, outa_t &OUTA)
Definition: aerodata.cc:271
doublereal force_position
Definition: aerodc81.h:111
C81MultipleAeroData(int i_p, int i_dim, AeroData::UnsteadyModel u, std::vector< unsigned > &p, std::vector< doublereal > &ub, std::vector< const c81_data * > &d, DriveCaller *ptime=0)
void SetSectionData(const doublereal &abscissa, const doublereal &chord, const doublereal &forcepoint, const doublereal &velocitypoint, const doublereal &twist, const doublereal &omega=0.)
virtual ~STAHRAeroData(void)
std::vector< const c81_data * > data
int GetForces(int i, const doublereal *W, doublereal *TNG, outa_t &OUTA)
doublereal chord
Definition: aerodc81.h:110
doublereal cl
Definition: aerodc81.h:65
doublereal * cfx_0
UnsteadyModel
Definition: aerodata.h:89
doublereal * prev_alpha_pivot
UnsteadyModel unsteadyflag
Definition: aerodata.h:116
AeroData * pAeroData
int c81_data_merge(unsigned ndata, const c81_data **data, const doublereal *upper_bounds, doublereal dCsi, doublereal dcltol, c81_data *i_data)
Definition: c81data.cc:396
integer iGetNumCols(void) const
Definition: matvec3n.h:298
std::ostream & Restart(std::ostream &out) const
virtual void SetAirData(const doublereal &rho, const doublereal &c)
#define ASSERT(expression)
Definition: colamd.c:977
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
virtual void PutCoef(integer iRow, const doublereal &dCoef)=0
int GetForces(int i, const doublereal *W, doublereal *TNG, outa_t &OUTA)
integer profile
Definition: aerodata_impl.h:62
virtual doublereal dGet(const doublereal &dVar) const =0
doublereal clalpha
Definition: aerodc81.h:69
static std::stack< cleanup * > c
Definition: cleanup.cc:59
doublereal * clalpha
doublereal dGetPnt(void) const
Definition: gauss.h:55
#define SAFENEWARR(pnt, item, sz)
Definition: mynewmem.h:701
static const doublereal a
Definition: hfluid_.h:289
std::vector< const c81_data * > data
Definition: aerodata_impl.h:92
virtual unsigned int iGetNumDof(void) const
doublereal cm
Definition: aerodc81.h:67
virtual void AssJac(FullSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr, integer iFirstIndex, integer iFirstSubIndex, const Mat3xN &vx, const Mat3xN &wx, Mat3xN &fq, Mat3xN &cq, int i, const doublereal *W, doublereal *TNG, Mat6x6 &J, outa_t &OUTA)
GradientExpression< BinaryExpr< FuncAtan2, LhsExpr, RhsExpr > > atan2(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2962
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51
GradientExpression< UnaryExpr< FuncTan, Expr > > tan(const GradientExpression< Expr > &u)
Definition: gradient.h:2979
integer profile
Definition: aerodata_impl.h:41
static const doublereal TheodorsenParams[2][4]
int aerod2(doublereal *w, doublereal *vam, doublereal *tng, doublereal *outa, integer *inst, doublereal *rspeed, integer *ipr)
TheodorsenAeroData(int i_p, int i_dim, AeroData *pa, DriveCaller *ptime=0)
#define SAFEDELETE(pnt)
Definition: mynewmem.h:710
Definition: gauss.h:50
virtual std::ostream & Restart(std::ostream &out) const