## 2 Vertical Acceleration and Turbulence Severity Estimation
### Vertical-Acceleration Based EDR Estimation
The time constant \(\tau_{1}\) is chosen such that \(h_{bp}(\tau)\to 0\) for \(\tau>\tau_{1}\). Dryden and von Karman turbulence model are commonly used in describing the stochastic behavior of high-altitude turbulence, while the latter can describe the turbulence characteristics with better accuracy [32]. The energy spectrum of von Karman model is built up according to large amounts of measuring and has a roll-off rate of -5/3 in high frequency section. As an approximation of von Karman model, the roll-off rate of Dryden model is -2 in high frequency section. 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. In addition, faster-moving aircraft will also be responsive to larger scales that are typically outside the inertial sub-range. [PERSON] turbulence model, representing both the inertial sub-range and the larger scales beyond it, has been widely used in both aerodynamics and meteorological field [10]. As a logical choice in describing the turbulence characteristics, von Karman model has been integrated into commonly used EDR estimation algorithm. In addition, according to the [PERSON] and [PERSON] energy spectrum theory, the PSD of vertical turbulence is approximated by [10]:  
\[\hat{\phi}_{W_{z}}(\omega)=0.7V^{2/3}\varepsilon^{2/3}\omega^{-5/3} \tag{23}\]  
The energy spectrum of turbulence depends on the parameter \(\varepsilon\). According to turbulence model, the theoretical EDR value, generally defined as \(\varepsilon^{1/3}\) and commonly used to verify the severity estimation algorithm, is obtained. By inserting Equation (23) into Equation (20),  
\[\hat{\sigma}_{\alpha_{n}}^{2}=0.7V^{2/3}\varepsilon^{2/3}\int_{0}^{\infty} \left|H_{bp}(\omega_{1},\omega_{2},\omega)\right|^{2}\left|\frac{a_{n}(i \omega)}{W_{z}(i\omega)}\right|^{2}\omega^{-5/3}d\omega=0.7V^{2/3}\varepsilon^ {2/3}I(\omega_{1},\omega_{2},t) \tag{24}\]If \(I(\omega_{1},\omega_{2},t)\) in Equation (24) is known, the EDR value \(\varepsilon^{1/3}\) can be estimated. To estimate the acceleration response, Fourier transform is performed on the vertical acceleration and vertical wind sequence in the period \(T.\) After that, the results are converted from time domain to frequency domain to obtain \(a_{n}(i\omega)\) and \(W_{z}(i\omega)\) respectively. The value of \(\frac{a_{n}(i\omega)}{W_{z}(i\omega)}\) forms a new sequence. According to Wiener-Khinchin theorem, the autocorrelation function \(R\) of the new sequence by IFT is obtained with \(R(0)\) as the estimation results. After this, the acceleration response \(I(\omega_{1},\omega_{2},t_{k})\) can be acquired. Finally, EDR value is estimated by the following equation:  
\[\varepsilon^{1/3}(t_{k})=\frac{\hat{\sigma}_{a_{n}}(t_{k})}{\sqrt{0.7V^{2/3}( t_{k})I(\omega_{1},\omega_{2},t_{k})}} \tag{25}\]  
To summarize, an integrated computation flow for vertical acceleration and turbulence severity estimation is shown in Figure 6. The grid on wing-tail combination is generated beforehand. In the algorithm loop, external turbulence and flight parameters are used to compute local velocity by Equation (4) firstly. After that, the aerodynamic force is computed by nonplanar and unsteady VRM which has been assembled in the dashed box. In each cycle, by toggling the Flag, total aerodynamic force \(\mathbf{F}_{i}\) and the force without turbulence effects \(\mathbf{F}_{i}^{0}\) are computed, respectively to obtain the force increment \(\Delta\mathbf{F}\). Based on the latest incremental derivatives, vertical acceleration is obtained by solving the FDM. Next, the time series of acceleration data are used to estimate \(\hat{\sigma}_{a_{n}}(t_{k})\) by Equation (22). In addition, after fast Fourier transform (FFT), the acceleration series in frequency domain are used to obtain the transfer function \(I(\omega_{1},\omega_{2},t_{k})\) together with the frequency-domain vertical turbulence. At the end of each iteration, the estimated EDR value is output by solving Equation (25).