1 /***************************************************************************
5 ----------------------------------------------------------------------------
7 FUNCTION: aerodynamics model based on constant stability derivatives
9 ----------------------------------------------------------------------------
11 MODULE STATUS: developmental
13 ----------------------------------------------------------------------------
15 GENEALOGY: Based on data from:
16 Part 1 of Roskam's S&C text
17 The FAA type certificate data sheet for the 172
18 Various sources on the net
19 John D. Anderson's Intro to Flight text (NACA 2412 data)
20 UIUC's airfoil data web site
22 ----------------------------------------------------------------------------
24 DESIGNED BY: Tony Peden
28 MAINTAINED BY: Tony Peden
30 ----------------------------------------------------------------------------
35 6/10/99 Initial test release
38 ----------------------------------------------------------------------------
48 Croll,Cl rolling moment (yeah, I know. Shoot me.)
51 o constant i.e. not a function of alpha or beta
60 de elevator deflection
67 ----------------------------------------------------------------------------
71 ----------------------------------------------------------------------------
75 ----------------------------------------------------------------------------
79 ----------------------------------------------------------------------------
83 --------------------------------------------------------------------------*/
87 #include "ls_generic.h"
88 #include "ls_cockpit.h"
89 #include "ls_constants.h"
91 #include "c172_aero.h"
99 #define DYN_ON_SPEED 33 /*20 knots*/
103 #define NZ generic_.n_cg_body_v[2]
109 extern COCKPIT cockpit_;
112 static SCALAR interp(SCALAR *y_table, SCALAR *x_table, int Ntable, SCALAR x)
119 /* if x is outside the table, return value at x[0] or x[Ntable-1]*/
123 /* printf("x smaller than x_table[0]: %g %g\n",x,x_table[0]); */
125 else if(x >= x_table[Ntable-1])
127 slope=(y_table[Ntable-1]-y_table[Ntable-2])/(x_table[Ntable-1]-x_table[Ntable-2]);
128 y=slope*(x-x_table[Ntable-1]) +y_table[Ntable-1];
130 /* printf("x larger than x_table[N]: %g %g %d\n",x,x_table[NCL-1],Ntable-1);
132 else /*x is within the table, interpolate linearly to find y value*/
135 while(x_table[i] <= x) {i++;}
136 slope=(y_table[i]-y_table[i-1])/(x_table[i]-x_table[i-1]);
137 /* printf("x: %g, i: %d, cl[i]: %g, cl[i-1]: %g, slope: %g\n",x,i,y_table[i],y_table[i-1],slope); */
138 y=slope*(x-x_table[i-1]) +y_table[i-1];
144 void c172_aero( SCALAR dt, int Initialize ) {
149 static SCALAR lastFlapHandle=0;
152 static SCALAR trim_inc = 0.0002;
154 static SCALAR alpha_ind[NCL]={-0.087,0,0.14,0.21,0.24,0.26,0.28,0.31,0.35};
155 static SCALAR CLtable[NCL]={-0.22,0.25,1.02,1.252,1.354,1.44,1.466,1.298,0.97};
157 static SCALAR flap_ind[Ndf]={0,10,20,30};
158 static SCALAR flap_times[Ndf]={0,4,2,2};
159 static SCALAR dCLf[Ndf]={0,0.20,0.30,0.35};
160 static SCALAR dCdf[Ndf]={0,0.0021,0.0085,0.0191};
161 static SCALAR dCmf[Ndf]={0,-0.0654,-0.0981,-0.114};
163 static SCALAR flap_transit_rate=2.5;
169 /* printf("Initialize= %d\n",Initialize); */
170 /* printf("Initializing aero model...Initialize= %d\n", Initialize);
172 /*nondimensionalization quantities*/
173 /*units here are ft and lbs */
174 cbar=4.9; /*mean aero chord ft*/
175 b=35.8; /*wing span ft */
176 Sw=174; /*wing planform surface area ft^2*/
177 rPiARe=0.054; /*reciprocal of Pi*AR*e*/
178 lbare=3.682; /*elevator moment arm / MAC*/
187 Cda=0.13; /*Not used*/
196 CLde=-Cmde / lbare; /* kinda backwards, huh? */
218 MaxTakeoffWeight=2550;
226 Cl > 0 => Right wing down
229 elevator > 0 => AND -- aircraft nose down
230 aileron > 0 => right wing up
234 /*do weight & balance here since there is no better place*/
239 else if(Weight < 1500)
245 else if(Dx_cg < -0.6)
248 Cg=0.25 - Dx_cg/cbar;
254 if(Flap_handle < flap_ind[0])
257 Flap_handle=flap_ind[0];
258 lastFlapHandle=Flap_handle;
259 Flap_Position=flap_ind[0];
261 else if(Flap_handle > flap_ind[Ndf-1])
264 Flap_handle=flap_ind[fi];
265 lastFlapHandle=Flap_handle;
266 Flap_Position=flap_ind[fi];
272 Flap_Position=Flap_handle;
275 if(Flap_handle != lastFlapHandle)
282 while(flap_ind[fi] < Flap_handle) { fi++; }
283 if(Flap_Position < Flap_handle)
285 if(flap_times[fi] > 0)
286 flap_transit_rate=(flap_ind[fi] - flap_ind[fi-1])/flap_times[fi];
288 flap_transit_rate=(flap_ind[fi] - flap_ind[fi-1])/5;
292 if(flap_times[fi+1] > 0)
293 flap_transit_rate=(flap_ind[fi] - flap_ind[fi+1])/flap_times[fi+1];
295 flap_transit_rate=(flap_ind[fi] - flap_ind[fi+1])/5;
297 if(fabs(Flap_Position - Flap_handle) > dt*flap_transit_rate)
298 Flap_Position+=flap_transit_rate*dt;
302 Flap_Position=Flap_handle;
306 lastFlapHandle=Flap_handle;
309 if(Aft_trim) long_trim = long_trim - trim_inc;
310 if(Fwd_trim) long_trim = long_trim + trim_inc;
312 /* printf("Long_control: %7.4f, long_trim: %7.4f,DEG_TO_RAD: %7.4f, RAD_TO_DEG: %7.4f\n",Long_control,long_trim,DEG_TO_RAD,RAD_TO_DEG);
313 */ /*scale pct control to degrees deflection*/
314 if ((Long_control+Long_trim) <= 0)
315 elevator=(Long_control+Long_trim)*28*DEG_TO_RAD;
317 elevator=(Long_control+Long_trim)*23*DEG_TO_RAD;
319 aileron = -1*Lat_control*17.5*DEG_TO_RAD;
320 rudder = -1*Rudder_pedal*16*DEG_TO_RAD;
322 The aileron travel limits are 20 deg. TEU and 15 deg TED
323 but since we don't distinguish between left and right we'll
324 use the average here (17.5 deg)
328 /*calculate rate derivative nondimensionalization (is that a word?) factors */
329 /*hack to avoid divide by zero*/
330 /*the dynamic terms are negligible at low ground speeds anyway*/
332 /* printf("Vinf: %g, Span: %g\n",V_rel_wind,b);
334 if(V_rel_wind > DYN_ON_SPEED)
336 cbar_2V=cbar/(2*V_rel_wind);
337 b_2V=b/(2*V_rel_wind);
346 /*calcuate the qS nondimensionalization factors*/
348 qS=Dynamic_pressure*Sw;
353 /* printf("aero: Wb: %7.4f, Ub: %7.4f, Alpha: %7.4f, elev: %7.4f, ail: %7.4f, rud: %7.4f, long_trim: %7.4f\n",W_body,U_body,Alpha*RAD_TO_DEG,elevator*RAD_TO_DEG,aileron*RAD_TO_DEG,rudder*RAD_TO_DEG,long_trim*RAD_TO_DEG);
354 printf("aero: Theta: %7.4f, Gamma: %7.4f, Beta: %7.4f, Phi: %7.4f, Psi: %7.4f\n",Theta*RAD_TO_DEG,Gamma_vert_rad*RAD_TO_DEG,Beta*RAD_TO_DEG,Phi*RAD_TO_DEG,Psi*RAD_TO_DEG);
359 /* sum coefficients */
360 CLwbh = interp(CLtable,alpha_ind,NCL,Alpha);
361 /* printf("CLwbh: %g\n",CLwbh);
363 CLo = CLob + interp(dCLf,flap_ind,Ndf,Flap_Position);
364 Cdo = Cdob + interp(dCdf,flap_ind,Ndf,Flap_Position);
365 Cmo = Cmob + interp(dCmf,flap_ind,Ndf,Flap_Position);
367 /* printf("FP: %g\n",Flap_Position);
368 printf("CLo: %g\n",CLo);
369 printf("Cdo: %g\n",Cdo);
370 printf("Cmo: %g\n",Cmo); */
376 CL = CLo + CLwbh + (CLadot*Alpha_dot + CLq*Theta_dot)*cbar_2V + CLde*elevator;
377 cd = Cdo + rPiARe*Ai*Ai*CL*CL + Cdde*elevator;
378 cy = Cybeta*Beta + (Cyp*P_body + Cyr*R_body)*b_2V + Cyda*aileron + Cydr*rudder;
380 cm = Cmo + Cma*Alpha + (Cmq*Q_body + Cmadot*Alpha_dot)*cbar_2V + Cmde*(elevator);
381 cn = Cnbeta*Beta + (Cnp*P_body + Cnr*R_body)*b_2V + Cnda*aileron + Cndr*rudder;
382 croll=Clbeta*Beta + (Clp*P_body + Clr*R_body)*b_2V + Clda*aileron + Cldr*rudder;
384 /* printf("aero: CL: %7.4f, Cd: %7.4f, Cm: %7.4f, Cy: %7.4f, Cn: %7.4f, Cl: %7.4f\n",CL,cd,cm,cy,cn,croll);
387 /*calculate wind axes forces*/
392 /* printf("V_rel_wind: %7.4f, Fxwind: %7.4f Fywind: %7.4f Fzwind: %7.4f\n",V_rel_wind,F_X_wind,F_Y_wind,F_Z_wind);
395 /*calculate moments and body axis forces */
399 /* requires ugly wind-axes to body-axes transform */
400 F_X_aero = F_X_wind*Cos_alpha*Cos_beta - F_Y_wind*Cos_alpha*Sin_beta - F_Z_wind*Sin_alpha;
401 F_Y_aero = F_X_wind*Sin_beta + F_Y_wind*Cos_beta;
402 F_Z_aero = F_X_wind*Sin_alpha*Cos_beta - F_Y_wind*Sin_alpha*Sin_beta + F_Z_wind*Cos_alpha;
404 /*no axes transform here */
405 M_l_aero = croll*qSb;
406 M_m_aero = cm*qScbar;
409 /* printf("I_yy: %7.4f, qScbar: %7.4f, qbar: %7.4f, Sw: %7.4f, cbar: %7.4f, 0.5*rho*V^2: %7.4f\n",I_yy,qScbar,Dynamic_pressure,Sw,cbar,0.5*0.0023081*V_rel_wind*V_rel_wind);
411 printf("Fxaero: %7.4f Fyaero: %7.4f Fzaero: %7.4f Weight: %7.4f\n",F_X_aero,F_Y_aero,F_Z_aero,Weight);
413 printf("Maero: %7.4f Naero: %7.4f Raero: %7.4f\n",M_m_aero,M_n_aero,M_l_aero);