1 Introduction

The wake vortex of heavy aircrafts poses danger to the following light planes. As the air traffic increases, this threat limits the efficiency and security of airports. It is meaningful to study the evolution of aircraft wake vortex for reducing the wake separations to increase the airport capacity. As the development of computer technology progresses, numerical simulation has become an effective method in wake vortex research. In the early works, the two-dimensional (2D) Navier-Stokes equations were solved. Robins and Delisi.^[1] used the 2D large eddy simulation (LES) to study the effects of vertical shear and stable stratification on a pair of vortices. The terminal area simulation system (TASS) model was used to simulate wake vortex evolution numerically with the 2D LES, and the model results were compared with field measurement in Idaho and Memphis^[2-3]. Later on, Proctor studied the sensitivity of vortex behavior to the atmospheric parameters, such as crosswind shear, stratification, and the 2D turbulence^[4].

Wake vortices will decay by not only diffusion, but also the linear instabilities such as Crow's long wave instability^[5]. The Crow instability is not the unique mechanism in distortion of wake vortex. The three-dimensional (3D) secondary vortex structures, induced by strong turbulence and stable stratification, also play an important role in vortex bursting, indicating that 3D numerical simulations are necessary. In 1999, Shen et al.^[6] simulated the behavior of wake vortex measured in 1995 Memphis field experiments using a 3D LES. The influence of initial atmospheric turbulence and stable stratification was studied by Proctor and Switzer^[7] with the 3D LES. German scientists have performed a series of 3D LES since 1999^[8-14], and systematically revealed the behavior of vortices in complex atmospheric environment. It is recognized that the wake vortex undergoes a two-phase decay: the diffusion phase and the rapid decay phase. In a nearly neutral atmosphere with weak turbulence conditions, the rapid decay phase is caused by Crow's long wave instability, while in strong turbulence or stable stratification, the secondary vortices formed by turbulence disturbance and baroclinity are the main cause of the rapid decay.

As is significant for numerical simulation in any wake vortex study, it is meaningful to improve the computation efficiency while maintaining a high-degree of accuracy. The large velocity gradient inside and near the vortex core region requires a high resolution grid in the plane perpendicular to the vortex axis in numerical simulation. As the vortex core is transported via mutual-induction, a fine uniform grid is needed for the whole computation domain if using a fixed mesh. In Ref. ^[6], a simulation domain of 104 m$\times$202 m$\times$200 m (axial, lateral, and vertical directions), or $2.6 b_0\times5.1 b_0\times5.1b_0$, was used with the grid points $52\times202\times200$^[6], where $b_0$ is the initial vortex separation. In order to recognize Crow's longwave instability, the axial domain length should be at least 8 times of the vortex separation $b_0$. Holzäpfel et al.^[11] made the axial domain length 408 m (8.7$b_0$) with the grid spacing 6.375 m. In later works, the uniform resolution of 1 meter in all three directions was used to increase the axial resolution for more accurate calculation of vorticity in non-axial direction, which made the total grid points number as large as 26 million^[13-14]. However, the need for a fine grid is not as important in regions far from the vortex core. There is much room for reducing the total grid number. In the current work, a self-adaptive-grid LES method is applied to the wake vortex numerical simulation. Self-adaptive-grid method can obtain fine grid in high velocity gradient regions such as vortex cores, and coarse grid points in low velocity gradient areas. As the vortex moves, the fine-grid region can follow closely around the vortex core. This method decreases the grid point number to about 6 million with a resolution of about 0.5 m in the core area, which considerably improves the efficiency of wake vortex simulation. The Idaho Falls (IDF) measurement of run 9 case^[3] is simulated numerically and compared with that of the field experiment data. In addition, six cases with different ambient atmospheric turbulence and Brunt Väisälä (BV) frequencies are investegated. In Section 2, the numerical method and the adaptive algorithm are described in detail. In Section 3, the simulation parameters are introduced. Simulation results are discussed in Section 4. Finally, the concluding remarks are displayed in Section 5.

2 Numerical method

The code of Tsinghua Turbulence Laboratory large eddy simulation (TTLES) solves the incompressible Navier-Stokes equations on a non-staggered grid. The Boussinesq approximation is used to simulate the buoyancy effect. The TTLES code is originally used in atmospheric environment research and performed well in micro-scale urban environments numerical simulation^[15-16]. By filtering the Navier-Stokes equations, we can obtain the LES equations as follows:

in which

where $Pr_{\rm t}$, the turbulent Prandtl number, is set to be 0.72 in this paper. $\theta_0$ is the reference temperature which is 300 K in computation. A fourth-order finite volume method is used in space, and the fourth-order Runge-Kutta integration is implemented in time advancement. The numerical formulation is described in detail by Shi et al.^[15] and Xu et al.^[17].

The code of Adaptive Tsinghua Turbulence Laboratory large eddy simulation (ATTLES) is developed from the TTLES. The major change is that the ATTLES is implemented on a moving cartesian grid frame self-adapted to the flow field. The transformation to the time-depend computation domain can be written as follows:

where $(\xi, \eta, \zeta, \tau)$ are the coordinates in the computation domain. The LES equations in the moving frame with an eddy viscosity model can be written as follows:

in which

The moving velocity of grid points $u_{\rm g}$, $v_{\rm g}$, and $w_{\rm g}$ is shown in Fig. 1. Similar to the TTLES, a finite volume method, a fourth-order compact scheme, and a fourth-order Runge-Kutta method are applied to Eq. (4) . A Lagrangian dynamic Smagorinsky model^[18] is applied, which can reduce the excessive eddy viscosity in the vortex core^[14]. The partial derivatives $x_\xi$, $y_\eta$, and $z_\varsigma$ can be obtained from the grid spacing in the corresponding directions. The grid plane speeds $u_{\mathrm g}$, $v_{\mathrm g}$, and $w_{\mathrm g}$ are determined by a self-adaptive algorithm. The adaptive algorithm is similar to the spring analogy method by Gnoffo^[19]. As shown in Fig. 1, the grid planes are supposed to be connected with springs. Assuming the quasi-equilibrium model that the elastic forces on the $i$th grid point reach balance in the $x$-direction, one may have the balance equation

where $K_{i+\frac{1}{2}}$ is the stiffness of the spring between the $i$th and ($i+1$)th grid point and is modeled as

In the above equations, $a_i$ refer to the maximum velocity gradient in the $i$th grid plane, and $a_{{\rm max}}$ and $a_{{\rm min}}$ are the maximum and minimum of $a_i$ ($i=1, 2, \cdots, N_x$), respectively. The constant $A$ is the ratio of maximum and minimum grid spacing in the $x$-direction. $F(f)$ in Eq. (7) is a function valued within [0, 1]. Nakahashi and Deiwert^[20] suggested $F(f)=f^B$, where $B$ is a positive constant. In this paper, $F(f)$ takes the form of $F(f) ={\min}\{1.0, 2{\sin}(f)\}$. This setting will enlarge the fine grid region near the vortex core to improve the numerical accuracy. With these settings, springs in high gradient regions have greater stiffness. When the system achieves balance, fine grid is obtained in these areas. The details can be found in Nakahashi and Deiwert^[20].

In this paper, the grid is updated each time step by

where $u_{\rm g}$, $v_{\rm g}$, and $w_{\rm g}$ are updated every $N_u$ steps by the adaptive algorithm. Take the $x$-direction for example, the adaptation algorithm is implemented as follows:

(i) Compute and store $K_{i+1/2}$ by Eq. (7) .

(ii) Solve a tri-diagonal linear system of equations of size $N_x$ for the equilibrium position $x'_i$ of the grid planes

(iii) Update $u_{\rm g}$,

where $x_i$ is the current position of the $i$th grid.

(iv) Update $v_{\rm g}$ and $w_{\rm g}$ in the similar way.

In this way, the grid points need a time of $N_u\Delta t$ to move to the equilibrium position $x'_i$, which means that $N_u\Delta t$ should be as short as possible to ensure the instantaneity. However, in consideration of numerical stability, the grid speed $u_{\rm g}$ needs to be limited, so too short a $N_u\Delta t$ is not feasible. In this paper, $N_u\Delta t$ is chosen to be 0.1 second, during which the vortex core undergoes a tiny displacement.

The computation time needed by adaptive grid is mainly consumed on updating the grid speed. As this updating is applied every $N_u$ time steps, the additional computation time is negligible when $N_u$ is set to 10 or larger. Therefore, this method is efficient when the total grid number is remarkably reduced.

3 Initial conditions

The velocity field is initialized with a pair of counter-rotating vortices superimposed upon a background field of turbulence. The Burnham Hallock model^[21] is used.

In Proctor's IDF B757 run 9 LES^[3], the core radius $r_{\rm c}$ is set to be 5% wingspan of the generating aircraft, corresponding to $r_{\rm c}=0.06b_0$ in this paper. The initial turbulence will be generated by the isotropic turbulence model with given turbulent kinetic energy dissipation rate $\varepsilon$ by Rogallo's method^[22]. The turbulent velocity is expanded in the spectral space, and the random isotropic turbulent velocity is taken the following form and determined by the prescribed kinetic energy spectrum:

In Eq. (12) , $\theta_1$, $\theta_2$, and $\phi$ are the random numbers in $[0, 2\pi]$ with the uniform probability density. A modified von Karman spectrum is used for the kinetic energy^[23] as follows:

The spectrum reaches its maximum value when $k=k_{\rm p}$. In the atmospheric boundary layer, the energy containing scale $2\pi/k_{\rm p}$ has the magnitude of altitude from ground^[24]. In this paper, $k_{\rm p}$ is valued $2\pi/90 {\rm m}^{-1}$, which indicates that the energy containing scale equals 90 m accepted in atmosphere science. $K_0$ is determined from the relation $\varepsilon=2\nu \int{{k^2}E( k )}{\rm d}k$. The turbulence field is at first produced on a uniform mesh of $256^3$ and then interpolated onto the computational grid. The energy density spectrum of initial turbulence is shown in Fig. 2.

The IDF B-757 run 9 case was simulated numerically by Proctor^[3] with the 2D LES, and the results are compared with the field experiment data from the IDF field study^[25] sponsored by the Federal Aviation Administration (FAA) to evaluate B-757 and B-767 aircraft wakes. The temporal evolution of vortex circulation and position are recorded, and the real-time profiles of crosswind and potential temperature are also measured. Because the real-time atmospheric condition is provided detailed in this case, it is selected to examine the model performance in complex atmospheric conditions with strong ground effect. In this paper, this case is simulated with a 3D LES. The crosswind (wind speed in lateral direction) and potential temperature profiles are initialized with the measured data^[4](see Fig. 3) , and a turbulence field with a dissipation rate $\varepsilon=7.5\times10^{-6} {\rm m}^2/{\rm s}^3$ is added to the initial field. The vortex parameters are $\Gamma_0=365 {\rm m} ^2/{\rm s}$ and $b_0=30$ m representing a landing Boeing B-757-200. The initial altitude of wake vortex is 70 m above the ground. The simulation domain is $L_x\times L_y\times L_z=20.0 b_0\times8.0 b_0\times4.83 b_0$ with the total grid points $N_x\times N_y\times N_z=244\times240\times74$, where $x$, $y$, and $z$ refer to the lateral, axial, and vertical directions, respectively. The periodic boundary conditions are applied to the axial and lateral directions. The top boundary is treated with symmetry conditions, which means that the normal gradient and normal velocity components are set to be zero. The ground surface is treated with no-slip condition. A fixed time step of $\Delta t=1/150 {\rm s}$ is used.

In addition, six cases with different ambient atmospheric turbulence and BV frequency are simulated numerically. The vortex parameters are $\Gamma_0=446 {\rm m}^2/{\rm s}$ and $b_0=47.4$ m with $w_0=1.50$ m/s and $t_0=31.7$ s. The atmos-pheric parameters are shown in Table 1. In Case N01, `N' means neutral atmosphere, and the number 01 stands for the normalized eddy dissipation $\varepsilon^*=0.01$ (normalized by $\varepsilon^* = {({\varepsilon {b_0}} )^{1/3}}/{w_0}$), while `S06' stands for the stable atmosphere with a uniform BV frequency $N^*=t_0N=0.6$. In the cases of neutral atmosphere, the temperature equation will not be solved. In the stable stratification cases, the temperature field is initialized with a constant vertical gradient ${\rm d}\theta/{\rm d}z$, which is determined by a prescribed BV frequency $N=((g/\theta_0) ({\rm d}\theta/{\rm d}z))^{1/2}$. All these cases are simulated in a domain size of $L_x\times L_y\times L_z=6.3 b_0\times 8.0 b_0\times 6.3 b_0$. All three directions are treated as periodic boundary except for cases with non-zero ${\rm d}\theta/{\rm d}z$, in which the symmetry condition is used for bottom and top boundaries. The self-adaptive-grid is applied to the lateral and vertical directions with the initial minimum grid spacing $\Delta x=\Delta z=0.5$ m in the core area. A uniform grid of $\Delta y=1.0$ m is used in the axial direction, which is the same as the work of Hennemann and Holzäpfel^[13] and Misaka et al.^[14]. The most important structure in the axial direction is Crow's instability, the wave length of which is in the order of the vortex separation $b_0$, much larger than the current axial resolution 1 m. Therefore, the current axial resolution is enough. The total grid number is $N_x\times N_y\times N_z=144\times 380\times 120$. A fixed time step of $\Delta t=1/150$ s is used.

In this paper, the vortex strength is estimated with the circulation, which is calculated by Stokes's theorem

The vortex center is determined by the extreme point of $\lambda_2$, which is the second eigenvalue of ${S}^2+{\Omega}^2$. Jeong and Hussain^[26] proved that $\lambda_2< 0$ indicates low pressure, which is thought to be inside a vortex core. Hence, the extreme of $\lambda_2$ is chosen to show the vortex center because it performs well in the vortex identification under strong external strain, especially for the secondary vortices at the edge of core.

The averaged circulation $\Gamma _{r_{1}-r_{2}}$ is calculated by

In this paper, the circulation $\Gamma_{r_{1}-r_{2}}$ is further averaged in the axial direction. In the IDF B-757 run 9 case, $\Gamma _{(0-10) {\rm m}}$ is used to allow a comparison to Proctor's results, while in other cases, $\Gamma _{(5-15) {\rm m}}$ (normalized by its initial value) is used since it is a widely used character to reveal the evolution of main vortex.

4 Results and discussion 4.1 Performance of self-adaptive-grid

Figure 4(a) shows the grid performance in the IDF B-757 run 9 case. The vortex cores are descending due to mutual-induction and undergo a lateral offset under the influence of crosswind. The wake vortex separation increases as a result of ground effect^{[12, 26]}. It can be seen that the grid is well adapted to the velocity field and fine grid region is obtained near the vortex core and the ground surface. Figures 4(b) and 4(c) show the evolution of computation grid in Cases N05 and S10. For Case N05, as the crow instability develops, the vortices deviate from the base line and links to form a vortex ring. For Case S10, the secondary structure appears and the vortices burst. The fine grid regions are expanded to ensure enough resolution in the near core areas. The ratio of the maximum and minimum grid spacings is chosen to be 6 at the beginning of simulation. As the expansion of fine grid region, the parameter needs increasing to ensure enough resolution in the core region. The minimum grid spacing increases from the original 0.5 m to about 0.8 m.

Figure 5 shows the adaptive mesh near the left vortex core in Case N05. As the function $F(f) = {\min}(1.0, 2{\sin}(f))$ is used (see Section 2) , the mesh in the area of radius $2r_{\rm c}$ around the core is practically uniform. This local-uniform mesh can reduce the error near the core area, including the numerical error and the error caused by the change in the LES filter scale when turbulence flows through mesh of different grid sizes.

4.2 Numerical simulation of real case: IDF B-757 run 9 case

The results of numerical simulation of IDF B-757 run 9 case are introduced in this section. Figure 6(a) shows the evolution of 0--10 meter averaged circulation $\Gamma_{(0-10) {\rm m}}$. The results are compared with the 2D LES by Proctor^[3] and the measured data by the FAA. It is shown that both the upstream and downstream vortices undergo a obvious rapid decay at time $t$ = 60 s and 25 s, which is in agreement with the measurement, and was not observed in Proctor's 2D LES results. It indicates that 3D structures dominate the rapid decay. In Fig. 6(b), the downstream vortex descends to an altitude of about 40 m at $t$ = 25 s, where a strong crosswind shear is located (see Fig. 3) . This time corresponds to the rapid decay time. The downstream vortex responds under the influence of shear zone leading to the ascending of vortex core (see Fig. 6(b)), which is also observed by Proctor et al.^[4]. For the upstream vortex, it descends through the shear zone without responding or rapid decay (see Fig. 6(b)) because the shear has the same direction of the vortex. Under the effect of ground, the circulation decay rate increases at \linebreak $t$ = 50 s^[12]. Figure 6(c) shows the temporal evolution of vortex lateral position. The results of present work are in agreement with Proctor' LES and the experiment data. It can be seen that the separation of vortices increases under the induction of image vortices below the ground^[12].

The ground effect vorticity is generated near the ground by the no-slip condition at the ground. In the present results, the ground effect vorticity can be observed obviously. Figure 7 shows the axial vorticity field near the upstream primary vortex in the $xz$-plane. In Fig. 7(a), a positive (anticlockwise) vorticity zone appears in the left of primary vortex. Under the induction of it, the primary vortex responds from the ground. The response of downstream vortex above the shear zone has the analogous mechanism. In Fig. 7(b), a zone of negative (clockwise) vorticity appears at $t$ = 25 s in the right of primary vortex, which is produced by the crosswind shear, rather than the ground effect.

4.3 Vortex structure and two-phase decay

The six designed cases with different turbulence intensities and stable stratifications are used to show the relationship between vortex structures and its decay. Figure 8(a) shows the decay of $\Gamma_{(5-15) {\rm m}}$ for all cases. The two-phase decay of circulation, i.e., the diffusion phase and the rapid decay phase, can be seen obviously in all cases except for Case N23 ($\varepsilon^{*}=0.23$, $N^{*}=0.0$), in which the strong turbulence reduces the circulation of primary vortex from the very beginning^[14]. In the diffusion phase, the decay rate mainly depends on $\varepsilon^{*}$. The cases with the same $\varepsilon^{*}$ (N05, S06, S10) have the similar decay rate and the circulation decays faster under stronger turbulence conditions, which is in agreement with other existing simulations^[13-14]. It is worth mentioning that the decay rate in the later period of diffusion phase is greater than that at the beginning in Cases N05 and N10. The reason for that is the development of turbulent eddies (see Figs. 10 and 11) . The comparison with the results of Hennemann and Holzäpfel^[13] is shown in Fig. 8(b). The decay rate of diffusion phase is smaller than that in Ref. ^[13]. The reason is that a Smagorinsky closure model is used in Ref. ^[13], which leads to a stronger dissipation. In addition, the grid spacing can also make this difference. In Ref. ^[13], the resolution in the core is 1 m, while it is 0.5 m in the present work. In Ref. ^[14], it is proved that a larger grid spacing leads to a greater decay in the diffusion phase.

Figure 8(c) shows the vortex descending normalized by $b_0$. The vortices descend with the speed $w_0$ at the beginning. In the rapid decay phase, the strength of primary vortex decreases and leads to a smaller descending rate. In the case S10 ($\varepsilon^{*}=0.05$, $N^{*}=1.0$), the vortex rebounds under the contribution of buoyancy. These results are in agreement with the existing studies^{[4, 27]}.

The onset time of rapid decay $t_2^*$ is an important reference value for the circulation decay. Figure 9 shows $t_2^*$ of the present cases, the LES results^{[7, 13]}, and the empirical formula by Sarpkaya^[28] and Holzäpfel et al.^[11] (the P2P model). For the cases in which the vortices link and form a ring, such as Cases N01, N05, and N10, $t_2^*$ is the link time. $t_2^*$ of this case is compared with the link time model^[27] in Fig. 9(a). Cases N01 and N05 are in agreement with Ref. ^[27], while in Case N10, the vortices link later than the model forecasts. The reason is that the strength of primary vortices decays under the influence of background turbulence, which slows down the linking process. For Case N23, the vortices have no chance to link (see Fig. 11(b)). Therefore, it is meaningless to compare it with Sarpkaya's model. $t_2^*$ in Case N23 is similar to Proctor's LES7, but earlier than the result in Ref. ^[13]. The cases with stratification are shown in Fig. 9(b) with the LES results and P2P model, and the present estimation of the onset time of rapid decay is generally in agreement with the previous studies. According to Hennemann and Holzäpfel^[13], the turbulent integral length scales also have influence on $t_2^*$. For this reason, $t_2^*$ from different researchers can be different.

According to Crow's theory^[5], the vortices undergo a nearly sinusoidal instability. As the instability develops, the two vortices links and at last form a train of vortex rings. Figure 10(a) shows this process of Case N05 by the iso-surface of $\lambda_2=-0.1$. At $t^{*}=5.7$, the vortices begin to link, which is corresponding to the rapid decay time $t_2^*$. The helical vorticity structure^{[14, 29]} could be seen at the link point and propagated along the vortex tube. At $t^{*}=7.0$, a bone-like vortex ring formed, which is a similar to the results of the LES^[13-14]. After that, the ring is deformed and breaks at last. The decay process of Case S10 is shown in Fig. 10(b). Different from N05, the vortices decay rapidly under the strong stable stratification and expand into pieces before the forming of a ring. A large amount of secondary vortices appears at time $t^{*}=1.5$, which is the onset of rapid decay. It can be seen the baroclinic vorticity is the main source of secondary vortices.

Figure 11 shows the evolution of vortices in Cases N10 and N23 with a stronger background turbulence $\varepsilon^*=0.10$ and $\varepsilon^*=0.23$. In Case N10 (see Fig. 11(a)), the vortices link at time $t^*=5.5$, corresponding to the rapid decay time (see Fig. 9) . Similar to Case N05, a vortex ring forms, but it turns into pieces soon as the development of secondary vortices, which comes from the stretched turbulent eddies^[14]. In this condition, the vortex linking and the small scale instabilities have equally importance in the circulation decay. In Case N23 (see Fig. 11(b)), the primary vortices decay rapidly in a much stronger ambient turbulence, and the vortex crushes at time $t^*=6.0$. For this reason, the vortices have no chance to link. This means that the main decay mechanism is the small scale instabilities rather than linking.

5 Conclusions

In this paper, a self-adaptive-grid LES is performed for the numerical simulation of aircraft tail vortex. This method greatly improves the computational efficiency by considerably reducing the total grid number. The core resolution is set to be $\Delta x\times\Delta y\times \Delta z =0.5 {\rm m}\times1.0 {\rm m}\times0.5 {\rm m}$. Only about six million grid points are needed in the whole domain, while this number will be as large as 34 million when using a uniform mesh with the same core resolution. In view of this, the improvement made by self-adaptive-grid in computational efficiency is quite considerable.

A real case numerical simulation (IDF B-757 run 9 case) is implemented. The evolution of circulation and the position of vortices fit the experimental data very well for both the upstream and downstream vortices. The influence of crosswind shear and ground effect is accurately simulated by the method. The vorticity zone produced by ground and crosswind shear zones can be obviously observed in the simulation results. The method is proved to be reliable in the LES of aircraft wake vortex in complex atmosphere conditions. In addition, six cases with different atmospheric turbulences and BV frequencies are implemented with the method. The two phase decay of circulation can be seen in the results. The global vortex characteristics, such as the rapid decay time and the descending rate, are in well agreement with the previous studies. Besides, the important features of vortex evolution, such as the secondary vortices, the form of vortex ring, and the helical structures, can be distinguished in the results.

In summary, the self-adaptive-grid method is proved to be practical in the numerical research of tail vortex. With this method, the resolution in the core can be further increased under the limited computation ability. It is noteworthy that the sensibility to the axial resolution is not analyzed in this paper. Besides, as the self-adaptive method is not applied to the axial direction, the simulation of non-axial vorticity may be less accurate, which needs further works.