Flow simulation around circular cylinder at low Reynolds numbers using vortex particle method

The algorithm of the Viscous Vortex Domain method implemented in the author’s code VM2D is investigated. This code allows simulating 2D viscous incompressible flows and solving fluid-structure interaction problems. The verification of the algorithm, namely, the correctness of taking into account the viscosity effect, is performed on well-studied test problems of the flow simulation around a circular cylinder for Reynolds numbers in the range 20 . . . 200: the problem of the impulse start of the circular cylinder is considered, and the dependency of the flow separation angle on the Reynolds number is investigated. The comparison shows that for low Reynolds numbers (up to 140), the flow separation angle obtained in numerical simulations by using the VM2D code is in good agreement with the results of other studies. For higher Reynolds numbers, the estimated separation angle is slightly overestimated.


Introduction
The problem of simulation of the viscous incompressible flow around bodies and determining the hydrodynamic loads is one of the most important problems of computational fluid dynamics (CFD) and structures design. Despite the long history of development of various methods for numerical solving of such problems, new approaches and computationally efficient algorithms are being developed up to the present day.
Lagrangian vortex particle methods of CFD have gained popularity in recent decades. Their relevance and higher efficiency in comparison to Eulerian mesh methods for solving some classes of problems are caused by several factors: • meshless nature, that is undoubtedly convenient when simulating the flow around moving bodies or solving coupled hydroelastic problems; • perturbations decay condition at infinity is satisfied automatically, so there is no need to artificially bound the flow domain when solving external flow problems; • computing resources are concentrated only in the domain with non-zero vorticity (usually it is localized around and behind the body); • since the vortex methods are particle methods, the algorithm of the vortex method can be effectively adapted for computations on graphics accelerators.
All these points raise the attention of researchers and engineers to vortex particle methods (VPM). New approaches to improving the VPM are developing. There are several different 2 modifications of the vortex methods, the key difference between them is the approaches of viscosity effect accounting. These approaches can be split into two groups: stochastic and deterministic. Among the stochastic approaches, the random walk method developed by Chorin [1] should be mentioned, where the vorticity diffusion is modelled as a stochastic process. Among the deterministic approaches, we note the Particle Strength Exchange method [2], Vorticity Redistribution Method [3], and the approach based on the introduction of the diffusive velocity [4,5].
Although a lot of researchers are working on the development of VPM, today it is still difficult to find a high-quality software implementation of VPM that would be available to a wide range of users, would allow solving a wide class of problems by using the capabilities of modern computing systems. Of course, each scientific group working in the field of vortex methods has its own codes, but in most cases, it is not open-source and can be used only by members of these groups.
The authors have developed a software implementation of VPM, called VM2D [6], that satisfies the above-listed requirements. As the basic modification of the VPM, the Viscous Vortex Domain (VVD) method is used, which makes it possible to simulate two-dimensional flows and is based on the introduction of the diffusive velocity [7]. Note that the authors have made some modifications to the VVD method, primarily concerning the satisfaction of the boundary condition. To date, the basic data structures and the computational kernel of the algorithm are developed, however, in order to serve this program as a reliable tool for flows simulating and coupled hydroelastic problems solving, it must be thoroughly verified and the range of model problems should be determined. In [8], the VM2D was verified on the test problem of flow modeling around two closely spaced airfoils and flow modelling in a channel with a backward-facing step. The results are in good agreement with the results of other researchers.
The accuracy of viscosity accounting makes a significant contribution to the accuracy of flow simulation at low Reynolds numbers. The goal of this work is to check the correctness of taking into account the viscosity in the algorithm of the Viscous Vortex Domains method, which is implemented in the VM2D code. So, in this paper the well-studied problem of the flow simulation around a circular cylinder is considered, and the position of the flow separation point for small Reynolds numbers is investigated. These tests make it possible to estimate the correctness of viscous terms modelling by using VVD method implemented in the VM2D code.

Viscous Vortex Domain Method and the VM2D code
The VVD method makes it possible to simulate two-dimensional flows of a viscous incompressible medium and to solve coupled hydroelastic problems. It is based on the introduction of the diffusive velocity [7]. The main calculated variable is vorticity Ω, its distribution in the flow domain is simulated with a large number of vortex elements with positions r i and circulations Γ i . The motion of vortex elements is carried out along the velocity field, which is a superposition of the convective velocity of the flow V and the diffusive velocity W . The convective velocity can be calculated according to the Biot -Savart law, and the diffusive velocity is determined as the ratio of the vorticity gradient to the vorticity where ν is kinematic viscosity coefficient. Formulae for calculating the diffusive velocity are described in [5,7]. At each time step of the calculation, as a result of satisfying the no-slip boundary condition, new vorticity is generated at the airfoil boundary, which then becomes part of the vortex wake. The boundary condition is written in the form of a boundary integral equation for the intensity of the vortex sheet. Two types of boundary integral equation are possible: expressing the equality of the normal and tangential components of the airfoil velocity and the medium velocity at the airfoil boundary. In this work, an approach based on the equality of tangent components is used, since it leads to a boundary Fredholm-type integral equation of the second kind with an integrable kernel (the approach with normal components leads to a singular equation of the first kind).
A detailed description of the approaches and algorithms implemented in the VM2D code can be found in [9]. In this work, we do not discuss the computational efficiency of the VM2D code; moreover, we note that for some computational blocks, algorithms that are more efficient from a computational point of view can be used, for example, to calculate the mutual influence of vortex elements in calculating convective velocities, one can use approximate fast methods: FFT-based method [10], Barnes -Hut algorithm adapted to vortex methods [11], fast multipole method. In this research, the acceleration of the used brute force algorithms (for example, direct calculation of the mutual influence of all vortex elements by Biot -Savart law) is carried out by parallelizing them with OpenMP and MPI technologies and adapting them to perform calculations on graphics accelerators using the Nvidia CUDA technology [12].

Unsteady drag force computation for impulsively started cylinder
To verify the code, the well-known test case of unsteady flow simulation for impulsively started cylinder was considered firstly. In [13] a lot of results of both numerical and physical experiments are presented, we use for verification such of them that correspond to low Reynolds number: Re = 40 and Re = 200. Note, that for these regimes pressure-based drag force and viscous friction force are of the same order.
In fig. 1 the results of unsteady drag coefficient are shown in comparison to [13]. One can see, that the results are in acceptable agreement.

Determination of separation angle for flow around a circular cylinder
In the case of flow simulation around an airfoil that does not have sharp edges or corner points, the position of a separation point depends not only on the shape of the airfoil but also on the Reynolds number. In numerical simulations, the obtained position of a separation point can say a lot about the correctness of the flow regime simulation on the whole and the accuracy of the hydrodynamic loads estimation.
In this paper, in continuation of a series of problems considered for verification of the VM2D code, we investigate the accuracy of determining of a separation point for a circular cylinder for low Reynolds numbers (Re = 20, . . . , 200). In this range, there are different flow regimes: steady regime, when the flow behind the airfoil is symmetric, and regime with vortex shedding when von Karman vortex street is observed. In the case of Reynolds numbers up to 40, the separation point maintains its position, while for regimes with vortex shedding, the separation point moves significantly during the vortex-shedding period.
Separation angle θ is measured as shown in the Fig. 2. The separation point at the airfoil surface line is determined as the point on the airfoil boundary where the shear stress changes its sign.  Table 1 shows the results of determining the separation angle in simulations with the mentioned discretization parameters for Reynolds numbers 20, 60, 120 and 200.  Table 1, one can see that not for all modes the obtained separation point is independent of discretization parameters. For further calculations, we will use the most detailed discretization (N = 1000 and τ = 0.002), since the results obtained by this discretization seem to be the most consistent with the results obtained by other researchers.
The Fig. 3 shows the vortex particles distribution around the airfoil and the positions of separation points for different time instants during one vortex-shedding period T for the case of flow simulation with Re = 200. To determine the position of the separation point at time instant, we find the point on the airfoil boundary, at which the share stresses vanish.
The Table 2 shows the separation angles for various Reynolds numbers. For the regimes without Karman vortex street, the separation angleθ is shown, which just slightly varies in time during the numerical simulation. For unsteady regimes with Karman vortex street formation, when the position of the separation point varies significantly during the vortex-shedding period, the average valueθ of the separation angle is given, as well as its oscillation range θ min . . . θ max .
A detailed review of the results obtained numerically and experimentally by different authors for the considered problem is given in [14]. The Fig. 4 shows a comparison of the results obtained in current study from numerical simulation using VM2D code, with experimental and numerical results of other authors: experimental investigations by Thom, 1933 [15]

Conclusions
Using the VM2D code, a two-dimensional flow around a circular airfoil for Reynolds numbers in the range of 20. . . 200 was simulated. Results of the drag force acting on the impulsively started cylinder obtained by using the VM2D code corresponds to the experimental data and simulations of the other researchers. For these low-Reynolds regimes, the angles of flow separation are determined and the results are compared to the results of other studies. For small Reynolds numbers (up to 140), the results obtained using VM2D are in good agreement with the results of other studies. For higher Reynolds numbers, the obtained separation angles are slightly overestimated, which creates a field for further research aimed at improving the VM2D code to obtain more correct simulation results.