Estimation of Aircraft-Dependent Bumpiness Severity in Turbulent Flight

Featured Application: This work proposes a severity estimation method of aircraft bumpiness in turbulence, which is fundamental to turbulent ﬂight safety of civil aviation aircraft. Abstract: Atmospheric turbulence threatens ﬂight safety of civil aviation aircraft by inducing aircraft bumpiness. A severity estimation method of aircraft bumpiness in turbulent ﬂight is explored according to in-situ Eddy Dissipation Rate (EDR) indicator. With the turbulence intensity derived from EDR value, a time series of longitudinal and vertical turbulence was generated according to von Karman turbulence model. In order to obtain the vertical acceleration response of aircraft, the continuous change of aerodynamic force on the assembly of wing and horizontal tail was computed by Unsteady Vortex Lattice Method (UVLM). The computing accuracy was improved by using semi-circle division and assigning the vortex rings on the mean camber surface. Furthermore, the adverse effects of control surface deﬂections on bumpiness severity estimation can be effectively removed by separating turbulence-induced and aircraft maneuvers-induced aerodynamic force change. After that, the variance of vertical acceleration, as the severity indicator of aircraft bumpiness, was obtained by Welch spectrum estimation. With the reﬁned grid level, the pitching moment change due to control surface deﬂections can be solved accurately by UVLM. The instantaneous acceleration change obtained by UVLM approximates recorded acceleration data with better accuracy than linear transfer function model. A further test with a set of ﬂight data on the same airway shows that compared with in-situ EDR indicator, the proposed method gives an aircraft-dependent estimation of bumpiness severity, which can not only be used to estimate in-situ bumpiness but also be applied to forecast the bumpiness severity of other different aircrafts. This paper is the ﬁrst attempt to apply the UVLM on the estimation of aircraft bumpiness severity in turbulence. With the proposed method, the bumpiness severity of the next aircraft ﬂying through the same turbulence ﬁeld can be predicted. In the future, the vortex ring model and related UVLM will be reﬁned to further improve the computing accuracy of vertical acceleration. In addition, by extending and applying the method in Chinese airlines, the practicability of the method can be improved by closed-loop veriﬁcation with ﬂight data.


Introduction
Inducing civil aviation aircraft to unexpected bumpiness, atmospheric turbulence is by far the leading cause of injuries [1]. The commonly used Eddy Dissipation Rate (EDR) value indicates turbulence intensity, which is unrelated to aircraft. However, to avoid or mitigate the adverse effects of turbulence, the bumpiness severity of specific aircraft needs to be determined according to the objective EDR indicator [2].
Since the rapid change of vertical acceleration induces aircraft bumpiness [3], pilots can percept aircraft bumpiness by the change of vertical acceleration. Pilot weather report (PIREP) is a traditional indicator of aircraft bumpiness. However, the indicated bumpiness severity is not the same experienced by the other aircraft because of different aircraft types and different flight states [4]. Furthermore, as a subjectively reported indicator, the reported value can easily be changed by different pilots. rapid-changing turbulent wind is different from low-frequency wind shear, the steady VLM lacks enough accuracy in computing aerodynamic force in turbulence. Instead, if the turbulence-effected aerodynamic performance is obtained by UVLM, a more accurate vertical acceleration response can be acquired, which forms a better basis to obtain an aircraft-dependent bumpiness severity estimation.
In this paper, a method of estimating aircraft-dependent bumpiness severity in turbulent flight is explored with in-situ EDR value provided. The intensity of turbulence is first obtained according to the von Karman and Kolmogorov turbulence theory. With the generated time series of turbulence as the excitation, the continuous change of vertical acceleration in turbulent flight is computed based on UVLM. After that, an aircraft-dependent bumpiness severity indicator is computed by Welch spectrum estimation. Compared with the recorded vertical acceleration in flight data, the acceleration response in turbulent flight is deeply analyzed. Furthermore, according to in-situ EDR value, the bumpiness severity of different aircraft types on the same airway is estimated.

Turbulence Generation Based on In-Situ EDR Value
With the isotropy and large Reynolds number assumption, there exists an inertial subrange of wave numbers in the turbulence field. The turbulence energy is neither generated nor dissipated. Instead, it is just transferred inertially from larger to smaller eddies. The energy spectrum of turbulence E(k) only depends on ε, the Eddy Dissipation Rate (EDR), by where A = 1.6 and k is the spatial frequency of turbulence [3]. The von Karman and Dryden turbulence model allows researchers to describe the random characteristics of isotropic and small-scale turbulence in high-altitude. The energy spectrum of von Karman model has a roll-off rate of -5/3 in high-frequency section, while the roll-rate of Dryden model is -2 because Dryden model is an approximation of von Karman model. Furthermore, faster-moving aircraft is also responsive to larger scales, which are typically outside the inertial sub-range. The von Karman turbulence model, representing both the inertial sub-range and the larger scales beyond it, has been widely used in both aerodynamics and meteorological field [3]. For large frequencies in the inertial subrange, the energy spectrum of the von Karman model is consistent with that of the Kolmogorov model [30]. The von Karman energy spectrum is described as where Γ(·) is the gamma function, σ 2 i refers to the intensity of three turbulence components. By equating Equations (1) and (2), the relationship between σ 2 i and ε is Based on the assumption of isotropy, the intensities of three turbulence components are identical, which can be described as σ Wx = σ Wy = σ Wz . The length scale in three directions is nominally set by L 0 = 669 m [31]. The EDR indicator is generally defined as ε 1 3 , the cubic root of eddy dissipation rate. According to Equation (3), the turbulence intensity is obtained from in-situ EDR indicator.
With the provided intensity and length scale, the longitudinal and vertical temporal spectra of von Karman model are described by where ω is the temporal frequency of turbulence. L Wx and L Wz stand for the longitudinal and vertical length scale of turbulence, respectively. In order to obtain the vertical acceleration response of aircraft, a sample of turbulence components needs to be generated beforehand by feeding unit-intensity white noise into forming filters. By decomposing Equation (4), the transfer functions of forming filter are described by By reducing to the one-order for rational approximation and further discretizing with the first-order backward differential method, the forming filters in Equation (5) are transformed into the following differential equations as [32] where T s stands for the sampling period. In this way, the longitudinal and vertical turbulence components are generated according to in-situ EDR value. After that, the time series of turbulence components is further used as the excitation of acceleration response model. It should be noted that a full turbulence model not only contains the turbulence components W x , W y and W z , but also contains three angular velocity components, which are q W = ∂W z ∂x , p W = − ∂W z ∂y and r W = − ∂W y ∂x [33]. Compared to three turbulence components, the angular velocity components have small effects on aircraft dynamics. Thus, the derivations of angular velocity components are not considered in this paper.

Vertical Acceleration Response by UVLM
In this paper, aircraft vertical acceleration response excited by turbulence is computed via non-planar Unsteady Vortex Lattice Method (UVLM). The Laplace's equation governs the principle of irrotational, inviscid, and incompressible air flow [15]. Based on linear Laplace's equation, the flow around complex geometries can be solved by superposing elementary solutions. To improve the computing accuracy, the vortex rings, each composed of four vortex filaments, are placed on the mean camber surface of the lifting body. Compared with the horseshoe vortices assigned on the flat surface of the lifting body, the accuracy can be improved by assigning vortex rings.
In high-altitude turbulent flight, the longitudinal and vertical turbulence components make the longitudinal aerodynamic force change rapidly, further leading to the fluctuation of vertical acceleration. In order to obtain the longitudinal aerodynamic performance with better computing efficiency, target aircraft is simplified to a lifting body of wing-horizontal tail assembly. As shown in Figure 1, vortex rings are assigned on the whole reference area of the assembly and distributed on the mean camber surface. The wing-tail assembly is further divided by dozens of vortex rings. Figure 2 shows the detailed grid arrangement and geometric parameters on the right-hand side. Semi-circle division method is used to design the neighboring lattices at the wingtip, wing root, and leading edge because the circulation distribution changes rapidly there [34]. In addition, the boundaries of each vortex ring are close to the geometric edge of control surfaces since elevators and spoilers can be deflected by pilot's manual operation or Automatic Flight Control System (AFCS). The arrangement of vortex ring is illustrated in Figure 3. According to the vortex lattice method, the leading edge of the vortex ring is placed at the 1/4 chord of the panel, while the Neumann boundary condition of no penetration is enforced by locating the collocation point at the 3/4 chord. In this way, the lift curve slope of the lifting body conforms to the thin airfoil theory [35]. As shown in Figure 3, the induced velocity of the vortex ring is described as the vector sum of the induced velocity superposed by four vortex filaments, where k is the coefficient vector defined according to the Biot-Savart theorem [15]: where Γ is the strength of the vortex filament, and r 0 is the unit vector of the vortex filament.  A time-marching computing process is used to perform the unsteady simulation. At each time step t k , the newly generated vortex rings at the trailing edge of the lifting body are released into the wake. As shown in Figure 4, the wing-tail assembly releases force-free wakes with no aerodynamic loads. The circulation Γ of the newly generated wake panel stays unchanged for the remainder of the simulation according to Helmholtz's theorem [15]. Since the scale of the high-altitude turbulence is large enough compared to the aircraft, it is feasible to analyze the boundary condition and further compute the aerodynamic force by UVLM according to the change of turbulence at each time step. At each time step t k , the local velocity induced by unsteady motion at an arbitrary point can be divided by the far-field free flow velocity, the effects of aircraft angular motion, and instantaneous turbulence. As a result, the Neumann boundary condition at any collocation points is described as An in-depth analysis of the boundary condition is shown in Figure 5, which shows the decomposition of boundary conditions along the wing chord. From Figure 5, the boundary condition at the collocation point is expressed as where V ∞ represents the far-field free flow, δ stands for the incretion of angle of attack due to control surface deflection. Besides, α i = tan −1 (− dz dx ) shows the camber effects.
[u, w] T represents the induced velocity decomposition along each axis, and [x, y, z] T stands for the coordinate of arbitrary point. A linear algebraic equation is derived by further expanding Equation (9). The circulation of each vortex ring is obtained by has major effects on local velocity. The local velocity of arbitrary point at any time step t k is described by The aerodynamic force is further computed by Kutta-Joukowski lift theorem based on the circulation of each vortex ring. Taking the right half wing as an example, any vortex and its neighboring vortices are superposed on the lattice panel to obtain total circulation distribution. As shown in Figure 6, the aerodynamic force at the collocation point is: where ρ is air density, V b,i is the total induced velocity of airfoil-attached vortex at the collocation point. In order to differentiate the aerodynamic force induced by aircraft maneuvers and external turbulence, Equation (12) is solved twice, with or without turbulence effects, respectively, which are recorded as F i and F 0 i . At each time step t k , total aerodynamic force F = [F x , F z ] T is the force sum of all vortex rings on the wing and horizontal tail. Besides, the pitching moment around the gravity center is obtained by where x cg , z cg represents the location of gravity center. According to the flight dynamics model, the differential equations of vertical acceleration and pitching moment induced by turbulence are where θ(t k ) is pitching angle at t k , m is the mass of aircraft, and I y stands for the inertia moment around pitching axis. In order to separate aircraft maneuvers effects from the acceleration response, Equation (14) also needs to be computed twice. Firstly, total aerodynamic force F(t k ) is used in Equation (14) to get a total acceleration response a z (t k ), induced by aircraft maneuvers and external turbulence. After that, the aerodynamic force only with aircraft maneuvers, F 0 (t k ) acts as the input to get the acceleration a 0 z (t k ). As a result, the incremental acceleration value only induced by turbulence is obtained by a z (t k ) − a 0 z (t k ).
There inevitably exists a gap of longitudinal aerodynamic performance between the whole aircraft and the wing-tail assembly. However, this paper makes a compromise between the computation efficiency and algorithm accuracy. Compared with the fuselage, aircraft wing and horizontal tail play a major role in pitching and plunge motion, while the integration of fuselage will greatly increase the algorithm complexity. Furthermore, only the acceleration increment caused by turbulence needs to be considered for aircraft bumpiness estimation. Thus, the error accumulation during the computation process can be reduced.
To summarize, with the time series of turbulence components as the input, the time series of vertical acceleration is obtained by UVLM. The next step is to get the estimation of bumpiness severity with the computed acceleration data.

Severity Estimation of Bumpiness
Compared to the RMS-g by directly obtaining the root mean square of acceleration data, this paper develops a more theoretically accurate method to estimate the variance of vertical acceleration, which functions as the indicator of aircraft-dependent bumpiness severity. In turbulent flight, the intensity and spatial frequency of turbulence are measured by the flying-through aircraft, inducing the change of vertical acceleration. The Welch spectrum estimation is used to obtain the power density spectrum of measured acceleration data. Furthermore, as the aircraft bumpiness indicator, the variance of vertical acceleration is obtained by Parseval theorem.
Firstly, a Tukey-Hanning window is used to decrease the spectrum leakage due to finite-length time series of vertical acceleration. The m-point Tukey-Hanning window is formulated as where M = f loor(0.1m − 0.2), which corresponds to a taper factor of 0.2 in traditional formulation of Tukey-Hanning window. The power normalized Tukey-Hanning window is Accordingly, the windowed vertical acceleration data series is Furthermore, by data segmenting and overlap setting, Welch spectrum estimation is used to obtain the power spectrum of a w z (t k ) [36]. Given that aircraft maneuvers and aero-elastic vibration would also induce the change of vertical acceleration, a band-pass filter is designed in EDR estimation to eliminate these adverse effects. Similarly, in the estimation of bumpiness severity, the high-end cutoff frequency f h lies below the aeroelastic response mode of target aircraft, while the lower cutoff frequency f l is set higher than the phugoid frequency of target aircraft. According to Parseval theorem, the variance of vertical acceleration in time domain equals the power spectrum in the band-pass of the frequency domainσ 2 To summarize, an integrated computation flow for the severity estimation of aircraft bumpiness is shown in Figure 7. The grid on the wing-tail assembly of target aircraft is generated beforehand. According to real-time in-situ EDR indicator, a time series of longitudinal and vertical turbulence components is generated by Equation (6). In the algorithm loop, the generated turbulence and flight parameters are used to compute local velocity by Equation (11) firstly. After that, the aerodynamic force in turbulence is obtained by UVLM. Vertical acceleration and related flight parameters necessary for computing local velocity are obtained by solving the flight dynamics model. For each EDR value, the time series of vertical acceleration data is generated based on the input of turbulence components. If each time step of turbulence components is processed, the estimated aircraft bumpiness indicatorσ 2 a z is obtained by Welch spectrum estimation. Section 2 provides a severity estimation method of aircraft bumpiness based on in-situ EDR indicator. The estimation method is dependent on the geometry of the wing-tail assembly and several flight parameters of specific aircraft. Compared to the aircraft-dependent RMS-g and DEVG indicators, which are computed according to recorded in-situ acceleration data, it provides an effective way of estimating the bumpiness severity in advance. Furthermore, since the adverse effects of aircraft maneuvers are eliminated in computing turbulence-induced acceleration response, the computing accuracy of acceleration response can be improved. On the contrary, the recorded in-situ acceleration data are blended with aircraft maneuvers, which leads to an in-accurate bumpiness severity estimation by RMS-g and DEVG.

Grid Refinement
The performance of UVLM was tested as a foundation for the acceleration response analysis. The non-planar UVLM requires a reasonable grid division (or grid resolution) to achieve aerodynamic performance that is not related to lattice number. If the grid resolution is low, the vortex rings may not fit the mean camber surface well, thus leading to low accuracy. On the other hand, too high grid resolution will greatly increase the computation cost. Too small vortex lattices would also make the collocation points too close to the vortex filament as shown in Figure 3, which also reduces the computation accuracy. Taking the Boeing 737-800 aircraft as an example in this paper, the grid resolution of UVLM was first studied to improve the accuracy on the premise of computing efficiency assurance.
The geometric and aerodynamic parameters of B737-800 aircraft are listed in references [37,38]. According to the geometry and airfoil parameters of the wing and horizontal tail, a non-planar vortex lattice model was built, in which different surface grids of increasing mesh density were generated as shown in Table 1. The changes of the lift, drag, and pitching moment coefficients were presented as a function of the grid refinement level. The experiments show that the results of the eighth refinement level are within 1% of the values obtained with the most refined grid. It means that a satisfactory accuracy was obtained by using a grid with 40 columns, 20 rows on the wing and 12 columns, 6 rows on the horizontal tail.

Verification of Pitching Moment
There inevitably exists a gap between the aerodynamic performance of the whole aircraft and that of the wing-tail assembly. Fortunately, only the acceleration increment induced by turbulence is considered for aircraft bumpiness estimation. The pitching moment change, which has major effects on instantaneous acceleration change as shown in Equation (14), is mainly caused by the deflection of elevator or horizontal stabilizer. As a result, the pitching moment change based on wing-tail assembly should approximate the authentic aerodynamic data. The longitudinal pitching moment changes of B737-800 aircraft caused by elevator and stabilizer deflections are provided in [37,38]. Figure 8 shows the results between UVLM results and the modeling data in three flight conditions, where c mδstab represents the pitching moment change due to stabilizer deflection, and c mδe stands for the moment change caused by elevator deflection. Taking Figure 8a as an example, it shows the comparison of c mδstab between UVLM and modeling data in mach number M = 0.2, α = 3 • at sea level (H = 0). It can be found that different flight conditions would change the dynamic response of the aircraft. The above results are obtained by steady vortex lattice method because the test conditions remain unchanged. Several improvements on traditional VLM make the results approximate modeling data. Firstly, vortex rings are assigned on the mean camber surface rather than using the flat surface assumption. Secondly, the semi-circle division method is adopted to refine the lattices neighboring to wingtip, wing root, and leading edge. Thirdly, the divided lattices are as close as possible to the structural boundary of the control surface. As a result, the non-planar UVLM based on wing-tail assembly has satisfactory computing accuracy, which is beneficial to further analysis of vertical acceleration response. In addition, the accurate computation of pitching moment coefficient is also significant to separate aircraft maneuvers effects on bumpiness estimation.

Instantaneous Acceleration Response
Since accurate vertical acceleration response is significant to severity estimation of aircraft bumpiness, the acceleration response computed by UVLM was further compared to measured acceleration data provided by Quick Access Recorder (QAR) flight data. The airborne flight data acquisition system provides hundreds of flight data in the form of time series, including flight dynamics parameters and control surface deflections. Supported by China Academy of Civil Aviation Science and Technology, QAR flight data of 880 flight segments on the same scheduled airway from Beijing Capital International Airport (ZBAA) to Guangzhou Baiyun International Airport (ZGGG) were collected. According to the QAR data, eight Wind Fields (WF) were selected and used to explore the acceleration response, covering the light, moderate, and severe turbulence. The maximum and minimum vertical acceleration data of the eight Wind Fields (WF1~WF8) are listed in Table 2. The acceleration response by UVLM is compared with the linear transfer function model and recorded acceleration data. It is worth noting that the acceleration data is recorded at the gravity center of aircraft.
A linear transfer function model of B737-800 aircraft was built according to the modeling method of Model-4 in ref [7]. To better understand the acceleration response mechanism in turbulence, four linear models of increasing gust excitation complexity were developed in ref [7]. Considering the gust pitching rate effects in pitching and plunge dynamics and further incorporating a quasi-steady aerodynamic model with finite lag value, the Model-4 is able to track instantaneous acceleration change with the best accuracy while the accuracy of peak response is slightly weak.
The longitudinal and lateral wind components are first derived by [39] where V Gx and V Gz are the ground speed components recorded in QAR flight data. According to the QAR flight data of WF5, the time series of W x and W z are computed and shown in Figure 9a. With the wind components as the input, both the linear model and UVLM give the response of vertical acceleration. Figure 9b shows the recorded vertical acceleration data (QAR Record), the computed vertical acceleration from the linear transfer function model (Linear Model) and UVLM, respectively. Before 90 s in WF5, the data show the aircraft was at a trimmed level flight condition. From 91.5 s, the aircraft began to experience initial turbulence, approaching the minimum acceleration near 96 s, and then reaching the maximum positive load at 96.8 s. After that, the acceleration showed a continuously decrease. In this short period, the linear model was momentarily out-of-phase once encountering sudden turbulence, while the tracking accuracy was still satisfactory. However, observations from Figure 9 suggest that both linear model and UVLM are able to follow the acceleration response in moderate turbulence, while the tracking accuracy of UVLM is better than that of the linear model. Figure 10 shows the acceleration response of the two methods in WF8. This flight data sample shows a severe aircraft bumping event, in which the vertical acceleration varied from +2.034 g to +0.413 g in 1 s. The response of the two methods tended to track the overall behavior of the recorded data. However, as far as the Linear Model is concerned, the sacrifice of peak acceleration response accuracy led to the serious deviation from recorded acceleration data. In contrast to the linear model, UVLM performed much more sensitively to severe turbulence. The tracking MSE was reduced from 0.0197 of Linear Model to 0.0081 of UVLM. Quantitative comparison leads to the conclusion that Linear Model was less sensitive and caused more error when encountering severe turbulence. In severe turbulence, the recorded acceleration series can be better recovered by UVLM than by the linear model.  In order to compare the acceleration response more comprehensively, the MSE and absolute peak acceleration in WF1~WF8 are computed and shown graphically in Figure 11. The Linear Model did a good prediction in light turbulence as UVLM. However, the tracking error increased in moderate and severe turbulence. The increase of turbulence intensity leads to bigger MSE in linear model. On the contrary, UVLM is able to track the acceleration response with better accuracy than the linear model. It is not unexpected that the accuracy of acceleration response by UVLM is better than by Linear Model. On the one hand, severe turbulence could easily change the governing aerodynamic coefficients of Linear Model, leading to a severe deviation from trim condition. Differing from Linear Model, UVLM is a kind of nonlinear algorithm in which the highorder response feature can be computed, rather than linear approximation. On the other hand, the recorded acceleration response is excited by both turbulence and control surface deflections. The recorded control deflections are input into the UVLM, changing the boundary condition of local vortex ring, while the control deflections are not considered in linear transfer function model.  Figure 12b,c, respectively. It can be found that both the elevator δ e and high-speed flight spoiler δ s HS experienced big deflections. The deflection of elevators had a greater impact on aerodynamics force than that of spoilers. Figure 12d shows the computation results of aerodynamic force with or without control deflections by using UVLM. By differentiating the aerodynamic force F from F 0 , the longitudinal aerodynamic force induced by control deflections or turbulence is able to be distinguished effectively. The effects of turbulence on vertical acceleration were separated from aircraft maneuvers effects. Therefore, the change of incremental vertical acceleration only induced by turbulence was obtained. Besides, it can be found that the sudden deflections of control surfaces may cause an increase in aerodynamic force and deteriorate the bumpiness severity. Figure 11e shows the EDR indicator by wind-based estimation algorithm given in [39]. The EDR value is commonly output once a minute. In this paper, in order to analyze the dynamic process of EDR and aircraft bumpiness estimation, the 90th percentage EDR value is given in a sub-window of 10 s with 50% overlap. Figure 11f shows the acceleration variation with or without aircraft maneuvers. The turbulence intensity was first obtained by Equation (3) according to the EDR value indicated once a minute. Then the longitudinal and lateral turbulence components were generated by Equation (6) in a sampling frequency of 8Hz. With the turbulence components as the excitation, the time series of vertical acceleration is computed by UVLM. After obtaining the acceleration response, the data length of one-minute vertical acceleration data is 480 points. In the process of bumpiness severity estimation, the band-pass frequency f l and f h are typically chosen as 0.1 and 1.0 Hz respectively for civil aviation aircraft. The one-minute acceleration data are further divided by six sections, while the length of Tukey-Hanning window is set as 160 and the overlapping is set by 50%. Figure 12f shows the total variance estimation,σ 2 a z , which is blended by aircraft maneuvers. In the computation process of incremental acceleration, the effects of control surface deflection are eliminated by σ 2 a z − σ 2 a z0 . It can be concluded that during the severity estimation of aircraft bumpiness, the acceleration response model based on UVLM is capable of effectively avoiding the adverse effects of aircraft maneuvers.

Severity Estimation of Aircraft Bumpiness
In this section, a special section of QAR data was chosen for bumpiness severity estimation of different aircrafts. Figure 13 shows the horizontal flight path from ZBAA to ZGGG. With the recorded QAR data, the wind-based EDR indicators are described by the color gradient. It can be found that in Zone 1, near 36.2 • N (North latitude), 114.7 • E (East longitude), there was moderate turbulence, while in Zone 2, near 30.1 • N,113.6 • E, the aircraft would experience severe turbulence. The one-minute EDR values are shown in the two sub-windows. Among the 880 flight segments, there are 12 different kinds of aircraft type, of which the Boeing 737, Boeing 787, and Airbus 330 aircraft are in the majority. The bumpiness severity of these three aircraft types was analyzed. The average Take-Off Weight (TOW), average Mach number, and grid refinement of the three aircraft types are shown in Table 3. According to the computing flow shown in Figure 7, the bumpiness severity of the above three aircraft types is shown in Figure 14. It can be found that different aircraft would experience different bumpiness severity in the same turbulence field. With the increase of aircraft TOW, the bumpiness severity decreases. In addition, further experiments show that with the increasing cruising airspeed, the bumpiness becomes more intense. According to estimated bumpiness severity indicator beforehand, the pilot can decide whether to fly through the turbulence field or to deviate from the severe turbulence. The bumpiness severity estimation is based on in-situ EDR indication, which is updated once a minute. Thus, the estimated severity would change with the EDR indication per minute. Benefitted by the generated isotropic turbulence components and rigorous algorithm derivations, the estimated value gives a credible bumpiness severity indication per minute.
It should be noted that if the computing results of target aircraft are saved into a table beforehand according to different aircraft mass and airspeed, the bumpiness severity of target aircraft can be estimated by looking up the table with in-situ EDR value. In addition, the bumpiness severity of the aircraft that is about to fly through the turbulence field can also be predicted.

Conclusions
To estimate and prevent severe aircraft bumpiness in advance is beneficial to the flight safety of civil aviation aircraft. This paper put forward a new severity estimation method of aircraft bumpiness based on in-situ EDR indicator. Some innovations are as follows: (1) Based on in-situ EDR indicator, the aircraft-dependent bumpiness severity is estimated. A turbulence time series according to the von Karman model is generated. Based on the assembly of wing and horizontal tail, the unsteady vortex lattice method is proposed to obtain the aerodynamic performance induced by turbulence. After obtaining the time series of vertical acceleration, the bumpiness severity is estimated by Welch spectrum estimation. (2) To improve the computing accuracy, there are several improvements on traditional UVLM. Vortex rings are assigned on the mean camber surface of lifting body, and the lattices neighboring to the wingtip, wing root, and leading edge are refined by semi-circle division method. In addition, when it comes to the control surfaces, the divided lattices are assigned as close as possible to the structural boundary. By the comparison of aircraft modeling data and recorded QAR flight data, it shows that the aerodynamic force change in turbulent flight can be computed with better accuracy by UVLM. (3) By separating turbulence-induced and aircraft maneuvers-induced aerodynamic force change, the adverse effects of control surface deflections on bumpiness severity estimation can be effectively removed. It is also beneficial to estimate in-situ aircraftdependent bumpiness.
This paper is the first attempt to apply the UVLM on the estimation of aircraft bumpiness severity in turbulence. With the proposed method, the bumpiness severity of the next aircraft flying through the same turbulence field can be predicted. In the future, the vortex ring model and related UVLM will be refined to further improve the computing accuracy of vertical acceleration. In addition, by extending and applying the method in Chinese airlines, the practicability of the method can be improved by closed-loop verification with flight data.