Vortex method with turbulence sub-grid scale modelling

Vortex method with turbulence sub-grid scale modelling

L. A. A. Pereira^I; M. H. Hirata^II; A. Silveira Neto^III

^IDepartamento de Mecânica, Universidade Federal de Itajubá - UNIFEI, 37500-903 Itajubá, MG. Brazil. [e-mail]

^IIDepartamento de Mecânica, Universidade Federal de Itajubá - UNIFEI, 37500-903 Itajubá, MG. Brazil. [e-mail]

^IIIDepartamento de Engenharia Mecânica, Universidade Federal de Uberlândia - UFU, 38401-136 Uberlândia, MG. Brazil. [e-mail]

ABSTRACT

In this paper a method for the simulation of convection and diffusion of the vorticity, generated on a body surface is presented. A purely Lagrangian scheme is used for the vorticity convection, and thereby avoiding mesh associated problems. The body surface is simulated by straight-line panels, with constant-strength vortex distribution. The strength of the discrete vortices is obtained directly without going through any additional calculation. Using a primary diffusion process this vorticity is replaced by Lamb vortices located nearby the body surface. The diffusion process of the vorticity is simulated using the random walk scheme. Turbulence sub-grid scale modelling employs a Second Order Velocity Structure Function model adapted to the Lagrangian scheme. With a cloud of discrete vortices this velocity structure function model is employed to simulate the micro structures of the flow. The flow over a circular cylinder and a NACA 0012 aerofoil are considered to evaluate the integrated aerodynamic loads. Comparisons are made with previous theoretical and experimental studies.

Keywords: Vortex method, panel method, turbulence model, sub-grid scale model

Introduction

The flow around a bluff body or even the flow around a slender body at an incidence is complex and many of its aspects can be considered as open questions deserving additional investigations. In order to overcome the difficulties of the analysis arising from the complexities of the flow, much effort is required from the investigators. The development of new techniques and a fresh approach to the solution of the problems are urgently needed.

Among others, two aspects of the flow contribute most of the difficulties found in the analysis. Firstly, the viscous wake, which develops downstream of the body, concentrates the vorticity generated on the body surface. As the important phenomena of the flow develop in the boundary layer and in the wake region, it is of fundamental importance, therefore, to understand and to be able to analyse the vortex dynamics. The second important aspect is related to the highly turbulent flow, particularly in the boundary layer for large Reynolds number, which dominates the flow behaviour. New techniques, which can handle such aspects together with the dynamics of the vorticity, must be developed.

In this paper a new approach to the problem is presented. Its main characteristic is the utilization of the Lagrangian description (instead of the Eulerian description, which constitutes the foundations of the most usual approaches such as the Finite Volume and the Finite Element Methods) in combination with a sub-grid scale modelling for turbulence.

A cloud of discrete vortices is used in order to simulate the vorticity, which is generated on the body surface and develops into the boundary layer and in the wake. Each individual free vortex of the cloud is followed during the numerical simulation in a typical Lagrangian fashion. This is in essence the foundations of the Vortex Method (Chorin, 1973; Sarpkaya, 1989; Kamemoto et al, 1995; Lewis, 1999). With the Lagrangian formulation a grid for the spatial discretization of the fluid region is not necessary. Thus, special care to handle numerical instabilities associated to high Reynolds numbers is not needed.

In the Lagrangian simulation, attention is focused on the region of high activities which are the regions containing vorticity. On the other hand, with the Eulerian schemes, the entire fluid domain is considered independent of the fact that there are regions less important in which any activity of the flow can be found.

With the Lagrangian tracking of the vortices, one need not take into account the far away boundary conditions. This is very important in the wake region, which is not negligible in the flows of present interest, see "Fig. 1", where turbulence activities are intense and unknown, a priori.

The operation counts with the Vortex Method is of N², N being the number of free vortices used in the numerical simulation. The so-called Multipole Expansions schemes allow the reduction of this operation count. The "particle-box algorithm" (Guedes et al., 1999) reduces the operation count to Nln(N) and the "box-box algorithm" (Greengard and Rokhlin, 1987), reduces the operation count to N. Nevertheless, the current computational resources make the direct numerical simulation of turbulence prohibitively expensive, due to the wide range between the largest scales of interest and the smallest dissipative scales.

It is appropriate to mention the existence of hybrid methods that, although keeping the Lagrangian tracking of the free vortices, utilise a grid, in order to calculate the velocity induced at the locations of the free vortices. These are the so-called Vortex-in-Cell Methods (Christiansen, 1973; Walther et al., 1999) which are not considered here.

It is, therefore, necessary to keep the number of free vortices within a manageable range. For this, an analysis of the flow phenomena, according to their characteristic scale, has been developed, since most the phenomena arise at larger scales of the flow, and yet, they are modified by the sub-grid scale phenomena (Lesieur, 1990). As a result, this leads to the analysis of the macro scale phenomena, which furnishes information about the large structures of the flow, using a manageable cloud of vortices. The micro scale phenomena are modelled as shown by Alcântara Pereira etal. (2000).

To separate the wave number scales an appropriately chosen low-pass filter, characterised by a function , is used. The filtered field is defined, for any quantity, as

Nomenclature

b = characteristic length

C_D= drag coefficient

C_k = Kolmogorov constant

C_L= lift coefficient

DRAG = drag coefficient

E = local kinetic energy spectrum

= local second order velocity structure function

= filter function

i = complex number

k_c = cut-off wave number

LIFT = lift coefficient

M = number of panels

N = number of free vortices

N_v = number of discrete vortices of the cloud found in the region defined by the distances (s₀–e) and (s₀+e), from the center of the reference vortex k

n = vector normal

p = pressure field

P = random number defined in the range 0.0 to1.0

Q = random number defined in the range 0.0 to1.0

Re = Reynolds number

Re_c = turbulent Reynolds number

r_p = position of a particle fluid

r_d = radial distance between the vortex center at the point in the flow field where the induced velocity is calculated

S_t = Strouhal number

S₁ = body surface

S₂ = far away boundary

S_ij = deformation tensor of the filtered field

T_ij = sub-grid scale tensor

U = incident flow

u = velocity vector

= filtered field or large scale field

u' = fluctuation field or the sub-grid scales

u = component of the velocity vector

u_q = tangential velocity

v = component of the velocity vector

z = complex plane

z_k= center of the reference vortex k

Greek Symbols

D = radius of sphere

Ds = panel length

Dt = time increment

G = vortice strength

a = incident angle

2e = width of the region taked as 2e =2/3s0

g = constant-strength vortex distribution

n = fluid kinematic viscosity

n_t = fluid turbulent viscosity

p = 3,14159265

r = fluid density

s₀ = molecular radius of the Lamb vortex

s_0c = corrected radius of the Lamb vortex

w = vorticity field

t = tangential vector

x = random displacement

Formulation of the Physical Problem

Consider the incompressible flow of a Newtonian fluid in a large two-dimensional domain. "Figure 1" shows the incident flow, defined by U, and the work domain with boundary S=S₁ÈS₂, S₁ being the body surface and S₂ the far away boundary.

The governing equations are given as follow

where the summation convention applies. For a complete definition of the problem, on S₁, the impenetrability and no-slip conditions are written as

n and t being, respectively, the unit normal and tangential vectors. One assumes that far away the perturbation caused by the body fades away as

According to Smagorinsky (1963) the above equations must be filtered using "Eq. (1)"

where H_Dx is the filter and (*) denotes the convolution product. The turbulence scales must be split in the large scales and in the sub-grid scales

where is the filtered field or large scale field and u' is the fluctuation field or the sub-grid scales, which are smaller than Dx.

It results that the filtered and the fluctuation fields are divergence-free and the filtered equations are written as (Lesieur, 1990)

where the generalised sub-grid scale tensor is defined as (Smagorinsk, 1963)

It is worth observing that "Eq. (10)" resembles the Reynolds equations for the mean flow, but the sub-grid scale tensor T_ij is different.

It is now proposed to numerically simulate the large structures, represented by , and to use appropriate models to represent the small-scale effects. With the use of the eddy-viscosity assumption (Boussinesq’s hypothesis; Boussinesq, 1877) to model the sub-grid scale tensor, the large structures are governed by

where n_t is the turbulent viscosity and

is the deformation tensor of the filtered field.

In order to take into account the local activity of turbulence, Métais and Lesieur (1992) considered that the small scales may not be too far from isotropy and proposed to use the local kinetic energy spectrum E(k_c) to define the eddy viscosity as

where k_c is the cut-off wave number.

Using a relation proposed by Batchelor (1967) the local kinetic energy spectrum at k_c is calculated with a local Second Order Velocity Structure Function (Lesieur and Métais, 1996)

From the Kolmogorov spectrum the eddy-viscosity can be written as a function of

where C_k =1.4 is the Kolmogorov constant.

The great computational advantage of this formulation over the Smagorinsky model, vis a vis the vortex method, is that in "Eq. (15)" the notion of velocity fluctuations (differences of velocity) is used instead of the rate of deformation (derivatives). The velocities u(x + r) are calculated over the surface of a sphere of radius D.

The filtered "Eq. (12)" are non-dimensionalized in terms of U and b (a characteristic length). The Reynolds number is defined as

The pressure terms in the non-dimensional Navier-Stokes equations can be eliminated (Batchelor, 1967) yielding the vorticity equation

One should observe that "Eq. (18)" is scalar equation whereas the vorticity vector in 2-D has only one component w. It is also worth to observe that the turbulence is essentially a 3-D phenomena and yet one is modelling it using a 2-D approach; obviously it is then assumed a 2-D turbulence.

With this procedure one are still left with important turbulence aspects and the final results are also improved. The use of 2-D turbulence may explain some numerical results that depart from the experimental values.

Numerical Simulation: The Vortex Method With a Sub-Grid Scale Modelling

The Vortex Method has been employed in this work. The method is modified to take into account the sub-grid scale phenomena, see "Eq. (18)". As already mentioned, the Vortex Method is based on a Lagrangian scheme used to describe the evolution of the vorticity in the fluid region. "Equation (18)" is then the starting point (Alcântara Pereira et al., 1999). Chorin (1973) proposed a viscous splitting algorithm, in which it is assumed that in a time increment the convection of the vorticity can be simulated independently of its diffusion. These effects are governed by

The vorticity, generated on the body surface, is simulated by a cloud of discrete vortices. The strength of individual vortices, generated in each time step, is determined using the no-slip boundary condition, see "Eq. (4)" and "Eq.(5)". Once generated, the new vortices are incorporated to the vortex cloud which, in turn, is subjected to the convection and diffusion process.

Convection is governed by "Eq. (19)", thus treated as a potential flow problem. In the complex z-plane the velocity field is given by (Sarpkaya, 1989)

Here, u and v are the x and y components of the velocity vector u and i = . The first term in the right hand side is the contribution of the incident flow; the summation of M integral terms comes from the M panels used to represent the body surface; on each panel a vortex distribution with constant density is used. The second summation is associated to the velocity induced by the cloud of N vortices; it represents the vortex-vortex interactions. In order to remove the singularities in the second summation, the potential vortices are substituted by Lamb vortices (Sherman, 1990). The circumferential velocity induced by k^th is (Mustto et al., 1998)

In this particular equation r_d is the radial distance between the vortex center and the point in the flow field where the induced velocity is calculated.

The core radius of a Lamb vortex originally recommended by Mustto et al. (1998) is modified to

The new position of each discrete vortex of the cloud is given by the second order Adams-Bashforth expression (Ferziger, 1981)

in which r_p is the position of a fluid particle; the diffusion of vorticity is governed by "Eq. (20)" and simulated using the random walk scheme (Lewis, 1991). In "Eq. (24)" the random displacement are given by

where P and Q are random numbers defined in the range 0.0 to 1.0.

The velocity structure function is now evaluated as

where N_v is the number of discrete vortices of the cloud found in the region defined by the distances (s₀–e) and (s₀+e) from the center of the reference vortex k. The velocities u_t(x + r_i) and v_t(x + r_i) are evaluated at the center of these N_v vortices. A correction is necessary due to the fact that the N_v vortices are not located at equal distances from the center of the reference vortex; see "Eq. (26)". The width of the region was taken as .

represents a local statistical average of square velocity differences between free vortices located in the region defined by the distances (s₀–e) and (s₀+e) from the center of the reference vortex. Fisically, this turbulence model represents the flow turbulent activities in the vicinity from the center of vortex located at x.

Finally the extended Blasius formulas (Sarpkaya, 1989) can be used to evaluate the integrated aerodynamic loads.

Results

The methodology presented was initially used to simulate the flow around a bluff body: a circular cylinder. "Figure 2" shows the flow pattern after 300 time steps of simulation. The pattern is illustrated by the instantaneous position of the vortices in the cloud.

"Figure 3" shows the vorticity distribution at the same instant.

The calculation of C_D (drag coefficient) and C_L (lift coefficient), as the simulation proceeds is illustrated in "Fig. 4". "Figure 4 (a)" shows the result without sub-grid scale modelling and "Fig. 4 (b)" shows the result with sub-grid scale modelling. As it can be seen, qualitatively, the behaviour of the results with sub-grid scale modelling is more regular, showing already the improvements obtained with turbulence modelling..

In the numerical simulation, the following parameters were assumed: Re=10⁵ (Reynolds number), M=32 (number of panels), Dt=0.1025 (time increment) and s₀ =0.0070 (radius of the vortex core).

The simulation was extended to other values of Re and the results are shown in Table 1. Experimental results are from Roshko (1961) and Blevins (1984) (with 10% uncertainty) and results referred to by Mustto et al. (1998) are derived using a Vortex Method without sub-grid scale modelling; these results are taken as a reference for numerical comparison since they were obtained using the Circle Theorem (Milne-Thompson, 1955), which ensure the exact satisfaction of the boundary condition on the cylinder surface.

As it can be seen from "Fig. 4" the lift coefficient oscillates around zero, as expected. However the mean value, calculated over the simulated period, is slightly different from zero. The values are, nevertheless, closer to zero when the sub-grid scale modelling is considered.

The drag coefficient shows a higher value as compared to the experimental results. This is also expected from a numerical simulation in 2-D. The numerical results are therefore acceptable for a wide range of Re. The methodology, however seems, at Re=10⁶, to be in the super critical region. More investigations are needed and one can imagine that with the use of more panels (and therefore more free vortices in the cloud) the results tend to be in closer agreement with the experiments.

The values of the Strouhal number are not in so good agreement with the experimental ones. This is probably due to two reasons; a small number of panels used to simulate the body surface and mainly due the time evolution of the loads. A new methodology for the calculation of the aerodynamics loads, which is based on an integral equation derived from a pressure Poisson equation are presently being tested and the preliminary results are much encouraging as far as the load and Strouhal number calculation.

The simulation of the flow around a slender body was also performed. The results of the numerical simulation of the flow around the NACA 0012 are presented below. The same simulation was used to derive the aerodynamic loads evolution, which is presented in "Fig. 5". From this figure the mean value, presented in Table 2, was evaluated.

In this numerical simulation, the following parameters were assumed: Re=10⁶ (molecular Reynolds number), M=32 (number of panels), Dt=0.0500 (time increment),s₀ =0.0030 (radius of the vortex core) and a =5° (incident flow).

As can be seen, the sub-grid scale modelling improved the results but the drag coefficient is still high. The lift coefficient, which is of importance, is better predicted with sub-grid scale model as compared with experimental results.

Conclusions

The main objective of the work with the implementation and initial test of a sub-grid scale model in connection with the vortex method, has been achieved. The results show that the Vortex Method with turbulence modelling, can improve previously obtained results without modelling, being therefore encouraging. Additional analysis of the influence of numerical parameters will have to be carried out. The differences encountered in the comparison of the numerical results with the experimental results are attributed mainly to the inherent three-dimensionality of the real flow for such a value of the Reynolds number, which is not modelled in the simulation. The results for the Strouhal number indicated that there was a lack of resolution, and new simulations with higher values of N are necessary. A new methodology to calculate the aerodynamics forces using an integral equation derived from the pressure Poisson equation is now under test and preliminary results are also encouraging.

Acknowledgments

This work was supported by FAPEMIG/FIEMG/IEL/EFEI that have supported the Training of Human Resources/Visitor Professor (to Dr. A. Silveira Neto) and by CNPq a Brazilian Government Agency (Proc. 300126/92–1).

Paper accepted March, 2003

Technical Editor: Clóvis Raimundo Maliska

• Alcântara Pereira, L.A., Hirata, M.H. and Oliveira, W., 1999, "Aerodynamics Loads (Numerical Simulation Using a Vortex Cloud)" (In Portuguese), 15^th Brazilian Congress of Mechanical Engineering, Águas de Lindóia, São Paulo, Brazil.
• Alcântara Pereira, L.A., Hirata, M.H. and Silveira Neto, A. 2000,"Numerical Simulation of the Newtonian Fluid Flow Using a Vortex Cloud and a Turbulence Sub-Grid Scale Modelling" (In Portuguese), IV SIMMEC, Uberlândia, Minas Gerais, Brazil.
• Batchelor, G.K. , 1967, "An Introduction to Fluid Dynamics", Cambridge University Press, Cambridge, UK.
• Blevins, R.D., 1984, "Applied fluid Dynamics Handbook", Van Nostrand Reinhold Co.
• Boussinesq, J., "Essai Sur La The Des Eaux Courantes", Memories Preesentes par Diver Savants & Agrav.
• Chorin, A.J., 1973, "Numerical Study of Slightly Viscous Flow", Journal of Fluid Mechanics, Vol. 57, pp. 785-796.
• Christiansen, J.P., 1973,"Numerical Simulation of Hydrodynamics by the Method of Point Vortices", J. Comp. Phys., Vol. 13, pp. 363-79.
• Ferziger, J.H., 1981, "Numerical Methods for Engineering Application", John Wiley & Sons, Inc.
• Guedes, V.G., Hirata, M.H. and Bodstein,G.C.R. , 1999, "Vortex Method Simulation of the Flow Around a Circular Cylinder Using the Multipole Expansion Algorithm", International Conference on Computational Heat and Mass Transfer, Gazimagusa, Cyprus of the North, Vol. 1, pp. 101-108.
• Greengard, L. and Rokhlin, V., 1987,"A Fast Algorithm for Particle Simulations", Journal of Computational Physics, Vol. 73, pp. 325-348.
• Kamemoto, K., Matsumoto, H.and Yokoi, Y., 1995, "On the Extensive Applicability of the Vortex Methods to the Prediction of Flow-Induced Vibration Problems", Yokohama National University, Yokohama, Japan.
• Lesieur, M., 1990, "Turbulence in Fluids", Kluwer, Dordrecht.
• Lesieur, M. and Métais, O., 1996, "New Trends in Large-Eddy Simulation of Turbulence", An. Review in Fluid Mechanics, Vol. 28, pp. 45-82.
• Lewis, R.I., 1991, "Vortex Element Methods for Fluid Dynamic Analysis of Engineering Systems", Cambridge University Press, Cambridge, UK.
• Lewis, R.I., 1999, "Vortex Element Methods, the Most Natural Approach to Flow Simulation - A Review of Methodology With Applications", Proc. of 1^st Int. Conf. on Vortex Methods, Kobe, Nov. 4-5, pp. 1-15.
• Métais, O. and Lesieur, M., 1992,"Spectral Large-Eddy Simulations of Isotropic and Stably-Stratified Turbulence", J. Fluid Mech. , Vol. 239, pp. 157-94.
• Miley, S.J., 1982, "A Catalogue of Low Reynolds Number Airfoil Data For Wind Turbine Applications", Department of Aerospace Engineering, Texas A&M, University, College Station, Texas, RFP-3387, UC-60.
• Milne-Thompson, L.M., 1955, "Theorical Hydrodynamics", MacMillan & Co.
• Mustto, A.A., Hirata, M.H. and Bodstein, G.C.R., 1998, "Discrete Vortex Method Simulation of the Flow Around a Circular Cylinder With and Without Rotation", AIAA Paper 98-2409, Proceedings of the 16^th AIAA Applied Aerodynamics Conference, Albuquerque, NM, USA, Vol. 1, pp. 59-69.
• Roshko, A., 1961, "Experiments on the Flow Past a Circular Cylinder at Very High Reynolds Number", Journal of Fluid Mechanics, Vol. 10, pp. 345-356.
• Sarpkaya, T., 1989, "Computational Methods With Vortices - The 1988 Freeman Scholar Lecture", Journal of Fluids Engineering, Vol. 111, pp. 5-52.
• Sherman, F.S., 1990, "Viscous flow", McGraw-Hill, International Editions Mechanical Engineering Series.
• Smagorinsky, J., 1963, "General Circulation Experiments With the Primitive Equations", Mon. Weather Rev., Vol. 91, pp. 99-164.
• Walther, J.H., Sagredo, J.T. and Koumoutsakos, P., 1999, "Simulation of Particle Flow Using Vortex Methods", Proc. of 1^st Int. Conf. on Vortex Methods, Kobe, Nov. 4-5, pp. 221-231.

Publication Dates

• Publication in this collection
  18 Mar 2004
• Date of issue
  Apr 2003

History

• Accepted
  Mar 2003