A Lagrangian Description of Buoyancy Effects on Aircraft Wake Vortices from Wing Tips near a Heated Ground Plane

Abstract

The present paper introduces the key ideas of a purely Lagrangian temperature particle method, which includes preheating effects on fluid flow nearest a heated wall. The numerical approach is then applied for the study of mixed heat transfer on aircraft wake vortices from wing tips in the vicinity of a heated ground plane, a situation commonly found during landing or takeoff operations at airports around the world. It was found in the literature experimental results of an investigation without the effects of heat transfer and crosswind, which were useful for a comparison with some present numerical results. Other numerical results are also discussed, focusing on the physics of the effects of mixed convection heat transfer and crosswind. As a contribution, the Richardson number is defined in terms of both aircraft wingspan and constant ground plane temperature, being the most important dimensionless group to capture the effects of laminar ascending mixed convection flow. The present methodology presents potentialities for predicting the transport and decay of primary vortical structures (under buoyancy forces), including their interaction with secondary vortical structures generated from a ground plane.

1. Introduction

The performance and reliability of airport traffic control around the world are directly related to the elapsed time reduction between landing or takeoff operations. A critical situation is found during simultaneous operations between parallel runways, where key factors, such as the ever-increasing airport traffic and aircraft sizes, cannot be disregarded. Further, tracking the aircraft viscous wake, which remains over the ground plane, taking some time to dissipate or to be removed by crosswind, directly affects the safety of landing and takeoff areas. It is still crucial to ensure that the flight of an airplane is not be interfered by the wake of another one [1]. According to Critchley and Foot [2], “accidents occur in subsequent operations, mainly in the 30~70 m range above ground level, when strong vortical structures are interacting with the runway ground”. In addition, Zheng and Ash [3] commented that “the free vortices, leaving the wing tips, have intensity proportional to the aircraft size and develop for considerable distances”.

In addition, Bobylev et al. [4] commented about the prohibitive flight problem near another one, provoking the risk of the vortex corridor and its consequences. With respect to aircraft wake vortices’ behavior in landing or takeoff operations, Vyshinsky [5] cited three characteristics to take into account: wake vortex dissipation, aerodynamic forces, and safe separation distance between flights near the airport. Vechtel et al. [6] reported that airports could use plate lines as a simple solution to reduce the vortex encounter severity, because plate lines have shown to accelerate the vortex decay and improve safety in airports.

To the best of our knowledge, aircraft wake vortices’ behavior near a ground plane was first studied assuming the incompressible inviscid flow model, where a fluid domain above the airport runway was adopted. In this model, the trajectory of two free vortices, shed from the wing tips, induces separation on the ground plane, but does not rebound [7]. Still, in the context of inviscid flow models, Donaldson and Bilanin [8] presented a thorough literature survey.

On the other hand, viscous flow models combined with crosswind effects enable the simulation of the boundary layer formed on the plane ground, being very useful because, substantially, they affect the aircraft wake behavior in accordance with the expected physics for the problem. Furthermore, the use of a pair of vortical structures shed from wing tips, instead of two free vortices, allows us to observe the primary vortical structure deformation and the development of secondary structures [9,10,11,12,13,14]. When analyzing the wake pattern near the ground, Qu et al. [15] reported that the viscous wake decay quickly responds because of the secondary vortical structures’ generation and the interaction between the primary and secondary vortical structures, being also slower under viscous effects.

Hallock and Holzäpfel [16] presented a review of wake vortex studies approaching the operational problem of spacing aircraft to improve airport capacity and operability. According to them, “although much has been learned about wake vortices and their behavior, there is still more to be learned about the phenomena of aircraft wake vortices”. The authors commented that the introduction of heavy aircraft, particularly B747-100 in 1969, led to aircraft wake separation standards and the beginning of extensive studies of vortical structure dissipation. The authors also observed that the use of extensive data and CFD has permitted the development of wake models to efficiently evaluate potential new procedures.

In the literature, there are studies discussing the problem of heat transfer in mixed convection linked with aircraft wake behavior near a ground plane [17,18,19,20,21,22,23,24,25]. In general, the effect of stable stratification on the descent of wake vortices is well known, since there is consensus about as primary vortical structures descend due to the mutual induction of their vorticity and temperature fields increasing buoyancy forces because of adiabatic compression. The descent rate with time of the pair of vortical structures is also reduced because of the positive buoyancy forces. The so-called vortex oval, i.e., the descending fluid volume involving the pair of vortical structures, also presents temperature differences along its surface, which generate counter-sign vorticity. However, there are unresolved issues, such as whether the spacing between the pair of vortical structures, shed from wing tips and rotating in opposite directions, increases or decreases with time, and whether there is a limiting descent altitude because of the stratification magnitude. Another issue related to atmospheric conditions is vertical shear of crosswind velocity profiles affecting the motion of wake vortices, where significant vertical changes of crosswind conditions can be found [26,27].

With respect to Lagrangian temperature particle methods, they present great potential to solve a variety of problems, such as shear layer, external flow, internal flow, and flow with multiple bodies. The key seed-corn paper in this field was presented by Smith and Stansby [28]. The authors used the vortex in a cell method in association with the random walk model [29] to introduce both discrete vortex and temperature particle elements according the similarity between vorticity transport equation and energy equation (forced convection heat transfer) in a two-dimensional flow. It is also important to mention the pioneering work of Ghoniem and Sherman [30], which investigated one- and two-dimensional diffusion heat transfer by using temperature particles associated with a random walk scheme to solve some specific problems. Kamemoto and Miyasaka [31] extended the discrete vortex element method based on the Biot–Savart law and core spreading scheme (diffusion solution) to study the time-dependent and forced-convective heat transfer around a bluff body in a uniform flow. The latter presented the basis to develop the Lagrangian approach used in this paper.

Ogami [32] originally proposed two models to create vorticity from temperature particles. First, he modeled the vorticity equation through of a natural extension of the approach presented by Ghoniem and Sherman [30]. In the second one (new), a discrete vortex element pair (one positive and one negative) was generated from each temperature particle. In the one-dimensional domain, both models were compared with an analytical solution, and the accuracy and validity were proven. In the two-dimensional domain, the application sample was also presented with the natural convection and the interaction between heat and vorticity. The first model proposed by Ogami [32] has been incorporated in the present Lagrangian approach because of low computational cost and very good efficacy.

The main contribution of this paper is to propose a viscous numerical model to investigate the problem of heat transfer in a mixed convection in two dimensions. The numerical method is based on the Lagrangian description. The focus is to investigate the influence of atmospheric effects on the development of aircraft wake vortices from wing tips near a ground plane. A model to evaluate diffusion heat transfer with a preheating effect between the wall plane and fluid is originally proposed here, aiming to reduce the final time of the numerical simulations. Literature records show that Lagrangian temperature particle methods with preheating models of fluid because of a heated wall effect are scarce. Therefore, the present paper contributes by proposing a preheating model to be incorporated into the Lagrangian approach, which simultaneously solves the vorticity transport equation and energy equation. The prediction of vortical structure trajectories is compared with experimental results [9], the latter with no heat transfer and no crosswind. Instantaneous vorticity contours and velocity field distributions are explored to show the formation and deformation of secondary structures. The influence of the Richardson number and crosswind conditions on primary vortical structure trajectories is investigated and properly discussed. The paper also discusses the potentialities of the present methodology for future applications focusing on surface roughness effects.

2. Problem Mathematical Formulation and the Vorticity and Temperature Fields

A schematic of a pair of vortical structures shedding from wing tips and rotating in opposite directions is shown in Figure 1, in which important quantities can be properly identified and defined. The primary vortical structures are characterized by the circulation strength Γ* = ±W*/(ρb*Ua*), W* being the aircraft weight, ρ the fluid density, b* the wingspan, and Ua* the aircraft approaching velocity. The fluid domain (Ω) includes the ground plane surface, that is, S = S₁∪S₂∪S₃, where the boundary conditions for the formulated problem must be imposed. The unsteady flow of a Newtonian fluid is assumed to be two dimensional, incompressible, and with moderate density differences (Boussinesq approximation). The crosswind is defined by the velocity U*_∞, the constant temperature far from the ground plane is defined by T*_∞, and the constant temperature at surface S is defined by T*_W. The runway length is B*, and the starting positions of the center of the vortical structures are defined by coordinates (B*/2 ± b*/2, h*/b*). The origin of the inertial frame of the reference is placed at x* = 0 and y* = 0.

In the general formulation of the problem, the quantities are nondimensionalized by the length scale b*. The dimensionless time is defined by t*U*_∞/b* (with crosswind condition) or t*Γ*/b*² (without crosswind condition). The temperature difference used as a temperature scale is ΔT* = T*_W − T*_∞. The symbol * identifies dimensional quantities.

The phenomenon to be investigated is governed by the continuity equation; the Navier–Stokes equations, where the Boussinesq approximation is assumed; and the energy equation. As the first methodology feature, the vorticity field is discretized and represented by an instantaneous vortex blob collection or Lamb discrete vortices [33]. The vorticity transport equation governs the motion of each vortex blob, and it is obtained by taking the curl of the Navier–Stokes equations, with additional support of the continuity equation [34]. Thus, the vorticity transport equation can be written as:

$$\frac{\partial \mathsf{\omega }}{\partial \mathrm{t}}+\left(\mathbf{u}\cdot \nabla \right)\mathsf{\omega }=\mathrm{R}\mathrm{i}\frac{\partial \mathsf{\theta }}{\partial \mathrm{x}}+\frac{1}{\mathrm{R}\mathrm{e}}{\nabla }^2\mathsf{\omega }$$

(1)

where ω is the only component of the vorticity vector in two dimensions, u is the velocity vector, Re is the Reynolds number, Ri is the Richardson number, and θ is the dimensionless temperature (see also Equation (5)).

The Reynolds number and the Richardson number are defined in the following forms, respectively:

$$\mathrm{R}\mathrm{e}=\frac{\mathsf{\rho }{\mathrm{U}}_{\infty }^{\ast }{\mathrm{b}}^{*}}{\mathsf{\mu }}=\frac{{\mathsf{\Gamma }}^{*}}{\mathsf{\nu }}$$

(2)

$$\mathrm{R}\mathrm{i}=\frac{\mathrm{G}\mathrm{r}}{\mathrm{R}{\mathrm{e}}^2}$$

(3)

where μ is the dynamic viscosity, ν is the kinematic viscosity, and Gr is the Grashof number.

The Grashof number represents the relationship between the buoyancy force and the viscous force, which can be represented in the following manner:

$$\mathrm{G}\mathrm{r}=\frac{\mathsf{\beta }{\mathrm{g}}^{\ast }\left({\mathrm{T}}_{\mathrm{W}}^{\ast }-{\mathrm{T}}_{\infty }^{\ast }\right){\mathrm{b}}^{{\ast }^3}}{{\mathsf{\nu }}^2}$$

(4)

where β is the thermal expansion coefficient and g is the gravity acceleration.

The dimensionless temperature is defined as:

$$\mathsf{\theta }=\frac{{\mathrm{T}}^{\ast }-{\mathrm{T}}_{\infty }^{\ast }}{{\mathrm{T}}_{\mathrm{W}}^{\ast }-{\mathrm{T}}_{\infty }^{\ast }}$$

(5)

where θ = 1 at the ground plane surface S and θ = 0 far from the ground plane (Figure 1).

As the second methodology feature, the heat is discretized and represented by an instantaneous temperature particle collection. The energy equation governs the motion of each temperature particle, and it assumes the following form:

$$\frac{\partial \mathsf{\theta }}{\partial \mathrm{t}}+\left(\mathbf{u}\cdot \nabla \right)\mathsf{\theta }=\frac{1}{\mathrm{R}\mathrm{e}\mathrm{P}\mathrm{r}}{\nabla }^2\mathsf{\theta }$$

(6)

where Pr is the Prandtl number defined such as:

$$\mathrm{P}\mathrm{r}=\frac{\mathsf{\nu }}{\alpha }$$

(7)

where α is the thermal diffusivity coefficient.

The boundary conditions of the mathematical problem must be satisfied at the ground plane surface S and far from the ground plane (Figure 1). More details about the boundary condition imposition will be discussed in Section 3.

3. The Temperature Particle Method with Buoyancy Forces

The present Lagrangian temperature particle method has been developed and applied to analyze unsteady, incompressible, and vortical flows with heat transfer [24,31,32]. The algorithm is simple and based on flow physics. As previously commented in Section 2, the vorticity-containing regions are discretized and represented by an instantaneous vortex blob distribution. In an analogy, the heat-containing regions are discretized and represented by an instantaneous temperature particle distribution. The vorticity and heat move around as the time evolutes. The heat is advected by the vorticity, and the heat interacts with the vorticity when the buoyancy effects are computed by assuming the Richardson number, Ri, different from zero (Equation (1)).

In an Eulerian manner, it can be more difficult numerically to solve the ever-changing vorticity and heat fields because of the use of a fixed grid. Therefore, the present Lagrangian approach follows individually throughout the numerical simulation both the vortex blobs and temperature particles. The Lagrangian manner advantage is that the attention is only directed to the regions of high activities, which are the regions containing vorticity and heat. The Lagrangian tracking of the particles avoids imposing the faraway boundary conditions; that is, they are automatically satisfied. In Figure 1, therefore, it is not necessary to impose boundary conditions far from the wall plane representing the runway airport.

A practical approach to approximate the ground plane surface, depicted in Figure 1 as S = S₁∪S₂∪S₃, is to break it down into a finite number M of straight-line elements (or flat panels) of the length ΔS_i. Each flat panel has a center point, named “pivotal point”, where the boundary conditions of the mathematical problem must be imposed. The no-slip boundary condition is satisfied on each pivotal point by choosing a number M of “shedding points” above the pivotal points, where M vortex blobs (Lamb vortices [34]) must be generated during each time step, Δt. The impermeability condition along all the ground plane length is automatically satisfied by the vortex mirror image method [35]. Each vortex blob is characterized by a distribution of vorticity, ${\mathsf{\varsigma }}_{{\sigma }_0}$ (commonly called the cut-off function); a circulation strength, ΔΓ; a core size, σ₀ (initially with the center at a “shedding point”); and a spatial position, x_i [33]. The solution of the system of linear algebraic equations, via the least squares method, is the circulation strength vector, {ΔΓ}. The vorticity generated from each flat panel in a time stepping remains concentrated into a viscous core, σ₀, and it is regarded as “free” to undergo advection and diffusion processes, satisfying the vorticity transport equation (see Equation (1)).

In Figure 2, the dimensionless heat flux can be represented by a heat quantity, Q, of the uniform density, q, and described in the Cartesian coordinate system (x, y), such that the y−axis is aligned with the thermal flux direction and the x−axis is parallel with the ground plane. The heat flux induces a time-dependent temperature distribution on the fluid near a flat panel, where T_W* is the wall temperature (isothermal condition) and T_∞* is the temperature far from the ground plane. The dimensionless heat quantity density is defined by q = q*/(B*²ΔT*). Thus, the ith dimensionless elemental heat quantity from the ith flat panel in a time stepping is represented in the following form:

$$\mathsf{\Delta }{\mathrm{Q}}_{\mathrm{i}}=\mathrm{q}\mathsf{\Delta }{\mathrm{S}}_{\mathrm{i}}$$

(8)

In order to handle temperature particles, Equation (8) gives us the basis to define the i^th temperature particle strength, ΔQ_i. The heat flux from the surface S to the fluid domain during each time stepping is determined by the well-known Fourier’s law [36], such as:

$$\frac{\partial \mathrm{q}}{\partial \mathrm{t}}=-\mathsf{\alpha }\frac{\mathrm{d}\mathsf{\theta }}{\mathrm{d}\mathrm{n}}$$

(9)

where n is the normal direction to the ground plane.

According to Equation (8), in the Lagrangian temperature particle method, as developed by Los Reis et al. [24] and Kamemoto and Miyasaka [31], the heat flux over the ground plane is represented by heat quantity density elements on flat panels, each producing one “free” temperature particle of the strength ΔQ_i at every time stepping; this conception is sketched in Figure 2. In order to simplify the numerical procedure in the heat calculation, the nascent rectangular heat element (Figure 2a) is replaced by an equivalent temperature particle, which has a thermal core, σ_i, with the center at a “shedding point” above the pivotal point (Figure 2b). A similar explanation can be presented for the discrete vortex elements.

The boundary condition of the temperature constant T_W* must be ensured on each pivotal point during each time step by using an equivalent temperature mirror image method [32]. Each nascent temperature particle has an instantaneous temperature defined by T_i* [24,31]. As a consequence, the temperature particle strength is numerically computed from Equation (9) in the following discretized form [24]:

$$\mathsf{\Delta }{\mathrm{Q}}_{\mathrm{i}}=\frac{1}{\mathrm{Re}\mathrm{Pr}}\frac{{\mathrm{T}}_{\mathrm{W}}^{\ast }-{\mathrm{T}}_{\mathrm{i}}^{\ast }}{{\mathrm{T}}_{\mathrm{W}}^{\ast }-{\mathrm{T}}_{\infty }^{\ast }}\frac{\mathsf{\Delta }{\mathrm{S}}_{\mathrm{i}}^{\ast }\mathsf{\Delta }{\mathrm{t}}_{}^{\ast }}{{\mathsf{\sigma }}_{\mathrm{i}}^{\ast }}$$

(10)

As any temperature particle is identified by a thermal core, σ_i, and strength, ΔQ_i (Equation (10)), a Gaussian distribution of the relative temperature around the center of a j^th vortex blob or at the point ${\mathrm{y}}_{\mathrm{j}}^{}$ (for y > 0 in Figure 1) is computed as [24,31]:

$${\mathrm{T}}_{\mathrm{j}}={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{Z}}\frac{\mathsf{\Delta }{\mathrm{Q}}_{\mathrm{i}}}{\mathsf{\pi }{\mathsf{\sigma }}_{\mathrm{i}}^2}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\mathrm{r}}_{\mathrm{j}\mathrm{i}}^2}{{\mathsf{\sigma }}_{\mathrm{i}}^2}\right)}$$

(11)

where Z is the total number of temperature particles, r_ji defines the radial distance between the points jth and ith, and Tj is the temperature at the jth point induced by the cluster of temperature particles.

As a contribution of this paper, the ground plane can also be represented by an infinite heat density sheet, in which the analytical solution is presented in accordance with Equation (12). For the initial sheet strength of ΔQ, the heat flux consists of densities ±q/2 above and below the heat density sheet. The key idea is inspirited in the modeling of inviscid potential flow presented by Martensen [37]. Thus, as time elapses, the heat will diffuse in the y direction, resulting in a temperature distribution that satisfies the heat diffusion equation [38]. Similar to the analytical solution of the vorticity distribution [34], the temperature induced at the point y_j (for y > 0 in Figure 1) can be obtained through the following analytical solution:

$${\mathrm{T}}_{\mathrm{j}}^{}=\frac{\mathrm{q}/2}{\sqrt{\mathsf{\pi }\mathsf{\tau }}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\mathrm{y}}_{\mathrm{j}}^2}{4\tau }\right)$$

(12)

where τ is defined as dimensionless preheating time, in which the temperature particle generation is started for t = τ.

According to Equation (12), a heat quantity density, q, can be estimated by adopting y_j = 0.0 and T_j = T_W, and also fixing a value for τ. Equation (12) refers to the proposed model to compute preheating effects on fluid near the ground plane during a preheating time of τ.

The total temperature induced at the point y_j (for y > 0 in Figure 1) is then computed, such as:

$${\mathrm{T}}_{\mathrm{j}}={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{Z}}\frac{\mathsf{\Delta }{\mathrm{Q}}_{\mathrm{i}}}{\mathsf{\pi }{\mathsf{\sigma }}_{\mathrm{i}}^2}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\mathrm{r}}_{\mathrm{j}\mathrm{i}}^2}{{\mathsf{\sigma }}_{\mathrm{i}}^2}\right)}-{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{Z}}\frac{\mathsf{\Delta }{\mathrm{Q}}_{\mathrm{i}}}{\mathsf{\pi }{\mathsf{\sigma }}_{\mathrm{i}}^2}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\mathrm{R}}_{\mathrm{j}\mathrm{i}\mathrm{m}}^2}{{\mathsf{\sigma }}_{\mathrm{i}}^2}\right)}+\frac{\mathrm{q}/2}{\sqrt{\mathsf{\pi }\mathsf{\tau }}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\mathrm{y}}_{\mathrm{j}}^2}{4\mathsf{\tau }}\right)$$

(13)

where the first term of the right-hand side represents the temperature particle contribution, the second term represents the mirror image temperature particle contribution because of the ground plane, and the third term establishes the wall temperature contribution for the preheating of the fluid near the ground plane.

Chorin [29] proposed an algorithm that splits the viscous term in the vorticity transport equation, Equation (1), to separately solve the advection and diffusion process. Because of the similarity between the vorticity transport equation and energy equation, Equation (6), the advective–diffusive operator in the latter equation is also split in two steps. Thus, the advection of each particle is calculated using the following second-order Adams–Bashforth time marching scheme [39]:

$${\mathrm{x}}_{\mathrm{i}}(\mathrm{t}+\mathsf{\Delta }\mathrm{t})={\mathrm{x}}_{\mathrm{i}}(\mathrm{t})+\left[1.5{\mathrm{u}}_{\mathrm{i}}(\mathrm{t})-0.5{\mathrm{u}}_{\mathrm{i}}(\mathrm{t}-\mathsf{\Delta }\mathrm{t})\right]\mathsf{\Delta }\mathrm{t}$$

(14)

where u_i represents the vector velocity calculated according to the Biot–Savart law (vortex–vortex interaction), image method (vortex–mirror image vortex interaction), and crosswind (vortex–lateral wind interaction) contributions at the point occupied by the ith vortex blob. Similarly, u_i represents the vector velocity calculated according to the Biot–Savart law (temperature–vortex interaction), image method (temperature–mirror image vortex interaction), and crosswind (temperature–lateral wind interaction) contributions at the point occupied by the ith temperature particle.

According to the random walk method [29], the vorticity and heat diffusion processes are calculated through a displacement of each particle in the following form:

$${\mathsf{\varsigma }}_{\mathrm{i}}\left(\mathrm{t}\right)=\sqrt{\frac{4\mathsf{\Delta }\mathrm{t}}{\mathrm{R}\mathrm{e}\mathrm{P}\mathrm{r}}\mathrm{l}\mathrm{n}\left(\frac{1}{\mathrm{P}}\right)}\left[\mathrm{c}\mathrm{o}\mathrm{s}{\left(2\mathsf{\pi }\mathrm{Q}\right)}_{\mathrm{x}}+\mathrm{s}\mathrm{i}\mathrm{n}{\left(2\mathsf{\pi }\mathrm{Q}\right)}_{\mathrm{y}}\right]$$

(15)

where P and Q are random numbers between 0 and 1, and the Prandtl number, Pr, is not considered when vortex blobs are diffused into the fluid domain.

In other cases, when vortex blobs or temperature particles migrate to the interior of the ground plane, they are reflected from their paths.

If the temperature gradient in the fluid domain causes moderate density differences, then the fluid motion is set up by buoyancy effects. Thus, the intensity of each vortex blob is increased by the amount ΔΓ⁺ because of direct interpretation from Equation (1), resulting in the following expression [32]:

$$\frac{\partial \mathsf{\omega }}{\partial \mathrm{t}}=\mathrm{R}\mathrm{i}\frac{\partial \mathsf{\theta }}{\partial \mathrm{x}}$$

(16)

Here, the circulation definition is used [34], that is, $\mathsf{\Gamma }={\displaystyle {\iint }_{\mathrm{A}}\mathsf{\omega }\cdot \mathbf{n}\mathrm{dA}}$, and Equation (16) can be easily discretized in accordance with the following equation:

$$\mathsf{\Delta }{\mathsf{\Gamma }}^{+}=\mathrm{R}\mathrm{i}\frac{\partial \mathsf{\theta }}{\partial \mathrm{x}}\mathsf{\Delta }\mathrm{t}\mathsf{\Delta }\mathrm{A}$$

(17)

where ΔA is the area occupied by the vortex blob under analysis, defined by a core radius of σ₀, during each time stepping, Δt. The derivative in Equation (17) is solved by the central difference method and its result is represented in the following manner [24]:

$$\frac{\mathrm{d}\mathsf{\theta }}{\mathrm{d}\mathrm{x}}\left|\begin{array}{l}\\ \\ \mathrm{i}\end{array}\right.=\frac{{}_{{\mathsf{\theta }}_{\mathrm{i}+1}-{\mathsf{\theta }}_{\mathrm{i}-1}}}{2\mathsf{\Delta }\mathrm{x}}-\frac{\mathsf{\Delta }{\mathrm{x}}^2}{6}\frac{{\mathrm{d}}^3\mathsf{\theta }}{\mathrm{d}{\mathrm{x}}^3}\left|\begin{array}{l}\\ \\ \mathrm{i}\end{array}\right.-\frac{\mathsf{\Delta }{\mathrm{x}}^4}{120}\frac{{\mathrm{d}}^5\mathsf{\theta }}{\mathrm{d}{\mathrm{x}}^5}\left|\begin{array}{l}\\ \\ \mathrm{i}\end{array}\right.-\dots$$

(18)

In the present numerical method, Equation (18) is truncated in the first term, where Δx is the ith vortex blob core radius. Therefore, the cluster of temperature particles, their images, and the preheating effects (Equation (13)) are required to compute the dimensionless temperatures θ_i+1 and θ_i−1 in the neighborhood of the ith vortex blob under analysis.

The interaction between the particles presents a high computational cost, and these calculations are efficiently made in parallel computing (OpenMP) in Fortran. The present work adopts the same methodology proposed by Los Reis et al. [24], where for each thread is provided each instantaneous particle position and its intensity. Therefore, the interaction of an arbitrary particle can be computed and, as consequence, each computational core evaluates a different particle, reducing the total time of a typical numerical simulation. The present strategy has proved to be extremely efficient, providing high CPU usage and near-null thread idle time.

Summarizing Section 3, this paper proposes a novel Lagrangian algorithm for the temperature particle method, aiming to calculate the buoyancy effects in a fluid caused by an isothermal ground plane. The vorticity is discretized and represented by a cloud of vortex blobs (Lamb vortices). The heat is discretized and represented by a cluster of temperature particles. The numerical approach exactly imposes the impermeability boundary condition on the ground plane surface using the image method and, simultaneously, conserves global circulation generating vortex blobs. The latter utilizes a linear system of algebraic equations to obtain the unknown circulation strength of the nascent vortex blobs. The boundary condition of constant temperature on the ground plane is satisfied through the Fourier law adapted to the thermal boundary layer for the unknown heat quantity strength of the temperature particles. The numerical model contributes by proposing an analytical solution for preheating effects to be combined with the mirror (image) temperature particle method.

The buoyancy effects are computed by changing the strength of each vortex blob during each time stepping of the simulation through a direct interpretation of the incompressible vorticity transport equation linked with the Boussinesq approximation. A second-order Adams–Bashforth time marching scheme is used to advance the position of the vortex blobs and temperature particles into the fluid domain because of advection. The random walk method is used to simulate the vorticity and heat (molecular) diffusion.

Finally, the developed algorithm is implemented using an in-house code via parallel computing (OpenMP) in Fortran through 12 steps in accordance with the following sequence:

(I).

Ground plane construction including the number M of “shedding points” above the pivotal points.

(II).

Creation of the pair of vortical structures shedding from wing tips.

(III).

Preheating effect computation on “shedding points” using Equation (12).

(IV).

Simultaneous generation of new vortex blobs and temperature particles (the buoyancy effects are included, when Ri ≠ 0.0).

(V).

Velocity field computation at vortex blobs and temperature particle locations into the fluid domain (i.e., real plan, as sketched in Figure 1).

(VI).

Advection of the vortex cloud and of the cluster of temperature particles using Equation (14).

(VII).

Vorticity and heat diffusion using Equation (15).

(VIII).

Reflection of vortex blobs and temperature particles that accidently migrate into the ground plane (i.e., the imaginary plan, as sketched in Figure 1).

(IX).

Buoyancy effect computation using Equations (17) and (18), if Ri ≠ 0.0.

(X).

Velocity field computation at each pivotal point to satisfy the no-slip condition and, consequently, to obtain a new generation of vortex blobs.

(XI).

Temperature computation at each pivotal point using Equation (13) for new generation of temperature particles.

(XII).

Advance by the time Δt.

4. Results and Discussion

4.1. Initial Conditions and Validation

In order to perform the simulations, it is necessary to estimate some numerical parameters as organized in

Table 1. Estimated dimensionless values of numerical simulation parameters.

Parameter	Value
Runway Length (B)	8.0
Release height of the primary vortical structures (h)	2.0
Center-to-center spacing of the vortical structure pair (b)	1.0
Circulation strength of the vortical structures (±Γ)	1.0
Core radius of each primary vortical structure (r)	0.1
Temperature far from the ground plane surface	0.0
Temperature at ground plane surface	1.0
Time stepping (Δt)	0.05
Number of time steps (NSteps)	1000
Number of flat panels (M)	40
Number of Lamb vortex elements of each primary vortical structure (N)	50
Core radius of each Lamb vortex (σ0)	0.001
Core radius of each temperature particle (σ)	0.001
Preheating time (τ)	1.5

.

First of all, the time stepping of Δt = 0.05 was appropriately chosen by solving two particular problems, which present an exact solution for purposes of comparison. These problems are related to Equations (14) and (15), as will be reported bellow.

Equation (14) numerically estimates the advective motion of vortex and temperature particles by the second-order Adams–Bashforth time marching scheme. In the present numerical approach, temperature particles are always advected by vortex blobs. Figure 3 shows a vortex blob pair of equal strength Γ = 1.0, distance AB apart. The vortex blob identified as #1 is initially positioned at (x, y) = (0.0, 1.0), while the vortex blob #2 is placed at (x, y) = (0.0, −1.0). A temperature particle pair identified as #1 and #2 is also considered for the present numerical simulation. These temperature particles are initially positioned at (x, y) = (0.0, 1.0) and (x, y) = (0.0, −1.0), respectively. All particles into the fluid domain have a core size of σ = 0.001. A vortex particle will experience an advection velocity induced by the other and computed by the Biot–Savart law in the following manner:

$$\overline{\mathrm{u}{\mathrm{v}}_{\mathrm{k}}^{\mathrm{n}}}\left({\mathbf{x}}_{\mathrm{k}},\mathrm{t}\right)={\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{N}\mathrm{V}}{\mathsf{\Gamma }}_{\mathrm{j}}{\mathrm{c}}_{\mathrm{k}\mathrm{j}}^{\mathrm{n}}\left[{\mathbf{x}}_{\mathrm{k}}\left(\mathrm{t}\right)-{\mathbf{x}}_{\mathrm{j}}\left(\mathrm{t}\right)\right]},\mathrm{n}=1,2;\mathrm{k}=1,\mathrm{N}\mathrm{V}$$

(19)

where, for a generic vortex cloud with NV particles, ${\mathsf{\Gamma }}_{\mathrm{j}}$ represents the strength of the jth vortex particle, and ${\mathrm{c}}_{\mathrm{k}\mathrm{j}}^{\mathrm{n}}\left[{\mathbf{x}}_{\mathrm{k}}\left(\mathrm{t}\right)-{\mathbf{x}}_{\mathrm{j}}\left(\mathrm{t}\right)\right]$ establishes the nth component of the velocity that the jth vortex particle (of unitary strength) induces at the location of the kth vortex particle.

Particularly, in Figure 3, the advective velocity induced on each vortex blob due the other will act in the direction normal to the line AB joining them.

As the time runs, the vortex blob pair will process in a circular path about the midpoint of AB. Equation (14) computes the advective motions along the time of both vortex and temperature particles. The vortex blobs of equal strength Γ = 1.0 carrying temperature particles move in a clockwise direction. The numerical simulation is stopped after 790 “Steps”, when the vortex blob #1 approximately attains the position (0.0, −1.0). Concurrently, the vortex blob #2 approximately attains the position (0.0, 1.0). Thus, an estimation of the relative error between the true circular drift path of the constant distance AB = 2.0b and the numerical solution obtained by Equation (14) is about 0.0011%. Therefore, the time stepping Δt = 0.05 is suitable to be used in the second-order Adams–Bashforth time marching scheme. As an advantage here, it can be commented that the Lagrangian approach to solve advection motions of vorticity and heat does not need any explicit treatment of advective derivatives, as requested by Eulerian schemes.

Equation (15) computes the diffusive motion of vortex and temperature particles by the random walk method. The diffusion method convergence is presented by analyzing the radial diffusion problem of a vortical structure of the strength Γ_v = 1.0, which is initially centered on the origin (0.0, 0.0), as shown in Figure 4. For that problem, there is an analytical solution for the vorticity distribution in space and time, such that [34]:

$$\mathsf{\omega }\left(\mathrm{r},\mathrm{t}\right)=\frac{{\mathsf{\Gamma }}_{\mathrm{v}}}{4\mathsf{\pi }\mathsf{\nu }\mathrm{t}}{\mathrm{e}}^{\left(-{\mathrm{r}}^2/4\mathsf{\nu }\mathrm{t}\right)}$$

(20)

Figure 5 exemplifies the test case for t = 1.0, Γ_v = 1.0, and Re = Γ_v/ν = 1.0. For purposes of comparison, the numerical solution utilizes a vortical structure Γ_v discretized and represented by Z = 5000 vortex blobs of strength Γ_v/Z and of core size σ = 0.001. All vortex particles were initially distributed on the origin (0.0, 0.0) at t = 0, such that, after 20 time steps of the magnitude Δt = 0.05, they diffuse (see in Figure 4), obeying Equation (15). Thus, the vorticity distribution can be statistically evaluated in accordance with the following formula:

$$\mathsf{\omega }\left(\mathrm{r}\right)=\frac{{\mathrm{z}}_{\mathrm{j}}}{\mathrm{Z}\mathsf{\pi }\left(r_{j+1}^2-r_j^2\right)}$$

(21)

In Figure 4, 20 annular bins of the same thickness were adopted, that is, Δr = 0.4, and center at the origin (0.0, 0.0), which are suitable to capture all z_j vortex blobs into the correspondent annular area lying between 0.0 ≤ r_j ≤ 7.6 and 0.4 ≤ r_{j + 1} ≤ 8.0 (Equation (21)). As shown in Figure 5, the numerical solution agrees very well with the exact solution. As additional information, in Equation (20), the exact solution of the radial vorticity distribution utilized the r.m.s. radius $\sqrt{\left[\frac{1}{2}\left({\mathrm{r}}_{\mathrm{j}}^2+{\mathrm{r}}_{\mathrm{j}+1}^2\right)\right]}$ for every strip j of an annular bin. Since a similar solution can be obtained for radial heat diffusion, the time stepping Δt = 0.05 is also suitable for the random walk method.

An amount of N = 50 vortex blobs (discrete vortices) is used to represent each of the wing tip counter-rotating vortical structures, rather than the option of using a pair of simple free vortices that cannot simulate the deformation of these primary structures. Initially concentrated in a single point, the N = 50 vortex blobs that represent a primary vortical structure in sequence are spread via random displacement by establishing a Gaussian distribution confined by a circumference of a dimensionless radius equal to 0.1 [24]. As consequence, the initial configuration to represent each of the wing tip counter-rotating vortical structures is very similar to the representation in Figure 4.

Therefore, in all the simulations, a pair of counter-rotating vortical structures sheds from the position (B/2 ± b/2, h) with the circulation Γ = ±1 (Figure 1). These wing tip vortices can be submitted to crosswind velocities of U_∞ = 0.0, 0.02 and 0.04, and the Reynolds number chosen for all test cases is fixed as Re = 7650. Half the wingspan of the aircraft and the initial shedding altitude are b/2 = 0.5 and h = 2.0, respectively, once that configuration allows a comparison with the experimental data given by Liu and Srnsky [9] using Re = 7650, the latter in the absence of crosswind and buoyancy effects. The relation h/b = 2.0 is typical for problems of landing or takeoff operations between parallel runways, when an aircraft develops velocity about approximately 240 km/h.

The length of the airport runway, sufficient to establish the interaction between the runway boundary layer and the primary vortical structures shedding from the aircraft, is defined by B = 8b; that length is also typical for problems of landing or takeoff operations between parallel runways.

In simulations with the presence of mixed convection heat transfer, the fluid temperature far from the ground plane is fixed at T_∞* = 20 °C, and two temperatures are, strategically, adopted for the airport runway, that is, T_W* = 40 °C and 60 °C, respectively, resulting in the respective Richardson numbers of Ri = 1.018 and 2.036. In accordance with the present methodology, an increase in Richardson number will produce unrealistic temperatures on a ground plane (Equations (2)–(4)), as will be discussed in Section 4.2.

Finally, the Prandtl number is properly chosen as Pr = 0.7 (fluid is the air) for all the simulations presented in this paper, including thermal effects. As already commented and for preheating effects, the density q/2 is computed using y_j = 0.0, T_j = T_W, and τ = 1.5 (Equation (12)).

Figure 6 presents a comparison between the experimental results [9] and the numerical results (particularly for Ri = 0.0 and without crosswind). These results refer to the predicted trajectory of the centroid of the right-side vortical structure from the wing tip (Figure 1). It was found that the predicted numerical results are consistent with the experiment data, and the deviation is less than 12% for a wide range of the trajectory. In other words, the trajectory obtained by the numerical method adjusts, in general, to the experimental data in its descending and ascending condition until the beginning of the first loop. After that, there is no information about the sequence of the data obtained experimentally for comparison, which did not allow the calculation of the average relative error.

4.2. Buoyancy Effects

When a flow is in the same direction and orientation of the buoyancy forces effect, it is called aided flow; on the other hand, a flow in the opposite direction is called concurrent flow. Accordingly, Hussein [40] reinforced that during the mixed convection heat transfer, it is necessary to add the contributions of the natural and forced convection heat transfer in the aided flows and to subtract them in opposite flows. That treatment must be made for each vortex blob present in the fluid domain in the context of the present numerical method. Therefore, the natural convection heat transfer is a computed vortex blob by a vortex blob, whereas the forced convection heat transfer is linked to the flow velocity induction effect. Figure 6 also shows the trajectories of the centroid of the right vortical structure from the wing tip for other two test cases without crosswind, that is, with thermal effects (buoyancy forces) for Ri = 1.018 and Ri = 2.036.

As already mentioned in Section 4.1, in the numerical simulations carried out, the condition of both absences of buoyancy effects and crosswind is the reference test case for comparisons. Therefore, in addition, the buoyancy effects without crosswind are considered in the numerical simulations by adopting a Richardson number from order of 1, which means that natural convection and forced convection are equivalent and a Richardson number from order of 2, which represents the preponderance of natural convection heat transfer.

The addition of the heat transfer from the ground plane surface using Ri = 0.0 shows that there is practically no difference in the downward path and during the recovery of the right vortical structure from the wing tip after the interference of the two boundary layers from the ground plane (preheating effects are included in all test cases). On the other hand, there is a subsequent increase in the recirculation of the same primary vortical structure for mixed convection heat transfer conditions when compared with the reference case (Ri = 0.0). That behavior is already verified in the first loop for the test case with Ri = 2.036 and in the second loop for the test case with Ri = 1.018 (Figure 6).

For a heated ground plane, therefore, the buoyancy effects (Ri > 0.0) are sensitive and interfere with the development of the hydrodynamic boundary layer producing favorable pressure gradients [41]; consequently, an increase in the temperature gradient near the ground plane, because of the Richardson number increase, can be identified. Figure 6 illustrates a strong recirculation of the same primary vortical structure increasing to the left and with its lateral propagation to the right being inhibited, when is chosen Ri = 1.018. It can be mentioned that the buoyancy torque, resulting from movement within a stratified fluid away from the ground plane, is also inhibited by that recirculation [3].

It is interesting to note that the action of the primary vortical structures, after approaching the ground plane and coupling with the secondary structures generated on the same ground plane, can reach altitudes equivalent to the shedding altitude, although their centroid are limited to lower altitudes. This can be compared in Figure 7, for the condition without buoyancy effects (Figure 7a) and for conditions with buoyancy effects (Figure 7b,c). In the present investigated phenomenon, the mixed convection heat transfer, interestingly, is more effective for Ri = 1.018 (Figure 7b). In general, natural convection heat transfer effects are negligible for Ri < 0.1, forced convection is negligible for Ri > 10, and neither is negligible in the range 0.1 < Ri < 10. The test case for the Richardson number of the order of unity (i.e., Ri = 1.018) indicates that the flow is likely buoyancy driven.

4.3. Combined Effects of Buoyancy Forces and Crosswind

In this section, the main results are presented and discussed for the condition of buoyancy effects on aircraft wake vortices under the action of different crosswind velocity values.

It is well known that there is a discordant decay of the vortical structures from the wing tip, left and right, while their trajectories are altered by the crosswind shear as they approach the ground plane. In that situation, such an asymmetric movement is consolidated, and the right vortical structure (downwind) grows and has a greater lateral propagation, while the left vortical structure (against the crosswind) loses altitude in the upward recovery movement, which, obviously, depends on the intensity of the crosswind.

It can be observed in Figure 8a that the greater the crosswind, the greater the tendency of the left vortical structure to remain rotating on a more restricted region close to the ground plane and to lose altitude, immediately, after the first approach to the ground plane. Figure 8b shows the right vortical structure trajectories for all the crosswind velocities tested without buoyancy effects (Ri = 0.0). It can be seen that, as the crosswind velocity increases, the vortical structure moves even more to the right, as expected, and it rolls up higher.

For configuration without both crosswind and buoyancy effects, after interacting with the hydrodynamic boundary layer formed on the ground plane, the centroid of the right vortical structure makes an upward trajectory up to a maximum altitude of approximately y = 1.32 (Figure 8b). On the other hand, under influence of different crosswind velocity values, maximum altitudes around y = 1.41 and 1.63 are identified for velocities of U_∞ = 0.02 and 0.04, respectively (Figure 8b). It is evident that there is coupling between these primary structures and the secondary structure generated on the ground plane, which induces a more pronounced recirculation for the primary vortical structures. There is also a more pronounced drop in the second approximation of the vortical structures of the ground plane (x = 6.50 and y = 0.72) for a greater crosswind velocity value of U_∞ = 0.04.

Figure 9 presents the instantaneous vorticity field behavior when the vortical structures from wing tips are closer to the ground plane. It is possible to identify how the positive (clockwise) and negative (anticlockwise) vortical structures induce the generation of secondary vortical structures of opposite signs from the ground plane; it can be noticed that, as the crosswind velocity increases, the magnitude of the vorticity developed from the ground increases. However, Zheng and Ash [3] reported that the secondary vortical structures dissipate more quickly, presenting little influence on the primary vortical structures for low Reynolds number values.

The mixed convection heat transfer is now considered in the presence of crosswind velocities of magnitude U_∞ = 0.02 and 0.04 to allow investigations of these combined effects. For this purpose, as performed in Section 4.2, as previously presented, the Richardson number also assumes the values of Ri = 1.018 and 2.036.

The predicted trajectory of the centroid of the right-side vortical structure from the wing tip (Figure 1) in the presence of different crosswind velocity values, and considering the Richardson number of Ri = 1.018, is shown in Figure 10. The condition without crosswind is also presented, as a reference case for comparisons. Thus, the presence of thermal effects on each situation shown in Figure 8 provokes different consequences on them (see Figure 10). When there is no crosswind, it is clear that the thermal effect makes the right vortical structure roll up more, and in addition, it goes higher.

On the other hand, it is remarkable that the crosswind attenuates the influence of mixed convection heat transfer, since the trajectories of the right vortical structure are almost the same as those shown in Figure 8; however, they do not reach the altitudes verified in the absence of thermal effects (compare the trajectories predicted by the numerical method with crosswind in Figure 8b and Figure 10, for Ri = 0.0 and Ri = 1.018, respectively). For a crosswind velocity of U_∞ = 0.04, for example, in the first case (Figure 8), there is a maximum altitude of y = 1.63, and in the second case, a maximum altitude of y = 1.57 is observed (Figure 10). It is also verified that the heat transfer effects are not able to inhibit, even if moderately, the lateral propagation of the vortical structure on the right side in the presence of different crosswind velocity values.

Figure 11 illustrates the approximated instant where secondary vortical structures are being generated below the primary vortical structures, when the latter are close to the ground plane. These secondary vortical structures are weak. Further, asymmetric decay of the primary vortical structures has begun, which is governed by the magnitude of the crosswind; and the interaction with the heated ground plane implies in the injection of momentum modifying the vorticity and recirculation of these structures. It can be also observed that as the crosswind increases, the left primary vortical structure moves more to right than the right primary vortical structure. The reason for that behavior is related to the barrier caused by a stronger vorticity of an opposite signal generated from the ground plane (under crosswind effect) constructing the right secondary vortical structure.

Figure 12 presents the same test cases illustrated in Figure 11, after 1000 time steps, that is, when the numerical simulations are stopped. Clearly, the influence of the higher crosswind velocity can be visualized (Figure 12c), which moves further the primary vortical structures from left to right when leaving the runway ground.

The trajectories of the vortical structures in the presence of different crosswind velocity values, now considering a Richardson number of Ri = 2.036, are shown in Figure 13. When comparing the results synthesized in Figure 13b with those in Figure 10, the Richardson number increase interferes on predicted trajectories for the right-side vortical structure under different crosswind magnitudes. The emphasis here is about the maximum altitudes, in general, moderately larger for an increase in the Richardson number and for a greater limitation of the lateral propagation only for the crosswind of greater magnitude, that is, U_∞ = 0.04.

For the test case without crosswind, the maximum altitude established in Figure 13b is around y = 1.34, being moderately lower than that established in Figure 10 for y = 1.43. When analyzing the crosswind of magnitude U_∞ = 0.02, the maximum altitude established in Figure 13b is around y =1.43, which is moderately higher than that established in Figure 10 for y = 1.38. Finally, for the crosswind of magnitude U_∞ = 0.04, the maximum altitude established in Figure 13b is around y =1.60; that value is moderately higher than that established for y = 1.57 in Figure 10. On the other hand, the lateral propagation inhibition of the vortical structure is observed only for U_∞ = 0.04, as well as a greater deepening in its second approach to the ground, for an increase in the Richardson number (see the red curve in Figure 10 and Figure 13b).

With regard to the left vortical structure, when comparing its recirculation trajectories in combined conditions of thermal effects (Ri = 2.036) and crosswind influences (Figure 13a) with the trajectories consolidated only under the influence of crosswind (Figure 8a), moderate changes are observed for U_∞ = 0.0 and U_∞ = 0.02. However, in the presence of the crosswind of greater magnitude, U_∞ = 0.04, the influence of buoyancy forces is more intense and acts to restrict the recirculation and the lateral propagation of the vortical structure against the crosswind. A maximum lateral propagation and a maximum altitude are observed equal to x = 4.1 and y = 1.11 in Figure 13a and equal to x = 4.38 and y = 1.26 in Figure 8a, respectively (compare the trajectory represented in red color in the mentioned figures).

In practice, this kind of study helps to reduce the time elapsed between subsequent operations at airports. That time directly affects the safety of a landing or takeoff operation, since the pair of counter-rotating vortical structures from wing tips remains over the plane ground and takes some time to dissipate or be removed by a typical crosswind.

Therefore, the results described in this section demonstrate that the dynamic of a pair of wing tip vortices close to a heated ground plane is a complex phenomenon, dependent on the coupling of the change in the momentum (because of the heat transfer) of the vortical structures with the action of different crosswind velocities, and also with the interaction of the hydrodynamic boundary layer formed on the ground plane, even for two-dimensional conditions without turbulence modeling.

4.4. Perspectives of the Present Methodology for New Applications

In the last years, our research group has made an effort to develop an in-house code aiming to demonstrate that the effects of two-dimensional roughness are much more sensitive than a single turbulence modeling [42,43,44,45].

The main results published have captured important phenomena of nonlinear hydrodynamics. The early proposed methodology specifically injects momentum in the hydrodynamic boundary layer to capture delay in the separation point of flows past a bluff body [43]. The problem of interference between a circular cylinder and a moving ground has been successively explored, including important physical interpretations [44]. The viscous flow past two circular cylinders is another problem that has been investigated too [45].

In closing, the numerical approach presented in this paper is very suitable to be incorporated into the problems of interference between solid boundaries mentioned above. Results available in the literature discussing combined effects of surface roughness and buoyancy forces are scarce. Therefore, our future investigation seems to be applicable for thermo-fluid engineering problems.

5. Conclusions

The numerical simulations have captured the expected physical behavior and identified in the experiments from Liu and Srnsky [9] Re = 7650 without heat transfer and crosswind. The main investigations were carried out through the centroid trajectories of a pair of vortical structures (or primary vortical structures) from wing tips and of instantaneous vorticity and velocity contours. The results obtained for the combined conditions of mixed convection heat transfer and different crosswind velocities have allowed us to conclude that the natural convection (buoyancy forces) interferes with the behavior of the primary vortical structures, both left and right, in different ways: the increase in the Richardson number culminates in increments of altitude in the rebound of the right vortical structure, which rotates downwind, and culminates in the restriction of the rebound of the left vortical structure, which rotates against the crosswind. In other words, an increase in the Richardson number tends to keep the left vortical structure turning closer to the airport runway. On the other hand, an increase in the Richardson number always acts, in both primary vortical structures, limiting their lateral propagation.

The Richardson number expresses the ratio of the buoyancy forces to the viscous forces (shear flow term). If the Richardson number is much less than unity, buoyancy forces are unimportant in the developed flow. On contrary, if the Richardson number is much greater than unity, buoyancy forces are dominant, and there is no necessary kinetic energy homogenizing the fluid. The present investigation has demonstrated that Richardson number defined in terms of aircraft wingspan and constant ground plane temperature establishes a flow that is likely buoyancy driven.

The present numerical method is effective in capturing the primary vortical structure dynamics and, satisfactorily, in predicting the main global quantities of this kind of unsteady flow. The evolution of the primary vortical structures is consistently simulated in accordance with the expected physics. It can be noticed that the thermal behavior in laminar mixed convection heat transfer of a fluid near a heated ground plane has been calculated with acceptable accuracy for the purposes of the present methodology starting physics reported by experimental results without the effects of buoyancy forces and crosswind using the Reynolds number given by Re = 7650.

In a future work, a large eddy simulation (LES) modeling will be implemented to capture surface roughness effects [43,44,45]. It is also planned to study the interference of higher relative roughness values (associated with the surfaces S₁ and S₃, see Figure 1) on the predicted trajectories of aircraft wake vortices near a heated ground plane. Unfortunately, there is a lack of results for airplane wakes (under buoyancy effects) nearest a heated ground plane, which has motivated the present study to contribute in the literature.