Vortex methods in CFD problems

Different mesh-using methods are widely applied to solve CFD problems. However, these methods are really hard to implement because of the problems associated with mesh construction, especially in complex-shaped domains. Using the mesh-free methods, e.g. vortex methods, could be the way to avoid these difficulties. This paper is focused on investigation the implementation features of vortex methods. Application of some vortex methods to modelling the incompressible flows around bodies is considered. Sometimes authors of vortex methods do not pay enough attention to details, which could complicate the practical application. These aspects of considered vortex methods are also discussed. In particular, it is shown that separate calculation of vortices with positive and negative circulations (the separate calculation approach) significantly improves the reliability of computations. Mathematical background of this approach is provided. The flow around the stable circular cylinder is considered as the test case. Several flow models are examined.

Here is a kinematic viscosity of the fluid. In two-dimensional case the vorticity is always directed orthogonal to the plane of the flow, so z  Ωe , where z e is the corresponding unit vector. Thus, equation (1) is equivalent to the following scalar equation for the  : (2) Only the case of a two-dimensional flow will be considered in this paper. It is important to mention here, that the VVD method could be also applied to modelling of axisymmetric three-dimensional flows with some insignificant changes [18].
In the VVD method equation (1) is solved in a following manner. The set of points called pointwise vortices or just vortices is introduced. Each vortex represents the fluid's domain and has its own position ( ( )) ( ( )), 1..
where d     V is called the diffusion velocity of vorticity [16]. Also N is the total number of vortices. One can show that if vortices move according to equations (3) Here  V is the flow velocity at the infinity. It could be shown [6] that this velocity field satisfies the continuity equation in the case of incompressible flow. Expressions for the diffusion velocity are more complicated (see [8] for details). Let us describe the fulfilment of the boundary conditions in details. The body's surface is approximated by a polygon with K vertices in the VVD method. The vortices (they are called surface ones) with the unknown circulations i  are generated at these vertices at each time-step. Surface vortices are used for modelling the forming of vorticity, which happens only on the body's surface due to the discontinuity of the fluid velocity [8]. Unknown circulations are obtained from a specific system of linear equations. It is achieved by writing the non-permeability condition 0  Vn in the control points on the body's surface and substituting the Biot-Savart law there ( n is outer normal to the body's surface).
Usually the control points are selected as the mid-points of the sides of the approximating polygon. One can recall, that more strict no-slip condition 0  V is required on the body's surface for the viscous flow, but it could be shown that this condition also holds [6].
After the unknown circulations are determined, surface vortices and other previously emerged vortices (free ones) are moved in the fluid according to equations (3). Some vortices (surface or free ones) could intersect the body's surface after their moving due to the approximation error. It is suggested in [6] to move such vortices anyway and then exclude them from the consideration. We faced significant numerical oscillations in vortices' forming using this method, and that's why we suggest to put these vortices on their places and not to move them instead. It is important to mention that they could be taken away into the flow after some time.
Also the VVD method allows the determination of some important characteristics of the flow, like the pressure in an arbitrary point using an analogue of the Cauchy-Lagrange integral, as well as the pressure force and the drag force acting on the body [6].
The VVD method has one serious problem. The diffusion velocity, as it could be seen from its definition above (see eq. (3)), is evaluated by dividing some expression by the local vorticity value. This local value could be close to zero in numerical calculations, which means that there is no vortex in this place and the diffusion velocity has to be equal to zero as well. But instead of it the unreasonably large value for diffusion velocity would be obtained after the division.
Ogami mentions this problem in [17]. He suggested to calculate the diffusion velocity separately for the vortices, which have positive and negative circulations (separate calculation approach). He also showed that this trick works in the test case, but mathematical background is required to use it confidently. It is done below for the case of smooth vorticity field. Denote , 22 After substituting this identity to equation (2), the following equation is obtained: Thus, equation (5) in this point has the following form: Region, where 0  , is considered in the same way. It follows that equation (5) could be transformed: Equation (7) describes the transport of positive and negative vorticity, but not the zero one. This region is a boundary between the regions of positive and negative vorticity, so it's moving could be also restored. Finally, the VVD method could be applied to the both equations in system (7). This procedure is found to be more stable compared to the original scheme of the VVD method. The fact that computed values of   and   are better separated from zero than the values of  could be the possible explanation for this phenomena.

The VD and VVHD methods
The case of an inviscid flow is studied in the same manner after 0   is assumed. The corresponding method is called the vortex domains (VD) method [5]. Only one aspect has to be considered in details. It is connected with an evaluation of the fluid velocity in the control points on the body's surface.
As it was mentioned previously, the velocity field calculated using the Biot-Savart law satisfies the no-slip condition on the body's surface, and this is not the case for the inviscid flow. Thus, some specific correction is required in this case. It could be found in [5]. The fluid velocity in the k-th control point may be evaluated from the formula: Here c k r is the position of the k-th control point, 0 () c k Vr is the fluid velocity in this point calculated with equation (4), is the local value of surface vorticity, which is determined from the circulations surface vortices placed in k-th and k+1-th polygon vertices respectively ( τ is a tangent vector to the body's surface in the k-th control point.
The case of non-isothermal flow also could be handled. The energy equation in the viscous incompressible fluid has the following form [6]: Here T is the fluid's temperature, V is the fluid's velocity again, a is the thermal diffusivity, diss N is the thermal energy, which is dissipated through the unit area during unit time in the fluid. Equation (9) should be supplied with some boundary conditions. The case, where the body has fixed temperature w T , is considered in this paper. It is also assumed that the fluid's temperature is zero at the infinity.
One can show that the energy dissipation could be neglected if where C is a fluid's thermal capacity [6]. This case is considered in this paper. Thus, equation (9) is simplified: Equation (10) has the same form as the equation (1) and could be solved in the same way after the heat diffusion velocity T T a T   V is defined [15].
The scheme of solving the equation (10)  generated in the thin layer near the body's surface and its heat energies are selected to satisfy the corresponding boundary condition for the temperature (see [15] for more details).
Notice that the separate calculation of heat vortices is not needed in the heat problem. Actually, all heat vortices have heat energies with the same sign equal to the sign of body's temperature w T . It is also worth mentioning that heat vortices' moving is partially determined from the fluid vortices' parameters, because the fluid velocity is required there. Thus, solving the heat problem requires the simultaneous solving of the fluid dynamics problem.

Numerical results and discussion
Specific software was developed in order to implement the considered vortex methods. The obtained results for modelling the flow over the stable circular cylinder is considered in this section. This problem is one of the test problems usually used in computational fluid dynamics in order to verify the implemented numerical methods.

Pressure coefficient for the inviscid flow
In the case of inviscid flow, we compute pressure coefficient on the surface of circular cylinder and compare it with analytic solution. A pressure coefficient is one of the important characteristics of the flow. It is defined by the formula: Here R is an arbitrary point in space,   p R is the pressure in this point, p  is the pressure at the infinity, ρ is the fluid density. The pressure coefficient is usually used to demonstrate the pressure distribution on the body's surface. There is an analytic expression for this coefficient [20]: Here  is a polar angle, which is calculated from the point on the cylinder's surface, where the corresponding outer normal is directed towards the flow. The considered problem was solved using the VD method. The body's surface was approximated by the regular polygon with K=50 vertices. The pressure coefficient was calculated in the control points using the Cauchy-Lagrange integral. It has the following form in this case: Notice that the fluid velocity in equation (13) is calculated using the correction to the Biot-Savart law (see eq. (8)). The calculated distribution of the pressure coefficient compared with the values given by equation (13) are given in

Modelling the viscous flow with Re = 250
A viscous flow around the circular cylinder has many differences from the inviscid flow. The Reynolds number

D Re
is the major characteristic in this case (here D is the diameter of the cylinder) [20]. The flow structure depends on the relationship between the viscous force and the inertia force, which is characterized by the Reynolds number.
The results for the Reynolds number Re=250 will be considered next, because according to the classification of flow regimes given in [21] the periodic vortex shedding is noticed in this case (so called the von Karman vortex street). This problem was solved using the VVD method. The cylinder's surface was approximated with the regular polygon like in the previous case. The number of vertices in the approximating polygon was increased to the value K=500 in order to resolve the complicated structure of the occurring flow. The vortex pictures [6] were obtained on each time-step to visualize the flow. They depict positions of all vortices in the flow. Thus, vortex pictures show the support of the vorticity field. The symmetric recirculation zone is formed after the cylinder instantly after the start of the simulation. Then it loses its stability with time, and the vortices, placed in the unstable part of the zone, move away forming the structure looking like the von Karman vortex street. It is shown in Fig.  2 (upper part).
Vortices in the street form compact separated structures (the vortex blobs). It is noticed that all vortices in each blob have one sign of circulation, which is an important property of the von Karman vortex street. Some vortices are placed outside of the street. It happens because the vortices' equations of motion are integrated with the explicit Euler method, which has only the first order of approximation. It was much more vortices outside of the vortex street before the separate calculation approach was applied (see subsec. 2.1). The visualized flow over stable circular cylinder from [22] is shown at the lower part of Fig. 2. Comparison of computed flow pattern with experimental one demonstrates that the VVD method describes the structure of the occurring flow correctly.

Fig. 2.
Comparison of the vortex picture calculated using the VVD method with experimental visualization from [22].

Separation point for different Reynolds numbers
Another test case is considered to confirm the fact that the VVD method allows to predict the characteristics of the occurring flow with high precision for different Reynolds numbers. The separation point of the flow is calculated for the Reynolds numbers Re=40,150,250 .
The separation point (where the boundary layer becomes detached from the body's surface) could be relatively easily found by the VVD method. The vertices of the approximating polygon are examined in order, where the polar angle is increased, and the first vertex, where the corresponding surface vortex's circulation has different sign with the previous one, is called the separation point. Actually, the tangential component of the fluid velocity is changing its sign near this point, so the flow is separated from the cylinder there.
The computed results are displayed in Fig. 3 compared to the experimental data from [23]. The dependency on the Reynolds number is considered. The experimental data are plotted with circles, and the computed data are drawn with crosses. Notice that the computed data are very close to the experiment, and the difference decreases with the Reynolds number increasing. It could be explained by the fact that smooth increasing of the fluid velocity at the beginning of the simulation was not provided. Actually, this increasing is not required at high Reynolds numbers, because in this case the flow is not smooth enough from the very beginning.

Heat transfer problem with Re=250, Pr = 5
The heat transfer problem was solved in order to test the VVHD method. We note that the heat transfer essentially depends also on the Prandtl number Pr a   [20]. The test problem with Re=250 and 5 Pr  is considered. The cylinder has fixed temperature 20 w TC  . Other parameters are selected the same as in the example from subsec. 3.2.
The temperature distribution could be obtained on each time-step using the positions and heat energies of the heat vortices (see [15] for details) and one of them is shown in Fig. 4. The temperature field follows the flow structure and the heat von Karman vortex street is also obtained. Notice that the heat vortex blobs are more compact than the corresponding vortex blobs in the fluid problem (see Fig. 2). It is explained by the fact that the Prandtl number is taken more than one. Because of that heat diffusion velocities have smaller values than the diffusion velocities (see eq. (3), (10)), and subsequently heat vortices move closer to each other than the fluid ones.

Conclusion
Using the vortex methods in CFD problems allows to avoid difficulties connected with the mesh-using methods that are usually applied in this case. Sometimes authors of vortex methods do not pay enough attention to details, which may be significant in practical applications. These aspects of the VD, VVD and VVHD methods are discussed in this paper. Note the separate calculation approach, which has to be applied in the VVD method in order to obtain the stable vortex structure. The mathematical background of this approach is provided in this paper for the first time. Some details of fulfillment of the boundary conditions on the body's surface are also discussed.
It is shown that considered methods allow simulation of important effects occurring in flows (like the von Karman vortex street) and prediction of important properties of the flow (like the position of the separation point). Good agreement with theoretical data is also presented in the case of inviscid flow.
The disadvantage of the VVD method is the fact that it could be applied only to the two-dimensional flows or to axisymmetric three-dimensional ones. The vorticity field becomes non solenoidal, when vortex methods are applied to general three-dimensional case, and mesh-using methods have to be employed to restore this property. However, there are some attempts to extend vortex methods to three-dimensional flows, among them the dipole domains method (DDM) is worth mentioning [6]. Presently the intensive research is conducted in this direction.