Modelling of Temporal‐Spatial Distribution of Airplane Wake Vortex for Scattering Analysis

Aircraft wake vortex is a pair of intensive counter-rotating airflow generated by a flying aircraft. Wake vortex is one of the most dangerous hazards in aviation because it may cause a following aircraft to roll out of control, particularly during the taking off and landing phases. The real-time detection of wake vortex is a frontier scientific problem emerging from many fields like aviation safety and atmospheric physics, and the dynamics and scattering characteristics of it remain as key problems to develop corresponding detection technologies. This chapter aims at presenting a simulation scheme for the dynamics of wake vortex under different weather conditions. For wake vortex generated in clear air, changes of the atmospheric dielectric constant produced by the density variation and water vapour variation are analysed; for wake vortex generated in rainy condition, the raindrop distribution in the wake vortex is also analysed. Both of them are essential for further analysing the scattering characteristics and developing new detection algorithms.


Introduction
Wake vortex is an inevitable physical phenomenon that exists in the rear zone of a flying aircraft, which rotates intensively and has a complex structure.The wake vortex generated by a large aircraft could be very hazardous to aviation safety since it might cause a following aircraft to roll out of control, particularly during the departure and landing phases.
In Air Traffic Management (ATM) field, International Civil Aviation Organization (ICAO) established a series of flight interval rules.These rules can ensure the flight safety in most time, but they are too conservative.In order to reach a good balance between avoiding the encountering hazard of wake vortex and increasing the transport capacity of airports, much attention has been paid on the real-time monitoring and detection of wake vortex in the past decades.Some major ATM programmes like Single European Sky ATM Research (SESAR) and Next Generation Air Transportation System, USA (NextGen) have also launched many projects on this topic, and the representative research institutes include Thales, Office National d'Etudes et de Recherches Aerospatiales (ONERA), Deutsches Zentrum für Luft-und Raumfahrt e.V. (DLR), Université catholique de Louvain (UCL), National Aeronautics and Space Administration (NASA), Federal Aviation Administration (FAA), Lincoln Lab, Boeing, and so on.In all these studies, the characteristics, detection technology and parameter retrieval are the key issues, and the characteristics of wake vortex serve as the basis for the rest studies.
In aviation safety, we mainly concern the clear air condition and wet weather condition.According to the scattering theory, the scattering of wake vortex in clear air is mainly determined by the fluctuation of dielectric constant inside the wake; while under wet condition, the key scattering factor becomes the massive number of precipitation droplets carried by the velocity of wake vortex.In this chapter, we present simulation schemes for the dielectric constant distribution and droplet distribution of wake vortex.The distributions are caused by the dynamics of wake vortex and serve as the physical basis for scattering analysis.
First, we study the dielectric constant distribution of wake vortex generated in clear air.

Two key parameters for determining the dielectric constant of wake vortex
The relative dielectric constant of atmosphere (ε r ) can be well depicted by the following expression [1]: where p,T,q are the pressure (pa), absolute temperature (K), and water vapour content (kg/kg), respectively.Generally, the second term in the square bracket is much smaller than 1, so the variation in dielectric constant between wake vortex and ambient air can be approximated as follows when the Taylor expansion is taken into account: In the expression, the parameters with subscript "a" refer to the ambient parameters and those without "a" refer to the wake vortex's parameters.
As is known, an isentropic process is a process in which there is neither heat exchange nor any friction effect [2].Typically, the wake vortex of a subsonic airplane can be assumed as an isentropic flow, and the thermodynamic parameters at different points along a streamline can be written as follows: where γ ¼ 1:4 is the adiabatic coefficient for air.Consequently, the variation in dielectric constant is transformed to Here we have denoted ξ ¼ ρ=ρ a , and the effects due to density and water vapour are mixed in the term ξ 2−γ q.In order to separate the two factors, this term is transformed to with Δq ¼ q−q a being the water vapour variation between the local wake and the ambient air.Since the variations in density and water vapour for wake vortex are very small, say ξ−1 << 1 and Δq << q a , the dielectric constant variation can be approximated as follows when the Taylor expansion is adopted: In this expression there are two undetermined parameters, ξ and Δq.They are separated into two different terms: The first term is determined by the density variation ξ: and the second term is determined by the water vapour variation Δq: In this manner, the key of modelling the dielectric constant is to determine the two parameters, ξ and Δq.

Velocity field of wake vortex
When the stationary phase of a subsonic wake vortex is taken into account, the dynamics can be well characterized by the steady Lamb momentum equation [2]: where V is the total velocity of wake vortex, and Ω ¼ ∇ • V is the vorticity.In this expression, the velocity and density are separated, which makes it convenient to work out the thermodynamics parameters according to the velocity field.
Typically, when no cross-wind is considered, a stable-stage wake is composed of two contour rotating vortices of the same strength, and they descend at a velocity V d .In this manner, in a coordinate system descending with the vortex cores, the wake vortex is steady, and the velocity for a given point P can be written as (see Figure 1): where V L and V R are the velocities deduced by the left and right vortices, and V d is the descending velocity.In Figure 1, the parameter Γ is the circulation which defines the strength of the wake vortex.
• In the expression, the deduced velocity of each vortex can be presented by existing velocity profile models, such as Rankine model, Lamb-Oseen model, and so on.Among them, the Rankine model is a widely used one, and the corresponding tangential velocity (V θ ) follows [3]: where r is the distance of a given point to the vortex centre, r c ≈0:052b 0 the vortex core radius, and b 0 the vortex spacing.As shown in Figure 2, the velocity outside the vortex core is irrotational, while the inside part is with a uniform vorticity: Vortex Structures in Fluid Dynamic Problems • The descending velocity can be derived from the deduced velocity of a vortex core:

Stream function of wake vortex
When the velocity components are determined, the vorticity is finally obtained as In this expression, we have considered ∇ • V d ¼ 0, and HðrÞ is an identification function for the left vortex core region (C L ) and right one (C R ): As a result, we have where u and v are the velocity components in x and y directions, and ψ is the stream function: In this manner, the Lamb momentum equation for wake vortex is rewritten as Furthermore, integrating the above equation from a point r to infinity gives where we have considered V ∞ ¼ V d ŷ, and the total stream function ψ is with the four terms on the right-hand side being the up-wash flow at infinity, the stream function due to the left vortex and right vortex, and a constant.
For a Rankine vortex, the stream function follows Then ψ L and ψ R can be obtained by replacing the circulation Γ with −Γ 0 and Γ 0 , respectively.
Also, the constant C is chosen to meet ψðrÞ¼0 when r is on the vortex core boundary: Due to the impact between two vortices, the constant C has a very slight variation along the vortex core boundary.In practice, the average of ψ 1 ðrÞ along one vortex core boundary is chosen as the constant.
As a result, combination of the total stream function and velocity gives the distribution of integral I as shown in Eq. ( 20), which is then used to calculate the density distribution of wake vortex.

Determination of dielectric constant due to density variation
Substituting the isentropic relationship (3) into the item on the right hand of Eq. ( 20) gives which can then be transformed as Vortex Structures in Fluid Dynamic Problems For a wake vortex, the integral I has the same magnitude as V 2 θ (generally not larger than 1000).On the other hand, ρ a and p a are, respectively, on the magnitude of 1 and 10 5 .Therefore, the second term in Eq. ( 25) is much smaller than 1, and the Taylor expansion gives: Consequently, the effects of dielectric constant can be rewritten as

Effect of water vapour variation on the dielectric constant
Generally, the atmospheric water vapour inside the wake vortex can be modelled as a passive scalar, which is convected by the wake vortex velocity field and is governed by the convectiondiffusion equation [4]: where q is the water vapour concentration, D is the diffusion coefficient for an air/water vapour system, and V is the velocity filed of a wake.
In fluid dynamics, the Péclet number is a dimensionless number indicating the rates of convection and diffusion of a flow [5]: where L is the characteristic length, U is the velocity, and D is the mass diffusion coefficient.Generally, a flow is convection-dominated, if the Péclet number is large.In this study, the diffusion coefficient for air (D ¼ 2:42 • 10 −5 ) is very small, which leads to a big Péclet number and the flow is convection dominated.In this manner, neglecting the impact of diffusion leads to a simplified governing equation: As is known, the initial water vapour gradient is very important to characterize the equation.
Here the stratified model is adopted in this work qðr, tÞj t¼0 ¼ q a ¼ m q ðy−y 0 Þþq 0 (31) with q 0 and m q being the offset and gradient of water vapour content, respectively.Substituting the initial condition into the governing equation gives ∂qðr, 0Þ ∂t þðV Á ∇Þqðr, 0Þ¼V y m q (32) Therefore, Eq. (30) minus Eq. (32) leads to a new equation: with Qðr, tÞ¼Δqðr, tÞ¼qðr, tÞ−qðr, 0Þ being the water vapour variation.At the same time, the initial condition of Eq. ( 33) becomes Qðr, tÞj t¼0 ¼ qðr, tÞj t¼0 −qðr, 0Þ¼0; this is coincident with the physical image that the water vapour variation is initially zero.
Eq. ( 33) is a hyperbola differential equation, which is often numerically difficult to solve.In the present work, the leapfrog scheme is adopted to solve the target equation, and good convergence and stability are achieved.The scheme is as follows: with u and v being the velocity components in x and y direction, respectively.In the process, the Von Neumann method leads to the following stable condition [6]: where Δt is the time step, Δ is the minimum grid spacing, and V max is the maximum velocity in the wake vortex.
In addition, non-uniform grids and symmetric condition are used to reduce computational cost.On the one hand, the velocity distribution shows that the flow in and around the vortex core is relatively complex and the flow far from the vortex core is slowly variational.Typically, sparse grids are adequate to characterize the slowly variational zones, but complex zones require dense grids.In this manner, non-uniform grid scheme can be adopted to reduce the total number of grids.On the other hand, the wake vortex is symmetric, so only half of wake vortex needs to be computed.Overall much computational cost can be saved through the above scheme.
If the water vapour variation, Qðr, tÞ is solved from Eq. (33), then Δε v r can be immediately obtained according to Eq. ( 9).
In this above simulation, the moving coordinate system is also used.The dielectric constant distribution can be transformed into the stationary coordinate system if the following transform is used: where V d is the descending velocity, t is the evolutional time, y and y ′ are the coordinates in the moving coordinate system and stationary coordinate system, respectively.

Numerical examples
Here we give several numerical examples for the dielectric constant distribution due to different effects.The parameters of airplane and atmosphere are as shown in Table 1.
According to the parameters, the distribution of IðrÞ can be worked out, and the intensity of dielectric constant due to density variation can be obtained as shown in Figure 3.It is observed that Δε d r in the vortex cores are much larger than that outside, so the vortex core could present a big contribution to the scattering of wake vortex.
The dielectric constant due to water vapour variation (Δε v r ) can be obtained when the given partial equation ( 33) is solved, and Figure 4 presents Δε v r at t ¼ 40s.It is observed that the convection effect of water vapour results in a non-uniform laminar structure in and around the wake vortex cores, and these structures could be good contributor to the scattering in high frequency.
Consequently, the sum of Figures 3 and 4 gives the total distribution of dielectric constant (see Figure 5).Comparing the magnitude of Δε d r and that of Δε v r shows that Δε v r is much less than Δε d r , and Δε d r dominates the magnitude of Δε r .However, Δε d r and Δε v r have different structures; they make different contribution to the scattering in different frequency bands.Typically, the scattering of clear air wake vortex at a frequency lower than 100MHz is mainly determined by the density variation Δε d r ; otherwise, the water vapour variation Δε v r makes the major contribution.

Extrapolation of dielectric constant distribution
The extrapolation includes two parts: one is relate to the density distribution and another is related to the water vapour distribution.Vortex Structures in Fluid Dynamic Problems

Extrapolation of dielectric constant related to density distribution
As shown in Eq. ( 27), the key of density distribution is the integral IðrÞ; this could be obtained with the normalized parameters.
First, if the space is normalized by vortex separation b 0 , say r ¼ r=b 0 , the Rankine model can be normalized as with Ṽ θ being the normalized velocity: Other velocity profile models have similar expressions.Substituting the normalized velocity into the integral Eq. ( 20) gives Where the normalized integral ĨðrÞ¼ ð Cð Ṽ Á ∇ ˜Þ Ṽ Á ds is only related to the velocity model, and ∇ ˜¼ð∂=∂x ∂=∂ỹÞ T ¼ b 0 ∇, s ¼ðx, yÞ.
For a stably flying airplane, the lift force balances the weight, which leads to the following initial circulation expression: Consequently, substituting the circulation into the integral gives The dielectric constant related to density distribution is then rewritten as Δε d r ðb 0 rÞ≈1:552 • 10 −6 R½T a ðb 0 rÞþ4668q a ðb 0 rÞ γp a ðb 0 rÞT a ðb 0 rÞ Since the normalized integral Ĩ ðrÞ¼ ð Cð Ṽ Á ∇ ˜Þ Ṽ Á ds is only related the velocity model, the following relationship can be obtained when the different airplane parameters and air conditions are introduced: with r 1 ¼ b 0, 1 r and r 2 ¼ b 0, 2 r.

Extrapolation of dielectric constant related to water vapour distribution
The key of water vapour distribution is to solve the partial differentiation equation (33).
Define some related normalized parameters as follows: where V 0 ¼ Γ 0 =b 0 and t 0 ¼ b 0 =V 0 .Consequently, the target partial differentiation equation is rewritten as If the variable Qðr, tÞ is solved from above equation, the dielectric constant related to water vapour distribution becomes In this manner, the dielectric constant for different airplane and air parameters is extrapolated as where For stably flying airplane, we have t 0 ¼ ρ a b 3 0 V a =ðMgÞ, and then the relationship becomes This is a very simple relationship.
With the combination of two extrapolation formulae, the dielectric constant distribution is then determined.This can save a lot of computation cost when different airplane wake vortices are to be analysed.
Another condition we always experience is wet weather condition (fog, rain, and snow).
Here we mainly concern the rainy condition since it is the most common situation around airports.

Wake vortex generated in rainy weather
In still air, the raindrops fall vertically towards the ground and reach a constant terminal falling velocity.When an aircraft takes off or lands in rainy weather, the raindrops will be inevitably involved in the aircraft wake vortices.Raindrops' motion in wake vortex is modified by the vortex flow.This modification of the trajectories of the raindrops may induce changes in raindrops' number concentration and velocity distribution in wake vortices, therefore results in changes in the recorded radar signal.This section presents a modelling scheme for raindrops' motion and distribution within the wake vortex.

Parameterization of raindrops
In still air, a falling raindrop reaches its terminal fall velocity V T with the equilibrium between the inertial force and the drag force acting on it [7].A widely used exponential expression between V T (m/s) and the diameter D(mm) is given by [8] where α 1 =9 .6 5m / s ,α 2 =1 0 .3m / s ,α 3 = 0.6 m/s, and ðρ 0 =ρÞ 0:4 is a density ratio correction factor adjusting deviation of the terminal fall velocity due to the air density change with the fall altitude.
Drop size distributions (DSD) have been widely used by radar meteorologists as they are directly related to radar reflectivity [9].In the following analysis, a suitable model to describe the size distribution of the rainfall in Europe is adopted [10] NðDÞ¼N 0 D 2 e −∧D (49) where N 0 =64500R −0:5 (m -3 mm -3 ) with R (mm/h) being the considered rain rate, ∧ ¼ 7:09R −0:27 (mm −1 ), NðDÞðm −3 mm −1 Þ represents the number of raindrops of the diameter D per unit volume per unit diameter class interval.

The motion equation of raindrops in wake vortices
Typically, the diameter of raindrops disperses between 0.5 and 4 mm.Usually their Stokes number in wake vortex flow is approximate to 1, which makes their motion trajectories significantly differ from the streamlines of total velocity field.To obtain those trajectories, the motion equation of the raindrops is studied [11].
When a raindrop enters into the wake vortex flow, its movement is governed by where t is the time, a is the acceleration of the raindrop, F d is the fluid drag force acting on the raindrop, m p is its mass, and g is the gravity acceleration.For a raindrop moving with velocity v p in the fluid whose velocity field is u½z p ðtÞ, if its diameter D ranges from 0.5 to 4 mm, the drag force F d can be approximately considered in the Newton regime [12] and given by where z p ðtÞ denotes raindrops' position, δv is the relative velocity between the vortex flow and the raindrop, and C d is the fluid drag coefficient.The impact of air density variations in the vortex flow on C d can be neglected because the density of raindrops is much larger than that of air.Thus, C d for a raindrop of diameter D is derived by the equilibrium equation of its weight and the drag force when falling at terminal falling velocity in still air: where ρ w is the density of raindrops.Substituting Eqs. ( 51) and (52) into Eq.( 50), the motion equation of raindrops within wake vortices can be further expressed as The instantaneous position and velocity of raindrops can be obtained from the above equations, but it is not easy to obtain a simple and closed expression.Here, a fourth-order four variables Runge-Kutta algorithm is proposed to solve the equation of motion [11].

Examples of trajectories of raindrops in wake vortices
In still air, the raindrop falls along a vertical trajectory to the ground.In presence of wake vortices, the trajectory of a raindrop is depending on its diameter and the location where it enters into the wake vortex flow.Circulation is a very important parameter to characterize wake vortex since it describes the strength of wake vortex.For raindrops moving in the vortex flow, their motion characteristics, that is, trajectory and velocity also largely depend on the vortex circulation.For simplicity, only the impacted part of the wake vortex region is shown in Figure 6, where four sets of trajectories are illustrated for raindrops of diameters 0.5, 1.0, 2.0, and 4.0 mm.Each set of trajectories corresponds to a given circulation value of the wake vortices.For a given circulation, the corresponding trajectory and velocity of raindrops in the wake vortices are unique.Comparisons between different trajectories give the following conclusions: (1) the smallest raindrops are much more sensitive to the vortex circulation and (2) the motion characteristic of raindrops in wake vortices is representative of the vortex strength.
Vortex Structures in Fluid Dynamic Problems

Raindrops' distribution in wake vortices
A Doppler radar will be very sensitive to raindrops' motion and possibly enable the detection of wake vortices in rain.To better understand the impact of wake vortices on the raindrops' motion, it is necessary to develop the methodology to quantitatively analyse the raindrops' distribution in wake vortices.

Raindrops' number concentration
The box counting method is adopted here to quantitatively compute raindrops' distribution in wake vortices.For simplification, we consider the situation where the raindrops are falling into the two dimensional rectangular region of wake vortex in stable phase.Before entering into the wake vortex region, the raindrops are falling in still air with the constant terminal falling velocity, and they are named as "initial raindrops".The number density of initial raindrops of diameter D (mm) is assumed to be N 0 ðDÞ (m −3 mm −1 ).In the wake vortex region, the raindrops' trajectory and velocity are changed and they are denoted as "disturbed raindrops".The number density of disturbed raindrops is assumed to be NðD, x, yÞ (m −3 mm −1 ), where (x, y) are the coordinates in the wake vortex region.Obviously, for a given wake vortex pair, NðD, x, yÞdepends on both the diameters of raindrops and their locations in wake vortex.In order to better illustrate the influence of wake vortices on the raindrops' distribution, the raindrops' relative number concentration at a point (x, y) is defined as η N ðD, x, yÞ¼ NðD, x, yÞ N 0 ðDÞ (54) where η N ðD, x, yÞ depicts the change in raindrops' concentration induced by the wake vortex.If η N > 1, the concentration of raindrops is enhanced, otherwise it is reduced.
In order to obtain the quantitative estimation of η N , the wake vortex region is divided into n x • n y grid boxes with equal size.The size of each grid box in xy plane are Δx and Δy, respectively.Above the wake vortex region, there are n x • 1 grid boxes where the initial raindrops of diameter D are homogeneously distributed, and the number of initial raindrops in each grid box is recorded as Num 0 ðDÞ.At each time step, their positions and velocities are updated by computing the equation of motion.If some of the initial raindrops enter into the wake vortex region, the same number of new initial raindrops is added to the first row of the n x • n y grid boxes, say, the n x • 1 grid boxes above the wake vortex region.When all the raindrops released at initial time arrive at the bottom of wake vortex region, the number of disturbed raindrops NumðD, x, yÞ in each grid box in wake vortex region is counted.Thus, η N ðD, x, yÞ can be approximated by the  In Figure 7, the raindrops' number concentration in wake vortices is illustrated.The simulation parameters are listed in Table 2.It can be noticed that in the wake vortex region, there are two columns where the raindrops' concentration is very small, even to zero.These two columns appear symmetrically below the two vortex cores and the distance between them in (a) is much wider than the others.Between these two columns, there are two narrow regions where the number concentration of raindrops is enhanced.For the raindrops of 1 and 2 mm of diameter, the number concentration value exceeds 8 in some grid boxes.

Raindrops' velocity distribution
Besides the number concentration, the raindrops' velocity distribution in each grid box is of great interest.If the grid box size used for the box counting method is small enough and the number of raindrops in each grid box is large enough, the velocity components of the raindrops in one grid box can be thought to obey Gaussian distributions.The mean value and variance of the velocity of the raindrops in each grid box are computed.If the number concentration of the grid box is zero, the raindrops' velocity in this grid box is set to 0.
In Figures 8 and 9, the raindrops' horizontal and vertical velocity distribution in wake vortices are illustrated, respectively.From Figure 8, it is interesting to find that the raindrops'

Figure 3 .
Figure 3. Dielectric constant distribution due to density variation.

Figure 4 .
Figure 4. Dielectric constant distribution due to water vapour variation.

Table 1 .
Parameters of airplane and atmosphere.