1 /***************************************************************************
5 ----------------------------------------------------------------------------
7 FUNCTION: aerodynamics model based on stability derivatives
8 + tweaks for nonlinear aero
10 ----------------------------------------------------------------------------
12 MODULE STATUS: developmental
14 ----------------------------------------------------------------------------
16 GENEALOGY: based on aero model from crrcsim code (Drela's aero model)
18 ----------------------------------------------------------------------------
22 CODED BY: Michael Selig
24 MAINTAINED BY: Michael Selig
26 ----------------------------------------------------------------------------
31 7/25/03 LaRCsim debugging purposes
33 ----------------------------------------------------------------------------
37 ----------------------------------------------------------------------------
41 ----------------------------------------------------------------------------
45 ----------------------------------------------------------------------------
49 ----------------------------------------------------------------------------
53 --------------------------------------------------------------------------*/
57 #include "ls_generic.h"
58 #include "ls_cockpit.h"
59 #include "ls_constants.h"
61 #include "basic_aero.h"
68 #define NZ generic_.n_cg_body_v[2]
74 extern COCKPIT cockpit_;
77 void basic_aero(SCALAR dt, int Initialize)
78 // Calculate forces and moments for the current time step. If Initialize is
79 // zero, then re-initialize coefficients by reading in the coefficient file.
82 static SCALAR elevator_drela, aileron_drela, rudder_drela;
136 SCALAR CD_stall = 0.05;
145 SCALAR Phat, Qhat, Rhat; // dimensionless rotation rates
146 SCALAR CL_left, CL_cent, CL_right; // CL values across the span
147 SCALAR dCL_left,dCL_cent,dCL_right; // stall-induced CL changes
148 SCALAR dCD_left,dCD_cent,dCD_right; // stall-induced CD changes
149 SCALAR dCl,dCn,dCl_stall,dCn_stall; // stall-induced roll,yaw moments
150 SCALAR dCm_stall; // Stall induced pitching moment.
151 SCALAR CL_wing, CD_all, CD_scaled, Cl_w;
152 SCALAR Cl_r_mod,Cn_p_mod;
153 SCALAR CL_drela,CD_drela,Cx_drela,Cy_drela,Cz_drela,Cl_drela,Cm_drela,Cn_drela;
155 SCALAR G_11,G_12,G_13;
156 SCALAR G_21,G_22,G_23;
157 SCALAR G_31,G_32,G_33;
158 SCALAR U_body_X,U_body_Y,U_body_Z;
159 SCALAR V_body_X,V_body_Y,V_body_Z;
160 SCALAR W_body_X,W_body_Y,W_body_Z;
161 SCALAR P_atmo,Q_atmo,R_atmo;
163 // set the parameters
169 Alpha_0 = 0.349066E-01;
170 Cm_0 = -0.112663E-01;
174 CD_prof = 0.200000E-01;
185 Cn_p = -0.740898E-01;
188 Cn_r = -0.687755E-01;
191 CY_dr = 0.000000E+00;
192 Cl_dr = 0.000000E+00;
193 Cn_dr = 0.000000E+00;
195 Cl_da = -0.307921E-02;
196 Cn_da = 0.527143E-01;
211 // jan's data goes -.5 to .5 while
212 // fgfs data goes +- 1.
213 // so I need to divide by 2 below.
214 elevator = Long_control + Long_trim;
215 aileron = Lat_control;
216 rudder = Rudder_pedal;
218 elevator = elevator * 0.5;
219 aileron = aileron * 0.5;
220 rudder = rudder * 0.5;
224 // printf("%f\n",V_rel_wind);
226 /* compute gradients of Local velocities w.r.t. Body coordinates
227 G_11 = dU_local/dx_body
228 G_12 = dU_local/dy_body etc. */
231 G_11 = U_atmo_X*T_local_to_body_11
232 + U_atmo_Y*T_local_to_body_12
233 + U_atmo_Z*T_local_to_body_13;
234 G_12 = U_atmo_X*T_local_to_body_21
235 + U_atmo_Y*T_local_to_body_22
236 + U_atmo_Z*T_local_to_body_23;
237 G_13 = U_atmo_X*T_local_to_body_31
238 + U_atmo_Y*T_local_to_body_32
239 + U_atmo_Z*T_local_to_body_33;
241 G_21 = V_atmo_X*T_local_to_body_11
242 + V_atmo_Y*T_local_to_body_12
243 + V_atmo_Z*T_local_to_body_13;
244 G_22 = V_atmo_X*T_local_to_body_21
245 + V_atmo_Y*T_local_to_body_22
246 + V_atmo_Z*T_local_to_body_23;
247 G_23 = V_atmo_X*T_local_to_body_31
248 + V_atmo_Y*T_local_to_body_32
249 + V_atmo_Z*T_local_to_body_33;
251 G_31 = W_atmo_X*T_local_to_body_11
252 + W_atmo_Y*T_local_to_body_12
253 + W_atmo_Z*T_local_to_body_13;
254 G_32 = W_atmo_X*T_local_to_body_21
255 + W_atmo_Y*T_local_to_body_22
256 + W_atmo_Z*T_local_to_body_23;
257 G_33 = W_atmo_X*T_local_to_body_31
258 + W_atmo_Y*T_local_to_body_32
259 + W_atmo_Z*T_local_to_body_33;
262 //printf("%f %f %f %f\n",W_atmo_X,W_atmo_Y,G_31,G_32);
264 /* now compute gradients of Body velocities w.r.t. Body coordinates */
265 /* U_body_x = dU_body/dx_body etc. */
268 U_body_X = T_local_to_body_11*G_11
269 + T_local_to_body_12*G_21
270 + T_local_to_body_13*G_31;
271 U_body_Y = T_local_to_body_11*G_12
272 + T_local_to_body_12*G_22
273 + T_local_to_body_13*G_32;
274 U_body_Z = T_local_to_body_11*G_13
275 + T_local_to_body_12*G_23
276 + T_local_to_body_13*G_33;
278 V_body_X = T_local_to_body_21*G_11
279 + T_local_to_body_22*G_21
280 + T_local_to_body_23*G_31;
281 V_body_Y = T_local_to_body_21*G_12
282 + T_local_to_body_22*G_22
283 + T_local_to_body_23*G_32;
284 V_body_Z = T_local_to_body_21*G_13
285 + T_local_to_body_22*G_23
286 + T_local_to_body_23*G_33;
288 W_body_X = T_local_to_body_31*G_11
289 + T_local_to_body_32*G_21
290 + T_local_to_body_33*G_31;
291 W_body_Y = T_local_to_body_31*G_12
292 + T_local_to_body_32*G_22
293 + T_local_to_body_33*G_32;
294 W_body_Z = T_local_to_body_31*G_13
295 + T_local_to_body_32*G_23
296 + T_local_to_body_33*G_33;
299 /* set rotation rates of airmass motion */
316 // printf("%f %f %f\n",P_atmo,Q_atmo,R_atmo);
319 /* set net effective dimensionless rotation rates */
320 Phat = (P_body - P_atmo) * B_ref / (2.0*V_rel_wind);
321 Qhat = (Q_body - Q_atmo) * C_ref / (2.0*V_rel_wind);
322 Rhat = (R_body - R_atmo) * B_ref / (2.0*V_rel_wind);
331 // printf("Phat: %f Qhat: %f Rhat: %f\n",Phat,Qhat,Rhat);
333 /* compute local CL at three spanwise locations */
334 CL_left = CL_0 + CL_a*(Std_Alpha - Alpha_0 - Phat*eta_loc);
335 CL_cent = CL_0 + CL_a*(Std_Alpha - Alpha_0 );
336 CL_right = CL_0 + CL_a*(Std_Alpha - Alpha_0 + Phat*eta_loc);
338 // printf("CL_left: %f CL_cent: %f CL_right: %f\n",CL_left,CL_cent,CL_right);
340 /* set CL-limit changes */
346 if (CL_left > CL_max)
348 dCL_left = CL_max-CL_left -CL_drop;
352 if (CL_cent > CL_max)
354 dCL_cent = CL_max-CL_cent -CL_drop;
358 if (CL_right > CL_max)
360 dCL_right = CL_max-CL_right -CL_drop;
364 if (CL_left < CL_min)
366 dCL_left = CL_min-CL_left -CL_drop;
370 if (CL_cent < CL_min)
372 dCL_cent = CL_min-CL_cent -CL_drop;
376 if (CL_right < CL_min)
378 dCL_right = CL_min-CL_right -CL_drop;
382 /* set average wing CL */
383 CL_wing = CL_0 + CL_a*(Std_Alpha-Alpha_0)
384 + 0.25*dCL_left + 0.5*dCL_cent + 0.25*dCL_right;
386 // printf("CL_wing: %f\n",CL_wing);
389 /* correct profile CD for CL dependence and aileron dependence */
391 + CD_CLsq*(CL_wing-CL_CD0)*(CL_wing-CL_CD0)
392 + CD_AIsq*aileron*aileron
393 + CD_ELsq*elevator*elevator;
395 // printf("CD_all:%lf\n",CD_all);
397 /* scale profile CD with Reynolds number via simple power law */
398 if (V_rel_wind > 0.1)
400 CD_scaled = CD_all*pow(((double)V_rel_wind/(double)U_ref),Uexp_CD);
407 // printf("CD_scaled:%lf\n",CD_scaled);
411 /* Scale lateral cross-coupling derivatives with wing CL */
412 Cl_r_mod = Cl_r*CL_wing/CL_0;
413 Cn_p_mod = Cn_p*CL_wing/CL_0;
415 // printf("Cl_r_mod: %f Cn_p_mod: %f\n",Cl_r_mod,Cn_p_mod);
417 /* total average CD with induced and stall contributions */
418 dCD_left = CD_stall*dCL_left *dCL_left ;
419 dCD_cent = CD_stall*dCL_cent *dCL_cent ;
420 dCD_right = CD_stall*dCL_right*dCL_right;
422 /* total CL, with pitch rate and elevator contributions */
423 CL_drela = (CL_wing + CL_q*Qhat + CL_de*elevator)*Cos_alpha;
425 // printf("CL:%f\n",CL);
427 /* assymetric stall will cause roll and yaw moments */
428 dCl = 0.45*-1*(dCL_right-dCL_left)*0.5*eta_loc;
429 dCn = 0.45*(dCD_right-dCD_left)*0.5*eta_loc;
430 dCm_stall = (0.25*dCL_left + 0.5*dCL_cent + 0.25*dCL_right)*CG_arm;
432 // printf("dCl: %f dCn:%f\n",dCl,dCn);
434 /* stall-caused moments in body axes */
435 dCl_stall = dCl*Cos_alpha - dCn*Sin_alpha;
436 dCn_stall = dCl*Sin_alpha + dCn*Cos_alpha;
438 /* total CD, with induced and stall contributions */
440 Cl_w = Cl_b*Std_Beta + Cl_p*Phat + Cl_r_mod*Rhat
441 + dCl_stall + Cl_da*aileron;
443 + (CL*CL + 32.0*Cl_w*Cl_w)*S_ref/(B_ref*B_ref*3.14159*span_eff)
444 + 0.25*dCD_left + 0.5*dCD_cent + 0.25*dCD_right;
446 //printf("CL: %f CD:%f L/D:%f\n",CL,CD,CL/CD);
448 /* total forces in body axes */
449 Cx_drela = -CD_drela*Cos_alpha + CL_drela*Sin_alpha*Cos_beta*Cos_beta;
450 Cz_drela = -CD_drela*Sin_alpha - CL_drela*Cos_alpha*Cos_beta*Cos_beta;
451 Cy_drela = CY_b*Std_Beta + CY_p*Phat + CY_r*Rhat + CY_dr*rudder;
453 /* total moments in body axes */
454 Cl_drela = Cl_b*Std_Beta + Cl_p*Phat + Cl_r_mod*Rhat + Cl_dr*rudder
455 + dCl_stall + Cl_da*aileron;
456 Cn_drela = Cn_b*Std_Beta + Cn_p_mod*Phat + Cn_r*Rhat + Cn_dr*rudder
457 + dCn_stall + Cn_da*aileron;
458 Cm_drela = Cm_0 + Cm_a*(Std_Alpha-Alpha_0) +dCm_stall
459 + Cm_q*Qhat + Cm_de*elevator;
461 /* set dimensional forces and moments */
462 QS = 0.5*Density*V_rel_wind*V_rel_wind * S_ref;
464 F_X_aero = Cx_drela * QS;
465 F_Y_aero = Cy_drela * QS;
466 F_Z_aero = Cz_drela * QS;
468 M_l_aero = Cl_drela * QS * B_ref;
469 M_m_aero = Cm_drela * QS * C_ref;
470 M_n_aero = Cn_drela * QS * B_ref;
471 // printf("%f %f %f %f %f %f\n",F_X_aero,F_Y_aero,F_Z_aero,M_l_aero,M_m_aero,M_n_aero);