' WIND TURBINE DESIGN CODE....by Joe Supercool September 2004 ' Code is Quick Basic. Run time of order 10 minutes ' Estimation of turbine characteristics . ' Code determines blade angle and radial chord distribution ' for minimum induced loss. Lightly loaded condition applies. ' Airfoil lift coefficient is set to 0.5, drag to .03 ' Airfoils may conveniently be NACA 4-digit sections ' Air flow vectors joined tip to tail for display purposes. ' Downwash vector drawn normal to net inflow of V and OMr ' Froude efficiency from "Wind Turbine Engineering Design" by ' David M. Eggleston and Forrest S. Stoddard, Von Nostrand Reinhold ' 1987 ISBN 0-442-22195-9 ' For a wind turbine, wind speed and rotor RPM being known/anticipated, ' it is possible to draw a vector diagram of the airflows at a ' a blade element. Missing is the downwash vector, because the operating ' airfoil angle-of-attack is unknown. However, the angle-of-attack is ' constrained by its relationship with the lift coefficient. By applying ' the theory given by Larrabee in the 1977 National Free-flight Symposium ' report "Analytic Design of Propellers Having Minimum Induced Loss", it ' is possible to compute the air motion with respect (wrt) to the blade. ' The values so obtained have an apparent agreement with observation; in ' particular, the rotor power absorption appears to match known available ' wind power. ' The independent blade element assumption is made: ie, the vorticity of ' adjacent blade elements is ignored. DIM epitch(25), chord(25), pitch(25) DIM phi(25), alpha(25), beta(25), alpha0(25) DIM Cl(25), Cd(25) DIM advance(25), Vtorq(25) DIM aa(25), aad(25) DIM nind(25), nprof(25), nn(25), Nf(25) DIM dCt(25), dCp(25) DIM xxx(50), u(50), yc(50), xu(14), yu(14), yl(14), xl(14) DIM T(25), P(25), M(25), vdV(25), station(25) SCREEN 9, , 1, 1: pi = 4 * ATN(1): rtd = 180 / pi: dtr = 1 / rtd exp1 = 2.7182818# CLS PRINT " Rotor analysis code, graphical presentation" PRINT : PRINT COLOR 14 PRINT " if nothing is happening, press Enter" PRINT PRINT " to terminate, hit Ctrl/Pause" PRINT : COLOR 2 INPUT vbv CLS PRINT " Code by Supercool December 17, 2003" start: PRINT " Start design for new blade": INPUT vbv ' Data inputs V = 20 ' wind-speed in m/s RPM = 1500 ' Revs per minute RPS = RPM / 60 ' Revs per second R = 1 ' rotor radius in metres rho = 1.23 ' air density, kg/m^3 B = 2 ' number of blades chord75 = .1 ' blade width at 75% radius, in metres ' chord 75 will be the blade width at 75% ' after radial chord distribution optimisation startt = 1 ' start getting blade element values, outside spinner nloops = 3 ' number of loops for chord optimisation INPUT "Input wind speed, m/s ", V: IF V = 0 THEN V = 20 INPUT "Input rotor RPM ", RPM: IF RPM = 0 THEN RPM = 1500 INPUT "Input rotor radius, m ", R: IF R = 0 THEN R = 1 INPUT "Input number of blades ", B: IF B = 0 THEN B = 2 INPUT " Blade chord at 75% radius, m", chord75: IF chord75 = 0 THEN chord75 = .1 PRINT : COLOR 14 PRINT "Wind speed, m/s ", V PRINT "Rotor RPM ", RPM: RPS = RPM / 60 PRINT "Rotor radius, m ", R PRINT "Number of blades ", B PRINT "Blade chord,75% radius, m", chord75 PRINT : COLOR 12: INPUT " hit return to continue"; vbv GOSUB turbinedata ' propeller initial geometry and airfoils ' Initialise displacement velocity ' "displacement velocity" is the axial component of the flow velocity ' downwind of the rotor. You could call it the slipstream velocity. FOR pp = 1 TO 25: vdV(pp) = 1: NEXT pp start1: FOR qqq = 1 TO nloops ' Iteration loop for finding ' best radial chord distribution ' adjust rotor chord to make radial distribution ' of displacement velocity constant...ie. get uniform slipstream velocity ' Reference point on the rotor is the chord at 75% radius, which thereby ' remains unchanged FOR i = startt TO 25: chord(i) = chord(i) * vdV(18) / vdV(i): NEXT i ' Power avaiable from wind L.L.Freris Wind enery Conversion Systems ' Prentice Hall ISBN 0-13-960527-4: 1990 p57 Actuator disc Ch4 D.J.Sharpe ad = pi * R ^ 2 ' area of rotor disc actualpower = .593 * .5 * rho * V ^ 3 * ad ' power available from wind ' prop speed in plane of rotation: m/s Vtorq = 2 * pi * R * RPS ' tip speed in torque plane OM = 2 * pi * RPS ' shaft speed, radians per second OMr = OM * R lam = V / Vtorq ' TAN of tip blade advance angle ' ****************************************************8 FOR i = startt TO 25 ' there are 25 blade elements in this model Nf.old = 0 ' old value of Froude efficiency Nf xzeta = i / 25 - .02 ' r/R, mid blade element fo = B / 2 * SQR(1 + lam ^ 2) / lam * (1 - xzeta) F = 2 / pi * ATN(SQR(exp1 ^ (2 * fo) - 1)) ' Prandtls velocity ' averaging fraction OMrr = 2 * pi * R * xzeta * RPS ' torque plane velocity ar r x = OMrr / V G = F * x ^ 2 / (1 + x ^ 2) ' Goldstein factor ' Get actual airflow angle .. blade motion wrt to air ' loop to apply lift-slope constraint to alpha, Cl ' downwash must be normal to inflow W Vtorq(i) = 2 * pi * xzeta * R * RPS ' tangential velocity of element i ' in the torq plane Zqq = SQR(V ^ 2 + Vtorq(i) ^ 2) ' resultant of V and Vtorq, Zqq is ' freestream velocity relative to blade thetad = ATN(V / Vtorq(i)) ' angle between V and Vtorq ' Get limiting values for alpha, the angle of attack to chordline alphamin = alpha0(i) ' lowest possible value of angle-of-attack ' Froude loop begins...depends on beta(i) only ' We are trying to get beta from the forces acting on a blade element ' The interference factors are known from the geometry. ' Try a range of beta's (blade angles) FOR beta = .01 * dtr TO pi / 2 STEP .1 * dtr ' beta is blade angle in radians beta(i) = beta ' For each trial beta, try a range of angle-of-attacks, starting at zero-lift angle FOR j = alpha0(i) + .0001 TO pi / 2 STEP .25 * dtr alpha(i) = j ' angle of attack of blade element to airflow phi(i) = beta + alpha(i) ' angle between plane of rotation and inflow W ' Get downwash velocity vector gamma = -phi(i) + thetad ' downwash angle ax = Zqq * SIN(gamma) ' downwash velocity downwash = ax ' drag = 0' .01 * Zqq ' vestigial code ' downwash = ax * (Zqq - drag) / Zqq ' downwash, viscous corrected ' W = Zqq * COS(gamma) ' inflow velocity..Liebeck and Adkins ' Get interference factors a' and a (aad and aa) from vector diagram epsilon = thetad - gamma aad(i) = downwash * SIN(epsilon) / Vtorq(i) aa(i) = downwash * COS(epsilon) / V ' Inflow Wair, as seen by the blade: ie freestream deflected by the downwash. Wair = SQR(((V * (1 - aa(i))) ^ 2) + ((Vtorq(i) * (1 + aad(i))) ^ 2)) vd = 2 * downwash / COS(phi(i)) ' vd is displacement velocity ' vd = 2 * V * (aa(i) + aad(i)) ' Larrabee ' equation from Larrabee, p93, fig.4 vdV(i) = vd / V ' fractional displacement velocity ' note: vd/V is constant radially for minimum induced loss ' light loading approximation applies here ' NB: I have a problem finding the correct equation for vd. The above ' equation vd = 2 * downwash / COS(phi(i)) is from Larrabee, p93, fig.4 ' vector diagram, which defines v' as the velocity far downstream ' after slipstream contraction. In turbines, the slipstream expands! ' Further, adding Larrabee's equations 7 an 8 yields ' vd = 2 * V * (aa(i) + aad(i)), which is not the same as above. ' Since this code is based on the vector diagram Fig. 4, I have ' chosen that. ' ' Perhaps for the turbine I could choose v' in the developed slipstream, ' in which case the velocity at the blade would be 2v' !!! ' Hey, why not? Then: ' vd = .5 * downwash / COS(phi(i)) ' vd is displacement velocity ' for expanding slipstream ! ' vdV(i) = vd / V ' fractional displacement velocity ' Liebeck and Adkins .... exact lift coefficient Cl(i) = 4 * pi * vd * G / (chord(i) * B * OM * (1 - aa(i))) * SIN(phi(i)) ' Larrabee light loading .... not the same as Liebeck and Adkins ' Cl(i) = 4 * pi * vd * G / (chord(i) * B * OM * SQR(1 + x ^ 2)) ' Get Cl from lift slope from thin airfoil theory: alpha is ' consistent with Cl Cldummy = 2 * pi * (alpha(i) - alpha0(i)) IF Cldummy > Cl(i) THEN GOTO 3 NEXT j 3 ': INPUT vbv GOSUB vectors IF Cl(i) > .5 THEN GOTO beta1 ' Power, thrust coeffs Cl(i) = -Cl(i) ' change sign to match propeller equations ' doubtful logic in this Cx = Cl(i) * SIN(phi(i)) + Cd(i) * COS(phi(i)) Cy = Cl(i) * COS(phi(i)) - Cd(i) * SIN(phi(i)) solid = B * chord(i) / (2 * pi * xzeta * R) ' solidity ' thrust and power coefficients ' note sign change for aad(i), opposite to propeller dCt(i) = pi ^ 3 / 4 * ((1 - aad(i)) / COS(phi(i))) ^ 2 * xzeta ^ 3 * solid * Cy dCp(i) = pi ^ 4 / 4 * ((1 - aad(i)) / COS(phi(i))) ^ 2 * xzeta ^ 4 * solid * Cx ' ************************************************************** ' Get Froude efficiency ee = -ATN(Cd(i) / Cl(i)) ' change sign back to turbine theory tp = TAN(phi(i)) a = aa(i) lamm = a * SIN(phi(i) / ((1 - a) * (1 / tp + ee))) ad = a * (tp - ee) dum = ((1 - a) * (1 / tp + ee) - a * (tp - ee)) ad = ad / dum x = (1 - a) / ((1 + ad) * tp) Cp = 4 * x * a * (1 - a) * (tp - ee) * (1 - ee * tp) Nf(i) = (27 / 16) * Cp ' Froude efficiency Nf.old = Nf(i) IF Cl(i) <= .5 THEN GOTO tt beta1: NEXT beta tt: NEXT i GOSUB integrate ' Get power absorption, thrust and efficiency NEXT qqq COLOR 2 PRINT " Iteration "; qqq: INPUT vbv PRINT " I chord alpha beta Cl Cd a a' Nf v'/V " COLOR 14 FOR i = startt TO 25 STEP 2 PRINT USING " ## ###.# ##.# ##.# ##.## ##.## #.## #.### #.## #.##"; i; chord(i) * 1000; alpha(i) * rtd; beta(i) * rtd; Cl(i); Cd(i); aa(i); aad(i); Nf(i); (aad(i) + aa(i)) * 2 NEXT i PRINT PRINT : COLOR 4 PRINT " For maximum induced efficiency propeller,"; : COLOR 2: PRINT " v'/V"; : COLOR 4: PRINT " is radially constant" PRINT " Nf is blade element Froude efficiency...use to optimise chord, blade amgle" COLOR 10: INPUT "Display rotor?"; vbv x1 = -.1 * R: x2 = 1.1 * R y1 = -.5 * R * .7: y2 = .5 * R * .7 WINDOW (x1, y1)-(x2, y2): CLS LINE (0, 0)-(R, 0) FOR i = 1 TO 25 station(i) = R / 25 * i NEXT i FOR i = startt TO 25 LINE (station(i - 1), .3 * chord(i - 1))-(station(i), .3 * chord(i)) LINE (station(i - 1), -.7 * chord(i - 1))-(station(i), -.7 * chord(i)) NEXT i: LINE (station(25), -.7 * chord(25))-(station(25), .3 * chord(25)) INPUT vbv GOTO start END rotate: betat = ATN((ya - yb) / (xa - xb)) xq = xc * COS(betat) + yc * SIN(betat) yq = -xc * SIN(betat) + yc * COS(betat) yq = 0 betat = -betat ' back transform xx = xq * COS(betat) + yq * SIN(betat) yx = -xq * SIN(betat) + yq * COS(betat) RETURN arrow: ax1 = -V * .05: ay1 = OMr * .003 ax2 = 0: ay2 = 0 ax3 = ax1: ay3 = -ay1 x11 = ax1 * COS(angle) + ay1 * SIN(angle) y11 = -ax1 * SIN(angle) + ay1 * COS(angle) x12 = ax3 * COS(angle) + ay3 * SIN(angle) y12 = -ax3 * SIN(angle) + ay3 * COS(angle) ax1 = x11 + tipx: ay1 = y11 + tipy ax2 = tipx: ay2 = tipy ax3 = x12 + tipx: ay3 = y12 + tipy LINE (ax1, ay1)-(ax2, ay2) LINE (ax3, ay3)-(ax2, ay2) LINE (ax1, ay1)-(ax3, ay3) RETURN integrate: dxzeta = 1 / 25 ' differential sumt = 0: sump = 0 FOR kk = startt TO 25 sumt = sumt + dCt(kk) ' integrate thrust sump = sump + dCp(kk) ' " power NEXT kk rho = 1.23 T = rho * RPS ^ 2 * (2 * R) ^ 4 * sumt * dxzeta ' thrust, Newtons P = rho * RPS ^ 3 * (2 * R) ^ 5 * sump * dxzeta ' power, Watts PRINT : COLOR 4: CLS PRINT " Iteration "; qqq PRINT USING " Thrust = ####.## lbs Power at shaft = ####.# HP "; T / 4.45; P / 746 PRINT USING " Actual power (recoverable from wind) = #####.# HP "; actualpower / 746 PRINT USING " Power coeff (conversion efficiency) = ###.### "; -P / actualpower PRINT RETURN naca: ' AIRFOIL COORDINATES FOR NACA 4-DIGIT FAMILY T = T(i) M = M(i) P = P(i) alpha0 = ATN(M / (1 - P)) ' t=thickness-chord ratio ' m=fractional camber ' p=fractional high point distance from L.E. ' DATA - CHORD FRACTION xx xxx(1) = 0: xxx(2) = .025: xxx(3) = .05: xxx(4) = .1: xxx(5) = .2: xxx(6) = .3 xxx(7) = .4: xxx(8) = .5: xxx(9) = .6: xxx(10) = .7: xxx(11) = .8: xxx(12) = .9 xxx(13) = .95: xxx(14) = 1 ' THICKNESS DISTRIBUTION FOR pp = 1 TO 14 Z = xxx(pp) u(pp) = T / .2 * (.2969 * Z ^ .5 - .126 * Z - .3516 * Z ^ 2 + .2843 * Z ^ 3 - .1015 * Z ^ 4) NEXT pp ' CAMBER MEAN LINE zs = INT(10 * P) + 3 ' Forward of high point FOR pp = 1 TO zs yc(pp) = M / P ^ 2 * (2 * P * xxx(pp) - xxx(pp) ^ 2) NEXT pp ' Rearward of high point FOR pp = zs + 1 TO 14 yc(pp) = M / (1 - P) ^ 2 * (1 - 2 * P + 2 * P * xxx(pp) - xxx(pp) ^ 2) NEXT pp FOR ik = 1 TO 14 xu(ik) = (xxx(14) - xxx(ik)) yu(ik) = (-yc(ik) + u(ik)) yl(ik) = (-yc(ik) - u(ik)) NEXT ik ' Scale to chord(18) dum = chord(i) / chord75 FOR ik = 1 TO 14 xu(ik) = xu(ik) * dum yu(ik) = yu(ik) * dum yl(ik) = yl(ik) * dum 'PRINT ik, xu(ik): INPUT vbv NEXT ik ' Rotate thru "angle" FOR ik = 1 TO 14 xudum = xu(ik) * COS(angle) + yu(ik) * SIN(angle) yudum = -xu(ik) * SIN(angle) + yu(ik) * COS(angle) xldum = xu(ik) * COS(angle) + yl(ik) * SIN(angle) yldum = -xu(ik) * SIN(angle) + yl(ik) * COS(angle) xu(ik) = xudum xl(ik) = xldum yu(ik) = yudum yl(ik) = yldum NEXT ik x0 = xu(1) * COS(alpha0) + yu(1) * SIN(alpha0) y0 = -xu(1) * SIN(alpha0) + yu(1) * COS(alpha(0)) RETURN ' 컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴� vectors: ' ***********************************************8 LOCATE 1, 1: COLOR 4: PRINT USING " rotational motion vector Omr ###.#"; Vtorq(i) LOCATE 2, 1: COLOR 5: PRINT USING " wind motion vector V ##.#"; V LOCATE 3, 1: COLOR 6: PRINT USING " interference factor, thrust a = #.### ###.#"; -aa(i); -aa(i) * V LOCATE 4, 1: COLOR 7: PRINT USING " interference factor, torque a' = #.### ###.#"; aad(i); aad(i) * Vtorq(i) LOCATE 5, 1: COLOR 8: PRINT USING " downwash vector ##.#"; downwash LOCATE 6, 1: COLOR 9: PRINT USING " blade motion vector wrt air W ###.#"; Wair LOCATE 13, 1: COLOR 6: PRINT USING " Zero-lift = ###.#"; alpha0(i) * rtd LOCATE 14, 1: COLOR 6: PRINT USING " Alpha = ###.#"; alpha(i) * rtd LOCATE 15, 1: COLOR 6: PRINT USING " Phi = ###.#"; phi(i) * rtd LOCATE 16, 1: COLOR 6: PRINT USING " Beta = ###.#"; beta(i) * rtd LOCATE 17, 1: PRINT USING " Cl = #.##"; Cl(i) ' Plot vectors y1 = -.1 * Vtorq: y2 = 1.05 * Vtorq x1 = 1.4 * y1 - Vtorq / 2: x2 = 1.4 * y2 - Vtorq / 2 WINDOW (x1, y1)-(x2, y2): CLS ' torque vector COLOR 4: LINE (0, 0)-(0, Vtorq(i)) angle = pi / 2: tipx = 0: tipy = 0: GOSUB arrow ' wind vector COLOR 5: LINE (0, 0)-(V, 0) angle = 0: tipx = V: tipy = 0: GOSUB arrow 'INPUT "wind vector", vbv ' axial interference vector COLOR 6: LINE (V * aa(i), Vtorq(i))-(0, Vtorq(i)) angle = pi' 0' pi - ATN((aad(i) * vtorq(i)) / (V * aa(i))) tipx = 0: tipy = Vtorq(i): GOSUB arrow ' torque interference vector COLOR 7: LINE (V * aa(i), Vtorq(i) * (1 + aad(i)))-(V * aa(i), Vtorq(i)) angle = pi / 2 tipx = V * aa(i): tipy = Vtorq(i): GOSUB arrow ' downwash COLOR 8 LINE (V * aa(i), Vtorq(i) * (1 + aad(i)))-(0, Vtorq(i)) angle = pi - ATN((aad(i) * Vtorq(i) / (aa(i) * V))) tipx = 0: tipy = Vtorq(i): GOSUB arrow ' relative wind (resultant) COLOR 9: LINE (V * aa(i), Vtorq(i) * (1 + aad(i)))-(V, 0) angle = pi + ATN(Vtorq(i) * (1 + aad(i)) / (V * (1 - aa(i)))) tipx = (V * aa(i)): tipy = Vtorq(i) * (1 + aad(i)): GOSUB arrow ' plot airfoil (not true chord) y1 = -.1 * Vtorq: y2 = 1.05 * Vtorq x1 = 1.4 * y1 - Vtorq / 2: x2 = 1.4 * y2 - Vtorq / 2 x1 = 0: x2 = 2 * chord75 y1 = x1: y2 = .7 * x2 WINDOW (x1 - V, y1)-(x2 - V, y2) ' WINDOW (x1, y1)-(x2, y2) ' WINDOW (0 + V, 0)-(2 * chord75 + V, 1.4 * chord75) COLOR 12 angle = -(beta(i) + pi / 2): GOSUB naca scl = 1 FOR pp = 2 TO 14 LINE (scl * xu(pp - 1), scl * yu(pp - 1))-(scl * xu(pp), scl * yu(pp)) LINE (scl * xl(pp - 1), scl * yl(pp - 1))-(scl * xl(pp), scl * yl(pp)) NEXT pp COLOR 10: LINE (0, 0)-(scl * xu(1), scl * yu(1)) ' chordline COLOR 2: LINE (0, 0)-(scl * x0, scl * y0) ' zero-lift line RETURN turbinedata: FOR i = 1 TO 25 ' parameters for 25 blade elements T(i) = .12 ' thickness-to-chord ratio of airfoil M(i) = .02 ' camber P(i) = .35 ' camber high point chord(i) = -chord75 * (i - 1) / 17 + 2 * chord75 ' initial values for blade chord alpha0(i) = -ATN(M(i) / (1 - P(i))) ' airfoil zero lift angle, radians Cd(i) = .03 ' drag coeff: put your guess here! Vtorq(i) = 2 * pi * (i / 25) * R * RPS ' blade element speeds,m/s in torque plane NEXT i RETURN SUB acos (x, cosine(), angle(), angle) ' returns the angle whose cosine is equal to x pi = 3.14159 i = 1 switch = 1 DO i = i + 1 IF i < 181 THEN IF x > ABS(cosine(i)) THEN switch = 0 ELSE switch = 0 END IF LOOP UNTIL switch = 0 ' loop until i>181 or x>abs(cosine(i)) IF i < 181 THEN angle = (x - cosine(i - 1)) / (cosine(i) - cosine(i - 1)) angle = angle * (angle(i) - angle(i - 1)) + angle(i - 1) angle = angle * 180 / pi ELSE PRINT "subscript of of range" END IF END SUB SUB plotext (c, x, y, x$, position$, script, x1, y1, x2, y2) ' This subprogram prints the text specified by the string, X$ centered on the ' co-ordinates X,Y in the colour specified by the variable COLOUR. ' The text can therefore be printed anywhere on the graghics screen rather ' than being limited to just the normal text rows and columns. DIM ARRAY(1 TO 34) Nspace = 8: Ysub = 0 COLOR c LGT = LEN(x$) ' LGT = length of X$ SELECT CASE position$ CASE "C" xc = x - (LGT * 8 * (x2 - x1) / 640) / 2' THIS CENTRES THE TEXT ABOUT X CASE "S" xc = x ' THIS STARTS THE TEXT AT X CASE "E" xc = x - LGT * 8 * (x2 - x1) / 640' THIS ENDS THE TEXT AT X CASE ELSE xc = x - (LGT * 8 * (x2 - x1) / 640) / 2' THIS CENTRES THE TEXT ABOUT X END SELECT SELECT CASE script CASE 0 Ysub = 0 ' PRINT NORMALLY CASE -1 Ysub = Ysub ' PRINT AS SUBSCRIPT CASE 1 Ysub = -Ysub ' PRINT AS SUPERSCRIPT END SELECT SCREEN 9, 2, 0, 1 ' make the invisible screen (0) active CLS 0 ' clear the invisible screen COLOR c ' select colour WINDOW (0, 349)-(639, 0) LOCATE 23, 3: PRINT x$ ' print FOR i = 0 TO LGT - 1 WINDOW (0, 349)-(639, 0) GET (15 + 8 * i, 31)-(22 + 8 * i, 39), ARRAY ' GET text letter by letter SCREEN 9, , 1, 1 ' MAKE VISUAL SCREEN ACTIVE WINDOW (x1, y1)-(x2, y2) Xplace = xc + 8 * (x2 - x1) / 640 * i Yplace = y IF Xplace < x2 AND Xplace > x1 THEN PUT (Xplace, Yplace), ARRAY'PUT letters of text to visual screen END IF SCREEN 9, 2, 0, 1 ' MAKE INVISIBLE SCREEN ACTIVE NEXT i SCREEN 9, , 1, 1 ' MAKE VISUAL SCREEN ACTIVE x = x + LEN(x$) * Nspace script = 0 END SUB 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