Prediction and analysis of 3D hydrofoil and propeller under potential flow using panel method

Potential flow over an airfoil plays an important historical role in the theory of airfoil. The governing equation for potential flow is Laplace’s equation. One of Green’s identities can be used to write a solution to Laplace’s equation as a boundary integral. Using distributions of singularity solutions and determining their strength via the boundary conditions is the essence of panel method. This paper introduces a quick prediction method of three-dimensional hydrofoil and propeller performance based on panel method. The surface of hydrofoil and propeller is divided into numbers of quadrilateral panels. Combined sources with doublets singularities will be distributed on the corners of panels. Calculated blade pressure distributions of hydrofoil and propeller agree well with experimental data. Several sample calculations have been included using panel method.


Introduction
A singularity panel method based on Green's function integral equation is used to calculate the motion characters. The motion and general behavior of a fluid is governed by the fundamental laws of classical mechanics and thermodynamics and plays an important role in many fields. An introductory text on fluid mechanics, such as the study of Sabersky [1], surveys the basic concepts of fluid dynamics and the various mathematical models used to describe fluid flow under different restrictive assumptions. After the first successful application of numerical approach to the solution of flow around an arbitrary body by Hess and Smith [2] in the 1960's the singularity panel method has developed quickly.
Panel method was mainly applied to solving lifting flows for thick airplane, but now panel code can be also used in shipbuilding field [3] and even in submarine [4]. Now panel method code is suitable for marine propellers [5] and numerical resistance calculation of hulls [6]. The recent trend is the application of low-order methods [7], and many different panel methods are used now. This led to some comparisons [8] about different panel methods and other methods. Apparently, low-order panel methods are cheaper and faster to use. The panel method which is proposed by this paper is able to predict the three dimensional hydrofoil and propeller performance without performing a computational simulation. Generally, for performing the computational simulations, a high performance computer and a lot of computational time would be needed [9], [10]. However, using this method can reduce the computational cost.
In this paper, panel method is applied to the prediction of three-dimensional hydrofoil and propeller performance. Surface pressure distributions of hydrofoil and propeller will be calculated to compare with experimental data.

Solution of potential flow using panel method
Panel methods are numerical models based on simplified assumptions about inviscid and incompressible flow problem. In principle, if a problem can be solved by distributing unknown quantities on the boundary surface surrounding the body under simplified potential flow, characteristic coefficients of foil will be obtained. Under these assumptions, the vector velocity describing the flow field can be represented as the gradient of a scalar velocity potential. Therefore, the equation is as following: (1) A statement of conservation of mass in the flow field leads to Laplace's equation as the governing equation for the velocity potential. If the flow in the fluid region is considered to be incompressible and irrotational then the continuity equation reduces to Eq. (2). 0 2 I (2) ∇Φ is indicated in the reference that is attached to the body. Apparently, the disturbance created by the motion must disappear far from the body as Eq. (3).   I  I  I  I  I  I  I  I This equation is one of Green's identities. The integral of surface is all the boundaries S. The body is enclosed in a volume V, and Φ 1 , Φ 2 are two scalar functions of position. It is well known that any element which satisfies Laplace's equation can be written as an integral over the boundary surface S. In Eq. (4) S is the surface integral, and S=S B +S W +S ∞ where S B means body model boundaries, S W means wake model boundaries and S ∞ means outer boundary.
One important observation is that Eq. (2) can be solved by distributing elementary solutions on the boundaries. Commonly singular solutions for panel methods includes vortex, source and doublet distributions. Now the solution is reduced to find the appropriate element strength over the known boundaries. When the potential is specified on the boundaries this mathematical problem is called the Dirichlet problem [11]. Another approach to the solution is to describe that the normal flow boundary condition on the solid boundaries is zero. This problem is named as the Neumann problem [11]. For a given boundary condition, the above solution can be used in specific condition respectively.

Two-dimensional foil performance in panel method
Based on the distribution of singularity elements, surface geometry, and the exact type of boundary conditions, computational methods can be constructed. What should be notified is by the combination of some types of singularities the panel method can be more accurate. Further improvement is obtained by changing spacing and density of grid, wake model, location of collocation points and singularity elements, method of enforcing the boundary conditions, and Kutta condition [12].
A step sequence ( Fig. 1) is specified when a numerical solution of panel method is established. The panel code is developed following these steps.

Figure 1 Steps toward constructing a numerical solution
In this paper we use the combination of source and doublet singularity elements on the surface of the foil. It means that every panel has a local source and doublet strength and either the source or the doublet strength values should be specified. Also, any doublet distribution can be replaced by a vortex distribution, and solutions based on the vortex can be used in this case too [13].
In the process of grid generation, the singularity elements' corner points and collocation points should be defined. Singularity elements are placed at corner points and the boundary conditions, such as zero normal flow on the solid surface, will be enforced at collocation points. Fig. 2 shows the distribution of collocation points and singularity elements on panels.

Figure 2 Distribution of singularity points and collocation points on foil
To represent a numerical solution the surface S is divided into N panels and the integration is applied to every panel, and the general solution to Eq. (2) can be obtained by several sources σ and doublets μ. Here we apply Dirichlet boundary condition to this case. Then we can obtain the integral equation : Here the equation is specified at each collocation point, so we obtain a linear algebraic equation for each point. As mentioned above, the steps of establishing such numerical solution are as follows.
The potential at point P can be included into two subroutines in MATLAB that calculate the potential at point (x, z) due to the doublet and source element j: As the N+1 panel corner points and N collocation points are generated, Dirichlet boundary condition will be used and thus the collocation points should be placed inside the body. Then the influence of the potential at collocation point i due to element j is obtainable.
The influence coefficients will result in a N×N matrix, but the number of unknown elements are N+1 because of the wake doublet μ w . Therefore additional equation is necessary by applying Kutta condition: 0 ) ( 1 w N P P P (10) In this case we replace μ w with μ N -μ 1 . Therefore the order will be reduced to N. Only the first and the Nth columns should change because of the term ±c iw . So we can replace the doublet influence like this: After the known matrix multiplication is moved to the right-hand side of this equation, we can obtain RHS vector.
Based on the theory of panel method above a MATLAB program is developed to calculate the strength of singularities, and then we can obtain the coefficient of pressure on the surface of the foil. The external potential Φ u can be represented as the internal potential Φ i plus the doublet strength μ. Then the local external tangential velocity component of each collocation point can be represented as Eq. (13) where l is the distance between the two adjacent collocation points. After that the pressure coefficient of jth panel can be calculated in Eq. (14).
For the lack of wind-tunnel data, the calculated result is compared with the pressure distribution derived from Profili. Profili is one of the most famous tools for aeromodel designers and it offers a large amount of airfoils database and explicit coefficients of pressure.   Figure 4 we can see the coefficients of pressure on upper surface and lower surface specifically. Although it shows that the present panel method is suitable for calculating the coefficient of pressure on the surface of foil, still unconformity exsits when x /c is at 0.4 and 0.6 in Fig. 3 and Fig. 4 because of divided panel size and viscosity condition. As a whole this method is suitable for calculating the coefficient of pressure on the surface of foil.
A further step can be taken to do some secondary computations. After the C p is obtained, coefficients of lift and drag are obtainable. The secondary computation can be illustrated by the case of NACA0012.
(a) (b) Figure 5. Coefficients of lift and drag on NACA0012 at different attack angles Fig. 5(a) represents the coefficients of lift on NACA0012 with the increase of attack angle, and it's obvious that the coefficients of lift on NACA0012 increases with the increase of attack angle. From Fig.8 we can see that the C l reaches peak at 38 degrees, and then it decreases. Figure 5(b) shows the coefficients of drag on NACA0012 with increase of attack angle. Therefore, this program based on panel method provides a solution to potential flow around the surface of foil, and it needs less computational time and provide a quick estimation of performance of foil. But it is a sensitive program. The mismatching between calculated data and Profili data depends on the grid, the number of panels and other various parameters. These deficiencies can be improved in the future.

Prediction of hydrofoil and propeller performance 4.1 Geometry and wake models
Three-dimensional numerical solution based on panel method is similar to two-dimensional solution. In attempting to construct a potential-flow problem of threedimensional case, only the wake and trailing conditions (Kutta condition) will need more additional attention [14]. In this case, the most difficult part is the modeling of geometry, and some solutions to arbitrary geometry need to be sought [15]. The geometries of hydrofoil and propeller are shown in Fig. 6.
The key point of this step is the definition of the propeller coordinate. The cartesian and cylindrical coordinates of a point on the camber surface with radial coordinate r and chordwise coordinate s can be expressed in terms of skew(θ ), rake(x m ), chord (c), and camber(f) [16]. With this definition, the coordinate of a point on the camber surface of the key blade may be written as°°°°°°®  13 13 Trailing vortex leaves the trailing edge of the blade and flows into the slipstream with the local veocity at that position. However, the wake surface is usually approximated by prescribed helical surface in order to avoid time consuming calculation [17]. In this paper, we will use the wake model based on QCM [18]. Namely, the wake surface leaves the trailing edge in the direction tangent to the mean camber surface in order to satisfy the Kutta condition. Then the pitch of wake changes linearly. In the down faraway stream, the wake becomes helical surface with the ultimate constant pitch.

Results and analysis
In this case, the basic panel singularity element has a constant strength source or doublet which is similar to two-dimensional case, and surface of panel is also planar. Following the previous equation, the Dirichlet boundary condition on a thick body can state that the perturbation potential inside is zero. Therefore this equation should be evaluated for every collocation point inside. l and k are the counters of panels. The influence coefficients C k , C l of the body and wake doublets and B k of the sources can be calculated. Here the zero normal velocity boundary condition should be applied, and the equations on the collocation points of panel i: In this three-dimensional panel method the shape of body is divided into panel elements. The grid is usually constructed of rectangular patches, and in this case the collocation points may be slightly under the surface. The counter k for every panel is assigned and the outward normal vector n k is identified in Figure 7.
Therefore, this equation can be specified to: Here A 1k =C 1k if no wake flows from the panel and A 1k =C 1k ±C 1p if it is shedding a wake panel. Where the N linear equations for the N unknown μ k and σ k is already known. In this case, note that a kk =1/2, except when the panel is at the trailing edge. Therefore, the right-handside matrix can be specified for the strengths of the sources are already known.
In this paper NACA0012 is selected as the root and tip geometry of the three-dimensional hydrofoil and the surface of the hydrofoil is divided into 72 panels. The influence of the wing's image can be added if needed.
Here we set that the length of chord is 1 m, the length of semispan is 10 m, and the angle of attack is 5 degrees, and then the pressure coefficient distribution of the hydrofoil is obtainable.  In this paper we will calculate DTRC4119 propeller as example. Principal characteristics of propeller 4119 are shown in Table 1. The propeller has 3 blades with no skew and rake. The blade sections are designed with NACA 66 modified profiles and a camber line a=0. 8. The advanced ratio is j=0.833. The comparison of chordwise pressure distributions at r/R=0.3 and r /R=0.7 between calculated data and experimental data is shown in Fig. 9.
Comparison of K T , 10*K Q and η between calculated data and experimental data is shown in Fig. 10.
(a) (b) Figure 9. Comparison of chordwise pressure distributions between calculated data and experimental data Figure 10. Comparison of K T , 10*K Q and η between calculated data and experimental data

Conclusion
In this paper, a potential-based panel method has been used to make prediction of hydrofoil and propeller performance. A panel method program was developed in MATLAB to calculate the pressure coefficients on the surface. This panel method program is certainly less time-consuming than CFD programs and even cheaper and easier to operate with more controlled variables. Fairly good agreements are obtained, which shows the present method is effective and accurate. However, some mismatching can be still found in some cases especially when the number of panels changes. Although many inputs can be set in panel code, the geometry of other creative shapes are still difficult to be modeled in panel code. Future work includes refinement of the wake model, generalization to the unsteady case, and interaction of propeller and rudder.