' PROGRAM NAME: DESIGN.BAS (C:\CNCMOULD\F2C\ANALYSE\DESIGN.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 ' ************************************************************* ' 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, compressibility ' effects of the formation of shock waves. ' 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. ' 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. ' Arrays use din 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) DIM xu(100), xl(100), yu(100), yl(100) DIM xur(100), yur(100), xlr(100), ylr(100) DIM yt(100), xord(100) SCREEN 9, , 1, 1: pi = 4 * ATN(1): rtd = 180 / pi: dtr = pi / 180 ' rtd converts radians to degrees, and dtr conversely. '--------------------------------------------------------------------- ' INPUT VARIABLES spin = 1.2 ' spinner diameter in inches P = .8: Pdes = P ' sea level Power (BHP) at 15 C, 101.325 kPa RPM = 26000: RPMdes = RPM ' RPM vV = 130: VVVdes = vV ' Flight speed of aircraft (MPH) Ddes = 6 ' Diameter of prop in inches R = Ddes / 2 ' Propeller tip radius (inches) R = R * .0254 ' " " " (metres) B = 2 ' Number of propeller blades ' Calculate initial solidity.. amount of the prop disc shadowed by blades FOR i = 1 TO 25 ' there are 25 equally spaced zeta = i / 25 - .02 ' blade elements ch(i) = 20 ' width of initial blade (mm) z = 2 * pi * R * zeta sol(i) = B * ch(i) * .001 / z ' solidity NEXT i Niterations = 6 ' 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 (not used) bbb = INT(spin / 2 / R * 25) ' Working part of blade, elements ' outside spinner ' ****************************************************************** ' Select airfoil data GOSUB nacafoil ' Choose camber, thickness and camberhigh-point ' ------------------------------------------------------------------ ' UNIT CONVERSIONS TO MKS spin = spin * .0254 ' spinner diameter, metres P = P * 746 ' BHP to WATTS om = RPM * 2 * pi / 60 ' RPM to radians/sec v = vV * 2.54 * 12 * 5280 / 60 / 60 / 100 ' airspeed, MPH to m/sec V bbb = 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 ' POWER MATCHING LOOP PPPq = Niterations ' Number of iterations on slip and power absorption FOR PPP = 1 TO PPPq ' PPP is iteration counter RPM = RPMdes: vV = VVVdes: COLOR 4 PRINT " "; prop$: PRINT : COLOR 2 PRINT " RPM = ", RPM PRINT " MPH = ", vV: COLOR 9: PRINT PRINT " "; PPP, PPPq: 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 ' set window for graphics plots GOSUB plotnaca CLS ' ------------------------------------------------------------------ ' 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 ' ------------------------------------------------------------------ ' Find zero lift angle Az0 at station J Az0 = ATN(M2(J) / (1 - P.p(J))) ' zero-lift angle relative to chordline ' 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) ' Now consider advance ratio at current blade station ' Set initial values 323 ll = 1 ' counter AOA = 7 ' chosen angle of attack (degrees from ' zero-lift to advance) alpha(ll) = .5 * dtr ' start just above zero lift beta(J) = advance(J) + AOA * dtr ' 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 ' First guess (iteration) Cl1 = .1 * alpha(ll) * rtd ' value .1 is lift slope Cd1 = .02 ' for infinite aspect ratio lift-slope 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 from airfoil z = sol(J) * Cy1 / (4 * F(J) * Sp ^ 2) aa1 = z / (1 - z) ' swirl interference .. tangemtial downwash 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 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. NEXT J ' continue to next blade element '------------------------------------------------------------ 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 IF PPP = PPPq 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 sol(i) = B * ch(i) * .001 / z ' new solidity 543 NEXT i GOSUB plotchord2 NEXT PPP ' End power loop, display results ' ---------------------------------------------------------- 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) - Az0) ' pitch to chordline fpitch(i) = z / .0254 * (TAN(beta(i) - Az0) - faceangle(i)) ' face epitch(i) = z * TAN(beta(i)) / .0254 ' experimental pitch NEXT i ' (measured to xero-lift) CLS : LOCATE 1, 1 COLOR 12: BEEP: PRINT : PRINT : PRINT " ", prop$; PRINT " DESIGN.BAS (no correction for compressibilty)" PRINT : COLOR 2 PRINT " RPM = "; : PRINT USING "#####"; RPMdes; PRINT USING " V = ### MPH"; VVVdes; PRINT USING " Diameter = ###.## inches"; Ddes PRINT " Blade number = "; : PRINT USING "##"; B; PRINT USING " No. of iterations ###"; Niterations PRINT PRINT " POWER (HP) = "; : PRINT USING "###.#"; P PRINT " THRUST (lbs) = "; : PRINT USING "###.#"; T / 4.45 PRINT "THRUST POWER (hp) = "; : PRINT USING "###.#"; P * n PRINT " Efficiency = "; : PRINT USING "####.#"; n * 100 PRINT GOSUB delay COLOR 10 PRINT "I Cl Cl/Cd PITCH ePITCH Ch BETA Advance zero-angle AOA " PRINT " (ins) (ins) mm deg deg deg deg " FOR i = bbb 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) - Az0) * rtd; : PRINT " "; PRINT USING " ##.#"; advance(i) * rtd; : PRINT " "; PRINT USING " ####.#"; -Az0 * rtd; PRINT USING " ####.#"; AOA - Az0 * rtd NEXT i COLOR 0 INPUT vbv CLS GOSUB plotchord2 COLOR 12 LOCATE 2, 2: PRINT " Chord radial distribution adjusted to desired power input" LOCATE 3, 2: PRINT " and also adjusted to give constant slip" COLOR 0 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 for 25 stations along propeller blade prop$ = "F2C" name$ = " CNCMOULD\F2C\ANALYSE\DESF2C.BAS" ' t=thickness-chord ratio ' m=fractional camber ' p=fractional high point from L.E. FOR i = 1 TO 25 T(i) = .1 M2(i) = .04 P.p(i) = .4 ' camber high point Re(i) = 180000 ' guess, Reynolds number NEXT i RETURN atmos: ' ATMOSPHERICS RHO1 = 1.225 ' kg/m^3 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 = 0 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 = 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 plotnaca: te(J) = ch(J) / 60 ' thicken trailing edge ' Thickness form distribution NACA 4-digit FOR pp = 0 TO 100 xord(pp) = pp / 100 z = xord(pp) yt(pp) = T(J) / .2 * (.2969 * z ^ .5 - .126 * z - .3516 * z ^ 2 + .2843 * z ^ 3 - .1015 * z ^ 4) NEXT pp ' modify forward part of thickness form to circle startn = 1 stopn = 30 x1 = xord(startn): y1 = yt(startn) x2 = xord(stopn): y2 = yt(stopn) GOSUB rotate2 ax1 = x10: ay1 = y10 am1 = TAN(-betak + pi / 2) b1 = ay1 - am1 * ax1 x1 = xord(startn): y1 = yt(startn) x2 = xord(stopn + 1): y2 = yt(stopn + 1) GOSUB rotate2 ax2 = x10: ay2 = y10 am2 = TAN(-betak + pi / 2) b2 = ay2 - am2 * ax2 x1 = ax1: y1 = ay1: m1 = am1 x2 = ax2: y2 = ay2: M2 = am2 x = (b2 - b1) / (m1 - M2) y = (M2 * b1 - m1 * b2) / (M2 - m1) Rx = SQR((x - xord(startn)) ^ 2 + (y - yt(startn)) ^ 2) FOR pp = startn TO stopn yt(pp) = SQR(Rx ^ 2 - (xord(pp) - x) ^ 2) + y NEXT pp ' Correct equation at trailing edge .. Abbott etc wrong shift = yt(100) / 100 FOR pp = 0 TO 100 yt(pp) = (yt(pp) - shift * pp) NEXT pp ' Correct thickness to chord ratio to true value "t" .. wrong again max = 0 ' find max thickness normal to camber line FOR pp = 1 TO 99 IF yt(pp) > yt(pp + 1) AND yt(pp) >= yt(pp - 1) THEN max = yt(pp) NEXT pp FOR pp = 0 TO 100 ' correct to desired thickness t yt(pp) = yt(pp) * T(J) / (2 * max) NEXT pp ' CAMBER MEAN LINE FOR NACA 4-DIGIT AIRFOIL FOR pp = 1 TO 100 IF xord(pp) > P.p(J) THEN 55 NEXT pp 55 zs = pp - 1 ' high point location on camber line FOR pp = 0 TO 100 z = xord(pp) yc(pp) = M2(J) / P.p(J) ^ 2 * (2 * P.p(J) * z - z ^ 2) IF pp > zs + 1 THEN yc(pp) = M2(J) / (1 - P.p(J)) ^ 2 * (1 - 2 * P.p(J) + 2 * P.p(J) * z - z ^ 2) NEXT pp ' DETERMINE forward ORDINATES NORMAL TO MEAN LINE FOR pp = 0 TO zs z = xord(pp) th = ATN(2 * M2(J) / P.p(J) * (1 - z / P.p(J))) xu(pp) = z - yt(pp) * SIN(th) yu(pp) = yc(pp) + yt(pp) * COS(th) xl(pp) = z + yt(pp) * SIN(th) yl(pp) = yc(pp) - yt(pp) * COS(th) NEXT pp ' DETERMINE rearward ORDINATES NORMAL TO MEAN LINE FOR pp = zs + 1 TO 100 z = xord(pp) th = ATN(2 * M2(J) / (1 - P.p(J)) ^ 2 * (P.p(J) - z)) xu(pp) = z - yt(pp) * SIN(th) yu(pp) = yc(pp) + yt(pp) * COS(th) xl(pp) = z + yt(pp) * SIN(th) yl(pp) = yc(pp) - yt(pp) * COS(th) NEXT pp ' Scale ordinates to true chord in mm FOR pp = 0 TO 100 xu(pp) = xu(pp) * ch(J) yu(pp) = yu(pp) * ch(J) xl(pp) = xl(pp) * ch(J) yl(pp) = yl(pp) * ch(J) xord(pp) = xord(pp) * ch(J) yc(pp) = yc(pp) * ch(J) NEXT pp ' Rotate upper surface to thicken TE for flicking FOR pp = 0 TO 100 theta = ATN(.5 * (-te(J)) / (xu(100) - xu(0))) xur(pp) = xu(pp) * COS(theta) + yu(pp) * SIN(theta) yur(pp) = -xu(pp) * SIN(theta) + yu(pp) * COS(theta) xu(pp) = xur(pp) yu(pp) = yur(pp) NEXT pp ' Rotate lower surface to thicken TE for flicking FOR pp = 0 TO 100 theta = ATN(.5 * te(J) / (xl(100) - xl(0))) xlr(pp) = xl(pp) * COS(theta) + yl(pp) * SIN(theta) ylr(pp) = -xl(pp) * SIN(theta) + yl(pp) * COS(theta) xl(pp) = xlr(pp) yl(pp) = ylr(pp) NEXT pp dum = xu(100): dum1 = xl(100) FOR pp = 0 TO 100 xu(pp) = ABS(xu(pp) - dum) xl(pp) = ABS(xl(pp) - dum1) x1 = xu(pp): y1 = yu(pp) x2 = xl(pp): y2 = yl(pp) GOSUB rotate NEXT pp scl = 1 x1 = -10 * scl: y1 = -10 * scl x2 = 140 * scl: y2 = 93 * scl IF J > 15 THEN x1 = 110 * scl IF J > 15 THEN x2 = 220 * scl WINDOW (x1, y1)-(x2, y2) FOR pp = 1 TO 100 LINE (xu(pp - 1), yu(pp - 1))-(xu(pp), yu(pp)), 14 LINE (xl(pp - 1), yl(pp - 1))-(xl(pp), yl(pp)), 14 NEXT pp LINE (0, 0)-(240, 0), 10: LINE (0, 0)-(0, 100), 10 'LINE (0, 0)-(50, 0), 14: LINE (0, 0)-(0, 50), 14 'INPUT vbv 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 betan = -(beta(J) - Az0): GOSUB plotnaca x = Vr / TAN(advance(J)): y = Vr ' advance x4 = ch(J) * 1.8: y4 = x4 * TAN(beta(J)) ' zero-lift line x4 = 40: y4 = x4 * TAN(beta(J)) x1 = Wx(ll) * COS(phi2(ll)): 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, 86)-(15, 86), 4: LOCATE 2, 15: COLOR 4: PRINT "advance" LINE (5, 82)-(15, 82), 2: LOCATE 3, 15: COLOR 2: PRINT "Vr" LINE (5, 78)-(15, 78), 13: LOCATE 4, 15: COLOR 13: PRINT "Beta" LINE (5, 74)-(15, 74), 3: LOCATE 5, 15: COLOR 3: PRINT "W" LOCATE 2, 40: PRINT "Iteration vectors" LOCATE 4, 30: PRINT USING " Blade element J = ## Iteration = ##"; J; PPP CIRCLE (0, 0), 30, 13, 0, beta(J), .78 CIRCLE (0, 0), 32, 3, 0, ATN(y1 / x1), .78 CIRCLE (0, 0), 34, 4, 0, ATN(y / x), .78 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 = 15 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 ': ' PLOT NEW CHORD DISTRIBUTION x1 = -10: x2 = 35: y1 = -100: y2 = 348 WINDOW (x1, y1)-(x2, y2) COLOR 12 FOR th = bbb + 1 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 GOSUB delay RETURN delay: timex = .5 starttime = TIMER WHILE (TIMER - starttime) < timex WEND 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 = bbb + 1 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 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