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"
98 #define DYN_ON_SPEED 33 /*20 knots*/
102 #define NZ generic_.n_cg_body_v[2]
108 extern COCKPIT cockpit_;
111 SCALAR interp(SCALAR *y_table, SCALAR *x_table, int Ntable, SCALAR x)
118 /* if x is outside the table, return value at x[0] or x[Ntable-1]*/
122 /* printf("x smaller than x_table[0]: %g %g\n",x,x_table[0]); */
124 else if(x >= x_table[Ntable-1])
127 /* printf("x larger than x_table[N]: %g %g %d\n",x,x_table[NCL-1],Ntable-1); */
129 else /*x is within the table, interpolate linearly to find y value*/
132 while(x_table[i] <= x) {i++;}
133 slope=(y_table[i]-y_table[i-1])/(x_table[i]-x_table[i-1]);
134 /* 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); */
135 y=slope*(x-x_table[i-1]) +y_table[i-1];
142 void aero( SCALAR dt, int Initialize ) {
148 static SCALAR trim_inc = 0.0002;
150 static SCALAR alpha_ind[NCL]={-0.087,0,0.175,0.209,0.24,0.262,0.278,0.303,0.314,0.332,0.367};
151 static SCALAR CLtable[NCL]={-0.14,0.31,1.21,1.376,1.51249,1.591,1.63,1.60878,1.53712,1.376,1.142};
153 /*Note that CLo,Cdo,Cmo will likely change with flap setting so
154 they may not be declared static in the future */
155 /* printf("Initialize= %d\n",Initialize); */
156 /* if (Initialize != 0)
158 /* printf("Initializing aero model...Initialize= %d\n", Initialize);
166 Cda=0.13; /*Not used*/
193 /*nondimensionalization quantities*/
194 /*units here are ft and lbs */
195 cbar=4.9; /*mean aero chord ft*/
196 b=35.8; /*wing span ft */
197 Sw=174; /*wing planform surface area ft^2*/
198 rPiARe=0.054; /*reciprocal of Pi*AR*e*/
204 Cl > 0 => Right wing down
207 elevator > 0 => AND -- aircraft nose down
208 aileron > 0 => right wing up
212 if(Aft_trim) long_trim = long_trim - trim_inc;
213 if(Fwd_trim) long_trim = long_trim + trim_inc;
215 /* 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);
216 */ /*scale pct control to degrees deflection*/
217 if ((Long_control+long_trim) <= 0)
218 elevator=(Long_control+long_trim)*28*DEG_TO_RAD;
220 elevator=(Long_control+long_trim)*23*DEG_TO_RAD;
222 aileron = -1*Lat_control*17.5*DEG_TO_RAD;
223 rudder = -1*Rudder_pedal*16*DEG_TO_RAD;
225 The aileron travel limits are 20 deg. TEU and 15 deg TED
226 but since we don't distinguish between left and right we'll
227 use the average here (17.5 deg)
231 /*calculate rate derivative nondimensionalization (is that a word?) factors */
232 /*hack to avoid divide by zero*/
233 /*the dynamic terms might be negligible at low ground speeds anyway*/
235 if(V_rel_wind > DYN_ON_SPEED)
237 cbar_2V=cbar/(2*V_rel_wind);
238 b_2V=b/(2*V_rel_wind);
247 /*calcuate the qS nondimensionalization factors*/
249 qS=Dynamic_pressure*Sw;
254 /* 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);
255 */ //printf("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);
258 /* sum coefficients */
259 CLwbh = interp(CLtable,alpha_ind,NCL,Alpha);
260 CL = CLo + CLwbh + (CLadot*Alpha_dot + CLq*Theta_dot)*cbar_2V + CLde*elevator;
261 cd = Cdo + rPiARe*CL*CL + Cdde*elevator;
262 cy = Cybeta*Beta + (Cyp*P_body + Cyr*R_body)*b_2V + Cyda*aileron + Cydr*rudder;
264 cm = Cmo + Cma*Alpha + (Cmq*Q_body + Cmadot*Alpha_dot)*cbar_2V + Cmde*(elevator+long_trim);
265 cn = Cnbeta*Beta + (Cnp*P_body + Cnr*R_body)*b_2V + Cnda*aileron + Cndr*rudder;
266 croll=Clbeta*Beta + (Clp*P_body + Clr*R_body)*b_2V + Clda*aileron + Cldr*rudder;
268 /* 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);
269 */ /*calculate wind axes forces*/
274 /* 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);
276 /*calculate moments and body axis forces */
280 /* requires ugly wind-axes to body-axes transform */
281 F_X_aero = F_X_wind*Cos_alpha*Cos_beta - F_Y_wind*Cos_alpha*Sin_beta - F_Z_wind*Sin_alpha;
282 F_Y_aero = F_X_wind*Sin_beta + F_Y_wind*Cos_beta;
283 F_Z_aero = F_X_wind*Sin_alpha*Cos_beta - F_Y_wind*Sin_alpha*Sin_beta + F_Z_wind*Cos_alpha;
285 /*no axes transform here */
286 M_l_aero = croll*qSb;
287 M_m_aero = cm*qScbar;
290 /* 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);
292 /* printf("Fxaero: %7.4f Fyaero: %7.4f Fzaero: %7.4f Weight: %7.4f\n",F_X_aero,F_Y_aero,F_Z_aero,W);
293 *//* printf("Maero: %7.4f Naero: %7.4f Raero: %7.4f\n",M_m_aero,M_n_aero,M_l_aero);