# Approach and Landing Aircraft Wake Encounter Risk Based on Reynolds-Averaged Navier-Stokes Numerical Simulation

Authors: Weijun Pan; Jingkai Wang; Yaxing Xu; Qianlan Jiang; Yuming Luo
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

In order to shorten the wake safe separation, numerical simulation technology with aerodynamic response models and strip models have been combined to calculate wake hazard zone. As a realistic case, a medium aircraft ARJ21 following a heavy aircraft A330-200 is considered, and the Reynolds-averaged Navier-Stokes (RANS) method is used to explore the wake vortex evolution process of the leading aircraft at the decision height. A strip algorithm is proposed to calculate the rolling moment coefficient and overload increment of the ARJ21 after encountering the wake of the front aircraft in the three-dimensional space. The proposed algorithm identifies the area where the wake of the front aircraft can cause risks to the following aircraft and analyzes the evolution process of the hazard zone of the section where the decision height is located. The minimum safe separation of the ARJ21 following the A330-200 is 1.32 nmile, which is 26.4% of the ICAO separation standard of 5 nmile. When the average runway occupancy time (ROT) is reduced to match the separation of 1.32 nmile, the capacity of runway 02R/20L in Tianfu International Airport could theoretically reach 102.37 sorties/h under this aircraft pair combination. Compared to original 27 sorties/h, the runway capacity can be improved up to 279.14%, which will increase the airport operation efficiency.


---

## Body

## 1. Introduction

During the period 1983–2002, accidents caused by wake during the approach and landing phases accounted for 70% of the total number of wake accidents, among which the medium aircraft followed the heavy aircraft on the approach in most cases [1]. To ensure safety, International Civil Aviation Organization (ICAO) has imposed strict regulations for wake separations during the take-off and landing phases [2]. However, due to the inconsistency in the wake intensity between different aircraft types, the wake intensity that can be tolerated also differs, so aircraft separation cannot be strictly defined. Therefore, to ensure flight safety, the take-off and landing efficiency of an airport is sacrificed. In the final approach stage, an aircraft is close to the ground, and the wake will undergo special evolution due to the influence of ground effects. Therefore, how to predict the wake structure in the final approach channel of an aircraft is of great significance in shortening the wake separation. An ARJ21 aircraft is a new medium- and short-range turbofan regional medium-size passenger aircraft. It has the layout of tail-mounted twin-engine, a T-shaped flat tail, and a low monoplane [3]. The aerodynamic force and torque an ARJ21 aircraft receives are different from most traditional jet aircraft. For lacking of research on the response of an ARJ21 aircraft to the front aircraft wake vortex, ARJ21 aircraft maintains a large wake separation with the front aircraft due to the lack of wake encounter risk assessment in flight operation. Therefore, it is necessary to shorten the wake separation while ensuring safety.Kovalev and Vanhoenacker-Janvier [4] detected the wake turbulence under clear atmospheric conditions and conducted the large-eddy simulation based on the wake vortices in turbulent stratification. It is shown that the high-energy dissipation rate within the WV oval provides higher values of the RCS relative to the clear atmosphere turbulence. Stephan et al. [5] conducted the mixed Reynolds-averaged–large-eddy simulation research on a landing aircraft. They analyzed the impact of the landing gear on the aerodynamics of an aircraft during the landing phase. The result shows that the landing gears have a notable effect on the aircraft’s drag, as well as on the wake footprint. However, the effect on the wake vortex circulation evolution is moderate. Cui et al. [6] used adaptive grid technology in the large-eddy simulation test of the aircraft wake vortex evolution. This approach reduced the number of grids and improved the calculation efficiency of numerical simulations significantly. At the same time, they also proposed a lift surface wake vortex generation method, which was used to study the evolution characteristics of aircraft wake vortex in the atmosphere. Based on this method, a rapid wake prediction system was constructed. Zhou et al. [7] used the Euler-Eulerian multiphase flow model to study the evolution characteristics of the wake vortex under rainfall conditions. Campo et al. [8] proposed a calculation formula of the safety separation between two aircraft and analyzed the roll response by calculating the rolling moment of the aircraft after encountering a wake. They also used the maximum roll angular velocity as an indicator of the wake encounter index severity. Gerben et al. [1] used the dimensionless RMC (roll moment coefficient) to measure the severity of wake encountered by various types of aircraft. Through simulation tests, it was verified that the RMC classification under the reasonable worst case of the current standard was acceptable. Holzäpfel et al. [9] collected meteorological and traffic data at Vienna Airport for 12 months. They used a wake prediction system to analyze the potential for separation reduction based on the collected data. Any combination of aircraft that requires wake separation has a certain potential to reduce the separation minimum safely. Pan et al. [10] considered factors such as runway configuration, crosswind, ground effect, and wake bearing capacity of the following aircraft and established the corresponding wake encounter response model. The research results that the current wake separation had a certain reduction potential.The contributions of this paper are mainly reflected in the following three aspects: First, the strip-like algorithm is designed to calculate the aerodynamic response characteristics of the following aircraft precisely based on the wake vortex-induced velocity distribution in this paper. The aircraft is subdivided into innumerable strips, the vertical resultant velocity of the wake vortex field acting on each strip is calculated, and the angle of attack caused by the wake on each strip is calculated according to the lift line theory. In this way, the additional aerodynamic torque generated by the following aircraft encountering the wake at any position can be accurately calculated. Second, the influence of the position change of the tail vortex of the front aircraft on the aerodynamic characteristics of the following aircraft has been considered adequately. Then, the hazard zone is divided according to the influence degree of value on the following aircraft. Combined with the approach segment of the following aircraft, a new minimum safety separation standard is proposed to improve the runway capacity. At the same time, the change in the wake hazard zone of the section at the DH (decision height) in the latest approach segment is analyzed to provide a reference for dynamic wake separation.


## 2. Wake Vortex Evolution at DH

### 2.1. Near-Field Evolution of Wake Vortex

#### 2.1.1. Governing Equations and Turbulence Model

In this paper, the RANS method decomposes the turbulent motion into two parts: pulsation and average motion. By solving the time-averaged N-S equation of turbulence, the turbulent motion can be rapidly simulated. Compared with direct numerical simulation, the large-eddy simulation (LES) can save more computing resources. The SST model [11] is selected as a turbulence model. Compared with the standard model, this model has higher calculation accuracy and is suitable for rotating flow and other situations.


#### 2.1.2. Computing Domains and Grids

The DH is the specific altitude at which a go-around should be initiated during a precision approach if the necessary visual reference to continue the approach cannot be established. Therefore, it is an important parameter of the precision approach stage and is a critical point for landing and overtaking. The study of the wake vortex evolution at the DH has important reference significance for preventing the following aircraft from entering the wake hazard zone and establishing the dynamic wake separation. The ILS (instrument landing system), which is considered a blind landing system, is the most widely used aircraft precision approach and landing guidance system. The ILSs can be roughly classified into three types based on the boot performance. The DH values of the ILS types are presented in Table1.Table 1
Classification and performance standards of the three ILS types.
Performance standardILS typeDH (m)IIIIIIABC6030Below 30 or 0Below 15 or 00Therefore, the wake vortices at different DHs corresponding to different types of precision approaches are selected for research in this study. The altitudes of an aircraft above the ground are set toH=60m, H=30m, H=10m, and H=5m. The aircraft model used in this study is A330-200, having a wingspan of 60.3 m, a wing area of 361.6 m2, and a fuselage length of 58.82 m. Since the horizontal tail vortex dissipates rapidly and has a slight effect on the safety of the following aircraft [12], only the wing part is retained in the calculation process. The computational domain meshing is completed using a structured mesh based on the finite volume method, and the wing surface is locally refined to meet the accuracy requirement of the calculation. The boundary layer mesh height y of the wing surface mesh is calculated by
(1)y=y+μVρ,where the wing surface y+ is one; μ is the dynamic viscosity of the turbulent flow; V is the fluid velocity; and ρ is the fluid density.The computational domain with an aircraft as the origin is presented in Figure1(a), where the x-direction denotes the approach direction of the aircraft, the y-direction is the wingspan direction, and the z-direction is the height direction. The wing is 3 b0 from the entrance, 2 b0 from the top surface, and 10 b0 from the exit, thus having enough spatial scales to distinguish between the long and short wave instability of crow [13]. By setting different distances between an aircraft and the lower wall, the ground effect at different heights can be studied. The structured grid of the fluid domain is presented in Figure 1(b). To capture the evolution and development of the wake vortex accurately, special grid refinement is performed in the wake region to ensure calculation accuracy. In addition, the mesh of the vortex area should have a good transition to ensure that the airflow boundary layer is accurately resolved. Therefore, the aspect ratio of the grid in the wake region is close to one, and the minimum grid size in the wake vortex area is 0.1m×0.1m.Figure 1
Computational domain and its grid. (a) Schematic diagram of computational domain and (b) computational domain grid.
(a)(b)The leading edge mesh of the wing is presented in Figure2(a). According to Equation (5), the boundary layer is densified to satisfy the condition of y+=1; the height of the first layer of the boundary layer is 6.14e-6 m, and the boundary layer has 40 layers. Figure 2(b) shows the surface mesh of the winglet. The surface mesh of the wingtip is refined to capture the evolution of the wingtip vortex better.Figure 2
The A330-200 wing partial grid. (a) Grid at the leading edge of the wing and (b) grid at the winglet.
(a)(b)


#### 2.1.3. Boundary Conditions and Calculation Parameters

To simulate the real scene of an aircraft flight, the fluid in the flow field is the ideal gas, and the fluid velocity is set to 72 m/s, which represents the approach speed of an aircraft. At the same time, considering the landing attitude of an aircraft, a 5° clamp between the fluid and the aircraft is set. The computational domain is a regular hexahedron configuration; the inlet is set as a velocity inlet; the outlet is set as a pressure outlet. The top surface and the left and right sides are set as symmetry planes. The bottom surface and the aircraft wing surface are set as nonslip solid wall surfaces, and the friction coefficient is 0.5. The fluid parameters are shown in Table2.Table 2
Fluid parameter table.
PressureTemperatureReynolds numberAerodynamic viscosityTurbulent intensityTurbulent kinetic energySpecific dissipation rate101.325 kPa288.15 K2.92e71.81e-51.86%2.70 m2/s28.76 1/s。In the discrete spatial method, the gradient adopts the least squares cell-based method. The pressure adopts the second order, and the momentum, energy, turbulent kinetic energy, and specific dissipation rate all adopt the second order upwind.


#### 2.1.4. Method Feasibility Verification

Reliable numerical simulations require a high degree of consistency between the predicted value and the experimental value. To verify the reliability of the numerical simulation results, this paper selects the NACA0012 airfoil, the SST turbulence model, and the incoming flow speed of 68 m/s. The above calculation method was used to calculate the static pressure coefficient at 72.5% of the wingspan, which was compared with the experiments of Chow et al. [14]. The comparison and verification results are shown in Figure 3.Figure 3
The pressure coefficient distribution at 72.5% span.As shown in Figure3, the simulation values are in good agreement with the experimental values, but the numerical simulation results reflect the pressure distribution near the airfoil in the flow field more accurately. The error of multiple calculations is within 2%. The flow structure is consistent with the reality is consistent with reality, so the numerical simulation method has high reliability in predicting the flow field around the wing.


#### 2.1.5. Wake Vortex Contour Results

The contour graphs of the wake vortex Q criterion at different positions on the trailing edge of the wing when an aircraft is at the height of 30 m from the ground is presented in Figure4. From the contours, the curling of the wake vortex and the change in the vortex structure, as well as the wake vortex intensity, can be observed.Figure 4
Distribution ofQ criterion at 0–10b0 behind the wing.
(a)
x/b0=0(b)
x/b0=2(c)
x/b0=4(d)
x/b0=6(e)
x/b0=8(f)
x/b0=10Figure4 shows that after the air flows through the wing surface, a pair of wingtip and flap vortices are formed. At the initial moment, the wingtip and flap vortices are two separate parts, and the wingtip vortex has a higher vortex degree. Over time, the flap vortices make the wingtip vortices roll up inward and merge with the flap vortices forming a pair of main vortices. Figure 4(d) shows that at the trailing edge 6 b0 of the wing, the two groups of vortex structures merge into the main vortex, and the position of the strongest vortex is from the tip vortex to the center of the main vortex. In addition, under the combined action of turbulent viscosity and mutual induction between the two vortices, the main vortex structure begins to spread outward and roll up, the vortex length gradually shortens, the vortex width gradually widens, and the vortex shape is rolled up from the vortex surface forming a ring-shape wake vortex.The vorticity distribution at different DHs is presented in Figure5. When the distance between the aircraft wake vortex and the ground is less than 1.5 times the initial vortex spacing, the wake vortex evolution process is significantly affected by the ground effect, and this is the near-earth evolution stage [15]. Compared to the cases of H=60m and H=30m, at H=10m and H=5m, the length of the flap vortex is significantly shortened, and the structural attenuation of the main vortex end is more obvious. A secondary vortex is induced near the ground, and it moves around the main vortex, making the main vortex unstable and decreasing the main vortex structure. Since the wake vortex induces a velocity field near the ground, it can be approximated to the pair of mirror vortices below the ground [16]. As shown in Figures 5(a) and 5(b), under the induction of mirror vortex, the closer the wake vortex is to the ground, the larger the distance between the vortex cores is. Due to the impenetrability of the ground, the wake bounces off the ground and increases in altitude.Figure 5
The vorticity distribution at different DHs.
(a)
H=5m(b)
H=10m(c)
H=30m(d)
H=60m


### 2.2. Far-Field Evolution of Wake Vortex

#### 2.2.1. Initial Conditions

A three-dimensional numerical simulation requires enormous computing resources. To study the long-term evolution process of the wake vortex at the DH, this study uses an ideal wake vortex tangential velocity model to initialize the convective field region and analyzes the stationary property relative to the ground and perpendicular to the flight path. For the evolution process of the wake vortex in cross section, the computational domain isLx×Ly=6.5b0×4.5b0, which is large enough to simulate the rebound and diffusion of the wake vortex under the influence of the ground effect. The grid size is a uniform grid size Δx=0.01b0 (0.5 m), and the bottom surface is a nonslip wall surface in the boundary conditions; the other three planes are symmetrical planes; the time step is 0.01 s. The remaining parameters and calculation methods are the same as in chapter 1. The Burnham-Hallock model is selected to initialize the wake vortex field as follows:
(2)vθ=Γ02πrr2r2+rc2,where Vθ is the tangential velocity of the wake vortex; Γ0 is the initial wake vortex circulation; and rc is the initial vortex core radius.The initial wake vortex circulation is calculated by(3)Γ0=MgρVb0,b0=π4B,rc=0.035B,where b0 is the initial vortex core distance; M is the maximum landing weight of an A330-200 aircraft, and it is 186 tons; B is the aircraft wingspan, and it is 60.3 m; ρ is the air density at the position of the wake vortex, which is 1.208 kg/m3; g is the local acceleration of gravity; and V is the A330-200 final approach speed, and it is 72 m/s.The characteristic speed, which is the initial descent speed of the wake vortex, is calculated by(4)ω0=Γ02πb0,t0=b0ω0,t∗=tt0,where t0 is the reference time, which represents the time required for the wake vortex to descend for distance b0 at the characteristic speed ω0, and t∗ is the dimensionless value of time. The wake vortex parameters are given in Table 3.Table 3
Partial parameters of the wake vortex.
Physical quantitiesValueΓ0m2/s442.75b0m47.36rcm2.11ω0m/s1.49t0s31.83The tangential velocity distribution of the tail vortex field initialized by the B-H model is presented in Figure6. The curve is the largest at the vortex center and gradually decreases along the radius.Figure 6
Tangential velocity distribution of the B-H model.To verify the authenticity of the numerical simulation, Figure7 compares radial velocity of the numerical simulation and that detected by radar under the same conditions.. The comparison shows that the two pairs of velocities have opposite directions, and the maximum radial velocity is approximately 5 m/s. Therefore, by compiling the trailing vortex with a tangential velocity model, the trailing vortex attenuation and its change in space position can be accurately simulated.Figure 7
Comparison of the initial state of the wake vortex obtained by the radar detection and numerical simulation. (a) Radial velocity of wake vortex detected by radar. (b) Numerical simulation of vertical velocity distribution of wake vortex.
(a)(b)


#### 2.2.2. Numerical Results

As shown in Figure8(a), due to the near-ground effect, after the wake vortex interacts with the ground, the ground will generate a boundary layer perpendicular to the axis of the main vortex. As the main vortex declines, a lateral pressure gradient is generated on the boundary layer, inducing the separation of the boundary layer from the ground (as shown in Figure 8(b)). Figure 8(c) shows that the separation layer continues to increase, forming a secondary vortex.Figure 8
Generation of the secondary vortex and rebound of the main vortex.
(a)
t∗ =0.1(b)
t∗ =0.19(c)
t∗ =0.4The rotation direction of the secondary vortex is opposite to that of the main vortex; it moves around the main vortex and reacts to the main vortex fast, causing the main vortex to deform and making the main vortex rebound and destabilize.Figure9 shows that different heights from the ground have different effects on the spatial position of the wake vortex. Due to the induction effect of the mirror vortex, the vortex spacing increases; the lower the distance from the ground is, the faster the vortex spacing will increase. The impermeability of the ground and the induction of secondary vortices cause the rebound. When H=60m or H=30m, the wake vortex will first descend and then rebound, and the trajectory of the wake vortex will be arc-shaped. At H=10m or H=5m, the height of the wake vortex will continue to increase, but the trajectory of the wake vortex will be linear.Figure 9
Changes in the wake vortex position at different heights from the ground. (a) The position of the vortex center varies with time; (b) the movement trajectory of the wake vortex center; and (c) the variation in vortex center spacing.
(a)(b)(c)Figure10 shows that different distances from the ground have different effects on the wake vortex intensity. The wake vortex intensity decay generally includes two stages, the diffusion stage (near-field vortex) and the fast decay stage (far-field vortex). In the diffusion stage, the wake vortex decay mainly depends on the radial diffusion of the vorticity, and the decay is relatively slow. When there is no ground effect, the wake vortex decays to 90% of its initial intensity [17]. After entering the fast decay stage, the wake vortex deforms, and the vortex intensity is fast. In the diffusion stage, due to the increase in the distance between the vortex cores, the two vortices develop into an isolated vortex stage, and the mutual-induced force between the two vortices weakens, preventing the development of crow instability and prolonging the life of the wake vortex [18]. When H=10m and H=5m, although the vortex spacing increases, the secondary vortex generated on the ground moves around the main vortex, thus making the main vortex rebound and destabilize, accelerating the dissipation of the wake vortex and shortening the time of the wake vortex diffusion stage. When H=5m, the wake vortex skips between the diffusion stages and enters the fast decay stage. Therefore, in the diffusion stage, the ground effect at the heights of 30 m and 60 m delays the dissipation of the wake vortex but accelerates the wake vortex decay at the heights of 5 m and 10 m. After entering the fast decay stage, the secondary vortex also becomes unstable, and the inductive effect on the main vortex weakens. However, when the vortex spacing increases, the isolated vortex stage is formed, and the induced force between the two vortices decreases gradually. The slowing effect of the vortex spacing increase on the wake vortex dissipation is stronger than the promotion effect of the secondary vortex on the dissipation of the wake vortex.Figure 10
Circulation changes at different distances from the ground.


## 2.1. Near-Field Evolution of Wake Vortex

### 2.1.1. Governing Equations and Turbulence Model

In this paper, the RANS method decomposes the turbulent motion into two parts: pulsation and average motion. By solving the time-averaged N-S equation of turbulence, the turbulent motion can be rapidly simulated. Compared with direct numerical simulation, the large-eddy simulation (LES) can save more computing resources. The SST model [11] is selected as a turbulence model. Compared with the standard model, this model has higher calculation accuracy and is suitable for rotating flow and other situations.


### 2.1.2. Computing Domains and Grids

The DH is the specific altitude at which a go-around should be initiated during a precision approach if the necessary visual reference to continue the approach cannot be established. Therefore, it is an important parameter of the precision approach stage and is a critical point for landing and overtaking. The study of the wake vortex evolution at the DH has important reference significance for preventing the following aircraft from entering the wake hazard zone and establishing the dynamic wake separation. The ILS (instrument landing system), which is considered a blind landing system, is the most widely used aircraft precision approach and landing guidance system. The ILSs can be roughly classified into three types based on the boot performance. The DH values of the ILS types are presented in Table1.Table 1
Classification and performance standards of the three ILS types.
Performance standardILS typeDH (m)IIIIIIABC6030Below 30 or 0Below 15 or 00Therefore, the wake vortices at different DHs corresponding to different types of precision approaches are selected for research in this study. The altitudes of an aircraft above the ground are set toH=60m, H=30m, H=10m, and H=5m. The aircraft model used in this study is A330-200, having a wingspan of 60.3 m, a wing area of 361.6 m2, and a fuselage length of 58.82 m. Since the horizontal tail vortex dissipates rapidly and has a slight effect on the safety of the following aircraft [12], only the wing part is retained in the calculation process. The computational domain meshing is completed using a structured mesh based on the finite volume method, and the wing surface is locally refined to meet the accuracy requirement of the calculation. The boundary layer mesh height y of the wing surface mesh is calculated by
(1)y=y+μVρ,where the wing surface y+ is one; μ is the dynamic viscosity of the turbulent flow; V is the fluid velocity; and ρ is the fluid density.The computational domain with an aircraft as the origin is presented in Figure1(a), where the x-direction denotes the approach direction of the aircraft, the y-direction is the wingspan direction, and the z-direction is the height direction. The wing is 3 b0 from the entrance, 2 b0 from the top surface, and 10 b0 from the exit, thus having enough spatial scales to distinguish between the long and short wave instability of crow [13]. By setting different distances between an aircraft and the lower wall, the ground effect at different heights can be studied. The structured grid of the fluid domain is presented in Figure 1(b). To capture the evolution and development of the wake vortex accurately, special grid refinement is performed in the wake region to ensure calculation accuracy. In addition, the mesh of the vortex area should have a good transition to ensure that the airflow boundary layer is accurately resolved. Therefore, the aspect ratio of the grid in the wake region is close to one, and the minimum grid size in the wake vortex area is 0.1m×0.1m.Figure 1
Computational domain and its grid. (a) Schematic diagram of computational domain and (b) computational domain grid.
(a)(b)The leading edge mesh of the wing is presented in Figure2(a). According to Equation (5), the boundary layer is densified to satisfy the condition of y+=1; the height of the first layer of the boundary layer is 6.14e-6 m, and the boundary layer has 40 layers. Figure 2(b) shows the surface mesh of the winglet. The surface mesh of the wingtip is refined to capture the evolution of the wingtip vortex better.Figure 2
The A330-200 wing partial grid. (a) Grid at the leading edge of the wing and (b) grid at the winglet.
(a)(b)


### 2.1.3. Boundary Conditions and Calculation Parameters

To simulate the real scene of an aircraft flight, the fluid in the flow field is the ideal gas, and the fluid velocity is set to 72 m/s, which represents the approach speed of an aircraft. At the same time, considering the landing attitude of an aircraft, a 5° clamp between the fluid and the aircraft is set. The computational domain is a regular hexahedron configuration; the inlet is set as a velocity inlet; the outlet is set as a pressure outlet. The top surface and the left and right sides are set as symmetry planes. The bottom surface and the aircraft wing surface are set as nonslip solid wall surfaces, and the friction coefficient is 0.5. The fluid parameters are shown in Table2.Table 2
Fluid parameter table.
PressureTemperatureReynolds numberAerodynamic viscosityTurbulent intensityTurbulent kinetic energySpecific dissipation rate101.325 kPa288.15 K2.92e71.81e-51.86%2.70 m2/s28.76 1/s。In the discrete spatial method, the gradient adopts the least squares cell-based method. The pressure adopts the second order, and the momentum, energy, turbulent kinetic energy, and specific dissipation rate all adopt the second order upwind.


### 2.1.4. Method Feasibility Verification

Reliable numerical simulations require a high degree of consistency between the predicted value and the experimental value. To verify the reliability of the numerical simulation results, this paper selects the NACA0012 airfoil, the SST turbulence model, and the incoming flow speed of 68 m/s. The above calculation method was used to calculate the static pressure coefficient at 72.5% of the wingspan, which was compared with the experiments of Chow et al. [14]. The comparison and verification results are shown in Figure 3.Figure 3
The pressure coefficient distribution at 72.5% span.As shown in Figure3, the simulation values are in good agreement with the experimental values, but the numerical simulation results reflect the pressure distribution near the airfoil in the flow field more accurately. The error of multiple calculations is within 2%. The flow structure is consistent with the reality is consistent with reality, so the numerical simulation method has high reliability in predicting the flow field around the wing.


### 2.1.5. Wake Vortex Contour Results

The contour graphs of the wake vortex Q criterion at different positions on the trailing edge of the wing when an aircraft is at the height of 30 m from the ground is presented in Figure4. From the contours, the curling of the wake vortex and the change in the vortex structure, as well as the wake vortex intensity, can be observed.Figure 4
Distribution ofQ criterion at 0–10b0 behind the wing.
(a)
x/b0=0(b)
x/b0=2(c)
x/b0=4(d)
x/b0=6(e)
x/b0=8(f)
x/b0=10Figure4 shows that after the air flows through the wing surface, a pair of wingtip and flap vortices are formed. At the initial moment, the wingtip and flap vortices are two separate parts, and the wingtip vortex has a higher vortex degree. Over time, the flap vortices make the wingtip vortices roll up inward and merge with the flap vortices forming a pair of main vortices. Figure 4(d) shows that at the trailing edge 6 b0 of the wing, the two groups of vortex structures merge into the main vortex, and the position of the strongest vortex is from the tip vortex to the center of the main vortex. In addition, under the combined action of turbulent viscosity and mutual induction between the two vortices, the main vortex structure begins to spread outward and roll up, the vortex length gradually shortens, the vortex width gradually widens, and the vortex shape is rolled up from the vortex surface forming a ring-shape wake vortex.The vorticity distribution at different DHs is presented in Figure5. When the distance between the aircraft wake vortex and the ground is less than 1.5 times the initial vortex spacing, the wake vortex evolution process is significantly affected by the ground effect, and this is the near-earth evolution stage [15]. Compared to the cases of H=60m and H=30m, at H=10m and H=5m, the length of the flap vortex is significantly shortened, and the structural attenuation of the main vortex end is more obvious. A secondary vortex is induced near the ground, and it moves around the main vortex, making the main vortex unstable and decreasing the main vortex structure. Since the wake vortex induces a velocity field near the ground, it can be approximated to the pair of mirror vortices below the ground [16]. As shown in Figures 5(a) and 5(b), under the induction of mirror vortex, the closer the wake vortex is to the ground, the larger the distance between the vortex cores is. Due to the impenetrability of the ground, the wake bounces off the ground and increases in altitude.Figure 5
The vorticity distribution at different DHs.
(a)
H=5m(b)
H=10m(c)
H=30m(d)
H=60m


## 2.1.1. Governing Equations and Turbulence Model

In this paper, the RANS method decomposes the turbulent motion into two parts: pulsation and average motion. By solving the time-averaged N-S equation of turbulence, the turbulent motion can be rapidly simulated. Compared with direct numerical simulation, the large-eddy simulation (LES) can save more computing resources. The SST model [11] is selected as a turbulence model. Compared with the standard model, this model has higher calculation accuracy and is suitable for rotating flow and other situations.


## 2.1.2. Computing Domains and Grids

The DH is the specific altitude at which a go-around should be initiated during a precision approach if the necessary visual reference to continue the approach cannot be established. Therefore, it is an important parameter of the precision approach stage and is a critical point for landing and overtaking. The study of the wake vortex evolution at the DH has important reference significance for preventing the following aircraft from entering the wake hazard zone and establishing the dynamic wake separation. The ILS (instrument landing system), which is considered a blind landing system, is the most widely used aircraft precision approach and landing guidance system. The ILSs can be roughly classified into three types based on the boot performance. The DH values of the ILS types are presented in Table1.Table 1
Classification and performance standards of the three ILS types.
Performance standardILS typeDH (m)IIIIIIABC6030Below 30 or 0Below 15 or 00Therefore, the wake vortices at different DHs corresponding to different types of precision approaches are selected for research in this study. The altitudes of an aircraft above the ground are set toH=60m, H=30m, H=10m, and H=5m. The aircraft model used in this study is A330-200, having a wingspan of 60.3 m, a wing area of 361.6 m2, and a fuselage length of 58.82 m. Since the horizontal tail vortex dissipates rapidly and has a slight effect on the safety of the following aircraft [12], only the wing part is retained in the calculation process. The computational domain meshing is completed using a structured mesh based on the finite volume method, and the wing surface is locally refined to meet the accuracy requirement of the calculation. The boundary layer mesh height y of the wing surface mesh is calculated by
(1)y=y+μVρ,where the wing surface y+ is one; μ is the dynamic viscosity of the turbulent flow; V is the fluid velocity; and ρ is the fluid density.The computational domain with an aircraft as the origin is presented in Figure1(a), where the x-direction denotes the approach direction of the aircraft, the y-direction is the wingspan direction, and the z-direction is the height direction. The wing is 3 b0 from the entrance, 2 b0 from the top surface, and 10 b0 from the exit, thus having enough spatial scales to distinguish between the long and short wave instability of crow [13]. By setting different distances between an aircraft and the lower wall, the ground effect at different heights can be studied. The structured grid of the fluid domain is presented in Figure 1(b). To capture the evolution and development of the wake vortex accurately, special grid refinement is performed in the wake region to ensure calculation accuracy. In addition, the mesh of the vortex area should have a good transition to ensure that the airflow boundary layer is accurately resolved. Therefore, the aspect ratio of the grid in the wake region is close to one, and the minimum grid size in the wake vortex area is 0.1m×0.1m.Figure 1
Computational domain and its grid. (a) Schematic diagram of computational domain and (b) computational domain grid.
(a)(b)The leading edge mesh of the wing is presented in Figure2(a). According to Equation (5), the boundary layer is densified to satisfy the condition of y+=1; the height of the first layer of the boundary layer is 6.14e-6 m, and the boundary layer has 40 layers. Figure 2(b) shows the surface mesh of the winglet. The surface mesh of the wingtip is refined to capture the evolution of the wingtip vortex better.Figure 2
The A330-200 wing partial grid. (a) Grid at the leading edge of the wing and (b) grid at the winglet.
(a)(b)


## 2.1.3. Boundary Conditions and Calculation Parameters

To simulate the real scene of an aircraft flight, the fluid in the flow field is the ideal gas, and the fluid velocity is set to 72 m/s, which represents the approach speed of an aircraft. At the same time, considering the landing attitude of an aircraft, a 5° clamp between the fluid and the aircraft is set. The computational domain is a regular hexahedron configuration; the inlet is set as a velocity inlet; the outlet is set as a pressure outlet. The top surface and the left and right sides are set as symmetry planes. The bottom surface and the aircraft wing surface are set as nonslip solid wall surfaces, and the friction coefficient is 0.5. The fluid parameters are shown in Table2.Table 2
Fluid parameter table.
PressureTemperatureReynolds numberAerodynamic viscosityTurbulent intensityTurbulent kinetic energySpecific dissipation rate101.325 kPa288.15 K2.92e71.81e-51.86%2.70 m2/s28.76 1/s。In the discrete spatial method, the gradient adopts the least squares cell-based method. The pressure adopts the second order, and the momentum, energy, turbulent kinetic energy, and specific dissipation rate all adopt the second order upwind.


## 2.1.4. Method Feasibility Verification

Reliable numerical simulations require a high degree of consistency between the predicted value and the experimental value. To verify the reliability of the numerical simulation results, this paper selects the NACA0012 airfoil, the SST turbulence model, and the incoming flow speed of 68 m/s. The above calculation method was used to calculate the static pressure coefficient at 72.5% of the wingspan, which was compared with the experiments of Chow et al. [14]. The comparison and verification results are shown in Figure 3.Figure 3
The pressure coefficient distribution at 72.5% span.As shown in Figure3, the simulation values are in good agreement with the experimental values, but the numerical simulation results reflect the pressure distribution near the airfoil in the flow field more accurately. The error of multiple calculations is within 2%. The flow structure is consistent with the reality is consistent with reality, so the numerical simulation method has high reliability in predicting the flow field around the wing.


## 2.1.5. Wake Vortex Contour Results

The contour graphs of the wake vortex Q criterion at different positions on the trailing edge of the wing when an aircraft is at the height of 30 m from the ground is presented in Figure4. From the contours, the curling of the wake vortex and the change in the vortex structure, as well as the wake vortex intensity, can be observed.Figure 4
Distribution ofQ criterion at 0–10b0 behind the wing.
(a)
x/b0=0(b)
x/b0=2(c)
x/b0=4(d)
x/b0=6(e)
x/b0=8(f)
x/b0=10Figure4 shows that after the air flows through the wing surface, a pair of wingtip and flap vortices are formed. At the initial moment, the wingtip and flap vortices are two separate parts, and the wingtip vortex has a higher vortex degree. Over time, the flap vortices make the wingtip vortices roll up inward and merge with the flap vortices forming a pair of main vortices. Figure 4(d) shows that at the trailing edge 6 b0 of the wing, the two groups of vortex structures merge into the main vortex, and the position of the strongest vortex is from the tip vortex to the center of the main vortex. In addition, under the combined action of turbulent viscosity and mutual induction between the two vortices, the main vortex structure begins to spread outward and roll up, the vortex length gradually shortens, the vortex width gradually widens, and the vortex shape is rolled up from the vortex surface forming a ring-shape wake vortex.The vorticity distribution at different DHs is presented in Figure5. When the distance between the aircraft wake vortex and the ground is less than 1.5 times the initial vortex spacing, the wake vortex evolution process is significantly affected by the ground effect, and this is the near-earth evolution stage [15]. Compared to the cases of H=60m and H=30m, at H=10m and H=5m, the length of the flap vortex is significantly shortened, and the structural attenuation of the main vortex end is more obvious. A secondary vortex is induced near the ground, and it moves around the main vortex, making the main vortex unstable and decreasing the main vortex structure. Since the wake vortex induces a velocity field near the ground, it can be approximated to the pair of mirror vortices below the ground [16]. As shown in Figures 5(a) and 5(b), under the induction of mirror vortex, the closer the wake vortex is to the ground, the larger the distance between the vortex cores is. Due to the impenetrability of the ground, the wake bounces off the ground and increases in altitude.Figure 5
The vorticity distribution at different DHs.
(a)
H=5m(b)
H=10m(c)
H=30m(d)
H=60m


## 2.2. Far-Field Evolution of Wake Vortex

### 2.2.1. Initial Conditions

A three-dimensional numerical simulation requires enormous computing resources. To study the long-term evolution process of the wake vortex at the DH, this study uses an ideal wake vortex tangential velocity model to initialize the convective field region and analyzes the stationary property relative to the ground and perpendicular to the flight path. For the evolution process of the wake vortex in cross section, the computational domain isLx×Ly=6.5b0×4.5b0, which is large enough to simulate the rebound and diffusion of the wake vortex under the influence of the ground effect. The grid size is a uniform grid size Δx=0.01b0 (0.5 m), and the bottom surface is a nonslip wall surface in the boundary conditions; the other three planes are symmetrical planes; the time step is 0.01 s. The remaining parameters and calculation methods are the same as in chapter 1. The Burnham-Hallock model is selected to initialize the wake vortex field as follows:
(2)vθ=Γ02πrr2r2+rc2,where Vθ is the tangential velocity of the wake vortex; Γ0 is the initial wake vortex circulation; and rc is the initial vortex core radius.The initial wake vortex circulation is calculated by(3)Γ0=MgρVb0,b0=π4B,rc=0.035B,where b0 is the initial vortex core distance; M is the maximum landing weight of an A330-200 aircraft, and it is 186 tons; B is the aircraft wingspan, and it is 60.3 m; ρ is the air density at the position of the wake vortex, which is 1.208 kg/m3; g is the local acceleration of gravity; and V is the A330-200 final approach speed, and it is 72 m/s.The characteristic speed, which is the initial descent speed of the wake vortex, is calculated by(4)ω0=Γ02πb0,t0=b0ω0,t∗=tt0,where t0 is the reference time, which represents the time required for the wake vortex to descend for distance b0 at the characteristic speed ω0, and t∗ is the dimensionless value of time. The wake vortex parameters are given in Table 3.Table 3
Partial parameters of the wake vortex.
Physical quantitiesValueΓ0m2/s442.75b0m47.36rcm2.11ω0m/s1.49t0s31.83The tangential velocity distribution of the tail vortex field initialized by the B-H model is presented in Figure6. The curve is the largest at the vortex center and gradually decreases along the radius.Figure 6
Tangential velocity distribution of the B-H model.To verify the authenticity of the numerical simulation, Figure7 compares radial velocity of the numerical simulation and that detected by radar under the same conditions.. The comparison shows that the two pairs of velocities have opposite directions, and the maximum radial velocity is approximately 5 m/s. Therefore, by compiling the trailing vortex with a tangential velocity model, the trailing vortex attenuation and its change in space position can be accurately simulated.Figure 7
Comparison of the initial state of the wake vortex obtained by the radar detection and numerical simulation. (a) Radial velocity of wake vortex detected by radar. (b) Numerical simulation of vertical velocity distribution of wake vortex.
(a)(b)


### 2.2.2. Numerical Results

As shown in Figure8(a), due to the near-ground effect, after the wake vortex interacts with the ground, the ground will generate a boundary layer perpendicular to the axis of the main vortex. As the main vortex declines, a lateral pressure gradient is generated on the boundary layer, inducing the separation of the boundary layer from the ground (as shown in Figure 8(b)). Figure 8(c) shows that the separation layer continues to increase, forming a secondary vortex.Figure 8
Generation of the secondary vortex and rebound of the main vortex.
(a)
t∗ =0.1(b)
t∗ =0.19(c)
t∗ =0.4The rotation direction of the secondary vortex is opposite to that of the main vortex; it moves around the main vortex and reacts to the main vortex fast, causing the main vortex to deform and making the main vortex rebound and destabilize.Figure9 shows that different heights from the ground have different effects on the spatial position of the wake vortex. Due to the induction effect of the mirror vortex, the vortex spacing increases; the lower the distance from the ground is, the faster the vortex spacing will increase. The impermeability of the ground and the induction of secondary vortices cause the rebound. When H=60m or H=30m, the wake vortex will first descend and then rebound, and the trajectory of the wake vortex will be arc-shaped. At H=10m or H=5m, the height of the wake vortex will continue to increase, but the trajectory of the wake vortex will be linear.Figure 9
Changes in the wake vortex position at different heights from the ground. (a) The position of the vortex center varies with time; (b) the movement trajectory of the wake vortex center; and (c) the variation in vortex center spacing.
(a)(b)(c)Figure10 shows that different distances from the ground have different effects on the wake vortex intensity. The wake vortex intensity decay generally includes two stages, the diffusion stage (near-field vortex) and the fast decay stage (far-field vortex). In the diffusion stage, the wake vortex decay mainly depends on the radial diffusion of the vorticity, and the decay is relatively slow. When there is no ground effect, the wake vortex decays to 90% of its initial intensity [17]. After entering the fast decay stage, the wake vortex deforms, and the vortex intensity is fast. In the diffusion stage, due to the increase in the distance between the vortex cores, the two vortices develop into an isolated vortex stage, and the mutual-induced force between the two vortices weakens, preventing the development of crow instability and prolonging the life of the wake vortex [18]. When H=10m and H=5m, although the vortex spacing increases, the secondary vortex generated on the ground moves around the main vortex, thus making the main vortex rebound and destabilize, accelerating the dissipation of the wake vortex and shortening the time of the wake vortex diffusion stage. When H=5m, the wake vortex skips between the diffusion stages and enters the fast decay stage. Therefore, in the diffusion stage, the ground effect at the heights of 30 m and 60 m delays the dissipation of the wake vortex but accelerates the wake vortex decay at the heights of 5 m and 10 m. After entering the fast decay stage, the secondary vortex also becomes unstable, and the inductive effect on the main vortex weakens. However, when the vortex spacing increases, the isolated vortex stage is formed, and the induced force between the two vortices decreases gradually. The slowing effect of the vortex spacing increase on the wake vortex dissipation is stronger than the promotion effect of the secondary vortex on the dissipation of the wake vortex.Figure 10
Circulation changes at different distances from the ground.


## 2.2.1. Initial Conditions

A three-dimensional numerical simulation requires enormous computing resources. To study the long-term evolution process of the wake vortex at the DH, this study uses an ideal wake vortex tangential velocity model to initialize the convective field region and analyzes the stationary property relative to the ground and perpendicular to the flight path. For the evolution process of the wake vortex in cross section, the computational domain isLx×Ly=6.5b0×4.5b0, which is large enough to simulate the rebound and diffusion of the wake vortex under the influence of the ground effect. The grid size is a uniform grid size Δx=0.01b0 (0.5 m), and the bottom surface is a nonslip wall surface in the boundary conditions; the other three planes are symmetrical planes; the time step is 0.01 s. The remaining parameters and calculation methods are the same as in chapter 1. The Burnham-Hallock model is selected to initialize the wake vortex field as follows:
(2)vθ=Γ02πrr2r2+rc2,where Vθ is the tangential velocity of the wake vortex; Γ0 is the initial wake vortex circulation; and rc is the initial vortex core radius.The initial wake vortex circulation is calculated by(3)Γ0=MgρVb0,b0=π4B,rc=0.035B,where b0 is the initial vortex core distance; M is the maximum landing weight of an A330-200 aircraft, and it is 186 tons; B is the aircraft wingspan, and it is 60.3 m; ρ is the air density at the position of the wake vortex, which is 1.208 kg/m3; g is the local acceleration of gravity; and V is the A330-200 final approach speed, and it is 72 m/s.The characteristic speed, which is the initial descent speed of the wake vortex, is calculated by(4)ω0=Γ02πb0,t0=b0ω0,t∗=tt0,where t0 is the reference time, which represents the time required for the wake vortex to descend for distance b0 at the characteristic speed ω0, and t∗ is the dimensionless value of time. The wake vortex parameters are given in Table 3.Table 3
Partial parameters of the wake vortex.
Physical quantitiesValueΓ0m2/s442.75b0m47.36rcm2.11ω0m/s1.49t0s31.83The tangential velocity distribution of the tail vortex field initialized by the B-H model is presented in Figure6. The curve is the largest at the vortex center and gradually decreases along the radius.Figure 6
Tangential velocity distribution of the B-H model.To verify the authenticity of the numerical simulation, Figure7 compares radial velocity of the numerical simulation and that detected by radar under the same conditions.. The comparison shows that the two pairs of velocities have opposite directions, and the maximum radial velocity is approximately 5 m/s. Therefore, by compiling the trailing vortex with a tangential velocity model, the trailing vortex attenuation and its change in space position can be accurately simulated.Figure 7
Comparison of the initial state of the wake vortex obtained by the radar detection and numerical simulation. (a) Radial velocity of wake vortex detected by radar. (b) Numerical simulation of vertical velocity distribution of wake vortex.
(a)(b)


## 2.2.2. Numerical Results

As shown in Figure8(a), due to the near-ground effect, after the wake vortex interacts with the ground, the ground will generate a boundary layer perpendicular to the axis of the main vortex. As the main vortex declines, a lateral pressure gradient is generated on the boundary layer, inducing the separation of the boundary layer from the ground (as shown in Figure 8(b)). Figure 8(c) shows that the separation layer continues to increase, forming a secondary vortex.Figure 8
Generation of the secondary vortex and rebound of the main vortex.
(a)
t∗ =0.1(b)
t∗ =0.19(c)
t∗ =0.4The rotation direction of the secondary vortex is opposite to that of the main vortex; it moves around the main vortex and reacts to the main vortex fast, causing the main vortex to deform and making the main vortex rebound and destabilize.Figure9 shows that different heights from the ground have different effects on the spatial position of the wake vortex. Due to the induction effect of the mirror vortex, the vortex spacing increases; the lower the distance from the ground is, the faster the vortex spacing will increase. The impermeability of the ground and the induction of secondary vortices cause the rebound. When H=60m or H=30m, the wake vortex will first descend and then rebound, and the trajectory of the wake vortex will be arc-shaped. At H=10m or H=5m, the height of the wake vortex will continue to increase, but the trajectory of the wake vortex will be linear.Figure 9
Changes in the wake vortex position at different heights from the ground. (a) The position of the vortex center varies with time; (b) the movement trajectory of the wake vortex center; and (c) the variation in vortex center spacing.
(a)(b)(c)Figure10 shows that different distances from the ground have different effects on the wake vortex intensity. The wake vortex intensity decay generally includes two stages, the diffusion stage (near-field vortex) and the fast decay stage (far-field vortex). In the diffusion stage, the wake vortex decay mainly depends on the radial diffusion of the vorticity, and the decay is relatively slow. When there is no ground effect, the wake vortex decays to 90% of its initial intensity [17]. After entering the fast decay stage, the wake vortex deforms, and the vortex intensity is fast. In the diffusion stage, due to the increase in the distance between the vortex cores, the two vortices develop into an isolated vortex stage, and the mutual-induced force between the two vortices weakens, preventing the development of crow instability and prolonging the life of the wake vortex [18]. When H=10m and H=5m, although the vortex spacing increases, the secondary vortex generated on the ground moves around the main vortex, thus making the main vortex rebound and destabilize, accelerating the dissipation of the wake vortex and shortening the time of the wake vortex diffusion stage. When H=5m, the wake vortex skips between the diffusion stages and enters the fast decay stage. Therefore, in the diffusion stage, the ground effect at the heights of 30 m and 60 m delays the dissipation of the wake vortex but accelerates the wake vortex decay at the heights of 5 m and 10 m. After entering the fast decay stage, the secondary vortex also becomes unstable, and the inductive effect on the main vortex weakens. However, when the vortex spacing increases, the isolated vortex stage is formed, and the induced force between the two vortices decreases gradually. The slowing effect of the vortex spacing increase on the wake vortex dissipation is stronger than the promotion effect of the secondary vortex on the dissipation of the wake vortex.Figure 10
Circulation changes at different distances from the ground.


## 3. ARJ21 Wake Encounter Model

When the following aircraft enters the wake area of the front aircraft, the induced velocity field of the wake will cause uneven force on the aircraft’s body, and the following aircraft will be subjected to the rolling torque and overload of different degrees. The rolling moment will make the aircraft roll, and the upwash and downwash forces will cause aircraft turbulence and damage to the aircraft’s body structure. In the landing stage, the aircraft is at a lower height from the ground, and the downwash force will decrease the aircraft’s height, resulting in a hard landing. The above situations are all dangerous because the following aircraft encounters the wake of the front aircraft. Therefore, this study analyzes the high accident rate of a medium-size aircraft following the approach and landing of a heavy aircraft and uses the aerodynamic response model based on the strip method for modeling. Figure11 shows two different modes of the following aircraft encountering the front aircraft wake. In Figure 11(a), the tail vortices of the following aircraft and the front aircraft are at the same height. However, usually due to the turbulent viscous force and mutual induction between the two vortices, the tail vortex of the front engine will sink, as shown in Figure 11(b). When the following aircraft and the front aircraft are not at the same level, the danger will be greatly reduced.Figure 11
Two scenarios of the following aircraft encountering the wake of the front aircraft.
(a)(b)As shown in Figure12(a), the coordinate system is constructed using the midpoint of the two vortices as the coordinate system origin. The aircraft nose position (x, y) is defined as the aircraft position. The vertical induced velocity of the wake vortex is displayed in Figure 12(b), and the vertical velocity near the center of the vortex core is the largest. The induced velocity decreases with the position. Therefore, the position of the aircraft relative to the wake vortex is different, and its aerodynamic force is also different., and the downwash force generated by the induced velocity of the front wake vortex at each point on the aircraft is also different. According to the Hallock-Burnham model, the vertical velocity of any point in space can be calculated by
(5)VZ=VZ1+VZ2=Γ1y−y12πy−y12+z−z12+rc2+Γ2y−y22πy−y22+z−z22+rc2,where VZ is the induced resultant velocity of the wake in the vertical direction; Γ1 is the left vortex circulation; Γ2 is the right vortex circulation; and rc is the radius of the vortex core.Figure 12
The force of the ARJ21 aircraft encountering the wake. (a) Two-dimensional plane coordinate system of the wake encountered by the following aircraft. (b) The vertical velocity distribution of the wake vortex of the B-H model in the spanwise direction.
(a)(b)To calculate the aerodynamic force of the following aircraft precisely, this study divides the following aircraft into four parts, namely, wing, fuselage, engine, and tail, according to the geometric configuration of the aircraft’s body and designs a strip algorithm. As shown in Figure13, the wing, fuselage, engine, and tail are divided into countless strips. The strip model is applied to the velocity field of the tail vortex of the front aircraft. The Runge-Kutta integral algorithm is used to calculate the aerodynamic force magnitude, and the moment is accumulated. Then, the rolling moment coefficient and overload increment of the aircraft are obtained at different positions to judge the safety of the aircraft at these positions.Figure 13
Schematic diagram of the strip model.When an aircraft enters the induced velocity field formed by the front wake vortex field, the lift of the aircraft will change. Due to the different geometric parameters of different parts of the aircraft’s body, the corresponding lift calculation methods also differ. The additional lift variation of the wing or tail caused by the wake vortex field can be obtained by(6)ΔF=12ρ∞V2∫−B/2B/2ΔCLcydy,where ρ∞ is the atmospheric density; V is the incoming flow velocity (i.e., the aircraft vacuum speed); ΔCL is the change in lift coefficient; and cy is the wing chord at the wing or tail wing wingspan coordinate. The lift coefficient change is calculated by
(7)ΔCL=CLαarctanVzyV≈CLαVzyVwhere Vzy is the induced velocity of the front tail vortex field on the following wing or tail profile; arctanVzy/V is the change in wing or tail section angle of attack of the following aircraft caused by the front aircraft wake vortex; the angle of attack is typically small, so it can be approximated as Vzy/V; and CLα is the lift line slope.For the commonly used swept-wing aircraft, the wing or tail section wing’s chord length can be approximated to [19]
(8)cy=cr1−21−λyB=2SB+Bλ−2y1−λB21+λ,where cr is the chord length at the wing root; λ is the tip-to-root ratio; and S is the wing or tail area.The moment can be obtained byLR=L×F, so the induced rolling force of the wing and tail is expressed by
(9)LR=12ρ∞VCLα∫−B/2B/2VZycyydy,where is the induced rolling moment and y is the distance between the strips of the wing and tail and the symmetry axis of the aircraft’s body.The fuselage can be considered a slender cylinder with a small angle of attack. According to the slender spin linearization theory, its lift change can be expressed by [20]
(10)ΔFbody=Ncosα−Asinα,where ΔFL is the fuselage-induced lift change; N is the normal force; A is the axial force; and α is the angle of attack.For a slender body with a small angle of attack, the normal and axial forces can be obtained using the potential flow theory as follows:(11)N=12ρ∞V∞2Slr2α,A=−12ρ∞V∞2Slrα2,α=arctanVzyV∞≈VzyV∞,where Slr is the integral area of the fuselage strip.The rolling moment of the fuselage is given by(12)LR′=ρ∞V∞2∫‐D/2D/2VZyV∞2+VZ2ylyydy+12ρ∞∫‐D/2D/2VZ3yV∞2+VZ2ylyydy.The lift variation, roll moment, and engine roll moment coefficient can be obtained using the vortex plate numerical method as follows [21]:
(13)ΔFengine=ρ∞V∑j=1nVjSj=ρ∞V∫−b/2b/2Vzysydy,(14)L′R′=ρ∞V∫−b/2b/2Vzysyydy,where Vj is the induced speed of the front vortex in the engine of the following airplane.Since in aircraft engines are relatively close to the wing, they can be obtained based on the induced velocity on the wing along they-axis direction. In Equations (13) and (14), Sj denotes the strip length; b is the engine span length; y is the distance from a strip on the engine to its centerline; and S′ is the engine area.The roll moment coefficient and the overload increment indicate the danger degree of an aircraft encountering a wake vortex. According to Steven Lang’s experimental results [22], the overall rolling moment coefficient of an aircraft is 0.05–0.07, which is the maximum value that the aircraft’s roll control authority can use in the aileron control. If this safety threshold is exceeded, an aircraft will lose both stability and control. Generally, the critical value of the rolling moment coefficient is set to 0.05, and the rolling moment coefficient is calculated by
(15)RMC=2LRρ∞V2SB.The overload increment is a standard measure of the aircraft turbulence intensity, and it is generally considered that the critical overload safety is moderate overload [23], with a value of 0.5. The overload increment is obtained as follows:
(16)Δn=ΔFallMg,where ΔFall is the external force exerted on the aircraft in the vertical direction.


## 4. ARJ21 Encounter Wake Risk Analysis

### 4.1. Risk Analysis of ARJ21 Encountering Wake in the Approach Section

To verify the safety of the current wake separation, a medium-size aircraft ARJ21 is selected to follow a heavy aircraft A330-200 to approach and land at the 5-nmile wake separation specified by the ICAO.When the ARJ21 and the front vortex field are at the same height (as shown in Figure11(a)), Figure 14 shows the relationship between the rolling moment coefficient of different parts and the change in overload increment as the ARJ21 enters the spanwise position of the front wake vortex field. Positive and negative values of the rolling torque coefficient represent counterclockwise and clockwise rolling types, respectively. Also, positive and negative overload values denote the upwash the downwash overloads, respectively. The absolute value indicates the degree of rolling and turbulence. The torque coefficient increases significantly near the vortex core, reaching the maximum at the center of the vortex core. Also, the rolling torque coefficient of the whole aircraft is 0.027, which is less than the threshold of 0.05. Different from the maximum value of the tail rolling moment coefficient of each part of the aircraft’s body, the maximum overload position of the wing is 10.6 m from both sides of the midpoint; the maximum overload position of the whole machine is 21 m from both sides of the midpoint. The maximum of RMC is 0.052, which is much smaller than the threshold value of 0.5. The rolling moment coefficient and the overload increment are both smaller than the corresponding critical values, which verifies the safety of the current wake separation.Figure 14
The coefficient of each part of the ARJ21 at different spanwise positions of the wake vortex field at the same height.
(a)(b)Because of the front vortex sinking, the ARJ21 and the front vortex field are typically at different heights, as shown in Figure11(b). The relationship between the rolling moment coefficient and the overload increment of the ARJ21 with the spanwise position of the ARJ21 entering the wake vortex field at different heights from the front vortex center is displayed in Figure 15. The spanwise position corresponding to the maximum rolling moment coefficient does not change with the height. In contrast, the spanwise position corresponding to the maximum overload gradually approaches the midpoint with the height. As the height increases, both the rolling moment coefficient and the overload increment decrease rapidly, and the rolling moment coefficient decreases significantly. This indicates that the judgment on whether the following aircraft is safe or not depends not only on the wake intensity of the front aircraft but also on the position of the following aircraft entering the wake vortex field. The existing wake separation still has a large space for reduction.Figure 15
The overall coefficient of the ARJ21 with the spanwise position of the wake vortex field at different heights.
(a)(b)To shorten the wake separation further, it is needed to delineate the wake danger zone of the front aircraft according to the movement trajectory and intensity attenuation changes of the front aircraft. Then, it is superimposed with the approach course of the following aircraft, and the length of the overlapping part is set as a minimum wake separation. The procedure of the final phase of the precision approach is illustrated in Figure16(a). When the front aircraft is along the precision approach glide path to the DH, the wake behind the wing will affect the aircraft approaching along the same glide path. Based on the wake encounter and aerodynamic response model, the wake vortex at the DH is used as a reference frame to analyze the force of the following aircraft. As shown in Figure 16(b), when the approach direction is the x-axis direction, the spanwise direction is the y-direction, and the vertical direction is the z-axis, the intersection of the section where the DH is located, and the extension line of the runway centerline is the origin. The longitudinal position of the initial wake vortex at the DH is set as the x-axis zero point, the center point of the two vortices is the y-axis zero point, and the airport surface height is the z-axis zero point. Further, the position where the nose of the following aircraft enters the wake field is denoted by (x, y, z). According to the wake movement trajectory of the front aircraft and its influence on the following aircraft, the areas with different wake risk levels are defined.Figure 16
The ARJ21 encounter wake coordinate system in the final approach segment. (a) Precision approach procedure and (b) the three-dimensional coordinate system of the wake encounter.
(a)(b)According to the maximum roll control power that can be controlled by the ailerons of an aircraft, the roll risk is categorized into three levels according to the roll moment coefficient. The roll risk of 0.025–0.05 is a mild danger risk. Namely, when aircraft rolls slightly, it is easier for a pilot to recover the aircraft. The roll risk of 0.05–0.08 is a moderate-danger risk, and in this case, the rolling slope angle increases, and the risk of stalling arises. Therefore, a pilot needs to put in larger efforts to recover and maintain the original flight path. When the roll risk is higher than 0.08, it is a severe danger risk. In this case, the aircraft’s body rolls and shakes obviously; the aircraft is facing a stall; and it is difficult for a pilot to recover.As illustrated in [24], the risk range of 0.15–0.5 corresponds to a mild turbulence area. In this case, the aircraft has a slight altitude change; brief yaw, pitch, and roll attitudes occur but have minimal effect. Since the aircraft overload is difficult to be larger than one during the approach phase, the risk range of 0.5–0.75 is defined as a moderate turbulence area in this study. In this area, the flight altitude and attitude of aircraft change significantly, but the aircraft is still in a controllable state. The risk of above 0.75 corresponds to a severe turbulence area, where flight attitude and altitude change drastically. In this area, the airspeed indicated by the instrument changes significantly, and the aircraft loses control for a short time. In severe cases, damage to the body structure can be caused, resulting in an aviation accident. According to the above hazard level categorization method, the wake hazard area on the trailing edge of the front wing can be visualized, boundaries of different hazard levels can be determined, and the safe flight range can be delineated for the approaching following aircraft.Figure17 shows the area where the wake behind the wing can cause the ARJ21-700 to roll when the front A330-200 is at the DH (i.e., at 30 m above the ground). The isosurface map of the rolling hazard zone is presented in Figure 17(a), where the yellow, blue, and red colors indicate the mild-, medium-, and severe-hazard zones, respectively. The roll hazard zone represents a symmetrical two-part area, and the hazard zone closer to the vortex core has a longer radiation range in the aircraft approach direction. Since the wake generated by the preceding aircraft sinks over time, the hazard zone does not completely coincide with the glide path. The top view of and horizontal boundary of the hazard zone are presented in Figure 17(b). As shown in Figure 17(b), when the distance between the front and following aircraft is at least 4.136 nmile, the following aircraft is not in the moderate-danger zone at any position in the wake vortex field, and it is in a safe state; the ICAO Recat separation is decreased by 17.28%. Figure 17(c) presents the side view of the hazard zone and the vertical boundary of the danger area. Figure 17(c) shows that the wake safety separation cannot be judged only based on the attenuation of the wake intensity because it also depends on the change in the position of the wake vortex. Setting the aircraft channel size in the vertical direction to 30 m, it is only necessary to ensure that the channel, and the dangerous wake area do not overlap. In Figure 17(c), it can be seen that the wake separation is inversely proportional to the wake vortex sinking rate. The faster the wake vortex sinks, the smaller the overlap between the channel and the wake danger area will be, and the smaller the wake separation will be. Figure 17(c) demonstrates that, according to this safety standard, in a standard atmospheric environment, the distance-based safety separation is 1.32 nm, which is 68% shorter than separation 1 and 73.56% shorter than the ICAO separation.Figure 17
Rollover hazard area. (a) Hazardous area isosurface; (b) the top view of the hazard zone; and (c) the side view of the hazard zone.
(a)(b)(c)The isosurface map of the overload hazard zone is presented in Figure18(a), where the yellow, blue, and red colors indicate the mild-, moderate-, and severe-turbulence areas, respectively. The areas on both sides of the overload risk are upwash overload, while the middle area denotes downwash overload. The upper- and downwash overload areas are divided by the center of the two vortex cores so that the radiation range of the downwash overload in the approach direction is longer. Figure 18(b) shows the top view of the overload hazard zone and the horizontal boundary of the turbulence area. Compared to the roll hazard zone, the overload hazard zone has a larger influence range in the horizontal direction but a smaller influence range in the approach direction; the moderate- and severe-bump areas are small. Further, Figure 18(b) demonstrates that even under mild turbulence as a boundary, the wake safety separation is decreased by 33.57% compared to the ICAO separation. Figure 18(c) presents the side view of the overload hazard zone and the vertical boundary of the turbulence area. The moderate- and severe-turbulence areas have a small vertical influence range. If mild turbulence is used as a judgment standard, the distance-based safety separation is 1.167 nmile, which is 64.88% shorter than separation 1 and 76.67% shorter than the ICAO separation.Figure 18
Overload hazard area. (a) Hazardous area isosurface; (b) the top view of the hazard zone; and (c) the side view of the hazard zone.
(a)(b)(c)The maximum capacity of an airport runway represents the reciprocal of the weighted average service time of all aircraft waiting in line for approach and landing, and it is calculated by(17)C=1∑Pi,jΔTij,where C is the maximum runway capacity; ΔTij represents the time separation for arriving at the runway threshold, where the front aircraft is denoted by i, and the following aircraft is denoted by j; and Pi,j represents the probability that the front aircraft is i and the following aircraft is j.Taking runway 02R/20L of Chengdu Tianfu International Airport as an example, and simulated the two dangerous situations of roll and overload, the minimum wake separation 1.32 nmile is determined. Taking the average ROT as 67.24 s [24], the maximum runway capacity under ROT constraints is 53.54 sorties/h in our simulation situation. The runway occupancy time could be shortened in the future, and when the average ROT is reduced to the minimum time separation, which is 1.32 nmile, the runway capacity under all this type combination would be 102.37 sorties/h. An increase of 75.37 sorties/h and 279.14% compared to the 27 sorties/h under the ICAO 5-nmile separation.


### 4.2. Evolution of Wake Vortex Hazard Zone at DH

However, when the wake vortex is close to the ground, due to the ground effect, the wake vortex will undergo elastic deformation during the evolution process. The front plane wake vortex will not continue to sink, and its height and vortex core spacing will change with time. The resulting danger zone also changes dynamically. In addition, crosswinds may keep the wake vortices above the runway for extended periods of time. This will make the following aircraft with a similar channel be likely to encounter a wake vortex. If the following aircraft encounters a strong wake in the near-ground phase, the plane can be easily stalled and hardly recover, causing severe aviation accidents.When the meteorological conditions of an airport are calm, the wake vortex develops symmetrically. The wake vortex at a DH of 30 m for class II precision approach is selected as a research object in this study. The wake hazard area is defined according to the wake movement trajectory of the front aircraft, the rolling moment coefficient, and the overload increment limit the following aircraft can bear. The evolution process of the hazard zone over time is analyzed to provide a reference for the following aircraft to avoid the wake hazard zone.The evolution process of the roll hazard zone is presented in Figure19, where it can be seen that the maximum rolling moment coefficient that a commercial aircraft can withstand is 0.05. To ensure the flight safety in the near-ground area, a certain safety margin needs to be reserved. In this study, the margin value is set to 0.5 times the safety factor, so the threshold value of the safe rolling torque factor is 0.025, and the dangerous area is defined by this value. The color in Figure 19 indicates the plane’s rolling direction; blue and red indicate the plane rolls clockwise and counterclockwise, respectively. The color depth indicates the magnitude of the rolling moment coefficient of the following aircraft, reflecting the severity of the following aircraft’s rolling. When the following aircraft is in the center of the vortex core of the front aircraft, the rolling moment is the largest; namely, this is the position with the highest rolling risk. Under the combined action of the ground effect and the secondary vortex, the strength of the wake vortex decreases with time. The wake vortex first sinks and then rebounds; the vortex spacing increases continuously; the roll hazard zone deforms and displaces. When t=40s, the wake vortex roll risk zone is far enough apart. The maximum width of the airport runway is 60 m. In this case, the roll hazard zone is completely separated from the runway and will not have an uncontrollable roll impact on the aircraft that is approaching and landing.Figure 19
Variation in the wake roll hazard zone with time at the 30-m DH.
(a)
t=10s(b)
t=20s(c)
t=30s(d)
t=40sFigure20 shows the evolution of the overload hazard zone for moderate turbulence as a limit. Because the aircraft is very close to the ground during the landing phase, to ensure the stability of the aircraft attitude, a safety margin of 0.25 needs to be added, so the safe overload increment threshold is 0.25, and this value is used to define the dangerous area. The color in Figure 20 indicates the direction of the plane’s force; blue and red colors indicate that the plane is under the action of the downwash and upwash forces, respectively. The shade of color indicates the magnitude of the aircraft’s overload. At the DH, because the aircraft is close to the ground and the aircraft engine is in an idling state, the downwash-induced speed of the wake will result in aircraft hard-landing and can cause a safety accident. In this case, the downwash load is more dangerous than the upwash load. As shown in Figure 20(a), when t=10s, the dangerous area is a whole, and the influence range of the downwash load is much larger than that of the upwash load. The downwash overload danger zone is a circular area with a radius of 25 m, and the maximum downwash load is greater than the maximum upwash load. At t=20s, the danger zone is extended to the spanwise range but shortened in the vertical range. In this case, the danger zone is deformed into an ellipse. At t=30s, with the increase in vortex spacing and circulation attenuation, the danger zone is divided into two symmetrical areas, and the area of the danger zone is obviously reduced. At t=40s, the two parts of the danger zone have left the runway, and the space between them is enough to allow the following aircraft to land on the runway safely.Figure 20
Variation in the wake overload hazard zone with time at a 30-m DH.
(a)
t=10s(b)
t=20s(c)
t=30s(d)
t=40sAs shown in Figure15, the color represents the absolute value. Figure 21 shows the evolution process of the roll hazard zone at different DHs. When H=5m or H=10m, the hazard zone moves to both sides faster than H=30m or H=60m. If there is a crosswind, it may cause a roll risk to the adjacent runway, and the danger area gradually moves upward due to the rebound of the wake vortex. The danger zone gradually moves upward due to the rebound of the wake vortex. When H=30m or H=60m, the roll hazard zone first drops and then rebounds and rises and also gradually moves to both sides. However, the speed of the hazard zone moving to both sides is significantly decreased. The area of the rolling danger zone gradually decreases with the attenuation of the circulation. In addition, the rolling moment coefficient in the center of the danger zone gradually decreases, but the danger zone does not show obvious deformation.Figure 21
Variation in the roll hazard zone with time at different DHs.
(a)
H=5m(b)
H=10m(c)
H=30m(d)
H=60mThe evolution process of the overload hazard zone is presented in Figure22, where it can be seen that the duration of the overload risk is less than that of the rollover risk. Generally, the overload risk is completely dissipated in the first 50 s–80 s. Further, the reduction rate of the overload hazard zone is significantly higher than that of the rollover hazard zone. Over time, the overload hazard zone is gradually elongated and divided into two areas and deformation occur. Namely, the middle area is compressed to the center in the longitudinal direction but elongated in the span, changing from a circle to an ellipse. The top and bottom are gradually recessed towards the center until they are completely divided into two danger zones. The separation and complete dissipation times of the danger zone at different DHs are also different. The higher the DH is, the longer the separation time of the danger zone will be.Figure 22
Variation in the roll hazard zone with time at different DHs.
(a)
H=5m(b)
H=10m(c)
H=30m(d)
H=60m


## 4.1. Risk Analysis of ARJ21 Encountering Wake in the Approach Section

To verify the safety of the current wake separation, a medium-size aircraft ARJ21 is selected to follow a heavy aircraft A330-200 to approach and land at the 5-nmile wake separation specified by the ICAO.When the ARJ21 and the front vortex field are at the same height (as shown in Figure11(a)), Figure 14 shows the relationship between the rolling moment coefficient of different parts and the change in overload increment as the ARJ21 enters the spanwise position of the front wake vortex field. Positive and negative values of the rolling torque coefficient represent counterclockwise and clockwise rolling types, respectively. Also, positive and negative overload values denote the upwash the downwash overloads, respectively. The absolute value indicates the degree of rolling and turbulence. The torque coefficient increases significantly near the vortex core, reaching the maximum at the center of the vortex core. Also, the rolling torque coefficient of the whole aircraft is 0.027, which is less than the threshold of 0.05. Different from the maximum value of the tail rolling moment coefficient of each part of the aircraft’s body, the maximum overload position of the wing is 10.6 m from both sides of the midpoint; the maximum overload position of the whole machine is 21 m from both sides of the midpoint. The maximum of RMC is 0.052, which is much smaller than the threshold value of 0.5. The rolling moment coefficient and the overload increment are both smaller than the corresponding critical values, which verifies the safety of the current wake separation.Figure 14
The coefficient of each part of the ARJ21 at different spanwise positions of the wake vortex field at the same height.
(a)(b)Because of the front vortex sinking, the ARJ21 and the front vortex field are typically at different heights, as shown in Figure11(b). The relationship between the rolling moment coefficient and the overload increment of the ARJ21 with the spanwise position of the ARJ21 entering the wake vortex field at different heights from the front vortex center is displayed in Figure 15. The spanwise position corresponding to the maximum rolling moment coefficient does not change with the height. In contrast, the spanwise position corresponding to the maximum overload gradually approaches the midpoint with the height. As the height increases, both the rolling moment coefficient and the overload increment decrease rapidly, and the rolling moment coefficient decreases significantly. This indicates that the judgment on whether the following aircraft is safe or not depends not only on the wake intensity of the front aircraft but also on the position of the following aircraft entering the wake vortex field. The existing wake separation still has a large space for reduction.Figure 15
The overall coefficient of the ARJ21 with the spanwise position of the wake vortex field at different heights.
(a)(b)To shorten the wake separation further, it is needed to delineate the wake danger zone of the front aircraft according to the movement trajectory and intensity attenuation changes of the front aircraft. Then, it is superimposed with the approach course of the following aircraft, and the length of the overlapping part is set as a minimum wake separation. The procedure of the final phase of the precision approach is illustrated in Figure16(a). When the front aircraft is along the precision approach glide path to the DH, the wake behind the wing will affect the aircraft approaching along the same glide path. Based on the wake encounter and aerodynamic response model, the wake vortex at the DH is used as a reference frame to analyze the force of the following aircraft. As shown in Figure 16(b), when the approach direction is the x-axis direction, the spanwise direction is the y-direction, and the vertical direction is the z-axis, the intersection of the section where the DH is located, and the extension line of the runway centerline is the origin. The longitudinal position of the initial wake vortex at the DH is set as the x-axis zero point, the center point of the two vortices is the y-axis zero point, and the airport surface height is the z-axis zero point. Further, the position where the nose of the following aircraft enters the wake field is denoted by (x, y, z). According to the wake movement trajectory of the front aircraft and its influence on the following aircraft, the areas with different wake risk levels are defined.Figure 16
The ARJ21 encounter wake coordinate system in the final approach segment. (a) Precision approach procedure and (b) the three-dimensional coordinate system of the wake encounter.
(a)(b)According to the maximum roll control power that can be controlled by the ailerons of an aircraft, the roll risk is categorized into three levels according to the roll moment coefficient. The roll risk of 0.025–0.05 is a mild danger risk. Namely, when aircraft rolls slightly, it is easier for a pilot to recover the aircraft. The roll risk of 0.05–0.08 is a moderate-danger risk, and in this case, the rolling slope angle increases, and the risk of stalling arises. Therefore, a pilot needs to put in larger efforts to recover and maintain the original flight path. When the roll risk is higher than 0.08, it is a severe danger risk. In this case, the aircraft’s body rolls and shakes obviously; the aircraft is facing a stall; and it is difficult for a pilot to recover.As illustrated in [24], the risk range of 0.15–0.5 corresponds to a mild turbulence area. In this case, the aircraft has a slight altitude change; brief yaw, pitch, and roll attitudes occur but have minimal effect. Since the aircraft overload is difficult to be larger than one during the approach phase, the risk range of 0.5–0.75 is defined as a moderate turbulence area in this study. In this area, the flight altitude and attitude of aircraft change significantly, but the aircraft is still in a controllable state. The risk of above 0.75 corresponds to a severe turbulence area, where flight attitude and altitude change drastically. In this area, the airspeed indicated by the instrument changes significantly, and the aircraft loses control for a short time. In severe cases, damage to the body structure can be caused, resulting in an aviation accident. According to the above hazard level categorization method, the wake hazard area on the trailing edge of the front wing can be visualized, boundaries of different hazard levels can be determined, and the safe flight range can be delineated for the approaching following aircraft.Figure17 shows the area where the wake behind the wing can cause the ARJ21-700 to roll when the front A330-200 is at the DH (i.e., at 30 m above the ground). The isosurface map of the rolling hazard zone is presented in Figure 17(a), where the yellow, blue, and red colors indicate the mild-, medium-, and severe-hazard zones, respectively. The roll hazard zone represents a symmetrical two-part area, and the hazard zone closer to the vortex core has a longer radiation range in the aircraft approach direction. Since the wake generated by the preceding aircraft sinks over time, the hazard zone does not completely coincide with the glide path. The top view of and horizontal boundary of the hazard zone are presented in Figure 17(b). As shown in Figure 17(b), when the distance between the front and following aircraft is at least 4.136 nmile, the following aircraft is not in the moderate-danger zone at any position in the wake vortex field, and it is in a safe state; the ICAO Recat separation is decreased by 17.28%. Figure 17(c) presents the side view of the hazard zone and the vertical boundary of the danger area. Figure 17(c) shows that the wake safety separation cannot be judged only based on the attenuation of the wake intensity because it also depends on the change in the position of the wake vortex. Setting the aircraft channel size in the vertical direction to 30 m, it is only necessary to ensure that the channel, and the dangerous wake area do not overlap. In Figure 17(c), it can be seen that the wake separation is inversely proportional to the wake vortex sinking rate. The faster the wake vortex sinks, the smaller the overlap between the channel and the wake danger area will be, and the smaller the wake separation will be. Figure 17(c) demonstrates that, according to this safety standard, in a standard atmospheric environment, the distance-based safety separation is 1.32 nm, which is 68% shorter than separation 1 and 73.56% shorter than the ICAO separation.Figure 17
Rollover hazard area. (a) Hazardous area isosurface; (b) the top view of the hazard zone; and (c) the side view of the hazard zone.
(a)(b)(c)The isosurface map of the overload hazard zone is presented in Figure18(a), where the yellow, blue, and red colors indicate the mild-, moderate-, and severe-turbulence areas, respectively. The areas on both sides of the overload risk are upwash overload, while the middle area denotes downwash overload. The upper- and downwash overload areas are divided by the center of the two vortex cores so that the radiation range of the downwash overload in the approach direction is longer. Figure 18(b) shows the top view of the overload hazard zone and the horizontal boundary of the turbulence area. Compared to the roll hazard zone, the overload hazard zone has a larger influence range in the horizontal direction but a smaller influence range in the approach direction; the moderate- and severe-bump areas are small. Further, Figure 18(b) demonstrates that even under mild turbulence as a boundary, the wake safety separation is decreased by 33.57% compared to the ICAO separation. Figure 18(c) presents the side view of the overload hazard zone and the vertical boundary of the turbulence area. The moderate- and severe-turbulence areas have a small vertical influence range. If mild turbulence is used as a judgment standard, the distance-based safety separation is 1.167 nmile, which is 64.88% shorter than separation 1 and 76.67% shorter than the ICAO separation.Figure 18
Overload hazard area. (a) Hazardous area isosurface; (b) the top view of the hazard zone; and (c) the side view of the hazard zone.
(a)(b)(c)The maximum capacity of an airport runway represents the reciprocal of the weighted average service time of all aircraft waiting in line for approach and landing, and it is calculated by(17)C=1∑Pi,jΔTij,where C is the maximum runway capacity; ΔTij represents the time separation for arriving at the runway threshold, where the front aircraft is denoted by i, and the following aircraft is denoted by j; and Pi,j represents the probability that the front aircraft is i and the following aircraft is j.Taking runway 02R/20L of Chengdu Tianfu International Airport as an example, and simulated the two dangerous situations of roll and overload, the minimum wake separation 1.32 nmile is determined. Taking the average ROT as 67.24 s [24], the maximum runway capacity under ROT constraints is 53.54 sorties/h in our simulation situation. The runway occupancy time could be shortened in the future, and when the average ROT is reduced to the minimum time separation, which is 1.32 nmile, the runway capacity under all this type combination would be 102.37 sorties/h. An increase of 75.37 sorties/h and 279.14% compared to the 27 sorties/h under the ICAO 5-nmile separation.


## 4.2. Evolution of Wake Vortex Hazard Zone at DH

However, when the wake vortex is close to the ground, due to the ground effect, the wake vortex will undergo elastic deformation during the evolution process. The front plane wake vortex will not continue to sink, and its height and vortex core spacing will change with time. The resulting danger zone also changes dynamically. In addition, crosswinds may keep the wake vortices above the runway for extended periods of time. This will make the following aircraft with a similar channel be likely to encounter a wake vortex. If the following aircraft encounters a strong wake in the near-ground phase, the plane can be easily stalled and hardly recover, causing severe aviation accidents.When the meteorological conditions of an airport are calm, the wake vortex develops symmetrically. The wake vortex at a DH of 30 m for class II precision approach is selected as a research object in this study. The wake hazard area is defined according to the wake movement trajectory of the front aircraft, the rolling moment coefficient, and the overload increment limit the following aircraft can bear. The evolution process of the hazard zone over time is analyzed to provide a reference for the following aircraft to avoid the wake hazard zone.The evolution process of the roll hazard zone is presented in Figure19, where it can be seen that the maximum rolling moment coefficient that a commercial aircraft can withstand is 0.05. To ensure the flight safety in the near-ground area, a certain safety margin needs to be reserved. In this study, the margin value is set to 0.5 times the safety factor, so the threshold value of the safe rolling torque factor is 0.025, and the dangerous area is defined by this value. The color in Figure 19 indicates the plane’s rolling direction; blue and red indicate the plane rolls clockwise and counterclockwise, respectively. The color depth indicates the magnitude of the rolling moment coefficient of the following aircraft, reflecting the severity of the following aircraft’s rolling. When the following aircraft is in the center of the vortex core of the front aircraft, the rolling moment is the largest; namely, this is the position with the highest rolling risk. Under the combined action of the ground effect and the secondary vortex, the strength of the wake vortex decreases with time. The wake vortex first sinks and then rebounds; the vortex spacing increases continuously; the roll hazard zone deforms and displaces. When t=40s, the wake vortex roll risk zone is far enough apart. The maximum width of the airport runway is 60 m. In this case, the roll hazard zone is completely separated from the runway and will not have an uncontrollable roll impact on the aircraft that is approaching and landing.Figure 19
Variation in the wake roll hazard zone with time at the 30-m DH.
(a)
t=10s(b)
t=20s(c)
t=30s(d)
t=40sFigure20 shows the evolution of the overload hazard zone for moderate turbulence as a limit. Because the aircraft is very close to the ground during the landing phase, to ensure the stability of the aircraft attitude, a safety margin of 0.25 needs to be added, so the safe overload increment threshold is 0.25, and this value is used to define the dangerous area. The color in Figure 20 indicates the direction of the plane’s force; blue and red colors indicate that the plane is under the action of the downwash and upwash forces, respectively. The shade of color indicates the magnitude of the aircraft’s overload. At the DH, because the aircraft is close to the ground and the aircraft engine is in an idling state, the downwash-induced speed of the wake will result in aircraft hard-landing and can cause a safety accident. In this case, the downwash load is more dangerous than the upwash load. As shown in Figure 20(a), when t=10s, the dangerous area is a whole, and the influence range of the downwash load is much larger than that of the upwash load. The downwash overload danger zone is a circular area with a radius of 25 m, and the maximum downwash load is greater than the maximum upwash load. At t=20s, the danger zone is extended to the spanwise range but shortened in the vertical range. In this case, the danger zone is deformed into an ellipse. At t=30s, with the increase in vortex spacing and circulation attenuation, the danger zone is divided into two symmetrical areas, and the area of the danger zone is obviously reduced. At t=40s, the two parts of the danger zone have left the runway, and the space between them is enough to allow the following aircraft to land on the runway safely.Figure 20
Variation in the wake overload hazard zone with time at a 30-m DH.
(a)
t=10s(b)
t=20s(c)
t=30s(d)
t=40sAs shown in Figure15, the color represents the absolute value. Figure 21 shows the evolution process of the roll hazard zone at different DHs. When H=5m or H=10m, the hazard zone moves to both sides faster than H=30m or H=60m. If there is a crosswind, it may cause a roll risk to the adjacent runway, and the danger area gradually moves upward due to the rebound of the wake vortex. The danger zone gradually moves upward due to the rebound of the wake vortex. When H=30m or H=60m, the roll hazard zone first drops and then rebounds and rises and also gradually moves to both sides. However, the speed of the hazard zone moving to both sides is significantly decreased. The area of the rolling danger zone gradually decreases with the attenuation of the circulation. In addition, the rolling moment coefficient in the center of the danger zone gradually decreases, but the danger zone does not show obvious deformation.Figure 21
Variation in the roll hazard zone with time at different DHs.
(a)
H=5m(b)
H=10m(c)
H=30m(d)
H=60mThe evolution process of the overload hazard zone is presented in Figure22, where it can be seen that the duration of the overload risk is less than that of the rollover risk. Generally, the overload risk is completely dissipated in the first 50 s–80 s. Further, the reduction rate of the overload hazard zone is significantly higher than that of the rollover hazard zone. Over time, the overload hazard zone is gradually elongated and divided into two areas and deformation occur. Namely, the middle area is compressed to the center in the longitudinal direction but elongated in the span, changing from a circle to an ellipse. The top and bottom are gradually recessed towards the center until they are completely divided into two danger zones. The separation and complete dissipation times of the danger zone at different DHs are also different. The higher the DH is, the longer the separation time of the danger zone will be.Figure 22
Variation in the roll hazard zone with time at different DHs.
(a)
H=5m(b)
H=10m(c)
H=30m(d)
H=60m


## 5. Conclusion

In this paper, a simulation of medium aircraft following a heavy aircraft during the approach and landing stages is conducted. A330-200 aircraft is selected as the front aircraft, and the RANS method is employed for simulating the wake vortex evolution process. By constructing the wing structured grid and designing the wake vortex tangential velocity model, the law of the near- and far-field evolution of the wake vortex at different DHs is analyzed. A medium-size aircraft ARJ21 is selected as the following aircraft to construct a wake encounter model. The strip method is used to examine the aerodynamic response characteristics of the wake encounter. The safety is verified based on the rolling moment factor and overload increment. According to the changes in the wake vortex intensity and front aircraft position, the front aircraft wake hazard zone is defined. Also, the minimum wake separation is determined according to the overlapping range of the hazard zone and the approach channel of the following aircraft. Finally, the dynamic change in the wake hazard zone at different DHs with time is studied.(1)
The ground effect can accelerate both the deformation of the wake vortex and the decay of the vortex structure. It can increase the vortex spacing, make the wake vortex rebound, and accelerate the wake vortex to enter the rapid decay stage. At the same time, the boundary layer will be induced to separate from the ground, and the separation layer will continue to increase to form a secondary vortex. The promotion effect of the secondary vortex on the wake vortex attenuation occurs together with the suppression effect of the increase in the vortex spacing on the wake vortex attenuation, but latter effect is more obvious. WhenH=30m, the decay of the wake vortex circulation is the slowest(2)
The wake safety separation cannot be judged only based on the wake intensity attenuation but also considering the change in the wake vortex position. Judging the wake safety separation only based on the wake intensity, the minimum safe separation for the ARJ21 to follow the A330-200 approach is 4.163 nmiles, which is 17.28% shorter than that of the ICAO standard of 5 nmile. When combining the sinking of the front vortex and the approach channel of the following aircraft to determine the hazard zone and the overlapping range of the channel, the minimum distance between the two aircraft is 1.32 nmile under ideal conditions, which is 73.56% shorter than that of the ICAO standard. After the average ROT shortening, there is plenty of room for runway capacity improvement. When the average ROT is reduced to the minimum time separation based on a separation of 1.32 nmile, the capacity of runway 02R/20 L in Chengdu Tianfu International Airport could theoretically reach 102.37 sorties/h under all this aircraft pair combination. Compared to 27 sorties/h of the ICAO separation standard, the runway capacity is improved up to 279.14%, increasing the airport operation efficiency effectively(3)
The wake at the DH is affected by the ground effect, and the shape and position of a hazard zone change rapidly. When a hazard zone is deformed, it will rebound and spread to both sides of the runway. The overload hazard zone is elongated over time in the spanwise direction, changing from an integral area to two independent areas. The influence range of the downwash load is much larger than that of the upwash load. The shape of the roll risk zone does not change with time significantly, but the roll risk lasts longer than the overload risk. The lower the DHs are, the faster the hazard zone leave off the runway, which provides a reference for the following aircraft to avoid the hazard zone and land on a runway safelyIn the future, combined with QAR data and wind field data detected by radar, and considering the change of aircraft speed, the wake hazard zone under the atmospheric wind field will be predicted, and the dynamic wake separation under different combinations will be provided. And a rapid wake dynamic separation prediction system based on numerical simulation will be developed.



---

*Source: 9126755-2022-09-05.xml*