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

Transcription

1 atmosphere Article Acceleration-Based In Situ Eddy Dissipation Rate Estimation with Flight Data Zhenxing Gao 1 * Haofeng Wang Kai Qi 3 Zhiwei Xiang 1 and Debao Wang 1 1 College of Civil Aviation Nanjing University of Aeronautics and Astronautics Nanjing 1116 China; xiang_zw@nuaa.edu.cn (Z.X.); wdb@nuaa.edu.cn (D.W.) Aviation Safety Institute China Academy of Civil Aviation Science and Technology Beijing 18 China; wanghf@mail.castc.org.cn 3 Department of Flight Shandong Airlines Jinan 514 China; qi@sda.cn * Correspondence: z.x.gao@nuaa.edu.cn; Tel.: Received: 11 October ; Accepted: 16 November ; Published: 19 November Abstract: Inducing civil aviation aircraft to bumpiness atmospheric turbulence is a typical ris that seriously threatens flight 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 flight data of civil aviation aircraft vertical turbulence component is obtained as input of both algorithms. A new method of computing vertical acceleration response in turbulence is put forward through Unsteady Vortex Lattice Method (UVLM). The lifting surface of target aircraft is assumed to be a combination of wing and horizontal tail in a turbulent flight scenario. Vortex rings are assigned on mean camber surface forming a non-planar UVLM to furr improve accuracy. Moreover neighboring vortex lattices are placed as close as possible to structural edge of control surfaces. Thereby a complete algorithm for estimating vertical acceleration and in situ EDR value from Quic Access Recorder (QAR) flight data is proposed. Experiments show that aerodynamic performance is computed accurately by non-planar UVLM. The acceleration response by non-planar UVLM is able to trac recorded acceleration data with higher accuracy than that of linear model. Different acceleration responses at different locations are also obtained. Furrmore because adverse effects of aircraft maneuvers are separated from turbulence-induced aircraft bumpiness new acceleration-based EDR algorithm shows better accuracy and stability. Keywords: turbulence; Eddy Dissipation Rate (EDR); aircraft acceleration; Unsteady Vortex Lattice Model (UVLM); vortex ring 1. Introduction As a chaotic motion added to constant wind atmospheric turbulence is by far leading cause of injuries leading civil aviation aircraft to unexpected bumpiness [1]. On one hand turbulence deteriorates riding comfort and results in manipulation difficulty for pilots. On or hand significant riss are associated with turbulence including fuel loss aircraft structure fatigue or even damage and human injuries []. A great number of studies have been performed to better understand atmospheric turbulence including research on turbulence phenomenology and forecast numerical modeling in situ or remote detection and response of aircraft flying through turbulence. Furrmore a number of efforts are made to reduce and prevent such riss. It is worth estimating turbulence severity through pilot s subjective report or objective measurement so as to provide a method for ris assessment in turbulent flights [3]. Atmosphere ; doi:1.339/atmos

2 Atmosphere of Airborne observations of turbulence can be classified into pilot wear reports (PIREPs) and aircraft-based measurements. As a ind of direct and subjective method PIREPs have been widely used because of ir wide coverage in space and time following pilots participation in reporting meteorological conditions during flight [4]. However subjective judgment of turbulence leads to a high degree of uncertainty in terms of intensity timing and location of turbulence. Therefore use of in situ airborne data has been strongly demanded to obtain more accurate observational information regarding turbulence [5]. As far as indicators of turbulence severity are concerned re are Vertical Acceleration (VA) [6] Derived Equivalent Vertical Gust velocity (DEVG) [7] and Eddy Dissipation Rate (EDR) [8]. The commonly used VA indicator is RMS-g which is Root Mean Square of vertical acceleration divided by 1-g [9]. Since different aircraft flying through same turbulence field would have different acceleration responses directly measured RMS-g cannot reflect objective turbulence severity. Although with airspeed and aircraft mass considered improved DEVG still cannot eliminate flight effects on severity estimation satisfactorily [7]. To address above deficiencies EDR estimation algorithm based on measured vertical wind or aircraft vertical acceleration gives an objective and aircraft-independent indication of turbulence severity oretically. The wind-based method proposed by National Center for Atmospheric Research (NCAR) estimates EDR value from vertical wind spectrum within a discrete frequency range and maximum lielihood estimation [1]. A modified method is to employ a running-mean standard deviation calculation of bandwidth-filtered vertical wind with different angle of attac calibration and sampling frequencies [11]. Although integrated in an aircraft meteorological data relay system [1] wind-based EDR estimation still needs to be improved because some parameters including band-pass frequency and bias-correction term are empirically chosen which maes estimated EDR value sensitive to aircraft maneuvers [13]. By recovering turbulence severity from aircraft vertical acceleration acceleration-based EDR estimation is put forward before wind-based method [8]. Similar to VA method vertical acceleration is used as input of EDR estimation. However turbulence severity is obtained from vertical acceleration response. The acceleration response initially described by a polynomial fitting function although lacing adequate accuracy is still applied to EDR estimation nowadays [14]. The plunge and pitch motion of an aircraft providing largest contribution to vertical acceleration are modeled by two linear transfer functions based on small perturbation approximation. In addition with wind effects as only excitation of interest a quasi-static aerodynamic model and controls-fixed assumptions are made in derivation of linear model. To furr explore accuracy of acceleration response four linear models were developed in ref [15]. These models were tested respectively and compared with recorded acceleration data on Boeing757 aircraft. In ref [16] real-time turbulence effects during aerial refueling were studied and turbulence parameters were identified based on Dryden turbulence model assumption. By feeding identified turbulence into a flight dynamics model simulation results were compared with recorded data. In ref [17] to study riding comfort in turbulence acceleration responses in different seating locations were explored. The above research studies show that re are major differences between measured and oretical acceleration response. The linear model or fitting function cannot accurately describe acceleration response in turbulence. Furrmore restricted by turbulence-only excitation se methods are unable to distinguish acceleration change induced by aircraft maneuvers or turbulence; thus EDR values obtained are easily be affected by control surface deflections [18]. In potential flow ory panel method is applied to obtain aerodynamic performance of lifting body. Different inds of singularities can be assigned to panel element and strength of each singularity is determined according to its boundary condition. After that pressure distribution and aerodynamic force are furr computed. As an effective way of aerodynamic analysis and optimization Vortex Lattice Method (VLM) is used to obtain acceleration response with turbulence effects. By assigning horseshoe vortices planar VLM was initially applied to lifting surface computation. In addition non-planar vortex rings are assigned on mean camber surface

3 Atmosphere of to get solutions with better accuracy [19]. Compared with non-real time Computational Fluid Dynamics (CFD) method VLM is a ind of rapid and medium-accurate algorithm that has been widely used in aero-elastic analysis [] wind turbine design [1] lift and induced drag computation [] lifting body design optimization [3] and gust wind response analysis [4]. In special applications such as aerodynamic analysis in aerial refueling formation flight [5 7] and innovative flap design [8] VLM possesses advantage of rapid trial and error. When it comes to aero-elasticity and gust wind analysis Unsteady Vortex Lattice Method (UVLM) is developed with influence of aircraft maneuvers and wae on aerodynamics considered [3]. Furrmore with no-sideslip assumption VLM shows satisfactory accuracy in developing wind shear coefficients by computing additional aerodynamic effects in wind shear [93]. However steady VLM is not suitable for computing aerodynamic force in high-frequency turbulence. If 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 Quic Access Recorder (QAR) flight data of civil aviation aircraft 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 furr acceleration-based EDR estimation is explored in this paper. According to recorded acceleration acceleration response in turbulent flight will be deeply analyzed in both time and frequency domain. Furrmore compared with wind-based EDR estimation performance and accuracy of new EDR estimation algorithm will be analyzed.. Turbulence Severity Estimation Based on In-Situ Flight Data.1. Derived Vertical Wind From QAR Flight Data To estimate EDR value vertical turbulence component must be first obtained as 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 8 aircraft collects and records a great number of flight parameters in form of data frames via avionic bus [31]. With support of Chinese domestic airlines QAR flight data of B737 8 aircraft are gared and used as algorithm input. Moreover recorded acceleration data are used for verification. Table 1 compares parameter requirement of both wind-based and acceleration-based EDR estimation. Table 1. Parameter requirement for both Eddy Dissipation Rate (EDR) algorithms. Flight Parameter Symbol Wind-Based Acceleration-Based True airspeed V T Angle of attac α Pitching angle θ Roll angle φ Ground speed V e Pitching angular rate q Rolling angular rate p Yawing angular rate r The vertical turbulence component W z not recorded in QAR is mandatory input for wind-based and acceleration-based EDR estimation. Spatial turbulence is determined as difference between ground speed and airspeed W e = V e C e b Cb av T (1) where W e = [W x W y W z ] T V e stands for ground speed vector C e b represents transfer matrix from body-axis to earth-axis and C b a shows transition from wind-axis to body-axis. In straight

4 Atmosphere of 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 W x = V e 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 attac recorded by QAR needs to be converted into angle of attac in body axis before being used to obtain vertical wind. In solving Equation () body-axis angle of attac α b is obtained from measured angle of attac α through: α b = a + a 1 α. (3) The calibration constants a and a 1 are estimated in a preprocessing step using a least-squares linear fit of θ = a + a 1 α. In straight and level flight α b should be approximately equal to pitching angle θ under nominal smooth conditions. The or parameters in Equation () are obtained directly from QAR. A major advantage of wind-based EDR estimation is that it requires fewer parameters than acceleration-based method. When it comes to acceleration-based EDR estimation proposed in this paper angular rate [p q r] T is required for obtaining acceleration response. Since flight parameters used in Equation () are recorded with different sampling rates from 4 to 16 Hz turbulence components are estimated at 4Hz with lowest sampling rate... The von Karman Turbulence Model The Dryden and von Karman turbulence models are commonly used to describe stochastic behavior of well-developed and small-scale high-altitude turbulence. However von Karman model can describe characteristics of turbulence with better accuracy. On one hand energy spectrum of von Karman model is built according to large amounts of statistics data and is consistent with energy spectrum of Kolmogorov model. The spectrum function of von Karman model has a roll-off rate of 5/3 in high-frequency section. However as an approximation of von Karman model roll-off rate of Dryden model is. On or hand when it comes to EDR estimation since most of energy responsible for aircraft bumpiness is in 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 inertial sub-range. The von Karman model representing both inertial sub-range and larger scales beyond it has been widely used in both aerodynamics and meteorological communities [3]. As a logical choice in describing turbulence characteristics 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]: [ ] B NN (r) = σ /3 ( ) r 1/3 [ ( ) r W z K Γ(1/3) L 1/3 1 ( ( )] r r K L L) /3 (4) L where L is length scale and σ W z is variance of vertical turbulence component W z. Γ is Gamma function and K v is modified Bessel function. Equation (4) is used to obtain oretical power spectrum of turbulence in wind-based EDR estimation. According to von Karman and Kolmogorov energy spectrum ory energy spectrum of turbulence only depends on eddy dissipation rate ε. The EDR value is generally defined as ε 1/3 which is cubic root of EDR. The relationship between σ W z and ε is given as follows [8]: σ W z = ε /3 Aπ 1/ L /3 9 Γ(1/3) 55 Γ(5/6) (5)

5 Atmosphere of where A has been estimated to be 1.6. For typical flight conditions of civil aviation aircraft frequency band within which majority of vertical acceleration response occurs corresponds to inertial sub-range of atmospheric turbulence. The velocity spectrum of vertical turbulence φ Wz (ω) can be approximated by: φ Wz (ω) =.7V /3 ε /3 ω 5/3 (6) where Equation (6) forms basis of acceleration-based EDR estimation..3. Vertical Wind-Based EDR Estimation With estimated vertical turbulence component as input a frequency-domain and single-parameter maximum lielihood algorithm are used in vertical-wind based EDR estimation. Firstly a Tuey Hanning window is employed to reduce spectrum leaage due to finite length data of vertical turbulence component W z. As a well-performed window function main-side lobe ratio of Tuey Hanning window is higher and side lobe converges quicly. The m-point Tuey Hanning window is formulated as: τ = ( ( 1 1 cos π M+1)) M 1 M < m M 1 (1 cos( (m 1)π M+1 )) m M m 1 where M = f loor(.1m.). In more traditional formulation of Tuey Hanning window this corresponds to a taper factor of.. The power normalized Tuey Hanning window is: (7) τ = τ m 1 1 m τ j j =. (8) Accordingly windowed vertical wind data series is: W w = τ W. (9) Secondly velocity spectrum of windowed vertical wind series is estimated by Fourier transform: ˆφ W z = m 1 W w f s m zj e iπj/m (1) j = where i is complex imaginary number f s is sampling frequency of vertical wind m = 1 f s corresponds to 1-second wind data and =... 5 f s. By dividing power spectrum of empirical vertical wind by oretical power spectrum over a certain range of frequencies estimated EDR value ˆε 1/3 is obtained by: ˆε 1/3 = γ h l + 1 h = l ˆφ Wz φ model 1/ (11) where γ is a bias-correction term and l and h are lower and upper index bounds corresponding to frequencies f l and f h over which average is taen. In practice lower cutoff frequency f l must be greater than aircraft phugoid response frequency while high-end cutoff frequency f h is made to lie below aero-elastic response mode of target aircraft. When it comes to a turbulence

6 Atmosphere 11 x FOR PEER REVIEW 6 of 1 1 m T = j j m ττ +. (1) j= Atmosphere of After that φ model is average periodogram of windowed von Karman autocorrelation 1 3 function (shown in model spectrum φ model Equation (4)) with unit ε : biased autocorrelation function T for normalized Tuey Hanning m 1 m 1 window τ must model be computed 1 firstly jv by: πij jv πij φ = TmBNN( )exp( ) = Re TjBNN( )exp( ) BNN() fs j ( m 1) fs m fs j fs m (13) = =. T = 1 m 1 τ j τ m j+. (1).4.A New Vertical Acceleration-Based EDR Estimation Algorithm j = To estimate After that φ model aircraft-independent EDR value acceleration response excited by turbulence is average periodogram of windowed von Karman autocorrelation is firstly computed by integrating non-planar function (shown in Equation (4)) with unit ε 1/3 UVLM. After that a new acceleration-based EDR : estimation algorithm is put forward. φ model = 1 m 1 f s T m B NN ( jv f s ) exp( πij m ) = m Non-Planar Unsteady Vortex Lattice Method f s Re T j B NN ( jv f s ) exp( πij m ) B NN (). (13) j = (m 1) j = Inviscid irrotational and incompressible flow is governed by Laplace s equation which is derived.4. A New from Vertical continuity Acceleration-Based equation: EDR Estimation Algorithm To estimate aircraft-independent EDR value φ= acceleration response excited by turbulence (14) is firstly computed by integrating non-planar UVLM. After that a new acceleration-based EDR where estimation φ is algorithm velocity is put potential. forward. Based on linear Laplace s equation elementary solutions can be superposed to solve flow around complex geometries. Vortex rings composed of several vortex segments.4.1. Non-Planar are placed Unsteady on mean Vortexcamber Latticesurface Methodto model lifting body. As shown in Figure 1 in order Inviscid to achieve irrotational more accurate and incompressible results vortex flowsections is governed of AE by and Laplace s BD are furr equation divided whichby is derived section AF fromand FE continuity BC and equation: CD respectively. Taing vortex segment ED as an example induced velocity of ED is obtained with Biot Savart equation φ = [19]: (14) where φ is velocity potential. Based Γ on r1 linear r Laplace s r1 r V = r equation elementary solutions can be = Γ 4π r1 r r r (15) superposed to solve flow around complex geometries. 1 Vortex. rings composed of several vortex segments are placed on mean camber surface to model lifting body. As shown in Figure 1 in order In Equation to achieve (15) V moreis accurate induced results velocity Γ vortexis sections strength of AEof and BD vortex are furr filament divided and r by is section vector AFof and FE vortex BC and filament. CD respectively. r 1 and r are Taing position vortex vectors segment of both ED as ends an example to an arbitrary induced point velocity of ED is obtained with Biot Savart equation [19]: in space. By defining coefficient vector induced velocity of a vortex ring to an arbitrary point in space can be expressed as V = Γ ( rvector 1 r sum of 4π r 1 r r r1 r induced ) velocity induced by six vortex filaments: = Γ. (15) r r 1 V =Γ ( ) =ΓK. (16) AB BC CD DE EF FA Figure 1. Vortex ring and induced velocity. Figure In Equation 1 also (15) illustrates V is induced boundary velocity condition Γ is for strength vortex of ring vortex in which filament bound and r vortex is is vector placed of on vortex 1/4 filament. chord of r 1 and panel r are and position Neumann vectors boundary of bothcondition ends to an of arbitrary no penetration point in is

7 Atmosphere of space. By defining coefficient vector induced velocity of a vortex ring to an arbitrary point in space can be expressed as vector sum of induced velocity induced by six vortex filaments: V = Γ( AB + BC + CD + DE + EF + FA ) = ΓK. (16) Atmosphere 11 x FOR PEER REVIEW 7 of Figure 1 also illustrates boundary condition for vortex ring in which bound vortex is placed on 1/4 chord of panel and Neumann boundary condition of no penetration is enforced at collocation point located at 3/4 chord. This is referred to as a fundamental concept enforced at collocation point located at 3/4 chord. This is referred to as fundamental concept for vortex lattice method by which section lift curve slope corresponds exactly to that of thin for vortex lattice method by which section lift curve slope corresponds exactly to that of thin airfoil ory [34]. airfoil ory [34]. In trimmed flight it is vertical turbulence component that leads to fluctuation of In trimmed flight it is vertical turbulence component that leads to fluctuation of longitudinal longitudinal aerodynamic force and furr change of vertical acceleration. Thus target aircraft aerodynamic force and furr change of vertical acceleration. Thus target aircraft is simplified by is simplified by a wing-horizontal tail combination and vortex rings are distributed on ir mean a wing-horizontal tail combination and vortex rings are distributed on ir mean camber surfaces. camber surfaces. As shown in Figure vortex rings are assigned on whole reference area rar As shown in Figure vortex rings are assigned on whole reference area rar than lifting surface than lifting surface only exposed outside fuselage. The wing and horizontal tail are divided only exposed outside fuselage. The wing and horizontal tail are divided by dozens of rectangular by dozens of rectangular vortex rings respectively. Moreover a detailed grid arrangement and vortex rings respectively. Moreover a detailed grid arrangement and furr geometric parameters furr geometric parameters computation are conducted as shown in supplementary materials. computation are conducted as shown in supplementary materials. Figure S1 shows grid Figure S1 shows grid distribution on wing and horizontal tail. Equation S1 gives definition of distribution on wing and horizontal tail. Equation S1 gives definition of each grid number. Equation each grid number. Equation S and S3 give specific coordinate of each grid on wing and S and S3 give specific coordinate of each grid on wing and horizontal tail respectively. Since horizontal tail respectively. Since circulation distribution changes rapidly at wingtip wing circulation distribution changes rapidly at wingtip wing root and leading edge neighboring root and leading edge neighboring lattices are re-designed according to semi-circle division lattices are re-designed according to semi-circle division method [35]. In addition boundaries method [35]. In addition boundaries of each vortex ring approximate geometric edge of of each vortex ring approximate geometric edge of control surfaces since elevators and spoilers control surfaces since elevators and spoilers can be deflected by pilot s manipulation or can be deflected by pilot s manipulation or Automatic Flight Control System (AFCS). Automatic Flight Control System (AFCS). Figure. Vortex rings on mean camber surface of wing tail combination. In unsteady vortex lattice method a time-marching scheme is applied to perform unsteady simulation. At each timestep stept t a row of newly vortex rings is released into wae with a row of newly generated vortex rings is released into wae with same same circulation circulation as as trailing-edge trailing-edge panel from panel which from it which was shed it was at shed previous at previous time step. time As step. shown As in shown Figure in Figure 3 3 circulation circulation Γ of Γ of newly newly shed shed wae wae panel panel stays stays unchanged unchanged for for remainder remainder of of simulation simulation as as dictated dictated by by Helmholtz s Helmholtz s orem orem [19]. [19]. The The wing tail wing tail combination combination releases releases force-free force-free waes waes with with no no aerodynamic aerodynamic loads. loads.

8 Atmosphere of Atmosphere 11 x FOR PEER REVIEW 8 of Figure 3. Non-planar unsteady vortex ring modeling. Although turbulence field is is rotational scale scale of of turbulence is large is large enough enough compared to to aircraft. aircraft. So it So isit feasible is feasible to analyze to analyze boundary boundary condition condition and furr and furr compute compute aerodynamic aerodynamic force by force UVLM by UVLM according according to turbulence to turbulence variation variation at each at time each step. time Under step. Under effects of effects airfoil-attached of attached wae vortices vortices wae andvortices local velocity and local Neumann velocity boundary Neumann condition boundary at any collocation condition point at any is: vortices collocation point is: N N j T T ( K j i Γ j(t )) n( i + K[V i Γ Wx j( ti(t )) n) i + V[ Wy W ( ) ( ) ( )] +[ ( ) ( ) ( )] = x ii(t t V) V Wy i Wz t i(t VWz )] i t T n i + i V[V x i xi t (t V ) y i tv yi V(t z i t) V zi n(t i )] T n i =. (17) (17) j= 1 j = 1. By By furr furr expanding expanding above above equation equation a linear linear algebraic algebraic equation equation is is obtained. obtained. The The circulation circulation of Γ each vortex ring is obtained by Γ 1( t ) ( t)... ΓN( t) 1 (t ) Γ (t )... Γ N (t. It ). should It should be noted be noted that that in in turbulent flight flight in in addition to far-field free flow V V ( t ) [ ( ( ) ( )] addition to far-field free flow (t T ) = [V= x V (t x t ) V y y t(tv ) z tv z (t and )] T and aircraft angular motion aircraft angular motion [p(t [ pt ( ) ) ( ) ( )] T qt q(t rt ) r(t )] T instantaneous turbulenceww(t ( t ) [ ( ) ( ) ( )] instantaneous turbulence T = ) W= x t[w W x (t y t ) WW z t y (t ) W z (t )] T has major effects on local velocity. The local velocity of an arbitrary point at any has major effects on local velocity. The local velocity an point at any time t time t is described by Equation (18) with [x y z] T coordinate of arbitrary point: is described by Equation (18) with [ x yz ] T coordinate of arbitrary point: V x (t ) V x (t ) cos α(t ) cos β(t ) q(t )z + r(t )y W x (t ) V y (t ) Vx( t) V x( t )cos α( t)cos β( t) q( t) z+ t ) y Wx( t ) = V y (t ) sin β(t ) y( t) = V y( t )sin β( t) + + r(t )x + p(t )z r ( t) x+ t) z Wy( t ) V z (t ). + W y (t ). (18) V z (t ) sin α(t ) cos β(t ) p(t )y + q(t )x W z (t ) (18) V ( t ) V ( t )sin α( t )cos β( t ) p( t ) y+ q( t ) x W ( t ) z z z After obtaining circulation of each vortex ring aerodynamic force is furr computed After obtaining circulation of each vortex ring aerodynamic force is furr computed according to Kutta Jouowsi orem. As shown in Figure 4 any vortex and its adjacent vortices according to Kutta Jouowsi orem. As shown in Figure 4 any vortex and its adjacent vortices need to be superposed on lattice panel to get actual circulation distribution. Taing right half need to be superposed on lattice panel to get actual circulation distribution. Taing right half wing as an example aerodynamic force at collocation point is: wing as an example aerodynamic force at collocation point is: F i (t ) = Fρ(Γ ( t (nrw ) = n ρ CW ( Γ)(t ) Γ ( (nrw 1 t ) Γ n CW )(t ))(V( t bi ))((t V ) ( + t [V ) + xi [ V(t ) ( tv) yi V(t () t V) zi V(t ( t)] T )] + + [V[ V Wxi(t( t ) V Wyi(t ( t ) ) VV Wzi(t ( t )] )] ) T ) r r AB T T i ( nrw ncw ) ( nrw 1 ncw ) b i x i y i z i W x i W y i W z i AB T T + ρ ( Γ ( n )( ) ( 1 )( ))( ( ) [ ( ) ( ) ( )] RW n t CW Γ nrw n t CW V b i t + Vx i t V y i t Vz i t + + [ VW ( ) ( ) ( )] ) x i t VW y i t VW z i t rbc T + ρ ( Γ ( n )( ) ( 1 )( RW n t CW Γ nrw n t CW + x y z ))( V ( t ) + [ V ( t ) V ( t ) V ( t )] + [ V ( t ) V ( t ) V ( t T )] ) rcd b i x i y i z i W i W i W i +ρ(γ (nrw n CW )(t ) Γ (nrw n CW+1 )(t ))(V bi (t ) + [V xi (t ) V yi (t ) V zi (t )] T + [V Wxi(t ) V Wyi(t ) V Wzi(t )] T ) r BC (19) +ρ(γ (nrw n CW )(t ) Γ (nrw n CW+1 )(t ))(V bi (t ) + [V xi (t ) V yi (t ) V zi (t )] T + [V Wxi(t ) V Wyi(t ) V Wzi(t )] T ) r CD V where ρ is ρ air is air density density and Vand bi bi is is total induced total induced velocity velocity of airfoil-attached of airfoil-attached vortex at vortex collocation point. collocation At each point. time At each step time t step aerodynamic t aerodynamic force withforce without or without turbulence turbulence effects are effects solved are respectively solved respectively which are which recorded are recorded as F i andas F i F. The and total i F. aerodynamic The total aerodynamic force F = force [F x F i y F F= z ][ T Fcontains ] T x Fy Fz sum of all vortex rings on wing and horizontal tail. The pitching moment around gravity contains sum of all vortex rings on wing and horizontal tail. The pitching moment around center is obtained by: gravity center is obtained by: N N N N M M ( t ) = ( x x) F ( t ) + ( z z ) F ( t ) y (t ) = (x cg x i )F zi (t ) + (z i z cg )F xi (t ) () () y cg i z i i cg x i i= 1 i= 1 i = 1 where ( x ) cg zcg represent location of gravity center. i = 1

9 Atmosphere of where (x cg z cg ) represent location of gravity center. Atmosphere 11 x FOR PEER REVIEW 9 of Figure 4. Neighboring rings for a vortex ring. The differential equationof of vertical accelerationand and pitching moment induced by by turbulence are: are: b b az( t) = q( t) VGx( t ) p( t ) VGy( t) + g cos θ ( t)cos φ( t) + Fz( t) m a z (t ) = q(t )V b I yq ( t) = Gx M (t ) p(t )V b y( t) + ( I z Gy I x) (t ) + g cos θ(t ) cos φ(t ) + F z (t )/m (1) p( t) r( t) + I xz( p ( t) r ( t)). (1) θ( t ) where and φ( I t ) y q(t ) = M y (t ) + (I z I x )p(t )r(t ) + I xz (p are pitching and rolling angle at t (t ) r (t )) m is mass of aircraft and Ix Ixz I x I xz where Iθ(t y ) and φ(t ) are pitching and rolling angle at t m is mass of aircraft and I y Ixz I I z xz I z stands for inertia stands matrix. for Ininertia order to matrix. separate In aircraft order to maneuvers separate effects aircraft from maneuvers acceleration effects response from acceleration Equation (1) response needs to be Equation computed (1) twice. needs Firstly to be computed total aerodynamic twice. Firstly force Ftotal is used aerodynamic in Equation force (1) a ( ) F tois get used a complete in Equation acceleration (1) to response get a complete a acceleration response z t z (t ) which is induced by aircraft maneuvers which is and induced external by aircraft turbulence. maneuvers Secondly and xternal aerodynamic turbulence. force without Secondly turbulence aerodynamic effects F force functions without as input turbulence to get acceleration only containing maneuver effects a effects F functions as input to get acceleration z(t only ). The incremental acceleration only containing maneuver effects a z ( tinduced ). The by turbulence is obtained by a z (t ) a z(t ). In addition vertical acceleration at any location on incremental acceleration only induced by turbulence is obtained by az ( t) az( t) longitudinal plane is obtained by:. In addition vertical acceleration at any location on longitudinal plane is obtained by: a z(t a ) ( t = ) = a z a (t ( t ) ) + ( x (x x x ) cg ) + ( + z (z z ) z q cg ( ) t q ). (t ). () () z z cg cg It It shouldbe notedthat that agap gapof of longitudinal aerodynamic performance exists inevitably between real real aircraft and wing tail combination. However this this paper maes a compromise between real-time performance and and algorithm algorithm accuracy. accuracy. On On one one hand hand compared compared with with fuselage fuselage aircraft aircraft wing and wing horizontal and horizontal tail play atail major play role a major in pitching role in andpitching plunge motion and plunge whilemotion integration while of integration fuselage will of greatly fuselage increase will greatly complexity increase of complexity algorithm. of On algorithm. or hand On onlyor acceleration hand only increment acceleration caused increment by vertical caused turbulence by vertical needs turbulence to be considered needs to for be considered EDR estimation. for EDR A satisfactory estimation. A accuracy satisfactory of furr accuracy EDR estimation of furr is EDR ableestimation to be guaranteed is able byto be incremental guaranteed acceleration by incremental response. acceleration response..4.. Acceleration-Based EDR Estimation.4.. In Acceleration-Based turbulent flight EDR spatial Estimation frequency and intensity of turbulence are actually obtained by measuring temporal frequency and aircraft acceleration response. Given frequency response In turbulent flight spatial frequency and intensity of turbulence are actually obtained by function from turbulence input W measuring temporal frequency z (iω) to vertical acceleration a and aircraft acceleration response. z (iω) spectrum of acceleration Given frequency response response is obtained by: function from turbulence input Wz ( iω ) to vertical acceleration az ( iω ) spectrum of acceleration response is obtained by: φ az (ω) = a z (iω) φ W z (iω) Wz (ω). (3) az ( iω) φa ( ω) = φ ( ) z W ω. (3) z W ( iω) z Equation (3) is built up in temporal-frequency domain since measurements are made by aircraft flying through turbulence with airspeed V. According to Taylor s frozen-field hyposis

10 Atmosphere of Equation (3) is built up in temporal-frequency domain since measurements are made by aircraft flying through turbulence with airspeed V. According to Taylor s frozen-field hyposis length scale of vertical turbulence component L can be obtained by L = ω/v. The estimation is conducted in frequency-domain and above equation is integrated between frequency segments ω l and ω f to obtain acceleration response energy: ˆσ a z = ωh ω l φ az (ω)dω = ωh ω l a z (iω) φ W z (iω) Wz (ω)dω. (4) Similar to 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 influence of aircraft maneuvers and high-frequency aero-elastic vibration with ω l = π f l and ω h = π f h. The above equation can be modified by: ˆσ a z H bp (ω l ω h ω) a z (iω) φ W z (iω) Wz (ω)dω. (5) To estimate response energy band-pass filter is transformed into time-domain firstly by Inverse Fourier transform (IFT) as follows: h bp (t) = 1 H π bp (ω l ω h ω)e iωt dω. (6) After that Parseval s orem (assuming local stationary) and convolution orem are adopted to yield approximation [36]: ˆσ a z (t ) 1 t +T τ1 [ h T bp (τ )a z (τ 1 τ )dτ ] dτ 1 (7) t T where a z (t) is acceleration data series obtained by Equation (1). It states that mean-square value of frequency-limited response is obtained by computing a running mean-square value of band-pass filtered acceleration data. A nonzero mean value must be removed in computation of outer integration in Equation (7). The averaging integral T is selected to satisfy local-stationarity assumption and provide enough data values for statistical stability. The time constant τ 1 is chosen such that h bp (τ) for τ > τ 1. By inserting Equation (6) into Equation (7) ˆσ a z =.7V /3 ε /3 H bp (ω l ω h ω) a z (iω) ω 5/3 dω =.7V /3 ε /3 I(ω W z (iω) l ω h t). (8) If I(ω l ω h t) in Equation (8) is nown EDR value ε 1/3 will be estimated. To estimate acceleration response Fourier transform is performed on vertical acceleration and vertical wind sequence in period T. After that it is converted from time-domain to frequency-domain to obtain a z (iω) and W z (iω) respectively. With new time series a z(iω) autocorrelation function R W z (iω) is obtained by IFT. According to Wiener Khinchin orem R() is taen as estimation results of I(ω l ω h t ). As a result EDR value is estimated by: ε 1/3 (t ) = ˆσ az (t ).7V /3 (t )I(ω l ω h t ) (9) To summarize an integrated algorithm flow for new acceleration-based EDR estimation is shown in Figure 5. The grid on wing tail combination is generated beforehand. In algorithm loop measured flight parameters are used to obtain turbulence components by Equation () and furr local velocity is obtained by Equation (18). The parameters shown in Table 1 function as data input.

11 Atmosphere of Atmosphere 11 x FOR PEER REVIEW 11 of After flag that total aerodynamic aerodynamic force force is acquired F and by non-planar force without UVLM. turbulence In each cycle effects byf toggling are obtained flag respectively. total aerodynamic The vertical forceacceleration F and force increment withoutinduced turbulence by effects turbulence F areis obtained respectively. by solving The vertical acceleration increment induced by turbulence is obtained by solving Equation ˆ σ a ( t ) Equation (1) twice. Next time series of acceleration data are used to estimate z by (1) Equation twice. Next (7). In addition time series after offast acceleration Fourier Transform data are used (FFT) to estimate acceleration ˆσ az (t ) series by Equation in frequency-domain (7). In addition after Fast Fourier Transform (FFT) acceleration I series in frequency-domain are used to obtain are used to obtain transfer transfer function I(ω l ω h t ) function ( ω ) l ωh t toger with frequency-domain vertical toger with frequency-domain vertical turbulence component. turbulence component. At end of each iteration estimated EDR value is output by solving At end of each iteration estimated EDR value is output by solving Equation (9). Equation (9) Experiments and Discussion Figure Figure Algorithm flow flow of of new new acceleration-based EDR EDR estimation Experiments on on Non-Planar UVLM Analysis of Grid Convergence Analysis of Grid Convergence Before acceleration response analysis performance of non-planar UVLM was tested. UVLM Before acceleration response analysis performance of non-planar UVLM was tested. requires a sufficiently refined grid to achieve results that are grid-independent. The grid resolution of UVLM requires a sufficiently refined grid to achieve results that are grid-independent. The grid UVLM was first studied to improve computing efficiency on premise of accuracy assurance. resolution of UVLM was first studied to improve computing efficiency on premise of accuracy As far as Boeing aircraft is concerned geometry and airfoil parameters of wing assurance. and horizontal tail are listed in ref [3738]. A non-planar vortex lattice model was built in which As far as Boeing aircraft is concerned geometry and airfoil parameters of wing different surface grids of increasing mesh density were generated. Vortex rings are assigned on and horizontal tail are listed in ref [3738]. A non-planar vortex lattice model was built in which mean camber surface. Each grid has a constant spacing in both chord-wise and span-wise directions. different surface grids of increasing mesh density were generated. Vortex rings are assigned on The changes of lift coefficient drag coefficient and pitching moment coefficient are presented as mean camber surface. Each grid has a constant spacing in both chord-wise and span-wise directions. a function of grid refinement level. As a result satisfactory accuracy is obtained by adopting a The changes of lift coefficient drag coefficient and pitching moment coefficient are presented as a function of grid refinement level. As a result satisfactory accuracy is obtained by adopting a

12 Atmosphere x 147 FOR PEER REVIEW 1 1of of grid with 4 columns rows on wing and 1 columns 6 rows on horizontal tail. Furr experiments grid with 4are columns based on rows optimized on wing grid and level. 1 columns 6 rows on horizontal tail. Furr experiments are based on optimized grid level Computation of Pitching Moment Coefficient Computation of Pitching Moment Coefficient There exists a gap of longitudinal aerodynamic characteristics between whole aircraft and wing tail There exists combination. a gap of longitudinal However since aerodynamic pitching characteristics moment change between is mainly whole caused aircraft by and deflection wing tail of elevator combination. or stabilizer However computing since results pitchingbased moment on change wing tail is mainly combination causedshould by be deflection comparable of elevator with auntic or stabilizer aerodynamic computing data. results The modeling based ondata of wing tail B737 8 combination aircraft provides should longitudinal be comparable pitching with auntic moment changes aerodynamic caused data. by elevator The modeling and stabilizer data ofdeflections B737 8 aircraft [3738]. provides Figure 6 longitudinal pitching moment changes caused by elevator and stabilizer deflections [3738]. Figure c 6 shows results between UVLM results and modeling data in three flight conditions. shows results between UVLM results and modeling data in three flight conditions. c mδstab stands m δ stab c stands for for pitching pitching moment moment change change due todue stabilizer to stabilizer deflection deflection and c mδe and m represents δ e represents moment moment change change due to elevator due elevator deflection. deflection..4. Modeling data UVLM stab (o ) (a) M=. =3 o H= stab (o ) (b) M=.5 =3 o H=5m stab (o ) (c) M=.76 = o H=1m 1.5 Modeling data UVLM e (o ) (d) M=. =3 o H= e (o ) (e) M=.5 =3 o H=5m e (o ) (f) M=.76 = o H=1m Figure6. 6. Comparison of of computation computation and and experimental experimental pitching pitching moment moment coefficients. coefficients. (a) M = (a). M α = 3. H α = ; 3 (b) H M = = ;.5 (b) α M = = 3.5 H = α 5 = 3 m; H (c) = 5 M =.76 m; (c) α M = =.76 H = 1; α = (d) H M = = 1;. α (d) = 3 M H = =. ; (e) α = M 3=.5 H = α ; = (e) 3 M H = 5;.5 α (f) = 3M = H.76 = 5; α = (f) M H = α = H = 1. Compared with modeling data mean square error (MSE) of c Compared with modeling data mean square error (MSE) of c mδstab is.45 and MSE m δ stab is.45 and MSE of c mδe is.55. It shows that pitching moment change caused by deflection of stabilizer or of elevator cm δ e is approximates.55. It shows that modeling pitching data under moment change above flight caused conditions. by deflection Since of stabilizer geometrical or elevator boundary approximates of control surface modeling approaches data under boundary above offlight each vortex conditions. ring non-planar Since geometrical UVLM can boundary be used toof accurately control describe surface approaches pitching moment boundary change. of each vortex ring non-planar UVLM can be used Theto experiments accurately describe prove that pitching based onmoment wing tail change. combination non-planar UVLM has satisfactory The experiments computing prove accuracy that based whichon provides wing tail a goodcombination basis for furr non-planar estimationuvlm of aircraft has satisfactory accelerationcomputing and turbulence accuracy severity. which Furrmore provides a good accurate basis computation for furr of estimation pitching of moment aircraft acceleration coefficient isand significant turbulence to separate severity. aircraft Furrmore maneuvers accurate effects oncomputation EDR estimation. of pitching moment coefficient is significant to separate aircraft maneuvers effects on EDR estimation. 3.. Acceleration Response Analysis 3.. Acceleration Response Analysis Since accurate acceleration response forms basis for EDR estimation acceleration response by non-planar Since accurate UVLMacceleration will be analyzed response andforms furr compared basis for with EDR measured estimation acceleration acceleration in QAR. response Supportedby bynon-planar Chinese airlines UVLM QAR will flight be analyzed data including and furr 4 flight compared segments with have measured been gared acceleration from in B737 8 QAR. Supported aircraft onby Chinese same scheduled airlines QAR air route. flight Five-minute data including cruising 4 flight data segments in each have segment been gared from B737 8 aircraft on same scheduled air route. Five-minute cruising flight data in

13 Atmosphere of with 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 comparative study of acceleration response. The turbulence components are computed by Equation () as input of following EDR estimation. It should be noted that re are acceleration data in QAR of B737 8 aircraft while not all aircraft have acceleration record. Atmosphere 11 x FOR PEER REVIEW 13 of The acceleration response computed by non-planar UVLM is compared with linear model and recorded each QAR segment data. with The pea same geographic vertical acceleration area are selected datafor ofedr estimation. six wind fields According (WF1 6) to pilot are shown in Table. reports It is worth six Wind noting Fields that (WF) covering acceleration light moderate data and is recorded severe turbulence at aircraft are selected gravity to conduct center. comparative study of acceleration response. The turbulence components are computed by Equation () as input of following EDR estimation. It should be noted that re are Table. Selected six wind fields. acceleration data in QAR of B737 8 aircraft while not all aircraft have acceleration record. Wind Fields Turbulence Severity Max/Min Acceleration The acceleration response computed by non-planar UVLM is compared with linear model and recorded WF1 QAR data. The pea vertical Light acceleration data of six wind fields g/+.98 (WF1 6) g are shown in Table. It is WF worth noting that acceleration Light data is recorded at +1.3 aircraft g/+.97 gravity g center. WF3 Moderate +1.6 g/+.655 g WF4 Table Moderate. Selected six wind fields g/+.67 g Wind WF5 Fields Turbulence Severe Severity Max/Min +.4 g/+.169 Acceleration g WF1 WF6 Light Severe g/+.175 g/+.98 g WF Light +1.3 g/+.97 g WF3 Moderate +1.6 g/+.655 g Comparison of Acceleration Response WF4 Moderate g/+.67 g Referring to model 4 WF5 in ref [15] a linear Severe model of B aircraft g/+.169 was g built based on WF6 Severe +.14 g/+.175 g modeling data. To better understand gust penetration effect and mechanism underlying gust pitch3..1. ratecomparison term fourof linear Acceleration models Response of increasing gust excitation complexity were developed in ref [15]. Model 1 is simplest model by using a gust point approximation. Model considers gust Referring to model 4 in ref [15] a linear model of B737 8 aircraft was built based on acceleration modeling and gust data. pitch To better rate understand effects in gust pitch penetration dynamics. effect and Furrmore mechanism by underlying accounting for full-body gust pitch penetration rate term four effects linear model 3 models of incorporates increasing gust excitation se gust complexity excitation were signals developed in in plunge dynamics. ref Model 4 [15]. Model 1 incorporating is simplest amodel quasi-steady by using a aerodynamic gust point approximation. model with Model finite considers lag value gust is able to acceleration and gust pitch rate effects in pitch dynamics. Furrmore by accounting for fullbody gust penetration effects model 3 incorporates se gust excitation signals in plunge follow acceleration transients with better accuracy while pea prediction accuracy is sacrificed. Figure 7a shows dynamics. Model 4 derived incorporating vertical wind a quasi-steady and Figure aerodynamic 7b shows model awith short finite sample lag value of recorded is able to vertical acceleration follow transients acceleration overlaidtransients with with computation better accuracy results while from linear pea model prediction andaccuracy non-planar is UVLM in WF3. Upsacrificed. to 166s Figure data 7a shows imply derived aircraft vertical was wind at aand trimmed Figure 7b flight shows condition. a short sample Atof recorded s aircraft vertical acceleration transients overlaid with computation results from linear model and nonplanar UVLM in WF3. Up to 166s data imply aircraft was at a trimmed flight condition. At reacted to initial disturbance reached maximum positive load near s and n started anor local maximum s aircraft positive reacted to load initial at 17s. disturbance After that reached it experienced maximum apositive continuous load near decrease In this short sequence s and n linear started model anor was local momentarily maximum positive out-of-phase load at 17s. onceafter encountering that it experienced sudden turbulence a while UVLM continuous showed decrease. better In this tracing short sequence accuracy. Although linear model not was a conclusive momentarily argument out-of-phase once observation encountering sudden turbulence while UVLM showed better tracing accuracy. Although not a from Figure 7 suggests that both linear model and UVLM are able to replicate acceleration conclusive argument observation from Figure 7 suggests that both linear model and UVLM response inare moderate able to replicate turbulence. acceleration However response tracing moderate accuracy turbulence. of However UVLM is better tracing than that of linear model. accuracy of UVLM is better than that of linear model. W z (m/s) a z (g) a z (g) Figure 7. Vertical wind and acceleration response in WF3. (a) W z ; (b) vertical acceleration.

14 Atmosphere of Atmosphere 11 x FOR PEER REVIEW 14 of Figure 7. Vertical wind and acceleration response in WF3. (a) Wz; (b) vertical acceleration. Figure 8 shows acceleration response of two models for WF6 along with recorded data. This sample Figure 8 covers shows a severe acceleration vertical response turbulence of event two models resulting for WF6 in an along amplitude with recorded swinging from +.14 gdata. to This sample g over covers 1s. Grouping a severe vertical toger turbulence when event responding resulting toin an amplitude turbulence swinging two from responses tended to +.14 trac g to overall g behavior 1s. Grouping of toger recorded when data. responding The sacrifice to turbulence of pea acceleration two responses response accuracytended brought to trac about overall serious behavior deviation of recorded in trace data. The of linear sacrifice model of pea from acceleration recorded response data. On accuracy brought about serious deviation in trace of linear model from recorded data. On contrary UVLM performed higher sensitivity to severe turbulence with sharp gradients resulting contrary UVLM performed higher sensitivity to severe turbulence with sharp gradients resulting in in an over-prediction of pea values. The tracing MSE was reduced from.47 (linear model) to an over-prediction of pea values. The tracing MSE was reduced from.47 (linear model) to.81 (UVLM)..81 (UVLM). Qualitative Qualitative observation observation leads leads to to conclusion that linear linear model model is less is less sensitive sensitive for moderate for and moderate severeand turbulence. severe turbulence. It can be It concluded can be concluded that that recorded recorded acceleration is better is better followed by UVLM followed than by UVLM linear than model. by linear model. Figure Figure 8. Vertical 8. Vertical wind wind and and acceleration response in WF6. (a) (a) Wz; W(b) vertical acceleration. z ; (b) vertical acceleration. The result in Figure 8 is not unexpected that higher-order feature in recorded acceleration data The result in Figure 8 is not unexpected that a higher-order feature in recorded acceleration data is reliably replicated by UVLM rar than linear model. There are significant differences between is reliably replicated by UVLM rar than linear model. There are significant differences between two models. Firstly severe turbulence could easily change governing aerodynamic twocoefficients models. Firstly of severe linear model turbulence from trim couldcondition. easily change In contrast governing UVLM is a aerodynamic ind of nonlinear coefficients of linear computation model from instead trimof condition. linear approximation. In contrast Secondly UVLM is a ind recorded of nonlinear acceleration computation response is instead of linearactually approximation. excited by external Secondly turbulence recorded and control acceleration surface deflections responseby ismanual actually manipulation excited byor external turbulence AFCS. and The control recorded surface control deflections function by manual as input manipulation in UVLM while oraircraft AFCS. maneuvers The recorded are not control considered in linear model. deflections function as input in UVLM while aircraft maneuvers are not considered in linear model. For a more comprehensive comparison MSE and absolute pea acceleration in WF1 WF6 Forare a more presented comprehensive graphically in Figure comparison 9. The linear MSE model and does absolute a good prediction pea acceleration in light turbulence in WF1 WF6 are presented while suffering graphically accuracy in in Figure moderate 9. The and severe linearturbulence. model does The agreater good prediction turbulence in intensity light turbulence while suffering bigger accuracy MSE will in be. moderate However UVLM and severe is capable turbulence. of predicting The greater acceleration turbulence response with intensity bigger better accuracy MSE will than be. However linear model. UVLM The satisfactory capable prediction of predicting of pea acceleration acceleration holds promise response with for furr EDR estimation. better accuracy than linear model. The satisfactory prediction of pea acceleration holds promise for furr Atmosphere EDR estimation. 11 x FOR PEER REVIEW 15 of Figure 9. Acceleration response comparison of linear model and UVLM. Figure 9. Acceleration response comparison of linear model and UVLM Spectrum Analysis at Different Locations Three typical locations are chosen to investigate acceleration response at different locations on aircraft. Location A is at gravity center and Location B is at cocpit position which is about 15 feet from point A. Location C is at tailgate of aircraft about 6 feet away from gravity center. One of prominent features of UVLM is that acceleration at different locations of

15 Atmosphere Figure 9. Acceleration response comparison of linear model and UVLM. 15 of 3... Spectrum Analysis at Different Locations 3... Spectrum Analysis at Different Locations Three typical locations are chosen to investigate acceleration response at different locations on aircraft. Three typical Location locations A is aregravity chosen center to investigate and Location acceleration B is at response cocpit position at different which locations is about on 15 aircraft. feet from Location point AA. is Location gravity C is center at tailgate and Location of aircraft B is atabout cocpit 6 feet position away from which isgravity about 15 center. feet from One point of A. prominent Location Cfeatures is of tailgate UVLM of is aircraft that about acceleration 6 feet away at different from locations gravity center. of fuselage One of can prominent be obtained. features Since of it UVLM is nearly is that impossible acceleration to find at different differences locations from oftime-domain fuselage acceleration can be obtained. series Since it is spectral nearly impossible density estimation to find is differences used to from obtain time-domain frequency-domain acceleration acceleration series spectral spectrum density at different estimation locations. is used to The obtain acceleration frequency-domain spectra of three acceleration locations spectrum in WF1 at WF6 different are locations. shown in Figure The acceleration 1. To obtain spectra of spectrum threedensities locationsof in WF1 WF6 vertical acceleration are shown inseries Figure 1. To Welch obtain spectral spectrum density densities estimation of with vertical Hanning acceleration window series function Welch is used spectral [35]. density The data estimation length with in time-domain Hanning window for each function curve is is56 used points [35]. The which datarepresents length in a length time-domain of 64 s acceleration for each curve data. is With 56 points sampled which acceleration representsdata a length with of 4Hz 64 s acceleration window length data. is With 1 fs = 4 sampled and acceleration overlapping data is with set by 4Hz 5%. window length is 1 fs = 4 and overlapping is set by 5%. PSD((m/s ) /Hz) PSD((m/s ) /Hz) PSD((m/s ) /Hz) PSD((m/s ) /Hz) PSD((m/s ) /Hz) PSD((m/s ) /Hz) Figure 1. Acceleration spectrum at at different locations. (a) (f) Wind Fields (WF)1 6. It can be found that bumpiness at tail of fuselage is larger than that at gravity center It can be found that bumpiness at tail of fuselage is larger than that at gravity while bumpiness at nose is least. It is consistent with analysis in ref [17]. In moderate center while bumpiness at nose is least. It is consistent with analysis in ref [17]. In turbulence re is little difference in magnitudes of vertical acceleration between locations A B and C. moderate turbulence re is little difference in magnitudes of vertical acceleration between locations However with increase of turbulence severity difference of acceleration magnitude increases. As a result UVLM-based acceleration response algorithm offers an effective way for estimating bumpiness at different locations of fuselage EDR Estimation Comparison After that 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 new acceleration response algorithm a furr 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 window to provide stable computational statistics. For cruising aircraft a nominal 1-minute reporting interval is generally used. For each of 5-minute flight time in 4 flight segments re are totally EDR values. Within reporting window a 1-second sub-window with 1/ overlap is selected for EDR estimation in both wind-based and acceleration-based method. Furrmore median and

16 Atmosphere of 9th percentage ε 1/3 values are computed respectively in 1-minute period. In order to mitigate potential erroneous EDR output due to QAR data quality issues average and pea EDR value are replaced by median and 9th percentage value. If turbulence is relatively continuous median and 9th percentage ε 1/3 values should be close to each or in magnitude. Whereas for a discrete sudden bumpiness 9th percentage value will typically be larger than median value. In addition band-pass frequencies ω l and ω h are typically chosen as.1 and 1. Hz respectively for civil aviation aircraft Application to EDR Estimation Figure 11 shows EDR estimation results in WF6 along with vertical turbulence component control surface deflections and aerodynamic force change. Taing QAR data in WF6 as an example vertical turbulence component obtained by Equation () is shown in Figure 11a. From 9s to 1s aircraft experienced a sudden downdraft. Meanwhile deflection of elevators and spoilers by pilot manipulation or AFCS are shown in Figure 11b. It can be found that low-speed flight spoiler δ sls and high-speed flight spoiler δ shs both experienced sudden deflections. Although deflection of elevator was smaller than that of spoilers it had greater impact on aerodynamics force. Atmosphere In non-planar 11 UVLM x FOR PEER since REVIEW boundaries of vortex ring are close to geometric edge17 of of control surface change of aerodynamic force due to control deflection can be computed accurately. F(N) 1/3 (m /3 s -1 ) W z (m/s) ( o ) Figure 11. EDR calculation under elevator deflection. (a) W z ; (b) control deflections; (c) aerodynamic Figure force; (d) 11. EDR. calculation under elevator deflection. (a) Wz; (b) control deflections; (c) aerodynamic force; (d) EDR. Figure 11c shows aerodynamic force with or without control deflections. With proposed method Comparison in this paper of EDR longitudinal Estimation aerodynamic force induced by aircraft maneuvers or vertical turbulence is able to be effectively distinguished. The effects of vertical turbulence on vertical Based on QAR data of 4 flight segments Figure 1 shows a scatterplot of total acceleration were separated from aerodynamic force F. Therefore change of incremental EDR values for wind-based versus acceleration-based EDR estimation. Figure 1a shows acceleration was only excited by vertical turbulence. It can be found that pilot manipulation may lead scatterplot of one-minute median EDR values whereas Figure 1b shows 9th percentile EDR to an increase in aerodynamic force and deteriorate aircraft bumpiness. values. For each point horizontal coordinate stands for new acceleration-based EDR value and vertical coordinate represents wind-based EDR value. If both estimation methods recover same aircraft-independent EDR value estimated EDR value points should distribute neighboring to one-to-one line which is clearly identified by y = x in Figure 1. However a small number of points are clustered above one-to-one line at region of small estimated EDR values (lower-left corners). These points are due to incomplete removal of all maneuvers-induced accelerations in wind-based method. The band-pass filter with empirically selected band-width

17 Atmosphere of 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 dynamic process of EDR estimation 9th percentage EDR value is given in a sub-window of 1 s with 1/ overlap. Through computation of incremental acceleration influence of control deflection-induced maneuvers is eliminated. It can be concluded that new EDR estimation is able to effectively eliminate adverse effects of aircraft maneuvers Comparison of EDR Estimation Based on QAR data of 4 flight segments Figure 1 shows a scatterplot of total EDR values for wind-based versus acceleration-based EDR estimation. Figure 1a shows scatterplot of one-minute median EDR values whereas Figure 1b shows 9th percentile EDR values. For each point horizontal coordinate stands for new acceleration-based EDR value and vertical coordinate represents wind-based EDR value. If both estimation methods recover same aircraft-independent EDR value estimated EDR value points should distribute neighboring to one-to-one line which is clearly identified by y = x in Figure 1. However a small number of points are clustered above one-to-one line at region of small estimated EDR values (lower-left corners). These points are due to incomplete removal of all maneuvers-induced accelerations in wind-based method. The band-pass filter with empirically selected band-width parameters f h could not eliminate maneuvers effects completely leading wind-based EDR estimates to be more sensitive to control deflections. A first-order function fitting was made over data points with fitting equation y = 1.719x and correlation coefficient It means that wind-based EDR value is statistically higher than acceleration-based EDR value because of aircraft maneuvers contamination. For 9th percentile values fitting function is y = x +.85 and correlation coefficient is It means that in wind-based estimation discrete sudden bumpiness Atmosphere maes 11 x FOR 9th PEER percentage REVIEW EDR value much higher than acceleration-based EDR 18 value. of Figure 1. Comparison of EDR estimation. (a) median EDR value; (b) 9th percentile EDR value. With increase above of non-planar turbulenceuvlm intensity configuration number of discrete CPU time sudden of bumpiness vertical acceleration increases. The computation variance of which EDR value mainly bycontributes wind-based to method complexity increases of algorithm and number costs 1.96s of outliers with an increases. i7 97 Benefitted processor under by accurate a visual vertical studio 17 acceleration development response environment. acceleration-based The whole CPU EDR time estimation of a one-time is still complete satisfactory. EDR estimation is.47s. As a result algorithm is capable of outputing EDR value in real With time. After above being non-planar furr UVLM optimized configuration and embedded CPU into time ofairborne vertical acceleration system computation in situ EDR which estimation mainly is realized. contributes to complexity of algorithm costs 1.96s with an i7 97 processor under a visual studio 17 development environment. The whole CPU time of a one-time complete EDR estimation 4. Conclusions is.47s. As a result algorithm is capable of outputing EDR value in real time. 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 non-planar Unsteady Vortex Lattice Method (UVLM) is proposed to obtain aerodynamic force change in turbulent flight.

18 Atmosphere of After being furr optimized and embedded into airborne system in situ EDR estimation is realized. 4. 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 non-planar Unsteady Vortex Lattice Method (UVLM) is proposed to obtain aerodynamic force change in turbulent flight. To improve computing accuracy of UVLM vortex rings are assigned on mean camber surface and semi-circle division method is adopted to refine lattices neighboring to wingtip wing root and leading edge. Anor major improvement is that lattices neighboring to control surfaces such as spoilers and elevators are divided as close as possible to structural edge which is beneficial for furr EDR estimation.. 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 re are improvements in both response accuracy and tracing performance in moderate and severe bumpiness. The vertical accelerations at different locations of fuselage are computed accurately. Compared with wind-based EDR estimation new acceleration-based EDR algorithm shows better accuracy and stability because adverse influence of aircraft maneuvers on EDR estimation is eliminated effectively. Supplementary Materials: The following are available online at Figure S1: Grid distribution on wing and horizontal tail Equation S1: definition of ith grid number Equaiton S: specific coordinate of each grid on wing Equation S3: specific coordinate of each grid on horizontal tail. Author Contributions: Conceptualization Z.G.; methodology Z.G. and H.W.; software D.W. and Z.X.; 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. All authors have read and agreed to published version of manuscript. Funding: This research was funded by Natural Science Foundation of China grant number [NSFC: U17331 U15331]. Conflicts of Interest: The authors declare no conflict of interest. References 1. Sharman R.D.; Lane T. Aviation Turbulence: Processes Detection Prediction; Springer: Cham Switzerland 16; pp Gao Z.X.; Fu J. Research on Robust LPV modeling and control of aircraft flying through wind disturbance. Chin. J. Aeronaut [CrossRef] 3. Victor N.; Vladimir L.; Eugene N.; Andrei T.; Alesandr B. Measurement of atmospheric turbulence characteristics by ultrasonic anemometers and calibration processes. Atmosphere Boeing Jeppesen Company. JetPlan User s Manual; Boeing Jeppesen Company: Englewood CO USA 14; pp Kim J.H.; Chun H.Y. Statistics and possible sources of aviation turbulence over South Korea. J. Appl. Meteor. Climatol [CrossRef] 6. Sherman D.J. The Australian Implementation of AMDAR/ACARS and Use of Derived Equivalent Gust Velocity as a Turbulence Indicator; AR-4-51; Department of Defence Defence Science and Technology Organisation Aeronautical Research Laboratories: Melbourne Victoria Australia 1985; p Gill P.G. Objective verification of World Area Forecast Centre clear air turbulence forecasts. Meteor. Appl [CrossRef]

19 Atmosphere of 8. Cornman L.B.; Corinne S.M.; Gary C. Real-time estimation of atmospheric turbulence severity from in-situ aircraft measurements. J. Aircr [CrossRef] 9. Michael J.E.; Meymaris G.; Cornman L.B.; Sherry J.; Mulally D. EDR/RMS-G Correlation Analysis. In Proceedings of American Meteorological Society 98th Annual Meeting Austin TX USA 8 January Sharman R.D.; Corman L.B.; Meymaris G. Description and derived climatologies of automated in situ eddy-dissipation rate reports of atmospheric turbulence. J. Appl. Meteor. Climatol [CrossRef] 11. Haverdings H.; Chan P.W. Quic access recorder data analysis for windshear and turbulence studies. J. Aircr [CrossRef] 1. WMO. Aircraft Meteorological Data Relay (AMDAR) Reference Manual; World Meteorological Organization: Geneva Switzerland 3; p Cornman L.B.; Meymaris G.; Limber M. An update on FAA aviation wear research program s in situ turbulence measurement and reporting system. In Proceedings of 11th Conference on Aviation Range and Aerospace Meteorology Hyannis MA USA 4 8 October Huang R.S.; Sun H.B.; Wu C. Estimating eddy dissipation rate with QAR flight big data. Appl. Sci [CrossRef] 15. Billy K.B.; Brett A.N. Aircraft Acceleration prediction due to atmospheric disturbances with flight data validation. J. Aircr Atilla D.; Timothy A.L. Flight data analysis and simulation of wind effects during aerial refueling. J. Aircr Bizinos N.; Redelinghuys C. Tentative study of passenger comfort during formation flight within atmospheric turbulence. J. Aircr [CrossRef] 18. Robinson P.; Buc B.K.; Bowles R.L.; Boyd D.L.B.; Corman L.B. Optimization of NCAR in situ turbulence measurement algorithm. In Proceedings of 38th Aerospace Sciences Meeting and Exhibit Reno NV USA 1 13 January ; AIAA: Reston VA USA; pp Katz J.; Plotin A. Low Speed Aerodynamics nd ed.; Cambridge Press: Cambridge UK 1; pp Joseba M.; Rafael P.; Michael R.G. Applications of unsteady vortex-lattice method in aircraft aeroelasticity and flight dynamics. Prog. Aerosp. Sci Xu B.F.; Wang T.G.; Yuan Y.; Zhao Z.Z.; Liu H.M. A Simplified Free Vortex Wae Model of Wind Turbines for Axial Steady Conditions. Appl. Sci [CrossRef]. Robert J.S.S.; Rafael P.; Joseba M. Induced-drag calculations in unsteady vortex lattice method. AIAA J Oliviu S.G.; Andreea K.; Ruxandra M.B. A new non-linear vortex lattice method: Applications to wing aerodynamic optimizations. Chin. J. Aeronaut Liu Y.; Xie C.C.; Yang C.; Cheng J.L. Gust response analysis and wind tunnel test for a high-aspect ratio wing. Chin. J. Aeronaut [CrossRef] 5. Lue H.P.; Ruben E.P.; Peter W.J. Wae modelling for aerial refueling using aerodynamic adjoints. In Proceedings of 18 AIAA Atmospheric Flight Mechanics Conference Kissimmee FL USA 8 1 January 18; pp Lue H.P.; Ruben E.P.; Peter W.J. Wae modelling with embedded lateral and directional stability sensitivity analysis for aerial refueling or formation flight. In Proceedings of AIAA Scitech 19 Forum San Diego CA USA 7 11 January 19; pp Jan L.; Eie S. Flight dynamics of CS-5 aircraft in formation flight with atmospheric disturbances. In Proceedings of 18 AIAA Aerospace Sciences Meeting Kissimmee FL USA 8 1 January 18; pp Ezra T.; Nhan N. Unsteady aeroservoelastic modeling of flexible wing generic transport aircraft with variable camber continuous trailing edge flap. In Proceedings of 33rd AIAA Applied Aerodynamics Conference Dallas TX USA 6 June 15; pp Dan D.V.; Roland L.B. Effect of spatial wind gradients on airplane aerodynamics. J. Aircr Daniel C.; Gustavo E.C.; Eric T.; Nhan N. Transonic and viscous models for vortex lattice method applied to transport aircraft. In Proceedings of 35th AIAA Applied Aerodynamics Conference Denver CO USA 5 9 June 17; pp

20 Atmosphere of 31. Burns A.A.; Chenovich G.A. Digital Flight Data Acquisition Unit 737 6/ 7/ 8/ 9 Data Frame Interface Control and Requirements Document; Report No.: D6A11 ; The Boeing Company: Seattle WA USA 3; pp Duran I.B.; Schmidli J.; Bhattacharya R. A Budget-Based Turbulence Length Scale Diagnostic. Atmosphere [CrossRef] 33. Frehlich R.; Cornman L.; Sharman R. Simulation of three dimensional turbulent velocity fields. J. Appl. Meteor [CrossRef] 34. Pistolesi E. Vortex-Lattice Utilization; Report No.: SP 45; NASA: Washington DC USA 1976; pp Brenden P.E.; David S.G. A quasi-continuous vortex lattice method for unsteady aerodynamics analysis. In Proceedings of 18 AIAA Aerospace Sciences Meeting Kissimmee FL USA 8 1 January 18; AIAA: Reston VA USA; pp Stoica P.; Moses R. Spectral Analysis of Signals; Prentice Hall: Upper Saddle River NJ USA 4; p Boeing Company. Aerodynamic Data and Flight Control System Description for 737 6/ 7/ 8/ 9 Training Simulator; D611A1 vol1; rev H; The Boeing Company: Seattle WA USA 3; pp Boeing Company. Aerodynamic Data and Flight Control System Description for 737 6/ 7/ 8/ 9 Training Simulator; D611A1 vol3; rev H; The Boeing Company: Seattle WA USA 3; pp Publisher s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. by authors. Licensee MDPI Basel Switzerland. This article is an open access article distributed under terms and conditions of Creative Commons Attribution (CC BY) license (