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;
156 SCALAR G_11,G_12,G_13;
157 SCALAR G_21,G_22,G_23;
158 SCALAR G_31,G_32,G_33;
159 SCALAR U_body_X,U_body_Y,U_body_Z;
160 SCALAR V_body_X,V_body_Y,V_body_Z;
161 SCALAR W_body_X,W_body_Y,W_body_Z;
163 SCALAR P_atmo,Q_atmo,R_atmo;
165 // set the parameters
171 Alpha_0 = 0.349066E-01;
172 Cm_0 = -0.112663E-01;
176 CD_prof = 0.200000E-01;
187 Cn_p = -0.740898E-01;
190 Cn_r = -0.687755E-01;
193 CY_dr = 0.000000E+00;
194 Cl_dr = 0.000000E+00;
195 Cn_dr = 0.000000E+00;
197 Cl_da = -0.307921E-02;
198 Cn_da = 0.527143E-01;
213 // jan's data goes -.5 to .5 while
214 // fgfs data goes +- 1.
215 // so I need to divide by 2 below.
216 elevator = Long_control + Long_trim;
217 aileron = Lat_control;
218 rudder = Rudder_pedal;
220 elevator = elevator * 0.5;
221 aileron = aileron * 0.5;
222 rudder = rudder * 0.5;
226 // printf("%f\n",V_rel_wind);
228 /* compute gradients of Local velocities w.r.t. Body coordinates
229 G_11 = dU_local/dx_body
230 G_12 = dU_local/dy_body etc. */
233 G_11 = U_atmo_X*T_local_to_body_11
234 + U_atmo_Y*T_local_to_body_12
235 + U_atmo_Z*T_local_to_body_13;
236 G_12 = U_atmo_X*T_local_to_body_21
237 + U_atmo_Y*T_local_to_body_22
238 + U_atmo_Z*T_local_to_body_23;
239 G_13 = U_atmo_X*T_local_to_body_31
240 + U_atmo_Y*T_local_to_body_32
241 + U_atmo_Z*T_local_to_body_33;
243 G_21 = V_atmo_X*T_local_to_body_11
244 + V_atmo_Y*T_local_to_body_12
245 + V_atmo_Z*T_local_to_body_13;
246 G_22 = V_atmo_X*T_local_to_body_21
247 + V_atmo_Y*T_local_to_body_22
248 + V_atmo_Z*T_local_to_body_23;
249 G_23 = V_atmo_X*T_local_to_body_31
250 + V_atmo_Y*T_local_to_body_32
251 + V_atmo_Z*T_local_to_body_33;
253 G_31 = W_atmo_X*T_local_to_body_11
254 + W_atmo_Y*T_local_to_body_12
255 + W_atmo_Z*T_local_to_body_13;
256 G_32 = W_atmo_X*T_local_to_body_21
257 + W_atmo_Y*T_local_to_body_22
258 + W_atmo_Z*T_local_to_body_23;
259 G_33 = W_atmo_X*T_local_to_body_31
260 + W_atmo_Y*T_local_to_body_32
261 + W_atmo_Z*T_local_to_body_33;
264 //printf("%f %f %f %f\n",W_atmo_X,W_atmo_Y,G_31,G_32);
266 /* now compute gradients of Body velocities w.r.t. Body coordinates */
267 /* U_body_x = dU_body/dx_body etc. */
270 U_body_X = T_local_to_body_11*G_11
271 + T_local_to_body_12*G_21
272 + T_local_to_body_13*G_31;
273 U_body_Y = T_local_to_body_11*G_12
274 + T_local_to_body_12*G_22
275 + T_local_to_body_13*G_32;
276 U_body_Z = T_local_to_body_11*G_13
277 + T_local_to_body_12*G_23
278 + T_local_to_body_13*G_33;
280 V_body_X = T_local_to_body_21*G_11
281 + T_local_to_body_22*G_21
282 + T_local_to_body_23*G_31;
283 V_body_Y = T_local_to_body_21*G_12
284 + T_local_to_body_22*G_22
285 + T_local_to_body_23*G_32;
286 V_body_Z = T_local_to_body_21*G_13
287 + T_local_to_body_22*G_23
288 + T_local_to_body_23*G_33;
290 W_body_X = T_local_to_body_31*G_11
291 + T_local_to_body_32*G_21
292 + T_local_to_body_33*G_31;
293 W_body_Y = T_local_to_body_31*G_12
294 + T_local_to_body_32*G_22
295 + T_local_to_body_33*G_32;
296 W_body_Z = T_local_to_body_31*G_13
297 + T_local_to_body_32*G_23
298 + T_local_to_body_33*G_33;
301 /* set rotation rates of airmass motion */
318 // printf("%f %f %f\n",P_atmo,Q_atmo,R_atmo);
321 /* set net effective dimensionless rotation rates */
322 Phat = (P_body - P_atmo) * B_ref / (2.0*V_rel_wind);
323 Qhat = (Q_body - Q_atmo) * C_ref / (2.0*V_rel_wind);
324 Rhat = (R_body - R_atmo) * B_ref / (2.0*V_rel_wind);
333 // printf("Phat: %f Qhat: %f Rhat: %f\n",Phat,Qhat,Rhat);
335 /* compute local CL at three spanwise locations */
336 CL_left = CL_0 + CL_a*(Std_Alpha - Alpha_0 - Phat*eta_loc);
337 CL_cent = CL_0 + CL_a*(Std_Alpha - Alpha_0 );
338 CL_right = CL_0 + CL_a*(Std_Alpha - Alpha_0 + Phat*eta_loc);
340 // printf("CL_left: %f CL_cent: %f CL_right: %f\n",CL_left,CL_cent,CL_right);
342 /* set CL-limit changes */
348 if (CL_left > CL_max)
350 dCL_left = CL_max-CL_left -CL_drop;
354 if (CL_cent > CL_max)
356 dCL_cent = CL_max-CL_cent -CL_drop;
360 if (CL_right > CL_max)
362 dCL_right = CL_max-CL_right -CL_drop;
366 if (CL_left < CL_min)
368 dCL_left = CL_min-CL_left -CL_drop;
372 if (CL_cent < CL_min)
374 dCL_cent = CL_min-CL_cent -CL_drop;
378 if (CL_right < CL_min)
380 dCL_right = CL_min-CL_right -CL_drop;
384 /* set average wing CL */
385 CL_wing = CL_0 + CL_a*(Std_Alpha-Alpha_0)
386 + 0.25*dCL_left + 0.5*dCL_cent + 0.25*dCL_right;
388 // printf("CL_wing: %f\n",CL_wing);
391 /* correct profile CD for CL dependence and aileron dependence */
393 + CD_CLsq*(CL_wing-CL_CD0)*(CL_wing-CL_CD0)
394 + CD_AIsq*aileron*aileron
395 + CD_ELsq*elevator*elevator;
397 // printf("CD_all:%lf\n",CD_all);
399 /* scale profile CD with Reynolds number via simple power law */
400 if (V_rel_wind > 0.1)
402 CD_scaled = CD_all*pow(((double)V_rel_wind/(double)U_ref),Uexp_CD);
409 // printf("CD_scaled:%lf\n",CD_scaled);
413 /* Scale lateral cross-coupling derivatives with wing CL */
414 Cl_r_mod = Cl_r*CL_wing/CL_0;
415 Cn_p_mod = Cn_p*CL_wing/CL_0;
417 // printf("Cl_r_mod: %f Cn_p_mod: %f\n",Cl_r_mod,Cn_p_mod);
419 /* total average CD with induced and stall contributions */
420 dCD_left = CD_stall*dCL_left *dCL_left ;
421 dCD_cent = CD_stall*dCL_cent *dCL_cent ;
422 dCD_right = CD_stall*dCL_right*dCL_right;
424 /* total CL, with pitch rate and elevator contributions */
425 CL_drela = (CL_wing + CL_q*Qhat + CL_de*elevator)*Cos_alpha;
427 // printf("CL:%f\n",CL);
429 /* assymetric stall will cause roll and yaw moments */
430 dCl = 0.45*-1*(dCL_right-dCL_left)*0.5*eta_loc;
431 dCn = 0.45*(dCD_right-dCD_left)*0.5*eta_loc;
432 dCm_stall = (0.25*dCL_left + 0.5*dCL_cent + 0.25*dCL_right)*CG_arm;
434 // printf("dCl: %f dCn:%f\n",dCl,dCn);
436 /* stall-caused moments in body axes */
437 dCl_stall = dCl*Cos_alpha - dCn*Sin_alpha;
438 dCn_stall = dCl*Sin_alpha + dCn*Cos_alpha;
440 /* total CD, with induced and stall contributions */
442 Cl_w = Cl_b*Std_Beta + Cl_p*Phat + Cl_r_mod*Rhat
443 + dCl_stall + Cl_da*aileron;
445 + (CL*CL + 32.0*Cl_w*Cl_w)*S_ref/(B_ref*B_ref*3.14159*span_eff)
446 + 0.25*dCD_left + 0.5*dCD_cent + 0.25*dCD_right;
448 //printf("CL: %f CD:%f L/D:%f\n",CL,CD,CL/CD);
450 /* total forces in body axes */
451 Cx_drela = -CD_drela*Cos_alpha + CL_drela*Sin_alpha*Cos_beta*Cos_beta;
452 Cz_drela = -CD_drela*Sin_alpha - CL_drela*Cos_alpha*Cos_beta*Cos_beta;
453 Cy_drela = CY_b*Std_Beta + CY_p*Phat + CY_r*Rhat + CY_dr*rudder;
455 /* total moments in body axes */
456 Cl_drela = Cl_b*Std_Beta + Cl_p*Phat + Cl_r_mod*Rhat + Cl_dr*rudder
457 + dCl_stall + Cl_da*aileron;
458 Cn_drela = Cn_b*Std_Beta + Cn_p_mod*Phat + Cn_r*Rhat + Cn_dr*rudder
459 + dCn_stall + Cn_da*aileron;
460 Cm_drela = Cm_0 + Cm_a*(Std_Alpha-Alpha_0) +dCm_stall
461 + Cm_q*Qhat + Cm_de*elevator;
463 /* set dimensional forces and moments */
464 QS = 0.5*Density*V_rel_wind*V_rel_wind * S_ref;
466 F_X_aero = Cx_drela * QS;
467 F_Y_aero = Cy_drela * QS;
468 F_Z_aero = Cz_drela * QS;
470 M_l_aero = Cl_drela * QS * B_ref;
471 M_m_aero = Cm_drela * QS * C_ref;
472 M_n_aero = Cn_drela * QS * B_ref;
473 // 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);