Identification and Modeling of the Airbrake of an Experimental Unmanned Aircraft

This paper presents the modeling, system identification, simulation and flight testing of the airbrake of an unmanned experimental aircraft in frame of the FLEXOP H2020 EU project. As the aircraft is equipped with a jet engine with slow response an airbrake is required to increase deceleration after speeding up the aircraft for flutter testing in order to remain inside the limited airspace granted by authorities for flight testing. The airbrake consists of a servo motor, an opening mechanism and the airbrake control surface itself. After briefly introducing the demonstrator aircraft, the airbrake design and the experimental test benches the article gives in depth description of the modeling and system identification referencing also previous work. System identification consists of the determination of the highly nonlinear (saturated and load dependent) servo actuator dynamics and the nonlinear aerodynamic and mechanical characteristics including stiffness and inertia effects. New contributions relative to the previous work are a unified servo angular velocity limit model considering opening against the load or closing with it, the detailed construction and evaluation of airbrake normal and drag force models considering the whole deflection and aircraft airspeed range, the presentation of a unified aerodynamic - mechanic nonlinearity model giving direct relation between airbrake angle, dynamic pressure and servo torque and the transfer function-based modeling of stiffness and inertial effects in the mechanism. The identified servo dynamical model includes system delay, inner saturation, the aforementioned load dependent angular velocity limit model and a transfer function model. The servo model was verified based-on test bench measurements considering the whole opening angle and dynamic load range of the airbrake. New, unpublished measurements with gradually increasing servo load as the servo moves are also considered to verify the model in more realistic circumstances. Then the full airbrake model is constructed and tested in simulation to check realistic behavior. In the next step the airbrake model integrated into the nonlinear simulation model of the FLEXOP aircraft is tested by flying simulated test trajectories with the baseline controller of the aircraft in software-in-the-loop (SIL) Matlab simulation. First, the standalone airbrake simulation is compared to the SIL results to verify flawless integration of airbrake model into the nonlinear aircraft simulation. Then deceleration times with and without airbrake are compared underlining the usefulness of the airbrake in the test mission. Finally, real flight data is used to verify and update the airbrake model and show the effectiveness of the airbrake.

During the planned test flights the aircraft should fly a racetrack pattern (see Fig. 3) with two turns and two straight sections. One of the straight sections is the test leg where the aircraft should speed up to a given higher reference speed and then also slow down. As it has a BF Turbines BF300 jet engine with slow dynamics (5s run-up time from idle to full power) a high bandwidth additional actuator is required to make higher decelerations possible. To satisfy this requirement an airbrake was designed [16]. As this will be a critical component of control during flutter tests the precise knowledge of its dynamics is crucial. That's why two test benches were developed by Technical University of Munich (TUM). One to test the whole airbrake mechanism under load including the servo motor, the opening linkages and the airbrake with artificial aerodynamical load [9]. The second is to test the dynamics of the servo motor with different load levels [4]. A preliminary identification of the system dynamics was done in [9] which is extended and refined with several new contributions in this article. Complete identification of the airbrake requires to mathematically model all of the static (such as the ratio from airbrake torque to servo torque) and dynamic characteristics (such as servo motor dynamics). Considering the literature there are limited sources about servo dynamics identification and neither of them considers full airbrake dynamics together with aerodynamic loads. [3] deals with the identification and control of a direct current (DC) motor arriving to a first order plus integrator plus delay transfer function model. Multi parametric optimization is used to obtain the unknown parameters and uncertainty ranges are also considered. [5] reduces the problem of alternating current (AC) motor system identification to the identification of two unknown parameters as others are known. They are determined through a genetic algorithm. [13] identifies first the linear characteristics of the related mechanism then it identifies a 4th order linear model and a nonlinear model considering servo speed limitation. It finally reduces the required five unknown parameters to three and identifies them with least squares optimization. The nonlinear model proves to better approximate the system. [8] deals with identification and control together. It identifies the DC servo dynamics with third order dynamics plus delay and then a controller is designed. Finally, [18] identifies the dynamics of a DC servo motor considering its second order dynamics and the possible inner controller structures also as usually off-the-self servo motors include control electronics too. Finally, it derives a two transfer functions (TFs) model relating reference and rotation angles to servo torque.
In our case several factors should be considered to have a realistic model of the airbrake and several improvements are offered relative to previous work of the FLEXOP team. This is summarized below. as the literature sources also show. First a servo test bench was developed and measurements with doublet series servo angle references on different frequencies and with different load torque levels were carried out (see [4] and [9]). In our case system identification is complicated by saturation in the inner servo controller (control saturation also considered in [13]) and load dependent maximum angular velocities of the servo which differ if servo travels against or with the load (similar load issue is considered in [18]). This issue is partly covered in [9] but a simpler and unified model is derived in the current work. Pulling out all the nonlinearities and system delay (also identified) finally, a transfer function model of the servo is identified using the output error method [11]. As the measurements show a load dependent steady-state servo angle change in the system output (missing I term in servo control) a separate transfer function model was identified for this and the need for this term verified based-on real flight behavior. This is also a new contribution. 2. The nonlinear aerodynamics is covered by [9] only for 60m/s maximum flight speed and only for the normal force. Regarding that model the addition of the drag coefficient and the drawing of the opening angleairspeed-servo load characteristic is the contribution of this paper. 3. Regarding the nonlinearity in the opening mechanism described in [9] (examined also in [13] for another system) the contribution of the current paper is the construction of a unified aerodynamic plus mechanic nonlinear model to determine a direct relation between airbrake angle, dynamic pressure and servo torque. The precision of this model is also analysed. 4. A final challenge is the stiffness of the mechanism well analysed in [9]. Our contribution here is the transfer function-based frequency domain analysis of system dynamics considering stiffness and airbrake inertia.
The complete servo model was verified based on the test bench measurements (adding new variable load measurements to the previous work) and a full airbrake model was created after including servo motor model, opening mechanism and aerodynamic nonlinearity, stiffness and load dependent steady-state servo angle change. This model was first tested for realistic behaviour in Matlab simulations then integrated into the nonlinear simulation model of the FLEXOP aircraft and tested with software-inthe-loop 'flights' of the test pattern with baseline trajectory tracking controller (see [14]) applying airbrake opening during aircraft deceleration after the test leg. Finally, real flight measured position data of the airbrake is compared to the simulation model (driven by in-flight logged dynamic pressure and deflection command) output to verify the identified model. The airbrake effect on aircraft deceleration was also examined based-on the change in specific energy of the system.
The structure of the current paper is as follows. Section 2 briefly introduces the test aircraft and the design of the airbrake then Section 2.1 summarizes the overall concept of airbrake test and identifiaction. After that Section 3 gives brief information about the test benches. Section 4 presents the identification of the nonlinear airbrake dynamics, Section 5 introduces the full airbrake model with SIL simulation results in Section 5.1 and real flight test results in Section 5.2. Finally, Section 6 concludes the paper.

Airbrake Conceptual Design
Before discussing airbrake conceptual design the FLEXOP demonstrator aircraft and its operational circumstances should be briefly introduced to clarify the need and the requirements related to the airbrake. The general configuration of the FLEXOP demonstrator is of a conventional design asserting close similarity to state-ofthe-art commercial aircraft. An illustrative three view of the actual demonstrator is depicted in Fig. 2 together with main technical data in Table 1 (regarding the flutter speed see [17]). The high aspect ratio wings feature a moderate leading edge sweep to account for the characteristic bendtwist-coupling of high-subsonic planform designs. The final goal with the demonstrator is to prove possibility of flying beyond flutter onset speed applying active control of the wing dynamics. To isolate wing aeroelastics from thrust interferences and aerodynamic disturbances, an engine installation separated from the wings (on top of the fuselage) was chosen. Also a V-Tail for low-interference drag, as well as for minimum wetted area fuselage was applied. The required disturbance tolerance covers tolerance of moderate turbulence and deterministic wind disturbances until 5-10m/s as the 65kg mass of the aircraft make it a small one relative to the usual GA or larger aircraft. In future flutter control tests as turbulence and wind free weather as possible is required to remove additional disturbances of the structural dynamics. As the test airfield for the demonstrator is München Oberpfaffenhofen all operations should obey German airspace regulations. According to these regulations flight testing within the visual line of sight of a human pilot is required. The flight altitude is limited to a maximum of 305m above ground level due to the procedures of the airfield operator. A maximum distance between pilot and vehicle in the order of 1000m was considered to provide the required visibility (considering the 7m wingspan of the demonstrator) and so a flight test pattern as illustrated in Fig. 3 is assumed. In the figure the position of the pilot is marked with a red dot, with the resulting straight leg for aeroelastic (flutter) test phases marked by a red arrow [15].
Assuming 38 − 40m/s cruise speed and considering the 1000m maximum distance first estimations indicate an available straight distance of 1500m for the acceleration to test speed, test of flutter controller and deceleration to a safe airspeed for a turn maneuver. According to preliminary studies on flutter onset, an airspeed range up to 60 m/s is required. To account for the resulting requirements on the acceleration and deceleration performance of the vehicle, propulsion alternatives and aerodynamic configurations were studied in a point mass simulation, as illustrated in Fig. 4.
Regarding the engine selection reciprocating, turbojet, turboprop and electric propulsion concepts were studied. Optimization towards minimum system weight (including fuel or battery), complexity and costs finally led to the selection of a micro turbojet engine. The BF Turbines BF300 engine features a nominal thrust of 300N at ISA sea level conditions. The run-up from idle to full throttle consumes 5s. Despite of the unfavorable time constant and the increased consumption, the excellent thrust-to-weight ratio of the engine itself compensates for the additional fuel mass and slow throttle-up however, this arose the need for an additional deceleration device.
The airbrake design was optimized based-on the simple point mass simulation of the demonstrator (see Fig. 4) using a thrust dynamics model based on actual engine measurement data. Historical, empirical data was used to determine preliminary size and estimate drag force for the airbrake. To avoid aerodynamic disturbances or structural buffeting response, wing mounted devices are not applicable to the present problem -although spoilers can provide twice the effectiveness of fuselage installations (see [6]). Thus feasible fuselage locations were evaluated and maximum available airbrake areas determined. An integral solution, merging ventral brake flap and the landing gear door, as well as a double flap system, arranged symmetrically to the fuselage aft of the wing were studied in detail targeting to have a drag force vector going through the center of gravity of the aircraft (see Fig. 5).
With preliminary, conservative assumptions for the brake drag and the actuation speed, a minimum required brake flap area was determined with the help of the point mass simulation of the vehicle. This area estimation indicated integration problems for the ventral solution, since no ground clearance margin was left for the landing gear configuration. So in the next step trade studies for flap aspect ratio, reference area and actuator type were evaluated  for the side-by-side flap configuration. A simple actuator model assuming either a linear or a rotary actuator was derived based on manufacturer specifications. Dead time and actuator dynamics were neglected. A simple transmission ratio rule was used for the preliminary models. The transmission concept is shown in Fig. 6.
Both linear and rotary models were integrated into the point mass vehicle simulation to evaluate the best fitting concept. At maximum airspeed, the total deflection time accounted to 3.3 s for the linear actuator and below 0.5 s for the rotary actuator. In consequence, only a rotary actuator could provide the required short deceleration distance as rendered previously (for details see [16]). So, finally the side-by-side flap configuration and a KST X30-12-150 rotary servo actuator were selected as the main components of the airbrake system. After implementation of the whole system in a mock-up (see Fig. 8) and satisfactory test results installation into the demonstrator fuselage was done. However, to study effects on system dynamics and to tune on-board controllers a mathematical simulation model of the airbrake was required the construction and identification of which is the main topic of this article as follows.

Airbrake System Identification Concept
The identification of the airbrake is a complex task as several parameters and characteristics should be determined. The final construction of the airbrake consists of the servo motor, the opening mechanism and the airbrake itself (see Fig. 37). Going from the inner part to the outer first, the dynamics of the servo motor was determined by building a servo motor test bench (for details see Section 3.2) with a load motor which can be configured to apply different torque loads (consant / varying) on the shaft of the servo. Different dynamic deflections of the servo with different loads were done to cover the whole operational range (see Section 4.1 for the identification of servo dynamics).
After identification the servo model was verified based-on test bench measurements. Second the aerodynamic characteristics of the airbrake were determined including the generated drag force which slows down the aircraft and the load torque on the airbrake shaft resulting from the aerodynamic load. As wind tunnel testing is time consuming and expensive and there is exhaustive literature about such characteristics [6,12] the aerodynamics was finally determined based-on published formulas (see Section 4.2). Finally, the static characteristics were determined such as the ratio of the mechanism from airbrake torque to servo torque and the stiffness of the whole mechanism as this can have a significant effect on the dynamics. The ratio was determined based-on geometrical measurements, however the stiffness required to build a complete airbrake mock-up including the servo motor, the opening mechanism and the airbrake itself with the possibility to add different weights to simulate the aerodynamic loads (for details see Section 3.1). Measuring the deformations with different loads gave the opportunity to determine the opening angle dependent stiffness of the mechanism (see Section 4.2 for the static characteristics). An additional task was to determine the dynamic effect of opening mechanism stiffness and airbrake flap inertia on the overall dynamics. After constructing and identifying all of these model components and building the whole airbrake model the next task was to verify it in Software-in-the-loop (SIL) simulation and on real flight data. This is summarized in Section 5.

Actuator Test Benches
To make determination of airbrake static and dynamic characteristics possible two test benches were built one is a full size mock-up which also serves to test the whole airbrake opening mechanism together with the servo motor before installing it on the aircraft. The second is a separate test bench for the servo motor applying an electric load motor to simulate airbrake loading on the servo axle.

Full Size Airbrake Mock-Up
A mass-loaded mock-up of the airbrake installation is used to statically determine the full system stiffness over the range of deflection angles, as well as the static electric power consumption. The mock-up consists of the actual levers and rods, the airbrake flap, the servo actuator as well as a mounting rig. The servo actuator is powered by a laboratory power supply, limited to 20A. The voltage is set to a constant value of 12V . The stiffness of the rig is designed such that the angular error induced by the deflection of the assembly under maximum load accumulates to 0.0023 • giving minimal distortion in the deflection measurement of the airbrake itself. A Digital Pitch Gauge (DPG), as well as a potentiometer is used for the determination of the deflection. As the maximum airspeed induced by the angular velocity of the airbrake flap is negligible relative to the 38-60 m/s speed of the aircraft no unsteady correction of the aerodynamic forces was required to be applied. So a simple mass carriage/swing was used to solely apply a force normal to the airbrake. Within the quasi-linear load range up to 40 • , a mass of 12.658kg is used to simulate the aerodynamic load. For the range up to 60 • , the mass ballast was adapted to match with the higher aerodynamic load (for details see [9]). The drawings of the mock-up are shown in Fig. 7 while the photos can be seen in Fig. 8.

Servo Actuator Test Bench
In order to assist the identification of the airbrake servo actuator performance and allow for the derivation of a mathematical model representing its operating characteristics, a dynamometer test bench was developed [4]. It offers two main testing modes: a force-free configuration (see Fig. 9) where actuator movement is monitored in the absence of an external load and a torque-loaded configuration where a counter load is applied in a controlled manner by a load machine (see Fig. 10). In the first mode the servo actuator output shaft is connected to an angular position encoder the latter possessing negligible friction and moment of inertia. In the second mode the servo actuator output shaft is connected to a load machine with a torque sensor placed inside the load path between the two units. The next part focuses on the measurement campaign with the latter configuration ( Fig. 10).
Linear systems are often represented by means of a transfer function or a state-space model and offer the convenience, that classical model-based control algorithms can be applied on them. However, in the present case, several nonlinearities are present in the servo actuator e.g.: the built-in angular position controller with current and therefore angular velocity limitation and friction in the airbrake actuator gearbox which dissipates kinetic energy nonlinearly. For this reason a simple linear model representation of the airbrake actuator dynamics should be avoided. The implemented test scenarios are therefore tailored to the actual operation profile of the airbrake, which is designed to deflect against the airstream around the fuselage. Combined with the propulsive engine, it is used to help aircraft deceleration during the mission therefore dynamic deflections are expected in flight. Three test campaigns are conducted with the loaded test bench ( Table 4).
The results of the first two campaigns (Test campaign I & II) and a detailed description of the test setup are found in [4] while results of Test campaign III are newly generated applying gradually increasing load as the servo  Table 2 for different loads and the same deflection command (97.5 • ). Later Test campaign III data is applied to verify the realistic behavior of the servo model.
While the described test bench offered useful experimental data satisfactory to identify the basic mathematical model of the airbrake servo certain limitations are to be pointed out and discussed: 1. Due to the hinge kinematics of the airbrake mechanism, radial forces are exerted against the actuator shaft, in addition to the torsional moment. However, at the test bench the artificial load applied is solely the torsional moment, the radial forces are thus neglected. Depending on the level of the latter, additional friction on the actuator shaft bearing will arise. Successful validation of the identified airbrake model in flight tests (see Section 5.2 and Fig. 55 for example) has shown that there is no significant such effect in the real system. 2. The load torque control algorithm of the load motor in the test bench is based on a feedforward path. Although actual torque is available in real time, addition of a torque feedback control path is avoided, due to instabilities encountered. They are probably attributed to a slight delay in the load motor torque signal and the high stiffness of the load path. As the identified servo model describes the dynamics of the servo motor not the load motor its not possible to examine this effect with the identified model. Another issue is that at the beginning of the actuator movement a peak is visible at the torque profile (see Fig. 31) and also the current profile. This should be caused by the feedforward control experiencing a step change in the reference signal and so commanding a very high input at the first time. This assumption is verified by the behavior for gradually increasing load torque command as there is no such peak in the measured torque as Fig. 34 shows. 3. Sampling is conducted with a frequency of 200 Hz.
The maximum variation that can be accurately captured is therefore less than 100 Hz. Current data however suggests that fluctuations of the actual signal might be more rapid. However, as the operating range of the airbrake will be at low frequency (bang-bang control with full opening or closing) the exact representation of very high frequency dynamics is not required. As Fig. 28 shows the servo motor is incapable to make full range deflection even at 5Hz without load and limited range deflections ( Fig. 29) are also reduced by applying load on the servo. So there is no possibility to drive the servo and so the airbrake with high frequency reference inputs and so the limitation of the represented measured frequency range to 0-100Hz is more then acceptable. This is also underlined by the Bode diagram ( Fig. 11) of the identified servo motor dynamics (series connection of G ref (z) reference signal dynamics from Eq. 1 considering the maximum A N (T L ) = 6.1319 gain if T L = 8Nm, G sys (z) servo dynamics and the integrator) also underlines this with 6.2rad/s = 0.987H z cut-off frequency and 8.85rad/s = 1.4H z bandwitdh. The maximum A N (T L ) gain is larger than the maximum A P (T L ) and the maximum gain gives the maximum bandwidth that's why this is considered.
System identification is discussed in the next sections including reproduction of measurement results of the servo test bench and SIL simulation of the whole nonlinear aircraft dynamics with flight test pattern tracking baseline control without and with airbrake application during  deceleration. Finally, flight testing will be applied for airbrake model verification.

Identification of Airbrake Characteristics
After the design and construction (in the mock-up test bench see Fig. 8) of the airbrake this section deals with the identification of its dynamics. This is started with the identification of the servo dynamics based-on measurements applying servo test bench with load (see Fig. 10) and the verification of the resulting model comparing model outputs with test bench measurements. This servo model can be considered as the inner core of the whole airbrake model. To construct the whole airbrake model the characteristics of the opening mechanism and the aerodynamic load effects should be determined (based-on [9] but extending and improving the results) and integrated considering also the effect of airbrake inertia and the stiffness of the mechanism on the dynamic behavior.

Identification of Servo Dynamics
The selected rotary servo actuator is the KST X30-12-150 servo (see [16] and [9]). Its factory characteristic can be found in [10] page 3 and is summarized in Table 3.
To identify the servo dynamics the test bench presented briefly in Section 3 and in detail in [4] was applied making measurements with different frequency doublet servo angle setpoint inputs and different loads. The parameters of the whole measurement campaign are summarized in Table 4. In the first column the pulse width modulated (PWM) and angular values (α) of servo deflection references are shown. The second column shows the related airbrake angle deflection (φ) and the further columns show the applied torque load. Note that the upper opening limit of the airbrake is 60 • so the 67.9 • maximum value is a bit higher. Because the servo motor is tested independently on a separate test bench there is no limitation of this angle. The header lines of the blocks show the test signal frequencies. 1 − 5H z with 1H z steps (Test campaign I) and 0H z static tests were done (Test campaign II). The third block (Test campaign III) shows a test with step change of the servo reference angle from 0 to 97.5 degrees with continuously changing load to different maximum values as if were in case of inflight opening of the airbrake. Considering Fig. 39 shows that not all measurements are relevant regarding airbrake dynamics as on low opening angles the servo torque is also low. Relevant cases are shown with orange color in the table.
Some illustrative measurement plots are shown in Figs Figures 13 and 15 show that the possible angular velocity is limited inside the servo depending on the load. This is in agreement with the servo datasheet [10] as the first row of Table 3 also shows. However, the real maximum opening velocities against the load are lower then the values from datasheet. The approximate upper bounds are read from the measurements and summarized in Table 3. In case if the servo moves with the load the angular velocities seem to be almost unlimited, they have much larger absolute minimum values than maximums (except for the no-load case when the min/max angular velocities are the same, see Fig. 13). This case only approximate readings can be done again summarized in Table 3 motion the angular velocity is in a transient at all the time (see Fig. 16). The inner limitation of angular velocity by the actuator electronics made it impossible to identify simple transfer function or state space system models from the reference angle and torque to the output as in [18] for example because there is characteristic difference between opening and closing speeds. That's why nonlinearities were pulledout and identified separately from the transfer function models of the servo. Examining the measured data in detail shows that there can be an inner reference angular velocity value depending on the load and the angle tracking error of the servo. The servo tracks this inner angular velocity reference (see Fig. 15 for example where the maximum angular velocity is almost constant during opening of the airbrake) with a given dynamics. So firstly this inner limitation of the angular velocity reference and its tracking error related dynamics is studied and identified in the next part.

Transfer Function Identification Between Tracking Error and Angular Velocity Reference
At first, plotting together the angle tracking error and the angular velocity (for the 1800 to 1500 PWM (30deg to 60deg) 1Hz, 0Nm case in Fig. 17) shows that the angular velocity follows a square reference signal while the tracking error continuously changes. This means that the tracking error should be saturated to get a limited inner angular velocity reference from it. Note that considering Figs. 12 and 14 there is a delay, an offset and a scale error between the reference servo angle and the real one. As possible negative angular values are present in the  Studying several cases and considering 3-4 steps delay from the start of saturated tracking error decrease to the start of real angular velocity decrease has shown that the saturation limit of the tracking error should be about 5.8 • = 0.10123rad.
An almost square reference can be easily generated by a simple transfer function from normalized saturated tracking error to angular velocity reference in the form: G(s) = A T s+1 however, for a positive tracking error the velocity reference limit should be different than for a negative and this difference should depend on the load torque (T L ) of the servo. The positive-negative difference can advantageously be described by the following system structure: where A P (T L ), A N (T L ) are load torque (T L ) dependent scalar coefficients, T = 0.003s is the time constant determined by trial and error to give a realistic reference considering the system answer (see Fig. 18   The curves related to these values are also plotted in the Figs. 19 and 20 ('Refined model'). Of course the formulae have to be converted to give rad/s output in the form: After giving a method to generate inner angular velocity references from angle tracking error and load torque the tracking dynamics of this reference should be identified.

Servo Dynamics Identification Between Angular Velocity Reference and Angular Velocity
The servo dynamics from delayed inner angular velocity reference to measured angular velocity was identified using the Matlab System Identification Toolbox output error function considering zero delay and varying numerator and denominator polynomial degrees in discrete time. Finally, zero numerator and 3 denominator degrees (similarly to [8]) These discrete time poles are stable and oscillatory as the measured system output is also oscillatory. The parameters of the transfer function are shown in the Appendix. The slight offset error in the figure is the result of non-perfect offset removal from reference -output signal pair (see Fig. 18 where the estimated velocity reference is only approximately zero). Integrating the real angular velocity of the servo model gives servo deflection angle so this way the identification of the servo dynamics is ready. However, a load dependent angular position offset was detected in the test bench measurements (see Fig. 14) which should be also identified. This is done in the next part.

Load Dependent Servo Angle Offset Dynamics
Examining the measured data shows that application of the load moves the system steady-state servo angle with some value. This can be caused by a missing I-term in servo control electronics. Identification of the load torque (T L ) output angle (α) offset dynamics can be done considering load application in the fix 1500 PWM measurement cases. The delay between load and angle was estimated as 4 steps and then a transfer function (G T L (z)) with zero numerator and one denominator degree was identified with 96.98% fit (parameters published in Appendix). The load change and angle change are shown in Fig. 22 and the transient of the identified transfer function is shown in Fig. 23. After identifying the reference and system dynamics a Matlab Simulink simulation for the servo model itself was constructed and the identified model verified driven by test bench data with load application. The block scheme of the Simulink model is shown in Fig. 24.
The block delay represents the 3 steps time delay of the reference. SAT and NORM represent saturation of the error signal to ±0.10123rad and normalization with 1 0.10123 = 9.8786 respectively. G ref (z) represents the discrete time equivalent of the identified inner angular velocity reference model (1) G sys (z) is the identified system model from angular velocity reference to angular velocity and G T L (z) is the identified system model from load to angle including also the 4 steps time delay.

Servo Model Verification Based-on Measured Test Bench Data
After the identification and Simulink model construction the dynamic model of the servo motor was verified basedon measured data from the servo test bench. As a first   Fig. 27. The first two figures show that the final system model acceptably models the servo dynamics both without and with load. In the third figure it can be seen that the original lower angular velocity bounds gave very large overshoot in case of downward moving servo that's why the limits were increased (decreased in absolute value) and the expressions for A P (T L ) and A N (T L ) were modified.
As Fig. 24 shows the required inputs of the servo model are reference deflection and load torque. The reference deflection is given for all test bench measurement (see Fig. 26 for example) and the real load torque of the servo is saved (see for example Fig. 31) so both can be given as simulation model input.
At first, the identified model was checked considering a measurement from Test campaign I with full range deflection 0−120 • at 5Hz square wave excitation with 8Nm load application through some time. This is the most critical test case with the highest frequency and load. As shows this range with this frequency is too challenging for the servo controller and it can not follow the reference signal even with zero load. Comparison with data from the identified model shows that there is a drift in the test bench behavior which is not considered in identification.
Checking the model for test bench measurements with a tractable range (30 − 60 • ) again from Test campaign I in Fig. 29 shows that in normal mode there is neither drift in test bench output nor in identified model. Figure 29 also shows that the identified model well follows the behavior of the real servo considering the angular deflections. Figure 30 shows the measured and simulated angular velocities and Fig. 31 shows the reference and the measured torque of the test bench for completeness.
The simulated angular velocities are close to the real ones except for the transient periods but the overall performance is acceptable. So finally, the largest difference between model and measurements is the downward overshoot of the angle and the different behavior in the transients (when torque is applied and removed). However, in these test bench measurements the maximum load torque was applied immediately which is not realistic considering the opening of the airbrake, there the maximum load torque appears only on maximum deflection and it increases gradually as the To verify the servo model with more realistic conditions test bench measurements with gradually increasing (following the increase in opening angle) load torque from Test campaign III are considered. Figures 32, 33 and 34 show the results for the maximum load case when servo deflection to 97.5 • is followed by a load increase to 8Nm. The figures show that both the simulated angular velocity and angle follow well the measurements so the description of the real dynamics is satisfactory. The only difference is a large angular velocity glitch in the measurements when the  load is suddenly removed (see Fig. 34). Carefully examining Fig. 32 shows that the identified model well follows the small angle change in this case (at about 12.5s) and so the lack of the glitch does not cause any discrepancy and so the identified model's behavior is acceptable.
After the verification of the servo model the static characteristics of the opening mechanism and the airbrake aerodynamic effects are examined and described in the next subsection together with the possible dynamic effects of opening mechanism stiffness and airbrake inertia.

Determination of Static Characteristics
The most important static characteristic of the airbrake is the aerodynamics. Considering the dynamic pressure p dyn = ρ 2 V 2 known (here ρ is air density and V is airspeed) the aerodynamics can be described by the drag force c D and normal force c N coefficients depending on the opening angle φ of the airbrake. The normal force coefficients are published in [9] while the drag force coefficients are first published here in Table 5.
Regarding the drag the calculations were done considering a flat plate with a given opening angle based on [12]. Plotting the data together (see Fig. 35) it seems to be     Regarding the normal force coefficient its application point also changes with the opening angle. This is characterized by the a(φ) arm length in Table 5. As the final goal is to get an airbrake model as compact as possible its worth to unite the normal force coefficient and the arm length as The required servo torque to move the airbrake can be calculated considering the geometry of the opening There is a nonlinear relation between the angles β, γ , δ, φ summarized in tabular form in Table 6. Note that for small β angles δ can be negative as shown in Fig. 37 and for large β angles γ becomes negative (it is shown in the positive range in the figure).
The servo torque for a given airbrake opening angle φ, airspeed V and air density ρ can be determined as: In the second part of the equation all the geometric parameters are summarized in a virtual servo arm length s a (φ). The curves for c N (φ) are determined before so a curve for the s a (φ) = e The tabular data and the fitted 5th degree curve are shown in Fig. 38. The curve parameters are given in the Appendix.
After fitting polynomial models to the normal force coefficient and the virtual servo arm it is advisable to check Fig. 37 The opening mechanism of the airbrake (β servo arm angle, φ airbrake opening angle the precision of fitting. The largest servo torque will be required on sea level so servo torque values are calculated on sea level for airspeed from 10m/s to 60m/s and for all possible φ opening angles from Table 5. Calculations were done both from the mechanism and tabular data and from the fitted curves. Figure 39 shows that the approximate (dashed) curves are close to the real ones at every point while the absolute difference increases with airspeed. The relative errors of the fits (percentage) are shown in Fig. 40. These percentages are independent from the air density and airspeed. The figure shows that the maximum error is about 6% at φ = 10 • and then the errors are decreasing. The relative error can not be calculated for zero angle with zero torque that's why it is not plotted. Extrapolating the other values the relative error around zero angle can be about 8 − 10%.
After fitting curves to the aerodynamic coefficients and moment arm fast conversion formulae from servo angle α to airbrake angle φ and vice versa are required. The servo angle α is different form the servo arm angle β. It is defined to be zero at the minimum servo arm angle: Regarding the curve fits α − φ resulted as a 3rd degree polynomial going through zero while φ − α as a 4th degree polynomial again going through zero. The polynomials are summarized in the Appendix while the fit qualities can be checked in Fig. 41 (they are superior).
Finally, the stiffness characteristics of the airbrake mechanism should be identified to increase precision of deflection model considering the flexibility of the mechanism. Theoretical airbrake deflection angle from commanded actuator position and transmission of the rigid kinematic chain and measured deflection angle from the airbrake mock-up (see Fig. 8) were compared for a given set of measurement points. From the angle deviation, the system stiffness was derived (for details see [9]). The goal of this work was to extend the covered range by the stiffness formula as only the φ opening range 20 − 60 • was covered by the model in [9]. Finally, the stiffness data was converted from degrees φ angle input and output to radians and a curve fit was done on the data resulting in a third degree polynomial (k φ ) which is presented in the Appendix. The plot of the stiffness against the opening angle can be seen in  Fig. 42. The fitted formula can cover the whole angle range 0 − 60 • with reasonable non-negative values. As opening mechanism flexibility is present in the system its effects on the dynamic behavior should be examined. This is the topic of the next part.

Effects of Airbrake Inertia and Mechanism Stiffness
The inertial moment of the flap was estimated by modelling the flap as a rectangular prism with an external axis of rotation in distance c = 51.5mm, as illustrated in Fig. 43. The length l = 218mm is assumed corresponding to the edge length of the actual flap, the width b according to according to the actual sandwich core thickness of the structure. An isotropic density distribution is assumed. To compensate for the external axis of rotation, the parallelaxes theorem is used. The inertial moment can be estimated by using the following equation (see [7]): With the second term being the compensation from parallel axes theorem and using the measured mass m = 0.115kg of the actual flap. As the formula shows the width (b) is not required to estimate the moment of inertia with regard to the given axis. Due to the slight trapezoidal planform of the flap, the actual center of gravity is closer to  the axis of rotation. In consequence, the inertia is estimated a bit conservatively. Figure 44 shows the simple dynamical model of the servo-airbrake setup considering the stiffness k φ of the mechanism as a torsional spring. α is the servo arm angle (β) shifted to be zero at minimum deflection as α = β − 18.612 • while φ is the deflection angle of the airbrake itself. For the relation between β and φ through the airbrake deflection mechanism see Fig. 37. To simplify the model a one-to-one mapping between α and φ is assumed as any static gain between them would not modify dynamic behavior. The servo torque T equals the elastic force in the mechanism at every time as: The dynamic equation of the airbrake considering its mass moment of inertia I = 2.127 · 10 −3 kgm 2 the torque of the aerial forces on the airbrake T L and the fact that it is driven through the flexible mechanism results as: Note that here the airbrake load torque is transformed to the servo shaft (T L ) and the stiffness of the mechanism is considered as a simple torsional spring between. Reorganizing the above expression and making Laplace transform of each side results in two transfer functions, one  from servo angle α to airbrake deflection φ: G α (s) and one from aerial load T L to airbrake deflection φ: G L (s) So finally, one gets a two input (α servo angle and T L torque of aerial forces) one output (φ airbrake angle) system. Note that the k φ term depends on the actual angle of the airbrake (see Fig. 42) and so these transfer functions should be examined for different airbrake angle values. Bode plots were plotted for angle values 0 : 5 • : 60 • covering the whole range and shown in Figs. 45 and 46.  They show that both term can be well approximated with a constant gain on a wide frequency range. They give undamped oscillations only on high frequencies (well above 70 rad/s=11Hz) on which the airbrake will not be operated as its bandwidth is about 1.4H z shown by Fig. 11. Another fact to be considered is that the mechanism should also have damping which is unknown but surely present and so could damp the oscillations to a favorable level even on high frequencies. If flight testing does not show oscillations even with high frequency excitation then a static model including only the effect of stiffness can be used.
Considering low frequency approximations of the transfer functions one gets: Taking the inverse Laplace transform and reorganizing the terms gives: Where Δφ is the angular deformation and can be obtained as Δφ = α − φ = T L k φ . In the airbrake simulation the servo deflection α and the momentarily load T L on the servo are known and so k φ and φ can be obtained as shown in the structure in Fig. 47.

Full Airbrake Simulation Model and Test Results
After identifying the model components the whole system with servo model, mechanism and aerodynamic effects was constructed and simulated in Matlab Simulink. The block scheme of the whole simulation structure is shown in Fig. 47. Here, Servo is the servo simulation model from Fig. 24, p dyn calculates the dynamic pressure from ρ air density and V airspeed, α(φ) and φ(α) are the polynomials transforming airbrake angle to servo angle and vice versa. k φ is the angle dependent stiffness coefficient, the load torque T L is divided by it to get the angle change (Δφ) of the airbrake as derived at the end of Section 4.2. To be a bit more realistic the whole angle change is subtracted from the airbrake position used to calculate air drag and normal force, but only half of it is subtracted to calculate the virtual servo arm s a (φ) assuming non-uniform deformation inside the mechanism. SAT represents the saturation of the airbrake angle output between 0rad and 1.0472rad (0 − 60 • ). The blocks with three × represent the multiplication of all inputs. D is the air drag force as the useful output of the system contributing to the deceleration of the aircraft.
After the construction the model was tested with step airbrake angle reference changes, airspeed and air density (flight altitude) sweeps and even with chirp signal input. The overall system was stable and realistic in all cases so it was built into the nonlinear simulation model of the FLEXOP aircraft. These test cases are not presented in detail here to hold the length of the paper acceptable, rather the final SIL simulation results of the full aircraft nonlinear model are presented.

Airbrake Model Verification in SIL Mission Simulation
The full airbrake model was integrated into the nonlinear simulation model of the FLEXOP aircraft by DLR (Deutsches Zentrum für Luft-und Raumfahrt e.V.) and then the baseline controller of the aircraft [14] was completed with a simple logic which fully opens and then closes the airbrake to help deceleration as the aircraft approaches the end of the test leg. The baseline controller is a cascade gain-scheduled PID control with several inner and outer loops which is detailed in [14]. Giving more details is not the topic of this article. It is possible that in the future the airbrake operation logic will be replaced by a more sophisticated controller but if flight test results with the logic will be satisfactory then this won't be required. At first, the standalone airbrake model deflection and load values were compared to the SIL results to check if the integration of the model into the nonlinear simulation is flawless. The standalone model is driven by deflection commands and dynamic pressure saved from the SIL simulation. Airbrake load, servo and airbrake deflections are compared in Figs. 48 and 49. In Fig. 48 'SIL' means the result from the nonlinear SIL simulation of the FLEXOP aircraft while 'Standalone' means other results from the standalone simulation. 'servo' means the deflection of the servo arm while 'airbrake' means deflection of the airbrake. The load curves are almost exactly the same while the deflections show different transient responses but otherwise they perfectly match. The SIL airbrake and servo deflections have some over-and undershoot in the transient while the standalone ones do not. This is caused by the different solver settings in the standalone (fixed step, discrete) and SIL (fixed step, ode3) Matlab Simulink simulations as the ode3 solver is iterative. In conclusion the airbrake model integration is flawless the nonlinear model can be used for controller tuning, SIL and hardware-in-the-loop (HIL) tests.
The next step after the verification of airbrake integration into the model was to test its performance accelerating the demonstrator to different flight speeds and then decelerating it with or without airbrake. The nominal cruise speed of the demonstrator is 38m/s and the airbrake is started to be closed as the speed decreases to 39m/s during The airspeed characteristics are plotted in Figs. 50, 51 and 52. The without airbrake airspeed curves are shifted to start deceleration at the same time as the curves with airbrake application. That's why they seem to react to the initial reference change too late. The figures show that with the application of the airbrake the airspeed can reach the 38 m/s cruise value earlier with smaller or no overshoot.  Table 7 summarizes the settling times of SIL simulation regarding ±2% tolerance around the 38 m/s cruise speed (37.24 − 38.76m/s, the usually applied ±5% range was too large with almost 40 m/s upper bound). The time gain means the time with which the deceleration with airbrake is shorter and so the flutter test can be longer. The distance gain is the time gain multiplied by 38 m/s which is the minimum gained distance on the test leg. The table shows that the flutter test can be from 1.5s to as large as almost 4s longer with airbrake application which can largely help to better evaluate the behavior. After SIL testing the final verification of the airbrake model is done in the next section considering flight test results.

Flight Test Results
Real flight test results were obtained in two manual flights in August 2019 at München Oberpfaffenhofen airport. The aircraft was manually controlled into almost straight and level flight and then applying constant throttle (see Fig. 60) and trying to hold the flight direction (see Fig. 53) the airbrake was manually fully opened and closed giving a staircase reference until maximum deflection (see e. g. Fig. 54) and after some time back. Straight and level flight was targeted to give an opportunity to realistically evaluate the effect of the airbrake in decreasing the airspeed of the aircraft. The tests are approximately repeatable as the similar altitude and airspeed ranges show in Table 8. Unfortunately manual control can not perfectly hold altitude (see Fig. 59) so finally the system specific energy was used for evaluation. In the future in frame of the next FLiPASED (see [1]) project airbrake tests are planned with autopilot control during 2020 which will better be repeatable. However, precise evaluation of airbrake effectiveness is not strictly required to evaluate the airbrake dynamic model itself that's why model identification results can be published now.
The ranges of barometric flight altitude (h) and indicated airspeed (V ) are shown in Table 8 and Fig. 59 shows their relative changes starting from the first value at activation of airbrake. The figure well shows that these ranges do not cover monotonically increasing or decreasing values, both parameters can decrease and increase even multiple times during one airbrake application test section. That's why later airbrake effectiveness is evaluated based-on the specific energy of the system because there is no definite change of the airspeed and/or altitude. Figure 53 shows the flight trajectories which are close to straight during airbrake application. In autopilot controlled flight the direction hold of the aircraft will be much better than here in manual flight. The airbrake was fully opened and then closed at constant throttle value (see Fig. 60). As the on-board system saves the RC (radio controller) commands and measures the deflections of the actuators its possible to plot the command and answer together and simulate the airbrake deflection  based-on the command and the also logged barometric altitude and airpseed. Basically the servo deflections are commanded and measured but the airbrake opening angle is considered in the calibration characteristics. Examining the in-flight opening characteristic (see Fig. 54) shows that the opening position seems to occur before the command which is physically impossible. Taking a closer look shows that the start and end of opening is after the command starts and ends so the overall too early opening should be virtual only because of the large quantization error of saved command. This is underlined by the closing characteristic (see Fig. 55) where the airbrake deflection is after the command. Note that the manual command is highly quantized while the airbrake angle output is a smooth signal so the command on the servo should be smoother than the logged signal. This means that the airbrake simulation results will be quantized which should be neglected in the model comparison. Also in the opening cases the lag of the simulated curve is realistic relative to the command while seems to be too late relative to the in-flight deflection. Figure 54 also shows that the rise time is about 0.3−0.35s which is a bit higher than the ≈ 0.25s value from Table 2 related to the about 3Nm load resulting (consider Fig. 39) from the about 40m/s airspeed in flight 2 during airbrake opening. The airbrake model was simulated considering the actual barometric flight altitude and indicated airspeed values and the manual deflection commands of the airbrake. Note that the maximum physical deflection of the airbrakes on the aircraft is about 45 • in contrast to the maximum 60 • test bench deflection. The steady state values of simulated and in-flight airbrake angles are very close to each other (max 0.6 • difference) so this does not require model refinement.
The comparison of simulated closing characteristic (as the opening can not be taken as a precise reference) with the in-flight one shows that the simulation with the three time steps delay in the servo dynamics lags more than the flight data. This means that the delay could be decreased finally to 1 step giving results better following the flight data (see Fig. 56). The figure also shows that the slope of the simulated signal is a bit larger than the in-flight one. This can be possibly solved by further limiting the maximum closing speed of the servo under load. Examining the closing speed limits and the load resulting from inflight dynamic pressure shows that the load is about 2Nm and so the current maximum closing speed is −372 • /s in the simulation model. The maximum limitation of the speed can be −340 • /s but this is very close to the current model (91.3%) and does not improve significantly the performance so finally the closing speed characteristic remained unchanged. As the airbrake closes only if the airspeed is decreased to a low value the load during the closing will always be small so the load dependent closing speed will be around the above mentioned values.
Considering the slope of the opening in Fig. 54 in the first case the simulated slope is a bit larger while in the second case its almost the same as the in-flight one so finally the opening angular velocity limits are unchanged. Figure 57 shows that both the in-flight and the simulated data has a change in airbrake position when the manual deflection command is constant. This should be caused by the decrease in airbrake load and is well followed by the simulation model. This verifies the approach to identify also the load to servo angle transfer function G T L (z) and excludes the possibility of having this effect only due to test bench measurement errors. Figure 58 shows the commanded and real deflections of the left airbrake in the 2nd flight. Similar half doublet deflections were applied in all test cases and the airspeed, altitude and specific energy results are obtained with these  Figure 59 shows the change in airspeed and altitude during airbrake application. All initial values were shifted to zero to show the decrease or increase of the values relative to the initial. The airspeed decreases with about 2 m/s in the first flight and more then 5 m/s in the second while the time horizons are about the same. This is because in the first flight the initial airspeeds are 34-36 m/s while in the second flight they are about 43 m/s (see Table 8). As the drag force of the airbrake scales quadratically with airspeed this is not surprising as higher drag means higher deceleration. Regarding the altitude it increases in one of the cases while decreases in all other three. As the altitude can increase / decrease during the maneuvers its hard to estimate if the decrease in airspeed is caused by altitude change or by the airbrake.
To decide about this the overall energy state of the aircraft should be considered. With constant throttle and small changes in altitude and airspeed the thrust of the jet engine can be assumed to be almost constant and so starting from a trimmed condition the energy content should not change significantly. However, plotting the specific energy u = gh + 1 2 V 2 (see Fig. 60, h is barometric altitude and g is the gravitational constant) shows that the energy content of the aircraft decreases so the airbrake removes energy from the system. As the mass of the aircraft slightly changes due to fuel consumption the specific energy is considered instead of the whole potential and kinetic energies.
As a conclusion it can be stated that the flight test data verifies the airbrake simulation model having similar results in simulation then in flight so this model can be applied for more sophisticated airbrake control design and for the HIL verification of any autopilot applying airbrake.

Conclusion
This paper presents the modeling and system identification of the airbrake of an unmanned experimental aircraft. This airbrake consists of a servo motor, a nonlinear mechanism and the aibrake flap itself. After briefly introducing the demonstrator aircraft and summarizing the design considerations the system identification concept is discussed. Then the applied test benches are briefly described which are a full scale mock-up and a servo test bench with load motor.
From this point the work focuses on the identification of the servo dynamics based on test bench measurements first pulling out the significant nonlinearities from the system. These include the characterization of the load dependent opening and closing angular velocities and the estimation of the saturation level of servo angular velocity reference input. The remaining linear part of the servo dynamics was identified as a transfer function plus delay term. As the load causes a steady-state servo angle change this load dependent angle change was also modeled as a transfer function. This effect was verified by flight test results so this transfer function is a required part of the airbrake model. After constructing the servo model its verification based-on test bench measurements was done. The next step was the description and identification of static characteristics. This includes the aerodynamics as drag and normal force coefficients. The former is modeled with a line fit against the airbrake deflection angle, the latter is normalized and modeled as a combination of a parabola and a line again against the airbrake angle. Then a virtual servo arm is defined and its characteristics are determined considering the nonlinearities in the airbrake opening mechanism. This, together with the normalized normal force coefficient gives the relation between the dynamic pressure, airbrake opening angle and servo load torque. After generating these the airbrake opening angle -airspeed -servo load torque characteristics are plotted for the whole operational range of the aircraft. The overall precision of this modeling resulted between 6 to 1 percent depending on the opening angle (compared to the calculated tabular data). These characteristics made it possible to select the required measurements for system identification which represent operational cases from a set of test bench measurements (opening angle range and load range sweep). Finally, the servo angle to airbrake angle and the reverse conversion characteristics were calculated by fitting polynomials and the angle dependent stiffness of the opening mechanism was also characterized. Before construction of the whole airbrake model the effect of opening mechanism stiffness and airbrake inertia on the dynamics was examined. It has turned out that this can cause oscillations only on very high frequencies which are outside the planned operating range and estimated bandwidth of the airbrake.
The whole airbrake simulation model was first checked in Monte-Carlo simulation with several angle references, airspeeds and flight altitudes and showed realistic results. Then, the airbrake model was integrated into the nonlinear simulation model of the FLEXOP demonstrator aircraft and its effect on the deceleration characteristics was examined considering software-in-the-loop trajectory tracking with a baseline controller. It has turned out that airbrake application can gain 1.5 to 4 seconds additional time to test the flutter controller of the aircraft.
Finally, in-flight manual test of the airbrake was done in August 2019. Comparing the simulation model with the inflight results has shown that the delay of the servo model can be decreased (from 3 time steps to 1) but otherwise the model well covers the opening, the closing and even the load change characteristics of the real airbrake. As in manual flight there is no precise altitude hold but the throttle and so the thrust was approximately constant the specific energy of the aircraft was examined during the application of the airbrake. This shows that the airbrake decreased overall system energy in all of the cases so its application in deceleration is advantageous.