Numerical study on aerodynamic performance during the process of entering the formation flight

. Two or more planes can fly in close formation similar to migratory birds, making use of the wing-tip vortex of the leader plane to increase lift and reduce drag, thereby effectively improving the flying range. By conducting wind tunnel tests and numerical simulations, the aerodynamic performance of formation flight at different relative positions can be obtained, the optimal formation position can thereby be solved. However, significant nonlinear and unsteady aerodynamic characteristics, induced by the interference between the follower plane and the wing-tip vortices of the leader plane, will affect the flight safety of the whole formation. At present, there is no effective prediction methods. Numerical simulations adopting adaptive grid refinement and dynamic overset grid were conducted for the dynamic entering process of the formation flight containing two Ty-154 planes. The aerodynamic characteristics and vortex interferences were analysed considering the effects of approaching direction and speed. The results indicate that there is deterioration of stability at the position where maximum lift gain is reached; Compared to the entering speed, the impact of directions on dynamic aerodynamic characteristics is more significant.


Introduction
In 2001, Nature published the research on the long-distance migration of Tang geese using formation flight.The study found that by means of formation flight 11~14% of energy can be saved, and Tang geese flying in formation reach farther than those flying alone [1].Inspired by it, studies have been made on the aerodynamic characteristics of aircraft formation flight and the concept of Surfing Aircraft Voices for Energy (SAVE) has been proposed.It specifically refers to the close formation flight of two or more planes similar to migratory birds, during which the rear plane "rides" on the vortex of the front plane, achieving significant lift-drag ration promotion.
In a formation flight, the airflow around the front plane receives strong disturbances, and forms a complicated wake field downstream of components such as the wing and fuselage, which is the main reason for affecting the aerodynamic characteristics of the rear plane.Through theoretical analysis [2,3], numerical calculations [4,5], wind tunnel tests [5,6], and flight tests [7,8], the aerodynamic characteristics of the rear aircraft at different relative positions in the formation flight can be obtained, thereby deriving the optimal position of the formation to maintain maximum aerodynamic gain.But there is another problem that should be concerned: how to arrive or leave the optimal position safely?When the rear plane is approaching, leaving, or passing through the wing-tip vortex of the front plane, the flow nearby changes significantly and rapidly, leading to nonlinear and unsteady aerodynamic characteristics, which can easily excite unexpected motion, reduce the reliability and effectiveness of flight control.In specific situations, maintaining a stable flying manually is quite difficult as fast and accurate operations are needed in an extremely short period of time, which will affect the flight safety and mission accomplishment.Therefore, it is necessary to conduct investigations on the dynamic effects during the entering process of the rear plane, which will support the verification for formation flight control strategies design, and contribute to the aviation safety.
The dynamic characteristics for an aircraft during certain flight task can be analyzed through aerodynamic modeling and simulation [9,10], based on the static aerodynamic coefficients obtained from engineering estimation or Computational Fluid Dynamics (CFD).However, unsteady, complicated vortices interaction is the dominant mechanism forming the aerodynamic interferences in a formation flight.The whole entering process should be simulated using time-accurate numerical method in order to take account of the unsteady flow phenomenon.
The use of CFD for the dynamic entering process of formation flight requires innovation of the dynamic grid method.Traditional schemes include the grid deformation, the grid reconstruction and the overlapping grid.Among which, the grid deformation method only suits for small-scale motion, while the dynamic entering usually covers distance many times larger than the wingspan of a plane; The grid reconstruction method is inefficient as the grid topology should be regenerated at every time step; The overlapping grid method requires cell scale consistency at the boundary between both set of grid to ensure a successful interpolation, hence the background grid should be refined throughout the entire path of the entering process, leading to a severely increase of the total number of grid cells (up to 100 million), which may be impossible to be solved using large clusters.
Therefore, a dynamic, adaptive overlapping grid method was developed to simulate the large-scale motion in a formation flight using limited server resources.Numerical simulations were conducted on the dynamic entering process of the formation combined by two transporters simplified from Ty-154.The effects of entering direction and relative speed on the aerodynamic characteristics were then analyzed.

Numerical method
An in-house CFD code was adopted, which is a three-dimensional, cell-centered Reynolds-averaged Navier-Stokes (RANS) solver based on unstructured grid.Some commonly used turbulence models such as S-A and SST models have been implemented.The multigrid scheme were used to accelerate convergence, and the parallel-running mode of the solver was available to improve the computational efficiency.The control equations, discretization schemes, turbulent models and some extra strategies on grid refinement and movement are illustrated in following subsections.

Control equations and discretization schemes
For the numerical simulation of unsteady flow, the RANS equations along the Cartesian coordinate system was considered.The equations take the following form: (1) where Q is the conserved flow variables, (E, F, G) are the inviscid convective flux, and (Ev, Fv, Gv) are the viscous flux: ( ) ( ) The spatial discretization adopted the Finite Volume Method (FVM) based on unstructured grids.The Roe's flux difference scheme was used for the inviscid convection term discretization.The thirdorder upwind biased upwind scheme for conservation laws (MUSCL) was used for flow variables interpolation at the interface.The Venkatakrishnan limiters were used to suppress numerical oscillations for the discontinuous problems.The viscous term of the equations was discretized by the second-order central difference scheme.
The boundaries of flow field can be divided into actual boundaries (fixed wall) and artificial boundaries (free stream, symmetry, and overset interpolation).The surface of planes was set to non-slip, adiabatic wall boundary condition, the free stream boundary was calculated using local one-dimensional Riemann invariants, the overset interpolation boundary should be generated dynamically by a series of operations including hole cutting, contributors searching and orphan removing.
Menter's Shear Stress Transport (SST) model was adopted to approximately simulate the turbulence.The model was first proposed by Menter in 1994 [11,12], which combines advantages of both k-ε and k-ω models.It acts equivalent to k-ω model near the wall, while gradually shifts to k-ε model away from the wall, showing high accuracy in general turbulent flow calculations.
The control equation of the SST turbulence model is as follows: 1 21 ( ) where ρ is density, k P is the production term, k is the turbulent kinetic energy, ω is the dissipation rate, μ is the laminar viscosity coefficient, μt is the turbulent viscosity coefficient.β * , β, γ, σk, σω, a1 are closure constants, and F1, F2 is weighting functions for adjusting constants in different flow regions.

Vorticity based adaptive mesh refinement method
The aerodynamic characteristics for the dynamic entering of the formation flight may undergo significant changes due to the dynamic interferences between flow around the rear plane and the wingtip vortex of the front plane, especially during the process when the rear plane is passing through the vortex.In order to get the detailed changing trends, the "small-scale" wake vortices should be simulated in high fidelity.
It is necessary to adopt a sufficiently fine grid where the wake flow of the front plane exists to improve the simulation accuracy.However, the fact that the distance between the front and rear planes is usually quite large (up to 3~10 times of the wing-span) in a formation makes it impossible to conduct the grid refinement globally.The grid refinement level that can be accepted by current clusters is insufficient to meet the high-fidelity simulations.Therefore, an adaptive mesh refinement method based on vorticity was adopted which could significantly reduce the grid cells to improve computational efficiency, when ensuring adequate simulation accuracy of for complicated wake vortex flow.
To achieve adaptive refinement, it is important to select appropriate flow field variables and establish the relationship between the variables and the grid cell scale.Then the bisection method was used for mesh refinement based on the scale distribution.After the alternating iteration of the flow control equations solving and the refinement, adaptive grid can be finally obtained for calculation.The selection of flow field variables is related to the physical problems to be simulated, for example when simulating mass separated flow of a plane, the turbulent kinetic energy can be selected, and when it comes to shock capture in supersonic flow, the Mach number gradient is a proper characteristic variable.For the simulation of formation flight, the simulation accuracy of wake vortices propagating downstream of the front plane is dominant, so the vorticity along the x-direction ωx was selected as the variable for adaptive mesh refinement, through which the grid cell scale L was defined as: Where Lmax and Lmin are the maximum and minimum cell scales, k is the tuning factor.A larger k produces finer grid, thus more grid cells.

Adaptive assembly for overlapping grid
Besides the high-fidelity simulation on "small-scale" wake vortices, it is also necessary to simulate the "large-scale" motion which covers the distance several times of wing-span, greatly beyond the characteristic length of the planes.The automatic overlapping grid method was adopted in the simulation.Generally, when variables are interpolated between different computational regions at boundaries, it is important to ensure that the cell scales on both sides of the boundary are equal, otherwise the inconsistency will lead to interpolation accuracy loss or contributor searching failure.While, the length of grid cells around the rear plane is much smaller than that of the background grid.One option is refining all cells of the background grid covered by the entering path of the rear plane manually, which will sharply increase the total number of grid cells and slow the computational jobs.Therefore, an adaptive assembly method for the overlapping grid was adopted, which automatically refines the adjacent background grid cells based on the overlapping boundary nodes distribution to ensure successful interpolation at each time step.
The kernel of the method is to how to refine the background grid robustly and efficiently.Two criteria can be chosen to decide whether the refinement should be accomplished: One is making an attempt to build the interpolation templates between overlapping grids, refining the adjacent cells of background grid when the attempt fails; The other is making a comparison between adjacent cells of both zones, the refinement finishes only when their lengths are equal.Conducting overlapping grid interpolation includes hole cutting, contributors searching and orphan removing, which is time consuming and unnecessary when cells of one side are much larger than that of the other side.So the adjacent background grid cells were bisected without the interpolation building until the cell lengths of both sides of grid zone are equal.Then the interpolation was tried and the cells were refined again when occurring failure.This process will be repeated until the grid assembly succeeds and available grid sets are formed.During the entire unsteady simulation, the above operation was conducted at each time step as the relative position of two planes may be changed.

Model introduction
The formation composed of two planes simplified from Ty-154 transporter was studied.The plane includes fuselage, wings, vertical tail, and high flat tail.Both planes are scaled by 1:22.Table 1 provides the geometric parameters.The appearance of the plane is shown in Figure 1.

Flow field region and computational grid
The simulation was carried out using a semi-zone grid, with symmetric boundary applied at the longitudinal symmetric plane of the front plane.The flow field was divided into two regions: the background zone and the rear plane zone, shown in Figure 3.The freestream vector was along the xaxis in grid coordinate system and parallel to the symmetric boundary condition.Thus the case simulated a Λ-shaped formation with one leader and two followers in fact.
The grid close to the surface of the plane has been refined along the normal direction to simulate the boundary layer, and the first grid spacing is about 0.0003 mm to ensure that y + =1. Figure 4(a) shows the surface grid distribution of the rear plane, with local refinement at surfaces with large curvature, to the minimum value of about 0.4 mm. Figure 4(b) shows the grid distribution at wall surfaces and longitudinal sections of the rear plane wing.It can be seen that the background grid has been adaptively refined at the overlapping boundary, ensuring grid cells consistent.The overlapping assembly result was shown in Figure 4(c).The minimum spatial grid spacing away from the solid wall is set to Lmin=10mm to achieve a balance between the simulation accuracy and the efficiency.

Adaptive grid independence
A grid independence study was conducted on the formation flight process with relative positions (x, y, z)/b = (3.0,-0.1, 0.0), to determine the appropriate value of the adaptive tuning factor k. Selecting k = 1.1, 1.3, 1.5 and 1.7 respectively, Figure 5 shows the trend of the total number of grid cells during the iteration.The initial grid has 6.07 million cells, and the converged number of grid cells reached 7.04 million, 9.94 million, 12.63 million, and 14.51 million respectively by different k.Table 2 shows the impact of k on the aerodynamic characteristics of the rear plane.Compared to results obtained by the finest grid when k = 1.7, the relative deviation of the longitudinal aerodynamics obtained by k = 1.5 is only about 0.2%, and the deviation of the lateral aerodynamics is about 1.2%~3.5%.While for k = 1.3, the deviations are relatively larger, with the longitudinal aerodynamics error of about 1.5%, and the lateral aerodynamics error of 10%.Therefore, it can be considered that when k = 1.3 the grid refinement is insufficient and the discretization error is non-negligible.When k ≥ 1.5, the calculations show results of general convergence.To improve computational efficiency, k = 1.5 was selected for adaptive grid refinement control in following simulations.

Numerical validation
A comparison with experiment result has been conducted for validation of the numerical method and the grid.The experiment which adopted a formation of two same aircrafts was made in a transonic wind tunnel with the test section of 2.4m×2.4m.The wingspans of the leading and the following models are 0.735m and 0.49m respectively, the angles of attack were 2° and 2.4°.The Mach number for the experiment was 0.76, and the Reynolds number based on the wingspan of the rear plane is 8.0×106.The relative streamwise location x/b = 1.5.More experimental details can be checked in studies by Tao Y [5].
Figure 6 shows the comparison between CFD and experiment for lift force coefficient CL, pitching moment coefficient Cm and rolling moment coefficient Cl.The computational and experimental results showed similar distribution along spanwise and vertical directions.The magnitudes of aerodynamic interferences obtained from both methods were consistent with each other, with some small errors that still existed because of the complexity of the vortex flow and the interferences of wind tunnel walls and struts.Since there were only distinctions of numerical parameters of test models and the incoming flow, between the validation case and the problems to be studied in present paper, the comparative study is believed to be appropriate to validate the numerical accuracy, and the computational method and grid have been proved to be valid.

Calculation results and analysis
The following numerical simulations were conducted at a freestream Mach number of 0.74.The Reynolds number based on the wingspan of the aircraft is 1.8×10 8 , thus the approximate flight altitude is 12km.The turbulence energy of freestream is 0.1%, and the turbulent-to-laminar viscosity ration is set to 50.The angles of attack are 2° for both planes.Different relative speeds have been analyzed including low, medium and high speed.

Longitudinal entering
The dynamic process of the rear plane entering the wake vortex influence zone of the front plane longitudinally is analysed in present section.Table 3 shows the initial and final formation location parameters, as well as the relative entering speed of the rear plane.Different levels of relative speed were derived from the moving distance and the accomplished time which is 1s for high speed, 2s for medium speed and 3s for high speed.Considering the aircrafts are scaled by 1:22, the specific value of speed for real flight should be 57.86 m/s, 86.68 m/s and 173.58, respectively, which varies in a considerably wide range.
Before the unsteady simulation, an initial flow field result should be calculated first.Thus the start time of the moving process of the rear plane was set to t = 0.5s, when the flow pattern and aerodynamic coefficients were converged.
Figure 7 shows the grid distribution on a longitudinal slice around wingtips at different time during the simulation.Before simulation, the background grid is in its initial state, only cells adjacent to the rear plane have been refined.At t = 0.5s, the background grid has been adaptively refined along the wake vortex propagation area of the front plane.The adaptive refinement area changes with the movement of the rear plane and the influence area of vortex wake flow.From the grid distribution changes, the adaptive refinement method developed in this paper can effectively simulate dynamic entering problem. Figure 8 shows the changes of the aerodynamic coefficients of the rear plane along with the moving distance when it enters the formation longitudinally.The impact process of the wing tip vortex on the rear plane can be observed though the change of lift coefficient.As the rear plane gradually moves forward and upwards, the lift coefficient gradually increases and reaches its maximum value when L = 5.12m ~ 5.32m, where the lift coefficient increases by about 0.058 compared to the initial state.It can be considered a zone where the wing tip vortex has the greatest effect on the rear plane due to the maximum lift, which can be named as the Theoretical Maximum Profit Zone (TMPZ).It can be used as a reference for the analysis of other aerodynamic coefficients and marked with orange coloured bands.The drag coefficient shows an increasing -decreasing -increasing trend, with TMPZ located in the reduction range.The pitch moment coefficient gradually decreases to a peak value near TMPZ and then rapidly increases, reaching its maximum value after crossing TMPZ, and then gradually decreases to the initial value.The vortices wake flow also has a significant impact on the lateral characteristics of the rear plane.When approaching TMPZ, the lateral force coefficient slowly increases to a maximum value and then rapidly decreases to a minimum value before returning to zero.At the same time, the trend of yaw moment coefficient is opposite to that of the lateral force coefficient, indicating that the lateral force mainly acts on the tail wing.The rolling moment coefficient increases rapidly to its maximum value before approaching TMPZ, and rapidly decreases to a smaller negative value at the beginning of TMPZ, followed by a slow increase.Generally, despite the maximum lift gain in TMPZ, there are significant changes in lateral force and aerodynamic moments, posing a huge challenge to the stability control of the plane.At the same time, the aerodynamic coefficients of the rear plane show same trend at different entering speed, implying that the aerodynamic hysteresis effect is not the dominant factor within the study range.Figure 9 shows the surface pressure distribution and wake interference at different time points.During the entering process, the wingtip vortex wake of the rear plane shows a downward pattern due to the relative velocity.When the wingtip vortex of front plane is far from the rear plane, there is no significant effect on the pressure distribution of the rear plane.When the frontal wingtip vortex sweeps over the wing of the rear plane, a downward wash flow is generated on the left side of the vortex core, and an upward wash flow is generated on the right side, resulting in a decrease in the local angle of attack on the outer side of the left wing, and an increase of angle of attack on the inner side.Therefore, the low-pressure area on the leeward side is enlarged on the inner side and reduced on the outer side.As most part of the rear plane is immerged in the upwash flow induced by the frontal wingtip vortex, the lift generated by the left wing increases at this time.When the rear plane is located below the wingtip vortex, the induced wash flow has a rightward velocity component, resulting in a positive lateral force coefficient mainly acting on the vertical tail, which then induces a negative yawing moment.However, when the rear plane is located above the wingtip vortex, as the induced wash flow has a leftward velocity component, a negative lateral force and a positive yawing moment are brought out.The induced rolling moment is also generated by the contribution of unbalanced aerodynamic forces.When the rear plane approaches the front wingtip vortex, the lift on its left wing gradually increases, and the vertical tail generates a lateral force to the right, causing a rapid increase in the rolling moment; As the rear plane gradually leaves the wingtip vortex, the lift on the left wing gradually decreases, while the lateral force of the vertical tail gradually turns negative, resulting in a rapid decrease in the rolling moment.

Lateral entering
Entering the formation from different directions implies completely different flow mechanisms and aerodynamic characteristics.The present section investigates the dynamic lateral entering process of the rear plane.Table 4 shows the initial and final formation parameters, as well as the corresponding relative speeds.Similar to the longitudinal process, the specific value of speed for real flight of the lateral entering process should be 53.02m/s, 79.42 m/s and 158.84.
Figure 10 gives the changes in aerodynamic coefficients during the entering process.From the trend of lift coefficient, the maximum value was reached in the range of L = 6.15m ~ 6.60m, which can be defined as TMPZ of the lateral entering.It can be seen that the change of drag coefficient during the lateral entering is greater than that during the longitudinal entering, and the pitching moment shows a decreasing trend at first and then rapidly increases.The lateral force gradually decreases during the movement, and correspondingly the yawing moment gradually increases.Before entering TMPZ, the rolling moment rapidly increases at first, similar to the trend of lift coefficient, and then rapidly decreases to a negative value.In both longitudinal and lateral entering process, although in TMPZ the maximum lift gain can be easily achieved, the unfavourable interferences of lateral forces and aerodynamic moments are significant, too. Figure 11 shows the flow details during the lateral entering process of the rear plane.From the local perspective of the rear plane, the frontal wingtip vortex gradually approaches the left wing tip and sweeps over the wing surface until it reaches the wing root.During the process, it is evident that the pressure on the leeward side of the left wing decreases globally at first, influenced by the induced wash flow, and then shows unbalanced distribution as the vortex sweeps over the wing.As the area affected by the vortex gradually moves to the right side, the lift on the left wing first increases and then decreases, while the lift on the right wing gradually increases at a slower speed, resulting in a trend of increasing and then decreasing of lift coefficient and rolling moment of the plane.The leftward washing flow induced by the vortex near the vertical tail of the rear plane also causes a gradually decreasing lateral force and correspondingly increasing yawing moment.Comparing the dynamic processes of different entering directions, it is concluded that the amount of vortex induced aerodynamic forces and moments of the rear plane are basically equivalent.However, during the longitudinal entering process, the directional stability of the rear plane changes severely, while during the lateral entering process, the longitudinal and lateral stability changes significantly.

Conclusion
The numerical study on the aerodynamic performance of the rear plane entering the formation of Ty-154 has been conducted.The vortices interference and its effect on aerodynamic coefficients were analysed.The results showed that: • (a) During the entering process, there exists a Theoretical Maximum Profit Zone (TMPZ), in which the maximum lift coefficient gain of about 0.05~0.06can be obtained; • (b) In the TMPZ, the lateral force and aerodynamic moments of the rear plane changes significantly, posing a huge challenge to the stability control, based on this fact, the TMPZ should not be recklessly pursued for safety in formation flight, strict security border should be set and certain buffer domains are equally needed because of the complex meteorological conditions such as the sudden gust; • (c) Entering the formation from longitudinal direction causes greater changes in directional stability, while from lateral direction causes greater changes in longitudinal and lateral stability, so the entering path should be designed based on the control capability of the following aircraft; • (d) Within the present study envelop, the aerodynamic hysteresis effect caused by relative speed has no significant impact on aerodynamic characteristics.

Figure 1 .
Figure 1.Geometric outline diagram of Ty154 (Unit: m).The formation location parameters are defined by the relative position between wing-tips of adjacent side of two planes, and are nondimensionalized using the wingspan of the plane.When the rear plane is located downstream, upper, and right of the front one, all the location parameters of x/b, y/b, z/b are positive.Figure 2 shows a schematic diagram of the relative position of the formation.

Figure 2 .
Figure 2. Schematic diagram of formation location parameters.

Figure 5 .
Figure 5.Total number of grid cells in adaptive refining process, with different k.

Figure 6 .
Figure 6.Aerodynamic coefficients distribution on the slice at x/b = 1.5, obtained by CFD (left) and wind tunnel experiment (right).

Figure 7 .
Figure 7. Grid distribution during the entering process, at different time step.

Figure 8 .
Figure 8. Aerodynamic coefficients of rear plane during longitudinal entering.

Figure 10 .
Figure 10.Aerodynamic coefficients of rear plane during lateral entering.

Table 2 .
Calculation results of different grid scales.

Table 3 .
Setups of longitudinal entering process.

Table 4 .
Setups of lateral entering process.