Acta Polytechnica https://doi.org/10.14311/AP.2021.61.0378 Acta Polytechnica 61(2):378–390, 2021 © 2021 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague SIMULATION OF A MANEUVERING AIRCRAFT USING A PANEL METHOD Pavel Schořa,∗, Martin Kouřila, Vladimír Daněkb a Aeromobil R&D, Pripojna 5, 821 06 Bratislava, Slovak Republic b Brno University of Technology, Faculty of Mechanical Engineering, Institute of Aerospace Engineering, Technická 2896/2, 616 69 Brno Czech Republic ∗ corresponding author: schor85@email.cz Abstract. We present a method for numerical simulations of a maneuvering aircraft, which uses a first-order unsteady panel method as the only source of aerodynamic forces and moments. By using the proposed method, it is possible to simulate a motion of an aircraft, while the only required inputs are geometry and inertia characteristics, which significantly reduces the time required to start the simulation. We validated the method by a comparison of recordings of flight parameters (position, velocities, accelerations) from an actual aerobatic flight of a glider and the results obtained from the simulations. The simulation was controlled by deflections of control surfaces recorded during the actual flight. We found a reasonable agreement between the experimental data and the simulation. The design of our method allows to evaluate not only the integral kinematic quantities but also instant local pressure and inertia loads. This makes our method useful also for a load evaluation of an aircraft. A significant advantage of the proposed method is that only an ordinary workstation computer is required to perform the simulation. Keywords: Flight mechanics, aerodynamics, simulation, panel method, maneuvering aircraft. 1. Introduction The numerical simulation of a maneuvering aircraft is not a computationally demanding problem, when one uses a concept of aerodynamic derivations [1, 2]. First computer programs solving the equations of motion were developed in 1970’s [3]. These programs used the concept of aerodynamic derivations, which provide a very fast computation of forces acting on the aircraft. However, to use the aerodynamic derivations, a signif- icant number of dependencies must be determined in advance, mostly by experiments in wind tunnel or by computational fluid dynamics methods (CFD) such as [4]. Instead of pre-computing number of tables by CFD, we propose to couple the flight simulation program directly to a first-order unsteady panel method. The main advantage of this approach is that the free stream and wake always exactly match the current situation, therefore, proper forces and moments can be com- puted. Also, only the geometry of the aircraft must be known. Motions of control surfaces are introduced by morphing of the mesh. 2. Methods 2.1. Flight mechanics In order to simulate a flight of an aircraft, we solve its equations of motion in a non-inertial reference frame located at the centre of mass of the aircraft. We use a flat-earth approximation, as the simulated flight is supposed to be close to the ground. There is an interaction between the flight simulator and the code for the aerodynamic analysis. The interaction is explained in Section 2.3 bellow. This work was developed earlier by the first author as a part of his PhD thesis [5]. Reference frames We first define reference frames which are used during the simulation, as shown in Figure 1. Three different reference frames are used: Earth-fixed - E is an inertial reference frame. The origin is located at an arbitrary point, preferably on the surface of the Earth. The Z-axis points upwards (outside of the earth), Body-fixed - Bf is used for solving of the Equations of motion - Equation 1 and Equation 2. Its origin is always located at the instantaneous centre of mass of the aircraft. Its axes are aligned to axes of Struc- tural reference frame S, as shown on Figure 1. When viewed from origin of the Earth fixed reference E , the Body-fixed reference frame bf translates its origin and rotates its axes as the aircraft moves. Body-fixed Structural - S is used to define the ge- ometry and to compute the aerodynamic forces. It’s origin is defined by a geometry master model, the X-axis points to tail, the Y-axis points to the right wing, and the Z-axis points upwards. Equations of motion The motion of the aircraft is described by its state vector [v, ω, v̇, ω̇]bf The equa- tions of motion are solved in a non-inertial reference frame bf . We adopted the derivation of the equa- tions of motion in a non-inertial reference frame bf by 378 https://doi.org/10.14311/AP.2021.61.0378 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en vol. 61 no. 2/2021 Simulation of a maneuvering aircraft using a panel method Figure 1. Reference frames [5] Stevens [1] and we only present the final Equations 1 and 2. The force vector Fbf in Equation 1 consists of aerodynamic and propulsion forces, which are related to the instantaneous centre of mass, and the gravity force transformed into the Body-fixed reference frame bf . The same applies for moment Mbf . The aero- dynamic forces and moments are evaluated at each time-step for given flow conditions and given control vector ∆C. This control vector contains the position of control surfaces, propulsion control and position of flaps and other means of control, if applicable . v̇bf = Fbf m − ω̃bfvbf (1) ω̇bf = ( Jbf)−1 [Mbf − ω̃bfJbfωbf ] (2) From Equations 1 and 2, we can compute velocities for the next time step by the explicit Euler’s method, see Equations 3 and 4. vbf t+∆t = ∫ t+δt t v̇bf (t)dt .= v̇bf t ·∆t (3) ωbf t+∆t = ∫ t+δt t ω̇bf (t)dt .= ω̇bf t ·∆t (4) As the new state vector [v, ω, v̇, ω̇]bf is known, the new position of the aircraft in the Earth-fixed reference frame E is computed by transforming the velocity vector vbf into earth coordinates E and performing the explicit Euler integration, as shown in Equation 5: pE t+∆t = pE t + ∫ t+δt t vE(t)dt .= pE t + vE t ·∆t (5) The Equation 5 gives a new position of the reference frame bf expressed in earth fixed coordinates E. Also, the reference frame must be rotated to its new angular position. This is done by explicit Euler integration in Equation 6 : φE t+∆t = φE t + ∫ t+δt t ωE(t)dt .= φE t + ωE t ·∆t (6) The Equation 6 gives a new angular position of the reference frame bf expressed in earth fixed coordinates E at time t+ ∆t. Therefore, a new vector of gravity force inbf may be expressed. The actual transformations of vectors in Equa- tions 3, 5 and 6 are done by a multiplication of the vectors by corresponding Directional cosine matrices (DCMs). 2.2. Aerodynamic model Following the requirements from Section 1 and 2.1, we had to implement a very fast solver of the aerody- namic problem, which is capable of computing, within seconds, different combinations of angle of attack, side slip, angular velocities and deflected control sur- faces. All of these requirements have to be fulfilled automatically, without the need to manually alter the fluid domain (e.g., manually changing the mesh). Our choice was a Boundary Element Method, namely Unsteady Panel Method. A detailed description of this method can be found in the literature, see refer- ence [6], or [7]. Therefore, we only recapitulate the fundamentals of this method. Flow governing equations The flow is evaluated in the non-inertial reference frame S, see Figure 1. The flow is assumed to be unsteady, incompressible and irrotational. Such a flow is called a Potential flow. This flow is governed only by continuity equation 7 [6] : ∇2Φ = 0 (7) Following boundary conditions are applied: On the solid boundary, the normal component of the velocity must be zero: (∇Φ + v) · n = 0 (8) In Equation 8, the kinematic velocity vector v rep- resents the velocity of the fluid when viewed from the non-inertial reference frame S. This velocity is decom- posed into components shown in Equation 9, where vrel represents time derivative of body deformation. v = −[V0 + vrel + Ω× r] (9) Far from the body, the disturbance created by the motion should decay: lim r→∞ (∇Φ− v) = 0 (10) Also, the velocity vector v at each point can be obtained by a derivation of the velocity potential Φ: v = ∇Φ (11) A general solution of Equation 7 can be obtained when the surface of the body SB is modelled by source and doublet elements on the boundary of the body. Consider a cross section of a wing as shown on Fig- ure 2. The surface S∞ encloses the problem at infinity, surface S represents a body (wing) and surface W represents its wake. The surface S + W divides the 379 P. Schoř, M. Kouřil, V. Daněk Acta Polytechnica Figure 2. Potential flow domain domain into two regions: the external region with the flow field of interest and velocity potential Φ and the internal region with a fictitious flow and velocity po- tential Φi. The surface S is modelled by a doublet and source singularities. The wake region W is modelled by doublet singularities. The following applies at the boundary between Φ and Φi: − µ = Φ− Φi (12) − σ = ∂Φ ∂n − ∂Φi ∂n (13) Using Equation 12 and 13, the velocity potential at point P can be expressed as: Φ(P ) = 1 4π ∫∫ S+W µn · ∇1 r dS − 1 4π ∫∫ S σ 1 r , dS (14) Where r is the distance from the point P to the element dS on the surface, and n is a unit normal vector of the element dS. A detailed derivation of Equation 14 can be found in reference [6] or[7]. Next, the internal Dirichlet boundary condition is introduced. This boundary condition sets a constant value of the internal potential Φi at the surface of the body: Φi = const = 0 (15) Numerical solution of flow equations The sur- face of the aircraft is divided intom body panels and n wake panels. Each body panel has a constant distribu- tion of source singularities and constant distribution of doublet singularities, both of unknown strength. Following the procedure from [6, 7], or [8], the surface is discretized into n surface panels and nw panels in the wake region with a constant distribution of singularities. The Equation 14 is also discretized and the Dirichlet boundary condition is evaluated at the centroid of each surface panel. These points are called Collocation points. The result of this procedure is a set of linear equation [6] : n∑ k=1 Ckµk + nw∑ l=1 Clµl + n∑ k=1 Bkσk = 0 (16) Where Ck, Cl are doublet influence coefficients of body and wake panels, Bk is source influence coef- ficient, µk is doublet strength of body panel, µl is doublet strength of wake panel, σk is a source strength of body panel. The σk is set as: σk = −nk · [V0 + vrel + Ω× r] (17) Where nk is the panel normal unit vector. The last step is to compute local panel forces: First, local panel velocities are computed as derivates of dou- blet strengths between neighbourhood panels. Then a total local panel velocity and panel pressure coefficient are computed. Finally a local panel force is obtained [6]. Wake Modelling In the panel method, the wake is an integral part of the solution of the flow field. The wake is used to satisfy the Kutta-Joukowski condition at trailing edges. The wake is usually modelled as infinite thin vortex sheets using doublets panels with a constant intensity. To obtain the correct solution of the circulation, it is necessary to know the strength of doublet wake wake panels and their spatial locations. There are two main issues related to the wake: First: If the wake panels are very close to other body control points (e.g. fuselage or tail units) or penetrate it, an unrealistic pressure distribution is obtained on the body. This is due to very high velocities induced by the singular wake panels. This consequently leads to unphysically strong doublets on the penetrated body. A usual remedy to this issue is to attach the wake to the fuselage so that the edges of wake panels are coincident with the edges of the body panels. Second: the wake itself affects the flow field by its downwash. For instance, a downwash due wing wake may have a significant effect on the lift of the horizontal tail and, consequently, also on the pitching moment of the aircraft. As the wake is located closer to the horizontal tail, the more it affects the lift of the horizontal tail. During manuevering of the aircraft, the spatial position of the wake changes significantly. For example, during the pull-up maneuvre, the wake passes above the horizontal tail, as shown in Figure 3. However, if the aircraft flies inverted at a negative angle of attack, the wake passes bellow the horizontal tail as shown in Figure 4. Based on the statements above, for a simulation of a manuevering aircraft, it is required to have a wake model, which can move freely in the flow field and can pass close to or penetrate another objects. So far, 380 vol. 61 no. 2/2021 Simulation of a maneuvering aircraft using a panel method Figure 3. Wake above the horizontal tail Figure 4. Wake bellow the horizontal tail a few methods, which address these issues have been developed. For example Willis [9] solved the wake-body im- pingement by computing intersections between the wake and the body. The intersections formed a closed regions inside the body, where an additional potential was introduced to cancel the initial potential jump at the body-wake intersection. Also, the wake may be modelled by vortex particles, which allows to obtain a detailed wake shape without impingement issues [10]. We chose a method developed by Gennaretti and Bernardini [11]. They split the wake into near-wake and far-wake regions. The near wake is formed by standard doublet panels and their influence is evalu- ated directly in Equation 16. The far-wake panels are treated as vortex rings elements (with Rankine vortex core model [6]) and their influence is evaluated by altering body-panel source strengths in Equation 17 so that the Equation 17 becomes Equation 18. σk = −nk · [V0 + vwakefar + vrel + Ω× r] (18) A setup for our case is shown in Figure 5: The near-field wake is shown as red panels, the far-filed wake panels are white. Using this setup, the wake may penetrate the empennage without introducing unphysical forces on the empennage panels. 2.3. Coupling scheme A coupling between the flight simulation solver and the aerodynamic solver is straightforward, with following steps: (1.) Transform velocities vbf and ωbf from bf to S (2.) Update the control vector ∆C (3.) Compute aerodynamic forces and moments and assemble Fbf and Mbf Figure 5. Wake model: near and far fields (4.) Update the state vector and compute new posi- tion of the aircraft in E (5.) Update wake in E and transform its shape to S Wake model Earlier formulations of the unsteady panel method such as in [6] suggested to solve the wake shape in an inertial reference frame. The body panels are displaced to the new position and already existing wake panels keep their positions. To fill the gap between the trailing edge and the first wake pan- els, new wake panels are added and wake relaxation routine is applied. We found this unsuitable for our application, as adding new panels would result in an increased computational time for long-lasting simu- lations. We developed a modification of the original algorithm by [6] with respect to the far-field wake described above. Our algorithm can be described as follows: (1.) The aerodynamic problem is solved in the struc- tural reference frame S. Coordinates of body panels are fixed, except of moving control surfaces. (2.) The position of the structural reference frame S is computed with respect to the inertial coordinates E (3.) The position of far-wake panels in the inertial coordinates E is computed (4.) Equations of motion are solved and a new posi- tion of the structural reference frame S in inertial coordinates E is computed (5.) For each row of far-wake panels: coordinates of i- th panel become coordinates of (i+1)-th panel from the previous step (6.) For each row of far-wake panels: coordinates of first panels are computed from the position of last near-field panels and second far-field panels (7.) Displacement vectors for each far-field panel are transformed from inertial coordinates E to struc- tural coordinates S (8.) The wake is deformed according to the computed displacement vectors and wake relaxation routine This loop is repeated at each time step. The main advantage is that the number of panels is constant between each time step, therefore, long-lasting ma- neuvers can be computed with a minimum amount of memory. Also, the control surfaces can be precisely 381 P. Schoř, M. Kouřil, V. Daněk Acta Polytechnica Figure 6. Mesh morphing - elevator deflected, there is no distortion due numerical errors as the number of iterations increases. The disadvantage of this is that the physical dimen- sions of far-wake panels change with the change of the airspeed. Therefore, the wake shape may become rough and inaccurate if the number of wake panels is low and airspeed is too high. As can be noted from Figure ??, the wake-update routine does not apply to the near wake, whose shape remains fixed. This is to avoid the body penetration by the near wake, which would result into unrealistic body forces. However, the fixed shape of the near wake reduces the accuracy of the computed forces. 2.4. Control surfaces To control the simulation, it is necessary to simulate deflections of all control surfaces. We used the mesh morphing approach shown in Figure 6. The coor- dinates of panels belonging to a control surface are rotated around its axis of rotation and correspond- ing wake panels are translated according to a new position of the trailing edge. As shown in Figure 6, we use dense panelling at leading edges and coarse panelling at trailing edges. This reduces the risk of panel overlaps at large deflections of control surfaces. Viscous drag The potential flow theory can prop- erly predict only induced drag, a viscous drag can not be predicted. With few modifications, the vis- cous effects can be introduced to the panel method. For example, see the work of Cebeci [12]. Using this iterative boundary layer approach, it is possible to compute the viscous drag and lift reduction correctly due to boundary layer separation. However, this also increases the computational time. If the aerodynamic polar is measured, it is possible to compute a difference CDvisc between the measured drag and drag computed by the panel method. Using this difference, an increment of viscous force Dvisc can be computed as shown in Equation 19. This force is then added to Equation 1. Dvisc = v |v| 1 2ρ|v| 2 SCDvisc (19) 3. Comparison of simulation to flight data 3.1. Overview Flight measurements of aerobatic manoeuvres were done on the two-seat, all-metal sailplane L-13 AC Blanik within the project at the Institute of Aerospace Engineering, Faculty of Mechanical Engineering, Brno University of Technology. 3.2. Numerical simulation To perform a numerical simulation of an aerobatic flight, we had to create a virtual model of L-13 AC sailplane. There are two main purposes of the virtual model. First, an external shape of the sailplane must be discretized into quadrilateral panel for the panel method as described in Section 2.2. Also, hinge axis lines for each rudder must be determined. Second, inertia matrix related to the actual centre of gravity is required by Equation 2. Based on the actual sailplane, we reconstructed its geometry in CAD software. We also modelled all spars, bulkheads and ribs. We also measured the mass of an empty sailplane and the mass of each wing and fuselage. Using this measurement, we corrected the average area density in the CAD. Next, we measured the mass of the pilot (including the parachute) and the data acquisition unit. Then we computed the actual position of the centre of gravity and the actual inertia matrix related to this point using the CAD. The computed inertia matrix is given in Equation 20: J = 1785.11 0 0 0 1327.30 0 0 0 3326.18  (20) It is also necessary to set a proper time step in Equation 5 and 6. A shorter time step would result in a long runtime of the simulation, a longer time step would result in decreased acuraccy and possi- bly non-physical results. We assume that during the sumulation, the sailplane should not travel more than one chord length per time step. With respect to the maximum allowed speed VNE = 64[m/s], we set the time step δt = 0.001[s]. 382 vol. 61 no. 2/2021 Simulation of a maneuvering aircraft using a panel method Figure 7. Pitot-static tube and wanes 3.3. Flight test arrangement We equipped the sailplane with a data-acquisition unit. This unit recorded the following data: • position of control surfaces • attitude angles • linear acceleration (load factors) • angular rates • angle of attack and sideslip angle • airspeed • altitude and outside air temperature The sampling rate of data acquisition was 100Hz for all recorded data. The pitot-static probe with wanes for the reading of the total and static pressure and incidence angles (angle of attack and sideslip angle) was fitted to the glider nose, as shown in Figure 7. The sensors of aileron, rudder and elevator de- flection were installed on longitudinal and lateral- directional control system in the rear pilot seat area. This arrangement may cause errors due to elastic- ity of the control cables, however, we neglected this fact since other arrangements would be much more complicated and heavier [13]. The data logging and control unit along with the static pressure transducer (flight altitude determination), differential pressure transducer (flight speed determination) and optical gyroscope (reading of load factor, attitude angles and angular rate) were installed behind the back pilot seat. The power supply of the optical gyroscope and data logging and control unit were attached in the glider baggage room. Both installations can be seen in Figure 8. The outside air temperature sensor was placed on the lower side of the left wing to eliminate sunlight influence on the sensor. The whole measurement set was activated, its functionality was checked and the calibration of some sensors was done (deflection of control surfaces, temperature sensor, etc.). Figure 8. Sensors board installation and power supply 3.4. Aerobatic flight The actual aerobatic flight took place at airfield Brno-Medlanky (ICAO code: LKCM) on the 22nd of September 2005. The sailplane was towed to 1800m above sea level by a tow airplane. Then, an aerobatic sequence was performed by the third author. This sequence included aerobatic figures according to the FAI (Federation Aeronautique Internationale) aero- batic catalogue, including loops, stall turns, inverted flight, rolls, spins, half-roll-inverted flight-half-loop. A barograph of whole flight is shown in Figure 9. 3.5. Comparison of measured data to simulation To compare our method to experimental data, we selected a sequence of three consecutive loops, which begins at time t = 1360s and lasts 60s. Please note that time t = 0s in Figures 10, 11, 13, 14, 15, 16, and 17 corresponds to time t = 1360s on the barograph in Figure 9. Initial conditions for the simulation were set ac- cording to time t = 1360s. This means that control surfaces were deflected according to t = 0s in Figure 10 as explained in Section 2.4. Inputs for aerodynamic solver were set according to the measured airspeed, angle of attack and sideslip angle and measured an- gular velocities. Also, measured linear and angular accelerations were passed to Equations 5 and 6. The simulation was driven only by control inputs shown in Figure10. As we neglected the stiffness of control cables, the deflections of control surfaces in the simulation exactly copied the control inputs from Figure 10. Figure 11 shows a comparison of computed and measured angle of attack related to XY plane in the structural reference system S. The physical sensor of angle of attack (vane) was calibrated, so that it reads zero when the vane is parallel to the plane. 383 P. Schoř, M. Kouřil, V. Daněk Acta Polytechnica Figure 9. Barograph Figure 10. Control inputs for aerobatic sequence 384 vol. 61 no. 2/2021 Simulation of a maneuvering aircraft using a panel method Figure 11. Comparison of flight data - angle of attack Figure 12 shows a comparison of computed and measured sideslip angle related to the plane in struc- tural reference system. The physical sensor of sideslip angle (vane) was calibrated, so that it reads zero when the vane is parallel to the plane. Figure 13 shows a comparison of magnitude of com- puted and measured airspeed. The computed airspeed comes directly from Equation 5. Figure 14 shows a comparison of computed and measured altitude. Please note that the reference altitude z = 0m corresponds to the altitude at time t = 1360s from Figure 9. Figure 15 shows the pitch angle measured by the optical gyroscope compared to pitch angle obtained from Equation 6. Note that this angle is related to the inertial reference frame E. Figure 16 shows the normal load factor nz compared to the computed load factor. This load factor is related to the structural reference system S. Figure 17 shows a comparison of measured and computed trajectories. The computed trajectory is obtained by integrating Equation 5 and a transfor- mation from the body fixed frame bf to the inertial frame E. The measured trajectory was computed by second and third author in their previous work [14]. 4. Discussion Background information The presented method has been originally developed for evaluation an of loads on Aeromobil 4.0 XPF experimental aircraft shown in Figure 18. During its design process, there was a strong need for a method capable of evaluat- ing local pressure loads over a whole surface of an arbitrary shaped aircraft. Consider the empennage of the Aeromobil 4.0 with it’s panelling representation as shown in Figure 19. There are two horizontal sta- bilizers and three vertical fins. As designed, three is a strong mutual influence of each empennage surface. To compute it’s aerodynamic loads, following options are available: • Develop, validate and solve a set of analytical equa- tions • Use an unsteady vortex lattice method • Use an unsteady 3D panel method • Use an unsteady RANS CFD solver Considering the number of load cases to be evaluated, computational costs and accuracy, it has been decided to develop an unsteady 3D panel method coupled to flight mechanics solver, as presented in this paper. Since the data from Aeromobil 4.0 are confidential at the moment, the authors decided to adopt an earlier work of second and third author [14] to validate the method. Obtained results From Figures 10, 11, 13, 14, 15, 16, and 17, one can see a possible limitation of our method. During the first loop, the sailplane actually stalled at the inverted phase of the loop. During this stall, the angle of attack was α = −24◦, while the predicted angle of attack was only −7.5◦. This dis- crepancy may also be caused by some atmospheric turbulence, however the flight was conducted very early in the morning. Therefore, a presence of a gust is highly unlikely, but it cannot be denied. The predicted trajectory of the first loop in Figure 17 cor- 385 P. Schoř, M. Kouřil, V. Daněk Acta Polytechnica Figure 12. Comparison of flight data - sideslip Figure 13. Comparison of flight data - airspeed 386 vol. 61 no. 2/2021 Simulation of a maneuvering aircraft using a panel method Figure 14. Comparison of flight data - altitude Figure 15. Comparison of flight data - pitch angle 387 P. Schoř, M. Kouřil, V. Daněk Acta Polytechnica Figure 16. Comparison of flight data - load factor nz Figure 17. Comparison of flight data - spatial position 388 vol. 61 no. 2/2021 Simulation of a maneuvering aircraft using a panel method Figure 18. Aeromobil 4.0 XPF during takeoff Figure 19. Aeromobil 4.0 XPF empennage panelling. Note the deflected elevator on top horizontal surface responds to the stall during the inverted phase, but this stall was not reconstructed by second and third author in their work [14]. However, panel methods give a realistic lift, pitching moment and induced drag only for such an angle of attack where no significant separation is present. The remaining two loops show good agreement between the measured and predicted data. There are discrepancies in side-slip angle β shown in Figure 12. The computational analysis pre- dicts more side-slip oscillations, while the physical data seem to be more damped. Figure 10 can shed some light on this discrepancy.. There are three signif- icant rudder inputs at approx time 17s, 33s and 47s. From Figure 12, it is obvious that the model reacts to these inputs in the same way as the physical sailplane. However, it must be pointed out that the sensor of the rudder deflection was located on the rudder’s control cables. The stiffness of control cables is neglected in the numerical model. Therefore, what is presented as a rudder control input may not represent the actual deflection of the rudder during the flight. A slightly higher airspeed was predicted by our method as can be noted from Figure 13. This discrep- ancy could be attributed to damping of the airspeed measuring apparatus, as this was calibrated only to static cases. 5. Conclusions In their paper, Johnson et al. [15] expect that the Navier-Stokes CFD codes may serve as a “routine tool for the load analysis, stability and control analysis and high-lift design process”. We have demonstrated that, under certain assumptions (incompressible, in- viscid fluid and negligible slipstream due propeller), an unsteady panel method can also serve as a reliable tool for stability and control analyse. To support this conclusion, we present a comparison of simulated and actual flight data in Section 3. The comparison of the predicted flight parameters and the measured data gives a very good agreement. The difference between the predicted height loss and the actual height loss is less than 20 meters after finishing the third loop. The predicted load factors are in good agreement with the measured data. The presented simulation lasted more than 60s, which is much longer than simulations required for certification purposes. Considering the load analysis, the presented method provides pressure distribution on whole surface of the aircraft. This gives two advantages. First, a detailed pressure dis- tribution and detailed information of inertia forces is available, making a room for a structural optimization and eventually lighter structure. Second, the pressure distribution can be mapped into a structural finite element analysis. It can be concluded that our presented method may be proposed mainly for load evaluation of virtual prototypes of newly designed aircraft. By using our method, the designers have a tool, which predicts a pressure distribution on a surface of a whole aircraft for a wide range of maneuvers. The applicability of the method is not limited to conventional designs, but it is also possible to evaluate non-conventional designs as shown in Figure 19 or as proposed by Kroo [16], or non-conventional wing tip devices [17]. 6. Future work As a future work, we would like to validate the method on a different type of sailplane, performing different complex maneuvers. Between each maneuver, there shall be a short period of steady level flight, so that we can compare a difference between the simulation and the measured flight data for each maneuver. This is also needed, if restarts of the simulation are required, because these periods of short level flight allow to set accurate initial conditions for the simulation. The method of solving equations of motion shall be im- proved by an integrator with adaptive time. Also, the aerodynamic solver should take into account a velocity of deflecting control surfaces in Equation 9. 389 P. Schoř, M. Kouřil, V. Daněk Acta Polytechnica 7. Acknowledgement This research was funded by Research contract Nr.0685/2015 by Ministry of Education, Science, Re- search and Sport of the Slovak Republic as part of state aid measures for research and development from the state budget under the Act No. 185/2008 Coll. on research and development incentives, supported research and experimental development projects of key components of innovative means of transport for road and air transport in Aeromobil R&D,s.r.o.. The authors are also thankful to Institute of Aerospace Engineering, Brno University of Technol- ogy for providing the flight measurement datasets recorded under Aerospace Research Centre, Grant project LN00B051 supported by Ministry of Educa- tion, Youth and Sports, Czech Republic List of symbols ∇ Gradient operator Φ Velocity potential φ Attitude angles p Position vector n Normal vector u Velocity vector ∆C Control vector E Earth-fixed reference frame BF Body-fixed reference frame S Structure-fixed reference frame References [1] B. L. Stevens, F. L. Lewis, E. N. Johnson. Aircraft control and simulation: dynamics, controls design, and autonomous systems. John Wiley & Sons, 2015. [2] J. Roskam. Airplane flight dynamics and automatic flight controls. DARcorporation, 2003. [3] L. T. Nguyen. Simulator study of stall/post-stall characteristics of a fighter airplane with relaxed longitudinal static stability, vol. 12854. National Aeronautics and Space Administration, 1979. [4] J. Haider, C. H. Lee, A. Gil, et al. A first order hyperbolic framework for large strain computational solid dynamics in openfoam. In World Congress on Computational Mechanics and The Asia-Pacific Congress on Computational Mechanics, pp. 688–688. 2016. [5] P. Schor. Load State of an Aircraft with an Elastic Wing. Ph.D. thesis, Brno University of Technology, Technicka 2, Brno, 2018. [6] J. Katz, A. Plotkin. Low-speed aerodynamics, vol. 13. Cambridge university press, 2001. https://doi.org/10.1017/CBO9780511810329. [7] B. Maskew. Program VSAERO Theory Document: A Computer Program for Calculating Nonlinear Aerodynamic Characteristics of Arbitrary Configurations. NASA contractor report. National Aeronautics and Space Administration, 1987. [8] D. Ashby. Potential Flow Theory and Operation Guide for the Panel Code PMARC. National Aeronautics and Space Administration, 1999. [9] D. J. Willis. Implementing automatic-wake body intersections for full aircraft congurations in dirichlet panel methods. In 31st AIAA Applied Aerodynamics Conference, p. 2650. 2013. [10] D. J. Willis, J. Peraire, J. K. White. A combined pfft-multipole tree code, unsteady panel method with vortex particle wakes. International Journal for numerical methods in fluids 53(8):1399–1422, 2007. [11] M. Gennaretti, G. Bemardini. Novel boundary integral formulation for blade-vortex interaction aerodynamics of helicopter rotors. AIAA journal 2007. https://doi.org/10.2514/1.18383. [12] T. Cebeci, J. Shao, F. Kafyeke, E. Laurendeau. Computational Fluid Dynamics for Engineers: From Panel to Navier-Stokes Methods with Computer Programs. Springer, 2005. https://doi.org/10.1007/3-540-27717-x. [13] J. Bakunowicz, R. Meyer. In-flight wing deformation measurements on a glider. The Aeronautical Journal 120(1234):1917–1931, 2016. https://doi.org/10.1017/aer.2016.98. [14] V. Danek, M. Kouril, R. Sosovicka. Measurement of aerobatic flight characteristics. vol. 3. Czech Aerospace Research Centre, 2006. [15] F. T. Johnson, E. N. Tinoco, N. J. Yu. Thirty years of development and application of cfd at boeing commercial airplanes, seattle. Computers & Fluids 34(10):1115–1151, 2005. [16] I. Kroo. Nonplanar wing concepts for increased aircraft efficiency. ’VKI Lecture Series’ 2005. [17] N. N. Gavrilovic, B. P. Rasuo, G. S. Dulikravich, V. B. Parezanovic. Commercial aircraft performance improvement using winglets. FME Transaction 43(1):1–8, 2015. https://doi.org/10.5937/fmet1501001g. 390 https://doi.org/10.1017/CBO9780511810329 https://doi.org/10.2514/1.18383 https://doi.org/10.1007/3-540-27717-x https://doi.org/10.1017/aer.2016.98 https://doi.org/10.5937/fmet1501001g Acta Polytechnica 61(2):378–390, 2021 1 Introduction 2 Methods 2.1 Flight mechanics 2.2 Aerodynamic model 2.3 Coupling scheme 2.4 Control surfaces 3 Comparison of simulation to flight data 3.1 Overview 3.2 Numerical simulation 3.3 Flight test arrangement 3.4 Aerobatic flight 3.5 Comparison of measured data to simulation 4 Discussion 5 Conclusions 6 Future work 7 Acknowledgement List of symbols References