A330-300 Wake Encounter by ARJ21 Aircraft

Abstract

Today, aviation has grown significantly in importance. However, the challenge of flight delays has become increasingly severe due to the need for safe separation between aircraft to mitigate wake turbulence effects. The primary emphasis of this investigation resides in elucidating the evolutionary attributes of wake vortices in homogeneous isotropy turbulence. The large eddy simulation (LES) method is used to scrutinize the dynamic evolution of wake vortices engendered by an A333 aircraft in the atmospheric milieu and assess its ramifications on the ARJ21 aircraft. The research endeavor commences by formulating an LES methodology for the evolution of aircraft wake vortices, integrating adaptive grid technology to reduce the necessary grid volume significantly. This approach enables the implementation of axial and non-axial grid adaptive refinement, leading to more accurate simulations of both axial and non-axial vortices. Numerical simulations are conducted using the LES approach to scrutinize three distinct rates of turbulence dissipation amidst the ambient atmospheric turbulence, and the results are juxtaposed with Lidar measurements (Wind3D 6000 LiDAR) of wake vortices acquired at Chengdu Shuangliu International Airport (CTU). Subsequently, the rolling moment of the following aircraft is calculated, and three-dimensional hazard zones are determined for the A333. It is found that during the approach phase, the wake turbulence separation minima for an ARJ21 (CAT-F) following an A333 (CAT-B) is 3.35 NM, which represents a reduction of approximately 33% compared to ICAO RECAT (Wake Turbulence Re-categorization). The findings validate the dependability of the fine-grained mesh used in the vortex core region, engendered through the adaptive grid method, which proficiently captures the Crow instability and the interconnected phenomena of vortices in the numerical examination of aircraft wake. The safety of wake encounters primarily depends on the magnitude of environmental turbulence and the development of structural instability in wake vortices.

1. Introduction

Aircraft wakes are one of the leading causes of capacity problems in the air transport industry, and the lift forces exerted on an aircraft’s wings can create long-lived vortices in the wake. In particular, these may endanger any aircraft following closely behind during the critical landing phase of the aircraft [1], causing the following aircraft to be susceptible to pitching, rolling, and other safety compromises when entering the wake area of the leading aircraft. Therefore, a minimum safety separation must be ensured for aircraft taking off and landing on the same runway. ICAO, Europe, and the United States have established corresponding wake separation standards (RECAT in the United States and RECAT-EU in Europe) [2] to ensure flight safety. However, in the rapidly growing aviation environment, these aircraft safety separations limit the capacity of congested airports. Many scholars have recently conducted numerical simulation studies to reduce the minimum wake separation. Computational fluid dynamics (CFD) has been widely used as the primary research tool in the study of aircraft wake dissipation [3,4,5]. Crow et al. [6] found the wake vortex long-wave instability in 1970 and predicted that atmospheric turbulence would affect the wake vortex dissipation. Jongil Han et al. [5] used LES to calculate the wake vortex at different atmospheric turbulence intensities. Similarly, Proctor et al. [7] used LES to study the evolution of the wake vortex in stable and neutral atmospheres and gave an empirical model of the wake vortex evolution with the atmospheric turbulence dissipation rate ε and Brunt–Väisälä (BV) frequency N (which is related to the vertical gradient of temperature and is an atmospheric stability parameter) as parameters of the empirical model for the evolution of the wake vortex. In addition, Proctor et al.’s study also shows that more substantial atmospheric turbulence dissipation rates and BV frequencies will cause the wake vortex to enter into rapid decay more quickly. Proctor et al. [7] studied the wake vortex evolution in the near-surface phase under side winds and headwinds. Numerous studies have shown that side winds lead to an asymmetric decay of the two tail eddies, with the downstream eddies decaying faster than the upstream eddies [8,9,10].

Computer arithmetic has seen rapid advancements in recent years, primarily focusing on the study of the evolution law of two-dimensional wake vortices influenced by side winds, ground effects, and BV frequency in the initial stages, progressively transitioning to three-dimensional calculations [11,12,13]. The 3D numerical simulation of wake turbulence encounters several challenges: Firstly, the vast computational load; the strong shear in the wake vortex core region necessitates a high grid resolution. The movement of the wake vortex implies a requirement for global refinement. Additionally, to differentiate between long- and short-wave instabilities, the axial direction of the wake vortex demands a lengthy computational domain, approximately eight times larger than the axial computational grid of the vortex spacing [2,14,15]. This results in a massive total grid number for the wake vortex numerical simulation, with Misaka et al. [16] and De Visscher et al. [17] accounting for 60 and 125 million, respectively. Another challenge lies in the initialization method of the wake vortex, with the use of the entire aircraft winding simulation being computationally exhaustive [18]. Moreover, a significant headwind encourages the movement of the aircraft wake, necessitating enhanced grid precision and computational domain. The numerical simulation offers inferior accuracy for long-distance simulations of the wake vortex, with a significant surge in grid numbers as the computational domain expands, substantially inflating the simulation’s computational volume. This is beyond the computational capacity of contemporary computer simulations, necessitating sophisticated adaptive grid technology [19,20]. This research used a dynamic adaptive mesh LES method, using the finite volume CFD code Fluent for calculating the wake vortex evolution. This method captures the grids near the vortex core, enriches the grids around the vortex core region, and coarsens the grids farther away from the vortex core. Note that compared to the approach in [19], which also employed an adaptive grid LES for aircraft wake vortices, the proposed method differs significantly. The previous study focused solely on refining the grid axially. However, due to its grid intensification being restricted to the axial direction, it becomes less accurate in capturing non-axial vortices. In contrast, our proposed method conducts an analysis of the sensitivity to both axial and non-axial resolution, hence significantly improving the accuracy of simulating non-axis vortices. Ultimately, the proposed method reduces the overall grid size to about 5 million, enhancing computational efficiency and cutting down on computational costs.

In studying aircraft wake, Friedrich Köpp et al. [21] conducted field experiments at Tubbs Airport, validating that Coherent Doppler Lidar (CDL) could accurately measure wake vortices from their inception to the decay stages over an extended period. Smalikho et al. [22] developed a method to measure wake vortices by a 1.5 µm pulsed Coherent Doppler Lidar “Streamline” and proposed a way to estimate the volume and vortex core position of the wake from Lidar data. Consequently, this study uses previously acquired Lidar data from Chengdu Shuangliu Airport (CTU) [23,24,25], such as crosswind and turbulence dissipation rates, among other real-time meteorological information, to research the transmission and decay of aircraft wake under various turbulence dissipation rates. CDL measurements were carried out from 6 August 2018 to 22 October 2018, at CTU, with environmental conditions recorded simultaneously. Multiple aspects of the transmission and decay of wake vortices, such as velocity, circulation distribution, displacement, decay, kinetic energy, and lifespan, could be inspected by extracting them from the measured database.

As of 24 November 2023, the regional jet ARJ21 has celebrated a milestone, surpassing 10 million passenger trips. Since its official commercial operation commenced in June 2016, the Commercial Aircraft Corporation of China (COMAC) has successfully delivered 117 ARJ21 aircraft to various airlines. The aircraft has seamlessly integrated into the aviation landscape, serving over 400 domestic and international routes that span more than 140 cities. With nearly 1800 flights operated weekly, the ARJ21 has become a key player in the aviation sector.

Following the introduction of the ARJ21 into service in China, concerns arose about the potential impact of wake turbulence from other aircraft on the ARJ21. This led airports to implement larger separations for landing and take-off in their air traffic control procedures. While this cautious approach aimed to address safety concerns, it inadvertently limited runway capacity and impacted operational efficiency.

As the operations of the ARJ21 passenger jet continue to mature, there is a growing recognition of the need for dedicated research into wake turbulence specific to this aircraft model. Currently, there is a scarcity of literature addressing the effects of wake turbulence from leading aircraft on the ARJ21. Consequently, this study on wake turbulence separation holds significant importance in establishing and optimizing a dedicated wake turbulence separation system for the ARJ21, ultimately contributing to enhanced safety and operational efficiency in air traffic control.

Assessing the impact of the lead aircraft wake on a following aircraft typically involves using the Roll Moment Coefficient (RMC) threshold. Currently, the calculation of RMC predominantly adopts the strip theory method. Zhou et al. [26] divided the A320’s wings and horizontal tail into fourteen and six strips, respectively, incorporating the strip model into the wake flow field to approximately calculate the additional aerodynamic forces and moments induced by wake velocities, thereby evaluating the boundary of wake hazard during the cruising phase. Using wake dissipation models and aerodynamic response models, the ARJ21 aircraft was subdivided into 22 quadrilateral blocks for calculating aerodynamic forces and moments and analyzing its safety when encountering wakes from various aircraft [27,28]. To better reflect the impact of the A333 wake on the ARJ21 aircraft, this study divides the ARJ21 wing into 55 strips, equal to the number of grid nodes covered by the wing in the numerical simulation domain. Although this increases the computational load, it results in a more accurate RMC value for assessing the range of the wake hazard zone.

This study aims to investigate the wake vortices of the A333 aircraft by using LES. The second section introduces a dynamic self-adaptive grid to enhance the computational efficiency of the numerical simulations. The third section presents the background turbulence field and the initial wake vortex field. A comparative analysis between the numerical simulation results obtained under different turbulence dissipation rates and the wake vortex Lidar detection data from CTU is conducted in the fourth section. Furthermore, the variation in the RMC, a crucial metric for evaluating the rolling moment experienced by the aircraft, is examined by using the numerical simulation results when encountering the wake of the A333 aircraft. The RMC is used to visualize the hazard zones of the wake under varying turbulence dissipation rates.

2. Numerical Method and Adaptive Grid

2.1. Numerical Method

The wake evolution simulation used Fluent while using the Boussinesq approximation of the Navier–Stokes equations for LES. The formulation of the equations is as follows:

$$\frac{\partial u_j}{\partial x_j}=0$$

(1)

$$\rho \frac{\partial u_i}{\partial t}+\rho u_j\frac{\partial u_i}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\rho \left(v+v_t\right)\frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}$$

(2)

The formulation incorporates the subgrid-scale kinematic viscosity $v$. A Lagrangian dynamic model closes the equations to alleviate the excessive eddy viscosity within the vortex core [16]. Spatial discretization uses second-order upwind schemes, while temporal integration uses bounded second-order implicit schemes.

2.2. Wake Vortex Model

The strength of the wake is manifested through its circulation. An averaging process calculates the circulation, integrating the vorticity over 11 circular planes centered around the vortex core, with radii ranging from 5 to 15 m. The integration is performed along the direction of the aircraft’s forward flight. Furthermore, the initial circulation at t = 0 is used as a reference to normalize the results into non-dimensional values [8].

$${\Gamma }_{5-15}^{\ast }=\frac{1}{11}{\displaystyle \sum _{r=5}^{15}({\displaystyle {\int }_0^r\omega dA})/{\Gamma }_0}$$

(3)

The level of turbulence is characterized in this study by its dimensionless intensity.

$$w_0={\Gamma }_0/2\pi b_0$$

(4)

$${\epsilon }^{\ast }={\left(\epsilon b_0\right)}^{1/3}/w_0$$

(5)

In the three-dimensional simulations, the theoretical value of the vortex link time, TL, is used as a benchmark for comparing simulation results. Under the Lamb–Oseen initialization scheme, the vortex linking occurred more rapidly, whereas the BH initialization predicted the TL values that are closer to the theoretical expectations. Moreover, Section 5 of this paper focuses on studying and discussing the impact of the A330-300 wake on the following ARJ21 aircraft. Opting for the Lamb–Oseen model, which predicts faster wake vortex dissipation, might lead to lower wake turbulence separation standards. Nevertheless, reduced safety margins might entail greater safety risks. A comparative analysis of the hazard circulation, Γ_5–15m, for two different vortex core sizes is therefore conducted. The Γ_5–15m values obtained using the Lamb–Oseen model are significantly lower compared to those from the BH model [29]. An encounter response model based on Lamb–Oseen wake encounters underestimates the associated hazards. Given these considerations, the BH model is chosen for wake vortex initialization in this study. It offers a closer approximation to theoretical values, ensuring a more realistic and secure evaluation of wake impacts.

When the aircraft achieves vertical equilibrium, the vertical momentum of the wake vortex is balanced by the gravitational force exerted on the aircraft. The initial circulation of the wake vortex can be mathematically described as follows [30]:

$${\Gamma }_0=\frac{4MLWg}{BV\rho \pi }$$

(6)

The initial circulation of the wake vortex, denoted as ${\Gamma }_0$, is calculated using the maximum landing weight MLW of the aircraft, the local acceleration due to gravity g, the wingspan B, and the aircraft speed V. Simulations of the initial flow field of the wake vortex are performed using the Burnham–Hallock (BH) model [31]. The initial circulation ${\Gamma }_0$ determines the tangential velocity $V_{\theta }$, $r_c$ is the vortex core radius, and $r$ is the distance from the vortex center. The expression for $V_{\theta }(r)$ is given as follows:

$$V_{\theta }\left(r\right)=\frac{{\Gamma }_0}{2\pi }\frac{r}{r_c^2+r^2}$$

(7)

The initial vortex core radius $r_{co}$ is defined by Thomas Gerz [1] as the initial spacing between vortex cores, denoted as $b_0$:$b_0=\pi B/4$, $r_{c0}=0.052b_0$; the time scale, marked as $t_0$, characterizes the duration required for the downward motion of a sinking wake vortex generated by an aircraft or aircraft model, to cover a distance equal to the initial spacing between vortex cores, $t_0=b_0/w_0$, $t^{\ast }=t/t_0$. The Reynolds number, defined based on circulation, is denoted as $Re={\Gamma }_0/v$ [32].

2.3. Ambient Turbulence Field by Rogallo Spectrum

The three-dimensional spectrum is an essential characteristic of homogeneous isotropic turbulence and is expressed as follows:

$${\displaystyle {\int }_0^{\infty }E(k)dk=K}$$

(8)

$$2v{\displaystyle {\int }_0^{\infty }{k}^{2}E(k)dk=\epsilon }$$

(9)

Moreover, the modified Von Karman spectrum simulates the complete energy spectrum,

$$E(k)=A\frac{2/3K}{k_e}\frac{{(k/k_e)}^4}{{[1+{(k/k_e)}^2]}^{17/6}}\exp [-2{(k/k_{kol})}^2]$$

(10)

In Equation (10), with $k_e=2\pi /100$ and $k_{kol}={(\epsilon /v^3)}^{1/4}$, the product of A and K can be determined using Equations (8) and (9).

The method introduced by Rogallo [33] was used in this study to generate velocity fields in the spectral domain,

$$\stackrel{⌢}{u}(k)=\frac{1}{k\sqrt{k_1^2+k_2^2}}\left(\begin{array}{c}\alpha \left(k\right)k{k}_2+\beta (k)k_1{k}_3\\ \beta (k)k_{2}k-\beta (k)k{k}_1\\ -\beta (k)(k_1^2+k_2^2)\end{array}\right)$$

(11)

$$\begin{array}{l}\alpha (k)=\sqrt{\frac{E(k)}{2\pi k^2}}\exp (i{\theta }_1)\cos \varphi \\ \beta (k)=\sqrt{\frac{E(k)}{2\pi k^2}}\exp (i{\theta }_2)\sin \varphi \end{array}$$

(12)

The variables ${\theta }_1$, ${\theta }_2$, and $\varphi$ in Equation (12) represent random numbers that conform to a uniform probability distribution within the interval $[0,2\pi ]$,

$$u(x)={\displaystyle \sum _{k_1}{\displaystyle \sum _{k_2}{\displaystyle \sum _{k_3}\stackrel{⌢}{u}(k)\exp (ik,x)}}}$$

(13)

The equation mentioned earlier generates the background turbulence field within a 256 × 256 × 256 grid. This generated flow field is then interpolated into a 400 × 256 × 256 grid, where the grid spacing remains uniform at Δ = 2 m. Precursor simulations are followed to obtain a fully developed turbulent flow.

Based on the energy and dissipation spectra acquired through LES, as shown in Figure 1, the spectrum of turbulent kinetic energy gradually declines as wavenumbers increase, indicating a power-law relationship that aligns with Kolmogorov’s proposed cascade mechanism [34]. Specifically, during the initial time point (t* = 1.56), the energy spectrum exhibits a prominent concentration of energy at larger scales. However, as time progresses, the decaying trend of the energy spectrum slows down, suggesting a transfer of energy toward smaller-scale structures. This observation provides supporting evidence for the existence of an energy cascade and the phenomenon of energy transfer from larger to smaller scales. Notably, Figure 1 shows a relatively modest rate of energy decay over time, which does not significantly undermine its suitability as a background turbulence field for modeling the dissipation of wake vortex. Thus, it is an appropriate choice for numerically simulating wake vortex dissipation.

In addition, the process of energy dissipation through viscosity in turbulence is described by the dissipation spectrum. Observations of the dissipation spectrum reveal a gradual increase with increasing wavenumbers, following a power-law relationship. Within the lower wavenumber range, the slope of the dissipation spectrum is relatively small, indicating a more efficient energy transfer at larger scales. However, as wavenumbers increase, the slope of the dissipation spectrum gradually steepens, highlighting the dominant role of viscosity in dissipating energy at smaller scales and converting it into thermal energy. Over time, the turbulent energy spectrum undergoes energy transfer. At the early point (t* = 1.56), the energy is predominantly concentrated within smaller scales, corresponding to the high-wavenumber region. As time progresses, the turbulent energy gradually shifts toward larger scales, corresponding to the low-wavenumber region. Consequently, at the following time points (t* = 4.69, t* = 6.25), the energy gradually increases in the low-wavenumber range while decreasing in the high-wavenumber range. Concurrently, with time, the dissipation spectrum exhibits distinct characteristics, with the dissipation rate progressively increasing, indicating a more pronounced dissipation of turbulent energy at different scales.

Figure 2 shows the iso-surface of the decaying three-dimensional turbulent field. All iso-surfaces are color-coded based on the magnitude of |ω|. The vortex-surface field (VSF) iso-surfaces with |ω| > 0.5|ω|_max are truncated for clarity. In addition, a series of distinct multi-scale patterns can be observed in Figure 2c, showing the presence of eddies of various sizes within the turbulence, ranging from larger vortices to smaller ones. This multi-scale presence signifies the sufficient development of turbulence across different spatial scales. Furthermore, these patterns exemplify the intermittent nature of turbulence, as they manifest fluctuations in eddy intensity between strong and weak phases. This intermittent behavior is a significant manifestation of the dynamic nature of turbulence, indicating its continuous evolution, reorganization, and transfer of energy within the fluid domain, ultimately leading to the formation of turbulent structures. Moreover, exceedingly thin filaments have also been observed within the turbulent flow field shown in Figure 2a,b. These filaments exhibit remarkably subtle motion characteristics, potentially representing smaller-scale turbulent structures or localized regions of concentrated turbulence energy. These delicate filamentary structures further solidify the extent of turbulence development within the fluid domain.

2.4. Grid Independence

During the grid independence verification, the BH model was used for simulating the velocity distribution of three-dimensional vortices. In order to investigate the impact of different grid resolutions on the simulation outcomes, four grid resolutions (r_c0/Δ = 2.5, 5, 8, 12) were selected for comparative analysis. Emphasis was placed on comparing the computed results with the theoretical curve of the BH model, mainly focusing on the velocity peak. Figure 3a depicts the velocity distribution curves for the four different grid resolutions alongside the theoretical curve of the BH model. The graph reveals a gradual convergence of the computed results toward the ideal curve as the grid resolution increases. Notably, bigger grid resolutions (such as r_c0/Δ = 8 and 12) exhibit enhanced capability in capturing velocity variations, particularly near the velocity peak. However, it is essential to consider computational resource limitations in practical applications. Hence, in this study, we selected a grid resolution of r_c0/Δ = 5 as the reference value. Despite its slightly reduced precision compared to bigger grid resolutions, it still yields relatively accurate simulation results while ensuring computational efficiency. Note that the BH model, chosen for wake vortex initialization in this paper, dictates the peak of tangential velocity. Consequently, this may lead to differences in grid sensitivity based on the selected model or its parameters. Further investigations are required on this topic.

To further quantify the impact of varying grid sizes on the velocity peak, we show Figure 3c, which provides two informative components. Firstly, the line graph shows the velocity peak’s ratio to the BH model’s ideal value. As the grid size decreases, the velocity peak approaches the theoretical value of 100%. Moreover, the grid size curve tends to smoothen for grid sizes smaller than 0.5 m, suggesting that further reduction in grid size has a limited effect on enhancing the simulation results within this range. Secondly, the bar chart shows the percentage difference between the actual velocity peak and the ideal value. The figure shows that smaller grid sizes exhibit a more negligible percentage difference at the velocity peak, indicating a closer match between the simulation results and the ideal value. By integrating the findings from Figure 3a,b, this study shows the influence of grid size on the simulation outcomes through a comparative analysis of velocity distribution curves and the percentage difference at the velocity peak across different grid sizes. In order to strike a balance between the use of computational resources and the accuracy of simulations, a grid size of 0.5 m was chosen as the reference standard. Future studies aim to establish a suitable benchmark for grid size.

2.5. Computational Domain and Boundary Conditions

Figure 4 shows the simulation domain, which has dimensions of Lx × Ly × Lz = 8.5 b₀ × 5.5 b₀ × 5.5 b₀, matching the size used by Hennemann [35]. The total number of grid points is Nx × Ny × Nz = 400 × 256 × 256, where x, y, and z correspond to the axial, lateral, and vertical directions, respectively. The axial dimension of 8.5 b₀ is chosen to be sufficiently large to capture potential Crow instability [6]. Periodic boundary conditions are implemented in all three directions: lateral, axial, and vertical. The simulation uses a fixed time step of Δt = 0.01 s.

2.6. Dynamic Self-Adaptive Grid

The cross-sectional views of the refined wake vortex mesh in the y-z and x-z planes after initialization are shown in the left and right sections of Figure 5, respectively. Dynamic adaptive mesh refinement technology has been implemented within the aircraft wake region, using Fluent to enhance the mesh and accurately represent its aerodynamic characteristics. This technique autonomously modifies the grid density by considering the local distribution of vorticity, a scalar quantity that characterizes the local rotation in fluid flow. The computation of vorticity follows the formula $\omega =\sqrt{{\omega }_x^2+{\omega }_y^2+{\omega }_z^2}$. Here, ω_x, ω_y, and ω_z represent the vorticity components along the x, y, and z directions, respectively.

In this study, Fluent was used to enhance the grid density in the aircraft wake region, aiming to improve the accuracy of aerodynamic simulations. Dynamic self-adaptive meshing technology, or Dynamic Adaptive Meshing, was used to adjust the grid density to achieve this objective automatically. The mesh refinement or coarsening was performed based on the local vorticity distribution, enabling precise control over the vortex core mesh structure. A vorticity threshold was set as a criterion to determine whether a grid cell required refinement or coarsening. This threshold ensured that the grid structure achieved sufficient refinement in the vortex region while avoiding excessive grid refinement in regions of low vorticity, thereby conserving computational resources. Figure 5 shows the execution of two rounds of adaptive refinement/coarsening processes. These multiple iterations enable a more accurate capture of details in the aircraft wake region, further enhancing the accuracy of the simulation results. An optimization strategy was also implemented to allocate computational resources effectively, reducing computational costs. This optimization approach facilitates the attainment of high-precision simulation results for the aerodynamic characteristics of the aircraft wake region within the constraints of limited computational resources. Through the application of dynamic adaptive meshing technology, the grid refinement of the aircraft wake region was achieved, leading to improved accuracy of the simulation results. Furthermore, optimizing computational resource allocation reduced computational costs, making it possible to obtain high-precision simulation results for the aerodynamic characteristics of the aircraft wake region using limited computational resources.

Figure 6 shows the temporal evolution of the vortex structure and its impact on the computational domain mesh. Specifically, the investigation focuses on the evolution of the computational domain mesh for the case of ε* = 0.01. The following key observations can be summarized: Firstly, at t* = 5, the vortex structure deviates from the reference state and transforms into a vortex ring due to the Crow instability. This vortex ring exerts a significant influence on the flow field. Secondly, a secondary structure emerges in the vortex flow field at t* = 6 as time advances. This occurrence indicates the burst of the vortex, further enhancing the development of the Crow instability and the formation of vortex rings. To accurately capture these intricate phenomena, adaptive meshing expands the region with refined mesh resolution. Moreover, the mesh refinement follows the displacement of the vortex core, enabling the gradual expansion of the refined region to accommodate the changing shape of the wake vortex.

3. Simulation Results

The descent of the wake vortex under mutual interaction is shown in Figure 7a. In all scenarios, the initial descent rate is approximately 1. As the circulation gradually diminishes, the descent speed progressively decreases. The vortex, characterized by a smaller circulation, undergoes a slower descent when it forms vortex chains through the process of rolling up. Notably, when ε* = 0.23, the decay in wake vortex intensity causes a reduction in the strength of the mutual induction effect, resulting in a decreased descent speed.

As shown in Figure 7b, the increase in spacing between wake vortex cores could be an additional factor contributing to this deceleration, as it weakens the mutual interaction between the two wake vortices. Specifically, in the case of ε* = 0.05, the formation of vortex linking results in a heightened interaction between the wake vortices, decreasing the spacing between their vortex cores.

In addition, Figure 7c shows the temporal decay of wake vortex circulation, with the circulation being normalized by its initial value at t = 0. The height is normalized using b₀, and time is normalized using t₀. Figure 7 shows the characteristic two-stage decay of wake vortex circulation. During the initial diffusion stage, the decay rates of circulation are nearly identical for all cases. The non-dimensional time signifies the transition to the second stage of rapid decay, which is influenced by meteorological conditions. A comparison of the cases with ε* values of 0.01, 0.05, and 0.23 reveals that they decrease as the turbulence dissipation rate increases, indicating that a higher dissipation rate accelerates the transition to rapid decay. Compared to Lin’s [19] use of the Gnoffo adaptive mesh algorithm, this study exhibits a slower dissipation rate under strong turbulence, primarily due to the larger grid size used for initializing the turbulent field, which inadequately captures strong turbulence eddies. However, during the initial diffusion stage, the decay rate in this study is slightly lower than the results obtained by Hennemann, who used an equidistant grid and a more potent standard Smagorinsky subgrid-scale model. The grid size near the vortex core in this study (0.5 m) is smaller than that of the equidistant grid (1 m) [35]. Previous research has indicated that a smaller grid size can induce numerical pseudo-diffusion near the vortex core, accelerating the decay rate during the diffusion stage. Nevertheless, this effect minimizes the wake vortex’s time to enter the rapid decay stage [16].

Figure 8 shows the development of wake vortex Crow instability in turbulent conditions (ε* = 0.05), revealing the formation of vortex linking at t* = 8, consistent with the findings of [17]. The decay of vortices in moderate turbulence and weak conditions is predominantly governed by Crow instability [6]. Furthermore, Figure 8 shows the limited occurrence of vortex linking in strong turbulence (ε* = 0.23) due to the presence of short-wave instability. This limited linking can be attributed to energy dissipation within turbulent vortices caused by the aforementioned short-wave instability, thereby hindering the vortices from maintaining sufficient strength for sustained linking. As the vortices approach the linking phase, the primary vortices undergo rapid elongation and transform into secondary vortices, analogous to the generation of secondary vortices through tilting, which experiences rapid dissipation during the later stages of their lifespan [16].

Figure 9 shows the self-induced interaction among wake vortices, resulting in the formation of a vortex chain. The vortex on the right displays greater vorticity and higher rotational intensity. The vortex on the left is attracted toward the right vortex due to the mutual interaction between the vortex tubes (Figure 9a). Vortices induce modifications in the surrounding flow, generating an attractive force between the vortex tubes. As the left vortex approaches the right vortex (see Figure 9b,c), their interaction strengthens, leading to the development of a tightly linked vortex chain structure. Vortex chains typically align themselves along their axis, creating continuous rotating formations. Despite the larger vorticity of the right vortex, the rotational direction of the left vortex remains unchanged. Still, it shifts from a horizontal position to the lower-left region relative to the right vortex (Figure 9d). The higher rotational intensity of the right vortex exerts a more pronounced influence on the surrounding fluid, causing the positional deflection of the left vortex under its effect. Eventually, the left vortex resides below and to the left of the right vortex. The formation of the vortex chain causes the interaction between the vortex tubes to pull the left vortex toward the lower region of the right vortex. In contrast, the influence of the right vortex further enhances the deflection of the left vortex, reducing the coverage area of the dense grid region.

4. Lidar Measurements

The parameter values of the A333 landing, retrieved through Lidar measurement at CTU, align closely with the numerical simulations assuming the A333 wake. The wake vortex, characterized by its rapid decay, may only remain detectable for a limited duration. Figure 10a compares the tangential velocity distribution of the initial wake vortex, computed using the BH model, and the results derived from Lidar.

Figure 10a shows the radial distribution of tangential velocity in the wake vortex, as measured by the Lidar. The centerline is the reference for the vortex core position, where positive values indicate radii above the horizontal line of the vortex core, and negative values indicate radii below it. The upper and lower curves represent the maximum and minimum velocities, respectively. Based on the legend, it can be observed that the left vortex rotates in a clockwise direction, with positive radial values in the upper half indicating positive tangential velocities. The upper curve exhibits a pattern of increasing and then decreasing tangential velocities within the 0–15 m range, showing a close agreement with the theoretical distribution. Similarly, the lower curve demonstrates a consistent pattern between the measured tangential velocities within the range of −15–0 m and the theoretical distribution.

Figure 10b compares the values obtained in this study and those predicted by the ref. [35] and Sarpkaya models [36,37]. The Sarpkaya model offers an empirical approach to estimation under neutral atmospheric conditions. However, the Sarpkaya model assumes wake vortex linking based on the theoretical premise of long-wave instability, and in the case of ε* = 0.23, vortex linking did not occur (as previously discussed), resulting in some deviation between the results of the ε* = 0.23 case and the Sarpkaya model. In conclusion, despite using a grid resolution that is only 11% of the one used in ref. [35], this study successfully achieves a high level of accuracy in simulating the decay of circulation.

Figure 11 shows the comparative analysis between the retrieved vorticity derived from Lidar measurements, the vorticity obtained through large-eddy simulations, and the fitting results of the BH model. The observed trends in these three data sets are consistent, indicating a good agreement. Hence, in scenarios where the vortex core is positioned at a lower location and the tangential velocity distribution is incomplete, the BH model can estimate the vorticity at the missing tangential velocity distribution location based on available vorticity observations. This approach provides a viable solution to address the challenge of vorticity inversion under varying environmental conditions.

Twenty cases of wake vortex Lidar measurement data for the A333 were extracted from CTU to optimize the Sarpkaya model, focusing on data collected under high-turbulence conditions. The data were fitted using an exponential decay function $\Gamma =a\cdot \exp (-b\cdot t)$. Among these parameters, the function’s magnitude is controlled by the fitting parameter $a$, while the fitting parameter governs the rate of vorticity decay $b$. Figure 12a shows the dissipation formula of circulation over time, denoted as $\Gamma =605.071\cdot \exp (-0.061\cdot t)$, obtained from the fitting curve, which closely aligns with the Sarpkaya model. Furthermore, Figure 12b shows a favorable level of agreement between the results and the numerical simulations.

5. Hazard Zones

5.1. Roll Moment Coefficient

During the evolution of wake vortices, they undergo deformation, giving rise to axial variations in their spatial structure. Consequently, the rolling moment fluctuates depending on the aircraft’s position along the flight path. The highest rolling moment, indicative of the most pronounced impact from the wake vortex field, is encountered when the aircraft aligns with the vortex core. As the wake vortices exhibit spatial variability along the flight path, the corresponding hazardous region within the wake vortex field also changes. Greater spatial extents of the wake vortices correspond to larger hazardous areas, while smaller spatial extents result in relatively smaller hazardous areas. Considering the potential encounters of the following aircraft with the wake vortex field at any given position, it is rational to position the aircraft at the location with the maximum spatial extent of the vortex to calculate the rolling moment and define the hazardous area, thereby ensuring the safety of the following aircraft. Figure 13 shows the simulated positioning of the following aircraft along the flight path within the computational domain.

In this study, the ARJ21 aircraft was chosen as the following aircraft, as shown in Figure 14. Upon entering the wake vortex field generated by the leading aircraft, the following aircraft experiences the influence of the downwash velocity field induced by the wake vortices. This influence encompasses alterations in lift and rolling moment induced by the wake vortices, with the induced rolling moment being the primary concern for flight safety. Hence, the RMC induced by the wake vortices is widely used as a safety criterion for the following aircraft [38]. This criterion guarantees that the decay of the wake vortices from the leading aircraft reduces the induced rolling moment coefficient to within the tolerable range for the following aircraft. To simplify the analysis, the following aircraft is represented as a flat-plate wing with an elliptical lift distribution, where the span of the following aircraft’s wing is denoted as B_F. The coordinates of the following aircraft crossing the predicted wake vortex plane are defined as (y₀, z₀). Assuming a lift line slope of $2\pi$ [39], the change in angle of attack $d\alpha$ and the change in induced lift dL caused by the downwash velocity field on a wing element can be calculated as follows:

$$d\alpha =u(y,z)/V_F$$

(14)

$$dL=2\pi c\frac{\rho V_F^2}{2}{d}_y{d}_{\alpha }$$

(15)

In Equation (15), the variable u(y, z) represents the downwash velocity field at a specific moment within the predicted wake vortex plane. In contrast, the variable V_F corresponds to the approach airspeed of the following aircraft. The induced rolling moment M(y, z) experienced by the following aircraft can be computed using the following equations:

$$S=\pi c_0{B}_F/4$$

(16)

$$M(y_0,z_0)={\displaystyle {\int }_{y_0-\frac{B_F}{2}}^{y_0+\frac{B_F}{2}}(y-y_0)}dL$$

(17)

The RMC [40] can be obtained from the above equation, using its inherent attributes of perceptual aircraft non-dimensional and discriminative properties.

$$\mathrm{RMC}=\frac{M(y_0,z_0)}{0.5\cdot \rho \cdot V_F^2\cdot S\cdot B_F}$$

(18)

The remaining parameters are defined as follows: the wingspan B_F = 27.3 m, the parameter c₀ = 2.95, the approach airspeed of the following aircraft is denoted by V_F = 70 m/s, and the wing area of the following aircraft is specified as S.

The distribution of the RMC under different turbulence dissipation rates is depicted in Figure 15. It is noted that the utmost rolling moment manifests in the proximity of the vortex cores, signifying the most substantial and hazardous impact of the wake vortex field on the following aircraft. The rolling moment converges to zero as the aircraft approaches the midpoint between the two vortices. When situated at the same altitude as the vortex core but at a distance from its lateral boundaries, the rolling moment initially diminishes to zero, subsequently undergoes a slight augmentation in the opposite direction, and then diminishes once more. To evaluate this risk, a criterion for the rolling moment coefficient ranging from 0.05 to 0.07 is adopted [40]. If the rolling moment coefficient surpasses 0.05, it has the potential to result in an unmanageable flight of the following aircraft. Therefore, 0.05 is deemed an impermissible threshold for assessing the safety of the aircraft in various turbulence conditions. Despite the irregular distribution of the rolling moment experienced by the following aircraft, the genuine hazardous boundary ought to exhibit a curved configuration, thereby augmenting the intricacy of flight path regulation. To streamline the flight path management, the hazardous boundary is approximated as a rectangular region demarcated by four vertical and horizontal constraints. Consequently, the aircraft must maintain a secure range of RMC within this rectangular region to avert the initiation of unmanageable circumstances near the vortex cores.

5.2. Hazard Zone

The temporal variations in the contour surfaces obtained from the computed RMC are shown in Figure 16. It is apparent that under low turbulence dissipation rates (ε* = 0.01), the wake generated by the leading aircraft spans a substantial area. Figure 17b shows the prolonged axial expansion of the wake, concomitant with an extended lifespan of the vortices. Moreover, the vertical profile shown in Figure 17a shows a similar downward behavior of the vortex core, as observed in Figure 7a, showing a reverse relationship between the rate of wake descent and the spacing between neighboring vortices.

In Figure 17a, a safety corridor spanning a width of 60 m in both positive and negative directions is shown to account for the vertical positional uncertainties of the leading aircraft, thus augmenting the safety margin for the following aircraft during wake turbulence encounters. Irrespective of the turbulence dissipation rates, the prescribed safety separation for the ARJ21 following the A333 remains stable. This consistency can be ascribed to the complete departure of the hazard zone from the A333 safety corridor during the wake vortex diffusion phase. The discrepancies in vortex intensity and dissipation trajectories among these three scenarios are negligible. Compared to the prevailing ICAO RECAT standards, the spacing requirements are diminished by approximately 33%. However, it is imperative to acknowledge that this reduction is solely predicated on considering turbulence dissipation rates, neglecting the impacts of BV frequency, crosswind, and headwind. Therefore, these findings can serve as a reference point for formulating dynamic separation guidelines for aircraft of the CAT-B and CAT-F categories (see

Table 1. Wake turbulence categories for leading and following aircraft.

Aircraft Type	ICAO	ICAO RECAT
ARJ21	Medium	CAT-F
A330-300	Heavy	CAT-B

) in atmospheric conditions.

6. Conclusions

The large eddy simulation (LES) was undertaken to examine the evolutionary evolution of wake vortices engendered by the A333 aircraft. The LES used dynamic adaptive mesh refinement techniques to enhance accuracy and efficiency. The mesh refinement strategy concentrated on augmenting the grid density near the vortex core region and integrating axial grid refinement. After two refinement iterations, the grid dimensions within the vortex core region were established as Δx × Δy × Δz = 0.5 m × 0.5 m × 0.5 m. Conversely, the grid resolution in the far-field region was diminished, substantially reducing the overall grid count. The equidistant grid (Δ = 1 m) encompassed around 24.4 million grid points, while the refined mesh approach reduced approximately 20 million grid points. This significant reduction considerably enhanced computational efficiency and mitigated computational costs.

The findings derived from three distinct turbulence dissipation rates offer compelling evidence of the primary diffusion phase and rapid dissipation phase of wake vortex circulation. The observations about the wake vortex linking time, descent velocity, and vortex spacing exhibit excellent concurrence with those acquired through a globally refined grid. Moreover, notable characteristics of vortex evolution, encompassing secondary vortices and helical structures, are prominently discernible in the findings.

The Lidar measurement data obtained at Chengdu Shuangliu International Airport (CTU) was extracted and juxtaposed with the results derived from LES regarding the vorticity distribution, velocity, and lifespan parameters of the wake vortex engendered by the A333 aircraft. The LES findings aligned remarkably with the Lidar measurement data, indicating qualitative and quantitative coherence concerning vortex decay, transportation, and topological structure. Moreover, the Lidar measurement data inferred a rapid dissipation formula for high-turbulence conditions for the A333, and the resulting curve showcased a favorable alignment with the Sarpkaya model.

The CAT-F aircraft ARJ21 was selected as the following aircraft to be paired with the CAT-B aircraft A333. The strip model was used to assess the aerodynamic response characteristics of wake encounters. The safety of the following aircraft was validated through the evaluation of the Roll Moment Coefficient (RMC). The wake evolution distribution was examined using LES to investigate the hazard zone and wake separation across three distinct turbulence dissipation rates. The findings unveiled a significant decrease in the safety margin by the Wake Turbulence Re-categorisation (RECAT) standards. Evaluating the risk of wake encounters for the ARJ21 aircraft can provide a foundation for the prompt anticipation of wake separation. However, the evaluation currently lacks encounters with other leading aircraft categories. Specifically, it is imperative to integrate BV frequency, crosswind, headwind, and other environmental factors to simulate the authentic evolution of wakes in atmospheric conditions, thus augmenting the dependability and veracity of numerical simulation outcomes.