THE UNIVERSIT Y OF MICHIGAN COMPUTING CENTER Tec-hnic al Report SIMU~LATION OF CONTROL OF MODEL POSITION AND WINDFLOW IN A TRANSONIC WIND TUNNEL John S. Tripp ORA Project 01646 suppore ted'by NATIONAL AERONAUTICS AND SPACE ADMINISTRATION GRANT NO, NGR-25-005-28`; WASHINGTONP Do Co administered thrrough: OFFICE OF RESEARCH ADMIN'TRATION ANN ARBOR February 1969

TABLE OF CONTENTS Page LIST OF FIGURES iv INTRODUCTION 1 MODEL ATTITUDE CON TROL 5 DYNAMIC WINDFLOW CONTROL 12 COMBENED WINDFLOW AND MODEL ATTITUDE CONTROL 17 APPENDIX —FLOW CHART AND LISTING OF SIULATION PROGRAM OF COMBINED ANGLE OF ATTACK AND MACH NUMBER CONTROL2:1 iii

LIST OF FIGURES Figure Page 1. Sting deflection versus true angle of attack. 2. Continuous halving algorithm. 3. Equal increment algorithm. 7 4. Quantized halving algorithm. 8 5. Single checkpoint algorithm. 6. Direct algorithm. 10 7. Data Set 1-Mach number versus spoiler setting. 14 8. Data Set 2'-variation in Mach number with angle of attack. 15 iv

INTRODUCTION Positioning of models and control of windflow in wind tunnels have normally been accomplished by means of human operators who must manipulate manual controls Several major advantages to be gained by automating these procedures are as follows: 1. Accurate model positioning, which is currently unfeasible to obtain because of support deflection due to wind forces. 2. Realization of a decrease in tunnel run time because of minimization of time required to assume test positions. 3. Automatic maintenance of windflow conditions under varying airflow blockage caused by changes in model attitude. 4. Automatic programming of aerodyna.mic tests under computer control. This report describes a study performed for the Eight-Foot Transonic Pressure Tunnel at the Langley Research Center of the National Aeronautics and Space Administration, in which control algorithms were developed for automatic control of model attitude and windflowo Included are descriptions of digital computer simulations of the performance of the various algorithms which were proposed and the ultimate combination of algorithms chosen to be implemented in the facility. In the 8-foot tunnel an aircraft or space vehicle model is attach ed tc a test balance, which contains strain gages for measuring six basic force and m.ment components; the balance is attached to the sting, a long tapered steel support located in the test section of the tunnel. The sting, in turn, is attached at the large end to a movable, massive support strut. During tests model angle of attack is varied by changing strut position by means of a drive motor, which is coupled through a magnetic clutch to a gear drive mechanism. Drive rate, which is regulated by the magnetic clutch, is determined by the setting of a control potentiometer. If the mechanism is driving and the drive control switch is turned off, a friction type drum brake engages, and the drive motor is uncoupled from the gear drive, but the strut assembly continues to coast for a variable distance, depending on drive rate and direction. The six basic force and moment components, which are measured by means of the test balance, are normal force, side force, axial force, pitching moment, rolling moment, and yawing moment. Model attitude is specified in terms of angle of attack (pitch angle), angle of yaw, and angle of roll. This study considered the case where only angle of attack was varied and where angle of roll and angle of yaw remained equal to 0~. 1

Wind forces acting upon the model cause sting-balance deflection, which results in an error between true angle of attack and indicated angle of attack. This error can be computed as a function of normal force and pitching moment. However, the nature of the test balance causes nonlinear interactions to appear between the indicated values of the basic forces and moments. These latter values must therefore be corrected for interactions, which is done by means of solving a set of six simultaneous nonlinear equations, before an accurate value of sting-balance deflection can be calculated. 2

MODEL ATTITUDE CONTROL The first phase of the study investigated the control and systems operations required to position the model precisely to a specified angle of attack with wind on. The required function of a control algorithm for angle of attack positioning is to provide a time sequence of control commands to the drive mechanism so as to cause the sting to drive to a specified true angle of attack, within a specified error limit, without overshooting, within minimal time. A control command specifies drive rate and drive direction. In order to accomplish its required function the controller must continually monitor indicated angle of attack and periodically calculate sting bending and true angle of attack. As the target is approached the drive rate may need to be reduced; eventually the brake is engaged at a point such that the sting will coast to the target position. This braking operation requires a knowledge of strut coast for all anticipated drive rates and directions. The skeleton control algorithm is as follows: 1. Read new target value for angle of attack, OT 2. Measure current true angle of attack, a. a. Read raw forces from data acquisition system. b. Correct raw forces for balance interactions. c. Compute sting deflection anid a. 3. Compute error between a and a. If target has been reached, stop; return to step 1. 4. Select next checkpoint, i.e., a at which sting deflection is to be calculated, where O = uncorrected a, 5. Predict sting deflection and a at checkpoint. If system is to stop at checkpoint, estimate coast after brake is applied. 6. Determine control command and drive to checkpoint. 7. Return to step 2. The operation of the model positioning mechanism under control of this algorithm was simulated on a digital computer, Some aspects of the simulation are now dis(cussed. 3

The six basic forces and moments were assumed to be functions of indicated angle of attack (aR) only, for constant Mach number and windflow parameters. Step 2a of the control algorithm was then simulated by obtaining the raw forces from static force data taken from aerodynamic tests of an aircraft model, which were run in the tunnel. This empirical data consisted of a set of the six forces and moments plotted versus indicated angle of attack, for constant Mach number. For each value of aR, the forces were determined by interpolation. Dynamic behavior of the strut drive and braking apparatus was simulated assuming constant acceleration to fixed drive velocities, constant fixed drive velocity dependent only on control potentiometer setting, and constant deceleration to zero velocity during braking. Experimental data was taken at the wind tunnel which showed these assumptions to be reasonable and from which values of acceleration, deceleration, and drive velocities were obtained. Using these values a dynamic routine was written for the simulation program which, when given the current time, sting position, drive velocity, and control command, calculates the time, sting position and drive velocity at the next checkpoint. The different versions of the control algorithm, which were tested using the simulation, differed at the following steps: 1. Sting deflection prediction. 2. Checkpoint selection. 3. Strut coast estimation. Description of each of these variations and its performance will now be given. The first version of the control algorithm utilized zero order sting deflection prediction which predicts the deflection at the current checkpoint to be equal to that at the previous checkpoint. Its performance is satisfactory only if the increments between checkpoints are kept small; this results in the disadvantage of having to perform frequent sting deflection calculations, each of which may require as much as 0.2 second of computer time. All later versions of the control algorithm utilized first order prediction (linear extrapolation)o With the empirical force data used in the simulation first-order predictions were made to an accuracy of ~0.020, or less, over checkpoint increments of 1~. Figure 1 shows a curve of sting deflection versus true angle of attack for a supersonic transport model which was tested in the 8-foot facility. The curve is nonlinear, but it can be seen that linear extrapolations over increments of 1~ can be made accurately. However, linear extrapolations over large intervals such as 10~ will be subject to considerable error. Hence checkpoint increments must be kept to the order of 1~. 4

3.2 30p SST - Mach 0.6 2.8 2.6 2.4/ 2.2/, 2.0 4) 1.8,/ o 1.6/ C 1.4 C. 1.2-.2.82 0.4/ 0 0 I 2 1 1 1 2,4 5 6 7 8 9 10 11 12 13 14 15 -0.2 -0.4 / True Angle of Attack in Degrees -0.6 Figure 1. Sting deflection versus true angle of attack. 5

The first checkpoint selection scheme to be tested selects each checkpoint increment to be equal to half of the remaining distance to the target until the latter value becomes less than some limit, E, whereupon the increment is made equal to the entire remaining distance to the target. In order to limit the number of sting deflection calculations, c is kept relatively large; the consequence is wide variability in the final drive increment, which results in widely varying final prediction errors, especially with zero order prediction, and a tendency for overshoot. The results of simulation tests of each variation of the control algorithm are presented graphically as a point plot of target error versus the number of sting bending computations per degree of change. Since for the optimal system one desires to minimize both of these parameters, a desirable performance characterization will consist of a compact cluster of points near the origin. The plot for the halving algorithm described above with first-order prediction appears in Figure 2. Unfortunately, too few points were available to obtain a meaningful plot..11t.10 noQ-.08 4, a.07 c.06 0 u.05 C.C.04 o, 03 0 F 1-, 0 i 11.LU 1.01 - 0 I I I U I 2 3 4 Number of Deflection Calculations per Degree of Displacement Figure 2. Continuous halving algorithm. In order to eliminate variability in the final drive increment, several versions of an increment selection algorithm were devised in which the value of the final increment is prechosen, and all previous increments, except the first, 6

are made equal to some multiple of this value. In the simplest version of this scheme, all increments are made equal to the final one, except the first, to which the remainder is added. Using 0.5~ as the increment value, and first order prediction, the positioning accuracy was excellent (less than 0.005~ error) for any amount of displacement between target values of a; however, over wide angle changes numerous computations of sting deflection were necessary. See Figure 3 for the performance plot..11.10.09.08 -.07 4)..05.04 LJ.03.02 ~ I Operoting Rnge o.01' pl 0n F Qs 0 I 2 3 4 Number of Deflection Calculations per Degree of Displacement Figure 3. Equal increment algorithm. The best features of the two previously described increment selection schemes were combined into a quantized halving increment selection algorithm in which increments are chosen backward as a sequence of multiples of the minimum increment, where the sequence of multiplying factors is an increasing sequence of powers of two. The first two increments at the front, which normally will not fit the sequence, are chosen to be half the distance to the nearest point of the sequence. The performance of this algorithm is superior to that of the others since it provides the accuracy and consistency of the equal increment scheme above, while requiring fewer calculations of sting deflection. Where the former requires n calculations this method requires log2 n + 1 calculations. If raw data can be acquired with no significant delay and the calculations can be performed in real time, the performance of the overall control algorithm is almost time optimal, with positioning errors falling within the allowable limits (see Figure 4). 7

.11 I. 10.09,.08 0).07.06 C I._ S.05 w.04 0 c.03..020 o *.01 - o 0 0 1 2 3 4 Number of Deflection Calculations per Degree of Displacement Figure 4. Quantized halving algorithm. Unfortunately the data acquisition system at the 8-foot tunnel has lowpass filters installed in the force data channels which have time delays of approximately 5 seconds. Their purpose is to average the variation in forces caused by model buffeting at high angles of attack. This delay makes realtime calculations of true angle of attack impossible. Hence an increment selection algorithm with a delay was devised which functions as follows: While the drive mechanism is stopped the current true angle of attack, a2, is computed. The difference I(T - o21, if less than 1.5~, is assigned to the increment, the checkpoint becomes the target, 1, and the system then drives directly to co. Otherwise the checkpoint is made 1 short of CT. The system drives to the checkpoint, stops, holds sufficiently long to allow for filter rise time, computes true a, predicts sting deflection at the target, and then drives on to aT. The positioning accuracy is inferior to that of the quantized halving increment selection algorithm, and the response is not close to being time optimal. However, in most cases tested in the simulation the prediction error was less than 0.02~. The error plot for this scheme appears in Figure 5. Al algorithm was tested briefly which utilized no intermediate checkpoints but drove directly from one target to the next. For small displacements between targets the performance was adequate, but over excursions of 100 prediction errors exceeded 0.35 which was not tolerable (see Figure 6). 8

.Ill 0.10.09 n.08 C.07 4) 0 c.06 2.05.O04 c C o.03 o 4.0 0 o.02 0 0.010 3 0 1 2 3 4 Number of Deflection Calculations per Degree of Displacement Figure 5. Single checkpoint algorithm. Comparing Figures 3 through 6 it is seen that the quantized halving algorithm of Figure 4 exhibits the best performance, since the plot most nearly approaches the origin. This algorithm is also insensitive to position displacements as the errors remain nearly constant. One would select this scheme over the others were it not for the delay considerations described above. The final characteristic of the control algorithms under discussion which was varied was the technique for estimating strut coast when braking to zero velocity. It can be calculated by the equation r2 d - + rt 2A D where d = distance coasted r = drive rate when brake command is given A = deceleration t = time delay before brake engages D Use of this equation provides an accurate estimate of coast. It provided an 9

1.00.90 o.80.70.60C.c 0 o...c.50 t C ~.40 0 a0.30- o.20.10 8 0 0 0g 1 2 3 4 Number of Deflection Calculations per Degree of Displacement Figure 6. Direct algorithm. 10

exact estimate of the simulated coast since the latter was computed using the same constants. The value of d was determined experimentally to be approximately 0.05~ +0.01~ for drive speeds corresponding to a control potentiometer setting of 2. It was not position dependent. A simple coast estimate technique, therefore consists of always estimating the coast to be 0.05~. Its use in the simulation resulted in satisfactory performance with total errors generally less than 0.02~. There is some tendency for coast to vary slowly with changes in ambient temperature. An automatic corrective scheme was added to the simple coast estimation technique described above. Separate coast estimates are maintained for each direction of travel. If the coast estimate error exceeds some limit, a correction is made to the estimate. The performance of the self-adjusting coast estimation method was tested using the simulation. The values of the corrections finally converged to 0.055~ and 0.046e for upward and downward directions, respectively.

DYNAMIC WINDFLOW CONTROL The requirements for windflow control are twofold: 1. To maintain a constant Mach number under changes in model attitude. 2. To automatically adjust Mach number to a new target value under program control, Requirement 1 arises because of the increasing throttling effect caused by the model upon the windflow with increasing angles of attack, which results in a decrease in Mach number. Compensation for this effect can be achieved with relative ease for all but larger models by means of the diffuser spoiler, which is a part of the windflow control mechanism. The simulation to be described here models a computer centered system for controlling diffuser spoiler settings to effect requirement 1. The diffuser spoiler consists of a pair of pivoted flaps installed downstream from the test section of the tunnel, which may be retracted flush with the tunnel walls or lowered into the windstream, to a maximum angle of 20~, to cause a choking effect. At the beginning of each test, with the model at 0~ angle of attack, the spoiler is set to some value, say 10~. Then, as the model assumes greater angles of attack, the spoiler is withdrawn sufficiently to restore Mach number to its target value. A simple control algorithm which compensates for Mach number variation due to angle of attack changes by means of the diffuser spoiler is now given. The basic element of the algorithm is a second order predictive technique which predicts the spoiler angle setting necessary to restore Mach number to the target value, based on previous Mach number values and corresponding spoiler angle settings. The sequence of operations is as follows: 1. Save old value of Mach number, M, and spoiler angle setting, 8. 2. Measure M at new angle of attack. 3. Is M sufficiently close to target value' M? If so, stop. T 4. Predict value of spoiler angle setting, %D, required to restore MtoMT. 5 Evaluate spoiler subtarget angle p = ( - )R + t where t = current spoiler angle setting, R = proportionality parameter. 6. Does b exceed allowable range for t? If so, stop and indicate error condition. 12

7. Drive to. 8. Determine M at s. Return to step 3. P The effectiveness of the algorithm was investigated by means of a computer simulation which used empirical data obtained from tests at the 8-foot tunnel. The independent variables affecting windflow are: 1, Fan velocity, N, rpm. 2. Free stream total pressure, H, lb/ft2. 35 Diffuser spoiler angle setting, 5, deg. 4o Reentry flap setting, y, counts. 5. Configuration of the model. 6. Angle of attack, a, deg. Two families of curves were obtained experimentally. The first family, to be denoted Data Set 1, consists of four curves of Mach number versus spoiler setting, all other independent variables held constant, with no model, for four values of fan velocity; the set is shown in Figure 7. The second family of curves, denoted Data Set 2, was obtained from aerodynamic test data of a research model, and consists of three curves of spoiler setting required to maintain constant Mach number versus angle of attack, for three values of Mach number. This presentation of the independent variable, 5, as a function of a and the dependent variable, M, was undesirable. Therefore, using Data Set 1, Data Set 2 was transformed into a set of curves of M versus a for constant 5; this new set, shown in Figure 8, will be denoted Data Set 2?. Data Sets 1 and 28 were used as the basis for simulation. The state of the system is represented by a five-element state vector, S = (a,5,7,)N,H) Now M is a function of S; however, for the case being simulated y, N, and H are held constant. Hence, using Data Set 2v, given the initial values ao, o0, and Mo, for a new value of a, z,, one can determine the corresponding value of Mach number, Mi. Likewise, using Data Set 1, given a new value of, 52, the corresponding value of Mach number M2 can be determined, assuming that the curves of Data Set 1 are not affected by the insertion of a model. Hence steps 2 and 8 13

1.2 1.1 I.0 0.9 0.8 0.7 1). z o (c C - - - - - Note: Values of N and ywere not Obtained N3, Y3 - - - 0.6 N4 \N 0.5 0.4 0.3 0.2 0.1 I I I I I I I I I I I 0 2 4 6 8 Spoiler Setting 10 12 14 in Degrees 16 18 20 Figure 7. Data Set 1-Mach number versus constant N, constant 7, and empty tunnel. spoiler setting for 14

*90r N=698 \ 8=7.0 Variation in Mach Angle of Attack Number with.85h Note: y=7261.80-.:3 E z (_u 0 2.75 70 on.65.60 N =538 ----— 8=96 I I I I I I I I I -4 -2 0 2 4 Angle of 6 8 10 12 14 16 Attack in Degrees Figure 8. Data Set 2'-variation in Mach number with angle of attack.

of the control algorithm are simulated by means of interpolating Data Sets 2' and 1, respectively. Step 4 of the control algorithm predicts the value of 6 which will restore M to MT. The prediction is effected by means of a second order extrapolation which uses the current and two previous values of M and the corresponding values of 8, Whenever angle of attack changes, causing an increment AM in Mach number, the two stored prior values of Mach number are incremented by AM, which assumes, in effect, that a change in a results in a translation of the & versus a curve left or right. The performance of the predictive scheme in the simulator indicated that such a translation of the 3 versus a curves is an effective technique. The simulation of the Mach control algorithm alone, then, consists of performing the eight steps of the control algorithm, iterating steps 3 through 8 until the target value is reached, for several angles of attack. It was evaluated for values of the proportionality parameter R, of 0.5, 0.8, 1.0, and 1.2, where R is the proportion of the distance to the predicted spoiler setting to be taken during each iteration. The performance of the control algorithm was dependent on the value of R, but in every case the resultant Mach number converged to the target value after a few iterations. For R = 0.5, six or more steps were required to reach the target, as M approached the final value similarly to the response of an overdamped servomechanism to a step input. For R = 1.0, the target was reached in two or three steps. For R = 1.2, the system overshot the target, but converged back to it in four or five steps in a damped oscillatory fashion. Presumably for sufficiently large values of R the system would be unstable. 16

COMBINED WINDFLOW AND MODEL ATTITUDE CONTROL Having devised separate control algorithms for model positioning to specified angles of attack at constant Mach number, and for maintaining constant Mach number under varying angles of attack, the two algorithms were combined. The algorithm for angle of attack control which introduces delays to allow for data acquisition filter delays, was selected because of its intended application in the 8-foot wind tunnel facility. System control proceeds as follows using this algorithm. Assume that the sting drive is stopped; the true angle of attack is computed, using force data from the.model and is compared to the target angle of attack, aT. If the difference is less than 1.50, the system drives directly to the target. Otherwise it drives to within 1~ of the target, and halts for a period long enough to obtain filtered measurements from which a can be computed. During this delay the current value of M is computed, and if unequal to the target value, MT, the Mach control routine is called to correct it. At the end of the delay true angle of attack is computed, and the system drives to aT where M is restored to the target value as above. The procedure repeats for each test point of the run, at MT. The combined algorithm functions as follows: 1. Read target values for Mach number and angle of attack, MT and T. 2. Measure current position. a. Read raw forces from data acquisition system. bo Correct raw forces for balance interactions. c. Compute sting deflection and a. 53 Measure Mach number. 4. Compute 1 T - a. If target reached, stop. 5. Select next checkpoint. a. If aT - al < 1.5 aT is checkpoint. b. Otherwise aT - 1 (T + 1 if driving in reverse direction) is checkpoint.T 6. Predict sting bending and OL at checkpoint. 17

7. Determine control strategy, drive to checkpoint and stop. 8. Measure current values of a and Mach number as in steps 2 and 3. 9. Invoke Mach control algorithm to restore Mach number to target value M. T 10. If checkpoint was target return to step 1 for next test point. Otherwise aT becomes new checkpoint; return to step 6 to continue. The performance of the combined algorithm was evaluated by means of a simulation program constructed from the previously described simulation programs of the two component control algorithms. See the Appendix for a listing and flow chart of the simulation program. A step-by-step outline of the simulation is now given. 1. Initialize program. 2. Read initial conditions. 35 Read target conditions. 4. Initialize test point. 5. Find current position. a. Obtain force vector as function of aR and M. b. Correct raw forces for balance interactions. c. Compute sting bending and a. 6. Compute l aT - al. If target reached, stop. 7. Find next checkpoint. a. If la - al < 1.5, aT is checkpoint. b. Otherwise aT - 1 (aT + 1 if driving in reverse direction) is checkpoint.T 8. Predict sting bending and raw a at checkpoint. 9. Determine control strategy. 10o Compute time and position at which sting drive stops after control stragety is applied to drive to checkpoint, 18

11. Find current value of a using step 5. 12. Compute Mach number as function of a. 13. Repeat steps 11 and 12 to eliminate interaction between a and M. 14. Call Mach control routine to restore M to target value M. Obtain new value of spoiler setting. 15. If checkpoint was target, return to step 3 for next test point. Otherwise aT becomes new checkpoint. Return to step 10 to continue. A discussion of the performance of the combined simulation is now given. Empirical data consisted of a family of static curves of the six basic forces as functions of angle of attack for three values of Mach number (Data Set 3), in conjunction with the family of three static curves of Mach number versus angle of attack (Data Set 22), which was used in the windflow control simulation previously described. The two data sets were taken, unfortunately, using two different models, A single set of correlated windflow and force data was unavailable from the tunnel. Interaction between Mach number and angle of attack was determined, as noted in step 7 above, by iterating the calculation of a as a function of M and the calculation of M as a function of a. Since the target values for a in the test data were separated by 2~, the checkpoint positions differed by 1~. Between the range of 0~ to 10~, the largest resultant change in Mach number, for 1~ increment in a, was 0.0080. The corresponding change in the calculated value of c, after one iteration was then 0.02140. Recalculating M yielded a change of less than 0.0002, indicating that one iteration was adequate to determine interaction. between Mach number and angle of attack for the test data. The performance of the angle of attack positioning algorithm and of the windflow control algorithm was not degraded by their combination into a single one, except, when more than one iteration is required, the time delay necessary for Mach stabilization (approximately 4 seconds per iteration) increases the amount of time required to reach a target state beyond that required for angle of attack control alone. (Note that the transient behavior of Mach number was not simulated; hence time delays were only estimated by constant values.) The accuracy of angle positioning was the same as with the parent algorithm. Likewise, the number of trials required by the windflow control routine to restore M to MT was unchanged; using a value of 1,0 for R, no more than two steps were required to restore Mach to the target value. See the appropriate sections of the paper for details about the parent algorithms. 19

APPENDIX FLOW CHART AND LISTING OF SIMULATION PROGRAM OF COMBINED ANGLE OF ATTACK AND MACH NUMBER CONTROL 21

ro ^'

REAL XSD(6) 0.,1.,2., 4.,8,,0./,YSO-(6y/0,.02,.25,.39.56,63/, 1SEL (6),K.YN( 12) YIN IT (6),YINROT ( 6 ), STNG/-100./, _..2'YNRDG'( 6), YNZERO(6), KNALF, KMALF,DC B RK (2)-/- 7,-, 9, 3ACCEL(2,2)/.83.81 - 54,- *6/ TDA/. 2/, TB/2/ COMMON AK(162),YNF(6),TARL(6),FITL(6),XNF(6),E(6),EPS, IRK I,RK2,TM,DELTADM,N,ALF2T,L ___L__ NAMELIST/NL A/ AK,KNALF,K MALF,KYN,S E LYINITYINROT,YNZiERO, -- lALFZRO,EPS9,RKiK K2/NLC/ALFCNTT,ALFTRM,TMDFLTA,MN __ RE AD(5,NLA) CALL TFRPIZ Dn 12 1=1,6 o iFY=YINIT( I)-YINROT( I ) IT=I..... IF(DIFY.LT.O,.)' IT=6 _+6. YNF(IT)=DIFY*SFL( I )*KYN(IT)/2. TARL ( )=O. —' —0................... 12 EITL(.I)=O. CALL INTR(O) 00 14 1=1,6 TARL( I )=XNF( I ) 14 EITL(I)=E(I) ALF2T=O. EPALFT=O. 200 READ(5,NLC) 202 WRITE(6,NLC) ITERCT=O CPS=O. ALFIND= (ALFCNT-ALFZRO) /STNG BDR IVE=O,. NSTOP=1 NSA=1 NCL= 240 IF(ITERCT.NE.0) EPALFT=EPALF IF(ITERCT.NE.O) ALF2T=ALF2 245 CALL TERP(ALFCN'T,TM,YNRDG, &200 ) DO 11 1=1,6 D._IFY=YNRDG( I )-YNZERO(I) ____....... IT=I IF(DIFY.LT.O.) IT=I+6 11 YNF(I)=DIFY*SEL(I)*KYNMIT) CALL INTR(O) TEPALF=KNALF*XNF( 1)+KMALF*XNF(3) T ALF2=ALF INNOTEPALF IF(L.EQ.1) EPALF=TEPALF GO TO (.310,43.0,430), NSTOP 430 CALL SPMACH(TALF2) WRITE(6,435) TALF2,TM..435 FORMAT('OPRIOR TO SPOILER ADJUST'/SX,t'ALPHA = ___ 1F 9.4,3X, MACH.=',F9,6) IF(L.NE.1) GO TO 440 25

L=L+1......" GO TO 245.......... 440 CALL SPSET(ALF2) ^...... -. -.-.-.. ~..~.....~.-..-..._...-.-.~.. —~.-..~.......... __ ^.................. L=l IF(NSTOP.EQ.3) GO TO 310 T=T+5..... BORV3=O. NCEL=l1 BDRIVE=O. GO TO 410........._ 310 WRITE(6,?65) ALFIND,ALF2,TM, T 265 FORMA T ('ORAW ALPHA =', F9.4, 3X,' TRUE ALPHA't=F.4-T3X.MACH..T 1F9.6,3X,'T=',F6.2) TS=ALFTRM-ALF2 IFINSTOP.EO.3) GO TO 280 IF(T-ERCT.GE.25) GO TO 500 I TERC T=I TERCT+l IF(ITFRCT.EQ.1.AND.TS.GE.O.) NSA=2 CPS=6. ATS=ABS (TS )............ NTS=ATS* 100. +5. I. "i(NTS.LT.150) GO'TO 410 DELALF=TS-SIGN( 1.,TS) NSTOP=2 GO TO 412 323 PALF2=ALF2+DELALF ALFCNT'=ALFCNT+DTHETA*S'TNG ALFIND=ALFIND+DTHFTA WRITE(6,325) CPS,ALFINO,PALF2 325 FORMAT('OORIVE CONTROL STRATEGY'/2X,'CPS =',F6.2,3X 1'DATA PT ALPHA =,F9.4,3X,,PRFD TRUF ALPHA =',F9.4) BORV2=BDR IVE NAC=1 CALL INTERP(XSD,YSD, 1,6,CPS, BORIVE'E,,6) IF (BDRV2-BDRIVE).GE.0.) NAC=2 MB=l ACL=ACCFL(NSA,NAC) 370 TD=TDB IF(BDPV2.EO.Q.) TD=TOA IF(NCEL.EQ.1) GO TO 375 DT3=( BDRV3-BDRV2)/ACL IF(DT3.LT.O.) GO TO 420 IF(DT3.GT.TD) GO TO 377 DX3=DT3*( BDR. V3+BRV2)/2.+BDRV3*( TD-D'T3).B DR V? =B DR V3 372 oT=(BDRIVE-BDV2 )/ACL DX=DT*(BDRIVE+RDRV2)/2.+OX3 DT=DT+TD GO TO (337,362),MB 375 DX3=TD*BDRV2 GO TO 372'"377' T.BD2= BDR V2+ACL*TD DX3=TD* (TBD?+BDR.V2 )/2. BDRV2=TB02 GO TO 372 337 IF(ADT.LT.DX) GO TO 340 T= T+DT+ I ADT-DX) /BDRR I V F NCEL=i GO TO 351 340 IF(ADT.LT.DX3) GO TO 342 26

BDRV3=BDPI VE..-.. TSOBD=B DRV2*2+2*, ACCEL(NSA, NAC )*( IAT —DX3- ) _ IF(TSQBD.LT.0.) GO TO 290 fBDR IVE=SQ RT(TSQ"B)D" T=T+TD+2* (ADT-DX3)/(BORIVE+-BDRV2) GO TO 350 342 BDRV3=BOR.IVF BDRIVE=BDRV2' T=T+ADT/BDR I VF __ 350 NCEL=2 351 IF(NSTOP.EQ.O) GO TO 240'360....MB=2 WRITE(6,620) TBDRIVF BDRV2 BDRV3 620 FORMAT('ODRIVE CONTROL STRATEGY'/5X,'BRAKE ON -lOXTh'T.. lF8.4/5X,'BDR.IVE=, FB,4,5X,I'BDRV2=',F8.4,.5X'BODRV3=,tF8.4) BDRV2=BDR IVE BDRIVE=0. ACL=DECBRK(NSA) GO'TO 370 362 T=T+DT.IF(NSA.EQ.1) OX=-DX ALFCNT=ALFCNT+DX *STNG... ~-...... — ALFIND=ALF'IND+DX GO TO 240 410 NSTOP=3 TS=ALFTRM-ALF2 DFLALF=TS....412 OEPALF= ( EPAL'F-EPALFT )*DELALF/( AIF2-'AL'F')2T~).. IF(ALF2.FQ.ALF2T) DEPALF=O. CALL INTERP(XSD,YSD,1,6,CPSTBD, 1,6) ACOAS'T=-TBD*TBD/( 2.*DECBRK( NSA) ) +TDB*TBD IF(NSA.EO.l) ACOAST=-ACOAST DTHETA=DELALF-DEPALF-ACOAST ADT=ABS (OTHETA) IF(OTHETA*DELALF) 360,'360,323 500 WRITF(6,501) 501 FOPMAT('O.TERATION LIMIT MAIN') GO TO 200 280 WRITF(6,281) TTS __ __ 281 FORMAT('OTARGET ALPHA REACHED T =',F7.2,' FRROR =.',F9.5) GO TO 200 290 WRITE(6,291) 291 FORMAT('OFATAL DRIVE ERROR') GO TO 200 420 WRITE(6,421) 421 FORMAT'DTE. LT. ) GO TO 290 END SUBROUJTINF TERP(A,TM,Y, *) RFAL F(6,7, 3),X(7,3),XM(3),Y{6),W(2) PD(R)=P+I PFJN( XI,X2,F',F2 ) =( ( XI-A) *F2-( X2-A) F 1 ) / (X -X?)'.... IF(TM.LT.XM(1.).OR.TM.GT.XM(3)) GO TO 290NR=O Jl=l IF(TM.GT.XM(2)) Jl=2 J2=J1+l DO 120 J=J1,J2 IF(A.LT.X(1,J).OP.,A.GT.X(7,J)) GO TO 2RO 120 CONTINUF 27

DVD= (TM-XM(J 1 ) ) / ( XM( J2)-XM(J 1)) DO 170 K=1,6 DO 150 LK=1,5 _____________ L=LK IF(A.LT.X(L+1,J)) GO TO 160 150 CONTINUE 160 DO 165 J=J1,J2 P1=PFUN(X(L,J),X(L+1,J),F(K,L,J),F(K,tL+I,J)) P2=PFUN(X(L,J),X(L+2 J ),F(K,L,J),F (K,L+2J)) W'J) PFUN(X{L+i,J)' XL+"-2JPI, P 2) 165 CONTINUE Y(K)=W(Jl)+(W(J2)-W(J)1 *DVO 170 CONTINUE 175 FORMAT('OTERP Y =1/6(2X,F8.2) RETURN 280 WRITE(6,281)A,TM 281 FORMAT('OTERP OVER RANGE A='',G15.85X,'TM =',G15'.'8) RETURN 1 ENTRY TERPI7 NAMEL IST/NLD/XM, X,F REA FD( 5,NLD). RETURN END SUBROUTINE INTR(NRI) DIMENSION El(6),AM(42), ICODE(6),EE(6)NXN12) COMMON AK(162 ),YNF (6),TAPL (6 ), EITL 6 ),XNF('6'T'r-"(~ --- NAMELIST/NL 1/XNFF DATA NXN/1,2,3,4,5,6,10000,1OOOC,10000, 10000,010000,10000/ NTY=l DO 14 I=1,6 EE(I)=0. E(I)=0. XNF(I)=YNF(I)+TARL(I) 14 ICODE(I)=1 16 DO 25 M=l,6 IF(NXN(M).EQ.0) ICODE(M)-= IF(ICODE(M).EQ.O) GO TO 25 K=6 DO 18 1=1,6 DO 17 J=I,6 AM( I )=XNFI) K=K+I 17 AM(K)=XNF(I)*XNF(J) 18 CON TINUE JK=NXN(M) J=(JK-1)*27 IF(ICODE(M)) 19,25,19 19 EF(JK)=O. DO 20 1=1,27 J=J+1 20 EE(JK)=EE(JK)+AK(J)*AM( I) IF(NRI) 24,21,24 21 CX=(EE(JK)-F(JK)-EI'TL(JK) )*.OE8 IX=ABS(CX) IF(IX-NXN(JK+6)) 24,22,22 24.. ICO DE(M)=O 22 E(JK)=EE(,JK)-FITL(JK) IFINRI) 25,26,25 28

26 XNF (JK)=YNF (JK)-E (JK) +TARL(JK 25 CONTINUE IF(NTY) 48,61,48 48' DO 49 1=1,6 49 EI(I)=E(I) NTY =O 61 CONTINUE DO 30I=1,6 28 IF(ICODE I ) 29,30,29 29 GO TO 16 30 CONTINUE DO 31 1=1,6 31'XNF( I)=YNF(I)-E(I) 63 CONTINUF r'101 FOR.MAT 0 INTR XNF ='/6( 2X, F9.3 ) RFTURN END 29

SUBROUTINE SPMACHALFP) RE AL YM-AC F-(3, 0 1),XA'L F(- 0,1,YM-ELct(4t'8-)-,XDE-10 ) X( 3 F(.3 FUN( X )=X X + 1..-...PFUN(Xl, X2,FI,F2 ) = (( XL-DM )*F2-' X2-DM ) F'l) /'X1-X?) COMMON WUI192),EPS,RKI,RK2,TM,DELTA,DM,N,ALF2T,L __ __ - IF(L.EQ.1) OTM=TM TM=FI NT( XALF,YMALF,ALF2T OTM, ALFP, KA, LA) 100 RETURN ENTRY SPSET(ALFP) NCT=1 DTM= TM-OTM X(2)=X(2) +DTM X(3)=X(3)+DTM 10 (-A-'5TDM-TMt r.LT.EPS) G GO TO 200 X(1)=TM F(1)=DELTA IF(N.NE.1) GO TO 120 PD=F(I)+(X(l)-DM)*RK1 GO TO 140'Tiy~ — =.....P'- OPF-UNIX( X- 2), X ( I ) F ( 2 )I F ) )I.. IF(N.EQ.2) Gn TO 140 P21=PFUN(X( 3),X( 1,F(3) F( 1 ) PD.=PFUN (X (3),X (2),P21,PO) 140 DELP'=(PD-DELTA)*RK2+DELTA WRITE(6,142) PD -14 —2 FORM.AT'('OPRFODCTED DELTA =',G1l2.6) IF(DELP.LT.O.OR.DELP.GT.20.) GO TO 220 PMINF=FINT(XDELYMDELDELTA,TM,DFLP,KD,LD) IF(N.LT.3) N=N+1 DO 160 I=1,2 X(4-I)=X(3-T) 160 F(4-I).=F(3-1 DELTA=DELP TM=PMINF WRITE(6,175) DELTA,TM 175 FORMAT(' SPOILER SETTING =',F10.4,3X,' MACH =',FIO.6) IF(NCT.GT.25) GO TO 210 NCT=NCT+ l GO TO 110 200 WRITE(6,201) 201 FORMAT('OTARGET MACH REACHED') GO TO 100 210 WRITE(6,211)'2.1'1[' FORMAT('OITFRATION LIMIT EXCEEDED)' GO TO 100 2.20 WRITE(6,221) DELP 221 FORMAT('ODELP OUT OF RANGE',Fl0.6) GO TO 100 240 DATA XALF/-I. 0.,2.,3.,4.,6.,8., 10.,I2.,14./, IYMALF/,599,.8,1,, 7 94599,,8,,99, 6.799, 998,598,.^Q^ ^ - - 2'.597,.796, 986,.595,.784,.964,.591,.775,.943,,5 R7, 763,.3 1I, 3.581,.749,.917,.576,. 731,.899/,KA,LA/3,10/, 4XDEL/O.,3.,6.,9.,12., 15.,18.,20. /, 5YMDEL/.615.905,,1.2,.605,.88,,985, 1. 195,.598,.835.945, 61.1 75,. 5175,, 79,.87, 1. 12,. 555, 74,. 795,. 975,.53,. 69,. 72 5, 84, 7.5.64,'.665,.755,.485,.6,.635,.715/,/ 8KD,LD,KA,LA/4,8, 3, 10/,F(2),F(3)/0.,0./ END

SUB ROUTINE INTERP( XT,YT,M,NX,YK,L) REAL XT(L),YT(K,L) F R)=R+1 _ PFUN X1,X2,FlF2)= (Xl-X)*F2-(X2-X )*F1)/( X-X? NN=N-2 IF(X.LT.XT(1).OR.X.GT.XT(N)) GO TO 170 DO 120 J=1,NN JK=J IF(X.LT.XT(J+1 ) GO TO 130 120 CONTINUE 130 P12=PFUN(XT(JK), XT JK+1 ),YT ( M JK),YT( M, JK+ 1) ) P13=PFUN(XT(JK), XT(JK+2),YT(M,JK),YT(M,JK+2)) Y=PFUN(XT(JK+1.),XT(JK+2),P12,P13) RETURN 170 WRITE(6,175) X 175. FORMA-T'OX OUT nF RANGE FOR INTFRP, X ='i,"-,g-' —1 " RETURN END FUNCTION FINT (XTYT,XA,YA,XP,K,L) REAL XT(L),YT(K,L),YK(20),YL(2) DO 10 J=1,K CALL INTFRP(XTYT,J,L, XAYK(J),K,L).K —10 CONTINUE IF(YA.LT.YK(1).OR.YA.GT.YK(K ) )GO 130 KK=K-l DO 20 J=1,KK JK=J _ IF(YA.LT.YK(J+) ) GO TO 21 20 CONTINUE 21 DO 22 J=1,2 22 CALL INTERP(XT,YT,JK+J- I,LXPYLJ), K,L) IF(YK(JK).EQ.YK(JK+1 ) GO TO 140 FINT=(YA-YK(JK))* (YL? )-YL( ) / (YK( JK+t)-YK(JK) )+YL(1) RETURN.. 130 WRITE(6,131) YA 131 FORMAT('OYA OUT OF RANGE FOR FINT, YA=',G15.8) F INT=YA RETURN 140 FINT =(YL(t1)+YL12))/2. RETURN. _END

UNIVERSITY OF MICHIGAN 3 9015 03526 7353