Acceleration-Based In Situ Eddy Dissipation Rate Estimation with Flight Data

: Inducing civil aviation aircraft to bumpiness, atmospheric turbulence is a typical risk that seriously threatens ﬂight safety. The Eddy Dissipation Rate (EDR) value, as an aircraft-independent turbulence severity indicator, is estimated by a vertical wind-based or aircraft vertical acceleration-based algorithm. Based on the ﬂight data of civil aviation aircraft, the vertical turbulence component is obtained as the input of both algorithms. A new method of computing vertical acceleration response in turbulence is put forward through the Unsteady Vortex Lattice Method (UVLM). The lifting surface of the target aircraft is assumed to be a combination of wing and horizontal tail in a turbulent ﬂight scenario. Vortex rings are assigned on the mean camber surface, forming a non-planar UVLM, to further improve the accuracy. Moreover, the neighboring vortex lattices are placed as close as possible to the structural edge of control surfaces. Thereby, a complete algorithm for estimating vertical acceleration and in situ EDR value from Quick Access Recorder (QAR) ﬂight data is proposed. Experiments show that the aerodynamic performance is computed accurately by non-planar UVLM. The acceleration response by non-planar UVLM is able to track the recorded acceleration data with higher accuracy than that of the linear model. Di ﬀ erent acceleration responses at di ﬀ erent locations are also obtained. Furthermore, because the adverse e ﬀ ects of aircraft maneuvers are separated from turbulence-induced aircraft bumpiness, the new acceleration-based EDR algorithm shows better accuracy and stability. validation, Z.G., H.W. and Z.X.; formal analysis, H.W.; investigation, H.W.; resources, K.Q.; data curation, K.Q.; writing—original draft preparation, Z.G.; writing—review and editing, Z.G., H.W. and Z.X.; visualization, Z.X.; supervision, Z.G.; project administration, Z.G.; funding acquisition, Z.G.


Introduction
As a chaotic motion added to constant wind, atmospheric turbulence is by far the leading cause of injuries, leading civil aviation aircraft to unexpected bumpiness [1]. On one hand, turbulence deteriorates the riding comfort and results in manipulation difficulty for pilots. On the other hand, significant risks are associated with turbulence, including fuel loss, aircraft structure fatigue or even damage, and human injuries [2]. A great number of studies have been performed to better understand atmospheric turbulence, including the research on turbulence phenomenology and forecast, numerical modeling, in situ or remote detection, and the response of aircraft flying through turbulence. Furthermore, a number of efforts are made to reduce and prevent such risks. It is worth estimating turbulence severity through pilot's subjective report or objective measurement so as to provide a method for risk assessment in turbulent flights [3].
to get solutions with better accuracy [19]. Compared with the non-real time Computational Fluid Dynamics (CFD) method, VLM is a kind of rapid and medium-accurate algorithm that has been widely used in aero-elastic analysis [20], wind turbine design [21], lift and induced drag computation [22], lifting body design optimization [23], and gust wind response analysis [24]. In special applications, such as aerodynamic analysis in aerial refueling formation flight [25][26][27] and innovative flap design [28], VLM possesses the advantage of rapid trial and error. When it comes to the aero-elasticity and gust wind analysis, the Unsteady Vortex Lattice Method (UVLM) is developed with the influence of aircraft maneuvers and wake on aerodynamics considered [20,22,23]. Furthermore, with no-sideslip assumption, VLM shows satisfactory accuracy in developing wind shear coefficients by computing additional aerodynamic effects in wind shear [29,30]. However, the steady VLM is not suitable for computing aerodynamic force in high-frequency turbulence. If the aerodynamic performance with turbulence effects is obtained by UVLM, a more accurate acceleration response can be acquired, and it should be beneficial for acceleration-based EDR estimation.
Based on the Quick Access Recorder (QAR) flight data of civil aviation aircraft, the vertical turbulence component is obtained as an EDR algorithm input. After that, a method for computing acceleration response based on non-planar UVLM and a further acceleration-based EDR estimation is explored in this paper. According to the recorded acceleration, the acceleration response in turbulent flight will be deeply analyzed in both time and frequency domain. Furthermore, compared with wind-based EDR estimation, the performance and accuracy of the new EDR estimation algorithm will be analyzed.

Derived Vertical Wind From QAR Flight Data
To estimate the EDR value, the vertical turbulence component must be first obtained as the input. Real-time flight data are used to obtain in-situ EDR estimation in this paper. The Digital Flight Data Acquisition Unit (DFDAU) in Boeing737-800 aircraft collects and records a great number of flight parameters in the form of data frames via avionic bus [31]. With the support of Chinese domestic airlines, the QAR flight data of the B737-800 aircraft are gathered and used as algorithm input. Moreover, the recorded acceleration data are used for verification. Table 1 compares the parameter requirement of both wind-based and acceleration-based EDR estimation. The vertical turbulence component W z , not recorded in QAR, is the mandatory input for wind-based and acceleration-based EDR estimation. Spatial turbulence is determined as the difference between ground speed and airspeed, where W e = [W x , W y , W z ] T , V e stands for the ground speed vector, C e b represents the transfer matrix from body-axis to earth-axis, and C b a shows the transition from wind-axis to body-axis. In straight and level flight, aircraft sideslip is relatively small and can be ignored for turbulence derivation. By expanding Equation (1), turbulence components are obtained as follows, z − V T (cos α b cos θ cos ψ + sin ψ sin φ sin α b + cos ψ sin θ cos φ sin α b ) W y = V e y − V T (cos α b sin ψ cos θ − cos ψ sin α b sin φ) W z = V e z − V T (sin α b cos θ cos φ − cos α b sin θ) . ( The angle of attack recorded by QAR needs to be converted into the angle of attack in the body axis before being used to obtain vertical wind. In solving Equation (2), the body-axis angle of attack α b is obtained from measured angle of attack α through: The calibration constants a 0 and a 1 are estimated in a preprocessing step using a least-squares linear fit of θ = a 0 + a 1 α. In straight and level flight, α b should be approximately equal to the pitching angle θ under nominal smooth conditions. The other parameters in Equation (2) are obtained directly from QAR. A major advantage of the wind-based EDR estimation is that it requires fewer parameters than the acceleration-based method. When it comes to the acceleration-based EDR estimation proposed in this paper, the angular rate [p, q, r] T is required for obtaining the acceleration response. Since the flight parameters used in Equation (2) are recorded with different sampling rates, from 4 to 16 Hz, the turbulence components are estimated at 4Hz, with the lowest sampling rate.

The von Karman Turbulence Model
The Dryden and von Karman turbulence models are commonly used to describe the stochastic behavior of well-developed and small-scale high-altitude turbulence. However, the von Karman model can describe the characteristics of turbulence with better accuracy. On one hand, the energy spectrum of the von Karman model is built according to large amounts of statistics data and is consistent with the energy spectrum of the Kolmogorov model. The spectrum function of the von Karman model has a roll-off rate of −5/3 in the high-frequency section. However, as an approximation of the von Karman model, the roll-off rate of the Dryden model is −2. On the other hand, when it comes to the EDR estimation, since most of the energy responsible for aircraft bumpiness is in the inertial sub-range, EDR estimation is generally based on an underlying turbulence model. Faster-moving aircraft will also be responsive to larger scales that are typically outside the inertial sub-range. The von Karman model, representing both the inertial sub-range and the larger scales beyond it, has been widely used in both aerodynamics and meteorological communities [32]. As a logical choice in describing turbulence characteristics, the von Karman model has been integrated into both wind-based and acceleration-based EDR estimation. The transverse autocorrelation function of von Karman model is [33]: where L is the length scale, and σ 2 W z is the variance of the vertical turbulence component W z . Γ is the Gamma function and K v is the modified Bessel function. Equation (4) is used to obtain the theoretical power spectrum of turbulence in wind-based EDR estimation. According to von Karman and Kolmogorov energy spectrum theory, the energy spectrum of turbulence only depends on the eddy dissipation rate ε. The EDR value is generally defined as ε 1/3 , which is the cubic root of EDR. The relationship between σ 2 W z and ε is given as follows [8]: where A has been estimated to be 1.6. For the typical flight conditions of civil aviation aircraft, the frequency band within which the majority of vertical acceleration response occurs corresponds to the inertial sub-range of atmospheric turbulence. The velocity spectrum of vertical turbulence φ W z (ω) can be approximated by: where Equation (6) forms the basis of acceleration-based EDR estimation.

Vertical Wind-Based EDR Estimation
With the estimated vertical turbulence component as the input, a frequency-domain and single-parameter maximum likelihood algorithm are used in vertical-wind based EDR estimation. Firstly, a Tukey-Hanning window is employed to reduce the spectrum leakage due to the finite length data of the vertical turbulence component W z . As a well-performed window function, the main-side lobe ratio of the Tukey-Hanning window is higher, and the side lobe converges quickly. The m-point Tukey-Hanning window is formulated as: where M = f loor(0.1m − 0.2). In the more traditional formulation of the Tukey-Hanning window, this corresponds to a taper factor of 0.2. The power normalized Tukey-Hanning window is: Accordingly, the windowed vertical wind data series is: Secondly, the velocity spectrum of the windowed vertical wind series is estimated by Fourier transform:φ where i is the complex imaginary number, f s is the sampling frequency of vertical wind, m = 10 f s corresponds to 10-second wind data, and k = 0, . . . , 5 f s . By dividing the power spectrum of the empirical vertical wind by the theoretical power spectrum over a certain range of frequencies, the estimated EDR valueε 1/3 is obtained by: where γ is a bias-correction term, and k l and k h are the lower and upper index bounds corresponding to the frequencies f l and f h over which the average is taken. In practice, the lower cutoff frequency f l must be greater than the aircraft phugoid response frequency, while the high-end cutoff frequency f h is made to lie below the aero-elastic response mode of the target aircraft. When it comes to a turbulence model spectrum φ model k , the biased autocorrelation function T k for the normalized Tukey-Hanning window τ must be computed firstly by: After that, φ model k is the average periodogram of the windowed von Karman autocorrelation function (shown in Equation (4)) with unit ε 1/3 :

A New Vertical Acceleration-Based EDR Estimation Algorithm
To estimate the aircraft-independent EDR value, the acceleration response excited by turbulence is firstly computed by integrating the non-planar UVLM. After that, a new acceleration-based EDR estimation algorithm is put forward.

Non-Planar Unsteady Vortex Lattice Method
Inviscid, irrotational, and incompressible flow is governed by Laplace's equation, which is derived from the continuity equation: where φ is the velocity potential. Based on linear Laplace's equation, elementary solutions can be superposed to solve the flow around complex geometries. Vortex rings, composed of several vortex segments, are placed on the mean camber surface to model the lifting body. As shown in Figure 1, in order to achieve more accurate results, the vortex sections of AE and BD are further divided by section AF and FE, BC and CD respectively. Taking the vortex segment ED as an example, the induced velocity of ED is obtained with the Biot-Savart equation [19]: Atmosphere 2020, 11, x FOR PEER REVIEW 6 of 20 ( 1) 0

2.4.A New Vertical Acceleration-Based EDR Estimation Algorithm
To estimate the aircraft-independent EDR value, the acceleration response excited by turbulence is firstly computed by integrating the non-planar UVLM. After that, a new acceleration-based EDR estimation algorithm is put forward.

Non-Planar Unsteady Vortex Lattice Method
Inviscid, irrotational, and incompressible flow is governed by Laplace's equation, which is derived from the continuity equation: where φ is the velocity potential. Based on linear Laplace's equation, elementary solutions can be superposed to solve the flow around complex geometries. Vortex rings, composed of several vortex segments, are placed on the mean camber surface to model the lifting body. As shown in Figure 1, in order to achieve more accurate results, the vortex sections of AE and BD are further divided by section AF and FE, BC and CD respectively. Taking the vortex segment ED as an example, the induced velocity of ED is obtained with the Biot-Savart equation [19]: In Equation (15), V is the induced velocity, Γ is the strength of the vortex filament, and 0 r is the vector of the vortex filament.    In Equation (15), V is the induced velocity, Γ is the strength of the vortex filament, and r 0 is the vector of the vortex filament. r 1 and r 2 are the position vectors of both ends to an arbitrary point in space. By defining the coefficient vector k, the induced velocity of a vortex ring to an arbitrary point in space can be expressed as the vector sum of the induced velocity induced by six vortex filaments: Figure 1 also illustrates the boundary condition for the vortex ring, in which the bound vortex is placed on the 1/4 chord of the panel and the Neumann boundary condition of no penetration is enforced at the collocation point located at the 3/4 chord. This is referred to as a fundamental concept for the vortex lattice method, by which the section lift curve slope corresponds exactly to that of thin airfoil theory [34].
In trimmed flight, it is the vertical turbulence component that leads to the fluctuation of longitudinal aerodynamic force and the further change of vertical acceleration. Thus, target aircraft is simplified by a wing-horizontal tail combination, and vortex rings are distributed on their mean camber surfaces. As shown in Figure 2, vortex rings are assigned on the whole reference area rather than the lifting surface only exposed outside the fuselage. The wing and horizontal tail are divided by dozens of rectangular vortex rings, respectively. Moreover, a detailed grid arrangement and further geometric parameters computation are conducted, as shown in the supplementary materials. Figure S1 shows the grid distribution on wing and horizontal tail. Equation S1 gives the definition of each grid number. Equation S2 and S3 give the specific coordinate of each grid on wing and horizontal tail respectively. Since the circulation distribution changes rapidly at the wingtip, wing root, and leading edge, the neighboring lattices are re-designed according to the semi-circle division method [35]. In addition, the boundaries of each vortex ring approximate the geometric edge of control surfaces, since elevators and spoilers can be deflected by the pilot's manipulation or the Automatic Flight Control System (AFCS).
Atmosphere 2020, 11, x FOR PEER REVIEW 7 of 20 enforced at the collocation point located at the 3/4 chord. This is referred to as a fundamental concept for the vortex lattice method, by which the section lift curve slope corresponds exactly to that of thin airfoil theory [34].
In trimmed flight, it is the vertical turbulence component that leads to the fluctuation of longitudinal aerodynamic force and the further change of vertical acceleration. Thus, target aircraft is simplified by a wing-horizontal tail combination, and vortex rings are distributed on their mean camber surfaces. As shown in Figure 2, vortex rings are assigned on the whole reference area rather than the lifting surface only exposed outside the fuselage. The wing and horizontal tail are divided by dozens of rectangular vortex rings, respectively. Moreover, a detailed grid arrangement and further geometric parameters computation are conducted, as shown in the supplementary materials. Figure S1 shows the grid distribution on wing and horizontal tail. Equation S1 gives the definition of each grid number. Equation S2 and S3 give the specific coordinate of each grid on wing and horizontal tail respectively. Since the circulation distribution changes rapidly at the wingtip, wing root, and leading edge, the neighboring lattices are re-designed according to the semi-circle division method [35]. In addition, the boundaries of each vortex ring approximate the geometric edge of control surfaces, since elevators and spoilers can be deflected by the pilot's manipulation or the Automatic Flight Control System (AFCS). In the unsteady vortex lattice method, a time-marching scheme is applied to perform unsteady simulation. At each time step k t , a row of newly generated vortex rings is released into the wake, with the same circulation as the trailing-edge panel from which it was shed at the previous time step. As shown in Figure 3, the circulation Γ of the newly shed wake panel stays unchanged for the remainder of the simulation as dictated by Helmholtz's theorem [19]. The wing-tail combination releases force-free wakes with no aerodynamic loads. In the unsteady vortex lattice method, a time-marching scheme is applied to perform unsteady simulation. At each time step t k , a row of newly generated vortex rings is released into the wake, with the same circulation as the trailing-edge panel from which it was shed at the previous time step. As shown in Figure 3, the circulation Γ of the newly shed wake panel stays unchanged for the remainder of the simulation as dictated by Helmholtz's theorem [19]. The wing-tail combination releases force-free wakes with no aerodynamic loads. Although the turbulence field is rotational, the scale of turbulence is large enough compared to the aircraft. So, it is feasible to analyze the boundary condition and further compute the aerodynamic force by UVLM according to the turbulence variation at each time step. Under the effects of airfoilattached vortices, wake vortices, and local velocity, the Neumann boundary condition at any collocation point is: By further expanding the above equation, a linear algebraic equation is obtained. The circulation of each vortex ring is obtained by 1 It should be noted that in turbulent flight, in addition to the far-field free flow and aircraft angular motion has major effects on local velocity. The local velocity of an arbitrary point at any time k t is described by Equation (18), with [ , , ] T x y z the coordinate of arbitrary point: After obtaining the circulation of each vortex ring, the aerodynamic force is further computed according to the Kutta-Joukowski theorem. As shown in Figure 4, any vortex and its adjacent vortices need to be superposed on the lattice panel to get actual circulation distribution. Taking the right half wing as an example, the aerodynamic force at the collocation point is: where ρ is air density, and , bi V is the total induced velocity of the airfoil-attached vortex at the collocation point. At each time step k t , the aerodynamic force with or without turbulence effects are solved respectively, which are recorded as i F and 0 i F . The total aerodynamic force contains the sum of all vortex rings on the wing and horizontal tail. The pitching moment around the gravity center is obtained by: where ( , ) cg cg x z represent the location of the gravity center. Although the turbulence field is rotational, the scale of turbulence is large enough compared to the aircraft. So, it is feasible to analyze the boundary condition and further compute the aerodynamic force by UVLM according to the turbulence variation at each time step. Under the effects of airfoil-attached vortices, wake vortices, and local velocity, the Neumann boundary condition at any collocation point is: By further expanding the above equation, a linear algebraic equation is obtained. The circulation of each vortex ring is obtained by ] T has major effects on local velocity. The local velocity of an arbitrary point at any time t k is described by Equation (18), with [x, y, z] T the coordinate of arbitrary point: After obtaining the circulation of each vortex ring, the aerodynamic force is further computed according to the Kutta-Joukowski theorem. As shown in Figure 4, any vortex and its adjacent vortices need to be superposed on the lattice panel to get actual circulation distribution. Taking the right half wing as an example, the aerodynamic force at the collocation point is: where ρ is air density, and V b,i is the total induced velocity of the airfoil-attached vortex at the collocation point. At each time step t k , the aerodynamic force with or without turbulence effects are solved respectively, which are recorded as F i and F 0 i . The total aerodynamic force F = [F x , F y , F z ] T contains the sum of all vortex rings on the wing and horizontal tail. The pitching moment around the gravity center is obtained by: where (x cg , z cg ) represent the location of the gravity center.
Atmosphere 2020, 11, x FOR PEER REVIEW 9 of 20 The differential equation of vertical acceleration and pitching moment induced by turbulence are: where ( ) . In addition, vertical acceleration at any location on the longitudinal plane is obtained by: It should be noted that a gap of longitudinal aerodynamic performance exists inevitably between real aircraft and the wing-tail combination. However, this paper makes a compromise between the real-time performance and the algorithm accuracy. On one hand, compared with the fuselage, the aircraft wing and horizontal tail play a major role in the pitching and plunge motion, while the integration of fuselage will greatly increase the complexity of the algorithm. On the other hand, only the acceleration increment caused by vertical turbulence needs to be considered for EDR estimation. A satisfactory accuracy of further EDR estimation is able to be guaranteed by the incremental acceleration response.

Acceleration-Based EDR Estimation
In turbulent flight, the spatial frequency and intensity of turbulence are actually obtained by measuring the temporal frequency and aircraft acceleration response. Given the frequency response function from turbulence input ( ) z W iω to vertical acceleration ( ) z a iω , the spectrum of the acceleration response is obtained by: Equation (23) is built up in the temporal-frequency domain, since the measurements are made by aircraft flying through turbulence with airspeed V. According to Taylor's frozen-field hypothesis, The differential equation of vertical acceleration and pitching moment induced by turbulence are: where θ(t k ) and φ(t k ) are the pitching and rolling angle at t k , m is the mass of aircraft, and stands for the inertia matrix. In order to separate aircraft maneuvers effects from acceleration response, Equation (21) needs to be computed twice. Firstly, the total aerodynamic force F is used in Equation (21) to get a complete acceleration response a z (t k ), which is induced by aircraft maneuvers and external turbulence. Secondly, the aerodynamic force without turbulence effects F 0 functions as input to get the acceleration only containing maneuver effects a 0 z (t k ). The incremental acceleration only induced by turbulence is obtained by a z (t k ) − a 0 z (t k ). In addition, vertical acceleration at any location on the longitudinal plane is obtained by: It should be noted that a gap of longitudinal aerodynamic performance exists inevitably between real aircraft and the wing-tail combination. However, this paper makes a compromise between the real-time performance and the algorithm accuracy. On one hand, compared with the fuselage, the aircraft wing and horizontal tail play a major role in the pitching and plunge motion, while the integration of fuselage will greatly increase the complexity of the algorithm. On the other hand, only the acceleration increment caused by vertical turbulence needs to be considered for EDR estimation. A satisfactory accuracy of further EDR estimation is able to be guaranteed by the incremental acceleration response.

Acceleration-Based EDR Estimation
In turbulent flight, the spatial frequency and intensity of turbulence are actually obtained by measuring the temporal frequency and aircraft acceleration response. Given the frequency response function from turbulence input W z (iω) to vertical acceleration a z (iω), the spectrum of the acceleration response is obtained by: Equation (23) is built up in the temporal-frequency domain, since the measurements are made by aircraft flying through turbulence with airspeed V. According to Taylor's frozen-field hypothesis, the length scale of the vertical turbulence component L can be obtained by L = ω/V. The estimation is conducted in the frequency-domain, and the above equation is integrated between the frequency segments ω l and ω f to obtain the acceleration response energy: Similar to the selection of lower and upper index bounds in wind-based EDR estimation, a band-pass filter H bp (ω l , ω h , ω) should be cascaded to eliminate the influence of aircraft maneuvers and high-frequency aero-elastic vibration, with ω l = 2π f l and ω h = 2π f h . The above equation can be modified by: To estimate the response energy, the band-pass filter is transformed into time-domain firstly by Inverse Fourier transform (IFT) as follows: After that, Parseval's theorem (assuming local stationary) and the convolution theorem are adopted to yield the approximation [36]: where a z (t) is the acceleration data series obtained by Equation (21). It states that the mean-square value of the frequency-limited response is obtained by computing a running mean-square value of the band-pass filtered acceleration data. A nonzero mean value must be removed in the computation of the outer integration in Equation (27). The averaging integral 2T is selected to satisfy the local-stationarity assumption and provide enough data values for statistical stability. The time constant τ 1 is chosen such that h bp (τ) → 0 for τ > τ 1 . By inserting Equation (6) into Equation (27), If I(ω l , ω h , t) in Equation (28) is known, the EDR value ε 1/3 will be estimated. To estimate the acceleration response, Fourier transform is performed on the vertical acceleration and the vertical wind sequence in the period T. After that, it is converted from the time-domain to the frequency-domain to obtain a z (iω) and W z (iω), respectively. With the new time series a z (iω) W z (iω) , the autocorrelation function R is obtained by IFT. According to Wiener-Khinchin theorem, R(0) is taken as the estimation results of I(ω l , ω h , t k ). As a result, the EDR value is estimated by: To summarize, an integrated algorithm flow for the new acceleration-based EDR estimation is shown in Figure 5. The grid on wing-tail combination is generated beforehand. In the algorithm loop, the measured flight parameters are used to obtain turbulence components by Equation (2) and further local velocity is obtained by Equation (18). The parameters shown in Table 1 function as the data input.
After that, the aerodynamic force is acquired by non-planar UVLM. In each cycle, by toggling the flag, the total aerodynamic force F and the force without turbulence effects F 0 are obtained, respectively. The vertical acceleration increment induced by turbulence is obtained by solving Equation (21) twice. Next, the time series of acceleration data are used to estimateσ a z (t k ) by Equation (27). In addition, after Fast Fourier Transform (FFT), the acceleration series in the frequency-domain are used to obtain the transfer function I(ω l , ω h , t k ) together with the frequency-domain vertical turbulence component. At the end of each iteration, the estimated EDR value is output by solving Equation (29).

Analysis of Grid Convergence
Before the acceleration response analysis, the performance of non-planar UVLM was tested. UVLM requires a sufficiently refined grid to achieve results that are grid-independent. The grid resolution of UVLM was first studied to improve the computing efficiency on the premise of accuracy assurance.
As far as Boeing 737-800 aircraft is concerned, the geometry and airfoil parameters of the wing and horizontal tail are listed in ref [37,38]. A non-planar vortex lattice model was built, in which different surface grids of increasing mesh density were generated. Vortex rings are assigned on the mean camber surface. Each grid has a constant spacing in both chord-wise and span-wise directions. The changes of the lift coefficient, drag coefficient, and pitching moment coefficient are presented as a function of the grid refinement level. As a result, satisfactory accuracy is obtained by adopting a

Analysis of Grid Convergence
Before the acceleration response analysis, the performance of non-planar UVLM was tested. UVLM requires a sufficiently refined grid to achieve results that are grid-independent. The grid resolution of UVLM was first studied to improve the computing efficiency on the premise of accuracy assurance.
As far as Boeing 737-800 aircraft is concerned, the geometry and airfoil parameters of the wing and horizontal tail are listed in ref [37,38]. A non-planar vortex lattice model was built, in which different surface grids of increasing mesh density were generated. Vortex rings are assigned on the mean camber surface. Each grid has a constant spacing in both chord-wise and span-wise directions. The changes of the lift coefficient, drag coefficient, and pitching moment coefficient are presented as a function of the grid refinement level. As a result, satisfactory accuracy is obtained by adopting a grid with 40 columns, 20 rows on the wing and 12 columns, 6 rows on the horizontal tail. Further experiments are based on the optimized grid level.

Computation of Pitching Moment Coefficient
There exists a gap of longitudinal aerodynamic characteristics between the whole aircraft and the wing-tail combination. However, since the pitching moment change is mainly caused by the deflection of elevator or stabilizer, the computing results based on the wing-tail combination should be comparable with authentic aerodynamic data. The modeling data of B737-800 aircraft provides longitudinal pitching moment changes caused by elevator and stabilizer deflections [37,38]. Figure 6 shows the results between UVLM results and the modeling data in three flight conditions. c mδstab stands for the pitching moment change due to stabilizer deflection, and c mδe represents the moment change due to elevator deflection. grid with 40 columns, 20 rows on the wing and 12 columns, 6 rows on the horizontal tail. Further experiments are based on the optimized grid level.

Computation of Pitching Moment Coefficient
There exists a gap of longitudinal aerodynamic characteristics between the whole aircraft and the wing-tail combination. However, since the pitching moment change is mainly caused by the deflection of elevator or stabilizer, the computing results based on the wing-tail combination should be comparable with authentic aerodynamic data. The modeling data of B737-800 aircraft provides longitudinal pitching moment changes caused by elevator and stabilizer deflections [37,38]. Figure 6 shows the results between UVLM results and the modeling data in three flight conditions. m stab c δ stands for the pitching moment change due to stabilizer deflection, and m e c δ represents the moment change due to elevator deflection. The experiments prove that based on the wing-tail combination, the non-planar UVLM has satisfactory computing accuracy, which provides a good basis for further estimation of aircraft acceleration and turbulence severity. Furthermore, the accurate computation of the pitching moment coefficient is significant to separate aircraft maneuvers effects on EDR estimation.

Acceleration Response Analysis
Since accurate acceleration response forms the basis for EDR estimation, the acceleration response by non-planar UVLM will be analyzed and further compared with measured acceleration in QAR. Supported by Chinese airlines, QAR flight data including 400 flight segments have been gathered from B737-800 aircraft on the same scheduled air route. Five-minute cruising flight data in Compared with modeling data, the mean square error (MSE) of c mδstab is 0.0045 and the MSE of c mδe is 0.0055. It shows that the pitching moment change caused by the deflection of stabilizer or elevator approximates the modeling data under the above flight conditions. Since the geometrical boundary of the control surface approaches the boundary of each vortex ring, non-planar UVLM can be used to accurately describe the pitching moment change.
The experiments prove that based on the wing-tail combination, the non-planar UVLM has satisfactory computing accuracy, which provides a good basis for further estimation of aircraft acceleration and turbulence severity. Furthermore, the accurate computation of the pitching moment coefficient is significant to separate aircraft maneuvers effects on EDR estimation.

Acceleration Response Analysis
Since accurate acceleration response forms the basis for EDR estimation, the acceleration response by non-planar UVLM will be analyzed and further compared with measured acceleration in QAR. Supported by Chinese airlines, QAR flight data including 400 flight segments have been gathered from B737-800 aircraft on the same scheduled air route. Five-minute cruising flight data in each segment with the same geographic area are selected for EDR estimation. According to pilot reports, six Wind Fields (WF) covering light, moderate, and severe turbulence are selected to conduct the comparative study of acceleration response. The turbulence components are computed by Equation (2) as the input of the following EDR estimation. It should be noted that there are acceleration data in the QAR of B737-800 aircraft, while not all the aircraft have the acceleration record.
The acceleration response computed by non-planar UVLM is compared with the linear model and recorded QAR data. The peak vertical acceleration data of the six wind fields (WF1-6) are shown in Table 2. It is worth noting that the acceleration data is recorded at the aircraft gravity center.

Comparison of Acceleration Response
Referring to model-4 in ref [15], a linear model of B737-800 aircraft was built based on the modeling data. To better understand the gust penetration effect and the mechanism underlying the gust pitch rate term, four linear models of increasing gust excitation complexity were developed in ref [15]. Model-1 is the simplest model by using a gust point approximation. Model-2 considers gust acceleration and gust pitch rate effects in the pitch dynamics. Furthermore, by accounting for the full-body gust penetration effects, model-3 incorporates these gust excitation signals in the plunge dynamics. Model-4, incorporating a quasi-steady aerodynamic model with finite lag value, is able to follow the acceleration transients with better accuracy, while the peak prediction accuracy is sacrificed. Figure 7a shows the derived vertical wind, and Figure 7b shows a short sample of recorded vertical acceleration transients overlaid with the computation results from linear model and non-planar UVLM in WF3. Up to 166s, the data imply the aircraft was at a trimmed flight condition. At 167.5 s, the aircraft reacted to the initial disturbance, reached the maximum positive load near 168.3 s, and then started another local maximum positive load at 170s. After that, it experienced a continuous decrease. In this short sequence, the linear model was momentarily out-of-phase once encountering sudden turbulence, while the UVLM showed better tracking accuracy. Although not a conclusive argument, the observation from Figure 7 suggests that both the linear model and UVLM are able to replicate the acceleration response in moderate turbulence. However, the tracking accuracy of the UVLM is better than that of the linear model. each segment with the same geographic area are selected for EDR estimation. According to pilot reports, six Wind Fields (WF) covering light, moderate, and severe turbulence are selected to conduct the comparative study of acceleration response. The turbulence components are computed by Equation (2) as the input of the following EDR estimation. It should be noted that there are acceleration data in the QAR of B737-800 aircraft, while not all the aircraft have the acceleration record. The acceleration response computed by non-planar UVLM is compared with the linear model and recorded QAR data. The peak vertical acceleration data of the six wind fields (WF1-6) are shown in Table 2. It is worth noting that the acceleration data is recorded at the aircraft gravity center. Referring to model-4 in ref [15], a linear model of B737-800 aircraft was built based on the modeling data. To better understand the gust penetration effect and the mechanism underlying the gust pitch rate term, four linear models of increasing gust excitation complexity were developed in ref [15]. Model-1 is the simplest model by using a gust point approximation. Model-2 considers gust acceleration and gust pitch rate effects in the pitch dynamics. Furthermore, by accounting for the fullbody gust penetration effects, model-3 incorporates these gust excitation signals in the plunge dynamics. Model-4, incorporating a quasi-steady aerodynamic model with finite lag value, is able to follow the acceleration transients with better accuracy, while the peak prediction accuracy is sacrificed. Figure 7a shows the derived vertical wind, and Figure 7b shows a short sample of recorded vertical acceleration transients overlaid with the computation results from linear model and nonplanar UVLM in WF3. Up to 166s, the data imply the aircraft was at a trimmed flight condition. At 167.5 s, the aircraft reacted to the initial disturbance, reached the maximum positive load near 168.3 s, and then started another local maximum positive load at 170s. After that, it experienced a continuous decrease. In this short sequence, the linear model was momentarily out-of-phase once encountering sudden turbulence, while the UVLM showed better tracking accuracy. Although not a conclusive argument, the observation from Figure 7 suggests that both the linear model and UVLM are able to replicate the acceleration response in moderate turbulence. However, the tracking accuracy of the UVLM is better than that of the linear model.  Figure 8 shows the acceleration response of the two models for WF6 along with the recorded data. This sample covers a severe vertical turbulence event resulting in an amplitude swinging from +2.124 g to +0.175 g over 1s. Grouping together when responding to the turbulence, the two responses tended to track the overall behavior of the recorded data. The sacrifice of peak acceleration response accuracy brought about the serious deviation in the trace of linear model from recorded data. On the contrary, UVLM performed higher sensitivity to severe turbulence with sharp gradients resulting in an over-prediction of peak values. The tracking MSE was reduced from 0.0247 (linear model) to 0.0081 (UVLM). Qualitative observation leads to the conclusion that the linear model is less sensitive for moderate and severe turbulence. It can be concluded that the recorded acceleration is better followed by UVLM than by the linear model.  Figure 8 shows the acceleration response of the two models for WF6 along with the recorded data. This sample covers a severe vertical turbulence event resulting in an amplitude swinging from +2.124 g to +0.175 g over 1s. Grouping together when responding to the turbulence, the two responses tended to track the overall behavior of the recorded data. The sacrifice of peak acceleration response accuracy brought about the serious deviation in the trace of linear model from recorded data. On the contrary, UVLM performed higher sensitivity to severe turbulence with sharp gradients resulting in an over-prediction of peak values. The tracking MSE was reduced from 0.0247 (linear model) to 0.0081 (UVLM). Qualitative observation leads to the conclusion that the linear model is less sensitive for moderate and severe turbulence. It can be concluded that the recorded acceleration is better followed by UVLM than by the linear model. The result in Figure 8 is not unexpected, that a higher-order feature in recorded acceleration data is reliably replicated by UVLM rather than the linear model. There are significant differences between the two models. Firstly, severe turbulence could easily change the governing aerodynamic coefficients of the linear model from trim condition. In contrast, UVLM is a kind of nonlinear computation, instead of linear approximation. Secondly, the recorded acceleration response is actually excited by external turbulence and control surface deflections by manual manipulation or AFCS. The recorded control deflections function as input in UVLM, while aircraft maneuvers are not considered in linear model.
For a more comprehensive comparison, the MSE and absolute peak acceleration in WF1-WF6 are presented graphically in Figure 9. The linear model does a good prediction in light turbulence while suffering accuracy in moderate and severe turbulence. The greater the turbulence intensity, the bigger the MSE will be. However, UVLM is capable of predicting the acceleration response with better accuracy than the linear model. The satisfactory prediction of peak acceleration holds promise for further EDR estimation. The result in Figure 8 is not unexpected, that a higher-order feature in recorded acceleration data is reliably replicated by UVLM rather than the linear model. There are significant differences between the two models. Firstly, severe turbulence could easily change the governing aerodynamic coefficients of the linear model from trim condition. In contrast, UVLM is a kind of nonlinear computation, instead of linear approximation. Secondly, the recorded acceleration response is actually excited by external turbulence and control surface deflections by manual manipulation or AFCS. The recorded control deflections function as input in UVLM, while aircraft maneuvers are not considered in linear model.
For a more comprehensive comparison, the MSE and absolute peak acceleration in WF1-WF6 are presented graphically in Figure 9. The linear model does a good prediction in light turbulence while suffering accuracy in moderate and severe turbulence. The greater the turbulence intensity, the bigger the MSE will be. However, UVLM is capable of predicting the acceleration response with better accuracy than the linear model. The satisfactory prediction of peak acceleration holds promise for further EDR estimation.

Spectrum Analysis at Different Locations
Three typical locations are chosen to investigate the acceleration response at different locations on aircraft. Location A is at the gravity center, and Location B is at the cockpit position, which is about 15 feet from point A. Location C is at the tailgate of aircraft, about 60 feet away from the gravity center. One of the prominent features of UVLM is that the acceleration at different locations of the fuselage can be obtained. Since it is nearly impossible to find the differences from time-domain acceleration series, the spectral density estimation is used to obtain the frequency-domain acceleration spectrum at different locations. The acceleration spectra of the three locations in WF1-WF6 are shown in Figure 10. To obtain the spectrum densities of the vertical acceleration series, the Welch spectral density estimation with the Hanning window function is used [35]. The data length in the time-domain for each curve is 2,560 points, which represents a length of 640 s acceleration data. With the sampled acceleration data with 4Hz, the window length is 10 × fs = 40, and the overlapping is set by 50%.

Spectrum Analysis at Different Locations
Three typical locations are chosen to investigate the acceleration response at different locations on aircraft. Location A is at the gravity center, and Location B is at the cockpit position, which is about 15 feet from point A. Location C is at the tailgate of aircraft, about 60 feet away from the gravity center. One of the prominent features of UVLM is that the acceleration at different locations of the fuselage can be obtained. Since it is nearly impossible to find the differences from time-domain acceleration series, the spectral density estimation is used to obtain the frequency-domain acceleration spectrum at different locations. The acceleration spectra of the three locations in WF1-WF6 are shown in Figure 10. To obtain the spectrum densities of the vertical acceleration series, the Welch spectral density estimation with the Hanning window function is used [35]. The data length in the time-domain for each curve is 2,560 points, which represents a length of 640 s acceleration data. With the sampled acceleration data with 4Hz, the window length is 10 × fs = 40, and the overlapping is set by 50%. It can be found that the bumpiness at the tail of the fuselage is larger than that at the gravity center, while the bumpiness at the nose is the least. It is consistent with the analysis in ref [17]. In moderate turbulence, there is little difference in magnitudes of vertical acceleration between locations It can be found that the bumpiness at the tail of the fuselage is larger than that at the gravity center, while the bumpiness at the nose is the least. It is consistent with the analysis in ref [17]. In moderate turbulence, there is little difference in magnitudes of vertical acceleration between locations A, B, and C. However, with the increase of turbulence severity, the difference of acceleration magnitude increases. As a result, the UVLM-based acceleration response algorithm offers an effective way for estimating the bumpiness at different locations of the fuselage.

EDR Estimation Comparison
After that, the new acceleration-based EDR estimation is tested with real flight data. The algorithm accuracy, especially after separating aircraft maneuvers effects, is studied and compared with wind-based estimation.
With the new acceleration response algorithm, a further acceleration-based EDR estimation is explored. The choice of window-length and EDR update rate is a compromise between having enough resolution to capture instantaneous aircraft bumpiness and having enough samples in the window to provide stable computational statistics. For cruising aircraft, a nominal 1-minute reporting interval is generally used. For each of the 5-minute flight time in 400 flight segments, there are totally 2000 EDR values. Within the reporting window, a 10-second sub-window with 1/2 overlap is selected for EDR estimation in both the wind-based and acceleration-based method. Furthermore, the median and 90th percentage ε 1/3 values are computed respectively in the 1-minute period. In order to mitigate potential erroneous EDR output due to QAR data quality issues, the average and peak EDR value are replaced by the median and 90th percentage value. If the turbulence is relatively continuous, the median and 90th percentage ε 1/3 values should be close to each other in magnitude. Whereas for a discrete sudden bumpiness, the 90th percentage value will typically be larger than the median value. In addition, the band-pass frequencies ω l and ω h are typically chosen as 0.1 and 1.0 Hz respectively for civil aviation aircraft. Figure 11 shows the EDR estimation results in WF6, along with the vertical turbulence component, control surface deflections, and aerodynamic force change. Taking the QAR data in WF6 as an example, the vertical turbulence component obtained by Equation (2) is shown in Figure 11a. From 90s to 100s, the aircraft experienced a sudden downdraft. Meanwhile, the deflection of elevators and spoilers by pilot manipulation or AFCS are shown in Figure 11b. It can be found that the low-speed flight spoiler δ s LS and high-speed flight spoiler δ s HS both experienced sudden deflections. Although the deflection of the elevator was smaller than that of spoilers, it had greater impact on the aerodynamics force. In non-planar UVLM, since the boundaries of the vortex ring are close to the geometric edge of the control surface, the change of aerodynamic force due to control deflection can be computed accurately.

Comparison of EDR Estimation
Based on the QAR data of 400 flight segments, Figure 12 shows a scatterplot of the total 2000 EDR values for wind-based versus acceleration-based EDR estimation. Figure 12a shows the scatterplot of one-minute median EDR values, whereas Figure 12b shows the 90th percentile EDR values. For each point, the horizontal coordinate stands for the new acceleration-based EDR value, and the vertical coordinate represents the wind-based EDR value. If both estimation methods recover the same aircraft-independent EDR value, the estimated EDR value points should distribute neighboring to the one-to-one line, which is clearly identified by y = x in Figure 12. However, a small number of points are clustered above the one-to-one line, at the region of small estimated EDR values (lower-left corners). These points are due to the incomplete removal of all the maneuvers-induced accelerations in the wind-based method. The band-pass filter with empirically selected band-width  Figure 11c shows the aerodynamic force with or without control deflections. With the proposed method in this paper, the longitudinal aerodynamic force induced by aircraft maneuvers or vertical turbulence is able to be effectively distinguished. The effects of vertical turbulence on vertical acceleration were separated from the aerodynamic force F 0 . Therefore, the change of incremental acceleration was only excited by vertical turbulence. It can be found that pilot manipulation may lead to an increase in aerodynamic force and deteriorate aircraft bumpiness.
Despite being filtered by a band-pass filter, aircraft maneuvers induced by control deflection still have adverse effects on traditional EDR estimation, as shown in Figure 11d. In order to show the dynamic process of EDR estimation, the 90th percentage EDR value is given in a sub-window of 10 s with 1/2 overlap. Through the computation of incremental acceleration, the influence of control deflection-induced maneuvers is eliminated. It can be concluded that the new EDR estimation is able to effectively eliminate the adverse effects of aircraft maneuvers.

Comparison of EDR Estimation
Based on the QAR data of 400 flight segments, Figure 12 shows a scatterplot of the total 2000 EDR values for wind-based versus acceleration-based EDR estimation. Figure 12a shows the scatterplot of one-minute median EDR values, whereas Figure 12b shows the 90th percentile EDR values. For each point, the horizontal coordinate stands for the new acceleration-based EDR value, and the vertical coordinate represents the wind-based EDR value. If both estimation methods recover the same aircraft-independent EDR value, the estimated EDR value points should distribute neighboring to the one-to-one line, which is clearly identified by y = x in Figure 12. However, a small number of points are clustered above the one-to-one line, at the region of small estimated EDR values (lower-left corners). These points are due to the incomplete removal of all the maneuvers-induced accelerations in the wind-based method. The band-pass filter with empirically selected band-width parameters f h could not eliminate maneuvers effects completely, leading the wind-based EDR estimates to be more sensitive to control deflections. A first-order function fitting was made over the data points, with the fitting equation y = 1.0719x + 0.0014395 and the correlation coefficient 0.95568. It means that the wind-based EDR value is statistically higher than the acceleration-based EDR value because of aircraft maneuvers contamination. For the 90th percentile values, the fitting function is y = 1.1854x + 0.002285, and the correlation coefficient is 0.97516. It means that in wind-based estimation, the discrete sudden bumpiness makes the 90th percentage EDR value much higher than the acceleration-based EDR value. With the above non-planar UVLM configuration, the CPU time of vertical acceleration computation, which mainly contributes to the complexity of algorithm, costs 1.96s with an i7-9700 processor under a visual studio 2017 development environment. The whole CPU time of a one-time complete EDR estimation is 2.47s. As a result, the algorithm is capable of outputing the EDR value in real time. After being further optimized and embedded into the airborne system, the in situ EDR estimation is realized.

Conclusions
The accurate estimation of atmospheric turbulence severity is of great significance to flight safety. The Eddy Dissipation Rate (EDR) indicator has been widely used in turbulence severity measurement. This paper put forward a new method for aircraft vertical acceleration and EDR estimation. Some innovations are as follows: 1. Based on wing and horizontal tail combination of target aircraft, the non-planar Unsteady Vortex Lattice Method (UVLM) is proposed to obtain the aerodynamic force change in turbulent flight. With the increase of turbulence intensity, the number of discrete sudden bumpiness increases. The variance of EDR value by the wind-based method increases, and the number of outliers increases. Benefitted by accurate vertical acceleration response, the acceleration-based EDR estimation is still satisfactory.
With the above non-planar UVLM configuration, the CPU time of vertical acceleration computation, which mainly contributes to the complexity of algorithm, costs 1.96s with an i7-9700 processor under a visual studio 2017 development environment. The whole CPU time of a one-time complete EDR estimation is 2.47s. As a result, the algorithm is capable of outputing the EDR value in real time.
After being further optimized and embedded into the airborne system, the in situ EDR estimation is realized.

Conclusions
The accurate estimation of atmospheric turbulence severity is of great significance to flight safety. The Eddy Dissipation Rate (EDR) indicator has been widely used in turbulence severity measurement. This paper put forward a new method for aircraft vertical acceleration and EDR estimation. Some innovations are as follows: 1.
Based on wing and horizontal tail combination of target aircraft, the non-planar Unsteady Vortex Lattice Method (UVLM) is proposed to obtain the aerodynamic force change in turbulent flight.
To improve the computing accuracy of UVLM, vortex rings are assigned on the mean camber surface, and the semi-circle division method is adopted to refine the lattices neighboring to wingtip, wing root, and leading edge. Another major improvement is that the lattices neighboring to the control surfaces, such as spoilers and elevators, are divided as close as possible to the structural edge, which is beneficial for further EDR estimation.

2.
A complete algorithm flow for estimating vertical acceleration and in situ EDR value from QAR flight data is proposed. Experiments on acceleration response show that compared with a linear model, there are improvements in both response accuracy and tracking performance in moderate and severe bumpiness. The vertical accelerations at different locations of the fuselage are computed accurately. Compared with wind-based EDR estimation, the new acceleration-based EDR algorithm shows better accuracy and stability because the adverse influence of aircraft maneuvers on EDR estimation is eliminated effectively.

Conflicts of Interest:
The authors declare no conflict of interest.