Vortex particle intensified large eddy simulation - VπLES

Vortex Particle Intensified Large Eddy SimulationN. Kornev and S. Samarbakhsh VI International Conference on Particle-based Methods – Fundamentals and Applications PARTICLES 2019 E. O˜nate, M. Bischoff, D.R.J. Owen, P. Wriggers & T. Zohdi (Eds) VORTEX PARTICLE INTENSIFIED LARGE EDDY SIMULATION - VπLES S. Samarbakhsh1and N. Kornev2 Chair of modeling and simulation University of Rostock Rostock, Germany 1e-mail: [email protected] 2e-mail: [email protected] Key words: LES, subgrid model, meshless method, free jet, hybrid Eulerian/Lagrangian method Abstract. This paper presents a novel Large Eddy Simulation approach with a direct resolution of the subgrid motion of fine concentrated vortices. The method, proposed first by [10], is based on combination of a grid based and the grid free computational vortex particle (VPM) methods. The large scale flow structures are simulated on the grid whereas the concentrated structures are modeled using VPM. Due to this combination the advantages of both methods are strengthened whereas the disadvantages are diminished. The procedure of the separation of small concentrated vortices from the large scale ones is based on LES filtering idea. The flow dynamics is governed by two coupled transport equations taking two-way interaction between large and fine structures into account. The fine structures are mapped back to the grid if their size grows due to diffusion. Algorithmic aspects specific for three dimensional flow simulations are discussed. Validity and advantages of the new approach are illustrated for a well tried benchmark test of the decaying homogeneous isotropic turbulence using the experimental data of [4] and free turbulent jet flow using experiments of [8, 2, 15]. 1 INTRODUCTION Insufficient resolution of fine vortex structures in turbulent flows is one of the key problems in Computational Fluid Dynamics (CFD). The most advanced and popular technique to resolve multi scale flow structures is the Large Eddy Simulation (LES) which is based on the idea of scale decomposition into large and small ones. While the large eddies are directly resolved on the grid, the effect of small vortices is taken into account through a subgrid stress (SGS) model. The subgrid motion is not resolved in LES but rather it is modelled using different functional and structural approaches. However, there are many problems which require direct representation of the subgrid motion to simulate, for instance, mixing or particle 1 174 S. Samarbakhsh and N. Kornev dynamics in turbulent flows. In our previous papers (see [9], [10], [11], [16] ]), we proposed a simulation technique resembling LES with an effort to directly reproduce the subgrid motion at least in the statistical sense. It is suggested to apply a hybrid grid and particle based method, utilizing a combination of the finite volume and computational vortex particle (VPM) (see [5]) methods. The large scale field is represented on the grid like in LES, whereas the small scale one (subgrid field) is calculated using the VPM. The new method called VπLES is a purely Lagrangian one for small structures and purely grid based one for large scale structures. The method is based on the decomposition of the velocity uand vorticity fields ωinto the distributed large scale (upper index ’g’) and concentrated small scale (upper index ’v’) fields: u(x,t)=ug(x,t)+uv(x,t),ω(x,t)=ωg(x,t)+ωv(x,t) (1) The fine vortex detection procedure utilizes the Large Eddy Simulation (LES) filtration applied to the grid based velocity field ug: ug(x,t)= ∞  −∞ ug(s,t)F(x−s)ds(2) where F(x−s) is a certain filter function. The small scale velocity field u�calculated as the difference between the original and filtered fields u�(x,t)=ug(x,t)−ug(x,t) (3) should be approximated by vortex particles in regions of concentrated vortices which are detected using any vortex identification criteria, for instance, λci ([1]). The cells with λi>λ ci,min contain the vortices which in principle can be converted to vortex particles. Such cells are marked as active ones using the λi,active field: λi,active =1,if λci >λ ci,min 0,otherwise (4) where λci,min is a certain small value introduced in order to limit the number of particles. 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. Neighboring cells which all have λi,active = 1 form large vortices. For them it is supposed that the larger vortices with scales of a few Δ can accurately be represented on the grid. Therefore the next task is to detect cells with fine vortices. According to the algorithm, all neighboring cells of the i−th cell are checked for the condition λci >λ ci,min. If all neighbors fulfill this condition, we identify a cell cluster which remains on the grid and all its cells become nonactive λi,active = 0. Only vortices in cells with λi,active = 1 are to be replaced by vortex particles. At each cell with λi,active = 1 the new vortex particle is introduced at the cell center if the permissible number of vortex particles per cell Npt is not exceeded. Otherwise, 2 175 S. Samarbakhsh and N. Kornev the new vortex replaces the cell’s weakest one. The number Npt was introduced to keep the total number at a reasonable level. This restriction is conform with the concept that the largest contribution to the subgrid kinetic energy is made by a small fraction of the strongest vortices. The radius of the new vortex is set as σ=βV ol1/3 i, where βis the overlapping ratio which is taken as β= 2 and V oliis the volume of the i-th cell. A thorough analysis of the influence of Npt and λci,min for the jet case is given in [11] in Sec. 4.3.5. The vortex particle strength is calculated as α=V oliωv=V oli(∇×u) (5) The velocity uv(x,t), induced by the vortex particles, is calculated at grid points xusing the Biot-Savart law uv(x,t)= 1 4π∇×FlowV olume ωv(ξ,t) |x−ξ|dV (ξ) and subtracted from the grid velocity ug,new =ug−uv. Thus, the total velocity at grid points ug,new +uv=ugremains constant after the vortex particle generation procedure. 2 GOVERNING EQUATION The evolution of vortex particles and large scale flow represented on the grid is described by a system of two coupled transport equations derived in [9, 10, 11] for incompressible isothermal flows: ∂ug ∂t +(ug·∇)ug=−1 ρ∇pg+νΔug+uv×ωg(6) dωv dt =(ωv·∇)(uv+ug)+νΔωv+∇×[uv×ωg−uv×ωg],(7) The sum of the curl of the first equation and the second equation retrieves the original Navier Stokes equation written in the form of the vorticity transport equation. The first equation (6) is coupled with the second one (7) through the additional term uv×ωg whereas the coupling of the second equation with the first one is due to the terms (ug· ∇)ωv,(ωv·∇)ugand ∇×[uv×ωg−uv×ωg. The equations (6) and (7) are solved sequentially. The first equation is solved on the grid whereas the second one uses the grid free Vortex Particle Method (VPM) [5]. The physical meaning of the coupling term uv×ωgis explained in [10]. The vortex particle displacement is calculated from the trajectory equation dri dt =ug i+uv i,(8) where iis the particle number. Computation of the velocity induced by vortex particles uvis performed with the direct summation of the Bio-Savart law taking into account one or two layers of neighboring cells. Only induction of the neighboring points is taken into account because the velocity uvis much less than ugand it is mostly determined 3 176 S. Samarbakhsh and N. Kornev by interaction of neighboring particles lying at a short distance. The justification of this simplification presented in [11] is that the correlation between neighboring small scale vortices is weak and they are well separated. This simplification can be considered as a kind of model which results in a very fast computational procedure of a local character suitable for parallel calculations. The velocity induced by a vortex particle can be calculated from the formula uv p=1 4π α×ξ ξ3(1 −e−ξ3/σ3) (9) proposed by [14]. Here αand σare the strength and the radius of a vortex particle, which are defined below. The velocity induced by a set of particles uvis calculated as the sum of uv p. It should be noted that the technique presented here is independent of any specific choice of vortex particles. Particularly, a set of functions introduced in [12] can be used within the present method. 2.1 Numerical solution of the equations (7) and (8) using the VPM Instability of numerical solution of the equation (7) caused by the stretching term (ω·∇)uis the most important problem of the VPM along with the computation of the velocity uv. In grid based methods with low and moderate order schemes, the action of the stretching is effectively counterbalanced by the numerical viscosity which is very low in Langrangian vortex particle methods. Theoretically, a stable VPM solution can be obtained by increasing the accuracy of the stretching and diffusion simulation which can be attained by a high number of vortex particles and high temporal resolution. Both make the method impractical at least for high Reynolds numbers. After many efforts the authors settled on the algorithm which was originally proposed by [7] and modified in [11]. This algorithm consists of the following substeps: •Calculation of the change of the vorticity strength magnitude d|ωv| dt =d√ωv·ωv dt =ωv |ωv|·dωv dt (10) where dωv dt is calculated from (7) without the viscous diffusion term. The term (ωv·∇)uv is calculated taking into account adjacent vortex particles located only within one or two layers of neighboring cells. •Calculation of the particle length from the equation of the elementary section dl transported in inviscid flow dl dt =l |ωv| d|ωv| dt (11) •Calculation of the particle core radius from the equation describing the transport of an elementary tube with length land radius σin inviscid flow dσ dt =−σ 2l dl dt (12) 4 177 S. Samarbakhsh and N. Kornev •Consideration of the viscosity influence using the core spreading method (CSM) (see [5]). The particle core radius is increased by Δσ: Δσ=√4νΔt(13) •Calculation of the new particle orientation ω∗=ωv(t)+dωv dt Δt(14) •Calculation of the particle strength magnitude from [7] |α(t+Δt)|=|α(t)|σ(t+Δt) σ(t)(15) •Calculation of the new strength vector α(t+Δt)=|α(t+Δt)|ω∗ |ω∗|(16) In the original version proposed by [7] the next step should be the redistribution of particles whose length has 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 [7] 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, as pointed out by [3], is an immanent part of every turbulence model. 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. If particles grow due to viscosity and flow stagnation and exceed some size they are mapped back to the grid. The simple Euler method is used for the integration of the differential equations. The flowchart of the whole algorithm is presented in [11]. The present method has the same error sources as every LES model [6]. Two comments should be made on the filtering errors. First, the applicants understand that the models relying on the small scales comparable with cell sizes can suffer from the filter aliasing errors inherently presented in each numerical method. For instance, such errors most strongly affect dynamic type models which rely heavily on the smallest scales to determine SGS properties [6].Second, since our algorithm does not use commutation of differencing and filtering operators which is the big difficulty in LES formalism and represents the second part of the filtering errors, the commutation error is not present. 5 178 S. Samarbakhsh and N. Kornev 3 RESULT AND DISCUSSION 3.1 Summary of previous result Validation and verification is performed for wall free flows including decaying isotropic turbulence (DIT) [4] and free turbulent jet [8, 2] test cases. Results revealed that the effect of the term uv×ωgin the equation (7) is similar to that of a LES subgrid model. The term uv×ωgbehaves as an energy drain transferring the energy of the grid based motion into the fine scale energy. At coarse resolutions, it acts as a diffusive Large Eddy Simulation subgrid model resulting in a LES-like behavior of the whole method. The additional term uv×ωgis automatically switched off when the resolution increases, the present method is consistent and converges to the Direct Numerical Simulation. The energy back scattering is also captured by the present method. As mentioned in [16], the intensification of the turbulent kinetic energy due to back scattering is proved to be very important to properly reproduce the jet breakdown and transition to turbulence close to the nozzle without any artificial turbulence forcing at the nozzle. The Reynolds stresses of the velocity field induced by particles possess the pronounced anisotropy which is space dependent. It was also shown that the model for the jet case can be sufficiently reduced by neglecting the inner interaction between particles. This results in a drastic reduction of the computational time. Some additional result on the model reduction and anisotropy Reynolds stresses are given in the next section. The result are obtain for the free jet at the Reynolds number Re = 104. More details about grid properties and set up of the simulation are thoroughly described in [16]. 3.2 Model reduction As shown in [11] the equations 7 and 8 can sufficiently be reduced by neglecting inner interaction between vortex particles without a significant loss of the simulation accuracy. The reduced equations take the form: ∂ωv ∂t +(ug·∇)ωg=(ωv·∇)ug+νΔωv,(17) and the r.h.s of the trajectory equation (18) contains only the grid based velocity dri dt =ug i(18) As seen in [11], the influence of inner interactions on spatially averaged kinetic energy and scalar dissipation rate is relatively weak and can be neglected in the calculation. Thereby the computations can be done sufficiently faster.This model is further referred to as the passive vortices model. Within the next simplification step the influence of the grid based solution on the evolution of vortex particles strengths is neglected. The equations describing the vortex particle evolution take the simplest form: ∂ωv ∂t +(ug·∇)ωg=νΔωv,(19) 6 179 S. Samarbakhsh and N. Kornev (a) influence of ignoring inner interaction between vortex particles (b) influence of considering vortex particles only in the curretn time step Figure 1: Influence of the interaction between vortex particles and grid on the evolution of urms dri dt =ug i(20) Figure 1-a demonstrates results for the r.m.s. of the axial velocity obtained using the full model equations 7 and 8, passive vortices model equations 17 and 18 and the model without influence of the grid based flow on vortex particles strengths equations 19 and 20. The difference between results is negligible pointing out that vortex particles serve just as triggers or intensifiers of turbulence and their inner interaction doesn’t contribute sufficiently to the flow evolution. Hence the name of the method is the LES intensified by the vortex particles or VπLES . In the next step of model reduction only the influence of vortex particles generated in the current time step are considered and vortex particles generated in the previous time step were mapped back to the grid. In this case the time of computations is less than full model since we only deal with a vortex particles in the current time step. As can be seen in Figure 1-b the result shows that the vortex particles trigger the turbulence and have a good agreement with experiment in the near jet exit region while in far field region decay of kinetic energy is not physical. It can be interpreted in this way that the energy drain of fine scale motion from large scale motion is not high enough and accumulation of kinetic energy on the grid flow motion happens. Concluding, not only new generated vortices but also the whole set of vortices including those generated upstream in previous time steps have a significant influence on the turbulence development downstream. With other words, the fine scale vortices model (7) and (8) can be reduces but not neglected. 3.3 Anisotropy of fine scale structures induced by vortex particles Since a deterministic prediction of a turbulent flow as mentioned in the definition of the turbulent motion by Lesieur [13] 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 7 180 S. Samarbakhsh and N. Kornev (a) total velocity (ug+uv)(b) fine scale velocity uv Figure 2: Distribution of the diagonal Reynolds stresses components for the total and fine scale velocity along the jet axis subgrid motion should be captured by a proper subgrid model: non-equilibrium effects including laminar-turbulent zone, energy backscatter and anisotropy of fine scale motion. In this subsection the anisotropy of velocities induced by fine vortices is discussed. The total flow shows a well pronounced anisotropy with the dominance of the axial fluctuations on the jet centerline (see Fig. 2-a). Reynolds stresses Rv ii of the velocity field uvshows also a clear anisotropy which is space-dependent. On the center line at x/D > 5 two diagonal stresses are equal to each other Rv 22 ≈Rv 33 and dominate over Rv 11(see Fig. 2-b). To explain this effect, we consider stochastic distribution of the statistically independent axis-symmetric vortex particles on the centerline with strengths aligned with the x-axis. They induce velocities uv x= 0 and uv y�= 0, uv z�= 0. Due to axis symmetrical character of each vortex the spatially averaged squares of velocities uv yand uv zare equal, i.e. (uv y)2= (uv z)2. Precession of vortices around their spins causes the appearance of the longitudinal velocities uv xwhich are much smaller than uv yand uv z. 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 r/D =0.25 the fine scale turbulence is nearly isotropic Rv 11 ≈Rv 22 ≈Rv 33 in the beginning of the jet development (see Fig. 3-a), i.e. this area is populated with vortex particles with orientations uniformly distributed around a sphere. Further downstream the same anisotropy takes place as that on the jet axis. At the jet boundary the fine scale turbulence becomes anisotropic with a clear dominance of the radial fluctuations Rv 22 >R v 11 ≈Rv 33 (see Fig. 3-b), i.e. the dominating fluctuations are in the direction of the dominating large scale entrainment motion. Figures 4 shows the p.d.f of vortex axes orientation where Lx, Lyand Lzare defined as Lx =αx |α|,Ly =αy |α|,Lz =αz |α|. 4 CONCLUSION The paper presents validation and verification study of a novel VπLES method which is based on the decomposition of the flow structures in large scale ones, resolved on the grid, 8 181 S. Samarbakhsh and N. Kornev (a) (b) Figure 3: Distribution of the diagonal Reynolds stresses components for the fine scale velocity uvalong the line r/D =0.25 and r/D =0.5 (a) (b) Figure 4: p.d.f of vortex axes orientation line r/D =0 9 182