- Email: [email protected]
Two-dimensional viscous vortex flow around a circular cylinder
Aerosp. Sci. Technol. 5 (2001) 85–94 2001 Éditions scientifiques et médicales Elsevier SAS. All rights reserved S1270-9638(00)01082-8/FLA
Two-dimensional viscous vortex flow around a circular cylinder Sébastien Rouvreau a,1 , Laurent Perault b a Laboratoire de Combustion et de Détonique (UPR CNRS 9028), ENSMA, BP 40109, 86961 Futuroscope-Chasseneuil cedex,
France
b Laboratoire d’Etudes Aérodynamiques (UMR CNRS 6609), ENSMA, BP 40109, 86961 Futuroscope-Chasseneuil cedex, France
Received 27 March 2000; revised 7 September 2000; accepted 25 October 2000
Abstract
A code based on a vortex method for simulating viscous vortex flow around a circular cylinder has been developed. The original method to obtain the no-through and no-slip condition introduced here is described in detail. The diffusive part of the vorticity transport equation is treated using a deterministic method based on the solution of the heat equation. Results have been obtained and are presented concerning the vorticity field, the velocity field, the evolution of drag and lift coefficients and the Strouhal number. 2001 Éditions scientifiques et médicales Elsevier SAS vortex method / viscous / cyclinder / no-slip
Résumé
Un code de calcul d’écoulement tourbillonnaire de fluide visqueux autour d’un cylindre, fondé sur une méthode particulaire a été développé. La méthode originale d’obtention de la condition d’adhérence mise en place ici est décrite en détail. Le traitement de la partie diffusive de l’équation de transport de la vorticité utilise, lui, une méthode déterministe fondée sur la solution de l’équation de la chaleur. Les résultats sont présentés pour l’évolution du champ de vorticité, du champ de vitesse, des efforts sur le profil et le nombre de Strouhal. 2001 Éditions scientifiques et médicales Elsevier SAS vorticité / visqueux / cylindre / adhérence
Nomenclature Cx Cz
drag coefficient lift coefficient
L n p
unit normal vector perpendicular to the planar domain reference length (diameter of the cylinder) vortices shedding frequency dimensionless pressure such that p = p∗2
→ ez
p∗ Re
ρU
pressure of the fluid Reynolds number defined as Re = UνL
St
Strouhal number, St = nL U
t t∗
dimensionless time L characteristic time such that t ∗ = U
u U v x
component of the velocity vector in x direction reference velocity (upwind flow velocity) component of the velocity vector in y direction longitudinal direction
→
x y ν
position vector transversal direction kinematic viscosity of the fluid
1 Correspondence and reprints. Present address: University of Maryland, Fire Protection Engineering Department, College Park, MD 20742-0001, USA. Email address: [email protected]
86
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
ρ
density of the fluid
ω
vorticity vector
Ω (t)
solution for the diffusive part of the transport equation
→ →
1. Introduction Many incompressible flows are characterised by rotational regions in a mainly irrotational flow. Vortex methods simulate such flows by discretising these concentrated rotational zones into vortex elements. These particle methods, where particles are advected in a Lagrangian way, and which are being developed for threedimensional studies, are particularly well adapted for incompressible two-dimensional instationnary flows dominated by convective effects. → → → The main variable is vorticity ω = × u , defined as the rotational of the velocity vector. For a two→ dimensional study, vorticity becomes a scalar ( ω = (0, 0, ω)), and the only non-zero component of the vector rotational can be written ω=
∂v ∂u − . ∂x ∂y
These methods, first used by Chorin [3], have been significantly developed during the last few years [7,9,12] and now constitute a family of methods with different approaches. The most used approach remains close to the Chorin’s exclusively Lagrangian method, where particles are displaced by convection and viscous diffusion in time (Random Walk Method). In 1986, Cottet and Gallic [4] and Choquin and Huberson [2] concurrently developed a method to deal with the diffusion term in a deterministic manner. For this method, when calculating the diffusion, it is not the particles’ position that is changed but their weight. Other widely used approaches are the Euler–Lagrange approach and the semi-Lagrangian approach (in the Vortex-In-Cell methods). Both of these approaches require a grid. 2. Mathematical formulation Viscous incompressible flows are governed by dimensionless Navier–Stokes equations (mass and momentum conservation) that can be written in 2D: ∂u ∂v + = 0, ∂x ∂y ∂u ∂u ∂p 1 ∂ 2u ∂ 2u ∂u +u +v =− + + , ∂t ∂x ∂y ∂x Re ∂x 2 ∂y 2 ∂v ∂v ∂p 1 ∂ 2v ∂ 2v ∂v +u +v =− + + 2 . ∂t ∂x ∂y ∂y Re ∂x2 ∂y
The boundary conditions associated with these equations include the no-slip condition on solid walls, which is the physical constraint that generates vorticity, and the no-penetration condition (i.e. zero normal speed on solid walls). Since it is impossible to get an analytical solution for this non-linear system, the use of a numerical scheme is necessary. Taking the rotational of the Navier–Stokes equations (divergence of the velocity vector being zero) and intro→ ducing vorticity ω leads to the two-dimensional vorticity transport equation: →
1 → ∂ω → → → + ( u · ∇) ω = ω . ∂t Re The main difficulty concerning time discretisation is to → treat the temporal derivative of ω . The method used here splits the equation and is known as the fractionated-stepmethod. The equation is split into two steps: one of advection and one of diffusion, solved simultaneously: →
∂ω → → → = −( u · ∇ ) ω , ∂t →
∂ω 1 → = ω . ∂t Re Proves of the convergence of these equations to the solution of the transport equation do exist [1,13,14]. The vorticity field is discretised into vortex elements that will be advected and diffused throughout time. Advective transport for an element is governed by a first order differential equation that uses the position of the element and its velocity induced by the vorticity field. The way this velocity is obtained depends on the method and an advective displacement is associated to it. Diffusive transport with the Random Walk method, is statistically simulated with Gaussian random displacements to give a diffusive displacement. Concerning the deterministic method used here, diffusive transport is simulated by varying particles weight, in other words their vortical intensity.
3. The deterministic diffusion model This method, developed firstly by Cottet and Gallic [4] and Choquin and Huberson [2] treats the diffusion equation in a deterministic way by using the exact solution of the heat equation. The conservative algorithm proposed by Choquin and Huberson, on which mathematical analysis has been made by Cottet and Gallic, takes into account the vorticity diffused from each particle. In this method, concerning the diffusion simulation, the vortex intensity varies, not the particle position. The solution for diffusive
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
part of the transport equation of the vorticity is then given by: →
Re Ω (t) = 4πt
87
With this equation, the vortex contained in a particle ‘P ’ can be calculated, and then, taking into account the vortex diffused out of each particle leads to: → → 2 → → Re |x − x | exp − Ω (δt) = Ω0 + 4πδt 4Re−1 δt
→ → 2 | x − x | → → →
ω0 (x ) d x . exp − 4Re−1 t
R2
→
→
P R 2 −P → →
→
R 2 −P P → →
× ω0 (x ) d x d x → → 2 Re |x − x | ···− exp − 4πδt 4Re−1 δt →
× ω0 (x ) d x d x . The second term of the right hand side of this equation represents vorticity entering the particle and the third one represents the vorticity leaving the particle. 4. Vorticity generation: no-slip condition 4.1. Conventional methods On solid walls, the velocity vector must be zero, which gives the no-slip condition. However, during a perfectfluid type calculation, this condition is not obtained, giving a tangential velocity called the slip-velocity ug . Thus, for each time step, an appropriate quantity of circulation must be generated on solid walls to obtain
Figure 1. Schematic of the method used to obtain the no slip condition on solid walls for one control point.
(a) Figure 2. Vorticity field, Re = 1000, black: clockwise, grey: counterclockwise: (a) t ∗ = 2.0; (b) t ∗ = 6.0; (c) t ∗ = 10.0; (d) t ∗ = 12.0; (e) t ∗ = 20.0.
88
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
(b)
(c) Figure 2. (Continued)
a zero slip-velocity, which then simulates the no-slip phenomenon due to fluid viscosity [5]. In the commonly used method to generate circulation, known as the Vortex Sheet Method, segment-vortices are emitted by solid walls and need to be transformed into vortex blobs when going out of a user-defined boundary layer called vortex sheet. The inverse transformation
also has to be performed when a blob enters the vortex sheet. This method relies on an arbitrary parameter (the vortex sheet thickness) which can be difficult to compute and increases the calculation cost, since at each time step the transformation of vortex elements (segment or blob) passing through the vortex sheet frontier must be computed.
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
89
(d)
(e) Figure 2. (Continued)
In other methods, that do not use the vortex sheet, two steps are commonly used to obtain the velocity boundary condition at the solid surface. One set of blobs is emitted to obtain the zero-normal velocity at the surface and then another one is created to obtain the zero-tangential velocity at the surface. It seems hazardous to use such a method without an appropriate iterative process since the creation of the second set of blobs must modify the
velocity field and then a normal velocity may appear at the solid surface. 4.2. A one-step method To avoid an iterative process or the use of the vortex sheet method that both induce time cost it is then of
90
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
(a)
(b) Figure 3. Velocity vectors of passive particles, Re = 1000: (a) t ∗ = 2.0; (b) t ∗ = 6.0; (c) t ∗ = 10.0; (d) t ∗ = 12.0.
interest to try to find a method to obtain this boundary condition in one step. In our method, solid walls are discretised into wallsegments with a control point in the middle of each one. In order to obtain the no-slip condition we use the association of one vortex-blob and one source-segment for each control point. The source is placed inside the profile close to the wall-segment and parallel to it
whereas the vortex blob is outside of the profile, straight above the control-point ( figure 1). In that position, the vortex blob and the source-segment induce tangential and normal velocities respectively at the control-point. This association enables a one step process to obtain the velocity boundary condition on a solid wall and avoid the use of an iterative process or the calculation of vortex elements passing through a vortex sheet frontier.
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
91
(c)
(d) Figure 3. (Continued)
Then, solving a single linear system of 2N equations at each time step gives directly the no-slip condition for an obstacle discretised into N wall-segments. Blobs are then advected, due to the velocity field and the new intensity of each blob due to the diffusion is calculated. Then, considering the resulting new velocity field, a new set of vortices is emitted and new intensities for the sources are calculated to keep the no-slip condition.
5. Drag and lift coefficients The force applied on the body by the fluid (with unit density) can be easily calculated in the case of a flow around a non-rotating circular cylinder [6]. It is given by: → Fb=−
d dt
→
→
ω ∧ x dσ.
fluid + body
92
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
(a)
(b)
(c)
(d)
Figure 4. Experimental visualisation of a flow around a circular cylinder, Re = 1000: (a) t ∗ = 2.0; (b) t ∗ = 6.0; (c) t ∗ = 10.0; (d) t ∗ = 12.0.
It is then easy to get the lift coefficient Cz and the drag coefficient Cx . 6. Results 6.1. Vorticity and velocity fields Cylinder has been chosen for this study in order to compare calculation results to some experiments conducted at Université de Poitiers in the last decade [8]. Simulation has been performed for cylinders discretised into 50 segments and for a Reynolds number of 1000 for the first twenty dimensionless seconds of the flow. The Reynolds number is defined as Re = UνL , where U is the velocity of the incoming flow, L is the diameter of the cylinder and ν is the kinematic velocity of the fluid. Concerning the velocity field and the vorticity field, a good accordance has been found between experimental and numerical results. The formation of two symmetrical vortices is observed ( figures 2(a) and 3(a) ) as in the experiment ( figure 4(a) ), creating a wake of about one diameter up to 3.0 < t < 4.0. Then, dissymmetry gradu-
ally takes place ( figures 2(b), 3(b) and 4(b) ), leading to the beginning of the periodic regime ( figures 2(c), 2(d), 3(c), 3(d), 4(c), 4(d) ). Observation of the trailing edge of the cylinder also shows an evolution in accordance with the literature ( figure 5 ). First, two primary vortex zones are created due to the impulsive start . Then contra-rotating layers appear between these vortex zones and the surface of the body, develop and make primary zones to separate from the surface. At the same time a third set of vortex zones appear, develop and join the primary vortex zone, making it grow until separation. This process gradually leads to the creation of the Von Kármán vortex street ( figure 2(e) ). 6.2. Evolution of lift and drag See figure 6. This evolution can be divided into two parts. The first is from t = 0 to t = 10. Here the drag obtained by calculation is characteristic of the flow under study. Initially, Cx is large, due to the impulsive start. Then the drag quickly decreases until t = 1.0 and then tends slowly towards Cx = 1.0, which is a
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
Figure 5. Velocity vectors of passive particles, t ∗ = 2.0, 3.0, 4.0 and 5.0, Re = 1000. Zoom on the trailing edge of the cylinder.
Figure 6. Evolution of drag Cx and lift Cz for 0.0 < t ∗ < 20.0 for a flow around a circular cylinder, Re = 1000.
93
94
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
characteristic value for a flow past a circular cylinder at Re ≈ 1000 [10]. As it can be expected from a flow around a circular cylinder, the mean lift is zero in the first period of this part. It then increases slowly and drops suddenly around t = 10, when the Von Kármán vortex street takes place. Beyond this point the drag is no longer steady and the lift oscillates widely and regularly. The frequency of this oscillations leads to St ≈ 0.22 which is a characteristic value of the Strouhal number for a flow around a circular cylinder with Re = 1000 [11]. 7. Conclusion In this paper a particle method, based on a very simple model for the no-slip condition, coupled with a well known deterministic diffusion model, is introduced to simulate a viscous flow past a circular cylinder. Results are in good correlation with both literature and experiments. Qualitatively, the vorticity field and the velocity field are in good accordance with experimental data. Moreover, concerning quantitative results, the body forces evolution and the Strouhal number arising from it correlate well with literature. The one step method used here to obtain the no-slip condition then proved correct. Since it is very simple and easily adaptable to more complicated geometries, this method should be of significant interest for particle simulations of flows around obstacles such as wings with flaps. Now that this method is proven to give qualitatively as well as quantitatively good results for the first 20 seconds of the flow, future work will include an evaluation of the advantage concerning time cost with using this method. Then, the following step will be to try to extend this method to 3d simulations.
References [1] Beale J.T., Majda A.J., Rates of convergence for viscous splitting for the Navier–Stokes equations, Math. Comput. 37 (1981) 243–259. [2] Choquin J.P., Huberson S., Particles simulation of viscous flow, Comput. Fluids 17 (2) (1989) 97–410. [3] Chorin A.J., Numerical study of slightly viscous flow, J. Fluid Mech. 57 (4) (1973) 785–796. [4] Cottet G.H., Gallic S., Rapport interne 115, CMAP, École Polytechnique, 1986. [5] Ghoniem A.F., Gagnon Y., Vortex simulation of laminar recirculating flow, J. Comput. Phys. 68 (2) (1987). [6] Koumoutsakos P., Leaonard A., High-resolution simulation of the flow around an impulsively started cylinder using vortex methods, J. Fluid Mech. 26 (1995) 1–38. [7] Leonard A., Vortex methods for flow simulation, J. Comput. Phys. 37 (1980) 289–335. [8] Pineau G., Mise en évidence et évaluation des effets tridimensionnels dans l’écoulement instationnaire autour d’un cylindre circulaire d’envergure finie, PhD thesis University of Poitiers, 1992. [9] Puckett E.G., Vortex methods: a survey of selected research topics, in: Nicolaides R.A., Gunzburger M.D. (Eds.), Incompressible Computational Fluid Dynamics – Trends and Advances, Cambridge University Press, 1991. [10] Roberson J.A., Crowe C.T., Engineering Fluid Mechanics, 6th edition, Wiley, 1997, pp. 429–436. [11] Roberson J.A., Crowe C.T., Engineering Fluid Mechanics, 6th edition, Wiley, 1997, pp. 436–437. [12] Sarpkaya T., Computational methods with vortices – The 1988 Freeman scholar lecture, J. Fluid. Eng.-T. ASME 111 (1989). [13] Strang G., Accurate partial difference method II. Nonlinear problems, Numer. Math. 6 (1964) 37–49.
Acknowledgements The authors would like to thank M. Gérard Pineau for the experimental data and for his precious contribution.
[14] Strang G., On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968) 506–517.
Two-dimensional viscous vortex flow around a circular cylinder Sébastien Rouvreau a,1 , Laurent Perault b a Laboratoire de Combustion et de Détonique (UPR CNRS 9028), ENSMA, BP 40109, 86961 Futuroscope-Chasseneuil cedex,
France
b Laboratoire d’Etudes Aérodynamiques (UMR CNRS 6609), ENSMA, BP 40109, 86961 Futuroscope-Chasseneuil cedex, France
Received 27 March 2000; revised 7 September 2000; accepted 25 October 2000
Abstract
A code based on a vortex method for simulating viscous vortex flow around a circular cylinder has been developed. The original method to obtain the no-through and no-slip condition introduced here is described in detail. The diffusive part of the vorticity transport equation is treated using a deterministic method based on the solution of the heat equation. Results have been obtained and are presented concerning the vorticity field, the velocity field, the evolution of drag and lift coefficients and the Strouhal number. 2001 Éditions scientifiques et médicales Elsevier SAS vortex method / viscous / cyclinder / no-slip
Résumé
Un code de calcul d’écoulement tourbillonnaire de fluide visqueux autour d’un cylindre, fondé sur une méthode particulaire a été développé. La méthode originale d’obtention de la condition d’adhérence mise en place ici est décrite en détail. Le traitement de la partie diffusive de l’équation de transport de la vorticité utilise, lui, une méthode déterministe fondée sur la solution de l’équation de la chaleur. Les résultats sont présentés pour l’évolution du champ de vorticité, du champ de vitesse, des efforts sur le profil et le nombre de Strouhal. 2001 Éditions scientifiques et médicales Elsevier SAS vorticité / visqueux / cylindre / adhérence
Nomenclature Cx Cz
drag coefficient lift coefficient
L n p
unit normal vector perpendicular to the planar domain reference length (diameter of the cylinder) vortices shedding frequency dimensionless pressure such that p = p∗2
→ ez
p∗ Re
ρU
pressure of the fluid Reynolds number defined as Re = UνL
St
Strouhal number, St = nL U
t t∗
dimensionless time L characteristic time such that t ∗ = U
u U v x
component of the velocity vector in x direction reference velocity (upwind flow velocity) component of the velocity vector in y direction longitudinal direction
→
x y ν
position vector transversal direction kinematic viscosity of the fluid
1 Correspondence and reprints. Present address: University of Maryland, Fire Protection Engineering Department, College Park, MD 20742-0001, USA. Email address: [email protected]
86
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
ρ
density of the fluid
ω
vorticity vector
Ω (t)
solution for the diffusive part of the transport equation
→ →
1. Introduction Many incompressible flows are characterised by rotational regions in a mainly irrotational flow. Vortex methods simulate such flows by discretising these concentrated rotational zones into vortex elements. These particle methods, where particles are advected in a Lagrangian way, and which are being developed for threedimensional studies, are particularly well adapted for incompressible two-dimensional instationnary flows dominated by convective effects. → → → The main variable is vorticity ω = × u , defined as the rotational of the velocity vector. For a two→ dimensional study, vorticity becomes a scalar ( ω = (0, 0, ω)), and the only non-zero component of the vector rotational can be written ω=
∂v ∂u − . ∂x ∂y
These methods, first used by Chorin [3], have been significantly developed during the last few years [7,9,12] and now constitute a family of methods with different approaches. The most used approach remains close to the Chorin’s exclusively Lagrangian method, where particles are displaced by convection and viscous diffusion in time (Random Walk Method). In 1986, Cottet and Gallic [4] and Choquin and Huberson [2] concurrently developed a method to deal with the diffusion term in a deterministic manner. For this method, when calculating the diffusion, it is not the particles’ position that is changed but their weight. Other widely used approaches are the Euler–Lagrange approach and the semi-Lagrangian approach (in the Vortex-In-Cell methods). Both of these approaches require a grid. 2. Mathematical formulation Viscous incompressible flows are governed by dimensionless Navier–Stokes equations (mass and momentum conservation) that can be written in 2D: ∂u ∂v + = 0, ∂x ∂y ∂u ∂u ∂p 1 ∂ 2u ∂ 2u ∂u +u +v =− + + , ∂t ∂x ∂y ∂x Re ∂x 2 ∂y 2 ∂v ∂v ∂p 1 ∂ 2v ∂ 2v ∂v +u +v =− + + 2 . ∂t ∂x ∂y ∂y Re ∂x2 ∂y
The boundary conditions associated with these equations include the no-slip condition on solid walls, which is the physical constraint that generates vorticity, and the no-penetration condition (i.e. zero normal speed on solid walls). Since it is impossible to get an analytical solution for this non-linear system, the use of a numerical scheme is necessary. Taking the rotational of the Navier–Stokes equations (divergence of the velocity vector being zero) and intro→ ducing vorticity ω leads to the two-dimensional vorticity transport equation: →
1 → ∂ω → → → + ( u · ∇) ω = ω . ∂t Re The main difficulty concerning time discretisation is to → treat the temporal derivative of ω . The method used here splits the equation and is known as the fractionated-stepmethod. The equation is split into two steps: one of advection and one of diffusion, solved simultaneously: →
∂ω → → → = −( u · ∇ ) ω , ∂t →
∂ω 1 → = ω . ∂t Re Proves of the convergence of these equations to the solution of the transport equation do exist [1,13,14]. The vorticity field is discretised into vortex elements that will be advected and diffused throughout time. Advective transport for an element is governed by a first order differential equation that uses the position of the element and its velocity induced by the vorticity field. The way this velocity is obtained depends on the method and an advective displacement is associated to it. Diffusive transport with the Random Walk method, is statistically simulated with Gaussian random displacements to give a diffusive displacement. Concerning the deterministic method used here, diffusive transport is simulated by varying particles weight, in other words their vortical intensity.
3. The deterministic diffusion model This method, developed firstly by Cottet and Gallic [4] and Choquin and Huberson [2] treats the diffusion equation in a deterministic way by using the exact solution of the heat equation. The conservative algorithm proposed by Choquin and Huberson, on which mathematical analysis has been made by Cottet and Gallic, takes into account the vorticity diffused from each particle. In this method, concerning the diffusion simulation, the vortex intensity varies, not the particle position. The solution for diffusive
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
part of the transport equation of the vorticity is then given by: →
Re Ω (t) = 4πt
87
With this equation, the vortex contained in a particle ‘P ’ can be calculated, and then, taking into account the vortex diffused out of each particle leads to: → → 2 → → Re |x − x | exp − Ω (δt) = Ω0 + 4πδt 4Re−1 δt
→ → 2 | x − x | → → →
ω0 (x ) d x . exp − 4Re−1 t
R2
→
→
P R 2 −P → →
→
R 2 −P P → →
× ω0 (x ) d x d x → → 2 Re |x − x | ···− exp − 4πδt 4Re−1 δt →
× ω0 (x ) d x d x . The second term of the right hand side of this equation represents vorticity entering the particle and the third one represents the vorticity leaving the particle. 4. Vorticity generation: no-slip condition 4.1. Conventional methods On solid walls, the velocity vector must be zero, which gives the no-slip condition. However, during a perfectfluid type calculation, this condition is not obtained, giving a tangential velocity called the slip-velocity ug . Thus, for each time step, an appropriate quantity of circulation must be generated on solid walls to obtain
Figure 1. Schematic of the method used to obtain the no slip condition on solid walls for one control point.
(a) Figure 2. Vorticity field, Re = 1000, black: clockwise, grey: counterclockwise: (a) t ∗ = 2.0; (b) t ∗ = 6.0; (c) t ∗ = 10.0; (d) t ∗ = 12.0; (e) t ∗ = 20.0.
88
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
(b)
(c) Figure 2. (Continued)
a zero slip-velocity, which then simulates the no-slip phenomenon due to fluid viscosity [5]. In the commonly used method to generate circulation, known as the Vortex Sheet Method, segment-vortices are emitted by solid walls and need to be transformed into vortex blobs when going out of a user-defined boundary layer called vortex sheet. The inverse transformation
also has to be performed when a blob enters the vortex sheet. This method relies on an arbitrary parameter (the vortex sheet thickness) which can be difficult to compute and increases the calculation cost, since at each time step the transformation of vortex elements (segment or blob) passing through the vortex sheet frontier must be computed.
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
89
(d)
(e) Figure 2. (Continued)
In other methods, that do not use the vortex sheet, two steps are commonly used to obtain the velocity boundary condition at the solid surface. One set of blobs is emitted to obtain the zero-normal velocity at the surface and then another one is created to obtain the zero-tangential velocity at the surface. It seems hazardous to use such a method without an appropriate iterative process since the creation of the second set of blobs must modify the
velocity field and then a normal velocity may appear at the solid surface. 4.2. A one-step method To avoid an iterative process or the use of the vortex sheet method that both induce time cost it is then of
90
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
(a)
(b) Figure 3. Velocity vectors of passive particles, Re = 1000: (a) t ∗ = 2.0; (b) t ∗ = 6.0; (c) t ∗ = 10.0; (d) t ∗ = 12.0.
interest to try to find a method to obtain this boundary condition in one step. In our method, solid walls are discretised into wallsegments with a control point in the middle of each one. In order to obtain the no-slip condition we use the association of one vortex-blob and one source-segment for each control point. The source is placed inside the profile close to the wall-segment and parallel to it
whereas the vortex blob is outside of the profile, straight above the control-point ( figure 1). In that position, the vortex blob and the source-segment induce tangential and normal velocities respectively at the control-point. This association enables a one step process to obtain the velocity boundary condition on a solid wall and avoid the use of an iterative process or the calculation of vortex elements passing through a vortex sheet frontier.
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
91
(c)
(d) Figure 3. (Continued)
Then, solving a single linear system of 2N equations at each time step gives directly the no-slip condition for an obstacle discretised into N wall-segments. Blobs are then advected, due to the velocity field and the new intensity of each blob due to the diffusion is calculated. Then, considering the resulting new velocity field, a new set of vortices is emitted and new intensities for the sources are calculated to keep the no-slip condition.
5. Drag and lift coefficients The force applied on the body by the fluid (with unit density) can be easily calculated in the case of a flow around a non-rotating circular cylinder [6]. It is given by: → Fb=−
d dt
→
→
ω ∧ x dσ.
fluid + body
92
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
(a)
(b)
(c)
(d)
Figure 4. Experimental visualisation of a flow around a circular cylinder, Re = 1000: (a) t ∗ = 2.0; (b) t ∗ = 6.0; (c) t ∗ = 10.0; (d) t ∗ = 12.0.
It is then easy to get the lift coefficient Cz and the drag coefficient Cx . 6. Results 6.1. Vorticity and velocity fields Cylinder has been chosen for this study in order to compare calculation results to some experiments conducted at Université de Poitiers in the last decade [8]. Simulation has been performed for cylinders discretised into 50 segments and for a Reynolds number of 1000 for the first twenty dimensionless seconds of the flow. The Reynolds number is defined as Re = UνL , where U is the velocity of the incoming flow, L is the diameter of the cylinder and ν is the kinematic velocity of the fluid. Concerning the velocity field and the vorticity field, a good accordance has been found between experimental and numerical results. The formation of two symmetrical vortices is observed ( figures 2(a) and 3(a) ) as in the experiment ( figure 4(a) ), creating a wake of about one diameter up to 3.0 < t < 4.0. Then, dissymmetry gradu-
ally takes place ( figures 2(b), 3(b) and 4(b) ), leading to the beginning of the periodic regime ( figures 2(c), 2(d), 3(c), 3(d), 4(c), 4(d) ). Observation of the trailing edge of the cylinder also shows an evolution in accordance with the literature ( figure 5 ). First, two primary vortex zones are created due to the impulsive start . Then contra-rotating layers appear between these vortex zones and the surface of the body, develop and make primary zones to separate from the surface. At the same time a third set of vortex zones appear, develop and join the primary vortex zone, making it grow until separation. This process gradually leads to the creation of the Von Kármán vortex street ( figure 2(e) ). 6.2. Evolution of lift and drag See figure 6. This evolution can be divided into two parts. The first is from t = 0 to t = 10. Here the drag obtained by calculation is characteristic of the flow under study. Initially, Cx is large, due to the impulsive start. Then the drag quickly decreases until t = 1.0 and then tends slowly towards Cx = 1.0, which is a
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
Figure 5. Velocity vectors of passive particles, t ∗ = 2.0, 3.0, 4.0 and 5.0, Re = 1000. Zoom on the trailing edge of the cylinder.
Figure 6. Evolution of drag Cx and lift Cz for 0.0 < t ∗ < 20.0 for a flow around a circular cylinder, Re = 1000.
93
94
S. Rouvreau, L. Perault / Aerosp. Sci. Technol. 5 (2001) 85–94
characteristic value for a flow past a circular cylinder at Re ≈ 1000 [10]. As it can be expected from a flow around a circular cylinder, the mean lift is zero in the first period of this part. It then increases slowly and drops suddenly around t = 10, when the Von Kármán vortex street takes place. Beyond this point the drag is no longer steady and the lift oscillates widely and regularly. The frequency of this oscillations leads to St ≈ 0.22 which is a characteristic value of the Strouhal number for a flow around a circular cylinder with Re = 1000 [11]. 7. Conclusion In this paper a particle method, based on a very simple model for the no-slip condition, coupled with a well known deterministic diffusion model, is introduced to simulate a viscous flow past a circular cylinder. Results are in good correlation with both literature and experiments. Qualitatively, the vorticity field and the velocity field are in good accordance with experimental data. Moreover, concerning quantitative results, the body forces evolution and the Strouhal number arising from it correlate well with literature. The one step method used here to obtain the no-slip condition then proved correct. Since it is very simple and easily adaptable to more complicated geometries, this method should be of significant interest for particle simulations of flows around obstacles such as wings with flaps. Now that this method is proven to give qualitatively as well as quantitatively good results for the first 20 seconds of the flow, future work will include an evaluation of the advantage concerning time cost with using this method. Then, the following step will be to try to extend this method to 3d simulations.
References [1] Beale J.T., Majda A.J., Rates of convergence for viscous splitting for the Navier–Stokes equations, Math. Comput. 37 (1981) 243–259. [2] Choquin J.P., Huberson S., Particles simulation of viscous flow, Comput. Fluids 17 (2) (1989) 97–410. [3] Chorin A.J., Numerical study of slightly viscous flow, J. Fluid Mech. 57 (4) (1973) 785–796. [4] Cottet G.H., Gallic S., Rapport interne 115, CMAP, École Polytechnique, 1986. [5] Ghoniem A.F., Gagnon Y., Vortex simulation of laminar recirculating flow, J. Comput. Phys. 68 (2) (1987). [6] Koumoutsakos P., Leaonard A., High-resolution simulation of the flow around an impulsively started cylinder using vortex methods, J. Fluid Mech. 26 (1995) 1–38. [7] Leonard A., Vortex methods for flow simulation, J. Comput. Phys. 37 (1980) 289–335. [8] Pineau G., Mise en évidence et évaluation des effets tridimensionnels dans l’écoulement instationnaire autour d’un cylindre circulaire d’envergure finie, PhD thesis University of Poitiers, 1992. [9] Puckett E.G., Vortex methods: a survey of selected research topics, in: Nicolaides R.A., Gunzburger M.D. (Eds.), Incompressible Computational Fluid Dynamics – Trends and Advances, Cambridge University Press, 1991. [10] Roberson J.A., Crowe C.T., Engineering Fluid Mechanics, 6th edition, Wiley, 1997, pp. 429–436. [11] Roberson J.A., Crowe C.T., Engineering Fluid Mechanics, 6th edition, Wiley, 1997, pp. 436–437. [12] Sarpkaya T., Computational methods with vortices – The 1988 Freeman scholar lecture, J. Fluid. Eng.-T. ASME 111 (1989). [13] Strang G., Accurate partial difference method II. Nonlinear problems, Numer. Math. 6 (1964) 37–49.
Acknowledgements The authors would like to thank M. Gérard Pineau for the experimental data and for his precious contribution.
[14] Strang G., On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968) 506–517.