Analysis of potential flow around two-dimensional body by finite element method

The paper presents a numerical method for analyzing the potential flow around two dimensional body such as single circular cylinder, NACA0012 hydrofoil and double circular cylinders by finite element method. The numerical technique is based upon a general formulation for the Laplace’s equation using Galerkin technique finite element approach. The solution of the systems of algebraic equations is approached by Gaussian elimination scheme. Laplace’s equation is expressed in terms of both steam function and velocity potential formulation. A finite element program is developed in order to analyze the result. The contours of stream and velocity potential function are drawn. The contour of stream function exhibits the characteristics of potential flow and does not intersect each other. The calculated pressure co-efficient shows the pressure decreasing around the forwarded face from the initial total pressure at the stagnation point and reaching a minimum pressure at the top of the cylinder.


INTRODUCTION
The flow past two dimensional body such as circular cylinders and hydrofoil has been the subject of numerous experimental and numerical studies because this type of flow exhibits the very fundamental mechanisms. The flow field over both the cylinder and hydrofoil is symmetric at low values of Reynolds number. As the Reynolds number increases, flow begins to separate behind the body causing vortex shedding which is an unsteady phenomenon. To achieve the goal of obtaining the detailed information of the flow field around two dimensional bodies, Finite Element Method (FEM) has been emerged as an attractive, powerful tool in many designing process.
The FEM was originated from the field of structural calculation in the beginning of the fifties and was introduced by Turner et al. (1956). The FEM was introduced into the field of computational fluid dynamics (CFD) by Chung (1977). The first study concerning the steady flow past a circular cylinder was reported by Thom (1933) for Reynolds number of 10 and 20. The works of Kawaguti (1953) and Payne (1958) were restricted to low Reynolds numbers (Re = 40) and relatively low Reynolds numbers (Re = 40~100) respectively. Oden (1969) has presented a theoretical finite element analogue for the Navier-Stokes equations, but without a practical numerical method. Dennis and Chung (1970) introduced finite element method into the field of computational fluid dynamics (CFD) by solving steady flow past a circular cylinder at Reynolds number (Re≤100).
` Tong (1971) presented results for steady flow using this method with pressure and velocities as dependent variables. Olson (1974) presented a numerical procedure to investigate steady incompressible flow problems using stream function formulation. Hafez (2004) simulates steady an inviscid flows over a cylinder using both potential and stream functions.
The objective of the present research is to analyze the potential flow around single circular cylinder, NACA 0012 hydrofoil and double circular cylinders by Galerkin technique of finite element method. Due to symmetry of circular cylinders and NACA 0012 hydrofoil, only the upper half portions have been considered as computational domain. Both stream function and velocity potential formulation have been used with definite boundary conditions. Contours of stream and velocity potential function, velocity above crest for both formulations and velocity and pressure distribution along the surface for various discretization are obtained which are compared with the analytical result available from the literature and shown graphically.

Potential flow around circular cylinder
Let us consider the potential flow of an ideal fluid around the circular cylinder placed with its axis perpendicular to the plane of the flow as shown in Figure 1. Now, the potential flow around circular cylinder can be represented by Laplace equation as: The velocity components u and v of the flow field in relation to stream function  or the velocity potential  are given by:

Boundary conditions for stream function
The half of the fluid domain is taken in the computations as shown in Figure

Boundary conditions for velocity potential
In case of velocity formulation, the boundary conditions that need to be satisfied in order to get the solution of Laplace equation: in Ω as shown in Figure 3:

Potential flow around a hydrofoil
Let us consider the flow of an ideal fluid around a hydrofoil placed with its axis perpendicular to the plane of the flow as shown Figure  4.

Boundary conditions for stream function
Now we need to solve the Laplace equation: in Ω as shown in Figure 5 with the following boundary conditions: Where a and b are the two constants. Now we need to solve Laplace's equation       ` (a) Ψ1 = U y on S1 (b) Ψ1 = 0 on S2 and S3 (c) Ψ2 = 0 on S1 and S3 (d) Ψ2 = 1 on S2 (e) Ψ3 = 0 on S1 and S2 (d) Ψ3 = 1 on S3

Numerical solution by stream function method
The stream function over the domain of interest is discretized into finite elements having M nodes: Using the Galerkin method, the element residual equations are: Where S represents the element boundary and (nx, ny) are the components of the outward unit vector normal to the boundary. Using Equation (4) in Equation (7) and substituting the velocity components into the boundary integrals, results in:

Numerical solution by velocity potential method
The finite element formulation of potential flow of an ideal fluid in terms of velocity potential is quite similar to that of the stream function approach, since the governing equation is Laplace's equation in both cases. By direct analogy with Equations (4) to (11) it is obtained as follows: Using the Galerkin method, the element residual equations are: Utilizing Equation (12) in the area integral of Equation (15)

RESULTS AND DISCUSSION
Based on the previous mathematical formulation as outlined as numerical solution by stream function method and numerical solution by velocity potential method a finite element program has been developed in FORTRAN 90 for calculating the potential flow around two dimensional bodies. For all finite element mesh configurations, nodes along the vertical line above the crest of the cylinder are numbered consecutively from top to bottom in order to be compatible with velocity calculations used in the program. The elements are taken in the form of triangle or quadrilateral for the convenience of discretization, thus the edge of the body may not be appeared as a circle or hydrofoil.

Single circular cylinder
Let us consider the flow around the circular cylinder of unit radius confined between two parallel plates having length of 7 m and height 4 m. A fluid of uniform velocity 1.0 m/s is assumed to be flowing from the left to the right of cylinder as shown in Figure 8. The choice of computational domain in the direction of flow is arbitrary and the free stream velocity is considered to prevail at distances sufficiently far from the cylinder. The upper half of the computational domain surrounded by the path (a-bc-d-e-f) is taken into account for numerical calculation due to symmetry of flow and is discretized by (24×5) triangular elements as shown in Figure 9 for stream function formulation. The contour of stream function has been obtained from stream function formulation and exhibits the characteristics of potential flow as shown in Figure 10. The stream lines have not intersected each other and mean the flow past the cylinder smoothly without any separation at the trailing edge. The upper half of the computational domain for velocity potential formulation is also discretized by (20×7) quadrilateral elements as shown in Figure 11. The contours of velocity potential have also been obtained from velocity potential formulation and exhibit the characteristic of potential flow i.e. no vortices exist at the trailing edge as shown in Figure 12.
The velocities along the vertical line above the crest of the cylinder are calculated and then compared with the analytical result in Figure 13. The average deviation for the velocity profiles between the two cases is less than one percent. Figure 14 is plotted by calculating the velocities above the crest at two points(x = 3.50, y = 1.00) and (x = 3.50, y = 2.00) against various number of nodes for stream function formulation which shows that computed velocities converges to the analytical solution as number of nodes increases. Figure 15 is obtained by plotting the velocity square along the surface of cylinder against the angular coordinates of nodal points. There are two types of curves of which first type shows a sinusoidal curve for the whole cylinder obtained theoretically and second type consists of four curves for four different finite element mesh configurations.
In Figure 16 the calculated pressure coefficient (C p ) is compared with the theoretical pressure distribution over the surface of the cylinder and the agreement is found to be quite satisfactory. The calculated results show the pressure decreasing around the forwarded face from the initial total pressure at the stagnation point and reaching a minimum pressure at the top of the cylinder.

NACA 0012 hydrofoil
Let us consider the flow around NACA 0012 hydrofoil confined between two parallel plates having length of 10 m and height 4m as shown in Figure 17. A fluid of uniform velocity 1 m/s is flowing from the left to right of the foil. The half of computational domain for stream function formulation is discretized by (16×3) elements as shown in Figure 18 and the contours of the stream lines are given in Figure 19.
Similarly, the half of computational domain for velocity potential formulation is discretized by (16×3) elements as shown in Figure 20 and the contours of the stream lines are given in Figure 21. Figure 22 depicts a comparison of pressure distribution pressure over the surface of the foil with the results obtained from constant strength source method 11 and shows very close agreement both at the leading and trailing edge of the foil.

Double circular cylinders
The flow around two circular cylinders of unit radius is confined between two parallel plates having length of 10 m and height 4m. The distance between the cylinders is unit length and the height above the cylinder is also unit length. A fluid of uniform velocity 1 m/s is flowing from left to right of cylinders as shown in Figure 23. The half of the computational domain is discretized by (16×3) elements for stream function formulation and (20×7) elements for velocity potential as shown in Figures 24 and 25, respectively. The contours obtained from these two formulations are drawn in Figures 26 and 27 respectively.

Conclusions
The paper presents a numerical method of calculating the potential flow around two dimensional bodies by finite element method. The following conclusions can be drawn from the present numerical analysis: (i) The present method can be an efficient tool for evaluating the potential flow characteristics of two dimensional body.                (iii) The calculated pressure co-efficient shows the pressure decreasing around the forwarded face from the initial total pressure at the stagnation point and reaching a minimum pressure at the top of the cylinder. (iv) The calculated results depend to a certain extent on the discretization of the computational domain and accuracy increases with increase of number of elements.