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"
96 #define DYN_ON_SPEED 33 /*20 knots*/
100 #define NZ generic_.n_cg_body_v[2]
106 extern COCKPIT cockpit_;
109 SCALAR interp(SCALAR *y_table, SCALAR *x_table, int Ntable, SCALAR x)
116 /* if x is outside the table, return value at x[0] or x[Ntable-1]*/
120 /* printf("x smaller than x_table[0]: %g %g\n",x,x_table[0]); */
122 else if(x >= x_table[Ntable-1])
125 /* printf("x larger than x_table[N]: %g %g %d\n",x,x_table[NCL-1],Ntable-1); */
127 else /*x is within the table, interpolate linearly to find y value*/
130 while(x_table[i] <= x) {i++;}
131 slope=(y_table[i]-y_table[i-1])/(x_table[i]-x_table[i-1]);
132 /* 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); */
133 y=slope*(x-x_table[i-1]) +y_table[i-1];
141 fprintf(out,"%g,%g,%g,%g,%g,%g,%g,%g,%g,",Long_control,Lat_control,Rudder_pedal,Aft_trim,Fwd_trim,V_rel_wind,Dynamic_pressure,P_body,R_body);
142 fprintf(out,"%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,",Alpha,Cos_alpha,Sin_alpha,Alpha_dot,Q_body,Theta_dot,Sin_theta,Cos_theta,Beta,Cos_beta,Sin_beta);
143 fprintf(out,"%g,%g,%g,%g,%g,%g,%g,%g\n",Sin_phi,Cos_phi,F_X_aero,F_Y_aero,F_Z_aero,M_l_aero,M_m_aero,M_n_aero);
147 void aero( SCALAR dt, int Initialize ) {
151 static SCALAR trim_inc = 0.0002;
155 SCALAR elevator, aileron, rudder;
157 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};
158 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};
163 /*Note that CLo,Cdo,Cmo will likely change with flap setting so
164 they may not be declared static in the future */
167 static SCALAR CLadot=1.7;
168 static SCALAR CLq=3.9;
169 static SCALAR CLde=0.43;
173 static SCALAR Cdo=0.031;
174 static SCALAR Cda=0.13; /*Not used*/
175 static SCALAR Cdde=0.06;
177 static SCALAR Cma=-0.89;
178 static SCALAR Cmadot=-5.2;
179 static SCALAR Cmq=-12.4;
180 static SCALAR Cmo=-0.062;
181 static SCALAR Cmde=-1.28;
183 static SCALAR Clbeta=-0.089;
184 static SCALAR Clp=-0.47;
185 static SCALAR Clr=0.096;
186 static SCALAR Clda=0.178;
187 static SCALAR Cldr=0.0147;
189 static SCALAR Cnbeta=0.065;
190 static SCALAR Cnp=-0.03;
191 static SCALAR Cnr=-0.099;
192 static SCALAR Cnda=-0.053;
193 static SCALAR Cndr=-0.0657;
195 static SCALAR Cybeta=-0.31;
196 static SCALAR Cyp=-0.037;
197 static SCALAR Cyr=0.21;
198 static SCALAR Cyda=0.0;
199 static SCALAR Cydr=0.187;
201 /*nondimensionalization quantities*/
202 /*units here are ft and lbs */
203 static SCALAR cbar=4.9; /*mean aero chord ft*/
204 static SCALAR b=35.8; /*wing span ft */
205 static SCALAR Sw=174; /*wing planform surface area ft^2*/
206 static SCALAR rPiARe=0.054; /*reciprocal of Pi*AR*e*/
210 SCALAR CLwbh,CL,cm,cd,cn,cy,croll,cbar_2V,b_2V,qS,qScbar,qSb,ps,rs;
212 SCALAR F_X_wind,F_Y_wind,F_Z_wind,W_X,W_Y,W_Z;
222 out=fopen("flight.csv","w");
223 /* Initialize aero coefficients */
233 Cl > 0 => Right wing down
236 elevator > 0 => AND -- aircraft nose down
237 aileron > 0 => right wing up
241 if(Aft_trim) long_trim = long_trim - trim_inc;
242 if(Fwd_trim) long_trim = long_trim + trim_inc;
244 /*scale pct control to degrees deflection*/
245 if ((Long_control+long_trim) <= 0)
246 elevator=(Long_control+long_trim)*-28*DEG_TO_RAD;
248 elevator=(Long_control+long_trim)*23*DEG_TO_RAD;
250 aileron = Lat_control*17.5*DEG_TO_RAD;
251 rudder = Rudder_pedal*16*DEG_TO_RAD;
257 /*check control surface travel limits*/
258 /* if((elevator+long_trim) > 23)
260 else if((elevator+long_trim) < -28)
265 The aileron travel limits are 20 deg. TEU and 15 deg TED
266 but since we don't distinguish between left and right we'll
267 use the average here (17.5 deg)
269 /* if(fabs(aileron) > 17.5)
271 if(fabs(rudder) > 16)
274 /*calculate rate derivative nondimensionalization (is that a word?) factors */
275 /*hack to avoid divide by zero*/
276 /*the dynamic terms might be negligible at low ground speeds anyway*/
278 if(V_rel_wind > DYN_ON_SPEED)
280 cbar_2V=cbar/(2*V_rel_wind);
281 b_2V=b/(2*V_rel_wind);
289 /*calcuate the qS nondimensionalization factors*/
291 qS=Dynamic_pressure*Sw;
295 /*transform the aircraft rotation rates*/
296 ps=-P_body*Cos_alpha + R_body*Sin_alpha;
297 rs=-P_body*Sin_alpha + R_body*Cos_alpha;
300 /* sum coefficients */
301 CLwbh = interp(CLtable,alpha_ind,NCL,Alpha);
302 CL = CLo + CLwbh + (CLadot*Alpha_dot + CLq*Theta_dot)*cbar_2V + CLde*elevator;
303 cd = Cdo + rPiARe*CL*CL + Cdde*elevator;
304 cy = Cybeta*Beta + (Cyp*ps + Cyr*rs)*b_2V + Cyda*aileron + Cydr*rudder;
306 cm = Cmo + Cma*Alpha + (Cmq*Theta_dot + Cmadot*Alpha_dot)*cbar_2V + Cmde*(elevator+long_trim);
307 cn = Cnbeta*Beta + (Cnp*ps + Cnr*rs)*b_2V + Cnda*aileron + Cndr*rudder;
308 croll=Clbeta*Beta + (Clp*ps + Clr*rs)*b_2V + Clda*aileron + Cldr*rudder;
310 /*calculate wind axes forces*/
315 /*calculate moments and body axis forces */
317 /*find body-axis components of weight*/
318 /*with earth axis to body axis transform */
320 W_Y=W*Sin_phi*Cos_theta;
321 W_Z=W*Cos_phi*Cos_theta;
323 /* requires ugly wind-axes to body-axes transform */
324 F_X_aero = W_X + F_X_wind*Cos_alpha*Cos_beta - F_Y_wind*Cos_alpha*Sin_beta - F_Z_wind*Sin_alpha;
325 F_Y_aero = W_Y + F_X_wind*Sin_beta + F_Z_wind*Cos_beta;
326 F_Z_aero = W_Z*NZ + F_X_wind*Sin_alpha*Cos_beta - F_Y_wind*Sin_alpha*Sin_beta + F_Z_wind*Cos_alpha;
328 /*no axes transform here */
329 M_l_aero = I_xx*croll*qSb;
330 M_m_aero = I_yy*cm*qScbar;
331 M_n_aero = I_zz*cn*qSb;