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 , 550, 3000, and 9500.
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 where the vorticity is defined as and the kinematic viscosity is assumed to be constant. This equation governs the evolution of vorticity in an incompressible viscous flow. The velocity is coupled with at all times in a self-consistent way through the relation of and . 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 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–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: and penalization. The penalization consists of adding a damping term on the velocity in the momentum equation, whereas the 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 penalization. Hence, most of the studies have practically adopted the penalization. In the 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.
2. 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 ) and diffused (according to ). 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 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 and its strength , that is, a circulation where is the area of fluid associated with the th particle, . The vorticity-carrying fluid element at a point 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 vector is expressed as where is the free stream velocity, 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 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 , , and . The velocity field on a grid is computed from the definition 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 where is the velocity of the solid body and denotes a mask function that yields in the fluid and in the solid [7, 12]. indicates the penalization parameter with dimensions , equivalent to an inverse permeability . In the case of explicit Euler time discretization, the penalization parameter must satisfy 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: where 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 .
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: 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, , 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 in a whole domain satisfies the incompressible mass conservation equation, [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 . Increasing leads to , equivalent to the parameter of an explicit scheme . Hence, the penalization parameter in the pressure Poisson equation is chosen to be 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, where is a fixed grid spacing, and the relation should be fulfilled to ensure numerical stability [27]. The time step size is also required to be for the diffusion. It is basically found that [7, 27]. This means that the penalization parameter is equivalent to . 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 , , , and . It is thought that is a good choice. This is numerically investigated in more detail in Section 4.
(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 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 () 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 time . In other words, an impulsive start of the cylinder is simulated by specifying uniform inflow velocities () in the whole computational domain. Boundary conditions at the body surface are imposed by the penalization terms.
In this study, the Reynolds number, , based on the free stream velocity, , and the diameter of the cylinder, , is selected to be , , , or . We only focus on the flow in the early stage of development, which is known to remain symmetric for up to . In our simulations, no symmetry constraint is imposed. The time () is nondimensionalized based on the radius () of the cylinder as . 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 for , , and and for . The penalization parameter for the diffusion is fixed to .
4.1. Steady Flow Past a Circular Cylinder at Re = 40
The numerical parameters are and , which are determined through the stability condition [7]. At , 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 and , respectively, in our simulation. Fornberg (1980) [28] made a similar remark and gave and experimentally. The drag coefficient which is defined as is computed as . This is close to the experimentally measured value of in [28].
Figure 2(a) shows the computed pressure coefficient () distribution for . Here, it is defined as . 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.
| (a) distribution |
| (b) along the body surface |
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, , it is deleted to avoid having too many particles. In the case of , 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 from the body surface is not removed.
4.2. Unsteady Flow Past a Circular Cylinder at Re = 550
The simulation parameters are , , and . The time step is determined by the condition . Simulation was carried out until 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 at 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 , , and . Figure 4 shows the time evolution of the drag coefficient for an impulsively started flow around a two-dimensional cylinder for 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 is very well simulated.
(a)
(b)
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 . 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 have the same upstream moving trends. As experimentally measured by Norberg (2002) [34], in the present simulation is located at the angular position .
(a)
(b)
(c)
(d)
4.3. Unsteady Flow Past a Circular Cylinder at Re = 3000
The simulation parameters are , , and . Simulation was carried out until 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 , the secondary vortices appear at an earlier time and grow larger. In the case of , 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 . Figure 7 shows that the streamline pattern computed for compares excellently with the flow visualization result of Loc and Bouard (1985) [21].
(a)
(b)
(c)
Pressure fields are explicitly computed through vorticity and velocity fields of well-simulated flows for . Figure 8 shows time evolutions of the pressure distribution around the cylinder body. As for the case of , 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 are also matched very well. These agreements help to demonstrate the quality of the numerical results obtained by the present algorithm.
(a)
(b)
(c)
(d)
4.4. Unsteady Flow Past a Circular Cylinder at Re = 9500
The simulation parameters are , , and . Flows were numerically simulated until , 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 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 .
(a)
(b)
(c)
(d)
(a)
(b)
Figure 12 shows the pressure distribution around the cylinder body, which is computed using the well-simulated flow for . 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.
(a)
(b)
(c)
(d)
4.5. 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 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 and the time step value is . We tested the FFT points of , , and . Actual computational domains correspond to , , and , respectively. Figure 14 confirms the convergence of numerical results for . The computational domain , 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.
(a) Comparison of pressure contours at different domain sizes
(b) Pressure coefficients on the body surface
| (a) |
| (b) |
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: , 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research was supported by the Korean Institute of Ocean Science & Technology sponsored by the Ministry of Trade, Industry & Energy (no. 10033668) and the Industrial Convergence Strategic Technology Development Program (no. 10044499) funded by the Ministry of Trade, Industry & Energy (MI, Korea).