c/////////////////////////////////////////////////////////////////////// c c Virginia Tech Department of Aerospace and Ocean Engineeirng c c Program LDstab estimates Lateral/Directional Stability and Control c derivatives for transonic tranport class aircraft. It consists of c a main program stabin that reads in the data, and a subroutine stab c that does the calculations. c c The program was written by Joel Grasmeyer, and is described in c VPI-AOE-254, which explains the input data required. The sign c conventions are also described there. c c program stabin c c This program reads the stability input file, and calls the stab c subroutine to calculate the maximum available yawing moment c coefficient. c c See stab.f for variable definitions c c Created by: Joel Grasmeyer c Last Modified: 10/28/97 c c Mods: c - Modified for use in Aircraft design class by W.H. Mason, 3/30/2000 c - Changed header, added pause at end, added an output file, c L.T.Leifsson March 2004 c - Added Cnbeta contribution from wing using DATCOM method, c this required adding two input variables: xbar and cbar, c by L.T. Leifsson, March 2004. c c/////////////////////////////////////////////////////////////////////// program stabin implicit none character*72 infile, outfile, header, title integer i, write_flag, unit_in, unit_outfile, new, nef real dihedral_wing, z_wing, dia_fuse, sref, hspan_vtail, & depth_fuse_vtail, c_vtail_root, c_vtail_tip, mach_eo, l_vtail, & sweep_wing_1_2, hspan_wing, cl, z_vtail, & length_fuse, rho_eo, a_eo, mu_eo, cn_avail, dia_nacelle, & ar_vtail_eff, dr_max, sh, sweep_vtail_1_4, xbar, cbar c Get input and output filenames from command line c write(6,10) write(6,3000) write(6,3001) write(6,3002) write(6,3003) write(6,3004) write(6,3005) read(5,20) infile write(6,3006) read(5,20) outfile c Read input that is common to both modes unit_in = 10 open(unit_in,file=infile) read(unit_in,"(a72)") header read(unit_in,"(a72)") title read(unit_in,*) write_flag read(unit_in,*) dihedral_wing read(unit_in,*) z_wing read(unit_in,*) dia_fuse read(unit_in,*) sref read(unit_in,*) hspan_vtail read(unit_in,*) depth_fuse_vtail read(unit_in,*) c_vtail_root read(unit_in,*) c_vtail_tip read(unit_in,*) mach_eo read(unit_in,*) sweep_vtail_1_4 read(unit_in,*) sweep_wing_1_2 read(unit_in,*) hspan_wing read(unit_in,*) cl read(unit_in,*) z_vtail read(unit_in,*) l_vtail read(unit_in,*) length_fuse read(unit_in,*) new read(unit_in,*) nef read(unit_in,*) dia_nacelle read(unit_in,*) sh read(unit_in,*) rho_eo read(unit_in,*) a_eo read(unit_in,*) mu_eo read(unit_in,*) dr_max read(unit_in,*) xbar read(unit_in,*) cbar c Close input file close(unit_in) c Open output file and write header unit_outfile = 20 open(unit_outfile,file=outfile) write(unit_outfile,3000) write(unit_outfile,3001) write(unit_outfile,3002) write(unit_outfile,3003) write(unit_outfile,3004) c Echo input to output file c Call stab subroutine call stab(outfile,title,write_flag,dihedral_wing,z_wing, & dia_fuse,sref,hspan_vtail,depth_fuse_vtail,c_vtail_root, & c_vtail_tip,mach_eo,sweep_vtail_1_4,sweep_wing_1_2,hspan_wing, & cl,z_vtail,l_vtail,length_fuse,new,nef,dia_nacelle, & sh,rho_eo,a_eo,mu_eo,dr_max,xbar,cbar,ar_vtail_eff,cn_avail) c Pause before closing program PAUSE 'Press ENTER to close program' c Write statements and formats 10 format(/2x, 'Lateral-Directional Stability and Control Derivative', 1 ' Estimation'/2x, 1 '-intended for transonic tranports, i.e., B747 -'//2x, 1 'Joel Grasmeyer'/2x, 2 'Aerospace and Ocean Engineering, Virginia Tech'/2x, 3 'Blacksburg, VA 24061'//2x, 4 'Aircraft Design Class Software'//2x, 5 'Enter the name of input data file') 20 format(a20) 3000 format(' **************************************************'/ 1 ' VT Aerospace Aircraft Design Software Series'/ 2 ' **************************************************'/) 3001 format(' Program LDstab'/ 1 ' The program estimates Lateral/Directional Stability'/ 2 ' and Control derivatives for transonic tranport class'/ 3 ' aircraft.'/) 3002 format(' Written by Joel Grasmeyer.'/ 1 ' This program modified by'/ 2 ' - W.H. Mason, March 2000,'/ 3 ' - L.T. Leifsson, March 2004.'/) 3003 format(' The Department of Aerospace and Ocean Engineering,'/ 1 ' Virginia Tech, Blacksburg, VA 24061'/ 2 ' http://www.aoe.vt.edu, email: whmason@vt.edu'/) 3004 format(' **************************************************'/) 3005 format(' Enter input filename:') 3006 format(' Enter output filename:') end ! End program stabin c/////////////////////////////////////////////////////////////////////// c c subroutine stab c c This subroutine calculates the maximum available yawing moment c coefficient. Note that right rudder deflection is defined as c positive, and right aileron up, left aileron down is defined as c positive. Both of these control deflections generate positive c moments about their respective axes. This is the convention used c by Roskam. c c Inputs c c outfile output filename c title title of aircraft configuration c write_flag write flag (0 = no output file, 1 = output file written) c dihedral_wing wing dihedral angle (deg) c z_wing distance from body centerline to quarter-chord point of c exposed wing root chord, positive for the quarter-chord c point below the body centerline (ft) c dia_fuse fuselage diameter (ft) c sref wing reference area (ft^2) c hspan_vtail vertical tail span (ft) c depth_fuse_vtail fuselage depth at the fuselage station of the c quarter-chord of the vertical tail (ft) c c_vtail_root root chord of vertical tail c c_vtail_tip tip chord of vertical tail c mach_eo mach number c sweep_vtail_1_4 vertical tail quarter-chord sweep angle (deg) c sweep_wing_1_2 average wing half-chord sweep angle (deg) c hspan_wing wing half-span (ft) c cl lift coefficient c z_vtail vertical distance from CG to AC of vertical tail (ft) c l_vtail horizontal distance from CG to AC of vertical tail (ft) c length_fuse fuselage length (ft) c new number of engines on the wing c nef number of engines on the fuselage c dia_nacelle nacelle diameter (ft) c rho_eo density at engine-out flight condition (slug/ft^3) c a_eo speed of sound at engine-out flight condition (ft/s) c mu_eo viscosity at engine-out flight condition (slug/(ft-s)) c dr_max maximum allowable steady-state rudder deflection (deg) c xbar longitudinal distance (positive rearward) from coordinate c origin (usually the center of gravity) to the wing aero- c dynamic center (ft), L.T. Leifsson March 2004. c cbar mean aerodynamic chord (ft), L.T. Leifsson March 2004. c c Outputs c c ar_vtail_eff effective aspect ratio of vertical tail c cn_avail maximum available yawing moment coefficient c c Internal Variables c c alpha angle of attack (rad) c alpha_d section lift effectiveness c alpha_d_cl section flap effectiveness (from Figure 10.2) c ar wing aspect ratio c ar_vtail actual aspect ratio of vertical tail c beta sideslip angle, positive from the right (rad) c beta_m square root of (1 - mach_eo)**2 c cf_c ratio of flap chord to wing or tail chord c cl_alpha lift-curve slope of entire aircraft (rad^-1) c cl_alpha_2d 2-dimensional lift-curve slope at MAC (rad^-1) c cl_alpha_vtail original lift-curve slope of vertical tail (rad^-1) c cl_alpha_vtail_eff effective lift-curve slope of vertical tail (rad^-1) c cl_beta variation of rolling moment coefficient with sideslip angle c cl_beta_htail horizontal tail contribution to cl_beta c cl_beta_vtail vertical tail contribution to cl_beta c cl_beta_wingbody wing-body contribution to cl_beta c cl_d rolling effectiveness of partial-chord controls c cl_da variation of rolling moment coefficient with aileron deflection c cl_dr variation of rolling moment coefficient with rudder deflection c cn_beta variation of yawing moment coefficient with sideslip angle c cn_beta_fuse fuselage contribution to cn_beta c cn_beta_vtail vertical tail contribution to cn_beta c cn_beta_wing wing contribution to cn_beta c cn_da variation of yawing moment coefficient with aileron deflection c cn_dr variation of yawing moment coefficient with rudder deflection c cy_beta variation of side force coefficient with sideslip angle c cy_beta_fuse fuselage contribution to cy_beta c cy_beta_vtail vertical tail contribution to cy_beta c cy_beta_wing wing contribution to cy_beta c cy_da variation of side force coefficient with aileron deflection c cy_dr variation of side force coefficient with rudder deflection c da aileron deflection, positive for right aileron up, left c aileron down (rad) c dr rudder deflection, positive for right deflection (rad) c eqn_7_5 factor estimated by Equation 7.5 c eqn_7_10 body-induced effect on wing height (from Equation 7.10) c eqn_7_11 d in Equation 7.10 (estimated from Equation 7.11) c eqn_7_12 another body-induced effect on wing height (from Equation 7.12) c eqn_11_2 rolling effectiveness of two full-chord ailerons (Eqn. 11.2) c f_cy_beta correction factor for cy_beta c f_cl_beta correction factor for cl_beta c f_cn_beta correction factor for cn_beta c f_cl_da correction factor for cl_da c f_cn_da correction factor for cn_da c f_cy_dr correction factor for cy_dr c f_cl_dr correction factor for cl_dr c f_cn_dr correction factor for cn_dr c fig_7_1 wing-body interference factor from Figure 7.1 c fig_7_3 empirical factor from Figure 7.3 c fig_7_5 ratio of the aspect ratio of the vertical panel in the presence c of the body to that of the isolated panel (from Figure 7.5) c fig_7_6 ratio of the vertical panel aspect ratio in the presence of c the horizontal tail and body to that of the panel in the c presence of the body alone (from Figure 7.6) c fig_7_7 factor accounting for relative size of horizontal and vertical c tails (from Figure 7.7) c fig_7_11 wing sweep contribution to cl_beta_wingbody (from Figure 7.11) c fig_7_12 compressibility correction to wing sweep (from Figure 7.12) c fig_7_13 fuselage correction factor (from Figure 7.13) c fig_7_14 aspect ratio contribution to cl_beta_wingbody (from Figure 7.14) c fig_7_15 dihedral effect on cl_beta (from Figure 7.15) c fig_7_16 compressibility correction to dihedral effect (from Figure 7.16) c fig_7_19 factor for body and body + wing effects (from Figure 7.19) c fig_7_20 Reynold's number factor for the fuselage (from Figure 7.20) c fig_10_2 flap chord factor (from Figure 10.2) c fig_10_3 span factor for plain flap (from Figure 10.3) c fig_10_5 theoretical lift effectiveness of plain TE flaps (Fig. 10.5) c fig_10_6 empirical correction for plain TE flaps (Fig. 10.6) c fig_10_7 empirical correction for lift effectiveness of plain flaps c at high flap deflections (from Figure 10.7) c fig_11_1 rolling moment effectiveness parameter (from Figure 11.1) c flap_eff_ratio flap effectiveness ratio (from Figure 10.2) c i index c k empirical factor for estimating the variation of yawing c moment coefficient with aileron deflection c kappa ratio of the actual lift-curve slope to 2*pi c phi bank angle, positive to the right (rad) c q dynamic pressure (lb/ft^2) c re_fuse fuselage Reynolds number c sbs body side area (ft^2) c sh area of horizontal tail (ft^2) c sv area of vertical tail (ft^2) c sweep_vtail_1_2 vertical tail half-chord sweep angle (deg) c x temporary variable for curve fits c c Created by: Joel Grasmeyer c Last Modified: 10/28/97 c c I/O changes by W.H. Mason, 3/30/2000 c c/////////////////////////////////////////////////////////////////////// subroutine stab(outfile,title,write_flag,dihedral_wing,z_wing, & dia_fuse,sref,hspan_vtail,depth_fuse_vtail,c_vtail_root, & c_vtail_tip,mach_eo,sweep_vtail_1_4,sweep_wing_1_2,hspan_wing, & cl,z_vtail,l_vtail,length_fuse,new,nef,dia_nacelle, & sh,rho_eo,a_eo,mu_eo,dr_max,xbar,cbar,ar_vtail_eff,cn_avail) implicit none character*72 outfile, title integer i, write_flag, unit_out, unit_outfile, new, nef real pi, dihedral_wing, z_wing, dia_fuse, sref, hspan_vtail, ar, & depth_fuse_vtail, c_vtail_root, c_vtail_tip, mach_eo, sv, sh, & sweep_wing_1_2, hspan_wing, cy_beta, ar_vtail, k, & cy_beta_wing, cy_beta_fuse, cy_beta_vtail, ar_vtail_eff, alpha, & cl_alpha_vtail, beta_m, eqn_7_5, kappa, cl_alpha_vtail_eff, q, & cl_beta, cl_beta_htail, cl_beta_vtail, cl_beta_wingbody, sbs, & cl, fig_7_1, fig_7_5, fig_7_6, fig_7_7, fig_7_11, cn_da, & fig_7_12, fig_7_13, fig_7_14, fig_7_15, fig_7_16, eqn_7_10, da, & eqn_7_11, eqn_7_12, cl_alpha, z_vtail, l_vtail, cn_beta_fuse, & cn_beta_vtail, cn_beta, cn_beta_wing, fig_7_19, fig_7_20, phi, & length_fuse, re_fuse, cl_da, cy_da, fig_11_1, eqn_11_2, cl_d, & alpha_d, fig_10_5, fig_10_6, cl_alpha_2d, cl_dr, dr, & cn_dr, cy_dr, fig_10_2, fig_10_7, alpha_d_cl, flap_eff_ratio, & fig_10_3, cf_c, beta, rho_eo, a_eo, mu_eo, cn_avail, fig_7_3, & dia_nacelle, dr_max, f_cy_beta, f_cl_beta, & f_cn_beta, f_cl_da, f_cn_da, f_cy_dr, f_cl_dr, f_cn_dr, x, & sweep_vtail_1_4, sweep_vtail_1_2, & cnbw1, cnbw2, cnbw3, cnbw4, cnbw5, cnbw6, cn_beta_wingterm, & cbar, xbar pi = acos(-1.) c Convert sweep angles from degrees to radians cwhm THIS IS A VERY DANGEROUS PRACTICE!! DON'T USE AS AN EXAMPLE cwhm of acceptable coding style sweep_wing_1_2 = sweep_wing_1_2*pi/180. sweep_vtail_1_4 = sweep_vtail_1_4*pi/180. c Write input data to screen for confirmation if (write_flag .eq. 1) then unit_out = 6 write(unit_out,90) title write(unit_out,*) write(unit_out,"(a25)") ' Input Data' write(unit_out,*) write(unit_out,101) write_flag, '= write_flag' write(unit_out,100) dihedral_wing, '= dihedral_wing' write(unit_out,100) z_wing, '= z_wing' write(unit_out,100) dia_fuse, '= dia_fuse' write(unit_out,100) sref, '= sref' write(unit_out,100) hspan_vtail, '= hspan_vtail' write(unit_out,100) depth_fuse_vtail, '= depth_fuse_vtail' write(unit_out,100) c_vtail_root, '= c_vtail_root' write(unit_out,100) c_vtail_tip, '= c_vtail_tip' write(unit_out,100) mach_eo, '= mach_eo' write(unit_out,100) sweep_vtail_1_4*180./pi, '= sweep_vtail_1_4' write(unit_out,100) sweep_wing_1_2*180./pi, '= sweep_wing_1_2' write(unit_out,100) hspan_wing, '= hspan_wing' write(unit_out,100) cl, '= cl' write(unit_out,100) z_vtail, '= z_vtail' write(unit_out,100) l_vtail, '= l_vtail' write(unit_out,100) length_fuse, '= length_fuse' write(unit_out,101) new, '= new' write(unit_out,101) nef, '= nef' write(unit_out,100) dia_nacelle, '= dia_nacelle' write(unit_out,100) sh, '= sh' write(unit_out,100) rho_eo, '= rho_eo' write(unit_out,100) a_eo, '= a_eo' write(unit_out,103) mu_eo, '= mu_eo' write(unit_out,100) dr_max, '= dr_max' write(unit_out,100) xbar, '= xbar' write(unit_out,100) cbar, '= cbar' end if c Echo input to output file unit_outfile = 20 write(unit_outfile,90) title write(unit_outfile,*) write(unit_outfile,"(a25)") ' Input Data' write(unit_outfile,*) write(unit_outfile,101) write_flag, '= write_flag' write(unit_outfile,100) dihedral_wing, '= dihedral_wing' write(unit_outfile,100) z_wing, '= z_wing' write(unit_outfile,100) dia_fuse, '= dia_fuse' write(unit_outfile,100) sref, '= sref' write(unit_outfile,100) hspan_vtail, '= hspan_vtail' write(unit_outfile,100) depth_fuse_vtail, '= depth_fuse_vtail' write(unit_outfile,100) c_vtail_root, '= c_vtail_root' write(unit_outfile,100) c_vtail_tip, '= c_vtail_tip' write(unit_outfile,100) mach_eo, '= mach_eo' write(unit_outfile,100) sweep_vtail_1_4*180./pi,'=sweep_vtail_1_4' write(unit_outfile,100) sweep_wing_1_2*180./pi, '= sweep_wing_1_2' write(unit_outfile,100) hspan_wing, '= hspan_wing' write(unit_outfile,100) cl, '= cl' write(unit_outfile,100) z_vtail, '= z_vtail' write(unit_outfile,100) l_vtail, '= l_vtail' write(unit_outfile,100) length_fuse, '= length_fuse' write(unit_outfile,101) new, '= new' write(unit_outfile,101) nef, '= nef' write(unit_outfile,100) dia_nacelle, '= dia_nacelle' write(unit_outfile,100) sh, '= sh' write(unit_outfile,100) rho_eo, '= rho_eo' write(unit_outfile,100) a_eo, '= a_eo' write(unit_outfile,103) mu_eo, '= mu_eo' write(unit_outfile,100) dr_max, '= dr_max' write(unit_outfile,100) xbar, '= xbar' write(unit_outfile,100) cbar, '= cbar' c Calculate stability and control derivatives via Roskam's methods c Sideslip angle derivatives cy_beta_wing = -0.0001*abs(dihedral_wing)*180./pi c Estimate fig_7_1 from Figure 7.1 (curve fit) if (z_wing/(dia_fuse/2.) .le. 0.) then fig_7_1 = 0.85*(-z_wing/(dia_fuse/2.)) + 1. elseif (z_wing/(dia_fuse/2.) .gt. 0.) then fig_7_1 = 0.5*z_wing/(dia_fuse/2.) + 1. end if c Estimate the side force coefficient due to the fuselage and nacelles cy_beta_fuse = -2.*fig_7_1*( pi*(dia_fuse/2.)**2 + & (new + nef)*pi*(dia_nacelle/2.)**2 )/sref c Estimate fig_7_3 from Figure 7.3 (curve fit) x = hspan_vtail/depth_fuse_vtail if (x .le. 2.) then fig_7_3 = 0.75 elseif (x .gt. 2. .and. x .lt. 3.5) then fig_7_3 = x/6. + 5./12. elseif (x .ge. 3.5) then fig_7_3 = 1. end if c Estimate fig_7_5 from Figure 7.5 (curve fit for taper ratio <= 0.6) x = hspan_vtail/depth_fuse_vtail fig_7_5 = 0.002*x**5 - 0.0464*x**4 + 0.404*x**3 - 1.6217*x**2 + & 2.7519*x + 0.0408 c Factor from Figure 7.6 is for zh/bv = 0. fig_7_6 = 1.1 c Estimate fig_7_7 from Figure 7.7 (curve fit) sv = hspan_vtail*(c_vtail_root + c_vtail_tip)/2. x = sh/sv fig_7_7 = -0.0328*x**4 + 0.2885*x**3 - 0.9888*x**2 + 1.6554*x - & 0.0067 c Estimate the effective aspect ratio for the vertical tail ar_vtail = hspan_vtail**2/sv ar_vtail_eff = fig_7_5*ar_vtail*(1. + fig_7_7*(fig_7_6 - 1.)) c Assume the section lift-curve slope is 2.*pi cl_alpha_vtail = 2.*pi c Estimate the effective lift-curve slope for the vertical tail kappa = cl_alpha_vtail/(2.*pi) beta_m = sqrt( 1. - mach_eo**2 ) sweep_vtail_1_2 = atan( (c_vtail_root/4. + hspan_vtail* & tan(sweep_vtail_1_4) + c_vtail_tip/4. - c_vtail_root/2.)/ & hspan_vtail ) cl_alpha_vtail_eff = 2*pi*ar_vtail_eff/( 2. + & sqrt( ar_vtail_eff**2*beta_m**2/kappa**2*( 1. + & tan(sweep_vtail_1_2)**2/beta_m**2 ) + 4. ) ) c Estimate the third term in eqn. 7.4 from eqn. 7.5 eqn_7_5 = 0.724 + 3.06*sv/sref/(1. + cos(sweep_vtail_1_4)) + & 0.4*z_wing/dia_fuse + 0.009*ar_vtail_eff cy_beta_vtail = -fig_7_3*cl_alpha_vtail_eff*eqn_7_5*sv/sref c Calculate total variation of side force coefficient with sideslip angle cy_beta = cy_beta_wing + cy_beta_fuse + cy_beta_vtail c Factor from Figure 7.11 is approximated by a curve fit for lambda = 0.5 fig_7_11 = -0.004/45*sweep_wing_1_2*180./pi c Factor from Figure 7.12 is approximated for 747 and 777 configurations c at low Mach numbers fig_7_12 = 1.0 c Factor from Figure 7.13 is approximated for 747 and 777 configurations fig_7_13 = 0.85 c Factor from Figure 7.14 is approximated for lambda = 0.5 and high AR fig_7_14 = 0.000 c Factor from Figure 7.15 is approximated by a linear curve fit for c lambda = 0.5, low sweep, and high AR ar = (2.*hspan_wing)**2/sref fig_7_15 = -0.00012 - 0.00013/10*ar c Factor from Figure 7.16 is approximated for 747 and 777 configurations c at low Mach numbers fig_7_16 = 1.0 c Estimate body-induced effect on wing height from eqns. 7.10, 7.11, and 7.12 eqn_7_11 = sqrt(pi*(dia_fuse/2.)**2/0.7854) eqn_7_10 = -0.0005*sqrt(ar)*(eqn_7_11/(2.*hspan_wing))**2 eqn_7_12 = -1.2*sqrt(ar)/(180./pi)*z_wing/(2.*hspan_wing)* & 2.*eqn_7_11/(2.*hspan_wing) c Wing-body contribution to cl_beta (wing twist effect is neglected) cl_beta_wingbody = ( cl*(fig_7_11*fig_7_12*fig_7_13 + & fig_7_14) + dihedral_wing*(fig_7_15*fig_7_16 + eqn_7_10) + & eqn_7_12 )*180./pi c Since the horizontal tail has a small lift coefficient, small dihedral, c and small area relative to the wing, it is negligible. cl_beta_htail = 0. c Estimate the angle of attack from the CL by assuming cl_alpha = 5.7 c based on the value given for the 747 in Nelson. The angle of attack c of the fuselage reference line is given as 5.7 deg for the M=0.25 c flight condition in NASA CR-2144. Therefore, the effective wing c incidence angle (with flaps at 20 deg) is 5.457 deg. cl_alpha = 5.7 alpha = cl/cl_alpha - 5.457*pi/180. c Estimate the vertical tail contribution to cl_beta cl_beta_vtail = cy_beta_vtail*( z_vtail*cos(alpha) - l_vtail* & sin(alpha) )/(2.*hspan_wing) c Calculate total variation of rolling moment coefficient with sideslip angle cl_beta = cl_beta_wingbody + cl_beta_htail + cl_beta_vtail c Wing contribution to cn_beta is negligible for small angles of attack. c cn_beta_wing = 0. c Calculate contribution of wing to cn_beta using DATCOM method. c Code added by L.T. Leifsson, code written by W.H. Mason, March 2004. cnbw1 = 1.0/(4.0*pi*ar) cnbw2 = tan(sweep_wing_1_2)/pi/ar/(ar + 4.0*cos(sweep_wing_1_2)) cnbw3 = cos(sweep_wing_1_2) cnbw4 = ar/2.0 cnbw5 = ar**2.0/8.0/cos(sweep_wing_1_2) cnbw6 = 6.0*xbar/cbar*sin(sweep_wing_1_2)/ar cn_beta_wingterm = cnbw1 - cnbw2*(cnbw3 - cnbw4 - cnbw5 + cnbw6) cn_beta_wing = cl**2.0*cn_beta_wingterm c Estimate empirical factor for body and body + wing effects from Figure 7.19 c Constant value assumed for 747 and 777-like configurations fig_7_19 = 0.0011 c Calculate fuselage Reynolds number at the engine-out flight condition re_fuse = rho_eo*mach_eo*a_eo*length_fuse/mu_eo c Estimate fuselage Reynolds number effect on wing-body from Figure 7.20 fig_7_20 = 1. + 1.2/log(350.)*log(re_fuse/1000000.) c Estimate fuselage contribution to cn_beta sbs = 0.83*dia_fuse*length_fuse cn_beta_fuse = -180./pi*fig_7_19*fig_7_20*sbs/sref* & length_fuse/(2.*hspan_wing) c Estimate vertical tail contribution to cn_beta cn_beta_vtail = -cy_beta_vtail*( l_vtail*cos(alpha) + & z_vtail*sin(alpha) )/(2.*hspan_wing) c Calculate total variation of rolling moment coefficient with sideslip angle cn_beta = cn_beta_wing + cn_beta_fuse + cn_beta_vtail c Assume variation of sideforce coefficient with aileron deflection is zero cy_da = 0. c Estimate the rolling moment effectiveness parameter from Figure 11.1 c for lambda = 0.5, and for 747 and 777-like ailerons at mach 0.2 fig_11_1 = 0.18 c Assume the 2D lift-curve slope is 2*pi/beta_m cl_alpha_2d = 2*pi/beta_m kappa = cl_alpha_2d/(2.*pi/beta_m) c Estimate the rolling effectiveness of two full-chord controls by Eqn. 11.2 eqn_11_2 = kappa/beta_m*fig_11_1 c Estimate aileron effectiveness by assuming cf/c = 0.20 and t/c = 0.08 fig_10_5 = 3.5 fig_10_6 = 1.0 cl_d = fig_10_6*fig_10_5 alpha_d = cl_d/cl_alpha_2d c Determine the rolling effectiveness of the partial-chord controls by c Eqn. 11.3. Note that this is the change in cl with respect to a change c in the sum of the left and right aileron deflections (d). cl_d = alpha_d*eqn_11_2 c Estimate variation of rolling moment coefficient with aileron deflection c by neglecting differential control effects. Since the aileron deflection c (da) is defined as half of the sum of the left and right deflections, cl_d c from the equation above must be divided by 2. cl_da = cl_d/2. c The method in Roskam for estimating cn_da does not account for the c effect of differential ailerons and the use of spoilers for roll control c on the yaw moment. Therefore, the factor k is estimated c based on the ratio of cn_da to cl_da from the 747 flight test data c presented in Nelson. Note that the effect of cl is absorbed into c the factor k. k = 0.0064/0.0461 c Estimate variation of yawing moment coefficient with aileron deflection cn_da = k*cl_da c Estimate the flap chord factor from Figure 10.2 for cf/c = 0.33 c The flap effectiveness ratio is estimated with a piecewise curve fit cf_c = 0.33 alpha_d_cl = -sqrt( 1. - (1. - cf_c)**2 ) if (alpha_d_cl .ge. -0.5) then flap_eff_ratio = 1.42 + 1.8*alpha_d_cl elseif (alpha_d_cl .ge. -0.6) then flap_eff_ratio = 1.32 + 1.6*alpha_d_cl elseif (alpha_d_cl .ge. -0.7) then flap_eff_ratio = 1.08 + 1.2*alpha_d_cl else flap_eff_ratio = 0.94 + alpha_d_cl end if flap_eff_ratio = 1. + flap_eff_ratio/( ar_vtail_eff - & 0.5*(-alpha_d_cl - 2.1) ) fig_10_2 = flap_eff_ratio*alpha_d_cl c Estimate empirical correction for lift effectiveness of plan flaps at c from Figure 10.7 for cf/c = 0.33. x = dr_max if (x .lt. 15.) then fig_10_7 = 1. else fig_10_7 = 4e-7*x**4 - 7e-5*x**3 + 0.0047*x**2 - 0.1453*x + & 2.3167 end if c Estimate span factor for plain flap from Figure 10.3 for delta eta = 0.85 fig_10_3 = 0.95 c Estimate variation of sideforce coefficient with rudder deflection cy_dr = cl_alpha_vtail_eff*fig_10_2*fig_10_7*fig_10_3*sv/sref c Estimate variation of rolling moment coefficient with rudder deflection cl_dr = cy_dr*( z_vtail*cos(alpha) - l_vtail*sin(alpha) )/ & (2.*hspan_wing) c Estimate variation of yawing moment coefficient with rudder deflection cn_dr = -cy_dr*( l_vtail*cos(alpha) + z_vtail*sin(alpha) )/ & (2.*hspan_wing) c Multiply empirical estimates by their respective correction factors c The correction factors are the ratio of the actual 747 derivatives to c the 747 derivatives predicted by the method above at the M=0.25 flight c condition defined in NASA CR-2144 and Nelson. The rudder deflection c was 15 deg for this calibration. cy_beta = 1.4068*cy_beta cl_beta = 0.7253*cl_beta cn_beta = 1.4232*cn_beta ! New correction factor, Cnbetawing now cl_da = 0.9202*cl_da ! estimated. Assuming xbar = 0.0 ft and cn_da = 0.9143*cn_da ! cbar = 27.31 ft, for 747,Leifsson, March '04. cy_dr = 0.6132*cy_dr cl_dr = 0.3004*cl_dr cn_dr = 0.7320*cn_dr c Set the rudder deflection to 20 deg, and the bank angle to 5 deg dr = dr_max*pi/180. phi = 5.*pi/180. c Solve for the sideslip angle and aileron deflection beta = -( cy_dr*dr + cl*sin(phi) )/cy_beta da = -( cl_dr*dr + cl_beta*beta)/cl_da c Calculate the maximum available yawing moment coefficient cn_avail = cn_da*da + cn_dr*dr + cn_beta*beta c Write output data if (write_flag .eq. 1) then write(unit_out,*) write(unit_out,"(a26)") ' Output Results' write(unit_out,96) write(unit_out,104) cy_beta_wing, '= cy_beta_wing' write(unit_out,104) cy_beta_fuse, '= cy_beta_fuse' write(unit_out,104) cy_beta_vtail, '= cy_beta_vtail' write(unit_out,*) write(unit_out,105) cy_beta, '= cy_beta' write(unit_out,*) write(unit_out,104) cl_beta_wingbody, '= cl_beta_wingbody' write(unit_out,104) cl_beta_htail, '= cl_beta_htail' write(unit_out,104) cl_beta_vtail, '= cl_beta_vtail' write(unit_out,*) write(unit_out,105) cl_beta, '= cl_beta' write(unit_out,*) write(unit_out,104) cn_beta_wing, '= cn_beta_wing' write(unit_out,104) cn_beta_fuse, '= cn_beta_fuse' write(unit_out,104) cn_beta_vtail, '= cn_beta_vtail' write(unit_out,*) write(unit_out,105) cn_beta, '= cn_beta' write(unit_out,97) write(unit_out,100) cy_da, '= cy_da' write(unit_out,100) cl_da, '= cl_da' write(unit_out,100) cn_da, '= cn_da' write(unit_out,*) write(unit_out,100) cy_dr, '= cy_dr' write(unit_out,100) cl_dr, '= cl_dr' write(unit_out,100) cn_dr, '= cn_dr' write(unit_out,98) write(unit_out,100) beta*180./pi, '= beta (deg)' write(unit_out,100) phi*180./pi, '= phi (deg)' write(unit_out,100) da*180./pi, '= da (deg)' write(unit_out,100) dr*180./pi, '= dr (deg)' write(unit_out,100) ar_vtail_eff, '= ar_vtail_eff' write(unit_out,100) cn_avail, '= cn_avail' write(unit_out,99) endif c Write results to output file write(unit_outfile,*) write(unit_outfile,"(a26)") ' Output Results' write(unit_outfile,96) write(unit_outfile,104) cy_beta_wing, '= cy_beta_wing' write(unit_outfile,104) cy_beta_fuse, '= cy_beta_fuse' write(unit_outfile,104) cy_beta_vtail, '= cy_beta_vtail' write(unit_outfile,*) write(unit_outfile,105) cy_beta, '= cy_beta' write(unit_outfile,*) write(unit_outfile,104) cl_beta_wingbody, '= cl_beta_wingbody' write(unit_outfile,104) cl_beta_htail, '= cl_beta_htail' write(unit_outfile,104) cl_beta_vtail, '= cl_beta_vtail' write(unit_outfile,*) write(unit_outfile,105) cl_beta, '= cl_beta' write(unit_outfile,*) write(unit_outfile,104) cn_beta_wing, '= cn_beta_wing' write(unit_outfile,104) cn_beta_fuse, '= cn_beta_fuse' write(unit_outfile,104) cn_beta_vtail, '= cn_beta_vtail' write(unit_outfile,*) write(unit_outfile,105) cn_beta, '= cn_beta' write(unit_outfile,97) write(unit_outfile,100) cy_da, '= cy_da' write(unit_outfile,100) cl_da, '= cl_da' write(unit_outfile,100) cn_da, '= cn_da' write(unit_outfile,*) write(unit_outfile,100) cy_dr, '= cy_dr' write(unit_outfile,100) cl_dr, '= cl_dr' write(unit_outfile,100) cn_dr, '= cn_dr' write(unit_outfile,98) write(unit_outfile,100) beta*180./pi, '= beta (deg)' write(unit_outfile,100) phi*180./pi, '= phi (deg)' write(unit_outfile,100) da*180./pi, '= da (deg)' write(unit_outfile,100) dr*180./pi, '= dr (deg)' write(unit_outfile,100) ar_vtail_eff, '= ar_vtail_eff' write(unit_outfile,100) cn_avail, '= cn_avail' write(unit_outfile,99) c Close output file close(unit_outfile) c Write statements and formats 90 format(/4x,'Data set title: ',a72) 96 format(/3x,'stability derivatives'/) 97 format(/3x,'control derivatives'/) 98 format(/3x,'engine out information - 5 deg bank and max rudder'/) 99 format(/3x,'Calculations completed') 100 format(3x,f11.4, 1x, a) 101 format(10x, i4, 1x, a) 102 format(3x,f11.0, 1x, a) 103 format(3x,g11.4, 1x, a) 104 format(22x,f11.4, 1x, a) 105 format(3x,f11.4, 1x, a, ' (adjusted)') c Convert sweep angles from radians back to degrees cwhm THIS IS A VERY DANGEROUS PRACTICE!! DON'T USE AS AN EXAMPLE cwhm of acceptable coding style sweep_wing_1_2 = sweep_wing_1_2*180./pi sweep_vtail_1_4 = sweep_vtail_1_4*180./pi return end