MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
TheodorsenAeroData Class Reference

#include <aerodata_impl.h>

Inheritance diagram for TheodorsenAeroData:
Collaboration diagram for TheodorsenAeroData:

Public Member Functions

 TheodorsenAeroData (int i_p, int i_dim, AeroData *pa, DriveCaller *ptime=0)
 
virtual ~TheodorsenAeroData (void)
 
virtual std::ostream & Restart (std::ostream &out) const
 
virtual void SetAirData (const doublereal &rho, const doublereal &c)
 
virtual void SetSectionData (const doublereal &abscissa, const doublereal &chord, const doublereal &forcepoint, const doublereal &velocitypoint, const doublereal &twist, const doublereal &omega=0.)
 
virtual unsigned int iGetNumDof (void) const
 
virtual DofOrder::Order GetDofType (unsigned int i) const
 
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)
 
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)
 
virtual void AfterConvergence (int i, const VectorHandler &X, const VectorHandler &XP)
 
- Public Member Functions inherited from AeroData
 AeroData (int i_p, int i_dim, UnsteadyModel u=STEADY, DriveCaller *pt=0)
 
virtual ~AeroData (void)
 
std::ostream & RestartUnsteady (std::ostream &out) const
 
virtual int GetForces (int i, const doublereal *W, doublereal *TNG, outa_t &OUTA)
 
virtual int GetForcesJac (int i, const doublereal *W, doublereal *TNG, Mat6x6 &J, outa_t &OUTA)
 
AeroData::UnsteadyModel Unsteady (void) const
 
- Public Member Functions inherited from AeroMemory
 AeroMemory (DriveCaller *pt)
 
virtual ~AeroMemory (void)
 
void Predict (int i, doublereal alpha, doublereal &alf1, doublereal &alf2)
 
void Update (int i)
 
void SetNumPoints (int i)
 
int GetNumPoints (void) const
 

Protected Attributes

integer iParam
 
doublereal d14
 
doublereal d34
 
doublereal chord
 
doublereal a
 
doublereal A1
 
doublereal A2
 
doublereal b1
 
doublereal b2
 
doublerealalpha_pivot
 
doublerealdot_alpha_pivot
 
doublerealdot_alpha
 
doublerealddot_alpha
 
doublerealcfx_0
 
doublerealcfy_0
 
doublerealcmz_0
 
doublerealclalpha
 
doublerealprev_alpha_pivot
 
doublerealprev_dot_alpha
 
doublereal prev_time
 
AeroDatapAeroData
 
- Protected Attributes inherited from AeroData
UnsteadyModel unsteadyflag
 
vam_t VAM
 
doublereal Omega
 
- Protected Attributes inherited from AeroMemory
DriveCallerpTime
 

Additional Inherited Members

- Public Types inherited from AeroData
enum  UnsteadyModel { STEADY = 0, HARRIS = 1, BIELAWA = 2, LAST }
 
enum  {
  VX = 0, VY = 1, VZ = 2, WX = 3,
  WY = 4, WZ = 5, FX = 0, FY = 1,
  FZ = 2, MX = 3, MY = 4, MZ = 5
}
 
- Protected Member Functions inherited from AeroData
int StorageSize (void) const
 
int GetForcesJacForwardDiff_int (int i, const doublereal *W, doublereal *TNG, Mat6x6 &J, outa_t &OUTA)
 
int GetForcesJacCenteredDiff_int (int i, const doublereal *W, doublereal *TNG, Mat6x6 &J, outa_t &OUTA)
 

Detailed Description

Definition at line 156 of file aerodata_impl.h.

Constructor & Destructor Documentation

TheodorsenAeroData::TheodorsenAeroData ( int  i_p,
int  i_dim,
AeroData pa,
DriveCaller ptime = 0 
)

Definition at line 363 of file aerodata_impl.cc.

References ASSERT, AeroData::iGetNumDof(), pAeroData, AeroData::STEADY, and AeroData::Unsteady().

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.),
374 cfx_0(0), cfy_0(0), cmz_0(0), clalpha(0),
376 prev_time(pTime->dGet()),
377 pAeroData(pa)
378 {
380  ASSERT(pAeroData->iGetNumDof() == 0);
381 }
doublereal * alpha_pivot
virtual unsigned int iGetNumDof(void) const
Definition: aerodata.cc:362
doublereal * cfy_0
doublereal * prev_dot_alpha
doublereal * dot_alpha_pivot
doublereal prev_time
doublereal * ddot_alpha
AeroData::UnsteadyModel Unsteady(void) const
Definition: aerodata.cc:208
doublereal * dot_alpha
DriveCaller * pTime
Definition: aerodata.h:67
doublereal * cmz_0
doublereal * cfx_0
doublereal * prev_alpha_pivot
AeroData * pAeroData
#define ASSERT(expression)
Definition: colamd.c:977
virtual doublereal dGet(const doublereal &dVar) const =0
doublereal * clalpha
AeroData(int i_p, int i_dim, UnsteadyModel u=STEADY, DriveCaller *pt=0)
Definition: aerodata.cc:181

Here is the call graph for this function:

TheodorsenAeroData::~TheodorsenAeroData ( void  )
virtual

Definition at line 383 of file aerodata_impl.cc.

References alpha_pivot, pAeroData, SAFEDELETE, and SAFEDELETEARR.

384 {
385  if (alpha_pivot) {
387  }
388 
389  if (pAeroData) {
391  }
392 }
doublereal * alpha_pivot
#define SAFEDELETEARR(pnt)
Definition: mynewmem.h:713
AeroData * pAeroData
#define SAFEDELETE(pnt)
Definition: mynewmem.h:710

Member Function Documentation

void TheodorsenAeroData::AfterConvergence ( int  i,
const VectorHandler X,
const VectorHandler XP 
)
virtual

Reimplemented from AeroData.

Definition at line 914 of file aerodata_impl.cc.

References alpha_pivot, DriveCaller::dGet(), dot_alpha, prev_alpha_pivot, prev_dot_alpha, prev_time, and AeroMemory::pTime.

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 }
doublereal * alpha_pivot
doublereal * prev_dot_alpha
doublereal prev_time
doublereal * dot_alpha
DriveCaller * pTime
Definition: aerodata.h:67
doublereal * prev_alpha_pivot
virtual doublereal dGet(const doublereal &dVar) const =0

Here is the call graph for this function:

void TheodorsenAeroData::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 
)
virtual

Reimplemented from AeroData.

Definition at line 599 of file aerodata_impl.cc.

References a, A1, A2, grad::atan2(), b1, b2, c81_aerod2_u(), outa_t::cd, cfx_0, cfy_0, chord, outa_t::cl, clalpha, outa_t::cm, cmz_0, copysign(), d14, d34, ddot_alpha, vam_t::density, dot_alpha, dot_alpha_pivot, AeroData::GetForcesJac(), Mat3xN::iGetNumCols(), FullSubMatrixHandler::IncCoef(), M_PI, pAeroData, grad::pow(), Mat3xN::Put(), grad::sqrt(), grad::tan(), AeroData::unsteadyflag, AeroData::VAM, AeroData::VX, AeroData::VY, AeroData::VZ, AeroData::WX, AeroData::WY, and AeroData::WZ.

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 }
doublereal cd
Definition: aerodc81.h:66
doublereal density
Definition: aerodc81.h:108
#define M_PI
Definition: gradienttest.cc:67
GradientExpression< BinaryExpr< FuncPow, LhsExpr, RhsExpr > > pow(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2961
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
vam_t VAM
Definition: aerodata.h:117
doublereal * cfy_0
int c81_aerod2_u(doublereal *W, const vam_t *VAM, doublereal *TNG, outa_t *OUTA, c81_data *data, long unsteadyflag)
Definition: aerodc81.c:243
void IncCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: submat.h:683
doublereal * dot_alpha_pivot
doublereal * ddot_alpha
doublereal * dot_alpha
doublereal * cmz_0
doublereal copysign(doublereal x, doublereal y)
Definition: gradient.h:97
doublereal cl
Definition: aerodc81.h:65
doublereal * cfx_0
UnsteadyModel unsteadyflag
Definition: aerodata.h:116
AeroData * pAeroData
integer iGetNumCols(void) const
Definition: matvec3n.h:298
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
doublereal * clalpha
doublereal cm
Definition: aerodc81.h:67
GradientExpression< BinaryExpr< FuncAtan2, LhsExpr, RhsExpr > > atan2(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2962
double doublereal
Definition: colamd.c:52
GradientExpression< UnaryExpr< FuncTan, Expr > > tan(const GradientExpression< Expr > &u)
Definition: gradient.h:2979

Here is the call graph for this function:

void TheodorsenAeroData::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 
)
virtual

Reimplemented from AeroData.

Definition at line 435 of file aerodata_impl.cc.

References a, A1, A2, alpha_pivot, grad::atan2(), b1, b2, vam_t::bc_position, outa_t::cd, cfx_0, cfy_0, vam_t::chord, chord, outa_t::cl, outa_t::clalpha, clalpha, outa_t::cm, cmz_0, copysign(), d14, d34, ddot_alpha, vam_t::density, DriveCaller::dGet(), dot_alpha, dot_alpha_pivot, vam_t::force_position, AeroData::GetForces(), AeroMemory::GetNumPoints(), iParam, MBDYN_EXCEPT_ARGS, pAeroData, grad::pow(), prev_alpha_pivot, prev_dot_alpha, prev_time, AeroMemory::pTime, VectorHandler::PutCoef(), SAFENEWARR, grad::sqrt(), grad::tan(), TheodorsenParams, AeroData::VAM, AeroData::VX, AeroData::VY, AeroData::VZ, AeroData::WX, AeroData::WY, and AeroData::WZ.

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 }
doublereal cd
Definition: aerodc81.h:66
doublereal density
Definition: aerodc81.h:108
GradientExpression< BinaryExpr< FuncPow, LhsExpr, RhsExpr > > pow(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2961
doublereal * alpha_pivot
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
vam_t VAM
Definition: aerodata.h:117
virtual int GetForces(int i, const doublereal *W, doublereal *TNG, outa_t &OUTA)
Definition: aerodata.cc:377
doublereal bc_position
Definition: aerodc81.h:112
doublereal * cfy_0
doublereal * prev_dot_alpha
doublereal * dot_alpha_pivot
doublereal prev_time
doublereal * ddot_alpha
doublereal * dot_alpha
int GetNumPoints(void) const
Definition: aerodata.cc:162
DriveCaller * pTime
Definition: aerodata.h:67
doublereal * cmz_0
doublereal copysign(doublereal x, doublereal y)
Definition: gradient.h:97
doublereal force_position
Definition: aerodc81.h:111
doublereal chord
Definition: aerodc81.h:110
doublereal cl
Definition: aerodc81.h:65
doublereal * cfx_0
doublereal * prev_alpha_pivot
AeroData * pAeroData
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
virtual void PutCoef(integer iRow, const doublereal &dCoef)=0
virtual doublereal dGet(const doublereal &dVar) const =0
doublereal clalpha
Definition: aerodc81.h:69
doublereal * clalpha
#define SAFENEWARR(pnt, item, sz)
Definition: mynewmem.h:701
doublereal cm
Definition: aerodc81.h:67
GradientExpression< BinaryExpr< FuncAtan2, LhsExpr, RhsExpr > > atan2(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2962
double doublereal
Definition: colamd.c:52
GradientExpression< UnaryExpr< FuncTan, Expr > > tan(const GradientExpression< Expr > &u)
Definition: gradient.h:2979
static const doublereal TheodorsenParams[2][4]

Here is the call graph for this function:

DofOrder::Order TheodorsenAeroData::GetDofType ( unsigned int  i) const
virtual

Reimplemented from AeroData.

Definition at line 429 of file aerodata_impl.cc.

References DofOrder::DIFFERENTIAL.

430 {
431  return DofOrder::DIFFERENTIAL;
432 }
unsigned int TheodorsenAeroData::iGetNumDof ( void  ) const
virtual

Reimplemented from AeroData.

Definition at line 423 of file aerodata_impl.cc.

424 {
425  return 2;
426 }
std::ostream & TheodorsenAeroData::Restart ( std::ostream &  out) const
virtual

Implements AeroData.

Definition at line 395 of file aerodata_impl.cc.

References pAeroData, and AeroData::Restart().

396 {
397  out << "theodorsen";
398 
399  return pAeroData->Restart(out);
400 }
virtual std::ostream & Restart(std::ostream &out) const =0
AeroData * pAeroData

Here is the call graph for this function:

void TheodorsenAeroData::SetAirData ( const doublereal rho,
const doublereal c 
)
virtual

Reimplemented from AeroData.

Definition at line 403 of file aerodata_impl.cc.

References pAeroData, and AeroData::SetAirData().

404 {
405  AeroData::SetAirData(rho, c);
406  pAeroData->SetAirData(rho, c);
407 }
virtual void SetAirData(const doublereal &rho, const doublereal &c)
Definition: aerodata.cc:214
AeroData * pAeroData
static std::stack< cleanup * > c
Definition: cleanup.cc:59

Here is the call graph for this function:

void TheodorsenAeroData::SetSectionData ( const doublereal abscissa,
const doublereal chord,
const doublereal forcepoint,
const doublereal velocitypoint,
const doublereal twist,
const doublereal omega = 0. 
)
virtual

Reimplemented from AeroData.

Definition at line 410 of file aerodata_impl.cc.

References pAeroData, and AeroData::SetSectionData().

416 {
417  AeroData::SetSectionData(abscissa, chord, forcepoint, velocitypoint, twist, omega);
418  pAeroData->SetSectionData(abscissa, chord, forcepoint, velocitypoint, twist, omega);
419 }
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
AeroData * pAeroData

Here is the call graph for this function:

Member Data Documentation

doublereal TheodorsenAeroData::a
protected

Definition at line 161 of file aerodata_impl.h.

Referenced by AssJac(), and AssRes().

doublereal TheodorsenAeroData::A1
protected

Definition at line 162 of file aerodata_impl.h.

Referenced by AssJac(), and AssRes().

doublereal TheodorsenAeroData::A2
protected

Definition at line 162 of file aerodata_impl.h.

Referenced by AssJac(), and AssRes().

doublereal* TheodorsenAeroData::alpha_pivot
protected

Definition at line 163 of file aerodata_impl.h.

Referenced by AfterConvergence(), AssRes(), and ~TheodorsenAeroData().

doublereal TheodorsenAeroData::b1
protected

Definition at line 162 of file aerodata_impl.h.

Referenced by AssJac(), and AssRes().

doublereal TheodorsenAeroData::b2
protected

Definition at line 162 of file aerodata_impl.h.

Referenced by AssJac(), and AssRes().

doublereal* TheodorsenAeroData::cfx_0
protected

Definition at line 164 of file aerodata_impl.h.

Referenced by AssJac(), and AssRes().

doublereal * TheodorsenAeroData::cfy_0
protected

Definition at line 164 of file aerodata_impl.h.

Referenced by AssJac(), and AssRes().

doublereal TheodorsenAeroData::chord
protected

Definition at line 160 of file aerodata_impl.h.

Referenced by AssJac(), and AssRes().

doublereal* TheodorsenAeroData::clalpha
protected

Definition at line 165 of file aerodata_impl.h.

Referenced by AssJac(), and AssRes().

doublereal * TheodorsenAeroData::cmz_0
protected

Definition at line 164 of file aerodata_impl.h.

Referenced by AssJac(), and AssRes().

doublereal TheodorsenAeroData::d14
protected

Definition at line 159 of file aerodata_impl.h.

Referenced by AssJac(), and AssRes().

doublereal TheodorsenAeroData::d34
protected

Definition at line 159 of file aerodata_impl.h.

Referenced by AssJac(), and AssRes().

doublereal * TheodorsenAeroData::ddot_alpha
protected

Definition at line 163 of file aerodata_impl.h.

Referenced by AssJac(), and AssRes().

doublereal * TheodorsenAeroData::dot_alpha
protected

Definition at line 163 of file aerodata_impl.h.

Referenced by AfterConvergence(), AssJac(), and AssRes().

doublereal * TheodorsenAeroData::dot_alpha_pivot
protected

Definition at line 163 of file aerodata_impl.h.

Referenced by AssJac(), and AssRes().

integer TheodorsenAeroData::iParam
protected

Definition at line 158 of file aerodata_impl.h.

Referenced by AssRes().

AeroData* TheodorsenAeroData::pAeroData
protected
doublereal* TheodorsenAeroData::prev_alpha_pivot
protected

Definition at line 166 of file aerodata_impl.h.

Referenced by AfterConvergence(), and AssRes().

doublereal * TheodorsenAeroData::prev_dot_alpha
protected

Definition at line 166 of file aerodata_impl.h.

Referenced by AfterConvergence(), and AssRes().

doublereal TheodorsenAeroData::prev_time
protected

Definition at line 167 of file aerodata_impl.h.

Referenced by AfterConvergence(), and AssRes().


The documentation for this class was generated from the following files: