Theoretical Background of the Hybrid V π LES Method for Flows with Variable Transport Properties

: The paper presents the theoretical basis for the extension of the V π LES method, originally developed in recent works of the authors for incompressible ﬂows, to ﬂows with variable density and transport properties but without chemical reactions. The method is based on the combination of grid based and grid free computational particle techniques. Large scale motions are modelled on the grid whereas the ﬁne scale ones are modelled by particles. The particles represent the ﬁne scale vorticity, and scalar quantities like e.g., temperature, mass fractions of species, density and mixture fraction. Coupled system of equations is derived for large and ﬁne scales transport.


Introduction
In our previous works [1][2][3] we propose a novel numerical technique, referred to as the VπLES approach, which utilizes the combination of grid based and grid free computational methods. The large scale motions are modeled on the grid whereas the fine scales are represented by vortex particles. Advantages of the new technique were clearly demonstrated in validation tests for canonical wall free flows. At present, the validation is continued for wall bounded flows, which results will soon be published. The idea of the method is relatively simple. Due to inherent instability of turbulent flows, the grid based solution shows the tendency to create flow structures with sizes comparable with the cell size. In grid based methods, these structures either disappear due to artificial viscosity in schemes with a low order accuracy or result in a high frequency instability which is well seen in the spectra as a pile up of the kinetic energy at high frequencies [2] when high order schemes are utilized. Within the framework of VπLES such structures are recognized and converted to vortex particles which don't depend on the grid. A fundamental assumption of the VπLES approach is the decomposition of any physical quantity, say the velocity u, into the grid based (large scale) u g and the fine scale u v parts, whereas u g is resolved on the grid and u v is represented by particles. Dynamics of large and fine scales is calculated from two coupled transport equations one of which is solved on the grid whereas the second one utilizes the Lagrangian grid free Vortex Particle Method (VPM). We classified the method as the Large Eddy Simulation (LES) with a direct resolution of the subgrid (SG) motion modeled by u v part. Numerical simulations [2,3] showed that vortex particles play a role as triggers or intensifiers of the turbulence in the grid based solution. Their inner interaction does not contribute sufficiently to the flow evolution. Hence the name of the method is the LES intensified by the vortex particles or VπLES .
The aim of the present paper is to develop a theoretical framework for extension of the VπLES method towards flows with variable density and transport properties. Variation of density and transport properties is assumed to be due to change of the temperature and chemical composition of the gas mixture. In this paper, such flows will further be referred to as the compressible flows at low Mach number or just compressible flows. The flow with a constant density is called incompressible in this paper. As the most general case with chemical reactions (e.g., combustion) would rise further questions and increase the volume of this work excessively, the present paper focuses on flows for which mixing and temperature gradients are important, but chemical reactions do not occur. Examples are the distribution of hazardous emissions from chimneys or vehicles in urban environments, flows in air-conditioned rooms and in large atria, exhaust gas recirculation or heat utilization.
Within the framework of VπLES (both compressible and incompressible) only the fine scales of vorticity field are represented by particles. The separation between the large and fine scales is done using the spatial LES filtration. Similarly to the incompressible case the fine scale velocity field u is calculated from the grid one u g using the LES filtration operation denoted by overlines: According to the Helmholtz decomposition every velocity field including the velocity pulsation u can be represented as the sum of three components: The first term is the velocity induced by vorticity, the second one is due to the flow expansion and the third one is the potential velocity arising due to influence of flow boundaries. We suppose that the fine vortices are generated at a some distance far from the boundaries and don't influence the boundary conditions on boundaries. In this case the boundary effect is entirely taken into account by the large scale velocity u g , i.e., u p (x, t) = 0. As follows from (2) the fine scale velocity should be represented not only by vorticity, but also by expansion sources. Therefore, in contrast to the incompressible case we have to introduce two sorts of particles: vortex particles and expansion particles. In principle they are different, because vortex particles arise in areas of high vorticity whereas the expansion particles should correspond to the areas of high expansion rate characterized by the material derivatives of the density. Although in the paper [4] the same set of particles is used for vortex particles and expansion particles, in this paper we will consider a more general case introducing different particle sets.
The algorithm for determination of the strengths of the vortex particles inducing u ω (x, t) (see [2,3]) is based on the identification of vortex structures in the field u (x, t) using the λ ci criterium [5] which is directly extendable to compressible flows [6]. At locations of large λ ci values we insert new vortex particles with the strength being equal to the local vorticity ω = ∇ × (u (x, t) − u v old (x, t)) multiplied with the cell volume, where u v old is the velocity induced by vortex particles already present in the flow. The vortex particles induce the velocity u v which, in the ideal case, should be equal to u − u v old . Because of discrete representation of ω by discrete vortex paticles, in the general case we have: old and a special rescaling procedure was introduced to match both velocities ( [2,3]). From the physical point of view u ε (x, t) arises due to the temporal change of the density fluctuation. Modeling of the expansion velocity is considered in details in Section 5.
Mathematical models of the flow with variable density and viscosity implies the solution of additional transport equations for scalar quantities like like, e.g., temperature and mass fraction of chemical species. Solution of the scalar transport equations within the VπLES framework is discussed in Section 8.

VπLES for Incompressible Flows
The purpose of this section is to shortly present the algorithm of the VπLES method for incompressible flows to make it easier to understand for compressible ones. The large scale vortices are reproduced on the grid. At each time step, the fine vortices consist of new generating vortices (filled circles) and "old" ones generated in the past (symbols * ) (see Figure 1). The main equation governing the large scale motion reads (see [2]): which is solved using the OpenFOAM. The motion of fine vortices is governed by the Equations (7) and (8). Once at any time instant, the grid solution step is completed, the evolution of fine vortices is carried out in the following sequence: generation, motion, feedback effect on large scales through the term u v × ω g , mapping to the grid or elimination. Below these processes are considered in more details.

The Algorithm of Creation of New Vortices
The new generated vortices are responsible for the fine scale velocity field where u v old is the velocity induced by vortices generated in the previous time steps. The vortex identification criterion λ ci [5] is calculated at each cell using the velocity field (4). The cells at which λ ci > λ ci,min contain the vortex candidates which in principle can be converted to vortex particles. They will be marked as active ones using the λ i,active field: where λ ci,min is a certain small value (a free parameter of the algorithm) introduced to sort out the weak vortices. To keep the required computational resources on an acceptable level, only small vortices with size proportional to the local cell size ∆ are to be converted to single vortex particles. It is supposed that the larger vortices with scales of a few ∆ can accurately be represented on the grid. Therefore, they should be identified and separated from small scale ones. For that, all neighboring cells of the i-th cell are checked for the condition λ i,active = 1. If all neighbors fulfill this condition we identify a cell cluster which remains on the grid and all their cells become non-active λ i,active = 0. Such algorithm recognizes A-cells in Figure 2 as those belonging to clusters whereas the cells B are cells with new generating vortices. Only activated vortices in cells with λ i,active = 1 are to be replaced by vortex particles. At each cell with λ i,active = 1 a single new vortex particle is introduced at the cell center. Thus, each cell contains old vortices and the new one ( Figure 1). The radius of the new vortex is set as σ = βVol 1/3 i , where β is the overlapping ratio which is taken as β = 2 and Vol i is the volume of the i-th cell. The vortex particle strength is calculated as α = Vol i ∇ × u . If the distance from the new particle to one already existing in the cell is smaller than some permissible distance ϕ∆, ϕ 1 they are merged (the cell B in Figure 1). The velocity u v new (x, t), induced by the new generated vortex particles, is calculated at grid points using the Biot-Savart law. The vorticity ∇ × u v new is not equal to the target vorticity ∇ × u because of discrete particle representation of the continuous vorticity field inducing u . Therefore the rescaling procedure should be introduced in the next step of the algorithm: The rescaled velocity field u v rescaled (x, t), induced by the new generated vortex particles, is calculated at grid points using the Bio-Savart law and subtracted from the grid velocity u g,new = u g − u v rescaled . Thus, the total velocity at grid points u g,new + u v rescaled = u g remains constant after the vortex particle generation procedure. To keep the number of vortices at any cell restricted, in case of N > N pt the weakest vortices populating this cell are deleted keeping the total number of vortices per cell less or equal to N pt . At N pt = 4 these are cells A in Figure 1. If the number of vortices is less than N pt no vortices are deleted.

Update Procedure for the Particle Positions and Strengths
The fine vortices (old and new generated ones) are moving and change their strengths. In our previous works [2,3] we showed that the vorticity transport equation can sufficiently be reduced by neglecting inner interaction between vortex particles. With the other words, the motion of vortices and the strength change can be calculated without consideration of the velocity induced by fine vortices. Therefore, in the next step the velocities u g and tensor ∂u g i /∂x j are calculated at vortex centers from the OpenFOAM (OF) simulations. The routines for the interpolation of u g , ∂u g i /∂x j and ω g = ∇ × u g from grid points to particles, particle tracking, filtering procedure to calculate u v × ω g were available in OpenFOAM and effectively used to develop a fast working code.
Vortex displacement ξ is found from the equation The vortex strength is changed according to the algorithm described in Section 6.2, in which the density should be constant, i.e., dρ dt = 0 and the temporal change of the fine vorticity is calculated as If during the motion the fine vortices become very small σ → 0 due to stretching, they are eliminated.

Calculation of the Influence of Fine Scales on the Grid Based Solution
The term u v × ω g in (3), which takes the influence of fine scales on the grid based solution, is treated in an explicit way. The velocities u v = u v rescaled are calculated at cell centers using the Biot-Savart law. The grid based vorticity ω g is taken from the OF numerical solution. The filtering is computed using the routines build in the OF code.

Mapping Big Vortex Particles Back onto the Grid
There is a permanent exchange between the small vortices and large scale ones represented on the grid. Large scale vortices become small due to stretching and are converted to particles according to the algorithm in Section 2.1. If particles grow due to viscosity and flow stagnation and exceed some size γ∆ they are mapped back to the grid. For that the velocities induced by this vortex particle are calculated at the grid points within adjacent cells and added to the grid based velocities u g . After that the vortex is deleted from the list of fine vortex particles. γ∆ can be taken equal to four according to the experience from [3].

Momentum and Vorticity Equations for Compressible Flows
Neglecting the volume viscosity, the Navier Stokes equation for compressible flows (see [7], pages 47-53) reads: where ρ, u, p, g, µ are, respectively, density, velocity vector, pressure, body force vector, dynamic viscosity and ε = ∇ · u. Applying curl operator to (9) we get consequently where The continuity equation imposes the mass conservation principle, which in case of variable density is: ∂ρ ∂t

Decomposition of Vorticity and Continuity Equations According to Scales
Following [1], the governing equations of the hybrid grid-free and grid-based method are derived using the decomposition of the velocity, vorticity, density and pressure fields into large scale, represented on the grid, and fine scale, represented by particles, components, i.e., Fine scale variations of both kinematic and dynamic viscosity coefficients as well as of diffusion coefficients are neglected as it is common in numerical methods for turbulent flows with mixing, buoyancy or combustion. The same assumption applies also for the pressure fluctuations [8]. Substitution of the decomposition (12) in the (10) gives The last equation is split in two equations for large and fines scales which are treated in a coupled manner. It should be noted that this splitting is ambiguous (see comments in Section 2 in [2]) and the best choice of it is the matter of experience. In this paper we suggest the following splitting: The sum of (14) and (15) retrieves the original Equation (13). The first equation is then rewritten in u − p variables to utilize any available CFD code based on a grid method. From (14) we obtain where P is a scalar function referred to as the grid based pressure and The Equation (16) is solved on the grid whereas the Equation (15) is treated using the VPM. The term (ρ g + ρ v )(u v × ω g ) takes into account the influence of fine scale motions on large ones. It has the same form as that derived in [1] for incompressible flows. Terms with the upper index "g" in the Equation (15) represent the influence of large scales on small ones.
The continuity equations can be split in a similar way: The Equations (16) and (17) involve terms of different spatial and time scales. Some terms can become invisible on the grid and need a special treatment. This can lead to the stiffness of these equations which is a common problem in modeling of multiscale physics. An usual way to overcome this problem is the spatial averaging of terms representing fine scales (see [9]) as done in our previous works for incompressible flows [1][2][3]. The spatially filtered equation takes the form: ∂ρ g ∂t Solution of the Equation (20) is performed using any grid based method and is further not discussed. Various LES filters build in a grid based code can be used to calculate the spatially averaged terms.

Solution of the Continuity Equation
Solution of the Equation (17) can be done in explicit, semi-implicit and fully implicit ways.
The presence of the source term − ∂ρ g u v ∂x i on the r.h.s of (17) doesn't cause principal difficulties. This is true for both explicit time-stepping methods and pressure-correction methods. We first note here, that this term is known at the time when the continuity equation is to be updated on the grid. Both explicit time-stepping and pressure-relaxation methods make use of source-terms, so that the r.h.s. of Equation (20) becomes simply a part of these source terms. Therefore, the computation of Equation (20) does not require any change in conventional algorithms for the solution of the continuity equation.
The remaining question is the handling of the second Equation (18). For that, two strategies are possible. First, one can assume that the whole density is represented only on the grid, i.e., we neglect the subgrid change of the density ρ v = 0 and the fine scale expansion velocity u ε (x, t) = 0. For instance, in the recent paper of Ghoniem's group on combustion [4] only the flow is calculated using the vortex particle method whereas all scalars including the density are computed on the grid. Then only the Equation (17) remains.
Within the second strategy the Equation (18) is taken into account in the computation of the velocity u ε (x, t) which is calculated from the expansion field ε v (see also the formula (8) in [4]): where Introduction of the expansion velocity is a common procedure developed for the application of vortex methods to compressible flows. The expansion velocity u v ε is expressed through ε v using the Helmholtz theorem: Discretization of (22) using the expansion particles is given in [4] (see formula (15)-(17) on the page 971 in [4]).
The last question to be discussed is the determination of the subgrid density ρ v . We consider the case that density is a function of the temperature T and concentration of chemical species Y. Therefore, it is quite reasonable to assume that the particles carrying the fine density fluctuations ρ v coincide with those carrying the fine scale fluctuations of the temperature T v and species Y v which will be introduced in Section 8. In this way we introduce three set of particles. The first one carries the vorticity field, the second one carries the temperature field and the third one describes the transfer of the chemical species. The last two sets carry also the density fluctuations ρ v . In case of single-component gases, only the fluctuations of temperature are considered for the density fluctuations. The movement of all particles is calculated from the trajectory equations where ξ is the particle center. The derivative Dρ v Dt describes the change of ρ v along the particle path:

Relations for Vortex Tubes in Compressible Flows
The algorithm of the solution is based on the Kelvin and Helmholtz theorems which state that for both compressible and incompressible smooth isentropic ideal flows, vortex tubes move with the flow, and their strength remains constant [10] |ω where σ is the radius of the vortex tube and |ω v | is the vorticity in the cross section of the vortex tube. Differentiating (25) on time gives Conservation of the mass written for the piece of the vortex tube of the length l reads Derivative of (27) on time provides the formula for the radius σ change and

Calculation of the Vortex Particle Transport
The formulas derived above are utilized to extend the algorithm [11], which we successfully used in our previous papers [2,3] to calculate the dynamics of vortex particles in incompressible flows, towards compressible flows. In what follows, under the designation ρ is used for the sum ρ g + ρ v . The algorithm consists of the following steps: • Calculation of the particle core radius from the equation describing the transport of an elementary tube with length l and radius σ in inviscid flow • Consideration of the viscosity influence using the core spreading method (CSM) (see [12]). The particle core radius is increased by ∆σ: • Calculation of the new particle orientation using the Euler method The last equation follows from (25) using the simple algebraic transformations •

Calculation of the new strength vector
In the original version proposed by [11] the next step should be the redistribution of particles whose length is doubled. According to our experience the redistribution results in an avalanche-like increase of the vortex particles number in areas of strong stretching. To prevent this [11] proposed a special elimination procedure based on a knowledge of a threshold for the dissipation rate which is difficult to set in a general flow case. To develop a robust code, to obtain a stable solution and to keep the particle number in a reasonable range we avoid the redistribution procedure in our computations. Thus, the smallest vortices are removed. This reduces the range of scales that must be resolved in a numerical calculation. Such a reduction is an immanent part of every turbulence model.

On the Model Reduction
A successful numerical method should provide the desired accuracy and be as simple as possible. At a glance the presented hybrid method may seem not attractive because of its seeming complexity. Especially the equation for the fine vorticity transport Equation (15) seems to be very complicated. However, as shown in [2,3] for the incompressible flow, the vorticity transport equation can sufficiently be reduced by neglecting inner interaction between vortex particles without a significant loss of the simulation accuracy. The following three levels of complexity have been considered for incompressible flows: • Full vorticity transport equation (full formulation): • Vorticity transport equation without stretching and convection caused by fine vortices (passive vorton): • Vorticity transport equation without stretching (designated as the woStrengthUpdating in Figure 3): Comparison between three models is presented in Figure 3 for the r.m.s of axial velocity along the axis of the free jet at the Reynolds number of 10 4 . Detailed description of the jet case and results are presented in [3]. The difference between results obtained using the full model (38), the passive vortices model (39) and the model without stretching term (38) is negligible pointing out that vortex particles serve just as triggers or intensifiers of turbulence and their inner interaction does not contribute sufficiently to the flow evolution. This conclusion raises the quite relevant hypothesis: perhaps the whole theory can radically be simplified, the moving vortex particles can be neglected, the trigger effect can artificially be simulated using any simple and computationally cheap procedure. To prove this hypothesis we left only new generated vortices, which induce perturbations of the grid solution by the term u v × ω g and then are deleted. No transport equation for the fine vorticity like (15) and no particle transport equation like (23) are solved. The particles are generated at each time instant, produce the perturbations, are deleted and new ones are generated again at the next time step, and so on. Such a trigger is not quite artificial and arbitrary because the perturbations caused by new generated vortices depend on the local flow parameters. The result of this proof are presented in Figure 4. The fluctuations obtained using the radically simplified model, designated as the reduced model, are very close to the target values in the initial stage of the jet development and substantially different at the distance of approximately eight nozzle diameters and downstream. The reason of this discrepancy is that the flow at large x/D has a lot of small vortices, generated upstream, which are neglected in the reduced model. These vortices make a great contribution to the dissipation of the turbulent energy, which in case of their absence is strongly overestimated. This analysis leads to the conclusion that the basic equations of fine scales transport can substantially simplified but not excluded at all. Bearing in mind these results, we expect that the fine vorticity transport equation for compressible flows (15) can also be strongly simplified without the loss of accuracy of simulations. The rate of simplification should be stated in serial methodological computations using the present model.

Transport of a Passive Scalar
In this paper we consider the case with the change of the density and viscosity happens due to mixing of different species without chemical reaction. The problems under consideration are formulated in terms of the mass fraction Y k of any of the k species, the temperature T, or, the mixture fraction f which transport is described by the equation: where Γ f is the diffusion coefficient that in general depends on both temperature and species composition. The passive scalar fluctuations f , are directly resolved in the extended VπLES approach, as shown below.
In the same manner as for the velocity field we decompose the scalar f into large and fine scales f = f g + f v using the LES filtering. The decomposed form of (41) reads The Equation (42) is solved on the grid, whereas the second Equation (43) is treated by the scalar particle method with the particles represented in the form where σ f is the scalar particle size, f 0 is the particle strength and ξ is the particle center. Numerical procedure for the calculation of the term u v j ∂ f g /∂x j is the same as that for the calculation of the term u v × ω g in VπLES.
The last Equation (43) is treated in the form: The first diffusion on the r.h.s of (45) can be treated using the core spreading method with f 0 remaining constant. The gradient ∂ f v ∂x j is calculated by the analytical differentiating (44). Trajectories of the scalar particles are calculated from (23).
In cases of no reaction, the temperature transport equation has the same form as the mixture fraction Equation (41) and can similarly be decomposed.
The mixture fraction formalism allows to reduce the number of particle sets to three sets: vortex particles carrying the vorticity, scalar particles carrying the mixture fraction, species concentration and ρ v as well as the scalar particles responsible for the temperature and ρ v .
The method can principally be extended to more complex chemical reactions.

Boundary Conditions
The boundary conditions (BC) are explicitly formulated only for the grid solution. There are no boundary conditions for the vortex method. For instance, the no slip BC for the velocity can be written in the form [1] u The necessary and sufficient boundary condition for the pressure in incompressible case, which is used in the Poisson equation within pressure correction approaches, is the Neumann BC obtained by the projection of the Equation (3) onto the normal direction n [1]: where the term proportional to the kinematic viscosity was traditionally neglected. Similar BC can be derived for the scalar field and for flows with variable transport properties. We don't derive these BC and use the traditional BC for the grid based solution from the following considerations: • Close to the wall the flow is rather smooth and the creation of energetically important concentrated fine structures is unlikely.

•
Representation of the continuously spatially distributed structures by discrete particles with a primitive spatial spherical distribution of vorticity ω = 3/(4πσ 3 )αe −ξ 3 /σ 3 results in a serious approximation error and cause numerical instability.

•
This approximation can be not physical if the vorticity induced by a particle is not zero beneath the wall ω(y < 0) = 0, where y is the distance from the wall (see Figure 5). This case takes place when the elongated cells with ∆x ∆y, where x is the coordinate along the wall, are usually used. In this case the vortex radius σ ≈ β∆x = 2∆x ∆y. If the vortex is located to the wall closer than 2∆x the particle approximation of vorticity is fully wrong.
To avoid these difficulties and taking into account the fact, that the influence radius of the vortex particle (area of a sufficient induced velocity ) is approximately 2σ we don't allow the particle generation at the distances y ≤ 4σ introducing the security zone schematically shown in Figure 1. The width of the zone depends on the grid and covers a few grid wall layers in vertical direction. Since the influence of particles at the boundary in such a way is avoided, we use traditional BC for the grid part of the solution. The same procedure is valid for the scalar particles. Figure 5. Wrong approximation of the vorticity by a vortex particle when the distance from the wall y is less than the vortex particle radius σ.

Implementation of the Method
The coupling of the Particle Method with the existing implicit codes is performed in an explicit way. The terms (ρ g + ρ v )(u v × ω g ), u v j ∂ f g /∂x j and ∂ρ g u v i ∂x i taking into account the influence of fine scales on large ones are taken from the previous time step. Following the particle method, the partial differential Equations (15) and (43) are reduced to the ordinary differential Equations (23), (24), (31), (32), (34), (35), (37) and (45) which are solved using the simple Euler approach, i.e., the r.h.s are taken from the previous time step. Higher order methods like Runge-Kutta or Bashforth-Adams are not necessary, because as mentioned in Section 7, the particles play the trigger role and the accuracy of their inner dynamics calculation is not important.
The full solution takes place in the following order employing the OpenFOAM (OF) code: • General structure of the time loop in OpenFOAM. Steps 2-5 and 6-11 are additional routines implemented into OpenFOAM (OF).
3.8 Calculation of the particle core spreading using Core Spreading Method (33).
3.12 At each cell N − N pt weakest vortices are deleted if N > N pt . • Fine Scalar Part calculates: 8.1 tensor ∂ f g i /∂x j at vortex centers from OF. 8.2 Particle displacement from equation dξ/dt = u v + u g .
8.3 Change of the scalar strength magnitude at particle centers from (45) without the diffusion term.
8.4 Calculation of the particle core spreading using Core Spreading Method (33). 8.5 At each cell N − N pt weakest particles are deleted if N > N pt . Their scalar value is redistributed among neighboring cell centers to ensure the conservation of the scalar.

Conclusions
The paper presents the theoretical background of the hybrid grid-based and grid-free method which is the extension of the VπLES method, developed originally for incompressible flows, to flows with variable density and transport properties. The method is based on the decomposition of flow parameters into large scale and fine scale parts. The governing equations are obtained through the substitution of this decomposition into the transport equations for vorticity, scalar dynamics and continuity equations. This results in the splitting of transport equations according to scales. Fine scale dynamics is simulated with Lagrangian grid free Particle Method. The numerical approach is designed so that the strengths and the weaknesses of grid-based and grid-free schemes can be seen as complementary. Particle methods are suitable for modeling of fine and fast vortices because they possess a low artificial viscosity, keep vortex particles from excessive diffusion and has no restrictions to Courant number.
The derivations are presented on the level of formulas and algorithms which can directly be implemented into the software code. The overall method is designed so that it can be incorporated into the existing grid based powerful codes like OpenFOAM which have extensive modeling capabilities of multi physical processes.
Author Contributions: Two authors N.K. and J.D. contributed equally to derivation and proof of the theoretical basis of the VπLES method extended to flows with variable density and viscosity. S.S. performed simulations presented in Section 7. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: