# Computation of Pressure Fields around a Two-Dimensional Circular Cylinder Using the Vortex-In-Cell and Penalization Methods

Authors: Seung-Jae Lee; Jun-Hyeok Lee; Jung-Chun Suh
Journal: Modelling and Simulation in Engineering
 (2014)

---

## Abstract

The vorticity-velocity formulation of the Navier-Stokes equations allows purely kinematical problems to
be decoupled from the pressure term, since the pressure is eliminated by applying the curl operator. The
Vortex-In-Cell (VIC) method, which is based on the vorticity-velocity formulation, offers particle-mesh
algorithms to numerically simulate flows past a solid body. The penalization method is used to enforce
boundary conditions at a body surface with a decoupling between body boundaries and computational grids.
Its main advantage is a highly efficient implementation for solid boundaries of arbitrary complexity on
Cartesian grids. We present an efficient algorithm to numerically implement the vorticity-velocity-pressure
formulation including a penalty term to simulate the pressure fields around a solid body. In vorticity-based
methods, pressure field can be independently computed from the solution procedure for vorticity. This
clearly simplifies the implementation and reduces the computational cost. Obtaining the pressure field at any
fixed time represents the most challenging goal of this study. We validate the implementation by numerical
simulations of an incompressible viscous flow around an impulsively started circular cylinder in a wide
range of Reynolds numbers: Re=40, 550, 3000, and 9500.


---

## Body

## 1. Introduction

Vortex methods are essentially a Lagrangian, grid-free approach in which fluid particles are used as basic computational elements to solve the Navier-Stokes equations in the vorticity-velocity formulation. The Navier-Stokes equations can be rewritten in terms of the vorticity. Taking the curl of the Navier-Stokes equation gives the vorticity transport equation (VTE), simplified as(1)∂ω∂t+u·∇ω=ω·∇u+ν∇2ω,
where the vorticity ω is defined as ω=∇×u and the kinematic viscosity ν is assumed to be constant. This equation governs the evolution of vorticity in an incompressible viscous flow. The velocity u is coupled with ω at all times in a self-consistent way through the relation of ∇·u=0 and ∇×u=ω. The position of the vorticity is determined from the velocity field and their strengths interact with each other. A fluid particle’s velocity can be evaluated through the Biot-Savart law. It is obvious that if the Biot-Savart law were used to calculate the velocity of each of the N particles in a simulation, then the direct calculation would involve O(N2) evaluations in a single time step. By contrast, Vortex-In-Cell (VIC) methods [1, 2] are used to determine the velocity field and place it onto a temporary grid, instead of direct integration using the Biot-Savart law, which is clearly inappropriate for large values of N. The VIC method has been improved with an immersed boundary method to handle more complex geometries (see [3–6]). The VIC method is a highly efficient hybrid particle-mesh algorithm to simulate an incompressible viscous flow past a solid body in an infinite domain [7]. Immersed boundary approaches, as reviewed in [8], are appropriate in free numerical computations of flows around complex geometries from technically difficult and time-consuming grid generation algorithms. Immersed boundary methods are used to enforce a no-slip condition at a body surface with decoupling between body boundaries and computational grids [3].Penalization methods are closely related to the immersed boundary approaches. Since penalty models tend to cancel vorticities inside a solid body, boundary conditions for both normal and tangential velocity components can be simultaneously imposed [9–11]. In penalization methods, boundary conditions are imposed by adding a penalization term to the momentum equations. There are two kinds of penalization modeling: L2 and H1 penalization. The L2 penalization consists of adding a damping term on the velocity in the momentum equation, whereas the H1 penalization additionally includes a perturbation of both the time derivative and a viscous term [12]. Angot et al. [12] have rigorously shown that the penalized equations can be used with confidence since penalty error is always negligible in the face of the error of approximation. They also pointed out that there is a small discrepancy induced by discretization of the penalized viscous term for the H1 penalization. Hence, most of the studies have practically adopted the L2 penalization. In the L2 penalty model, a solid body is regarded as a porous medium, where the permeability is infinite in the fluid part and tends to zero in the solid part [11, 13, 14]. The damping term has permeability η and a mask function χ, which describes the shape of the solid walls. The main advantage of penalization methods, which allows both normal and tangential velocity components at a solid body surface to be cancelled at once, is its simplicity. Penalization methods and VIC methods have been separately developed to simulate incompressible viscous flows around obstacles, and their combination has been achieved in recent years [7, 10, 15, 16].The vorticity-velocity formulation of the Navier-Stokes equations allows a purely kinematical problem to be decoupled from the pressure term. Clearly, the pressure term is eliminated by applying the curl operator. It is possible to obtain vorticity and velocity fields prior to any pressure calculation, and then the pressure field is explicitly treated by solving a suitable Poisson equation. In our previous studies [23, 24], we proposed an integral approach based on the vorticity-velocity-pressure formulation for obtaining the pressure field at any fixed time. A mathematical identity for a vector or scalar field is used to define field values of a quantity of interest, which involves an integral of singularities distributed over the surface and throughout the field. This integral approach has been successfully established in three different numerical schemes (refer to [4, 24]); an Eulerian finite volume method, a Lagrangian vortex method, and a hybrid Eulerian-Lagrangian method (VIC method). However, the boundary integral approach has disadvantages such as the presence of singular Green’s kernels. Special attention is needed to the accurate calculation of boundary integrals around a singular point. The biggest disadvantage in practice is the higher mathematical complexity needed to get a usable computational formulation. The matrices that result from the integral method are asymmetrical, and they are not easy to solve.In this study, we propose a simple way to numerically implement the vorticity-velocity-pressure formulation including a penalty term. To validate the numerical model, flows around an impulsively started cylinder are simulated for Reynolds numbers of 40, 550, 3000, and 9500. Most importantly, we attempt to solve the pressure field to extend the application of the penalization method. Governing equations for incompressible viscous flows are described below in Section2. Section 2.1 presents the conventional VIC method, and Section 2.2 describes the penalization method. Section 3 describes the vorticity-velocity-pressure formulation including a penalty term to simulate pressure fields of two-dimensional incompressible viscous flows around a solid body, which represents the most challenging goal of this study. This study is limited to two dimensions but can easily be extended to three dimensions. Flow simulations external to two-dimensional bodies are presented in Section 4.


## 2. Governing Equations

Two-dimensional incompressible unsteady flow of a viscous fluid can be rewritten in a domainΩ with boundary ∂Ω as follows:
(2)∂ω∂t+u·∇ω=ν∇2ωinΩ,
where ω is the scalar plane component of the vorticity vector (ω=ωe^z). The evolution of a flow is considered in discrete time steps. In each time step, vorticity is convected (according to u·∇ω) and diffused (according to ν∇2ω). The algorithm of viscous splitting consists of substeps in which the convective and diffusive effects are considered. Equivalently, the viscous splitting algorithm is expressed in a Lagrangian frame as (3a)DxDt=u,(3b)DωDt=ν∇2ω, where D/Dt=∂/∂t+u·∇ is the material derivative. N discrete Lagrangian fluid elements are linearly superposed to approximate the vorticity field as
(4)ω(x,t)=∑pN‍Γpζp(x-xp),
where ζi is a mollification of the Dirac-Delta function. Each particle is characterized by its position xp and its strength Γp, that is, a circulation Γp=∫‍ωdS≈ωpSp where Sp is the area of fluid associated with the pth particle, p=1,…,N. The vorticity-carrying fluid element at a point x=(xi,yi) is advanced with the local flow velocity in the convection step (see (3a)), and diffusion acts at these new locations to modify the vorticity field of a flow (see (3b)).


### 2.1. Vortex-In-Cell (VIC) Method

The convective part is solved through the VIC method. The velocity vectoru is expressed as
(5)u=U∞+uω=U∞+∇×ψ,
where U∞ is the free stream velocity, uω represents the rotational velocity, and ψ is a stream function. At infinity, the velocity field recovers the free stream velocity.  Taking the curl of (5) yields the Poisson equation
(6)∇2ψ=-ω,
linking vorticity to the stream function and in turn to the velocity. To solve (6), vorticity particles are remeshed onto an equally spaced Cartesian grid through
(7)ωm=1h2∑pN‍ΓpW(xp-xmh),
where h is the mesh spacing, using the following third order interpolation kernel:
(8)W(x)=M4′(x)={0for|x|>212(2-|x|)2(1-|x|)for1<|x|≤21-52|x|2+32|x|3for|x|≤1.
The stream function can be computed on a temporary grid from the Poisson equation with a homogeneous or nonhomogeneous Dirichlet boundary condition [4], using a fast Poisson solver based on the FFT (fast Fourier transform). Grids to solve the stream function can typically define a compact computational domain with nonhomogeneous Dirichlet boundary conditions which are imposed everywhere on the boundary. Unknown stream functions at domain boundaries are directly obtained from the Green’s function solution. It follows that
(9)ψ(xm)e^z=-12π∑p=1N‍Γpln⁡(|xp-xm|),
where xp∈Ω, xm∈∂Ω, and xm≠xp. The velocity field on a grid is computed from the definition uω=∇×ψ through the fourth order centered difference scheme.


### 2.2.  Penalization Method

The penalization method was initially designed to take into account solid obstacles in fluid flows. The main benefit of the penalty term is to replace the usual vorticity creation algorithm for enforcing a no-slip condition at a solid body. By adding the penalization term to the Navier-Stokes equation, the vorticity transport equation becomes(10)DωDt=ν∇2ω+∇×[λχ(us-u)],
where us is the velocity of the solid body and χ denotes a mask function that yields 0 in the fluid and 1 in the solid [7, 12]. λ indicates the penalization parameter with dimensions [s-1], equivalent to an inverse permeability η. In the case of explicit Euler time discretization, the penalization parameter must satisfy λΔt<O(1) to ensure stability. It is important to note that λ is an arbitrary parameter, independent of the spatial or temporal discretization, and thus the boundary conditions can be enforced to any desired accuracy by choosing λ appropriately [12, 25]. The parameter λ may be as large as necessary to accurately represent the no-slip boundary condition. The larger λ means that this term is stiff and must be solved implicitly, since the time step Δt should not be too small for long time-scale simulations of unsteady flows [7]. The penalization term can be evaluated implicitly by the following expression:
(11)ΔωΔt=∇×[λχ(usn-u~n+1)],
where
(12)u~n+1=un+λΔtχusn1+λΔtχ.
This scheme is unconditionally stable [7, 10]. Finally, (10) is rewritten as
(13)ωn+1-ωnΔt=ν∇2ωn+∇×[λχ1+λΔtχ(usn-un)],
by replacing (11) and (12). The diffusion and penalization terms in the penalized vorticity transport equation (see (13)) are solved separately. The diffusion term is evaluated onto a grid with a classical 9-point finite difference method. To solve the contribution of the penalization term, the derivative is computed through a centered difference approximation with fourth order error. Vorticity values at all grid points are updated from Δω.


## 2.1. Vortex-In-Cell (VIC) Method

The convective part is solved through the VIC method. The velocity vectoru is expressed as
(5)u=U∞+uω=U∞+∇×ψ,
where U∞ is the free stream velocity, uω represents the rotational velocity, and ψ is a stream function. At infinity, the velocity field recovers the free stream velocity.  Taking the curl of (5) yields the Poisson equation
(6)∇2ψ=-ω,
linking vorticity to the stream function and in turn to the velocity. To solve (6), vorticity particles are remeshed onto an equally spaced Cartesian grid through
(7)ωm=1h2∑pN‍ΓpW(xp-xmh),
where h is the mesh spacing, using the following third order interpolation kernel:
(8)W(x)=M4′(x)={0for|x|>212(2-|x|)2(1-|x|)for1<|x|≤21-52|x|2+32|x|3for|x|≤1.
The stream function can be computed on a temporary grid from the Poisson equation with a homogeneous or nonhomogeneous Dirichlet boundary condition [4], using a fast Poisson solver based on the FFT (fast Fourier transform). Grids to solve the stream function can typically define a compact computational domain with nonhomogeneous Dirichlet boundary conditions which are imposed everywhere on the boundary. Unknown stream functions at domain boundaries are directly obtained from the Green’s function solution. It follows that
(9)ψ(xm)e^z=-12π∑p=1N‍Γpln⁡(|xp-xm|),
where xp∈Ω, xm∈∂Ω, and xm≠xp. The velocity field on a grid is computed from the definition uω=∇×ψ through the fourth order centered difference scheme.


## 2.2.  Penalization Method

The penalization method was initially designed to take into account solid obstacles in fluid flows. The main benefit of the penalty term is to replace the usual vorticity creation algorithm for enforcing a no-slip condition at a solid body. By adding the penalization term to the Navier-Stokes equation, the vorticity transport equation becomes(10)DωDt=ν∇2ω+∇×[λχ(us-u)],
where us is the velocity of the solid body and χ denotes a mask function that yields 0 in the fluid and 1 in the solid [7, 12]. λ indicates the penalization parameter with dimensions [s-1], equivalent to an inverse permeability η. In the case of explicit Euler time discretization, the penalization parameter must satisfy λΔt<O(1) to ensure stability. It is important to note that λ is an arbitrary parameter, independent of the spatial or temporal discretization, and thus the boundary conditions can be enforced to any desired accuracy by choosing λ appropriately [12, 25]. The parameter λ may be as large as necessary to accurately represent the no-slip boundary condition. The larger λ means that this term is stiff and must be solved implicitly, since the time step Δt should not be too small for long time-scale simulations of unsteady flows [7]. The penalization term can be evaluated implicitly by the following expression:
(11)ΔωΔt=∇×[λχ(usn-u~n+1)],
where
(12)u~n+1=un+λΔtχusn1+λΔtχ.
This scheme is unconditionally stable [7, 10]. Finally, (10) is rewritten as
(13)ωn+1-ωnΔt=ν∇2ωn+∇×[λχ1+λΔtχ(usn-un)],
by replacing (11) and (12). The diffusion and penalization terms in the penalized vorticity transport equation (see (13)) are solved separately. The diffusion term is evaluated onto a grid with a classical 9-point finite difference method. To solve the contribution of the penalization term, the derivative is computed through a centered difference approximation with fourth order error. Vorticity values at all grid points are updated from Δω.


## 3. Pressure Poisson Equation

In vorticity-based methods, the pressure field is independently computed from the solution procedure for vorticity. The divergence operation on the penalized Navier-Stokes equation leads to the Poisson equation for pressure:(14)∇2H=∇·(u×ω)+∇·[λχ(us-u)],
where H is the total pressure including the static and dynamic pressure, which can be defined as
(15)H=p-P∞ρ+12(|u|2-|U∞|2),
where P∞ and U∞ represent the pressure and velocity of an undisturbed flow, respectively. The right hand side in (14) is approximated with the velocity and vorticity values on a temporary grid. The total pressure can be computed on a temporary grid using an FFT-based Poisson solver. In this case, pseudovelocities in a solid body become important because they affect pressure distributions through velocity divergence across solid boundaries. The pressure Poisson equation modified by the penalty term is similar to the conventional velocity-pressure formulation with an additional divergence damping term, α∇·u, to suppress the spurious divergence. This technique of adding a damping term (see [26]) or an appropriate stabilizing term by altering the pressure boundary condition (see [9]) is well known and has been used previously by a number of researchers in the field of incompressible flows. In penalization methods, however, the flow velocity field u in a whole domain Ω satisfies the incompressible mass conservation equation, ∇·u=0 [12].For simplicity, let us consider an incompressible viscous flow past a fixed solid obstacle. The pressure Poisson equation is thus rewritten as(16)∇2H=∇·(u×ω)-λ′(u·∇χ)=∇·(u×ω)-∇·(λ′χu).
The penalization parameter in (16) can be expressed as λ′ to avoid confusion with the penalization parameter λ in (13). Here, we address the issue of how large λ′ should be for numerical purposes. It is thought that no precision is needed for the selection of the parameter λ′, just an order of magnitude. Recall that in an implicit scheme for the vorticity transport equation, the penalization parameter (see (13)) actually becomes λχ/(1+λΔtχ). Increasing λ leads to λχ/(1+λΔtχ)≈1/Δt, equivalent to the parameter of an explicit scheme λ=1/Δt. Hence, the penalization parameter λ′ in the pressure Poisson equation is chosen to be λ′=1/Δt in this study. For a well-resolved simulation, vortex methods are not constrained with the usual CFL condition. Instead, it is natural to use the Reynolds number based on vorticity, Reh=|ω|h2/ν where h is a fixed grid spacing, and the relation Reh=O(1) should be fulfilled to ensure numerical stability [27]. The time step size Δt is also required to be νΔt/h2=O(1) for the diffusion. It is basically found that |ω|Δt=O(1) [7, 27]. This means that the penalization parameter λ′ is equivalent to λ′=1/Δt=O(|ω|). Finally, the two terms in the right hand side of (14) have the same order of magnitude inside a body. Figure 1 shows pressure fields calculated by using different orders of magnitude for the penalization parameter; consider λ′=0, 0.1/Δt, 1/Δt, and 10/Δt. It is thought that λ′=1/Δt is a good choice. This is numerically investigated in more detail in Section 4.Instantaneous pressure contours at different pressure penalty parameters: (a)λ′=0, (b) λ′=0.1/Δt, (c) λ′=1/Δt, and (d) λ′=10/Δt. Note that a flow over an impulsively started circular cylinder for Re=550 at T=3 is used to compute pressure fields.
(a)
(b)
(c)
(d)Dirichlet boundary conditions are also required to solve the pressure Poisson equations (see (16)). It is obvious that the total pressure H becomes zero at infinity. Although the total pressure H is directly undetermined at boundaries of a computational domain, this does not necessitate choosing excessively large computational domains. If domain boundaries are far enough from a solid body, the homogeneous Dirichlet boundary conditions (H=0) may be used. The computational domain size will be further discussed in Section 4.5.


## 4. Numerical Simulations

We begin with a series of flow simulations around an impulsively started circular cylinder. The impulsively started cylinder problem has been studied as a good prototype of unsteady separated flows and is considered to validate our numerical method. Initially, the flow is at rest and then impulsively started from rest at timet=0. In other words, an impulsive start of the cylinder is simulated by specifying uniform inflow velocities (u=U∞) in the whole computational domain. Boundary conditions at the body surface are imposed by the penalization terms.In this study, the Reynolds number,Re=U∞D/ν, based on the free stream velocity, U∞, and the diameter of the cylinder, D, is selected to be 40, 550, 3000, or 9500. We only focus on the flow in the early stage of development, which is known to remain symmetric for Re up to 9500. In our simulations, no symmetry constraint is imposed. The time (t) is nondimensionalized based on the radius (R=0.5D) of the cylinder as T=Ut/R. Vorticity particles are convected by the local flow velocity for each time increment. Vorticity values on the grid are transferred back to the particles before the positions of fluid particles are advanced in time using a second order Runge-Kutta method. Both Poisson equations for stream-function and pressure are solved through the FFT-based Poisson solver based on an open-source library called FFTW (fastest Fourier transform in the west).The initial computational domain is chosen asΩ=[-1.50D,3.62D]×[-2.56D,2.56D] for Re=40, 550, and 3000 and Ω=[-1.5D,2.6D]×[-2.05D,2.05D] for Re=9500. The penalization parameter λ for the diffusion is fixed to 108.


### 4.1. Steady Flow Past a Circular Cylinder at Re = 40

The numerical parameters areh/D=0.04 and Δt=0.016, which are determined through the stability condition Δt=h2/4ν [7]. At Re=40, it is well known that the flow reaches a steady state. A pair of stationary recirculating wakes develops behind the cylinder. The wake length L/D and separation angles θs is 4.16 and 51.7∘, respectively, in our simulation. Fornberg (1980) [28] made a similar remark and gave L/D=4.48 and θs=51.5∘  experimentally. The drag coefficient which is defined as CD=2Fd/(ρU∞2D) is computed as CD=1.483. This is close to the experimentally measured value of  1.498 in [28].Figure2(a) shows the computed pressure coefficient (Cp) distribution for Re=40. Here, it is defined as Cp=2(p-P∞)/ρU∞2. As Angot et al. [12] pointed out, the pressure is continuous through the cylinder body. There is no Gibb’s oscillation associated with the discontinuity at the body surface. In Figure 2(b), the pressure coefficient along the body surface is plotted together with the experimental data of Grove et al. [29]. The pressure coefficients along the body surface were obtained through the kriging interpolation algorithm in Tecplot. Comparing the pressure coefficient at the front and rear stagnation points with previous studies listed in Table 1, it is obvious that good agreement is achieved between our results and those of the previous studies.Table 1
Comparison of pressure coefficient at the front and rear stagnation points with previous studies.
Author(s)
C
p (Front)
-
C
p (Rear)
Grove et al. (1964) [29], experimental
1.196
0.509
Dennis and Chang (1970) [30], computational
1.133
0.487
Fornberg (1980) [28], experimental
1.14
0.46
He and Doolen (1997) [31], computational
1.133
0.487
Pontaza and Reddy (2003) [32], computational
1.197
0.484
Present
1.169
0.537Pressure contours (right) and pressure coefficients at the body surface (left) for steady flow past a circular cylinder atRe=40.
(a)
C
p distribution
(b)
C
p along the body surfaceIn vortex methods, redistribution is typically done every five time steps to maintain spatial uniformity of the particle distribution. If a new particle has a small vorticity magnitude, that is,|Γ|<0.0001|Γ|max⁡, it is deleted to avoid having too many particles. In the case of Re=40, however, it is noted that a solid body should be covered with fluid particles to obtain a stable solution. This is very critical at lower Reynolds numbers. The presence of particles around the cylinder body was therefore enforced after each redistribution. Actually, any fluid particle within the distance of 5h from the body surface is not removed.


### 4.2.  Unsteady Flow Past a Circular Cylinder at Re = 550

The simulation parameters areΔt=0.002, h/D=0.005, and λ=108. The time step is determined by the condition Δt<h2/4ν≈3.438×10-3. Simulation was carried out until T=10 to validate the present formulation in the early time stage after the impulsive start. The number of fluid particles ranged from approximately 30,000 to 150,000, and the total computation time is approximately 4 hours on 8 Intel Xeon64 3.3 GHz nodes. The code was implemented by adopting the Open-MPI (Message Passing Interface) library in C. A spatial computational domain for parallel computing is decomposed by equally splitting it into several subdomains depending on the number of processors used. Each processor is responsible for one subdomain in real space [33].A pair of secondary symmetric vortices appears and becomes stronger as time increases. The so-called bulge phenomena observed experimentally by Bouard and Coutanceau (1980) is well captured by our numerical simulation. Figure3 shows that the streamline pattern computed for Re=550 at T=5 compares very well with the flow visualization result of Bouard and Coutanceau (1980) [17]. The vortex core position indicated by the coordinates a and b is investigated. The nondimensionalized abscissa a/D and ordinate b/D are 0.34 and 0.27, respectively, and the wake length L/D is 0.82 in our simulation. Bouard and Coutanceau made a similar remark and gave a/D=0.36, b/D=0.28, and L/D=0.85. Figure 4 shows the time evolution of the drag coefficient for an impulsively started flow around a two-dimensional cylinder for Re=550 calculated with the present method, compared with results from the short time asymptotic solution of Bar-Lev and Yang (1975) [18], the vortex method result of Koumoutsakos and Leonard (1995) [19], and the vortex-in-cell method result of Kudela and Kozlowski (2009) [20]. The good agreements demonstrate that our penalized VIC method gives reliable force results and that the flow for Re=550 is very well simulated.Comparison of the streamline patterns forRe=550 at T=5.0: (a) flow visualization by Bouard and Coutancean (1980) [17] (left) and (b) present method (right).
(a)
(b)Figure 4
Time evolution of the drag coefficient for the impulsively started cylinder atRe=550 compared with the short-time theoretical result of Bar-Lev and Yang [18], the vortex method result of Koumoutsakos and Leonard [19], and the Vortex-In-Cell method result of Kudela and Kozlowski [20].Figure5 shows time evolutions of the pressure distribution. A pair of lower pressure regions moving downstream results from the strong vortical flow and the vorticity strength becomes weaker by viscous diffusion as time increases. The pressure distribution at the cylinder surface is shown in Figure 6. At the initial time in the numerical simulation, the surface pressure distribution is quite close to that obtained from the potential flow analysis. The pressure coefficient at the front stagnation point of the cylinder is almost constant, its value being 1.0. The pressure at the rear of the cylinder is more sensitive than at the front. The minimum pressure coefficient gradually moves upstream with increasing time. Meanwhile, the corresponding locations of Cp=0 have the same upstream moving trends. As experimentally measured by Norberg (2002) [34], Cp=0 in the present simulation is located at the angular position θ≈35∘.Pressure coefficient distribution around a circular cylinder forRe=550: (a) T=1, (b) T=3, (c) T=5, and (d) T=7.
(a)
(b)
(c)
(d)Figure 6
Surface pressure distribution at several instants forRe=550.


### 4.3.  Unsteady Flow Past a Circular Cylinder at Re = 3000

The simulation parameters areΔt=0.002, h/D=0.0025, and λ=108. Simulation was carried out until T=10 and the number of fluid particles ranged from approximately 130,000 to 420,000. The total run time is about 9 hours on 16 Intel Xeon64 3.3 GHz nodes. Compared with flow at Re=550, the secondary vortices appear at an earlier time and grow larger. In the case of Re=3000, the two secondary vortices formed are equivalent in size and in strength. It is the so-called α phenomena (see Bouard and Coutanceau (1980) [17]). The present numerical simulation is able to correctly capture the expected physics of the flow for Re=3000. Figure 7 shows that the streamline pattern computed for Re=3000 compares excellently with the flow visualization result of Loc and Bouard (1985) [21].Instantaneous streamlines for an impulsively started circular cylinder atRe=3000: (a) T=1, (b) T=3, and (c) T=5, compared with flow visualization results of Loc and Bouard (1985) [21].
(a)
(b)
(c)Pressure fields are explicitly computed through vorticity and velocity fields of well-simulated flows forRe=3000. Figure 8 shows time evolutions of the pressure distribution around the cylinder body. As for the case of Re=550, a pair of lower pressure regions moving downstream results from the strong vortical flow as time increases. In Figure 9, the surface pressure distributions at three different times are compared with the numerical results of Chang and Chern (1991) [22]. The present results are qualitatively similar to theirs, although an offset of similar order of magnitude remains. Therefore, it appears that the flow is developed slightly faster in our simulation. We conjecture that this offset may be due to the difference in primary variables for the numerical simulation. The pressure at the front and rear stagnation points of the cylinder agrees reasonably well with their results, and the corresponding locations of Cp=0 are also matched very well. These agreements help to demonstrate the quality of the numerical results obtained by the present algorithm.Pressure coefficient distribution around a circular cylinder forRe=3000: (a) T=1, (b) T=3, (c) T=5, and (d) T=7.
(a)
(b)
(c)
(d)Figure 9
Surface pressure distribution at several instants forRe=3000 compared with numerical results of Chang and Chern (1991) [22].


### 4.4.  Unsteady Flow Past a Circular Cylinder at Re = 9500

The simulation parameters areΔt=0.002, h/D=0.001, and λ=108. Flows were numerically simulated until T=4, and the number of fluid particles ranged from approximately 800,000 to 1,200,000. The total run time is about 25 hours on 24 Intel Xeon64 3.3 GHz nodes. Numerical results are presented for vorticity isocontours in Figure 10. The vorticity isocontour patterns are in good agreement with the vortex immersed boundary method results of Cottet et al. (2010) [35], although their results are not presented here. In Figure 11(a), it is seen that the so-called α and β phenomena observed experimentally by Bouard and Coutanceau (1980) are well captured by our numerical simulation. The streamline patterns computed for Re=9500 compare quite well with the flow visualization results of Loc and Bouard (1985) [21]. These agreements mean that vorticity and velocity fields computed using our penalized VIC method correctly capture the expected physics of the flow at Re=9500.Instantaneous vorticity isocontours for an impulsively started circular cylinder atRe=9500: (a) T=1, (b) T=2, (c) T=3, and (d) T=4.
(a)
(b)
(c)
(d)Instantaneous streamlines for an impulsively started circular cylinder atRe=9500: (a) T=2 and (b) T=4 compared with flow visualization results of Loc and Bouard (1985) [21].
(a)
(b)Figure12 shows the pressure distribution around the cylinder body, which is computed using the well-simulated flow for Re=9500. The pressure coefficients along the body surface are extracted from the pressure isocontour. In Figure 13, the coefficients at three different times are compared with the numerical results of Suh and Kim (1999) [23] and remarkable agreement is obtained.Pressure coefficient distribution around a circular cylinder forRe=9500: (a) T=1, (b) T=2, (c) T=3, and (d) T=4.
(a)
(b)
(c)
(d)Figure 13
Surface pressure distribution at several instants forRe=9500 compared with numerical results of Suh and Kim (1999) [23].


### 4.5.  Effect of Domain Size

The numerical parameters studied here are the outer boundary location, grid spacing,h, and time step value, Δt. The two first parameters are used for the discretization of the Poisson equations. The last one is the time discretization parameter and is closely related to the determination of the penalization parameter for the pressure Poisson equation. We studied the influence of these numerical parameters on the solution of the pressure Poisson equation for Re=550 as an example.To demonstrate that excessively large computational domains are not required, we carried out numerical simulations for flows with different outer boundaries. The grid spacing is fixed to0.005 and the time step value is 0.002. We tested the FFT points of 1024×1024, 2048×2048, and 4096×4096. Actual computational domains correspond to Ω=[-1.5D,3.62D]×[-2.56D,2.56D], [-4.0D,6.24D]×[-5.12D,5.12D], and [-9.0D,11.48D]×[-10.24D,10.24D], respectively. Figure 14 confirms the convergence of numerical results for Re=550. The computational domain Ω=[-1.50D,3.62D]×[-2.56D,2.56D], chosen at an early stage of flow development, is large enough to obtain a satisfactory solution in terms of pressure distributions near the body as shown in Figure 14(b). Enlarging the outer boundary barely has any effect on the results of vorticity as shown in Figure 15.Effect of the outer boundary forRe=550 at T=3.
(a)
Comparison of pressure contours at different domain sizes
(b)
Pressure coefficients on the body surfaceComparison of vorticity isocontours on different computational domain sizes forRe=550: FFT points of 1024×1024 (solid line) versus4096×4096 (dashed line).
(a)
T
=
3
(b)
T
=
7


## 4.1. Steady Flow Past a Circular Cylinder at Re = 40

The numerical parameters areh/D=0.04 and Δt=0.016, which are determined through the stability condition Δt=h2/4ν [7]. At Re=40, it is well known that the flow reaches a steady state. A pair of stationary recirculating wakes develops behind the cylinder. The wake length L/D and separation angles θs is 4.16 and 51.7∘, respectively, in our simulation. Fornberg (1980) [28] made a similar remark and gave L/D=4.48 and θs=51.5∘  experimentally. The drag coefficient which is defined as CD=2Fd/(ρU∞2D) is computed as CD=1.483. This is close to the experimentally measured value of  1.498 in [28].Figure2(a) shows the computed pressure coefficient (Cp) distribution for Re=40. Here, it is defined as Cp=2(p-P∞)/ρU∞2. As Angot et al. [12] pointed out, the pressure is continuous through the cylinder body. There is no Gibb’s oscillation associated with the discontinuity at the body surface. In Figure 2(b), the pressure coefficient along the body surface is plotted together with the experimental data of Grove et al. [29]. The pressure coefficients along the body surface were obtained through the kriging interpolation algorithm in Tecplot. Comparing the pressure coefficient at the front and rear stagnation points with previous studies listed in Table 1, it is obvious that good agreement is achieved between our results and those of the previous studies.Table 1
Comparison of pressure coefficient at the front and rear stagnation points with previous studies.
Author(s)
C
p (Front)
-
C
p (Rear)
Grove et al. (1964) [29], experimental
1.196
0.509
Dennis and Chang (1970) [30], computational
1.133
0.487
Fornberg (1980) [28], experimental
1.14
0.46
He and Doolen (1997) [31], computational
1.133
0.487
Pontaza and Reddy (2003) [32], computational
1.197
0.484
Present
1.169
0.537Pressure contours (right) and pressure coefficients at the body surface (left) for steady flow past a circular cylinder atRe=40.
(a)
C
p distribution
(b)
C
p along the body surfaceIn vortex methods, redistribution is typically done every five time steps to maintain spatial uniformity of the particle distribution. If a new particle has a small vorticity magnitude, that is,|Γ|<0.0001|Γ|max⁡, it is deleted to avoid having too many particles. In the case of Re=40, however, it is noted that a solid body should be covered with fluid particles to obtain a stable solution. This is very critical at lower Reynolds numbers. The presence of particles around the cylinder body was therefore enforced after each redistribution. Actually, any fluid particle within the distance of 5h from the body surface is not removed.


## 4.2.  Unsteady Flow Past a Circular Cylinder at Re = 550

The simulation parameters areΔt=0.002, h/D=0.005, and λ=108. The time step is determined by the condition Δt<h2/4ν≈3.438×10-3. Simulation was carried out until T=10 to validate the present formulation in the early time stage after the impulsive start. The number of fluid particles ranged from approximately 30,000 to 150,000, and the total computation time is approximately 4 hours on 8 Intel Xeon64 3.3 GHz nodes. The code was implemented by adopting the Open-MPI (Message Passing Interface) library in C. A spatial computational domain for parallel computing is decomposed by equally splitting it into several subdomains depending on the number of processors used. Each processor is responsible for one subdomain in real space [33].A pair of secondary symmetric vortices appears and becomes stronger as time increases. The so-called bulge phenomena observed experimentally by Bouard and Coutanceau (1980) is well captured by our numerical simulation. Figure3 shows that the streamline pattern computed for Re=550 at T=5 compares very well with the flow visualization result of Bouard and Coutanceau (1980) [17]. The vortex core position indicated by the coordinates a and b is investigated. The nondimensionalized abscissa a/D and ordinate b/D are 0.34 and 0.27, respectively, and the wake length L/D is 0.82 in our simulation. Bouard and Coutanceau made a similar remark and gave a/D=0.36, b/D=0.28, and L/D=0.85. Figure 4 shows the time evolution of the drag coefficient for an impulsively started flow around a two-dimensional cylinder for Re=550 calculated with the present method, compared with results from the short time asymptotic solution of Bar-Lev and Yang (1975) [18], the vortex method result of Koumoutsakos and Leonard (1995) [19], and the vortex-in-cell method result of Kudela and Kozlowski (2009) [20]. The good agreements demonstrate that our penalized VIC method gives reliable force results and that the flow for Re=550 is very well simulated.Comparison of the streamline patterns forRe=550 at T=5.0: (a) flow visualization by Bouard and Coutancean (1980) [17] (left) and (b) present method (right).
(a)
(b)Figure 4
Time evolution of the drag coefficient for the impulsively started cylinder atRe=550 compared with the short-time theoretical result of Bar-Lev and Yang [18], the vortex method result of Koumoutsakos and Leonard [19], and the Vortex-In-Cell method result of Kudela and Kozlowski [20].Figure5 shows time evolutions of the pressure distribution. A pair of lower pressure regions moving downstream results from the strong vortical flow and the vorticity strength becomes weaker by viscous diffusion as time increases. The pressure distribution at the cylinder surface is shown in Figure 6. At the initial time in the numerical simulation, the surface pressure distribution is quite close to that obtained from the potential flow analysis. The pressure coefficient at the front stagnation point of the cylinder is almost constant, its value being 1.0. The pressure at the rear of the cylinder is more sensitive than at the front. The minimum pressure coefficient gradually moves upstream with increasing time. Meanwhile, the corresponding locations of Cp=0 have the same upstream moving trends. As experimentally measured by Norberg (2002) [34], Cp=0 in the present simulation is located at the angular position θ≈35∘.Pressure coefficient distribution around a circular cylinder forRe=550: (a) T=1, (b) T=3, (c) T=5, and (d) T=7.
(a)
(b)
(c)
(d)Figure 6
Surface pressure distribution at several instants forRe=550.


## 4.3.  Unsteady Flow Past a Circular Cylinder at Re = 3000

The simulation parameters areΔt=0.002, h/D=0.0025, and λ=108. Simulation was carried out until T=10 and the number of fluid particles ranged from approximately 130,000 to 420,000. The total run time is about 9 hours on 16 Intel Xeon64 3.3 GHz nodes. Compared with flow at Re=550, the secondary vortices appear at an earlier time and grow larger. In the case of Re=3000, the two secondary vortices formed are equivalent in size and in strength. It is the so-called α phenomena (see Bouard and Coutanceau (1980) [17]). The present numerical simulation is able to correctly capture the expected physics of the flow for Re=3000. Figure 7 shows that the streamline pattern computed for Re=3000 compares excellently with the flow visualization result of Loc and Bouard (1985) [21].Instantaneous streamlines for an impulsively started circular cylinder atRe=3000: (a) T=1, (b) T=3, and (c) T=5, compared with flow visualization results of Loc and Bouard (1985) [21].
(a)
(b)
(c)Pressure fields are explicitly computed through vorticity and velocity fields of well-simulated flows forRe=3000. Figure 8 shows time evolutions of the pressure distribution around the cylinder body. As for the case of Re=550, a pair of lower pressure regions moving downstream results from the strong vortical flow as time increases. In Figure 9, the surface pressure distributions at three different times are compared with the numerical results of Chang and Chern (1991) [22]. The present results are qualitatively similar to theirs, although an offset of similar order of magnitude remains. Therefore, it appears that the flow is developed slightly faster in our simulation. We conjecture that this offset may be due to the difference in primary variables for the numerical simulation. The pressure at the front and rear stagnation points of the cylinder agrees reasonably well with their results, and the corresponding locations of Cp=0 are also matched very well. These agreements help to demonstrate the quality of the numerical results obtained by the present algorithm.Pressure coefficient distribution around a circular cylinder forRe=3000: (a) T=1, (b) T=3, (c) T=5, and (d) T=7.
(a)
(b)
(c)
(d)Figure 9
Surface pressure distribution at several instants forRe=3000 compared with numerical results of Chang and Chern (1991) [22].


## 4.4.  Unsteady Flow Past a Circular Cylinder at Re = 9500

The simulation parameters areΔt=0.002, h/D=0.001, and λ=108. Flows were numerically simulated until T=4, and the number of fluid particles ranged from approximately 800,000 to 1,200,000. The total run time is about 25 hours on 24 Intel Xeon64 3.3 GHz nodes. Numerical results are presented for vorticity isocontours in Figure 10. The vorticity isocontour patterns are in good agreement with the vortex immersed boundary method results of Cottet et al. (2010) [35], although their results are not presented here. In Figure 11(a), it is seen that the so-called α and β phenomena observed experimentally by Bouard and Coutanceau (1980) are well captured by our numerical simulation. The streamline patterns computed for Re=9500 compare quite well with the flow visualization results of Loc and Bouard (1985) [21]. These agreements mean that vorticity and velocity fields computed using our penalized VIC method correctly capture the expected physics of the flow at Re=9500.Instantaneous vorticity isocontours for an impulsively started circular cylinder atRe=9500: (a) T=1, (b) T=2, (c) T=3, and (d) T=4.
(a)
(b)
(c)
(d)Instantaneous streamlines for an impulsively started circular cylinder atRe=9500: (a) T=2 and (b) T=4 compared with flow visualization results of Loc and Bouard (1985) [21].
(a)
(b)Figure12 shows the pressure distribution around the cylinder body, which is computed using the well-simulated flow for Re=9500. The pressure coefficients along the body surface are extracted from the pressure isocontour. In Figure 13, the coefficients at three different times are compared with the numerical results of Suh and Kim (1999) [23] and remarkable agreement is obtained.Pressure coefficient distribution around a circular cylinder forRe=9500: (a) T=1, (b) T=2, (c) T=3, and (d) T=4.
(a)
(b)
(c)
(d)Figure 13
Surface pressure distribution at several instants forRe=9500 compared with numerical results of Suh and Kim (1999) [23].


## 4.5.  Effect of Domain Size

The numerical parameters studied here are the outer boundary location, grid spacing,h, and time step value, Δt. The two first parameters are used for the discretization of the Poisson equations. The last one is the time discretization parameter and is closely related to the determination of the penalization parameter for the pressure Poisson equation. We studied the influence of these numerical parameters on the solution of the pressure Poisson equation for Re=550 as an example.To demonstrate that excessively large computational domains are not required, we carried out numerical simulations for flows with different outer boundaries. The grid spacing is fixed to0.005 and the time step value is 0.002. We tested the FFT points of 1024×1024, 2048×2048, and 4096×4096. Actual computational domains correspond to Ω=[-1.5D,3.62D]×[-2.56D,2.56D], [-4.0D,6.24D]×[-5.12D,5.12D], and [-9.0D,11.48D]×[-10.24D,10.24D], respectively. Figure 14 confirms the convergence of numerical results for Re=550. The computational domain Ω=[-1.50D,3.62D]×[-2.56D,2.56D], chosen at an early stage of flow development, is large enough to obtain a satisfactory solution in terms of pressure distributions near the body as shown in Figure 14(b). Enlarging the outer boundary barely has any effect on the results of vorticity as shown in Figure 15.Effect of the outer boundary forRe=550 at T=3.
(a)
Comparison of pressure contours at different domain sizes
(b)
Pressure coefficients on the body surfaceComparison of vorticity isocontours on different computational domain sizes forRe=550: FFT points of 1024×1024 (solid line) versus4096×4096 (dashed line).
(a)
T
=
3
(b)
T
=
7


## 5. Conclusion

The hybrid VIC-penalization method was applied in this study to simulate flows around an impulsively started circular cylinder in a wide range of Reynolds numbers:Re=40, 550, 3000, and 9500. The penalization method was used to enforce solid boundary conditions. The main issue of this study is the simulation of pressure fields around a solid body though the vorticity-velocity-pressure formulation including a penalty term. The pressure fields were independently computed from the solution procedure for vorticity, and it maintains the decoupling between a solid body and uniform rectangular meshes. For validation, the results of steady state flow at low Reynolds numbers and the early stage of unsteady flow at moderate and high Reynolds numbers were obtained and compared well with previous experimental and computational data available in the literature. The flows demonstrate the ability of our numerical algorithm to obtain the pressure field on a temporary grid without solving the integral equation for the total pressure. It is highly efficient and numerical implementation is easy.



---

*Source: 708372-2014-05-25.xml*