Development of a hybrid grid-and particle-based numerical method for resolution of fine vortex structures in fluid mechanics

458 III International Conference on Particle-based Methods – Fundamentals and Applications PARTICLES 2013 M. Bischoff, E. O˜nate, D.R.J. Owen, E. Ramm & P. Wriggers (Eds) DEVELOPMENT OF A HYBRID GRID-AND PARTICLE-BASED NUMERICAL METHOD FOR RESOLUTION OF FINE VORTEX STRUCTURES IN FLUID MECHANICS NIKOLAI KORNEV1, VALERY ZHDANOV2, GUNNAR JACOBI1and IRINA CHERUNOVA3 1Chair of Modeling and Simulation in Mechanical Engineering and Marine Technology University of Rostock, 18059 Rostock, Germany e-mail: [email protected], web page: http://www.lemos.uni-rostock.de/ 2Chair of Technical Thermodynamics University of Rostock, 18059 Rostock, Germany e-mail: [email protected], web page: http://www.ltt.uni-rostock.de/ 3Chair of the Modelling and Design Don State Technical University (Schakhty branch) Shevchenko Str. 147, 346500 Shakhty, Russia e-mail: i sc[email protected] Key words: grid based method, vortex method, vortex dynamics Abstract. The paper presents a novel hybrid approach developed to improve the resolution of concentrated vortices in computational fluid mechanics. The method is based on combination of a grid based and the grid free computational vortex (CVM) methods. The large scale flow structures are simulated on the grid whereas the concentrated structures are modeled using CVM. Due to this combination the advantages of both methods are strengthened whereas the disadvantages are diminished. The procedures of the separation of small concentrated vortices from the large scale vortices 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 represented large structures if their size grows due to diffusion. Algorithmic aspects of the hybrid method are discussed. 1 INTRODUCTION Insufficient resolution of vortex structures is one of the key problems in Computational Fluid Dynamics (CFD). In our recent paper [1] we propose the hybrid gridand particle based method, based on a combination of the finite volume and computational vortex 1 Development of a hybrid gridand particlebased numerical method for resolution of fine vortex structures in fluid mechanics 459 N. Kornev, V. Zhdanov, G. Jacobi and I. Cherunova element [2] methods. In this paper we continue to describe the main principles of the new method not presented in [1]. In fact, the idea to combine the grid based and vortex methods is not quite new. As noted in [2]: ”The motivation for such techniques stems from the observation that the strengths and the weaknesses of grid-based and vortex schemes can be seen as complementary, depending on the physical problem”. Guermond et al. [3] proposed domaindecomposition technique based on subdividing the computational domain into two overlapping regions. The grid based method is applied in the region close to the body whereas the vortex method is utilized in the wake region. The matching of the solution between two domains is attained by Schwarz alternating method. The domaindecomposition technique is implemented in pure velocityvorticity and velocity-pressure and velocityvorticity [4] formulations. This idea is also well known for experts working in industrial aerodynamics. For instance, the grid based Navier Stokes solvers, applied for the calculation of flow around airfoils, are coupled with potential vortex method, used for simulation of tip vortex dynamics in the far field. All these works use decomposition of the domain in near and far fields. The wall bounded flows around the bodies are smooth and don’t contain large scale vortex structures. To reproduce accurately the flow gradients in the boundary layer the resolution should be very high. The wall bounded flows are modeled using the grid based techniques due to following reasons. Though the procedures for fulfillment of boundary conditions on solid walls within vortex schemes based on original works of Wu [5] has been developed for a long time (see [6]), there are still remaining many problems of algorithmic and principal characters. Among them are high costs of boundary element procedures necessary to calculate the vortex sheet on solid walls, artificial noise caused by discretization of continuous vortex fields through elements which results in the spurious turbulence close to the wall, etc. This is the reason why the most of vortex method computations especially in three dimensions have been performed for boundary free flows. A few impressive three dimensional calculations of complex geometries can be found in papers of Bernard et al. [7] and Kamemoto et al. [8]. Proper resolution of the concentrated vortex structures on grid is a big challenge for grid based method and can be easily done using their explicit representation within the vortex schemes. Within the vortexincell (VIC) method referred also to as the hybrid method, the grid is utilized for fast calculation of the velocities and for remeshing. First, due to application of the Poisson equation and Fast Fourier Transformation (FFT) instead of direct summation using the Biot-Savart integral the computations are sufficiently accelerated, especially when combined with Fast Multipole Method (FMM) for determination of boundary conditions. Second, and it is probably more important, the instability of the numerical simulation is sufficiently damped by using the remeshing procedure resulting in the redistribution of irregularly located vortex elements onto regular grid. In both cases the grid caused numerical diffusion, which absence is considered as the main advantage of vortex methods, is involved to reduce the stochastization of numerical solution. In 2 460 N. Kornev, V. Zhdanov, G. Jacobi and I. Cherunova some VIC versions the grid is also used for calculation of convection and diffusion operators. Although the grid introduction the VIC method is classified as Lagrangian or semi-Lagrangian approach since the vortices are tracked in Lagrangian way. However, the loss of the most important advantages of pure Lagrangian methods, i.e. grid independency, raises the big question about the efficiency and competitiveness of VIC with respect to common grid based methods. The present method differs principally from all vortex and hybrid methods mentioned above. It is based on the decomposition of the velocity and vorticity fields into the distributed large scale and concentrated small scale fields. The large scale field is represented on the grid, whereas the small scale one is calculated using the grid-free computational vortex method. The domain decomposition is not applied. The method is pure Lagrangian one for small structures and pure grid based one for large scale structures. The simulation with CVM is embedded into the grid simulation. There exist a permanent exchange between grid and particle represented vortices. Since we use the formalism different from the classical CVM method many of its weaknesses become irrelevant. Generation of the set of vortex elements from the grid distributed vortices is described in sections 2.1 and 2.2 in [1]. The governing equations are derived in sec. 2.3 [1] and described shortly below along with boundary condtions for the grid based solution. In this paper we present experimental data supporting the idea of the method and address a very important question concerning the interaction between scales. 2 SOME RECENT EXPERIMENTAL RESULTS SUPPORTING CONCEPT OF THE HYBRID METHOD Main concept of the hybrid method is that the fine scale vorticity in full developed turbulent flows is concentrated in a finite number of vortices which can be represented by single axisymmetric vortex elements. In this subsection we prove this concept using high resolved PIV measurements data. The flow under consideration is the turbulent axisymmetric jet developing in a coflow confined by a pipe of diameter D= 50mm and length 5000mm schematically given in Fig. 1. Medium in both flows is water. The inner tube had diameter d= 10mm and the length 600mm chosen from the condition that perturbations caused by the knee bend are suppressed near the nozzle exit. The test section of the mixer was installed in a Perspex rectangular box filled with water to reduce refraction effects. More detailed information about the hydrodynamic channel can be found in [9]. Since the Reynolds number based on the jet exit velocity Udis Red=dUd/ν = 104the jet can be considered as a fullydeveloped turbulent jet. PIV measurements were performed within the window 3.232mm×2.407mm with pixel distance of ∆ = 68.8µm. The laser thickness estimated as ∼40µm is very thin. The measurement window was located on the centerline of the jet mixer at the distances x/D = 1 and 7 from the nozzle. The vorticity was calculated using the central differential scheme (CDS). The snapshots of the vorticity component squared ω2 z/<ω 2 z>, where <> stands for quantity averaged over the window, is shown in Fig.2. Strong uneven distribution of ω2 zpointed 3 461 N. Kornev, V. Zhdanov, G. Jacobi and I. Cherunova Figure 1: Sketch of the flow. 1knee bend of nozzle, 2plate for damping of vortices shed from knee bend 1, 3outer tube, 4support plates, 5nozzle, 6test section, 7water box. Figure 2: Snapshot of the field ω2 z/<ω 2 z>within the measurement window in jet mixer. The averaged <ω 2 z>was 1.19s−2and 0.459s−2at, respectively, x/D = 1 and 7. clearly out, that the vorticity is concentrated in a relatively small number of spots or vortices. This character of the distribution is typical for both initial development of the confined jet at x/D =1.0 with weak anisotropy (R11 =0.092,R22 =0.0732) and in the region of its strong decay at x/D =7.0 where the flow is almost isotropic (R11 =0.0362, R22 =0.0352). Strong concentration of vorticity is especially obvious in Fig. 3. The cells are sorted in order of descend of ω2 z, i.e. the first cell has the maximum value of ω2 zand the last one with the number Nhas the minimum value. The ratio εk= k  i=1 ω2 zi/ N  i=1 ω2 zi shows the contribution of kcells to the total amount Ω = N  i=1 ω2 zi. The ratio along the horizontal axis shows the fraction of cells containing εk. As seen from this figure, the dependence εk(k/N ) is strongly nonlinear and reaches the saturation very quickly. Five percent of cells contains more than fifty five percent of the total Ω, twenty percent of cells contains more 4 462 N. Kornev, V. Zhdanov, G. Jacobi and I. Cherunova Figure 3: Ratios εkand Ekdepending on k/N . than eighty percent, sixty percent of cells contains less than five percent of Ω. Therefore, the number of active cells and, respectively, number of active vortices are very small. Note that the background vorticity obtained as the average over the whole measurement window is of order of ∼10−3although the maximum and minimum values of a few dozens. The distribution of Ek= k  i=1 u2 i/ N  i=1 u2 iis less nonlinear indicating the fact that the distribution of the energy is more uniform than that of ω2 z. The reason of more uniform distribution is that the energy is an integral quantity, whereas ωzis the local one. The contribution to Eis carried out not only by vortices located at adjacent cells within the measurement plane but also by all vortices of the volume including vortices located outside of the measurement window. The vortices are inclined to the measurement window at different angles β. The trace of vortices on the measurement plane is ωsin β. One can assume that the maximum ω2 z corresponds to vortices which are perpendicular to the measurement plane (β=π/2). Already visual analysis of Fig. 2 suggests that the strongest vortices are approximately axisymmetric. We apply the algorithm proposed below in the section 2.2 in [1] to detect the vortex structures in the field of ω2 zusing two dimensional linear approximation of ω2 z. Note that the linear approximation is consistent with CDS applied for the calculation of ωz. Fig. 4 shows the probability density function of the structures of the field ω2 z. The most frequent structures have radius around ∼ 2.5∆. 5 463 N. Kornev, V. Zhdanov, G. Jacobi and I. Cherunova Figure 4: Probability density functions of radius of structures of the field ω2 z(left) and of the axis ratio a−b √ab of structures of the field ω2 z(right) The p.d.f of the ratio a−b √ab indicating the circularity of the vortex cross section is given in Fig. 4, right. As seen, the circular cross section corresponding to a=bis the most frequent case. A large fraction of vortices is inclined to the measurement plane. Even if they are axisymmetric their intersection with measurement plane is not circular. It means that the true number of circular vortices is sufficiently larger than these corresponding to a−b √ab = 0 in Fig.4, right. Analysis of the circularity should be done with care because the peak of a−b √ab = 0 can be just due to a low resolution. Indeed, if the real size of vortex is smaller than the cell size ∆, being identified at any node of the grid, it occupies four adjacent cells. In our algorithm such a vortex is identified as the circle with the radius of R= ∆. The circularity of vortices with the radius equal to ∆ is indefinable. To exclude their influence we calculated conditioned p.d.f. of a−b √ab at R>m∆ shown in Fig. 5. As clearly seen the peaklike character of p.d.f. in vicinity of a−b √ab = 0 is kept even at m= 4. Taking the fact into account, that the most frequent vortices have according to Fig. 4 m≈2.5, and results in Fig. 5, one can conclude that the axisymmetric approximation of fine vortices can be considered as quite appropriate. Increase of the order of spline approximation of ω2 zfield up to three doesn’t change the qualitative conclusions drawn from the bilinear approximation. These results agree with these obtained in [10] for the turbulent boundary layer. Even in the shear flow the number of active vortices is small. The vorticity of the concentrated structures is one or two orders higher (see Fig.30 in [10]) than the vorticity of the background computed by differentiation of the averaged velocity field. 6 464 N. Kornev, V. Zhdanov, G. Jacobi and I. Cherunova Figure 5: Condtioned probability density function of the axis ratio a−b √ab (R>m∆) of structures of the field ω2 z. 3 EQUATIONS AND BOUNDARY CONDITIONS (BC) OF HYBRID GRID AND PARTICLE BASED METHOD Two governing coupled equations were derived in [1] using the splitting procedure applied to the NavierStokes equation: ∂ug ∂t +(ug∇)ug=∇P+ν∆ug+(uv×ωg) (1) dωv dt =(ωv∇)(uv+ug)+ν∆ωv(2) The first equation describes the flow of the background ug iwhereas the second one the flow induced by concentrated vortex structures uv i. The equations (1) and (2) are solved sequentially. The first equation is solved on the grid with finite volume method whereas the second one using the grid free computational vortex method. According to Gresho and Sani [11] the boundary conditions for the velocity are sufficient to allow the determination of both velocity and pressure from the NS equation. The no slip condition at the wall reads uv+ug=0→ug=−uv(3) The necessary and sufficient boundary condition for the pressure, which is used in Poisson equation, is the Neumann BC obtained by the projection of the Navier Stokes equation onto the normal direction: ∂P ∂n =ν∆un−  ∂un ∂t +(u∇)un  (4) 7 465 N. Kornev, V. Zhdanov, G. Jacobi and I. Cherunova Commonly, for high Reynolds numbers the first term on r.h.s. of (4) is neglected. The Neumann BC for grid based part reads ∂P ∂n =−∂ug n ∂t +(ug∇)ug n+(uv×ωg)n(5) There are no explicit boundary conditions for the vortex part. The interaction of fine vortices with boundaries is considered in boundary conditions for the grid based solution. 4 INTERACTION BETWEEN SCALES 4.1 Influence of large grid based vortices on small vortices This influence is taken by terms (ωv∇)ugand (ug∇)ωvin Eq. (2). The large vortices represented on grid contributes to the small vortex convection, rotation and amplification. Another mechanism of the interaction between scales is the transition of vortices. The small vortices are generated from grid based vortices using the algorithm described in Sec. 2 of [1]. The small vortices can become larger due to diffusion and be mapped back to the grid. 4.2 Impact of small scales on grid based solution The impact of small structures on gird based solution is taken by the term uv×ωg into account. The physical meaning of this term can easily be explained when applying the curl operator ∇×(uv×ωg)=−(uv∇)ωg+(ωg∇)uv(6) The first term on the r.h.s. (6) describes the transport of the grid based vorticity ωgby the velocity induced by concentrated vortices uvwhereas the second term is responsible for the rotation and amplification of the grid based vorticity in field of uv. Since the vortices are getting small due to stretching of vortex lines and become invisible on the grid, two following principal questions should be addressed below: •Is the term uv×ωgnegligible? Is its influence sporadic? •The impact of this term is local. How to not lose this local impact on the grid with the mesh size ∆ being much larger than the vortex size? Let us start with the first question. The dimensionless Navier Stokes equation reads: L UT ∂ug ∂t +( ug∇)ug=∇p+1 Re∆u+L U2(uv×ωg) (7) Let us evaluate the last term. The vortex induced velocity is proportional to the vortex circulation ωvσ2and inversely proportional to the distance from the vortex center ∼∆−1 outside of the vortex core: uv∼ωvσ2∆−1(8) 8 466 N. Kornev, V. Zhdanov, G. Jacobi and I. Cherunova Assuming that the grid based vorticity is finite L Uωg∼∂u ∂x ∼O(1), the term has the following order L U2(uv×ωg)∼κ∂u ∂x2σ ∆2∆ L∼κσ ∆2∆ L The parameter κ=ωv/ωgaccording to [10] (see Fig. 30) can be of the order of ten far from the wall. Close to the wall κ≈1.0. The ratio σ/∆ can be of the order of O(1) for big variety of vortices. The ratio ∆/L is the minimum for DNS simulations ∆/L ∼Re−3/4 t. Finally we have L U2(uv×ωg)∼κRe−3/4 t(9) Since Re Retand κ1 the last term in the equation (7) is not smaller than the diffusion term and can not be neglected even if the vortices are much smaller than mesh size ∆. Another support to keep this term in Eq.(2) is the fact that the vortices create clusters in turbulent flows which influence area can be sufficiently larger than that of a single vortex. Now we turn to the second question. Analysis of terms (uv∇)ωgand (ωg∇)uvshows that the second term is much larger: −(uv∇)ωg∼uv max ∂ωg ∂x (ωg∇)uv∼uv max σωg If the linear approximation is used within the grid based method, the first term is zero (uv∇)ωg= 0. Anyway it is much smaller than the second term, since the latter is proportional to σ−1. Therefore, the local interaction between the vortices and grid based flow can be captured using the simplification ∇×(uv×ωg)∼(ωg∇)uv(10) If Kis the smoothing function of the vortex element, the following formula are valid: uv=(γ×x)K(x) (ωg∇)uv=(ωgx)(γ×x)∂K ∂r +(ωg×γ)K(x) (11) where r=√xixi. The vector (ωg∇)uvdescribes the local impact of fine vortices on grid based vorticity. It is the local change of the grid based vorticity ωgat the place of the fine vortex with the strength γ. ∂ωg ∂t =(ωgx)(γ×x)∂K ∂r +(ωg×γ)K(x) (12) 9