Hybrid grid-free and grid-based method for simulation of turbulent flows Open Access

Computational vortex method (CVM) for simulations of free vortex fields was proposed first in 1932 by Rosenhead in a two-dimensional formulation. Later on, the CVM has been extensively used for simulations of vortex wakes behind aerodynamic configurations using potential flow models. The most prominent method among these is the nonlinear vortex lattice method, which is a very efficient tool to calculate lifting surfaces in a wide range of angles of attack. In the eighties, the CVM was extended to viscous flows using Navier–Stokes equation¹ and boundary layer theory.² An excellent review of vortex method is given in the book by Cottet and Koumoutsakos.³ The most modern state of the art in CVM or vortex particle method (VPM) is presented in a recent review by Mimeau and Mortazavi.⁴ 

The vortex particle intensified large eddy simulation (VπLES) method was proposed in our earlier works.^5–8 In contrast to the classic VPM, in VπLES, only the fine scales of vorticity field are represented by particles, whereas the large scales are represented on the grid. Due to this combination, the advantages of grid-based methods and VPM are strengthened, whereas their disadvantages are diminished. So, VπLES is a kind of hybrid grid-free and grid-based method. Hybrid methods based on the decomposition of the computational domain were proposed first by Cottet et al.⁹ and then, more recently, by Pasolari et al.¹⁰ Modern versions of the VPM (see, for instance, Ref. 11) also belong to hybrid methods. They utilize high order finite difference methods (FDMs), written on simple Cartesian grids, for all steps excluding the pure convection of vortex particles, which is treated with VPM. Unlike all hybrid methods, VπLES uses decomposition according to scales at which only the subgrid motion is modeled by the vortex method.

By reducing the contribution of the vortex method to a global solution, we have significantly reduced the problem of numerical instability of the VPM. The experience shows that it is quite a challenging procedure to get a stable solution for 3D problems, especially at high Re numbers using the vortex method in a pure Lagrangian formulation. The vortex amplification results in a strong increase in the vortex particle strength and vortex tube stretching. In inviscid fluid, the vorticity can tend to infinity. Chorin¹² proposed the folding mechanism explaining the conservation of the kinetic energy in ideal fluid under the condition of strong stretching when the vorticity tends to infinity. In the viscous flows, the folded vortex tubes experience mutual viscous diffusion and annihilate each other like in the vortex reconnection process. CVM, being just one of the numerical methods for solution of viscous flow equations, can, in principle, model such processes. However, this needs resolutions comparable to direct numerical simulation (DNS) and handling billions of vortex particles, which is a big computational problem typical of all particle methods. Another difficulty is the need for a high-precision viscosity modeling procedure, which is problematic within the framework of CVM/VPM. Therefore, the common way is to cut fine scales replacing their effect by a suitable subgrid stress (SGS) model. Successful SGS models for the VPM were developed by Winckelmans.¹³ A comparative study of different SGS models in the VPM is presented in Ref. 11.

VπLES as a kind of LES is based on a rather artificial separation of vortex structures into small-scale or subgrid and large-scale ones. Visualization of vortex structures from DNS indicates a rather complex structure of vortices, topologically similar to the structure of a tree in which, at the same time, there is an interweaving of vortices of different sizes (see, for instance, Ref. 14). Their subdivision into large and small structures is an artificial, purely mathematical procedure. This approach would be perfect if there were vortices with sizes of the order of, say, several cell sizes and vortices with sizes much smaller than the size of the grid, and there would be no vortices of other sizes intermediate between these two scale ranges. Then the scale decomposition would be absolutely natural and mathematically consistent. The size spectrum of vortices in a turbulent flow is, however, continuous. It is not entirely clear to which group the vortices with dimensions of the order of the grid size should be attributed, especially considering that their sizes are dynamically changing. It is assumed that vortices are small, and the error of such an artificial classification does not significantly reduce the accuracy of modeling large-scale flows. It is important that the subgrid models, modeling the collective action of small vortices in the statistical sense, correctly describe the main function of such vortices—to dissipate the energy of large-scale flow. In the VπLES method, we followed exactly this concept.

VπLES refers to models in which an attempt is made to restore subgrid fluctuations using a solution obtained for large resolved scales. In other words, VπLES is a closing model for fluctuations based on modeling of subgrid motion by vortex particles. The method is a kind of reduced-order model because the subgrid vortices are replaced by a finite number of discrete vortices that induce velocities determined from the Biot–Savart law. There are four categories of alternative approaches¹⁵ allowing one to restore the fluctuations with scales near the grid cutoff: stochastic, approximate–deconvolution (AD), blending combinations of AD and stochastic methods and kinematic-simulation models. Relevant literature references and short explanation for each category can be found in the article.¹⁵ Another promising approach of fluctuation resolution, which has been intensively developed recently, is the reconstruction of turbulent fields using machine learning methods.^16–19 This reconstruction is usually used for particle-laden flows when the particle motion is strongly influenced by small-scale vortices. The reconstruction is performed on the basis of data obtained for large-scale flows using the training of neural networks. The disadvantage of this approach is that it requires training data obtained from high-resolution DNS calculations. The universality of the training data with respect to the Re number and geometry of the flow remains an open question. One of the most innovative methods proposed recently to represent all scales is the pseudo-direct numerical simulation (P-DNS) (see, for instance, Ref. 20). The fine-scale motion is solved at the stage of preliminary calculation in simplified representative volume elements (RVEs) with various boundary conditions (BCs) written in a dimensionless form. A key assumption of the P-DNS technology is that various boundary conditions can be described by a finite number of dimensionless velocity gradients called the instability index tensor Id_ij. Depending on Id_ij, the fine-scale solutions are stored in a database, which according to the authors is problem-independent. The stored fine-scale flow parameters are utilized in coarse simulations. To do this, the index tensor Id_ij is calculated from the coarse solution and used as input to find the corresponding fine solution stored in the database. The various applications presented in Ref. 20 demonstrate the high performance of P-DNS.

The remainder of this paper is organized as follows: the governing equations and design of the VπLES method are described in Sec. II. Results of VπLES application to homogeneous decaying isotropic turbulence (HDIT), free jet case with passive scalar mixing, and channel flow are given in Sec. III. In addition to the results published in our previous papers,^5–7 we presented some new results, including an explanation of the cutoff of the vortex area influence using the Biot–Savart law, the anisotropy of small-scale motion, various passive mixing results and the results of applying VπLES to wall bounded flows (channel). The conclusions of this research are presented in Sec. IV.

Analyzing alternative approaches to determine the particle strengths, we considered methods of turbulence synthesis using the vortex particle method in the style of work.²⁴ In our works,^25–28 we managed to find vortices that provide the specified Reynolds stress, integral lengths, and spectra. Due to the fact that the velocity field is induced by vortices, it is automatically divergence free. An obstacle to using this technology in the VπLES method was the fact that it assumes the hypothesis of the constancy of Reynolds stresses over time. Strictly speaking, the definition of statistical parameters is based on the hypothesis of stationarity in time averaging or in ensemble averaging and the introduction of lengths presumes spatial homogeneity. This is quite acceptable in CFD to generate an inlet velocity where turbulence is assumed to be stationary. In the general case of an intensive change in flow parameters, in non-equilibrium flows, the use of the stationarity hypothesis in the formalism is not justified. Implementation of such a technique would face the following additional serious problems:

• It requires the additional simulation, which provides time averaged Re stress and Kolmogorov scales (see, for instance, Ref. 24).

• Not the entire Re stress could be ascribed to small-scale motion. There is also the problem of separating that part of the Reynolds stress that corresponds to small-scale motion. For these purposes, it will be inevitable to use models of subgrid motion with all their disadvantages.

Therefore, we determine vortices in VπLES dynamically in time and space as described earlier.

An important problem is that the generation of small-scale vortices is carried out on the basis of a large-scale solution in a region that, when represented in Fourier space, corresponds to high wave numbers near the grid cutoff. In this region, the large-scale flow energy is underestimated due to implicit numerical (or explicit in common LES) filtering. To overcome this disadvantage, the enrichment methods^29,30 can be applied. A simpler and more pragmatic way is to generate small-scale vortices in the region of smaller wave numbers. The newly identified vortex can be converted to the vortex particle when its size approaches to two or three or even more cell sizes. Both ways are utilized in VπLES. Due to the reduced viscosity of vortex methods, the vortex is given the opportunity to increase its intensity and reduce its size as a result of vortex filament stretching.

The algorithm for the particle strength update $α ( t + Δ t )$ is taken from Fukuda and Kamemoto,³¹ with a modification consisting in exclusion of the redistribution step, which results in an avalanche-like increase in the vortex particles' number in areas of strong stretching. Fukuda and Kamemoto³¹ prevented this growth by introduction of the threshold for the dissipation rate, which is difficult to introduce for the general case. This can be seen as a kind of model that removes the finest vortices and drains energy at high frequencies.

If particles grow due to viscosity and flow stagnation and exceed some size $γ Δ$⁠, they are mapped back to the grid. The value γ = 3 was used in the present simulations. The velocities induced by big vortex particles with $σ > 3 Δ$ 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.

If the Reynolds number decreases, and the flow becomes laminar, the velocity field becomes smooth, which leads again to the conditions (19) and (20). Therefore, the VπLES model tends to the laminar flow model. This result is clearly demonstrated in Fig. 17 of Ref. 7 for the free jet with Re number of 100.

Ideally, the scalar should be defined as the sum of $φ g$ and $φ v$⁠. We suppose that $φ ≈ φ g$ due to the following considerations. VπLES on a coarse grid should ensure the accuracy of LES or DNS results (including reproducing small-scale flow details), obtained at higher resolutions, with less computational resources. If this is not the case, the development of alternative approaches does not make much sense. VπLES uses the particle method, which is very laborious and poorly parallelized because it utilizes the integral law of Biot–Savart. For the dynamics of the velocity field, as already shown in Refs. 6 and 7 and will be shown herein, it was possible to demonstrate the advantages of VπLES. For mixing, the introduction of additional scalar particles led to a significant increase in computational time and memory. Therefore, the idea is to improve the accuracy of modeling large-scale mixing by using small-scale vortices, i.e., by the second term on the r.h.s. in Eq. (22).

The idea of enrichment was used in calculating the evolution of the scalar in HDIT. Indeed, in HDIT, the VπLES model leads to a significant decrease in the kinetic energy of a large-scale flow, which is the correct effect of the VπLES model (see Fig. 5), and therefore to a significant decrease in $u v$ pulsation and the strength of the formed vortex particles $α$⁠. As a result, vortex particles are proved to be weak and not able to enhance mixing, which turns out to be even lower than in simulations without turbulence models. When the resolution increases, the correct behavior of mixing is restored. To get correct results at coarse resolution, the de-filtering procedure (25) is applied to small-scale velocity $u v$ (instead of $u$⁠, we use $u v$⁠) in Sec. III A. In Refs. 15 and 29, the dynamic procedure is applied to find B. In this article, B is selected from the match of the scalar simulation results on a coarse grid and the DNS results. In the case of a free jet III B, enrichment was not required due to the intense generation of turbulence and small vortices that already occurred near the nozzle.

The VπLES method was implemented in an in-house code mentioned earlier as well as into OpenFOAM in a MPI parallel mode. The results presented in Fig. 2 demonstrate a very good scalability for number of processors till 128 and cell numbers 64³ and 128³. The computational time is reduced as $∼ 1 / N$⁠, where N is the core numbers.

The first verification study of VπLES was conducted for several two-dimensional cases in Ref. 5. Validations for three-dimensional cases are presented in Refs. 6 and 7. A thorough verification and validation was performed in Ref. 6 for the case of the homogeneous decaying isotropic turbulence in a box with periodic boundary conditions.

The computational domain in the form of a cube with the size $L = 0.508 m$ was meshed with uniform Cartesian grids with 32³,64³, 128³, and 256³ cells. The periodic boundary condition was enforced on each face of the cubical box. The Reynolds number based on the Taylor microscale, air kinematic viscosity $ν = 1.5 × 10 − 5 m 2 / s$⁠, and the root mean square of the longitudinal velocity $u rms = 0.222 m / s$ is $R e λ = u rms λ / ν = 71.6$⁠. The results are evaluated at dimensionless times of $τ = t U 0 M$⁠, where $M = 0.0508 m$ is the size of the grid, which was used for generation of turbulence in the experiments,³³  $U 0 = 10 m / s$ is the air stream speed, and t is the elapsed time in traveling at the mean flow velocity from the grid.

For time discretization, a second-order Crank–Nicolson scheme was used. The spatial discretization of the convective term was performed using the second-order central difference scheme (CDS). The CFL (Courant–Friedrichs–Lewy) criterion was kept less than 0.1. Simulations were carried out using the OpenFOAM code.³⁴ The simulation of Eqs. (3) and (4) is performed with the same time step although the large and fine scales should have different time scales. Since the interaction between vortex particles, which is expected to be fast, plays a negligible role at least for the test case under consideration, the single time step is utilized both for large-scale simulation and VPM. The initial stochastic velocity field was artificially generated according to the procedure proposed in Ref. 35.

As seen in Fig. 3 at coarse resolutions, the implemented VPM acts as a diffusive large eddy simulation subgrid model resulting in a LES-like behavior of the whole method. When the resolution grows at 128³ cells, the effect of the VPM subgrid disappears, and VπLES tends to DNS. A recent analysis has shown that the greatest contribution to energy drain, which is the goal of all turbulent models, is the extraction of the velocity of small vortices from the velocity based on the grid (see Sec. II B).

It was shown in Ref. 6 that VπLES is capable of resolving scales much less than the grid cell size. In such a way, very fine scales can be resolved on coarse grids. Probability density function of the vortex sizes σ is presented in Fig. 4 at different time steps. As seen, there is a wide range of fine vortices with sizes ranging from $0.1 Δ$ to Δ. This is an illustration of the fact that the present method is capable of resolving scales an order less than the grid cell size. Resolution of the vortices with $σ ∼ 0.1 Δ$ in pure grid methods would require grids with number of cells at least 10⁴ larger than in VπLES provided each vortex would be represented only by three cells. The plausibility of simulations at small scales were successfully checked for the second-order statistics, i.e., for the spectral density E and the longitudinal autocorrelation coefficient (see Figs. 12 and 13 in Ref. 6).

The performance of new methods often depends on the solver in which they are implemented. To prove the reliability of the method for other solvers and numerical schemes, VπLES has also been implemented in a completely different, simple in-house code, which is a semi-implicit finite difference code (finite difference method or FDM solver) based on a fractional step algorithm with splitting according to physical processes. The convection and pressure terms are handled explicitly, whereas the diffusion step is implicit. A simple Euler method is used for time discretization, while spatial discretization was performed using a second-order central difference scheme. The time step was chosen in such a way that the Courant criterion was less than 0.1. The same initial conditions are used as in OpenFOAM. The initial distribution of the passive scalar was set to that of the u_x velocity [ $φ ( t = 0 ) = u x ( t = 0 )$] with zero mean value $φ ¯ = 0$⁠. To increase the accuracy, the invariant conditions obtained in Sec. II F 3 were enforced.

As can be seen from Fig. 5, left, and Fig. 3, the performance of the in-house code is slightly worse than that of OpenFOAM, obviously due to the use of less advanced numerical schemes. However, the effect of VπLES on large-scale flow is the same in the FDM solver and OpenFOAM. Both turbulent velocity and scalar fields decay in time (see Fig. 5). The pure grid simulation with 64³ cells underpredicts the decay of the velocity and scalar fields. Deployment of VπLES substantially improves the accuracy of the velocity field evolution, as has already been demonstrated in Fig. 3 using OpenFOAM. However, the scalar dynamics prediction remains not satisfactory. The reason might be the underprediction of small-scale vortices' strength. Application of the enrichment (25) at $B ≠ 0$ allows one to substantially improve the result of simulation. The best agreement with the pure grid simulation with 256³ cells, which can be considered as DNS, is achieved at B = 2. The results are presented till τ = 80. At large τ, the decay of the scalar variance $φ 2$ is overestimated in VπLES because the viscous decay of vortex particle strengths is underestimated by the core spreading model. Enrichment leads to increase in the vorticity of the turbulent field, which is illustrated in Fig. 6. Stronger vortices increase mixing and cause faster attenuation of the scalar variance (see Fig. 7).

The second case of free turbulence utilized for the validation of the VπLES method was a free jet⁷ at the Reynolds number $R e = U j D / ν = 10 4$⁠, where U_j is the jet inflow velocity, and D is the jet diameter. The experimental data^36–38 for $R e = 2 × 10 4$ and $R e = 10 4$ were used for comparison. The momentum thickness was set as $D / θ = 180$ in all simulations.

For the convective term discretization, the second-order central difference scheme (CDS) was utilized, whereas the time discretization was based on a mixed scheme using the Euler (10%) and Crank–Nicolson (90%) schemes. The Courant number was kept less than 0.3 in all simulations. The time averaged values are calculated from 5 to 25 s.

Following Ref. 39, the mean velocity at the inlet is approximated by a constant velocity U_j in the jet core region and a laminar Blasius profile near the wall $r / D = 0.5$⁠. All LES simulations were performed with disturbances in the jet core region and inside the laminar boundary layer near the wall at the inlet according to Ref. 39. The grids used for validations are presented in Table I. The pure grid LES simulations with grid III are used further to evaluate the efficiency of the VπLES method. Based on methodical study, four particles per cell, N_pt = 4, were used in VπLES simulations.

A detailed analysis of results is presented in our paper.⁷  Figure 9 illustrates the evolution of axial velocity fluctuation along the jet centerline obtained using different grids and types of simulations.

The fine vortices in VπLES provokes a more earlier jet breakdown and a correct formation of fluctuations starting directly from the nozzle. The turbulence intensification by vortex particles is proved to be sufficient to trigger the turbulence already at small x/D without inlet perturbations. Up to $x / D = 4.5$⁠, the accuracy of VπLES coarse simulations is satisfactory. Especially at the nozzle exit in the near field (⁠ $x < 4 D$⁠), the coarse LES (grid I) with the dynamic Smagorinsky model shows a strong underestimation of the fluctuation and the discrepancy with reference results, while VπLES simulation shows very good performance. The grid II is derived from the grid I by refining the shear layer region. By increasing the mesh resolution in the shear layer, VπLES results for the axial velocity fluctuation are drastically improved, while the results of LES still shows a strong deviation from the experiment specifically in the region near the nozzle exit (⁠ $x < 4 D$⁠) and in the far field (⁠ $x > 7.5 D$⁠).

Figure 9 shows that VπLES with grid II predicts both near and far fields of the jet with the similar accuracy as the pure grid LES with 6.04 M cells. Figure 9 confirms the grid convergence of VπLES simulations already at $1.5 × 10 5$ cell grid, which is clear from the comparison of the turbulence kinetic energy calculated for $1.5 × 10 5$ and $3.0 × 10 5$ grids.

The positive effect of the coupling term $( u v × ω g )$ in Eq. (3) on the accuracy of the predicting velocity fluctuation is illustrated in Fig. 10. This effect is moderate, since this term describes only part of the interaction between small and large vortices in VπLES. The main effect, which consists in the transfer of energy from a large-scale flow to a small one, occurs due to the generation of vortex particles described in Sec. II B.

Since a deterministic prediction of a turbulent flow is practically impossible, the task of every SGS model is to reproduce the subgrid motion only in the statistical sense. The following features of the subgrid motion should be captured by a proper subgrid model: non-equilibrium effects including laminar-turbulent zones, energy backscatter and anisotropy of fine-scale motion. In this subsection, the anisotropy of velocities induced by fine vortices in VπLES is discussed.

The total flow shows a well-pronounced anisotropy with the dominance of the axial fluctuations on the jet centerline (see Fig. 11). Reynolds stresses $R i i v$ of the velocity field $u v$ also shows a clear anisotropy, which is space-dependent. On the centerline at $x / D > 5$⁠, two diagonal stresses are equal to each other, $R 22 v ≈ R 33 v$⁠, and dominate over $R 11 v$ (see Fig. 12). To explain this effect, we consider stochastic distribution of the statistically independent axisymmetric vortex particles on the centerline with strengths aligned with the x-axis. They induce velocities $u x v = 0$ and $u y v ≠ 0 , u z v ≠ 0$⁠. Due to the axisymmetric character of each vortex, the spatially averaged squares of velocities $u y v$ and $u z v$ are equal, i.e., $( u y v ) 2 ¯ = ( u z v ) 2 ¯$⁠. Precession of vortices around their spins causes the appearance of the longitudinal velocities $u x v$⁠, which are much smaller than $u y v$ and $u z v$⁠. Thus, the jet axis area at $x / D > 5$ is populated by vortex particles with axes predominantly oriented along the jet propagation or mean flow direction. At $x / D < 5$ on the centerline and at $r / D = 0.25$⁠, the fine-scale turbulence is nearly isotropic $R 11 v ≈ R 22 v ≈ R 33 v$ in the beginning of the jet development [see Fig. 13(a)], i.e., this area is populated with vortex particles with orientations uniformly distributed around a sphere. Further downstream, the same anisotropy takes place at $r / D = 0.25$ and on the jet axis. At the jet boundary (⁠ $r / D = 0.5$⁠), the fine-scale turbulence becomes anisotropic with a clear dominance of the radial fluctuations $R 22 v > R 11 v ≈ R 33 v$ [see Fig. 13(b)], i.e., the dominating fluctuations are in the direction of the dominating large-scale entrainment motion. This effect is difficult to explain using the concept of statistically independent distribution of axisymmetric vortex particles, which results in the equality of two stresses dominating over the third one. Most probably, such a distribution of diagonal Re stresses $R i i v$ is due to synchronization of vortex particles' positions and their orientations so that velocities induced by any particle are compensated by its neighbors.

This section presents the results of numerical simulation of scalar propagation in a jet. In the jet nozzle, the scalar value is equal to one and then varies in the range between 0 and 1 due to laminar and turbulent mixing. The profiles of the time averaged scalar $φ$ are shown in Fig. 14 for grid I (left) and II (right). The LES solution for grid III is presented as a reference. Due to the fact that there are no initial disturbances in the VπLES method, slow mixing occurs in the initial sections at $x / D = 2$⁠, which, however, accelerates and reaches acceptable values already at $x / D = 4$⁠. The coarse LES solution on grid I, on the contrary, shows excessive mixing at the beginning of the jet due to increased diffusion caused by low grid resolution. At $x / D > 4$⁠, the accuracy of reproducing the mean scalar on grid I becomes acceptable by both VπLES and LES. With an increase in resolution, there is a significant improvement in the accuracy of modeling over the entire length of the jet propagation (see Fig. 14 for grid II).

The fluctuation $φ ′ = ( φ ( x , t ) − φ ¯ ( x ) ¯ ) 2$ reproduction accuracy is, as expected, lower than for the average scalar $φ ¯ ( x )$ (Figs. 15 and 16). This is especially noticeable in the results obtained on a coarse grid I as shown in Fig. 15 at $x / D = 2$⁠. The simulations carried out with the help of VπLES are much closer to the results of modeling on a fine grid III using LES and without any turbulence model. LES on a coarse grid I significantly underestimates the development of mixing due to overestimated turbulent diffusion on a coarse grid. If we leave only numerical and laminar diffusion as in the simulation without turbulence model, then the results depend significantly on the presence of initial disturbances. With initial disturbances, no model simulation gives results close to VπLES; without disturbances, the fluctuation is underestimated like in a coarse LES. As in the previous case, the accuracy is sufficiently improved on finer grid II. Reliable prediction of pulsation evolution at the very beginning of mixing development is of great importance for the correct modeling of chemical reactions, especially fast chemical reactions that occur even at the initial stages of mixing. In this sense, the VπLES method has shown its advantages by accurately simulating the pulsation at the beginning of the mixing development on a rather coarse grid, that is, with less computer resources than LES with grid III. At large x/D, there is an overestimation of the fluctuation, which, as expected, decreases with increasing resolution (compare the left and right Fig. 16). The VπLES overestimates fluctuations, apparently due to insufficient decay of fine-scale vortices described by the core spreading method.

In the previous chapters, the potential of VπLES method has been evaluated for the simulation of free shear flows. The purpose of this chapter is to evaluate the VπLES method capability of treating near-wall flows. The test case is a fully developed turbulent channel flow through a rectangular duct of height $2 δ$⁠. The size of the computational domain was chosen as $6 δ × 2 δ × 3 δ$ in the streamwise, wall normal, and spanwise directions, respectively. Periodic boundaries have been applied in streamwise and spanwise directions. All simulations are performed for Reynolds number of $R e τ = u τ δ / ν = 395$⁠, which is based on the friction velocity $u τ$ and channel semi-height δ. DNS data taken from Ref. 40 is used as a reference. Parameters of three grids, C1, C2, and C3, used for simulations are presented in Table II.

When using the vortex particle method (VPM) near the wall, a number of problems arise, both common to particle methods and specific to VπLES. The first common problem is the non-physical wall crossing. Due to the temporal discretization of the vortex dynamics, vortices can cross the wall and penetrate into space outside the flow. This problem is usually solved by simply returning the vortex either to its original location or it is mirrored relative to the wall into the flow. Second, the representation of the structures continuously distributed in space by discrete particles with a simple spatial axisymmetric distribution of vorticity, which is usually used in VPM, leads to a serious approximation error and can cause numerical instability (numerical noise). As a rule, the task of synthesizing the near wall turbulence using a set of stochastic axisymmetric vortices faces theoretical difficulties when the anisotropy is strong at low y⁺. For instance, such a synthesis does not work when $R 11 ≫ R 22 + R 33$ in the method⁴¹ and when $1 − [ ( R 11 − R 22 ) 2 + 4 R 12 2 ] 1 / 2 > 0$ in the technique.²⁴ However, for small-scale vortices that are used in VπLES, it is assumed that this anisotropy will be weaker, and the use of axisymmetric vortices is practiced further.

A problem that is specific to VπLES is that the cells are usually very elongated in the streamwise direction $Δ x > Δ z ≫ Δ y$ (see Table II). This causes two difficulties. First, an axisymmetric vortex like a blob is geometrically not a proper model for the vorticity distributed within pancake-like cell. Second, the radius of the vortex particle determined from equality of volumes occupied by cell and particle $Δ x Δ z Δ y = 4 π σ 3 / 3$ can be larger than the distance to the wall $σ ≫ Δ y$⁠. As result, a part of new generated particle could be located beneath the wall, as shown in Fig. 17. If the vortex is located to the wall closer than $2 Δ y$⁠, the particle approximation of vorticity is fully wrong. Formally, the no slip boundary condition (BC) as well as the pressure BC should be modified to take into account the normal and tangential velocities induced by vortex particles at the wall.⁸ 

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 do not allow the particle generation at the distances $y ≤ 4 σ$ introducing the security zone. The width of the zone depends on the grid and covers a few grid wall layers in the 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.

Figure 18 demonstrates results for mean streamwise velocity (left) and turbulent kinetic energy depending on y⁺. Results for LES with dynamic Smagorinsky model and VπLES for a relatively coarse resolution (60 K cells) are very similar and show good agreement with DNS for $⟨ u x ⟩$⁠. For TKE k, both simulations show a significant excess in the range of $10 < y + < 40$⁠. Using $⟨ u x ⟩$ and k as method performance indicators, we have not documented advantages of VπLES with respect to LES with the same resolution. More positive conclusions can be drawn when analyzing the spectrum.

VπLES allows one to reproduce vortices smaller than the cell size. Therefore, it is expected to reproduce high frequencies even with the almost same temporal resolution $Δ t$ as in LES. This effect is illustrated in Fig. 19. The results are presented for the one-dimensional spectrum E_uu calculated using the streamwise velocity $u x v$⁠. Qualitatively, it looks similarly for the transversal velocity components E_vv and E_ww. For the grid C1, the LES with the dynamic Smagorinsky model (dash line) shows a sharp drop in energy content already at relatively small frequencies, whereas the VπLES spectrum (red line) is extended to substantially higher frequencies. The question arises whether the content of high-frequency energy in VπLES is physical or not. To prove it, the LES calculations on the fine grid C3 (solid black line) were performed, which show a relatively good agreement in the area of high frequencies with VπLES. At a certain frequency, the VπLES curve detaches and starts the spectrum part $k − γ$ with $γ ∼ 2$⁠. According to our experience, this behavior is typical for statistics of uncorrelated smoothed vortex particles.

The paper presents a novel hybrid grid-based and grid-free method, VπLES, which is based on the decomposition of scales into large and fine scales. Large scales are simulated on the grid, whereas the fine-scale motion is directly resolved by the vortex particle method. There is a strong interaction between both scales, which consists in (a) generation of fine-scale vortices from large-scale flow simulated on grid, (b) additional influence term $( u v × ω g )$ in Eq. (3), (c) mapping of growing small vortices onto grid solution, and (d) the convection, rotation, and amplification of small vortices caused by large-scale motion. A very time-consuming summation of velocities induced by vortices according to the Biot–Savart law is sufficiently reduced using the cutoff of the influence area, which is possible because of stochastic distribution of fine vortex orientations.

In the application examples described in this article and in our previous works, a good accuracy of VπLES is achieved with grid resolution and calculation time, which are lower than in purely grid simulations using DNS and highly resolved LES. Therefore, the new method can be an alternative to the existing purely grid techniques. The method is validated for the passive mixing using a reduced form of the scalar equation without introduction of scalar particles, the handling of which could be very laborious. This results in a certain underestimation of mixing. However, a good agreement with DNS was obtained for the passive mixing in homogeneous decaying isotropic box turbulence using the enrichment procedure, which enhances the fine scales close to the grid cutoff.

For the free jet, the convincing results were obtained for the turbulent velocity field. The comparable accuracy with highly resolved LES was obtained at substantially lower grid resolution and computational costs. Consideration of the additional coupling term $( u v × ω g )$ sufficiently improves the reproduction of the axial velocity fluctuations. It is proved that the velocity field induced by small-scale vortices is sufficiently anisotropic. For the passive mixing in the jet, the VπLES method has shown its advantages by accurately simulating the scalar pulsation at the beginning of the mixing development on a rather coarse grid, that is, with less computer resources than LES.

The application of VπLES to the flow in the channel shows a relatively weak effect of vortex particles on mean velocities and velocity fluctuations. First, this is due to elongated cells, which are usually used near the wall and which are not suitable for implementing VπLES utilizing axisymmetric vortex particles. Second, it is obvious that flows confined by walls require a minimum high resolution near the wall, and the results of coarse resolution cannot be sufficiently improved using any subgrid model. However, it was possible to better capture small scales and improve the reproduction of the spectrum in the high-frequency region.

The support of the German Research Foundation (Deutsche Forschungsgemeinschaft) under Grant Nos. AB 112/10-1 and INST 264/113-1 FUGG is gratefully acknowledged. Special thanks to Hesam Tofighian for parallelizing VπLES in OpenFOAM and Dr. Jordan Denev (both from Karlsruhe Institute of Technology) for numerous discussions.

The authors have no conflicts to disclose.

N. Kornev: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). S. Samarbakhsh: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal). J. Darji: Data curation (equal); Investigation (equal); Validation (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.