' PROGRAM NAME: DESIGN.BAS (C:\F2A\DESIGN/F2AGP.BAS ) ' ************************************************************* ' Design a propeller entirely from blade element analysis algorithms. ' Chord distribution must satisfy constant slip. ' Blade angle-of-attack arbitrarily chosen, but somewhere near best L/D ' Corrected by Prandtl-Glauert scaling rule ' 25 airfoils stored for use in C:\FITZGF2C\CUTMETAL\NACAF2C.BAS ' ************************************************************* ' Method is based on work of the late E. Eugene Larrabee, published ' in the 1979 NFFS symposium record, under the title "Propeller ' Design and Analysis for Modelers" ' In particular, the equations are from the section "Arbitrary Propeller ' Performance". This section was intended to permit the performance ' analysis of a propeller of known geometry. That is to say, where the ' blade angle, chord and airfoil parameters are known radially. ' However, in the following code, an inverse treatment is used, which ' yields a propeller design of known parameters by an iterative process. ' The resulting propeller has a known power absorption and is of minimum ' induced loss. ' The present code makes no allowance for Reynolds number or the effects ' of the formation of shock waves. ' However, it does allow for compressibilty, by using data corrected ' by the Prandtl-Glauert scaling rule. This data come from the airfoil ' at station 18 of F2C23, corrected back to zero mach. The remaining ' airfoils are then generated in CONVERT.BAS ' The process is as follows: ' 1. Choose a constant-chord blade with width somewhere near that desired ' 2. Assign a value of angle of attack for each blade element which gives ' good L/D. ' 3. Calculate the power absorption and slip for each blade element of ' the constant-chord prop. ' 4. Adjust the blade width and chord to reduce slip variation and ' make the power absorption closer to that desired. ' 5. Iterate, using the new design.ie, go back to step 3 above ' 6. Keep iterating until the slip is constant and the power absorption ' as desired. ' 7. Save the compressibilty corrected airfoils ' Some aspects of the coding are novel. Airfoil parameters are generally ' derived from the chordline location. eg blade angle and camber. ' Pitch in particular is measured to chordline using a pitch gauge ' instrument. Regretably, such a pitch quantity has no fundamental value ' as it tells us nothing about the aerodynamic qualities of the section. ' ' A better quantity, which simplifies the appearance of the maths, is ' the zero-lift angle. Unfortunately this angle is tricky to find, unless ' perhaps you have recourse to a wind tunnel or airfoil data tables. ' Luckily there is a very simple relationship between the zero-lift angle ' Az0 and the airfoil camber M and camber-line high point location P. ' This relationship is Az0 = ATN(M/(1-P)) , radians ' The accuracy is quite good, at least as good as more complex methods ' given by Abbott and Von Doenhoff in the classic "Theory of Wing ' Sections" NACA derived text. ' Accordingly, in the following code I have changed Larrabee's definition ' of the blade angle "beta", from being measured from the plane of rotation ' to the chordline, to measurement to the zero lift angle. This also ' permits the use of lift slope data where the origin at zero-lift is also ' at zero blade angle, as measured to the zero-lift line. ' The code is heavily commented, so hope fully it can be followed more ' readily than Larabees text. ' The reader may not have come across the concept of "interference factors" ' previously. They arise naturally from the way airfoil lift-slope data ' are presented. The plots of lift coefficient versus angle-of-attack ' are for airfoils of infinite aspect ratio. Clearly that is not the case ' we have with propellers or even "real" wings. So a correction for ' aspect-ratio is necessary. But what is the aspect ratio of a propeller? ' With infinite aspect ratio wings, the lift vector is perpendicular to ' the direction of motion. However, as aspect ratio is reduced, the lift ' vector angles back, so that the vector has a component retarding the ' motion. This is the quantity referred to as "induced drag". The ' components of the angled vector are the "interference factors" ' The lift of the section results from downwash of the air by the section. ' The components of downwash are axial and tangential, the latter resulting ' in "swirl" of the slipstream about the rotation axis. These quantites ' are, somewhat surprisingly, amenable to calculation. The trick is to match ' the downwash momentum change to the aerodynamic forces. Thanks to some ' clever fellows, in particular Joukowski, Prandtl, Glauert, Goldstein ' and not least, Larrabee, this has all been done for us. ' Now read on. ' Nota bene: The variables below let you get any answer you want. ' As a guide, apply Rothwell's law, by setting the chordline ' parallel to the induced inflow W. ' Variables: ' Enviromental (subroutine atmos) Temperature Ta Centigrade ' Pressure Pa kPa ' Altitude altitude feet ' Propeller Spinner diameter spin inches ' Power absorption P hp ' RPM RPM ' Diameter Ddes inches ' Blade number B ' Coding Angle of attack AA degrees ' (advance angle to zero-lift, at Mach = 0) ' Number of iterations Niterations ' drive$ = "C:\FITZGF2C\DESIGN\ ' ' Subroutines: ' nacafoil: sets trila value for chords (all equal) ' atmos: set air density ' rotate2: rotational tranformation ' rotate: " " ' vectors: plot vector diagram of airflow ' chordplot: plot chord radial distribution ' slipplot: plot slip ' plotchord2: plot chord radial distribution ' delay: wait for display ' applyPG: apply Prandtl-Glauert scaling to generic airfoil ' area: get airfoil C/G coordinates ' adder: optional, add chord at root, using scaling factor G ' generic: load generic airfoil coords. (Mach = 0) ' Must get this from other programs ' eg, F2C23.BAS and CONVERT.BAS ' printdata: store airfoil parameters for printing as text file ' saveairfoil: stores 25 airfoils corrected for Prandtl-Glauert ' rule. Airfoils at stations 1,2,3 identical ' Warning: Scaling of figures is set for F2C prop: may be no good for ' other sizes, but program should still run ' Scaling is done at the start of subprogram "vectors:" ' Arrays used in this code: DIM mach(25), beta(25), pitch(25), fpitch(25) DIM Re(25), ch(25), F(25) DIM phi1(450), phi2(450), Cl1(450), Cd1(450), Cy1(450), Cx1(450) DIM aa1(450), a11(450), dCt(450), dCp(450) DIM M2(25), P.p(25), yc(100), T(25) DIM faceangle(25), BCl0(25), face(25) DIM epitch(25), sol(25), slip(25), te(25) DIM Wx(450), alpha(450), phi(25) DIM advance(25), GP(25), mm(25) DIM xu(100), xl(100), yu(100), yl(100) DIM xur(100), yur(100), xlr(100), ylr(100) DIM yt(100), xord(100) DIM yua(25, 101), zua(25, 101), yla(25, 101), zla(25, 101) DIM station(30), tchord(30), beta1(30) DIM xurt(101), yurt(101), xlrt(101), ylrt(101) DIM xus(25, 100), yus(25, 100), zus(25, 100) DIM xls(25, 100), yls(25, 100), zls(25, 100) DIM PGrule(25) DIM suby(25), subz(25), angle(25), AA(25) DIM zur(101), zlr(101) DIM ylrnew(101), zlrnew(101), yurnew(101), zurnew(101) DIM yusg(101), zusg(101), ylsg(101), zlsg(101) DIM yurnext(101), ylrnext(101), zurnext(101), zlrnext(101) DIM dum(25), tc(25) DIM tcGP(25), zerolift(25) DIM betax(25), P(25), area(210), areatotal(210), cgx(25), cgy(25), cg(25) DIM g(25), camber(25), xcamber(101), ycamber(101) DEFSNG A-Z SCREEN 9, , 1, 1: pi = 4 * ATN(1): rtd = 180 / pi: dtr = pi / 180 ' rtd converts radians to degrees, and dtr conversely. flag = 1 ' plot all steps flag = 0 ' don't plot drive$ = "C:\FITZGF2C\DESIGN\" '--------------------------------------------------------------------- ' INPUT VARIABLES spin = .7 ' spinner diameter in inches P = 1.5: Pdes = P ' sea level Power (BHP) at 15 C, 101.325 kPa RPM = 40000: RPMdes = RPM ' RPM vV = 180: VVVdes = vV ' Flight speed of aircraft (MPH) Ddes = 75 * 2 / 25.4 ' Diameter of prop in inches R = Ddes / 2 ' Propeller tip radius (inches) R = R * .0254 ' " " " (metres) B = 1 ' Number of propeller blades Niterations = 4 ' Number of iterations ' --------------------------------------------------------------------- GOSUB atmos ' air density correction P = P * rhoratio ' Adjust engine power for air density Mo = .594 * Taa + 325.56 ' Speed of sound BBB = INT(spin / 2 / R * 25 / 25.4) ' Working part of blade, elements ' outside spinner PRINT BBB: INPUT vbv ' ****************************************************************** ' Select airfoil data GOSUB nacafoil ' Choose intial chord distribution ' ------------------------------------------------------------------ spin = spin * .0254 ' spinner diameter, metres P = P * 746 ' BHP to WATTS Vdes = vV * 2.54 * 12 * 5280 / 60 / 60 / 100 ' airspeed, MPH to m/sec V FOR i = 1 TO 25 station(i) = i * R / 25 * 1000 ' metres Vr = Vdes ' m/s mach(i) = SQR((i * R * RPMdes * .0010472 * 4) ^ 2 + Vr ^ 2) / Mo IF mach(i) > 1 THEN PRINT "Tip speed exceeds Mach 1" IF mach(i) > 1 THEN STOP advance = Vdes / (RPMdes / 60) ' metres per revolution angle(i) = ATN((advance / (2 * pi * station(i)))) * rtd PGrule(i) = SQR(1 / (1 - mach(i) ^ 2)) NEXT i ' INT(spin / 2 / R * 25) ' Calculations outside spinner FOR i = 1 TO 25 face(i) = ATN(1.1019 * T(i) ^ 2) ' angle between chordline and faceangle(i) = face(i) ' and airfoil face (bottom) NEXT i 'ÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ ' CALCULATION BEGINS ' Load chosen airfoil for Mach = 0 GOSUB generic INPUT "generic", vbv ' POWER MATCHING LOOP OPEN drive$ + "airfoil.dat" FOR OUTPUT AS #4 Endloop = Niterations ' Number of iterations ' on slip and power absorption FOR Powerloop = 1 TO Endloop ' Powerloop is iteration counter RPM = RPMdes: vV = VVVdes: COLOR 4 PRINT " "; prop$: PRINT : COLOR 2 PRINT " RPM = ", RPM PRINT " MPH = ", vV: COLOR 9: PRINT PRINT " "; Powerloop, Endloop: PRINT : COLOR 6 ' ADVANCE RATIO (lam) v = vV * 2.54 * 12 * 5280 / 60 / 60 / 100 ' m.s-1 om = RPM * 2 * pi / 60 lam = v / (om * R) ' ADVANCE RATIO '--------------------------------------------------------------- ' CALCULATE AT EACH STATION, THRUST, SLIP AND POWER ABSORPTION '--------------------------------------------------------------- ' MOVE ALONG BLADE ELEMENTS FOR j = BBB TO 25 Vr = v ' lamr = Vr / (om * R) ' advance ratio at element J zet = j / 25 - .02 ' offset to blade element centre (bin) bin = .04 ' fractional width of blade element OMr = om * zet * R ' VEL. AT POINT R IN PLANE OF ROTATION dum = (B / 2 * (1 - zet) / lamr * SQR(lamr ^ 2 + 1)) F(j) = 2 / pi * ATN(SQR(EXP(2 * dum) - 1)) ' Prandtl correction ' Calculate initial solidity.. amount of the prop disc shadowed by blades zeta = j / 25 - .02 ' blade elements z = 2 * pi * R * zeta sol(j) = B * ch(j) * .001 / z ' solidity ' free-stream inflow angle is advance(J) in radians. Angle is unaffected ' by the approaching airfoil. Corresponds to infinite aspect ratio wing advance(j) = ATN(Vr / OMr) ' ------------------------------------------------------------------ ' The following routine matches the radial and axial momentum change ' to the corresponding thrust and torque forces. To do this, an angle ' of attack is chosen (at zero lift), interference factors determined ' and balance with momentum change established ' ------------------------------------------------------------------ AA = 6 / (PGrule(j) ^ .5) ' Chosen angle of attack (degrees from AOA = AA + zeroliftz ' zero-lift to advance) ' Get airfoil parameters corrected for compressibilty GOSUB applyPG ' apply Prandtl-Glauert to generate ' airfoil J, then get reduction in ' angle of attack and zero-lift angle Az0 = zerolift(j) ' zero-lift angle relative to chordline ' from P-G subroutine applyPG ' Now consider advance ratio at current blade station ' Set initial values ll = 1 ' counter AOA = AA(j) ' angle between zero-lift and advance angle alpha(ll) = .5 * dtr ' start just above zero lift beta(j) = advance(j) + AOA ' Fixd value for beta ' NB: beta is blade angle measured ' to airfoil zero-lift line,NOT chordline ' Enter loop to get interference values 1930 phi1(ll) = beta(j) - alpha(ll) ' phi..angle to stream inflow for ' infinite aspect ratio data Cl1 = .1 * AA ' value .1 is lift slope deg^-1 ' for infinite aspect ratio lift-slope. Fixed ' value of Cl applies as GP rule normalises to mach = 0 Cd1 = .02 ' Choose drag coeff. Cp = COS(phi1(ll)): Sp = SIN(phi1(ll)) Cy1 = Cl1 * Cp - Cd1 * Sp ' Thrust load on blade element Cx1 = Cl1 * Sp + Cd1 * Cp ' Torque " " " " ' axial interference .. axial downwash fraction from airfoil z = sol(j) * Cy1 / (4 * F(j) * Sp ^ 2) aa1 = z / (1 - z) ' swirl interference .. tangemtial downwash fraction from airfoil z = sol(j) * Cx1 / (4 * F(j) * Sp * Cp) a11 = z / (1 + z) phi2(ll) = ATN(lamr / zet * (1 + aa1) / (1 - a11)) Wx(ll) = OMr * (1 - a11) / COS(phi1(ll)) ' inflow velocity for ' infinite aspect ratio IF ABS(phi1(ll) - phi2(ll)) < .0001 THEN GOTO 2390 ll = ll + 1 alpha(ll) = alpha(ll - 1) + .1 * (phi1(ll - 1) - phi2(ll - 1)) GOTO 1930 ' continue interference loop ' ------------------------------------------------------- 2390 ' Calculate power, thrust, efficiency and plot vectors phi(j) = phi2(ll) GOSUB vectors GOSUB delay Cl1(j) = Cl1: Cd1(j) = Cd1 OMr = om * zet * R ' VEL. AT POINT R IN PLANE OF ROTATION aV = aa1 * Vr ' AXIAL INDUCED VELOCITY aOMr = a11 * OMr ' SWIRL INDUCED VELOCITY slip(j) = aOMr * TAN(phi(j)) + aV ' lost motion from advance to W ' Slip .. Larrabee slip(j) = OMr * (TAN(phi(j)) - TAN(advance(j)))' lost motion from advance to W dum = ((1 - a11) / COS(phi(j))) ^ 2 dCt(j) = pi ^ 3 / 4 * dum * zet ^ 3 * sol(j) * Cy1 ' thrust coeff. dCp(j) = pi ^ 4 / 4 * dum * zet ^ 4 * sol(j) * Cx1 ' power coeff. ' Test for last iteration, then save airfoils IF Powerloop = Endloop THEN GOSUB saveairfoil NEXT j ' continue to next blade element in power loop '------------------------------------------------------------ GOSUB chordplot ' plot chord distribution GOSUB slipplot ' " slip " ' INTEGRATE TO GET THRUST AND POWER COEFFICIENT s = 0: FOR j = BBB TO 25: s = s + dCt(j): NEXT j: Ct = bin * s s = 0: FOR j = BBB TO 25: s = s + dCp(j): NEXT j: Cp = bin * s ' n is prop efficiency n = pi * lam * Ct / Cp ' prop efficiency ' Power absorption P = Cp * rho * (RPM / 60) ^ 3 * (2 * R) ^ 5 / 746 ' Horsepower ' Thrust T = Ct * rho * (RPM / 60) ^ 2 * (2 * R) ^ 4 ' Newtons Torq = 1 / (2 * pi) * rho * (RPM / 60) ^ 2 * (2 * R) ^ 5 * Cp' Newton metres THP = P * n ' Thrust-power: watts PRINT PRINT "POWER (HP)= ", P: PRINT "THRUST (lbs)=", T / 4.45 PRINT "THRUST POWER (lbs)=", P * n PRINT "Ntotal =", n * 100 PRINT " N slip = ", (1 - slip(18) / (slip(18) + Vr)) * 100 GOSUB delay IF Powerloop = Endloop THEN 2345 ELSE 2346 ' exit power loop ' otherwise continue iteration for power/slip ' ADJUST CHORD TO MATCH DESIGN POWER ABSORPTION AND SLIP ' The object is to set the radial chord disritbution to ' achieve the required power absorption, and get the slip radially ' constant for minimum induced loss 2346 FOR i = BBB TO 25 zeta = i / 25 - .02: z = 2 * pi * R * zeta ' ch(i) = ch(i) * (Pdes / P) ^ .5 * slip(18) / slip(i) ' new chord ' To match power only, use line below ch(i) = ch(i) * (Pdes / P) ^ .5 ' new chord sol(i) = B * ch(i) * .001 / z ' new solidity 543 NEXT i GOSUB plotchord2 NEXT Powerloop ' End power loop, display results CLOSE #4 ' ---------------------------------------------------------- 2345 ' Calculate pitch to chordline FOR i = BBB TO 25 zeta = i / 25 - .02 z = 2 * pi * R * zeta pitch(i) = z / .0254 * TAN(beta(i) + zerolift(i)) ' pitch to chordline fpitch(i) = z / .0254 * (TAN(beta(i) + zerolift(i)) - faceangle(i)) ' face epitch(i) = z * TAN(beta(i)) / .0254 ' experimental pitch NEXT i ' (measured to zero-lift) CLS : LOCATE 1, 1 COLOR 12: BEEP: PRINT : PRINT : PRINT " ", prop$; PRINT " DESF2CGP.BAS (with correction for compressibilty)" PRINT : COLOR 2 PRINT " RPM = "; : PRINT USING "#####"; RPMdes; PRINT USING " V = ### MPH"; VVVdes; PRINT USING " Diameter = ###.## inches (##.#)"; Ddes; Ddes * 25.4 / 2 PRINT " Blade number = "; : PRINT USING "##"; B; PRINT USING " No. of iterations ###"; Niterations PRINT PRINT USING " POWER = ###.## (hp)"; P PRINT USING " THRUST = ###.## (lbs)"; T / 4.45 PRINT USING " THRUST POWER = ###.## (hp) "; P * n PRINT USING " Efficiency = ####.# %"; n * 100: PRINT ' GOSUB delay COLOR 10 PRINT "I Cl Cl/Cd PITCH ePITCH Ch BETA Advance zero-angle AOA Mach" PRINT " (ins) (ins) mm deg deg deg deg " ' FOR i = 3 TO 5 ' dum = ch(i) * SIN(beta(i)) - ch(6) * SIN(beta(6)) ' sub1 = dum / SIN(beta(i)) ' ch(i) = ch(i) - sub1 ' NEXT i ' ch(1) = ch(3): beta(1) = beta(3) ' ch(2) = ch(3): beta(2) = beta(3) FOR i = BBB + 1 TO 25 STEP 2: COLOR 12 PRINT USING "##"; i; PRINT USING " #.###"; Cl1(i); : PRINT " "; PRINT USING "##.##"; Cl1(i) / Cd1(i); : PRINT " "; PRINT USING "##.#"; pitch(i); : PRINT " "; PRINT USING "##.#"; epitch(i); : PRINT " "; PRINT USING "###.#"; ch(i); : PRINT " "; PRINT USING " ##.#"; beta(i) * rtd; : PRINT " "; PRINT USING " ##.#"; advance(i) * rtd; : PRINT " "; PRINT USING " ####.#"; zerolift(i) * rtd; PRINT USING " ####.#"; (AA(i) + zerolift(i)) * rtd; PRINT USING " #.##"; mach(i) NEXT i INPUT vbv ' GOSUB printdata ' file .txt for hardcopy ' COLOR 0 CLS : LOCATE 2, 1: COLOR 12 PRINT : PRINT : PRINT PRINT " I Station Ch m P beta t/c Mach" PRINT " (mm) " FOR i = 3 TO 25 STEP 2: COLOR 12 PRINT USING " ##"; i; PRINT USING " ###.#"; station(i); : PRINT " "; PRINT USING "##.##"; ch(i); : PRINT " "; PRINT USING "#.####"; mm(i); : PRINT " "; PRINT USING " #.###"; P(i); : PRINT " "; PRINT USING " ##.#"; beta(i) * rtd; : PRINT USING " ##.###"; T(i); : PRINT " "; PRINT USING " #.##"; mach(i) NEXT i INPUT vbv ' Save parameters to file OPEN drive$ + "propdata.dat" FOR OUTPUT AS #7 FOR i = 1 TO 25 PRINT #7, station(i), ch(i), tcGP(i), camber(i), P(i) PRINT #7, pitch(i), pitch(i), epitch(i), beta(i) * rtd, mach(i) PRINT #7, cgx(i), cgy(i), cg(i) NEXT i CLOSE #7 PRINT "data saved" INPUT vbv CLS GOSUB plotchord2 INPUT vbv ' Plot blade angle to chordline against advance WINDOW (-5, -2)-(27, 80): CLS LINE (0, 0)-(0, 10), 2 LINE (0, 0)-(25, 0), 2: COLOR 14 LOCATE 2, 2: PRINT "Radial distn of blade angle and advance" FOR i = BBB + 1 TO 25 LINE (i - 1, (beta(i - 1) - Az0) * rtd)-(i, (beta(i) - Az0) * rtd), 14 LINE (i - 1, advance(i - 1) * rtd)-(i, advance(i) * rtd), 12 NEXT i INPUT vbv STOP END 'ÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ ' SUBROUTINES AND DATA 'ÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ nacafoil: ' Select airfoils initial trial value for chord at 25 stations ' along propeller blade prop$ = "F2C" name$ = " CNCMOULD\F2C\ANALYSE\DESF2CGP.BAS" FOR i = 1 TO 25: ch(i) = 20: NEXT i ' initial trial values RETURN ' for chord..parallel chord atmos: ' ATMOSPHERICS RHO1 = 1.225 ' kg/m^3, air density, standard conditions Ps = 101.325 ' International standard pressure Ts = 15 ' " " temperature Rg = Ps / (RHO1 * (Ts + 273.16)) ' gas constant CpCv = 1.4: q = 1: sig = 1 PRINT " ? Standard Atmosphere , enter 1": INPUT q: PRINT q PRINT " Enter altitude (ft) for analysis": INPUT altitude IF altitude = 0 THEN altitude = 500 PRINT " "; altitude sig = (1 - 6.879E-06 * altitude) ^ 4.258 IF q = 1 THEN rho = RHO1 * sig IF q = 1 THEN Ta = 15 IF q = 1 THEN 6 PRINT " Enter temp. (C) for analysis": INPUT Ta IF Ta = 0 THEN Ta = 30'15 PRINT " "; Ta Pa1 = RHO1 * sig * Rg * (Ta + 273.16) ' Standard atmosphere adjust PRINT " Enter pressure (kPa) for analysis": INPUT Pa: PRINT Pa IF Pa = 0 THEN Pa = Pa1 PRINT Pa rho = Pa / (Rg * (Ta + 273.16)) 6 rhoratio = rho / RHO1 RETURN rotate2: betak = ATN((y2 - y1) / (x2 - x1)) x11 = x1 * COS(betak) + y1 * SIN(betak) x21 = x2 * COS(betak) + y2 * SIN(betak) y11 = -x1 * SIN(betak) + y1 * COS(betak) y21 = -x2 * SIN(betak) + y2 * COS(betak) y = 0 x = (x11 + x21) / 2 betak = -betak ' back transform x10 = x * COS(betak) + y * SIN(betak) y10 = -x * SIN(betak) + y * COS(betak) RETURN rotate: xu(pp) = x1 * COS(betan) + y1 * SIN(betan) xl(pp) = x2 * COS(betan) + y2 * SIN(betan) yu(pp) = -x1 * SIN(betan) + y1 * COS(betan) yl(pp) = -x2 * SIN(betan) + y2 * COS(betan) RETURN vectors: ' Plot vectors CLS scl = 1 x1 = -10 * scl: y1 = -10 * scl x2 = 140 * scl: y2 = 85 * scl del = 0: IF j > 15 THEN del = 120 IF j > 15 THEN x1 = (x1 + del) * scl IF j > 15 THEN x2 = (x2 + del) * scl WINDOW (x1, y1)-(x2, y2): scl = 3 FOR pp = 1 TO 100 LINE (-yurnext(pp - 1) * scl, zurnext(pp - 1) * scl)-(-yurnext(pp) * scl, zurnext(pp) * scl), 14 LINE (-ylrnext(pp - 1) * scl, zlrnext(pp - 1) * scl)-(-ylrnext(pp) * scl, zlrnext(pp) * scl), 14 NEXT pp LINE (0, 0)-(0, 60), 11 LINE (0, 0)-(220, 0), 11 x = Vr / TAN(advance(j)): y = Vr ' advance x4 = ch(10) * 3 * COS(beta(j)) ' zero-lift line y4 = ch(10) * 3 * SIN(beta(j)) ' " " ' x1 = Wx(ll) * COS(phi2(ll)) ' induced flow angle y1 = Wx(ll) * SIN(phi2(ll)) ' " " " LINE (0, 0)-(x4, y4), 13 ' zero-lift line LINE (0, 0)-(x, y), 4 ' advance LINE (x, 0)-(x, y), 2 ' Vr LINE (0, 0)-(x1, y1), 3 ' W LINE (x1, y1)-(x, y), 12 ' interference resultant LINE (x1, y1)-(x1, y1 - aa1 * Vr), 12 ' vector a = aa1 LINE (x, y)-(x1, y1 - aa1 * Vr), 12 ' vector a'= a11 LINE (5 + del, 79)-(15 + del, 79), 4: LOCATE 2, 15: COLOR 4: PRINT "advance" LINE (5 + del, 76)-(15 + del, 76), 2: LOCATE 3, 15: COLOR 2: PRINT "Vr" LINE (5 + del, 72)-(15 + del, 72), 13: LOCATE 4, 15: COLOR 13: PRINT "Beta" LINE (5 + del, 68)-(15 + del, 68), 3: LOCATE 5, 15: COLOR 3: PRINT "W inflow vel." LOCATE 2, 40: COLOR 4: PRINT "Iteration vectors" LOCATE 4, 30: COLOR 12: PRINT USING " Blade element J = ## Iteration = ## (##) "; j; Powerloop; Niterations CIRCLE (0, 0), 30, 13, 0, beta(j), .87 CIRCLE (0, 0), 32, 3, 0, ATN(y1 / x1), .87 CIRCLE (0, 0), 34, 4, 0, ATN(y / x), .87 RETURN chordplot: ' PLOT CHORD DISTRIBUTION x1 = -20: x2 = 35: y1 = -100: y2 = 348 WINDOW (x1, y1)-(x2, y2) FOR th = BBB + 1 TO 25 LINE (th - 1, ch(th - 1) / 3)-(th, ch(th) / 3) LINE (th - 1, -2 * ch(th - 1) / 3)-(th, -2 * ch(th) / 3) NEXT th RETURN slipplot: ' PLOT SLIP DISTRIBUTION x1 = -2: x2 = 27: y1 = -1: y2 = 25 WINDOW (x1, y1)-(x2, y2): CLS LINE (0, 0)-(0, 10): LINE (0, 0)-(25, 0) LOCATE 2, 5: PRINT "Radial distribution of slip " LOCATE 3, 5: PRINT " Constant slip implies minimum induced loss" LOCATE 8, 8: PRINT "Slip velocity" LOCATE 23, 60: PRINT "station" FOR th = BBB + 1 TO 25 LINE (th - 1, slip(th - 1))-(th, slip(th)) NEXT th GOSUB delay RETURN plotchord2: ' PLOT NEW CHORD DISTRIBUTION x1 = -10: x2 = 35: y1 = -50: y2 = 100 WINDOW (x1, y1)-(x2, y2): COLOR 4 LINE (-2, 0)-(2, 0): LINE (0, -8)-(0, 8): COLOR 2 FOR th = 2 TO 25 LINE (th - 1, ch(th - 1) * .42)-(th, ch(th) * .42) LINE (th - 1, -.58 * ch(th - 1))-(th, -.58 * ch(th)) NEXT th LINE (25, ch(25) * .42)-(25, -ch(25) * .58) GOSUB delay GOSUB delay RETURN delay: RETURN timex = .5 ' display time, seconds starttime = TIMER WHILE (TIMER - starttime) < timex WEND RETURN applyPG: ' Correct airfoil at station J for compressibilty ' Get thickness, camber, high point, zero-lift angle ' Also get airfoil C/G coordinates WINDOW (-35, -8)-(15, 28): CLS ' LINE (0, 0)-(-10, 0): LINE (0, 0)-(0, 10) ' ================================================================ ' Airfoils are selected for the design point below *************** ' GOSUB adder ' root thickening factor ' Set new blade parameters. Only use airfoil from GENERIC at Mach = 0 ' Determine Mach number, advance per revolution first, then: ' Stretch generic airfoil in direction of motion ' First rotate generic airfoil to angle attack AOA of chordline IF flag THEN 3 ELSE 4 3 CLS FOR pp = 1 TO 100 LINE (yusg(pp - 1), zusg(pp - 1))-(yusg(pp), zusg(pp)), 10 LINE (ylsg(pp - 1), zlsg(pp - 1))-(ylsg(pp), zlsg(pp)), 10 NEXT pp 4 LOCATE 1, 1 IF flag THEN INPUT "generic", vbv theta = AOA * dtr FOR pp = 0 TO 100 yurnext(pp) = yusg(pp) * COS(theta) + zusg(pp) * SIN(theta) zurnext(pp) = -yusg(pp) * SIN(theta) + zusg(pp) * COS(theta) ylrnext(pp) = ylsg(pp) * COS(theta) + zlsg(pp) * SIN(theta) zlrnext(pp) = -ylsg(pp) * SIN(theta) + zlsg(pp) * COS(theta) NEXT pp ' Plot rotated airfoil at angle of attack AA IF flag = 1 THEN 51 ELSE 61 51 FOR pp = 1 TO 100 LINE (yurnext(pp - 1), zurnext(pp - 1))-(yurnext(pp), zurnext(pp)), 14 LINE (ylrnext(pp - 1), zlrnext(pp - 1))-(ylrnext(pp), zlrnext(pp)), 14 NEXT pp LINE (0, 0)-(-15, 0), 12 61 ' now stretch along abcissa FOR pp = 0 TO 100 yurnext(pp) = yurnext(pp) * PGrule(j) ylrnext(pp) = ylrnext(pp) * PGrule(j) NEXT pp IF flag THEN 71 ELSE 81 71 FOR pp = 1 TO 100 LINE (yurnext(pp - 1), zurnext(pp - 1))-(yurnext(pp), zurnext(pp)), 14 LINE (ylrnext(pp - 1), zlrnext(pp - 1))-(ylrnext(pp), zlrnext(pp)), 14 NEXT pp 81 IF flag THEN INPUT vbv ' Get angle of attack AA(J) for mach(J) z1 = zurnext(0) - zurnext(100) z2 = yurnext(0) - yurnext(100) AA(j) = -ATN(z1 / z2) ' Rotate to zero angle before re-scaling theta = -AA(j) FOR pp = 0 TO 100 yurnew(pp) = yurnext(pp) * COS(theta) + zurnext(pp) * SIN(theta) zurnew(pp) = -yurnext(pp) * SIN(theta) + zurnext(pp) * COS(theta) ylrnew(pp) = ylrnext(pp) * COS(theta) + zlrnext(pp) * SIN(theta) zlrnew(pp) = -ylrnext(pp) * SIN(theta) + zlrnext(pp) * COS(theta) NEXT pp ' Plot zero-angle of attack AA(J) airfoil IF flag THEN 91 ELSE 101 91 FOR pp = 1 TO 100 LINE (yurnew(pp - 1), zurnew(pp - 1))-(yurnew(pp), zurnew(pp)), 14 LINE (ylrnew(pp - 1), zlrnew(pp - 1))-(ylrnew(pp), zlrnew(pp)), 14 NEXT pp LINE (0, 0)-(-15, 0), 12 101 ' re-scale to correct chord length normalisechord = ABS(yurnew(0) - yurnew(100)) scale = ch(j) / normalisechord ' scale = scale * widen FOR pp = 0 TO 100 yurnew(pp) = yurnew(pp) * scale ylrnew(pp) = ylrnew(pp) * scale zurnew(pp) = zurnew(pp) * scale zlrnew(pp) = zlrnew(pp) * scale NEXT pp ' Thicken up airfoil root ' Do this by adding a proportion of upper camber to lower camber! ' FOR pp = 0 TO 100 ' zlrnew(pp) = zlrnew(pp) - G(ii) * zurnew(pp) ' NEXT pp CLS IF flag THEN 111 ELSE 121 111 FOR pp = 1 TO 100 LINE (yurnew(pp - 1), zurnew(pp - 1))-(yurnew(pp), zurnew(pp)) LINE (ylrnew(pp - 1), zlrnew(pp - 1))-(ylrnew(pp), zlrnew(pp)) NEXT pp 121 ' Get t/c of G-P airfoil with airfoil at zero angle max = -10 FOR pp = 10 TO 100 IF zurnew(pp) > max THEN max = zurnew(pp) IF zurnew(pp + 1) < max THEN GOTO 241 NEXT pp 241 tcGP(j) = -(zurnew(pp) - zlrnew(pp)) / (yurnew(0) - yurnew(100)) T(j) = tcGP(j) ' Get camber hi-point location ' First get camberline FOR pp = 0 TO 100 xcamber(pp) = (yurnew(pp) + ylrnew(pp)) / 2 ycamber(pp) = (zurnew(pp) + zlrnew(pp)) / 2 NEXT pp FOR pp = 10 TO 100 IF ycamber(pp) < ycamber(pp - 1) THEN GOTO 231 NEXT pp 231 dum = ABS(xcamber(0)) - ABS(xcamber(100)) P(j) = (ABS(xcamber(0)) - ABS(xcamber(pp - 1))) / dum ' get camber camber(j) = ycamber(pp - 1) / dum mm(j) = camber(j) ' Get zero lift angle zerolift(j) = -ATN(camber(j) * 1 / (1 - P(j))) * .7 ' Rotate to blade angle ... increase from zero angle of chordline betax = advance(j) + AA(j) + zerolift(j) FOR pp = 0 TO 100 yurnext(pp) = yurnew(pp) * COS(betax) + zurnew(pp) * SIN(betax) zurnext(pp) = -yurnew(pp) * SIN(betax) + zurnew(pp) * COS(betax) ylrnext(pp) = ylrnew(pp) * COS(betax) + zlrnew(pp) * SIN(betax) zlrnext(pp) = -ylrnew(pp) * SIN(betax) + zlrnew(pp) * COS(betax) NEXT pp ' Plot airfoil for Mach(j) CLS FOR pp = 1 TO 100 LINE (yurnext(pp - 1), zurnext(pp - 1))-(yurnext(pp), zurnext(pp)) LINE (ylrnext(pp - 1), zlrnext(pp - 1))-(ylrnext(pp), zlrnext(pp)) NEXT pp 'INPUT "new airfoil", vbv ' Get area of airfoil using Heron's formula for area of triangles ' and C/G of airfoil GOSUB area RETURN area: ss = 0 FOR pp = 0 TO 99 STEP 1 a = SQR((ylrnext(pp) - ylrnext(pp + 1)) ^ 2 + (zlrnext(pp) - zlrnext(pp + 1)) ^ 2) b1 = SQR((ylrnext(pp) - yurnext(pp + 1)) ^ 2 + (zlrnext(pp) - zurnext(pp + 1)) ^ 2) c = SQR((yurnext(pp + 1) - ylrnext(pp + 1)) ^ 2 + (zurnext(pp + 1) - zlrnext(pp + 1)) ^ 2) ss = ss + 1 s = (a + b1 + c) / 2 area(ss) = SQR(s * (s - a) * (s - b1) * (s - c)) a = SQR((yurnext(pp + 1) - yurnext(pp + 2)) ^ 2 + (zurnext(pp + 1) - zurnext(pp + 2)) ^ 2) b1 = SQR((yurnext(pp + 1) - ylrnext(pp + 1)) ^ 2 + (zurnext(pp + 1) - zlrnext(pp + 1)) ^ 2) c = SQR((ylrnext(pp + 1) - yurnext(pp + 2)) ^ 2 + (zlrnext(pp + 1) - zurnext(pp + 2)) ^ 2) ss = ss + 1 s = (a + b1 + c) / 2 area(ss) = SQR(s * (s - a) * (s - b1) * (s - c)) NEXT pp sum = 0 FOR jj = 1 TO ss sum = sum + area(jj) NEXT jj areatotal = sum ' Find mid area index to get C/G location sum = 0 FOR jj = 1 TO ss sum = sum + area(jj) IF sum > areatotal / 2 THEN GOTO 576 NEXT jj 576 sum = sum - area(jj) cg(j) = (sum / areatotal) * 100 ' now get exact C/G coordinates a = (tcGP(j) * ch(j) / 2) * SIN(betax) ' x-displacement b2 = (tcGP(j) * ch(j) / 2) * COS(betax)' y-disp. cgx(j) = ylrnext(jj / 2) + a cgy(j) = zlrnext(jj / 2) + b2 LINE (cgx(j) - 1, cgy(j))-(cgx(j) + 1, cgy(j)) LINE (cgx(j), cgy(j) - 1)-(cgx(j), cgy(j) + 1) IF j = 3 THEN GOTO fgh ELSE RETURN fgh: cgx(1) = cgx(j): cgy(1) = cgy(j) cgx(2) = cgx(j): cgy(2) = cgy(j) RETURN adder: g(1) = 1: g(2) = .78: g(3) = .58: g(4) = .41: g(5) = .26 g(6) = .14: g(7) = .06: g(8) = .02 FOR j = 9 TO 25: g(j) = 0: NEXT j RETURN RETURN ' outer nest generic: ' Load generic airfoil data, station 18 of F2C23 ' Coordinates for Mach = 0 ' Set chordline to zero angle 'flag = 1 OPEN drive$ + "GENERIC.DAT" FOR INPUT AS #7 FOR pp = 0 TO 100 INPUT #7, yusg(pp), zusg(pp), ylsg(pp), zlsg(pp) NEXT pp FOR i = 1 TO 25 INPUT #7, dum, tchord(i) NEXT i INPUT #7, s, beta, angle, dum INPUT #7, dum, dum, dum, dum CLOSE #7 tchord = tchord(s) ' generic airfoil chord ' Clear data FOR i = 1 TO 25: tchord(i) = 0: NEXT i WINDOW (-25, -15)-(25, 15): CLS ' LINE (0, 0)-(-10, 0): LINE (0, 0)-(0, 10) ' Plot rotated generic airfoil at angle beta(18) IF flag = 1 THEN 1 ELSE 2 1 FOR pp = 1 TO 100 LINE (yusg(pp - 1), zusg(pp - 1))-(yusg(pp), zusg(pp)), 10 LINE (ylsg(pp - 1), zlsg(pp - 1))-(ylsg(pp), zlsg(pp)), 10 NEXT pp 2 theta = -beta * dtr FOR pp = 0 TO 100 yurnext(pp) = yusg(pp) * COS(theta) + zusg(pp) * SIN(theta) zurnext(pp) = -yusg(pp) * SIN(theta) + zusg(pp) * COS(theta) ylrnext(pp) = ylsg(pp) * COS(theta) + zlsg(pp) * SIN(theta) zlrnext(pp) = -ylsg(pp) * SIN(theta) + zlsg(pp) * COS(theta) NEXT pp FOR pp = 0 TO 100 yusg(pp) = yurnext(pp) zusg(pp) = zurnext(pp) ylsg(pp) = ylrnext(pp) zlsg(pp) = zlrnext(pp) NEXT pp ' Plot rotated airfoil at angle of attack zero to chordline ' CLS IF flag = 1 THEN 511 ELSE 611 511 FOR pp = 1 TO 100 LINE (yusg(pp - 1), zusg(pp - 1))-(yusg(pp), zusg(pp)), 14 LINE (ylsg(pp - 1), zlsg(pp - 1))-(ylsg(pp), zlsg(pp)), 14 NEXT pp LINE (0, 0)-(-15, 0), 12 611 ' Get t/c of P-G airfoil with airfoil at zero angle max = -10 FOR pp = 10 TO 100 IF zusg(pp) > max THEN max = zusg(pp) IF zusg(pp + 1) < max THEN GOTO 2411 NEXT pp 2411 tcGPz = -(zusg(pp) - zlsg(pp)) / (yusg(0) - yusg(100)) ' Get camber hi-point location ' First get camberline FOR pp = 0 TO 100 xcamber(pp) = (yusg(pp) + ylsg(pp)) / 2 ycamber(pp) = (zusg(pp) + zlsg(pp)) / 2 NEXT pp FOR pp = 10 TO 100 IF ycamber(pp) < ycamber(pp - 1) THEN GOTO 2311 NEXT pp 2311 dum = ABS(xcamber(0)) - ABS(xcamber(100)) Pz = (ABS(xcamber(0)) - ABS(xcamber(pp - 1))) / dum ' get camber camberz = ycamber(pp - 1) / dum ' Get zero lift angle zeroliftz = -ATN(camberz * 1 / (1 - Pz)) RETURN printdata: RETURN OPEN drive$ + "design.dat" FOR OUTPUT AS #11 PRINT #11, "I Cl Cl/Cd PITCH ePITCH Ch BETA Advance zero-angle AOA Mach" PRINT #11, " (ins) (ins) mm deg deg deg deg " FOR i = BBB TO 25 STEP 2 COLOR 12: PRINT #11, USING "##"; i; PRINT #11, USING " #.##"; Cl1(i); : PRINT #11, " "; PRINT #11, USING "##.##"; Cl1(i) / Cd1(i); : PRINT #11, " "; PRINT #11, USING "##.#"; pitch(i); : PRINT #11, " "; PRINT #11, USING " ##.#"; epitch(i); : PRINT #11, " "; PRINT #11, USING " ###.#"; ch(i); : PRINT #11, " "; PRINT #11, USING " ##.#"; beta(i) * rtd; : PRINT #11, " "; PRINT #11, USING " ##.#"; advance(i) * rtd; : PRINT #11, " "; PRINT #11, USING " ####.#"; zerolift(i) * rtd; PRINT #11, USING " ####.#"; (AA(i) + zerolift(i)) * rtd; PRINT #11, USING " #.##"; mach(i) NEXT i CLOSE #11 RETURN saveairfoil: ' Store airfoil in drive$ + "airfoil.dat" FOR pp = 0 TO 100 yua(jdum, pp) = yurnext(pp): zua(jdum, pp) = zurnext(pp) yla(jdum, pp) = ylrnext(pp): zla(jdum, pp) = zlrnext(pp) PRINT #4, yua(jdum, pp), zua(jdum, pp), yla(jdum, pp), zla(jdum, pp) NEXT pp RETURN SUB prompt (x) ' The subprogram PROMPT displays "Press any key to continue" and loops until ' a key is pressed. X specifies the row on which the text appears. LOCATE 23, 15: c = 15 DO COLOR c: LOCATE x, 65 PRINT "Press any key"; FOR T = 1 TO 1000: NEXT IF c = 15 THEN c = 2 ELSE c = 15 LOOP UNTIL INKEY$ <> "" CLS : COLOR 15 END SUB