Study of the stall-spin phenomena using analysis and interactive 3-D graphics Citation Dubois, Pascale C. (1986) Study of the stall-spin phenomena using analysis and interactive 3-D graphics. Engineer's thesis, California Institute of Technology. doi:10.7907/SFMF-Z114. https://resolver.caltech.edu/CaltechETD:etd-03242008-113141 Abstract The purpose of this study is to gain a better understanding of the nonlinear stall-spin phenomenon through numerical analysis and interactive 3-D graphics. The linear aerodynamic range was thoroughly examined for the NAVION, a light aviation aircraft. Nonlinear aerodynamic behavior was modeled by adding nonlinearities to the lift, pitching and rolling moments. The results of this analysis are promising; however, a more sophisticated model is needed to fully simulate the stall-spin phenomenon. A graphic tool is described which allows the user to interact with the simulation process. This gives the user a "feel" for the dynamics of aircraft and effectively displays the characteristic features of the dynamic model. Item Type: Thesis (Engineer's thesis) Subject Keywords: Stall-Spin ; 3-D Degree Grantor: California Institute of Technology Division: Engineering and Applied Science Major Option: Aeronautics Thesis Availability: Public (worldwide access) Research Advisor(s): Culick, Fred E. C. Group: GALCIT Thesis Committee: Culick, Fred E. C. (chair) Roshko, Anatol Liepmann, Hans Wolfgang Defense Date: 1 December 1985 Record Number: CaltechETD:etd-03242008-113141 Persistent URL: https://resolver.caltech.edu/CaltechETD:etd-03242008-113141 DOI: 10.7907/SFMF-Z114 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 1102 Collection: CaltechTHESIS Deposited By: Imported from ETD-db Deposited On: 04 Apr 2008 Last Modified: 21 Dec 2019 03:57 Thesis Files Preview PDF (Dubois_pc_1986.pdf) - Final Version See Usage Policy. 1MB Repository Staff Only: item control page

STUDY OF THE STALL-SPIN PHENOMENA USING ANALYSIS AND INTERACTIVE 3-D GRAPHICS Thesis by Pascale C. Dubois In Partial Fulfillment of the Requirements for the Degree of Aeronautical Engineer June 1986 (submitted in December 1985) California Institute of Technology Pasadena, California Li ACKNOWLEDGEMENTS I would like to thank Professor Fred E.C. Culick for suggesting the topic treated in this thesis and for his constant encouragement and guidance. My gratitude also goes to the members of my committee, Professor H. Liepmann and Professor A. Roshko, for their advice and concern and to Professor Antonsson for his help with the communications between the host computer and the graphic device. Financial support during the course of this research was provided by the California Institute of Technology. Special thanks go to my friends without whom this work would never have _ been accomplished. Finally, I would like to express my most sincere thanks to my friend, Janet Campagna, for editing this thesis. iii ABSTRACT The purpose of this study is to gain a better understanding of the nonlinear stall-spin phenomenon through numerical analysis and interactive 3-D graphics. The linear aerodynamic range was thoroughly examined for the NAVION, a light aviation aircraft. Nonlinear aerodynamic behavior was modeled by adding nonlinearities to the lift, pitching and rolling moments. The results of this analysis are promising; however, a more sophisticated model is needed to fully simulate the Stall-spin phenomenon. A graphic tool is described which allows the user to interact with the simulation process. This gives the user a "feel" for the dynamics of aircraft and effectively displays the characteristic features of the dynamic model. iv NOMENCLATURE Diweeneeee Wing reference span Coceeseeee Wing mean aerodynamic chord Chie eeeeee Dimensionless drag coefficient CG... eens Center of gravity Ci .+e++++-Dimensionless lift coefficient Cy. ee eees -Dimensionless rolling moment coefficient Cate aeeeee Dimensionless pitching moment coefficient Ci..++...-Dimensionless yawing moment coefficient Cy......--Dimensionless side force coefficient Deve ceeee Drag force P s+++++--Aerodynamic and propulsive force Beeweeee . Acceleration due to gravity Gores eeeee Gravity force nh...... .. Altitude tly I,..Moments of inertia referred to body axis ’ xzereeres Product of inertia referred to body axis Levsccceee Rolling moment about the x-body axis due to aerodynamic torque (positive right wing down) | .. Lift M.........Mass rs -Pitching moment about the y-body axis due to aerodynamic torque (positive nose up) 1 ren Yawing moment about the z-body axis due to aerodynamic torque ( positive nose right) Pewaccceee Roll rate, angular velocity about x-body axis (positive right wing down) Qeececsaee piteh rate, angular velocity about y-body axis (positive nose up) GQeeeeceeee Dynamic pressure, 3 o V2 Povese «...Yaw rate, angular velocity about z-body axis (positive nose right) S.........Reference wing area U..eesees.Forward velocity along x-body axis Up+.s-....0teady~state forward velocity Vesesees ..Velocity along the y-body axis \ ren Magnitude of the velocity vector Vgwseeveee Reference velocity: 2mg=oS58 V2 Weeeeeeee Velocity along the z-body axis Wo sccesece Steady-state downward velocity W.....--.--Weignt of the aircraft X,Y,Z.....Inertial axis (Z pointing downward) X,Y,Z.....Body axis (x pointing forward, z downward) Aesseeee .-Angle of attack Aoeeeees. oLeady-state (trim) angle of attack Ao...e-..-Stall angle of attack Birveees ..-cideslip angle Y..e---.--Flignt path angle atte cseee Aileron control surface deflection (positive for positive rolling moment) 6 geese eee sElevator control surface deflection ( positive for nose-down pitching moment) vil Se eeeeeee Rudder control surface deflection ( positive for nose-left yawing moment) Note positive for negative yawing moment. Bi... eee eee Pitch angle, positive nose up Boeeeeceee Steady-state pitch angle Dosecevees Mass density of air Paveneces .-Roll angle, positive right wing down W...eee...Yaw angle, positive for right nose vii NONDIMENSIONAL DERIVATIVE DEFINITIONS LONGITUDINAL AXIS L eae C = oo L qs’ positive up Cc - 2 positive aft C - _ m qSe C Cc = a-b La oa C Cc _ ab Lé a6 ce Do da Cc . acd DS 36 c . ava Mma = 8a C = acm m6 36 2V Cc “T, g°m Cc mq fe) dq LATERAL AXIS y8 1g ng ys6 16 né lp np lr a w viii ix SIGN CONVENTION TABLE OF CONTENTS Acknowledgments..... AbDStract ccc. c cee c we cveee cece were sce e croc eecebil Nomenclature......... i i iran a a Definitions... creer e creer ccesccccesesvevves Vil Sign Convention... cscccccccccccccsssvcesesseseeiX Table of Contents... .. cee eee cece eer eee eee ene X List of Figures... cc eee e cer cveee see ec eees xii Chapter I : Introduction. ..ccccscenscsesccveel Chapter II: Equations of Motion...... Chapter III: Interactive Display Chapter IV: Linear Aerodynamics.............16 Chapter V : Nonlinear Aerodynamics.......... 2) CONCLUSIONS... ccc eee c ewer ccc enccrece cece cence 29 RePerences.. cc ccc ccc cece eee eee ene c ene eeseees 31 Appendix I : The (u,a,8) Form of the Equations of Motion... . cee ee eee eee 33 Appendix II: Study of the Singularity in Pitch Angle... ..ccccccccccece 200 0 39 Appendix III : Presentation of the Function Network Appendix IV : Linear Aerodynamic coefficients for the NAVION and the CHEROOKEE. 49 Appendix V: Linearized Longitudinal EqQuationS....cccccacecccccccseseeeedl Appendix VI: Nondimensional longitudinal Equations....... Appendix VII:Aproximation to Short Period and Phugoid Modes ee neers oer eeereaerace Fig.1.a Fig.2.a Fig.3.a Fig.3.b Fig.3.c¢c Fig.4.a Fig.4.b Fig.4.c Fig.4.d Fig.4.e Fig.5.a Fig.5.b Fig.5.c Fig.5.d Fig.5.e Fig.5.f Fig.5.g Fig.5.i Fig.5.j Fig.II.a xil LIST OF FIGURES An Example of Motion Representation..........ee06- 69 Euler Angle Representation... ...csseaes eee ewenee 70 Graphic Setup.........200- cece eee eee eee eee eens 71 Looping Maneuver....... pee ene eee eee e eres ce neces --72 Display Organization... ....2eeeeeeeces Short Period.........eeeaee Prhugoid Mode...secvenacccevcecee sete eens sec e eens 75 Looping and Phugoid...cccccsceccrvcccsceccsereeeece lO Dutch Roll Mode...... iri rie Spiral mode.......... a) Autorotation.....eees rare 6 Nonlinear Lift Curve...... cece eee ec wee cece ee eees 80 Erratic motion........6. cece reece eee eee ener eee e Ol Nonlinear Angle of Attack ResponSe.........cceeeee 82 Phase-Space Representation of Limit Cycle......... 83 Typical Phase-Space Configuration..........eee eee 84 Sensitivity to Control History......... cee eee eee 85 Variable Nonlinear Lift Curve......... cece ee ee eens 86 Nonlinear Rolling Moment Coefficient............6. 87 Study of the Singularity in Pitch........ eee eee .. 88 CHAPTER I INTRODUCTION The potential hazard associated with accidental stalling and Spinning of aircraft has received different attention depending on the type of aircraft. For a military aircraft, the high angle of attack range is part of the maneuver domain, but for a general aviation aircraft the high angle of attack range is certainly not a part of normal operation. However, the safety concern, for all types of aircraft, motivates interest in this area due to the high fatality rate associated with such accidents. In an effort to improve the safety record, which is of prime concern for general aviation aireraft, several controversial issues have arisen, including pilot training; e.g., should spin recovery be part of pilot testing and aircraft certification? Should an aircraft be able to recover from a fully developed spin? Although a lot of work has been done in this area, the stall-spin phenomenon is still poorly understood and seems to be dependent on many parameters [1], [2]. In contrast, the theory of flight dynamics in the low angle of attack range is well developed, and in most cases the calculations are sufficient to give confidence in a specific design. In the late 70's a significant research effort concentrated on that area [3], [4]. An extensive experimental program was launched, which included several standard testing methods ranging from easy to complex. Wind tunnel tests are usually the first to be performed. In conjunction with water channel flow visualizations, they provide an understanding of the general structure of the flow, as well as knowledge of the static stability of the airplane under study. The second stage is usually dynamic testing (forced oscillations){5]. Combined with data obtained with a rotary balance, this set of experiments produces the information necessary for an analytical prediction method. However, the most reliable source of information on stall-spin characteristics, prior to actual flight tests, which are done last, is obtained by testing dynamically scaled airplane models. These models may be dropped from a helicopter or an aircraft , powered and radio-controlled, or flown in a spin tunnel. The most noticeable fruits of this work are the observation of the high nonlinearity of the phenomenon and the strong coupling between the different modes of motion. The inherent nonlinearities of stall behavior prevent the generalization of tests from one configuration to another, and even a minor change in the geometry can induce some drastic effects , e¢.g., pitch-up moment at stall instead of pitch-down moment. No characteristic trends have been displayed by the studies other than the high costs of the experiments and the high sensitivity to many factors. The decreasing cost of numerical studies relative to experimental studies can only strengthen the interest in analytical tools. Two main approaches are taken in the literature to study this problem numerically. In the first approach [6], [7], the forces and moments of the airplane are computed using a table lookup method which yields coefficients based on test data (rotary balance and dynamic testings). The results are reliable, being in close agreement with experiments, and the computations are not complex. The second, and by far the most ambitious approach , aims for a thorough analysis through the use of lifting surface theory [8]. Inputs to the program are typically geometry specifications and a e-D lift curve. This method [9],[10] is complex and requires extensive computation. It is still in a development phase and the assumptions made are very restrictive. However, an obvious advantage is that an airplane need not actually be built before uSing this method. Therefore, it can be used in the design phase. It should be emphasized that the reliability of this analytical tool has not been determined. PRESENTATION OF THE SUBJECT The purpose of this study is to gain a better understanding of the stall-spin phenomenon and its associated nonlinear behavior , by using analysis and interactive 3-D graphics. The originality of this work resides in its use of graphics. Computational output is usually displayed by using plots of the several dependent variables as functions of time. However, with plots, even a slightly complex motion may become obscure and difficult to visualize. An example is given in Figure (1.a): The fuselage of the airplane is, in fact, decribing a cone while the airplane is rolling around its x body axis. The graphic output designed for this study displays a plane "flying" in real time, some flight instrumentation and control input devices. Output files are generated for further study and/or playback . This arrangement allows the user to see the motion as it is computed and to influence it through the controls. The interaction feature is very important in testing a numerical model of the aerodynamics and to quickly point out the interesting characteristics of the flight dynamics associated with it. The code developed here is intended to be more a learning device than a simulator and no special efforts were made to have the model respond as a specific airplane. The underlying idea was inspired by linear aerodynamics for which computations can be made easier by omitting those parameters which had little or no influence. Exploring the relative importance of the parameters through this interactive tool is the main goal of this work. CHAPTER II EQUATIONS OF MOTION The motion of an aircraft is a 6-degree of freedom problem and therefore can be fully described by a set of six nonlinear coupled differential equations of the second order, representing the translational and rotational accelerations of the airplane in a body-fixed coordinate system [11]. du Fs — + wr-oryv ); x 7m (gr 4 ) Fo Svepu-pw); y Ce u- pw) F =~m(a¥+pv-qu); (2.1) d d = —? ar - . Neel, Gtiew a’ Gz Tart Pai y x dt (I, -1,) r p+Iy, p a = dp ar - - No Tye dt * I dt * (Ty I)p q Tye qr. In most airplanes, the y plane is a plane of symmetry; therefore, Tay = Lyx 7 0 The last two terms of each of the equations are the kinematic coupling terms due to the rotation of the axis of the aircraft. These terms represent the inertial nonlinearities in the system. Other sources of nonlinearities are the aerodynamic forces and moments which depend on angle of attack, angle of sideslip, velocity and rotational rates. In order to minimize the complexity of the equations (expressions of aerodynamic forces and moments), another set of variables was chosen. The equations of motion were transformed to (u,8,a) form. (Appendix I) This body-fixed reference system is convenient for describing the aerodynamic forces and moments. However, a representation of the gravity field through Euler angles is then also necessary [12]. The Euler angle representation is as follows (see figure (2.a)). To transform the inertial axis into body axis, the inertial system is first rotated with respect to its Z-axis with an angle jy, yaw angle, then by an angle 6, pitch angle, with respect to the new y-axis and finally, by an angle 6, roll angle with respect to the x-body axis. The transformation matrix is consequently: cosycosé sinycos6é -sing R(v,9,¢) = 4 -sinvcosotcospsinesing cospcos¢+sinysinésing cosésing sinysind+cosycoseosineg -cosysing+sinycosgsiné cosscosq. The order in which the rotations are performed is important. Any vector V expressed in the inertial reference system can be resolved in the body axis system : Vy Vx vy Vo This representation is valid and unique for “3< as 3 O<p Sm and 0K oS The gravity force is along the Z-inertial axis, so G. = (m g cose sing , ~m g sine , m g cosé cos¢ ), (2.3) in body axes. The equations to be solved are : u = x -qutana+ru ae, Q ase [Fcosa - F sina] - p tgBcosa + q ~ Sinatgp r, 2 = Scenes [ cos Fy - singcosa F, - singsina F] - r cosa + p sina, NY =I, oe + 1. s + (I, - 2) qr+l, pq, (2.4) d No= 1, ot CU, -1z)r pty pa No Tyz dt . a dt . (Ty Typ q Tx qr $ =p + tgé (q sind + r cosd), 8 = q cos¢ - r sing, b = sing + cose r cose cose ’ with: Fl=-m g sing - D sina cosg + L sina + T, F y= ™ g cosé sing - D sing, (2.5) Fo= m g cos@ cosg ~- D sina cosg - L cosa. Note here that 6 = q; y r, > =p is true only in the linear range. The equations (2.5) are a complete representation of the mechanics of the system. They have several singularities. Some of them, at 8 = + 90°, u=0Q and a= + 90°, were ignored because of their unusual occurrence. The singularities at 9 = + 90° could not be discarded as they can occur in maneuvers which are of interest to our study, i.e., looping and terminal phase of a spin. This singularity in the pitch angle is not in the physics of the problem but is introduced by the Euler representation: at @= +90°, both roll and yaw angle are referenced with the same axis (the Z-inertial axis is then aligned with the x-body axis) so that the system can detect only o7v. This can also be seen in the transformation matrix R(¢,0=90°,w) for 6 = 90°. ¢ and wy cannot be determined individually; thus, the equations have to be singular. 0 0 cam | R(o,9=90°,¥) = sin(o-w) cos (¢-y) 0 (2.6) cos(o-p) -sin(o-v) 0 An analysis around these singularities was performed and is presented in Appendix II. However,because of a variable time step feature, which checks’ the value of the error, no further modifications of the program are necessary. The differential equation system solver, MODDEQ [13], sucessfully handles those singularities. This software subroutine uses a Runge~Kutta-Gill scheme featuring automatic control of truncation error and variable step size. The code is reasonnably fast ( 1000 time steps require about 15 CPUs). Results are presented and discussed in Chapter IV and Chapter V. CHAPTER III INTERACTIVE DISPLAY 3.1 General description of the PS300 The Evans and Sutherland PS300 is a graphic device which, through a Motorola 6800 microprocessor, can handle all graphic transformations independently from the host computer, in this case, a VAX-VMS; the interface is then best suited for infrequent communications of small amounts of data. The PS300 display rate is about 60 Hertz for a structure of up to 50,000 vectors; thus, the animation of the picture is smooth. A more technical description follows. Nr 2 GCP l=} MM lenl DP rf HOST VAX PS300 - 13-7 GCP : Graphics Control Processor *controls communications with the host *processes commands and creates data structures in the Mass Memor y (MM ) *performs memory management. MM : Mass Memory of 1 Megabyte DP : Display Processor *generates a picture on the screen. Internal communications *Interface 1 is a 16 bit-path. *Interface 2 is an 8 bit-path. External communications *Interface 3 is an asynchroneous line of 19.2 K baud. *Interface 4 is a 16-bit parallel direct memory access 1 mbyte/s 3.2 Operation First an ASCII file is downloaded through the serial interface to the PS300. This file will create all the display structure. for example, it will include a list of vectors representing the airplane as well as a function network. The function network establishes connections between transformations, a rotation, for example, and an object to be transformed (Appendix III). When all the display - 74 - Structures are in the PS300 memory, a FORTRAN program is run on the VAX, which determines the airplane motion by solving the equations. The translations and rotations are then sent down to the nodes of the function network in the PS300 memory, and information such as control surface deflections and thrust setting are requested from the graphic device where they can be set interactivally by turning the control dials. The flowchart below summarizes the organization of the program. At each step of time, the VAX requests the elevator, the rudder, the aileron deflections, as well as the thrust setting from the PS300 using fast subroutines developped at CALTECH by Professor Antonsson; the program calculates then the new position in space and velocities of the airplane and sends translations, rotations and values to the PS300 functions network. t = 0 Request &, 8a, Se, Tr —O elevator from the PS300 <<) aileron Mass —O rudder Memory h ti Compute the motion - ) thrust Send translations & rotations to PS300 re t= t+dt | All the communications between the VAX and the PS300 are through the parallel interface. The VAX program writes directly into the PS300 memory, bypassing the graphics control processor that usually deals with communications with the host. The GCP can then process the display commands faster. Figure (3.a) shows the physical graphic setup. Figure (3.b) represents a looping maneuver through a series of three pictures. On Figure (3.c), a typical display is represented. The biplane is an EVANS and SUTHERLAND design (it was in a demonstration program). On the top left, the control surfaces panel is showing the deflections of the surfaces. On the bottom left , the angle of attack indicator is next to the speed indicator. The attitude instrument on the right shows the pitch angle (distance from the horizontal line), the sideslip angle (distance from the vertical line), and the roll angle (angle between the horizontal line and the symbolized airplane). The number on the bottom far left is the iteration number, giving to the display a time reference. ~ 16 - CHAPTER IV LINEAR AERODYNAMICS 4.1 Linear aerodynamic model The first step of this study is to reproduce the well-known motion of an airplane [11] corresponding to small angles of attack . In this range of angle of attack, the aerodynamic is linear; i.e., the aerodynamic coefficients are expressed as linear functions of the different angles and velocities. The expansions of these coefficients are chosen as follows: L= + Sug (4 7 G0) * “Lge Se, Cy = Cy, + Chy Ca ~ do), Cy B+ Cyep Sry (4.1) Cp = Cyg B+ Cy Pp typ r+ Cygn 6, + Cage ba, m= Smg (4 ~to) + ng A+ Cmge Ser - 17 - néa Sa: Cc = The derivatives are assumed constant. The value of these derivatives for the NAVION and the CHEROKEE 180 are shown in Appendix IV. This representation is very close to reality in the small angle of attack range. When the angle of attack is close to A>, the Stall angle of attack, C;. changes drastically and this a model is no longer valid. 4.2 Longitudinal motion In the equations of motion described in Chapter II , the longitudinal mode can be uncoupled from the lateral modes; i.e., a longitudinal motion stays longitudinal as long as no lateral perturbation occurs [14]. x dt Foemt 2¥-qul] (4.2) 2 dt N =I 9&9, y dt . = de >. Fr z where: q ae and F Fo + GG. These equations are the longitudinal equations of motion. No assumptions have yet been made concerning the mechanics of the system, but since the aerodynamic model limits the problem to small angles of attack , the equations of motion can be linearized. First, they are nondimensionalized; then each variable is expressed as the sum of an equilibrium value and a perturbation part Finally, the new expressions are introduced in the equations where the second-order terms are neglected (see Appendix V for this derivation). The resulting linear system is: du . -_ 2 at Cou U CL cOs8, 8 d(a-8) _ at Coy ur CL sind, 8 (4.3) K 2 _y 2 da? u Co) ge28 = Cp, where: = e psec’ k I e mc? * This motion can be studied analytically . Two characteristic motions are found: the phugoid mode, which is a long period, typically 20 seconds, lightly damped oscillation of the pitch angle and the velocity; and the short period, which is a rapid, heavily damped oscillation acting mostly on the angle of attack. When the short period is excited, an oscillation dies out in about 2 periods, typically 2 seconds. Consequently, this mode is difficult to display. In Figure (4.a) the short period mode is noticeable during the first few seconds, after which the phnugoid mode takes over. However, crude estimates of its period and its damping are 1.4 - 19 - seconds and -2.6, respectively. The period of the phugoid displayed in figure (4.b) is about 30.8 seconds and the damping -0.0155. Another example of a phugoid mode is given Figure 4.c, where the airplane describes a looping before the long period oscillation. The analytical values (Appendix VII) are, for the NAVION [15]: Phugoid: T=30.1 seconds and damping=-0.0149 Short period: T=0.99 seconds and damping=~-2.13. The results from the numerical analysis are close to those of the theoretical. 4,3 Complete motion The full set of equations (2.4) is then used in order to obtain the lateral modes of motion. The numerical model used to get the results described in this paragraph includes’ the complete nonlinearized , 6 degrees-of-freedom equations of motion as well as the linearized aerodynamic model discussed in paragraph (4.1). The program was first sucessfully tested in the case of longitudinal motion, and similar results to those in paragraph 4.2 were found. The lateral motion is characterized by three modes, the roll subsidence, the spiral mode and the dutch roll. The roll subsidence is a stable mode, predominantly a heavily damped rolling motion. The spiral mode consists mostly of a heading change with small roll and sideslip angles. The dutch roll mode is a damped oscillation - 20 - where the three angles 96, 8, w have approximately the same magnitudes. Typically, the time it takes to damp to half-amplitude is approximately one period, In the dutch roll mode, as the airplane yaws to the right, it slips to the left and rolls to the right; then the motion reverses; i.e, it yaws to the left, slips to the right and rolls to the left in a continuous process. For the NAVION, both roll subsidence and spiral mode are stable [ii]. The roll mode is very heavily damped (o=-8.435). Consequently, both of these modes are almost impossible to notice, whereas the dutch roll has a strong influence and is easily characterized (see Figure 4.d). According to this plot, the period of the dutch roll is about 2.75 seconds and its damping is o=-0.422. By linearizing the equations around an equilibrium position and by looking for a nontrivial solution of the form exp(A t), a period of 2.69 seconds and a damping of -0.46 can be found analytically 11,15]. These results are, again, very close. For the Cherokee 180, the roll subsidence is also damped but the spiral mode is a diverging mode [11] and the characteristic change of heading, wy, is shown in Figure (4.e). This finishes a complete analysis of the linear aerodynamic case; the agreement of the numerical and analytical results validates the mechanics of the model. - 24 - CHAPTER V NONLINEAR AERODYNAMICS 5.1 Nonlinear aerodynamic model Near stall, the general picture of the flow becomes more complex with cells of separation, stalled surfaces and vortex arrangement [11]. These phenomena interact in a manner which makes theoretical analysis difficult and the aerodynamic behavior highly nonlinear. A few degrees of change in the angle of attack can trigger drastic changes in the aerodynamic coefficients of the aircraft. Autorotation phenomenon, an interesting behavior associated with the nonlinearities, was first explained by Glauert [16] as follows: For an aircraft flying at an angle of attack a,, a roll rate p can be damped or amplified depending upon the initial angle of attack a,. The angle of attack on the right wing, for a positive rolling moment, increases as the angle of attack of the left wing decreases. If the initial angle of attack is below stall, then the lift curve slope is positive and the difference in local angle of - 29 - attack creates a damping moment; i.e., the lift increases on the right wing and decreases on the left wing. If a, is in the range of the negative lift curve slope, as the lift decreases on the right wing and increases on the left wing, the resulting positive aerodynamic moment amplifies the rolling motion (see Figure 5.a). This can explain how the slightest nonsymmetrical perturbation encountered near stall transforms a symmetrical departure into a spin departure; i.e, the nonsymmetry is amplified by the same mechanism. However, experiments show that the autorotation leads to a steady roll rate [17] [20]. The experimental setup used to get these results restrains the degrees of freedom of the plane to a rotation. In the case of the actual airplane, the strong coupling between the degrees of freedom can be one of the reasons a steady autorotation rate is not reached.Since the spin entry appears to be related to the autorotation phenomenon, an aerodynamic model designed to produce such behavior should have the same instabilities, i.e., a negative lift curve slope and uncoupled wings, allowing the rolling amplification. A simple model gives insight into the driving mechanisms involved. The analysis done with the first model computes the angle of attack on each wing from the average angle of attack and the rolling rate (see Figure 5.a): Qnw = Arctan ( tana + BS ) (5.1) -~ 23 - S1y = Arctan ( tana - a ). The lift and rolling moment are then computed : If -0.2 $a < 0.2 , then Cr = 0.406 + 4.44 (a-ag) + 0.2 64; if 0.2 $a < 0.362, then Ce = (8472.32)? - 8a +3.26 + 0.26, 3 if 0.3628 a , then Ce = 0.56 + 0.2 64; (5.2) if -0.28S a < -0.2, then Cr =(8a+2.32)* - 8a - 2.46 + 0.2 S, i if a < -0.28, then Cc. = -0.219 + 0.2 6. Cy = -0.074 8 + 0.5 (Cr nw - Chiw + 0.01 r+ 0.01 on + 0.13 54 ’ : _- 8$ and finally, L=q 3 ( Cr nw + Chiw): Figure (5.b) represents the lift curve. As long as the system is experiencing a longitudinal motion, the behavior of the airplane is coherent. Some interesting features of this regime are reviewed in the next paragraph. When the angle of attack is greater than ap, i.e., deep-stall range, the airplane is very stable and responds accurately to the controls. On the contrary, a small lateral perturbation occurring in the negative lift-curve range (a, < a < ap) is drastically amplified and induces some violent variations in the angles and velocities (see Figure 5.c). The solutions seem to be reaching a chaotic stage and the - 24 - airplane displayed on the interactive graphics device goes through erratic maneuvers. 5.2 Longitudinal motion The 6-degree of freedom system has some strong instabilities that the longitudinal system does not have. The longitudinal motion computed with this model has some interesting features studied in furtner detail with a simplified 3-degree-of-freedom system including only the longitudinal variables. The system to be solved is then yth order. It is interesting to compare the linearized longitudinal nondimensional equations (2.4) to this set of nonlinear longitudinal nondimensional equations. (Appendix VI contains the derivation. ) Cc ° ; D sina . u = -sine - 2 2 to. u tana; cosa * cos*a Gur 8 ° g = cos(8-a) cosa _ Us sett cosa sina + 8 (5.3) u cos*a a o aaa om u2, Mo V2 = T . . where A =———2 and T=—— and C, , C, as defined in paragraph (5.1). -~ 25 - A typical solution of these equations is shown in Figure (5.d), a plot of the angle of attack variation versus time . The angle of attack diverges first and then reaches a stable oscillatory motion. This behavior is called a limit cycle; A representation in phase-space is given in Figure (5.e), where the angle of attack is plotted versus the pitch angle. For a positive or null lift slope at equilibrium (a = 6 =u =0) the motion is stable; l.e., a trajectory starting in the neighborhood of the equilibrium point will eventually spiral in to this point, but if the starting point is too far from the equilibrium point, the trajectory will reach a limit cycle. Figure (5.f) shows a typical phase-space plot for an equilibrium angle of attack of 12°. The equilibrium point and the unstable limit cycle (dashed line) as well as the stable limit cycle are represented here. When the trajectory starts inside the first curve (unstable limit cycle), it converges to the equilibrium point. When it starts outside this curve it tends to reach the stable limit cycle (the second curve) and stays there. The sequence of stable equilibrium point, unstable limit cycle and stable limit cycle is a typical configuration in a nonlinear dynamical system [18]. It should be noted that for a on the negative lift slope, no unstable limit cycle can be found; i.e, it does not matter how close to the equilibrium point the trajectory starts; it always goes back to the limit cycle. - 26 - Another major observation is that, depending on the control path, the results can be very different even if the same final control setting is reached. In Figure (5.g), two elevator setting histories are displayed as well as the associated angle of attack variations. In the plots on the left, a small positive elevator deflection provides, at first, the system with an excess in forward speed. Then the elevator is set to its final negative deflection; the airplane goes through a series of loopings and seems to reach a steady state in this looping mode. On the contrary, when the elevator is set directly to its final negative value (the same as before), the airplane encounters an oscillation in angle of attack which is amplified until it reaches a steady state. Prior to pursuing this course of study , it is necessary to address the question: " Is this phenomenon relevant to the stall behavior or is it a parasite solution occurring because of an inaccurate model?". No definitive answer can be given. This limit eycle phenomenon, while possibly part of the driving mechanism of stall-spin departure, may not be observed in real tests because the strong lateral instabilities such as autorotation cause a spin departure. However, the model is very simple and consequently may include some major flaws. In particular, the straight Ch curve is an oversimplification that may have led to inaccurate results. Therefore, more complex models have been tested, and these results are presented below. - 27 - 5.3 Other nonlinear aerodynamic models The negative lift slope range is of prime importance to the stall phenomenon; i.e, it triggers the strong lateral instabilities responsible for spin departure, A variable lift curve = as represented in Figure (5.h) has been incorporated in the previous aerodynamic model. The parameter A is the slope of the curve in the stall range. Surprisingly, the results show little sensitivity to this parameter. The limit cycle behavior is not changed. The period of the oscillations changes very slowly. The biggest flaw of the model is then the linear variation of the pitching moment versus a. Some airplanes may experience a pitch-up moment at stall. The loss of efficiency in the controls is also usual because the control surfaces stall when the angle of attack is too high. Both effects are integrated into the model through the pitching moment shown in Figure (5.j). In this figure, the solid line represents the pitching moment coefficient when the elevator setting is zero. The new feature here is mainly the change of slope of this curve around stall. The slope is negative in the linear range and in the deep-stall range but positive for the stall-development range ( a, <a <ap ). The dotted line represents the pitching moment coefficient for an elevator setting of -0.05 radian. For large angles of attack, the dotted line lies on top of the solid line because, at these angles of attack, the elevator surface is stalled and has lost its efficiency. In the linear - 28 - range, the distance between the two lines is constant as predicted by the linear model (see Chapter IV). In the stall-development range, the efficiency of the control surface is decreasing exponentially. With this model, the airplane experiences a stall from which recovery is difficult (the elevator loses almost all its efficiency). The solution still has limit cycles, but the results look quite different. Introduction into the system of some further damping by the means of Ch creates drastic changes; i.e., the limit cycles disappear. As the longitudinal equations of motion are back to a more stable behavior, we can try the 6-degree-of-freedom model. This last nonlinear aerodynamic model was tried on the complete 6-degree of freedom system. The results are promising. At stall, a step-rudder input triggers a spin, but this motion is quickly damped and the airplane reaches a Straight flight path after about half a turn. This behavior shows a lack in the aerodynamic model , more specifically, in the yawing moment. No nonlinearities are included in its expression . The steady developed spin requires a difficult balance between aerodynamic and kinematic moments [1], and an accurate description of these through the numerical method is necessary. CON CLUS TONS This study has shown a high sensitivity of the results to the chosen aerodynamic model, e.g, damping through Coat? ete. This fact seems to imply that a small change in the model might produce a completely different solution structure, and the validity of the idea that a simple aerodynamic model can give insight in the full phenomenon could be questioned. Furthermore, the stall-spin behavior and more specifically the steady developed spin involve a precise balance between kinematic and aerodynamic moments [1]. The kinematics of the system are exactly reproduced in the numerical analysis, but the aerodynamic model remains crude. Its major flaw is the linear yawing moment coefficient. Changing this part of the model is certainly the next logical step in reproducing stall-spin departure and subsequent motion. The numerical analysis developed in this study does not yet provide a satisfactory simulation of the phenomenon, but much progress has been made in this direction. The coding is complete and the aerodynamic model can be changed without changing the main program. The communication problem between the graphic device and the VAX host computer is solved, and a fast information exchange rate is now possible through the parallel interface. The graphic tool is complete and provides the user with a rapid "feel" of the - 30 - characteristic features of the computed motion as well as its validity. This is difficult to do by just using plots. This tool was heavily used throughout the study. The pictures 3.a, 3.b, 3.c¢ are the only output from the PS300 shown in this thesis, not because they are the only relevant things from the graphic display, but because only a movie could reproduce accurately the help gained through graphics and unfortunately a movie cannot be inserted into a book. 10. REFERENCES James S. Bowman Jr. : Summary of Spin Technology as Related to Light General Aviation Airplanes. NASA TN D-6575 = [1971] NASA : 1980 General Aviation Stall-spin Workshop. [1980] J. R. Chambers and S. Grafton: Aerodynamics Characteristics of Airplanes at High Angles of Attack. NASA TM 740097. [1977] Josepn R. Chambers : Overview of Stall/Spin technology. AIAA paper 80-1580. [1980] K.J. Orlick-Ruckemann : On Aerodynamic Coupling Between Lateral and Longitudinal Degrees of Freedom AIAA 15th Aerospace Sciences Meeting.[1977] Joseph R. Chambers, William P. Gilbert : Results of piloted Simulator Studies of Fighter Aircraft at High Angles of Attack AGARD CP-235. [1978] D. Johnston, D. Mitchell, T. Myers : Investigation of High Angle of Attack Maneuver Limiting Factors. TR AFWAL-tr-80-3141. [£1980] W. Bollay : A Nonlinear Wing Theory and its Application to Rectangular Wings of Small Aspect Ratio. Z. angew. Math. Mech. Bd.19 Nr.1. [1939] E. S. Levinsky : Theory of Wing Span Loading Instabilities near stall. AGARD CP-204. [1976] M.A.Hreha and F.H.Lutze : A Dynamic Model for Aircraft 1. 12. 13. 14. 16. 18. 20. - 32 - Poststall Departure. Journal of Aireraft, Vol.21,No.1. [1983] Barnes W. McCormick: Aerodynamics, Aeronautics, and Flight Mechanics. New York, Wiley,{1979]. James Wertz: Spacecraft Attitude Determination and Control. Amsterdam, Reidel,[1980]. Kiko Matsumoto: MODDEQ/ Differential Equation Solver. Mathematical Subroutine Library, California Institute of Technology. [1984] Perkins and Hage: Aircraft Performance Stability and Control. New York, Willey. [1949] G. Teper: Aircraft Stability and Control Data. NASA STI TR 176-1. [1976] H.Glauert: The Rotation of an Airfoil about a Fixed Axis. Advisory Committee for Aeronautics, Rep. Memo. No. 595. [1919] Montgomery Knight : Wind Tunnel Tests on Autorotation and "Flat Spin." NACA Report No. 273. [1927] Guckenheimer and Holmes: Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields. Springer, Verlag, New-York. [1983] H. Lugt : Autorotation. Ann. Rev. of Fluid Mech. 15:123-47. [1983] - 33 - APPENDIX I The (u,a,8) Form of the Equations of Motion We already know that the equations which express the linear velocities are: Foam [ Ms qgw-rv ] x dt Foe m [ Weru-pw] (I-1) y dt Fo =m [ = +pv-qu |. If V is the magnitude of the velocity vector, a the angle of attack and g the sideslip angle, then: u = V cosa cosg v= V sing (1-2) x it V sina cosg. tang cosa Therefore: v= w= tana u. This representation reveals a singularity at an 25 which will not be considered because light aviation aircraft usually do not reach such extreme angles of attack. - 34 - We can then write: d . <y . tang du, 8 #4 tang tana u dt cosa dt cosa cos72 4 cosa dw du G — = tana — + -—=> - dt dt cos?a 7° (T-4) Combining those equations with the above (I-1)equations we get Px _ tang u m q cosa ce i cos : . a= — [ Fo cosa - Fy sina ] - tangcosa p + q - Ssinatang r Wwe cosgpcosa [cosg F - singScosa F, - singsina F - cosa r+ Sina p. mu y 8 x 8 2] P -~ 35 - APPENDIX II Study of the Singularity in Pitch Angle The equations which govern the variation of Euler angles are the following: 8 = q cos¢ ~ r sind ¢ =p + tane (sing q + r cos) (II-1) = sind q+ cos r, cose cosé The equations have several singularities. In order to gain a better understanding of them, we must first remember that they were artificially introduced by the choice of Euler representation and then we must analyze the behavior of the equations around these singularities. Let us consider a pure longitudinal motion. We have p =r =0, and the equations become: 6 = q cos¢ -~ r sing > = tane sing q (II-2) - 36 - p= sind cosé The equilibrium position for » or p is @ = 0 or 3 points where $ = 0 andp~=0. When o =z 0, sind behaves like » and i.e., we get 6 = q tane 6. If both q and 6 were constants, then 9 = exp(q tane t) for small enough $, and » would be a stable equilibrium position as long as 9 <0; i.e., tane < 0 and unstable for @> 0. When @ approaches T 3 a small perturbation around = 0 would have an infinite exponential growth, i.e., taneé + * . In this case, the variation of q and 6 can be considered slow in comparison to ¢. $9 is expected to jump close to the stable equilibrium position, $ =m, when 6 is positive . This, in fact, is the result when the program is run. When 4 is exactly zero, then the pitch angle goes through the value 3. if @ takes on a small nonzero value, the pitch angle remains less than 3» and the roll and yaw angles jump by approximatively mn. Figure (II-a) presents some typical results. - 37 - APPENDIX III Presentation of the Function Network Before any program is run on the VAX, a file containing PS300 commands is downloaded to the PS300. This file defines, in the PS300 memory, objects and transformations. As an example, all the PS300 commands affecting the speed indicator are pointed out in the command file. The Vax is sending data directly into the PS300 memory to speed up the communication process. In order to avoid that a same buffer be accessed by both the PS300 (reading process) and the VAX (writing process), a double buffered structure is adopted. For example , when the VAX writes into the nodes called K***, the PS300 GCP traverses the nodes L*** and the process reverses; i.e., the VAX writes into the nodes L*** and the PS300 reads the nodes K***, The PS300 command file containing all the functions network is shown below: line:=vector list n=2 ~1,0 -0.75,0; line100:=rotate in line200:=rotate in line300:=srotate in -77.14 then line301; line400:=rotate in z ~102.86 then line401; line500:srotate in z “128.57 then line; line600:=rotate in z -154.28 then line601; lLine700:=rotate in z -180 then line; -25.7 then linel01; ~51.43 then line; NNN WN tchl:=begin_ structure character size 0.02 ; translate by -1.2,-0.05; set char screen_ ; character '100'; end_ structure; line101:=instance of line,tchi; tch3:=begin structure - 38 - character size 0.02 ; translate by -1.1,-0.2; set char screen ; character '300'; end_ structure; line301:=instance of line,tch3; tch4:=begin_ structure character size 0.02 ; translate by ~1.075,0; set char screen ; character '400'; end structure; line401:=instance of line,tch4; tch6:=sbegin_ structure character size 0.02 ; translate by -1.075,0 ; set char screen ; character '600'; end_ structure; line601:=instance of line,tch6; linem:=begin structure scale by 0.2; vector list n=2 0,0 -1,0; end_ structure; { The nodes receiving information from the vax are then created } { The double buffered feature can be noticed } k100:=srotate in z 0 then linen; 1100:=rotate in z0O then linem; velol:=instance of line,line100,line200, line300, line400,line500,1ine600,line700; velos:=scale by 0.2 then velol; k101:=instance of velos,k100; 1101:=instance of velos,1100; - 39 - k102:=translate by 0,-0.95 then k101; 1102:=translate by 0,-0.95 then 1101; KOOO:=SET CONDITIONAL BIT 1 ON APPLIED TO ABC; ABC:=BEGIN STRUCTURE IF CONDITIONAL BIT 1 ON THEN KO21; IF CONDITIONAL BIT 1 OFF THEN L021; END STRUCTURE ; KOO1:=ROTATE IN Y O THEN k002; LOO1:=ROTATE IN Y 0 THEN 1002; KO02:=ROTATE IN Z 0 THEN KO003; LO02:=ROTATE IN Z 0 THEN L003; KOO3:=ROTATE IN X 0 THEN fu3; LO03:=ROTATE IN X 0 THEN fu3; KO31:=TRANSLATE BY ~1,0.3,0 THEN KOO1; LO31:=TRANSLATE BY 7~1,0.3,0 THEN LOO1; k008:=instance of kO004,k007,k205; 1008:=instance of 1004,1007,1205; KOO5:=characters 'wait'; LOO5:=characters ‘wait'; K006:=scale by 0.02 then KO005; LO006:=scale by 0.02 then L005; KOO7:=translate by QO. 9, ,O then KOO06; LOO7:=translate by 0.9 ,0 then LOO6; Were) -0 ,70. kOO4:=begin structure a:=field of view 30 front boundary = 2 back boundary = 4; e:=set intensity on 0.75:1; instance of k031; - 4o - end_ structure; 1004:=begin_s a:=field_ of view 30 front boundary = 2 back boundary = 4; e:=set intensity on 0.75:1; instance of 1031; end_s; DISPLAY KOOO; sa oe ae ne a a ee es a ae as ee ee a a ee a i a ee a ww oe ee pri:=f: print; prif:=f: print; send 'reset cs' to <1>flabell; send 'reset th' to <1>flabel2; strin:=characters ~15,0 'AILERON stri2:=characters -15,0 'ELEVATOR stri3:=characters -15,0 'RUDDER stri4:=characters -15,0 'THRUST te ’ . t 1 ° , . ? {THIS ONE IS THE AILERON DEFLECTION} roll:=characters '0'; {THIS ONE IS THE ELEVATOR DEFLECTION} ecora:=characters '0'; {THIS ONE IS THE RUDDER DEFLECTION} ampl:=characters '0'; thr:=characters '315'; st2:=translate by 0,3,0 then str2; 8t3:=translate by 0,6,0 then str3; st4:=translate by 0,9,0 then str4; str2:=sinstance of stri2,cora; str3:=instance of stri3,ampl; stri:=instance of strin,roll; str4:=instance of stri4,thr; String:=instance of stri1,st2,st3,st4; Strinb:=scale by 0.02 then string; strine:=trans by ~0.5,0.75,0 then strinb; display strine; variable rudder ,aileron,elevon, ITERATION, thrust; rot:=f:dxrotate; roy:=f: dxrotate; rou:=f:dxrotate; roth:=f: dxrotate; send 10 to <3>rot; send 10 to <3>roy; send 10 to <3>rou; send 50 to <3>roth; printr:=f: print; printy:=f: print; printp:=f: print; printth:=f: print; COCA1 :=F:CONCATENATEC; COCA2:=F:CONCATENATEC; COCA3:=F:CONCATENATEC; coca4:=f: concatenate; limr:=f:limit; lima:=f: limit; lime:=f: limit; limth:=f: limit; send -15.0 to <3>limr; send -15.0 to <3>lima; send -15.0 to <3>lime; send 0 to <3>limth; send 15.00 to <2>limr; send 15.00 to <2>lima; send 15.00 to <2>lime; send 800 to <2>limth; pri:=f: print; send 0.0 to <1>rudder; send 0.0 to <1>aileron; send 0.0 to <1>elevon; send 315.0 to <1>thrust; SEND 0 TO <1>ITERATION; eonn dials<1>:<1>rot; conn dials<2>:<1>roy; conn dials<3>:<1>rou; conn dials<4>:<1>roth; conn rot<2>:<1>limr; conn limr<1>:<1>rudder; conn roy<2>:<1>lima; conn lima<1>:<1>aileron; - 4o - conn rou<2>:<1>lime; conn lime<1>:<1>elevon; conn roth<2>:<1>limth; conn Llimth<1>:<1>thrust; conn printr<1>:<1>COCA3; conn conn conn SEND SEND SEND SEND CONN CONN CONN conn conn conn conn conn send send send send brou: printp<1>:<1>COCAz2; printy<1>:<1>COCA1; printth<1>:<1>eoca4; ' ' TO <2>COCA1; ' ' TO <2>COCA2; ' ' TO <2>COCA3; ' " TO <2>COCA4; COCA1<1>:<1>AMPL; COCA2<1>:<1>CORA; COCA3<1>:<1>ROLL; cocad<1>:<1>thr; limr<1>:<i>printy; lima<1>:<1>printr; Lime<1>:<l>printp; limtn<i>:<1>printth; "AILERON' to <1>dlabel2; 'ELEVator'to <i>dlabel3; "RUDDER' to <1>dlabel1; 'thrust' to <1>dlabel4; =f: broutec; eqe:=f:eqe; broth:=f: broutec; eqcet: conn conn send send send conn send conn conn eonn send conn =f:eqe; fkeys<1>:<1>eqe; fkeys<1>:<1>peqct; fix(2) to <2>eqet; fix(1) to <2>eqc; 0.00 to <2>brou; eqe<1>:<1>brou; 315.0 to <2>broth; eqet<1>:<1>dbroth; brou<1>:< roth; broth<1>:<2>roth; 315 to <2>roth; brou<1>:<1>rou; - 43 - conn brou<1>:<2>rou; conn brou<1>:<2>rot; conn brou<1>:<1>rot; conn brou<1>:<2>roy; conn brou<1>:<1>roy; ee ee ne ee me ee ee ee me ee ey oe a ee me ee ae ee me me ee re ee ee ee ee ee oe ee ee oe ee ee ee ee ee Oe ee ee ee a on me ye ae ee one ae me ae ae ne me ee ae me ne ee a re ee ee ee ee ee ee ee ee te ee re re ee ee ee Ser:=vector_ list n=10 “0.1,0.2 -0.2,0.1 -0.2,-0.1 -0.1,-.2 0.1,-.2 -2,70.1 -2,0.1 0.1,.2 -0.1,.2;3 gage:=vector list itemized n=2 p .2,0 1 -.2,0 p 0,.2 1 0,-7.2; temoin:=vector list n=4 0.04,0 0.004,0 0,0.01 -0.004,0 -0.04,0; k020:=stranslate by 0,0,0 then k0O10; 1020:=translate by 0,0,0 then 1010; kO22:=translate by 0.75,-0.75,0 then k023; 1022:=translate by 0.75,-0.75,0 then 1023; k023:=instance of scr,gage,k020; 1023:=instance of ser,gage,1020; kO21:=instance of k022,k008,k102; 1021:=instance of 1022,1008,1102; KO10:=ROTATE IN Z O THEN TEMOIN; LO10:=ROTATE IN Z O THEN TEMOIN; at_lin:=vector_list itemized n=6 ~ 0,-1 10,1 “0.1,0.2 1 0.1,0.2 -0.05,0.362 1 0.05,0.362 -0.05,0 1 0.05,0 -0.1,70.26 1 0.1,-0.26; at_rw:=vector list n=1 0.5,0 0,0; at_lw:=vector list n=1 -0.5,0 0,0 ; k206:=translate by 0,0,0 then at_rw; 000 UU - yy - 1206:=translate by 0,0,0 then at_rw; k207:=translate by 0,0 then at_1lw; 1207:=translate by 0,0 then at_1lw; k205:=begin structure translate by -0.8,-0.5,0; scale by 0.2; instance of k207,k206,at_lin; end_ structure; 1205:=begin_structure translate by -0.8,-0.5,0; scale by 0.2; instance of 1207,1206,at_lin; end_ structure; PLANE := BEGIN S VECTOR _LIST BLOCK ITEMIZED N=157 P 0.08936,-0.07446, 0.0 L 0.08936,-0.08191, 0.07446 L -0.05957,-0.08191, 0.07446 L -0.05957,-0.07446, 0.0 L -0.05957,-0.08191,-0.07446 L 0.08936,-0.08191,-0.07446 L 0.08936,-0.07446, 0.0 L 0.08936,-0.08936, 0.14893 L -0.05957,-0.08936, 0.14893 L -0.05957,70.07446, 0.0 L 70.05957,-0.08936,-0.14893 L 0.08936,-0.08936,-0.14893 L 0.08936,-0.07446, 0.0 L 0.08936,-0.09680, 0.22339 L -0.05957,-0.09680, 0.22339 L -0.05957,-0.07446, 0.0 L -0.05957,-0.09680,-0. 22339 L 0.08936,-0.09680,-0.22339 L 0.08936,-0.07446, 0.0 L 0.08936,-0.10425, 0.29785 L -0.05957,-0.10425, 0.29785 L -0.05957,-0.07446, 0.0 L 70.05957,-0.10425,-0.29785 L 0.08936,-0.10425,-0.29785 L 0,.08936,-0.07446, 0.0 L 0.08936,-0.11169, 0.37231 L -0.05957,70.11169, 0.37231 L -0.05957,70.11169,-0. 37231 L 0.08936,-0.11169,-0. 37231 a ae a oe oe ol oft oft oll oft oll oll coll oll oll ll oll oll oll “oll oll dl a a 0 A a 2 - 45 - -08936,-0.07446, 0.0 -08936,-0.11914, 0.44678 -05957,-0.11914, ©. 44678 -05957,-0.07446, 0.0 -05957,70.11914,-0. 44678 -08936,70.11914,-0. 44678 -08936,-0.07446, 0.0 -05957,-0.07446, 0.0 01489, 0.05212, 0.28296 -02979,-0.10425, 0.29785 -O7446, 0.05212, 0.28296 -02979,-0.07446, 0.0 -O7446, 0.05212,-0. 28296 -02979,7-0.10425,-0.29785 01489, 0.04468,-0.28296 -05957,-0.07446, 0.0 10425, 0.07446, 0.05957 -O4468, 0.07446, 0.05957 ~O4468, 0.07446,-0.05957 10425, 0.07446,-0.05957 10425, 0.07446, 0.05957 10425, 0.06702, 0.13403 -O4468, 0.06702, 0.13403 -O4468, 0.07446, 0.05957 ~O4468, 0.07446,-0.05957 04468, 0.06702,-0.13403 10425, 0.06702,-0.13403 10425, 0.07446,-0.05957 ~ 10425, 0.07446, 0.05957 -10425, 0.05957, 0.20850 -O4468, 0.05957, 0.20850 -O4468, 0.07446, 0.05957 ~O4468, 0.07446,-0.05957 -O4468, 0.05957,-0.20850 ~10425, 0.05957,-0.20850 10425, 0.07446,-0.05957 10425, 0.07446, 0.05957 10425, 0.05212, 0.28296 -O4468, 0.05212, 0.28296 -O4468, 0.07446, 0.05957 -O4468, 0.07446,-0.05957 -O4468, 0.05212,-0. 28296 10425, 0.05212,-0.28296 10425, 0.07446,-0.05957 10425, 0.07446, 0.05957 10425, 0.04468, 0.35742 04468, 0.04468, 0.35742 -O4468, 0.07446, 0.05957 ~O4468, 0.07446,-0.05957 -O4468, 0.04468,-0.35742 P L L L L L L L L L L L L L L L L L L L L L L L L L L L L L L L L L L L P L L L - 46 - 0.10425, 0.04468,-0.35742 —0.10425, 0.07446,-0.05957 “0.10425, 0.07446, 0.05957 -0.10425, 0.03723, 0.43188 0.04468, 0.03723, 0.43188 0.04468, 0.07446, 0.05957 0.04468, 0.07446,-0.05957 0.04468, 0.03723,-0. 43188 -0.10425, 0.03723,-0. 43188 0.10425, 0.07446,-0.05957 0.23828,-0.02979, 0.02979 ~0.13403,-0.02979, 0.02979 ~0.46167,-0.02979, 0.0 -0.13403,7-0.02979,-0.02979 -0.05957,70.02979,-0.02979 -0.05957,-0.07446, 0.0 -0.05957,70.02979, 0.02979 -0.05957,-0.02979,-0.02979 0.08936,-0.02979,-0.02979 0.08936,-0.07446, 0.0 0.08936,-0.02979, 0.02979 0.08936,-0.02979,-0.02979 0.23828 ,-0.02979,-0.02979 0.23828,-0.02979, 0.02979 0.23828, 0.0, 0.02979 0.08936, 0.017489, 0.04468 0.13403, 0.01489, 0.04468 ~0.46167, 0.0, 0.0 “0.13403, 0.01489,-0.04468 0.08936, 0.01489,-0.04468 0.23828, 0.0,-0.02979 0.23828,-0.02979,-0.02979 0.23828, 0.04468,-0.02979 0.08936, 0.07446,-0.02979 0.04468, 0.07446,-0.02979 0.06702, 0.14893,-0.06702 0.08936, 0.07446,-0.02979 0.13403, 0.07446,-0.02979 —0.46167, 0.02979, 0.0 -0.13403, 0.07446, 0.02979 0.08936, 0.07446, 0.02979 0.06702, 0.14893, 0.06702 0.04468, 0.07446, 0.02979 0.08936, 0.07446, 0.02979 0.23828, 0.04468, 0.02979 0.23828,-0.02979, 0.02979 —0.13403,-0.04468, 0.0 —0.13403,-0.02979, 0.02979 -0.13403,7-0.02979,-0.02979 ~0.13403,-0.04468, 0.0 rererrerrerr errr rere errr errr eee INST PLANE_1; END S; PLANE 1 := ROTATE IN X 0. THEN PLANE 23 -0.40210,-0. ~0.43188,-0. -0.46167,-0. -0.46167, 0. -0.49890, 0. -0.49890,-0. ~0.46167,-0. -0.46167, ~0.46167, -0.43188, -0.40210, ~0.43188, -0.49890, ~0.49890, ~0.46167, ~0.49890, -0.49890, -0.46167, “0.07446, -0.02979,-0. 0.01489, 0. 0.05957,-0. 0.01489, 0. 0.01489, 0. 0.05957,70. 0.01489, 0. -0.02979,-0. L -0.07446, 0. 0 «8 « = 6 © « 6 oOCo0oO0 0000 0 0 0 0 0 0 Q. 0 0 i) 0 - 47 - 02979, 13403, 13403, 02979, 02979, 13403, 13403, 0.0 0.0 a ee ee ee ee) oOooo°oo°0o eo 2 «© © © O7446,-0. 10425,-0. 05212,70. 10425,70. O7446,-0. O7446, 0. 10425, 0 05212, O74H6, oooo0o°o 0 -0,70.14893 0.14893 0.14893 0.14893 -0.14893 -0.01489 0.01489 0.14893 -O, 0.14893 OQ. 02979 29785 28296 29785 02979 02979 - 29785 0.28296 10425, 0. 0 29785 02979 ; PLANE 2 := VECTOR_LIST BLOCK ITEMIZED N=8 creer oOqooqoo0o 00 0 fu3:=begin_ structure rotate in x 180; scale by 0.4; instance of plane; end_structure; PLANE 3 := F:CLFRAMES; PLANE 4 := F:MODC; PLANE 5 := F:XROTATE; - 23828, 0. -26807, 0. 27551, .26062,-0. - 23828, 0. .26062, 0. -27551,70. -26807, 0. Q. 0, 0.0 0, 0.0 01489, 0.09680 01489, 0.09680 0, 0.0 01489,-0.09680 01489 ,-0.09680 0, 0.0 ; CONNECT PLANE 5<1>:<1>PLANE_1; CONNECT PLANE _ 4<1>2<1>PLANE_ 53 CONNECT PLANE _3<2>: <1>PLANE_ 4; CONNECT PLANE 3<2>: <5>PLANE_ 33 SEND F1X(360.) TO <2>PLANE _ q; SEND FIX(2.) TO <1>PLANE 33 SEND FIX(1.) TO <2>PLANE 3; SEND FALSE TO <3>PLANE 3; SEND FIX(-34.) TO <4>PLANE_ 3; SEND FIX(0.) TO <5>PLANE 3; -~ 4g - - 49 - Appendix IV Linear Aerodynamic Coefficients for the NAVION and the CHEROKEE 180 n (ft?) oC (ft) e (ft) = (1b) m (slugs) I. (slug-ft?) Ty (slug-ft?) I, (slug-ft?) 1 (slug-ft?) Ve, (ft/s) ao (deg) 180 33.4 5.7 2,750 85.4 1,048 3,000 3,530 176 0.6 NAVION [15] mée y8 yéa yér ma =0.406 ah yy =0.355 =0.05 =0.33 =-9.96 ="1.74 =70.564 =0 =0.157 =-0.683 Cay ir léa Cisr np Cp Cc néa Cher =70.074 =-0.410 =0.107 =0,1342 =0.0118 =0.0701 =70.0575 =70.125 =-0.0035 =70.0717 w (ft?) oO (ft) e (ft) = (1b) m (slugs) IY (slug-ft?) IY (slug-ft?) I. (slug-ft?) I, (slug-ft?) Veo (ft/s) a, (deg) 160 25.77 5.25 2,400 TH 170 1,249 1,312 164 - 50 - CHEROKEE 180 [11] mée y8 ydéa yér ma =0.543 =4.68 =0.934 =0.06 20.44 =-7 42 =-2.40 =-0.396 =0.117 =7~0.741 ~ 51 - APPENDIX V Linearized Longitudinal Equations From the general longitudinal equations of motion : Frem{ Seqw], Fy-m({ S-qu], a linearization around a trim position is done . Each term is expressed as the sum of an equilibrium value and a perturbation part. The equilibrium values are noted ( ), as the varying parts are noted ( }'. u = Ug (1+ ul) t w= w' uy Wy=0 Pe =qs ( Cy. + Cy + Cxg, + Cg): Because the linearization is done around an equilibrium position, Cc = Kor Cg, 0. Fo = qs (Cc, + Cog)» where Cy and Cc, are the nondimensional forms of Fy and f.5, and Cus Cog are the dimensionless forms of Gp, and - 52 - miu t Geo. t' = = where t= oq 8° Note that the ( )' terms are dimensionless. This gives the nondimensional longitudinal linearized equations of motion. dt’ ~ “x xg d(a-6) 1 dt! * 2 ( Ce Cog) Estimates of the derivatives. T - D cosa + sina F = 0 ~L cosa - D sina where T is the thrust, L is the lift and D the drag. For a=0 we get: F -L . Therefore, for a propeller aircraft where T u =const, we have: Cy = Cy, ut Cy. a, where C. = -3. Cp and Cy = CG, -Cyy. The gravity field is: -m g siné Ge = fe) -m g cos6 : -~ 53 - Therefore: mg cos6é, 6' ey = 0 -m g sing, 8' . With a small, we can assume that -L, = mg; therefore, Cog = “Cy 6 cos8, @' and Cc, = “Cio sing, 6'. 8 APPENDIX VI Derivation of the Nondimensional Longitudinal Equations. The equations of motion (2.4) are simplified to get the dimensional longitudinal equations by letting 8 = ¢6=yp=pe=re=0. ay _ Fx at a7 au tana j d —a = COSA ([F, cosa - F, sina} + q; dt mu M =I 49, x dt dé at qs where FY = -m g sine - D cosa + L sina + T and F= mg cosé ~ D sina ~ L cosa. d : i “Uo = -g sino -p SOS%, , Sima, i - q u tana; dt m m dq L cosa sina de Gq _ cosa oy cosa | a de. dt u g cos(8-a) mu T mu * dt’ d79 M dt 2 I ~ 55 - In order to nondimensionalize these equations, some reference velocity and time have to be chosen. The reference velocity, V,, is defined by : mg = 5 pS V2, The reference time, t,, is defined by : t, = “ 2 therefore, the reference length has to be: 1, = Mo g * Let ue = th ' t” - © and T = 1, writing again the above ° t, mg equations we get: * Cy ; * sina oa ee at u = ~sing - u?+ u*+T -6u~ tana; cosa cos “a G, * ¥ . g = cos(87a) cosa _ u_ 7" cosa sina + 6 ¥ cosa u * 9=AC_y? 2 where A = moVg y g Dropping the * we recognize equation (5.3). - 56 - APPENDIX VII Approximation to Phugoid and Short Period Modes PHUGOID MODE. Assume that a = 0 and ignore rotational energy and damping in pitch. Equation (4.2) becomes then: Looking for a nontrivial solution of the form exp (\ t) , we find that A has to satisfy the following equation. 2 -~ = 4) 2c. CL Ci = 0 The solutions in i are: J, c2. XU zu? 4 c pe EL 4 ~ hole where Cu = 73 Cy and C,, = -2 Cy. Note that in the preceding equations everything is dimensionless. To get the damping and the frequency in hertz we need to divide \ by 3 Cy g -2 5 AS t, the time scale. We get: 0 = Tare LY ~ 597 - SHORT PERIOD MODE It is mostly a rapid incidence adjustment. Assume a constant forward speed; i.e., u= 0. This implies clearly that 6 =a. Therefore, the angular momentum equation gives: = C @ +C dt 2 Mo, me 9° Looking for a nontrivial solution in A as before gives: ; c C = -_| of Sea. 208 9, A= "TK ome + I sK (Gy): - 58 - APPENDIX VIII Program Listing INCLUDE 'PROCONST.FOR/NOLIST! COMMON /PSBUF/JCT,JBUF(1000) INTEGER*2 JCT, IRATE LOGICAL*1 JBUF,CHAR(8) DATA JCT/1000/ integer*4 IRETRY,ICHAN, IVALI1,IVALI2, IVALI3,ivali4,i INTEGER*2 IADRLO1, IADRLO2, IADRLO3, IADRHI1,IADRHI2, IADRHI3 integer*2 iadrlo4,iadrhni4,iwrite real*4 X(3),x2(3),atr,atl,t REAL Y(9),DETE,TE,yd,ydo,u(3) ,ui(3),a(3,3),xfov,angle REAL YDOT(9),RUD,AIL,ELE,ite,dins(3),thr,clal,x3 CHARACTER*80 DEVSPEC,label,string COMMON rud, ele,ail,thr,clal 0$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$3$$$5$$$5555555$$$$$$355$ $$ EXTERNAL cherokee EXTERNAL FUlila OPEN (UNIT=07,FILE='DEF.IN'’,STATUS='OLD') OPEN (UNIT=23,FILE='command.in' ,STATUS='OLD' ) OPEN (UNIT=22,FILE='input.in',STATUS='OLD! ) read(23,*)string e If idata=0, no output files are generated, if idata=1 then output files c are created. For example attac.out will contain the angle of attack e versus time variations. read(23,*)idata if (idata.eq.1) then eall open end if x2(3)=2 5 PESSSSSSSSSSESSSSSSSSSSSSSSSSSSSSSSSSSSSSS SESSSSSSSSSSSSSSSSSSSSSSS e If the motion is not interactive then the user has decided which - 59 - ¢@ motion to run. According to this choice, the program must run different ¢ sequences. CESSSSSSSSSSSSSSSSSSSS$HSSSSSSSSSSSSHSSS$SH$S$SSS$$$SSSH$SSSS$$SS5S3$$ $ READ(23,*)iinter IF (IINTER.EQ.1) THEN OPEN (UNIT=02,FILE='PLAYBACK.IN' ,STATUS='NEW' ) END IF IF (ITINTER.EQ.10) THEN OPEN (UNIT=02,FILE='playback.in',STATUS='old' ) END IF IF (TINTER.EQ.15) THEN OPEN (UNIT=02,FILE='looppb.in',STATUS='old') END IF READ(23,*)idisplay read(22,*)string READ(22,*)DETE write(*,*)'DETE=', dete READ(22,*)ITER read(22,*)devspec READ(22,*)Y(1) READ(22,¥*)¥(2) READ (22, *)¥(3) READ(22, *)¥(4) READ(22,¥*)Y(5) READ(22,*)¥(6) READ(22,*)¥(7) READ(22,¥*)Y¥(8) read(22,*)y(9) read(22,*)devspec read(22,*)clal read(22,*)ail read(22,*)ele read(22,*)rud thr=315.00 CHSS$SS$SSHSSHSSSSSSSSHSSSSSSSSH$SSSSSSSSSSSSSSSSSSS$SSS$SSSSSSSSSSSS TE=0 K=1 0$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$S$$$S$S5SS$$SSS$$$SS$$S$$$SSH5$$5$$ e If the program is to be run interactivaly, then the VAX and the PS300 ec are linked to allow fast communication rates. ~ 60 - if (iinter.eq.1) then CALL PS3_PHY_ATTACH(ICHAN) CESSSSSSSSSSHSSSH$SH$SH$SSSSSSSSSSSSSSSSHSS$SSSSSHSSHS$SSSSHSHSSHSHHHS$HSSS IRETRY=10 ec In order to achieve a fast communication rate, the GCP is bypassed c and the VAX program will request directly from a location in the PS300 ec memory a value needed in the computation. For that, the VAX has to know e the location in memory of some variables, i.e., rudder deflection. ec Those subroutines were developped at CALTECH by Doctor Antonsson label='rudder' CALL PS3_PHY_GET_ADDR(label,IRETRY,ICHAN, IADRHI1, IADRLO1) label='aileron' CALL PS3 PHY_GET_ADDR(label,IRETRY,ICHAN, IADRHI2, IADRLO2) label='elevon! CALL PS3_PHY_GET_ADDR(label,ITRETRY, CHAN, IADRHI3, IADRLO3) label='thrust' CALL PS3_ PHY GET_ADDR(label, IRETRY, ICHAN, IADRHI4, IADRLO4 ) end if CEHS$S$S$SSSSHSSS$SSSSSSSSSSSSSHSSSHSSS$HSSSSHSSSSHSSSSSSSSSSHSSHESH$SSSSSSSS c Then to send and display information on the screen, some PS300 subroutines e were used . if (idisplay.eq.1) then CALL PSINIT('K', IRATE) end if ifnum=10 if (linter.eq.1) then write(2,10)t,ail,ele,rud,thr end if if (i.ge.ifnum) then ele=0 ail=0 rud=0 c Now that we have the address, it is fast to get a value from ce the PS300 memory. Those subbroutines were also developped by Doctor Antonsson. if (iinter.eq.1) then CALL PS3 PHY READ VAR(ICHAN, IADRHI1,IADRLO1, IRETRY,RUD,IVALI1) - 61 - CALL PS3_PHY READ VAR(ICHAN, IADRHI2, IADRLO2,IRETRY,AIL, IVALI2) CALL PS3_ PHY READ VAR(ICHAN, IADRHI3, IADRLO3, IRETRY,ELE, IVALI3) CALL PS3_PHY_ READ VAR(ICHAN, IADRHI4, IADRLO¥, IRETRY,thr,IVALI4) 10 RUD=RUD*3.1416/180 AIL=AIL*3.1416/180 ELE=ELE* 3.1416/180 write(2,10)te,ail,ele,rud,thr end if if (iinter.eq.10.or.iinter.eq.15) then read(02,10)t,ail,ele,rud,thr end if format (5(x,e12.5)) ifnum=ifnum+10 end if if (iinter.ne.14) then CALL MODDEQ(FUlila,K,9,TE,Y,YDOT,DETE,1.E-4) goto 01 end if CALL MODDEQ(cherokee,K,9,TE,Y,YDOT,DETE,1.E7-4) CESSSSSS$SSSSSSSSSSSSSSSSSSSSSSSSSSSSSIS$SSSSSSSSSSSSSSSSSS$$SSSSSSSS $$$$$$3$ e COMPUTATION OF THE SPACE POSITION OF THE AIRPLANE C$$$$$$$$$$$$$$$$$$$$$$$$3$$$$$$$$$$$$$$$$$$$$$$$$55$$$$$$$$$$$55$$5$ $$3333$$ 01 if (idisplay.eq.1) then YD=-Y¥(4)+TAN(Y(8))*(Y(5) *SIN(Y(7))+Y¥(6) *COS(Y(7))) YDO=¥ (5) *COS(Y¥(7))-Y(6)*SIN(Y(7)) u(1)=y(1) u(2)=y(1)*tan(y (3) )/cos(y(2)) u(3)=y(1)*tan(y(2)) a(1,1)=cos(y(9))*cos(y(8)) a(1,2)=-sin(y(9))*cos(y(7))+cos(y(9))*sin(y (8) )*sinly(7)) a(1,3)=sin(y(9))*sin(y(7))+cos(y(9))*sin(y(8))*cos(y(7)) a(2,1)=sin(y(9))*cos(y(8)) a(2,2)=cos(y(9))*cos(y(7))+sin(y(9))*sin(y(8))*sin(ly(7)) a(2,3)=-cos(y(9))*sin(y(7))+sin(y(9))*sin(y(8))*cos(y(7)) a(3,1)=-sin(y(8)) a(3,2)=cos(y(8))*sinly(7)) a(3,3)=cos(y(8))*cos(y(7)) do 30 j=1,3 ui(j)=#a(j,1)*u(1)+a(j,2)*u(2)+a(j,3) *u(3) x(j)=x(j)+ui(j)*dete/2000 -~ 62 - 30 CONTINUE X3=x3+ui(3)*dete/2000 end if Seaeaeg EFS SS ISI SS SS ISSES SSSSESSSESSESSSSSSSSSSESSSSSSS C WRITE THE RESULTS SeSeSSg a SERS SSS SSS SESS SS SSSSSS SSSESSSSSSSSSS SSIES atr=atan(tan(y(2))+y(4)*8.35/(¥(1)/ (COS (Y¥ (2) ) *COS(Y¥(3))))) atl=atan(tan(y(2))-y(4)*8.35/(¥(1)/ (COS(Y(2))*COS(Y(3))))) if (idata.eq.1) then call write(y,te,x) end if CESSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSHSSSSSSSSHSSSSSSSSSSS $$$$$$$$ C SEND DOWN ALL THE INFORMATION TO THE PS300 CESS$S$SSSS$SHSSSSS$SSSSSSSSSSSSSSSSSSSSSSSSSSHSSSSSSSSS$SSSSSSSSSSESSSS $$$$$$$$$ if (idisplay.eq.1) then x2(1)=x(1) x2(2)=-x(3) x2(3)=-x(2) if (x2(3).gt.4) then x2(3)=1 end if if (x2(3).1t.1) then x2(3)=4 end if xf ov=x2(3)*0.2588 if (x2(1).gt.xfov) then x2(1)#-xfov X2(2)=-x2(2)/2 end if if (x2(1).1t.-xfov) then x2(1)=xfov x2(2)=-x2(2)/2 end if if (x2(2).gt.xfov) then x2(2)=-xfov X2(1)=-x2(1)/2 end if if (x2(2).1t.-xfov) then x2(2)=xfov x2(1)#-x2(1)/2 - 63 - end if x(1)=x2(1) X(2)=-x2(3) X(3)=-x2(2) dins(1)=y(3)/3 if (dins(1).ge.0.08) then dins(1)=0.08 end if if (dins(1).le.-0.08) then dins(1)=-0.08 end if dins(2)=(mod(y(8)+3.1416, 2*3.1416)-3.1416)/5 if (dins(2).ge.0.19) then dins(2)=0.19 end if if (dins(2).le.-0.19) then dins(2)=-0.19 end if dins(3)=0 encode(4,126,char)i 126 format (i4) call tran(206,0,atr,0) call tran(207,0,atl,0) CALL TRAN(031,X2(1),X2(2) ,X2(3)) CALL ROT(001,Y(9),2) CALL ROT(002,Y¥(8) ,3) CALL ROT(003,-Y(7),1) CALL TEXT(005,4,CHAR) eall tran(020,dins(1),dins(2),dins(3)) CALL ROT(010,Y(7),3) call rot(100,-y(1)*3.1416/700, 3) call rot(133,x3*3.1416,3) CALL NUFRAM end if 20 CONT INUE V4 FORMAT(E20.5) 0$$3$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$55$$5$55 $$$$$$$3$$ if (iinter.eq.1) then CALL PS3 PHY detACH(ICHAN) 0a0o0g0ngn0 AacCmNUCAON NO 10 end if END subroutine open ~ 64 - OPEN (UNIT=11,FILE='VELO.OUT',STATUS='NEW' ) OPEN (UNIT=12,FILE='ATTAC.OUT', STATUS='NEW' ) OPEN (UNIT=13,FILE='beta.OUT', STATUS='NEW' ) OPEN (UNIT=14,FILE='rora.OUT', STATUS='NEW’ ) OPEN (UNIT=15,FILE='PIra.OUT',STATUS='NEW' ) OPEN (UNIT=16,FILE='yara.OUT’,STATUS='NEW' ) OPEN (UNIT=17,FILE='phi.OUT', STATUS='NEW' ) OPEN (UNIT=18,FILE='teta.OUT', STATUS='NEW' ) OPEN (UNIT=19,FILE='psi.OUT',STATUS='NEW' ) OPEN (UNIT=20,FILE='verti.out',STATUS='new' ) WRITE(11,*) ‘TITLE WRITE(12,*) ‘TITLE WRITE(13,¥*) ‘TITLE WRITE(14,*) 'TITLE WRITE(15,*) 'TITLE WRITE(16,*) ‘TITLE WRITE(17,*) 'TITLE WRITE (19,*)'TITLE WRITE(18,*)'TITLE WRITE(22,*) 'TITLE return end VELOCITY VERSUS TIME' ANGLE OF ATTACK' ANGLE OF sideslip' roll rate’ PITCH RATE! yaw rate! phi' psi' teta! vertical trajectory' subroutine write(y,te,x) real y(9),te,x(3) WRITE (11,10)TE,Y(1) WRITE (12,10)TE,Y(2) WRITE (13,10)TE,Y(3) WRITE (14,10)TE,Y(4) WRITE (15,10)TE,Y(5) WRITE (16,10)TE,Y(6) WRITE (17,10)TE,Y(7) WRITE (19,10)TE,Y(9) WRITE (18,10)TE,Y¥(8) write (20,10)te,x(3) FORMAT(2(2X,E12.5)) ~ 65 - return end C THIS SUBROUTINE DETERMINES THE AERODYNAMIC OF THE PROBLEM CREE R ERE E REE EEK ER KE I TEE ER IE EER EER KEK Y(1) X(2) ¥(3) Y(4) X(5) X(6) X(7) X(8) X(9) gaARMRAANAANAAA SUBROUTINE FUlila(NN,TE,Y,YDOT) REAL CL,CD,CM,CY,CLL,CN,Q,L,D,SF,M,N,LL,FX,FY,FZ REAL Y(9),YDOT(9),TE,A,B,MU real RUD,AIL,ELE,thr,clal,effi INTEGER NN common rud,ele,ail,thr,clal IS IS IS IS IS IS Is Is Is effis1 THE THE THE THE THE THE THE THE THE FORWARD VELOCITY ANGLE OF ATTACK SIDE SLIP ANGLE ROLL RATE PITCH RATE YAW RATE ROLL ANGLE, PHI PITCH ANGLE, TETA YAW ANGLE, PSI C CL IS THE LIFT COEFFICIENT C HERE FOR EXAMPLE CL ALPHA IS EQUAL TO 4.44 CL*#2*(0.203+2.22*(y(2)-0.0105)+0.1*ELE) IF (y(2).GT.0.2.AND.y(2).LE.0.22) THEN CL=1.247+0. 2*effitele END IF IF (y(2).GT.0.22.and.y(2) .1t.(0.22+0.7323/clal)) then CL=1.247-clal*(y(2)-0.22) END IF if (¥(2).gt.(0.22+0.697/clal)) then cl=0.55+0.2*effitele end if IF (y(2).LT.-0.2.AND.y(2) .GE.-0.22) THEN CL=-0.4529 - 66 - END IF IF (y(2).1T.-0.22.and.y(2).ge.(-0.22-0.2673/clal)) THEN CL=-0.4529-clal*(y(2)+0.22) END IF if (y(2).1t.(-0.22-0.2429/clal)) then el=-0.21 end if 03$$$$$$3$$$$$$$$$$$$5$$$$$$$$5$$$$$$$$$$$3$$$$$$$$$$$$$$5$$$$5$5$5$$ $$$$$$$$$$$ C CD IS THE DRAG COEFFICIENT 333 CD=0.05 + 0.33*(¥(2)-0.00105) C CY IS THE SIDE FORCE COEFFICIENT CY¥=-0.564*Y(3) + 0.157*RUD C CLL IS THE ROLL MOMENT COEFFICIENT CLL=-0.074*Y¥(3)-0.039*%y (4) +0.0101 ¥¥(6)+0.0118*RUD+ 1 0O.1342*AIL C CM IS THE PITCHING MOMENT COEFFICIENT CM=~-0 .683*(Y¥(2)-0.0105)-0.161¥*Y(5)-1.77662*ELE C CN IS THE YAWING MOMENT COEFFICIENT CN=0.0701¥*Y¥(3)-0.00546*Y¥(4)-0.0118*Y(6)-0.0717*RUD-O. O0346*AIL C NOW THAT WE HAVE ALL THE COEFFICIENTS, LET'S COMPUTE THE FORCES. Q=0.5*0.002378*(Y(1)/ (COS (Y(2)) *COS(¥(3))))¥**2 C THE LIFT: L=Q*184*CL C THE DRAG: D=Q*184*CD C THE SIDE FORCE: SF=Q*184¥*CY C THE YAWING MOMENT; N= Q*CN*184*33.4 C THE ROLLING MOMENT: LL=Q* CLL*184*33.4 C THE FORCE ALONG THE X-BODY AXIS: F X=-2750¥*S IN(Y(8))-D*COS(Y(2))*COS(Y(3))+L*SIN(Y(2))+thr C THE FORCE ALONG THE Y~BODY AXIS: FY=2750*COS(Y¥(8))*SIN(Y(7) )-D*SIN(Y(3)) C THE FORCE ALONG THE Z-BODY AXIS: FZ=2750*COS(Y(7))*COS(Y(8)) ~D¥SIN(Y(2))*COS(Y¥(3))-L*COS(Y(2)) - 67 - C NOW, WE WRITE THE EQUATIONS OF MOTION AS DESCRIBED IN APPENDIX. YDOT(1)=Y%(6)*Y(1) *TAN(Y (3) )/COS(Y¥ (2) )-Y¥(5) *Y(1) *TANCY (2) ) i + FX/85.4 CG HERE THE PITCHING MOMENT IS COMPUTED M= Q* 1 84 *CM*5 .7 YDOT (2) =-Y¥(6)*SIN(Y(2)) *TANCY (3) )+¥(5)-Y( 4) *TAN(Y (3) ) *COS(Y¥ (2)) 1 +COS(Y(2))/(Y (1) 85.4) *(-FX*#SIN(Y(2))+FZ*COS(Y(2))) YDOT (3) = cos (y(2))*cos(y(3))*(cos(y(3)) *fy-sin(y(3))*(cos(y(2)) 1 *fx+ sin(y(2))*fz))/85.4/y(1)-y(6) *cos(y (2) )+y (4) *sinly(2)) A= (3000-3530) *Y(5)*Y(6)+LL B=(1048-3000) *Y(4) *Y(5)+N YDOT (4) =A/1048 YDOT (5) =M/ 3000-( 1048-3530) *Y(6)*¥(4)/3000 YDOT (6) =B/3530 YDOT (7)=¥(4)+TAN(Y(8))*(¥(5)*SINCY (7) )+Y%(6) *COS(¥(7))) YDOT (8) =¥(5) *COS(Y¥ (7) )-¥(6)*SIN(Y(7)) ydot (9) =y(5)*sin(y(7))/cos(y (8) )+y(6)*cos(y(7))/cos(y(8)) RETURN END - 68 - FIGURES {rad/sec) roll rate (rad/sec) ~ 69 - 50 ~ yaw rate Qo 4 NS pitch rate 4 4 -50 + 0 10 20 t (seconds) The airplane its describing a coning motion. 200 tag 0 0 10 20 t (seconds) FIG.1.a AN EXAMPLE OF MOTION REPRESENTATION -~ 7O - X,y,2 Inertial axis ——e Xx ” | No angle z y Pe ¥ 4 x @: pitch a 2:31 fe AiVaZ body fixed axis FIG.2q EULER ANGLE REPRESENTATION GRAPHIC SETUP FIG.3.a. FIG.3.b LOOPING MANEUVER tawer wens tsteormm. anenew « On the top left, the control panel indicating the elevator, rudder, aileron deflections as well as the thrust setting. On the bottom left, the angle of attack indicator. On the bottom center a speed indicator in feet per second. On the bottom right an attitude indicator . The symbolic airplane position displays the roll angle (angle with the norizontal line), the pitch angle (distance from the horizontal line) and the sideslip angle (distance from the vertical line). FIG.3.c¢ DISPLAY ORGANIZATION 0.15 0.10 0.05 alpha (rad) 0.00 -Q.05 220 200 180 velolf t/sec) 160 140 -~ 74 - ANGLE OF ATTACK 0.0 1.0 2.0 3.0 4.0 5.0 t (seconds) VELOCITY + { 0.0 1.0 2.0 3.0 4.0 5.0 t (seconds) FIG.4.a SHORT PERIOD teto (rad) velolf t/sec) 0.2 -0.0 -0.2 -0.4 220 200 180 160 140 - 75 - PITCH ANGLE FIG.4.b 30 100 t (seconds) VELOCITY 50 100 t (seconds!) PHUGOID MODE -~ 76 - PITCH ANGLE FIG.4.0a 100 time (sec} LOOPING AND PHUGOID betalft/sec} roll angle (rad) 0.04 0.02 0.00 “0.02 -0.04 0.08 0.06 0.04 0.02 9.90 Angle of Sideslip vom 0 5 10 15 20 25 t (seconds) Roll Angle 0 5 16 15 20 25 t {seconds} FIG.4.d DUTCH ROLL MODE psi (rad) beta (rad) 1.0 0.5 0.0 0.10 0.05 0.00 - 78 - YAW ANGLE 5 10 15 20 35 t (seconds) SPIRAL MODE ANGLE OF SIDESLIP 5 10 15 30 36 t (seconds) FIG.4.e SPIRAL MODE - 79 - AIRPLANE ROLLING AT A RATE P CL TYPICAL LIFT CURVE ee ee 3 > s ANGLE OF ATTACK EFFECTIVE ANGLE OF ATTACK ON RIGHT WING FIG.5.a AUTOROTATION C1 - 80 - 1.5 1.0 0.5 Q\-- -------------_2 BL------ -0.5 -1.0 -0.5 0.0 0.5 alpha (rad) FIG.5.b NONLINEAR LIFT CURVE 250 velocity T u(ft/3s) 200 150 100 50 100 angle of sideslip 10 8 y 8(rad) 0.0 100 time (sec) FIG.5.¢ - 81 - 1.0 0.0 “1.0 ERRATIC MOTION angle of attack a(rad) 0 {0g pitch angle r e(rad) 0 100 time (sec} CRAD) ANGLE OF ATTACK 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 - 8? - JUUUU FIG.5.d 10 20 . 30 40 TIME (SEC) NONLINEAR ANGLE OF ATTACK RESPONSE 50 ~ 83 - 1.5 PITCH ANGLE (CRAD) 1.0 0.5 0.0 -0.5 ANGLE OF ATTACK CRAD) Q.C 0.5 1.0 1.5 FIG.5.e PHASE-SPACE REPRESENTATION OF LIMIT CYCLE - 84 - 1.0 IMIT CYCLE og 4+ STABLE L cyc 0.6 t Q.4 UNSTABLE, LIMIT CYCLE 0.2 + 0.0 pitch anglje(rad) -0.2 + -0.4 if -0.6 8 =-0.0775 0.4 0.6 0.8 1.0 1.2 1.4 velocity : usy FIG.5.f TYPICAL PHASE~SPACE CONFIGURATION 1.6 - 85 - CONTROL INPUT 0.05 1 0.05 T 5, (rad) 8, (rad) 0.00 4 o.c0) «(+ -O.05 + “0.05 G.10 T “0.10 “0.15 “0.15 9 50 0 Time(sec) 50 ANGLE OF ATTACK 0.5 , 1.5 T a (rad) a (rad) 00 tA ‘| 1.0 “0.5 + 0.5 + 1S 0.0 0 Time(sec) 50 0 Time(sec) 50 TWO DIFFERENT HISTORY IN CONTROL SETTING BUT SAME INITIAL CONDITIONS FIG.5.¢ SENSITIVITY TO CONTROL HISTORY 1.5 on 1.0 0.5 0.0 -0.5 -1.0 ~ 86 - -0.2 0.0 0.2 0.4 angle of attack (rad) FIG.5.4 VARIABLE NONLINEAR LIFT CURVE - 87 - 0.05 0.00 5, 270.05 -0.05 “0.15 -0.20 “0.25 0.0 ANGLE OF ATTACK CRAD) 0-8 FIG.5.j NONLINEAR ROLLING MOMENT COEFFICIENT ~ 98 - 4 ROLL ANGLE YAW ANGLE "| . ; 2 4 time (sec) 0 time (sec) | 100 0 100 In plot 3, the motion in purely longitudinal and the pitch angle goes through the singularity. In plot 1,2,3, a small lateral perturbation was input in the system. As expegted, the yaw and roll anglesare jumping of about nr. 8 eC PITCH ANGLE 6 i PITCH ANGLE 4 T 3 2 4 is} 4 -2 0 100 0 100 time (sec) time (sec) FIG.II.a STUDY OF THE SINGULARITY IN PITCH