Discrete Vortex Prediction of Flow around Two Cylinders in Side by Side using Simple Grid System

Email: wisnu@oe.its.ac.id Abstract: Modelling of flows around two cylinders in side by side arrangement using a simple overlapping grid system is carried out. The use of such grid system is intended to give a considerable reduction in terms of the CPU time, especially during the calculation of the vortex velocity. It has been shown that this method is not only time efficient, but also gives a better distribution of surface vorticity as the scattered vortices around the cylinder are now concentrated on grid point located at uniform distance from the cylinder. The engineering applications of this topic is to simulate the loading on structural elements due to the presence of anodes. The in-line and transverse force coefficients and the flow patterns obtained are presented in order to provide more detail description of the flow phenomena and interaction involved. The comparison of the results with both experimental and numerical evidence is also presented and the range within which the algorithm produces good results is identified.


Introduction
During the last five decades, various authors have attempted to model the flow around cylinders. Stansby and Slaouti (1993) provided a numerical model using a vortex in cell and random walk technique to simulate the convection and diffusion process of the flow around two cylinders by solving the Poisson's equation which relates vorticity ω to the stream function ψ through 2 ψ ω ∇ = − . Three types of overlapping mesh were adopted so that a very fine mesh was used to give definition in the boundary layers, an intermediate mesh size was used for computing vortex roll up in the near wake region and a coarse mesh was used to transport vortex structures downstream. Meneghini et al. (1997;Williamson and Govardhan, 2004) applied similar method for predicting the flow around an oscillating circular cylinder.
Recently, Kostecki (2014) used random vortex method for a single cylinder involving turbulent models for high Reynolds number while Laroussi (2015) investigated the flows around two cylinders in close proximity under the influence of initial conditions using a finite element based program, at low Reynolds Numbers. Gordo (2011) on the other side, applied the method to explore the flow around airfoils using meshless methodology called full cloud vortex method.
In this study, a simple overlapping polar grid system in which the grid has 'square' elements which increase in area in the radial direction, is used for each individual cylinder to give good definition of the flow close to the cylinder surface. The grid node on the cylinder surface located at the center of the surface element provides the control points at which the Martensen equation is solved to give zero tangential velocity and so satisfy the Dirichlet boundary condition, Wardhana (1995).

Basic Formulation
As there is more than one cylinder in the flow domain, the complex potential w(z) can be written as follows: The solution procedure for the induced velocity at an element S n must include the influence of those two cylinders in the fluid domain. Mathematically, this can be expressed in the following Martensen equation: The introduction of vortices with circulation satisfying the boundary condition of zero tangential velocity is carried out by releasing a ring of vortices on the second ring from the cylinder surface. The identification of vortices shed from any given cylinder is maintained over the whole process. The implementation of the Biot-Savart velocity calculation over the polar grid nodes is done only among those vortices shed from the same cylinder. The area of proportions and the vortex strength distribution have to be calculated in either its own polar or rectangular grid system at every time step.

Methodology
The structure of the Discrete Vortex Algorithm is displayed in Fig. 1 above. The main calculations consist of the followings.

Surface Velocity Calculation
Equation 2 must satisfy the Dirichlet boundary condition at element m of cylinder p as follow: By absorbing the first term into the second term of coupling coefficient pq mn k , the equation may be represented as follow: In matrix from the two-cylinder problem can be written as: The RHS means the right-hand side of Equation 5 for the first and second cylinder respectively.

Segmentation of the Domain
A simple overlapping grid system to represent the fluid domain in the case of two-cylinder problems can now be carried out. Each cylinder has its own polar grid system and forms a local domain. The center of one of the cylinders in chosen as the global reference of the domain in which the coordinate of the rectangular grid nodes are measured and stored.
The size of grid segment in the polar grid system is expanding linearly from the surface of the cylinder. In the two-cylinder problem, there will be two overlapping polar grid systems expanding linearly from each cylinder surface. A rectangular grid system which overlaps the two polar grid systems is also incorporated into the system. This rectangular grid has a uniform square element through-out the whole domain. The structure of the grid system form twocylinder problem is displayed in Fig. 2 above. Any vortex that falls into a position where two grid elements of the polar and rectangular system overlap each other is located with reference to two base nodes for its relative position. The node points of the rectangular coordinate system are measured from the center of the referenced cylinder. The coordinates of each node were defined starting from the top-left corner to right-bottom corner to form a rectangular domain. The distance is approximately equal to u ∞ × t, where u ∞ is the free stream speed and t is the total time used, plus k x diameter, in which the value of k is approximated from experiments.

Distribution of Circulation to the Grid
In this multi-cylinder problems, a vortex of strength Γ v is not only distributed onto its four surrounding nodes in the polar grid, using the common bi-linear interpolation, but also onto nodes of the overlapping rectangular grid system. This means that each vortex has two base nodes i.e., the polar grid where it is shed from and the rectangular grid (Fig. 3).
The rectangular grid is used in the calculation of a vortex velocity due to the contribution from all other vortices shed from cylinders other than the one from which it originated. This implies that in the same segment, only vortices shed from the originating cylinder will have their strengths are summed on the surrounding grid nodes will be stored in the same array. Contribution from vortices shed from other cylinders will be summed in different arrays. This can also mean that the same rectangular grid nodes could be active with reference to one cylinder but be inactive with reference to the other cylinders.

Calculation of Velocity
Since each cylinder has its own polar grid which overlap with each other and the rectangular grid, the calculation of the velocities of a vortex shed from one cylinder due to the other vortices shed from the same cylinder is carried out with reference to their own individual polar grid. In addition, there is a contribution from these shed from the other cylinder and this is computed with reference to the rectangular grid system: • Calculate the velocity u pp (N v ) at the nodes surrounding a vortex z v , shed by the cylinder p and located at a polar grid segment with base node N' 1 , due to the other active nodes of the polar grid system of the same cylinder p using: where, µ = 1-4. • Using the rectangular grid system, calculate the velocity u pq (N µ '), at the nodes surrounding a vortex at z v with base node N 1 ', due to the other active nodes of the other cylinder q using: • Use bi-linear interpolations to find the vortex velocity due to both polar and rectangular active grid nodes by using the area of proportion P µ (v) and Q µ (v) using: As implied in the above procedures, there is no direct interaction between the nodes of the polar grids associated with each cylinder. Instead, the cylinder to cylinder node interactions are computed using the overlapping rectangular grid nodes. This strategy is chosen in an attempt to achieve a more uniform distribution of vorticity in the area in which the wakes of the cylinders interact.

Time Integration
A first order accurate Euler scheme is used to find the new locations of the vortices, that is: Since an overlapped grid system is used in the scheme, the new position of each vortex is then referenced to both polar and rectangular grid systems. In other words, each vortex has two base nodes from which its relative position at every time step is measured and renewed.

Force and Pressure Calculations
The force calculation can be carried out after the convection and diffusion processes by solving the following equation: The other force that contributes to drag and lift forces is the one due to the skin friction (viscous drag) on the surface of the cylinder. This force comes from the shearing stress at the surface: where, µ is the dynamic viscosity. The form drag, lift and skin friction coefficients can be calculated as follows: where, d is the diameter of the cylinder and β the tangent angle of the element. The basic procedure is based on the integration of the elemental pressure around a cylinder. The pressure around the cylinder can then be integrated numerically to get the value of the force coefficients.

Method of Enhancements Correction for Close Proximity
This method is for finding the influence coefficient of the elements located in close proximity. Similar technique can be adopted when two or more cylinders are close together. In these circumstances the influence coefficient can be written: which can be interpreted as the velocity field of a vortex whose strength decays with time because of the diffusion process. In order to model those effects in his discrete vortex analyses, Naylor (1982) used on an empirical reduction scheme for the individual vortex strengths similar to the above equation and given by: The method used the reduction of vorticity as follow: where, C r is a constant less than 1 and n is the number of iterations. With the time step ∆t = 0.15, the effect of variation of C r on the vortex strength Γ(t) is displayed in Fig. 4 above.

Fig. 4. Vorticity reduction
The implementation of the vortex reduction scheme produces some improvements in maintaining the symmetrical properties of the wake shed from both cylinders. Numerically, this is achieved because of the significant reduction in the interaction of the two vortex streets shed by both cylinders when they are close together.

Curvature Corrections
The correction for the curvature of the elements and the asymmetric enforcement when the arrangement of the cylinders produces a symmetrical flow pattern are also implemented in the present case study. This effect is a result of an increase in the local velocity in the gap region which convects the vortex sheets in this inner region more strongly than those located in the outer one.

Results and Discussion
At Re = 25,000, as the flow develops earlier at ˆ1 t = and ˆ5 t = , due to strong asymmetrical properties of the flow field Fig. 5, the rolling up of the vortices appears immediately while the formation region is building up behind the cylinders.
By comparing with the case of an isolated cylinder Downie (1981;Wardhana, 1995), it is shown, Fig. 6 and 7 above, that the blockage effect increase the length of the formation region by about 10-20% even though the width is relatively unchanged. This means that the influence of the increase in the local velocity in the gap region is quite significant.
The comparison of the flow pattern simulated numerically with the experimental flow visualizations produced by Bearman and Wadcock (1973) at R e = 25,000 as shown in Fig. 8 above shows good agreement between the two and suggest that the grid sizes and configuration used in the model are appropriate.  The increase from the value of the isolated cylinder is close to the values deducted from the experiments at Reynolds number of 25,000. By using fast Fourier transform to evaluate the frequency f, the value of Strouhal Number can then be calculated and it settles to a value of around 0.2, similar to that of the isolated cylinder, as seen in Fig. 9 below.
The general trend shows that the interaction between the two cylinder is weakened as the gap becomes wider. Each individual cylinder is behaving increasingly like an isolated cylinder. The drag coefficient also approaches the value of around 1.14 or, in the other words, the interference drag coefficient becomes equal to zero Fig. 10.
The result shown in this study have been obtained by implementing a vortex reduction scheme, similar to those of Naylor (1982), in which the strength of the vortices is reduced in such manner that the effect on the force coefficient and the flow pattern is as close as possible to that of the experimental results.
The flow pattern created behind the cylinders are seen to be more regular compared to the same configuration in the convective flow shown in Fig. 11. This is due to the implementation of the vortex reduction scheme simulating turbulent dissipation and reducing the level of irregularity occurring in the group of vortices with high strengths created in the early stages of the calculation. The values of the Strouhal numbers and the drag and lift coefficients are plotted onto the original graph of the experimental results as displayed in Fig. 9 and 10. The out-of-phase asymmetric Karman vortex street behind the cylinders can be achieved for the convection only cases in which neither the Random Walk, simulating the diffusion process, nor the vortex strength reduction scheme are yet included. This can be considered to model the flow at very high Reynolds number. After the incorporation of the Random Walk scheme, for computing intermediate Reynolds number flows, the asymmetry of the flow could not be maintained without the involvement of the Vortex Strength Reduction scheme. It appears that vortices with large strengths generated by each cylinder in the early stages of the flow interact to cause instabilities in its subsequent development, Fig. 12-14.
The effect of the interaction of the two cylinders further apart than 0.5 D, can be reasonably well represented with the present model even though the magnitudes of the force coefficients are slightly on the low side. However, when two equal cylinders are closer than 0.5 D, the complicated effect of the boundary layer interactions could not be approached realistically using the present model.

Conclusion
One of the main difficulties in the flow around two cylinders in an infinite fluid implementing the present model has been in achieving results within practical time limits. As has been mentioned already, the algorithm does not include explicitly a turbulence model. Results resembling experimental results have been achieved over a wide range, but inevitably there are flow configurations for which turbulence effects will not allow representation of the flow in this manner. By using certain mathematical transformation, this method can also theoretically be developed further for predicting flows about any two-dimensional shapes.