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

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.


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 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 selfconsistent 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 particles in a simulation, then the direct calculation would involve O( 2 ) 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 . The VIC method has been improved with an immersed boundary method to handle more complex geometries (see [3][4][5][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][10][11]. In penalization methods, boundary conditions are imposed by adding a penalization term to the momentum equations. There are two kinds of penalization modeling: 2 and 1 penalization. The 2 penalization consists of adding a damping term on the velocity in the momentum equation, whereas the 1 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 1 penalization. Hence, most of the studies have practically adopted the 2 penalization. In the 2 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 Section 2. 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.

Governing Equations
Two-dimensional incompressible unsteady flow of a viscous fluid can be rewritten in a domain Ω with boundary Ω as follows: where is the scalar plane component of the vorticity vector ( =̂). 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 where / = / + u ⋅ ∇ is the material derivative. discrete Lagrangian fluid elements are linearly superposed to approximate the vorticity field as where is a mollification of the Dirac-Delta function. Each particle is characterized by its position x and its strength Γ , that is, a circulation Γ = ∫ ≈ where is the area of fluid associated with the th particle, = 1, . . . , . The vorticity-carrying fluid element at a point x = ( , ) 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)).

Vortex-In-Cell (VIC) Method.
The convective part is solved through the VIC method. The velocity vector u is expressed as   where ∞ 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 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  [18], the vortex method result of Koumoutsakos and Leonard [19], and the Vortex-In-Cell method result of Kudela and Kozlowski [20].
where ℎ is the mesh spacing, using the following third order interpolation kernel: 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 where x ∈ Ω, x ∈ Ω, and x ̸ = x . The velocity field on a grid is computed from the definition u = ∇ × through the fourth order centered difference scheme.

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 where u 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 [ −1 ], equivalent to an inverse permeability . In the case of explicit Euler time discretization, the penalization parameter must satisfy Δ < 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 Δ 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: This scheme is unconditionally stable [7,10]. Finally, (10) is rewritten as 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 Δ .

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: where is the total pressure including the static and dynamic pressure, which can be defined as where ∞ and ∞ 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 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 + Δ ). Increasing leads to /(1 + Δ ) ≈ 1/Δ , equivalent to the parameter of an explicit scheme = 1/Δ . Hence, the penalization parameter in the pressure Poisson equation is chosen to be = 1/Δ 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, Re ℎ = | |ℎ 2 /] where ℎ is a fixed grid spacing, and the relation Re ℎ = O(1) should be fulfilled to ensure numerical stability [27]. The time step size Δ is also required to be ]Δ /ℎ 2 = O(1) for the diffusion. It is basically found that | |Δ = O(1) [7,27]. This means that the penalization parameter is equivalent to = 1/Δ = 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/Δ , 1/Δ , and 10/Δ . It is thought that = 1/Δ is a good choice. This is numerically investigated in more detail in Section 4.
Dirichlet boundary conditions are also required to solve the pressure Poisson equations (see (16)). It is obvious that the total pressure becomes zero at infinity. Although the total pressure 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 ( = 0) may be used. The computational domain size will be further discussed in Section 4.5.

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 time = 0. In other words, an impulsive start of the cylinder is simulated by specifying uniform inflow velocities (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 = ∞ /], based on the free stream velocity, ∞ , and the diameter of the cylinder, , 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. are determined through the stability condition Δ = ℎ 2 /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 / and separation angles is 4.16 and 51.7 ∘ , respectively, in our simulation. Fornberg (1980) [28] made a similar remark and gave / = 4.48 and = 51.5 ∘ experimentally. The drag coefficient which is defined as = 2 /( 2 ∞ ) is computed as = 1.483. This is close to the experimentally measured value of 1.498 in [28]. Figure 2(a) shows the computed pressure coefficient ( ) distribution for Re = 40. Here, it is defined as = 2( − ∞ )/ 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. In 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 5ℎ from the body surface is not removed.

Unsteady Flow Past a Circular
Cylinder at Re = 550. The simulation parameters are Δ = 0.002, ℎ/ = 0.005, and = 10 8 . The time step is determined by the condition Δ < ℎ 2 /4] ≈ 3.438 × 10 −3 . Simulation was carried out until = 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. Figure 3 shows that the streamline pattern computed for Re = 550 at = 5 compares very well with the flow visualization result of Bouard and Coutanceau (1980) [17]. The vortex core position indicated by the coordinates and is investigated. The nondimensionalized abscissa / and ordinate / are 0.34 and 0.27, respectively, and the wake length / is 0.82 in our simulation. Bouard and Coutanceau made a similar remark and gave / = 0.36, / = 0.28, and / = 0.85. Figure 4 shows the time evolution of the drag coefficient for an impulsively started flow around a twodimensional 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. Figure 5 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 = 0 have the same upstream moving trends. As experimentally measured by Norberg (2002) [34], = 0 in the present simulation is located at the angular position ≈ 35 ∘ .

Unsteady Flow Past a Circular
Cylinder at Re = 3000. The simulation parameters are Δ = 0.002, ℎ/ = 0.0025, and = 10 8 . Simulation was carried out until = 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].
Pressure fields are explicitly computed through vorticity and velocity fields of well-simulated flows for Re = 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  [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. Figure 12 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.

Effect of Domain
Size. The numerical parameters studied here are the outer boundary location, grid spacing, ℎ, and time step value, Δ . 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 to 0.005 and the time step value is 0.002.  Figure 14(b). Enlarging the outer boundary barely has any effect on the results of vorticity as shown in Figure 15.

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.