Improved Hydrodynamic Analysis of 3-D Hydrofoil and Marine Propeller Using the Potential Panel Method Based on B-Spline Scheme

In this paper, the hydrodynamic performance of lift-body marine propellers and hydrofoils is analyzed using a B-spline potential-based panel method. The potential panel method, based on a combination of two singularity elements, is proposed, and a B-spline curve interpolation method is integrated with the interpolation of the corner points and collocation points to ensure accuracy and continuity of the interpolation points. The B-spline interpolation is used for the distribution of the singularity elements on a complex surface to ensure continuity of the results for the intensity of the singular points and to reduce the possibility of abrupt changes in the surface velocity potential to a certain extent. A conventional cubic spline method is also implemented as a comparison of the proposed method. The surface pressure coefficient and lift the performance of 2-D and 3-D hydrofoils of sweepback and dihedral type with different aspect ratios are analyzed to verify the rationality and feasibility of the present method. The surface pressure distribution and coefficients of thrust and torque are calculated for different marine propellers and compared with the experimental data. A parametric study on the propeller wake model was carried out. The validated results show that it is practical to improve the accuracy of hydrodynamic performance prediction using the improved potential panel method proposed.


Introduction
The panel method has been widely used as a powerful tool in the aerodynamic and hydrodynamic fields since the successful work presented by Hess and Smith in 1967 [1]. The velocity potential method transforms the Laplace equation using the Green formula, which transforms an integral problem into a surface integral problem of the object, and the velocity potential is used as the unknown quantity. Brandner [2] and Baltazar et al. [3] first proposed this basic principle and compared the results obtained using the panel method and the RANS equation method. Moreover, applications of Green's function method to the frontier science field and engineering problems include the analysis of the stress gradient of nanobeams in nanotechnology with a possible loss of symmetry of Green's function for nonlocal boundary conditions [4], to predict accurately dynamic behavior of vibroacoustic systems [5], and to solve transmission loss problems [5,6], radiation, diffraction problems [7,8], etc.
For three-dimensional objects that generate lift, such as hydrofoils and propeller blades, apart from the basic objective of providing lift, we considered that there should be some attached vortices on the blades and that the additional wake vortices are assumed at the wake line, hence, all the vortices form a closed interval or the equal strength of an infinitely long vortex line. The panel method based on the velocity potential normally requires the Dirichlet boundary condition (B.C.) to solve the governing equation matrix, and the Dirichlet boundary condition requires the virtual velocity potential in the interior region of the object to be constant. Based on the results reported by Belibassakis et al. [9], it was assumed in the present research that the internal velocity potential of a three-dimensional object is equal to the perturbation velocity potential at infinity, hence, the intensity distribution of the source on the surface of the object can be determined, and then the governing equation, i.e., the Laplace equation, for each collocation point can be solved.
For hydrofoils with a continuous slope, there is no suitable additional condition. Because the slope at the trailing edge of a two-dimensional hydrofoil is discontinuous, it is possible to find a definite loop in which the velocity at the trailing edge of the hydrofoil is infinite. This additional constraint is the Kutta condition [10], i.e., zero velocity at the trailing edge. However, this condition is not directly suitable for performing numerical calculations. To meet the Kutta condition, Kinnas and Hsin [11] applied a simplex iterative method. This method can meet the equal pressure at the trailing edge of the hydrofoil, but it requires several iterative operations to reach convergence. Another approach, the Newton-Raphson method [12], has been proposed to satisfy the Kutta condition; however, it is difficult to reach convergence by this method when the pressures at the trailing edge are equal. Therefore, Srivastava and Roychowdhury [13] presented an improved Newton-Raphson iteration method for enhanced stability and reliability of the Kutta condition in numerical calculations. With this method, the computation efficiency is greatly improved.
For a specific problem, the panel method should be used according to the specific conditions. Mantia and Dabnichki [14] and Eça and Vaz [15] studied two-dimensional wings using the panel method based on velocity and based on velocity potential, respectively. They found that the panel method based on velocity and/or velocity potential can be applied to estimate the irrotational potential flow around 3-D hydrodynamic shapes with a certain thickness. However, for thin objects, the velocity potential panel method can provide more accurate results. The influence of the source and doublet combination elements on the velocity potential of the collocation point leads to a lower-order singularity than the velocity influence coefficient.
For a continuous geometric surface, Kanemaru and Ando [16] applied flat rectangular panels that are simple and easy to calculate; however, there will always be gaps between adjacent rectangular panels. Tarafder and Suzuki [17] applied surface modelling with discretized hyperboloid panels, where the gap between the panels can be narrow. The accuracy of the discretized form directly influences the accuracy of the calculations. If the difference in the strengths of the adjacent source and doublet is relatively large, a discontinuous situation will arise, leading to inaccuracy in the velocity potential derivative.
Based on this consideration, the grid meshing method based on the B-spline scheme was adopted in this study to make the spatial interpolation of the singular point of the panel elements more accurate and generate a reasonable surface mesh. The B-spline scheme based on the Bezier scheme [18] uses fewer control vertices to generate the high-order curves and surfaces and enables curve control over local deformation without affecting the global deformation. Hence, it is considered one of the most important geometric modelling methods [19]. The B-spline scheme is widely applied in vehicle design [20,21], aerodynamic optimization [22], submarine guidance systems with path planning [23], etc.
Prediction of hydrodynamic performance of hydrofoil and marine propeller using B-spline high-order panel method has been studied in the literature [24][25][26][27]. The lower number of panels and high computational efficiency were experienced in the high-order panel methods. However, the effect of the singularity war [28,29] on numerical instability is still difficult to overcome in the high-order methods whereas the low-order panel method [2,9,11,28,30,31], i.e. the constant potential strength with an equal singularity density on a panel, is more stable to numerical computation for the analysis of a smooth hydrodynamic shape [25]. To improve the accuracy of high-or low-order panel method, the integrated trailing wake surface shall be an improved wake modelling [24,30,32].
In this study, we propose a method by integrating the B-spline geometry representation with the low-order panel method so as to obtain a smoother and continuous discretized surface integrated with an analyzed trailing wake model for the purpose of enhancing the accuracy and numerical stability in the hydrodynamic analysis. Extending the earlier research work [21], the surface pressure coefficient and lift forces on a hydrofoil and propeller were calculated and analyzed using the cubic B-spline low-order panel method.
The potential panel method based on a combination of two singularity elements was proposed and a B-spline curve interpolation method was integrated with the interpolation of the corner points and collocation points to ensure the accuracy and continuity of the interpolation points. The B-spline interpolation was used for distribution of the singularity elements on a complex surface to ensure continuity of the results of the intensity of the singular points and to reduce the possibility of abrupt changes in the surface velocity potential to a certain extent.
Control polygon and control points of the B-spline interpolation were inversely calculated by the basis function for fitting and smoothing 3-D hydrofoil and propeller discreted panels integrated with a low-order panel method for efficient and effective hydrodynamic prediction improved in this study. The hydrodynamic performance of the 3-D hydrofoil, sweepback hydrofoil, and dihedral hydrofoil with different aspect ratios was analyzed to verify the rationality and feasibility of the method. Finally, the results of pressure distribution and the coefficients of thrust and torque for different 3-D screw propellers were obtained and compared with the experimental data. The comparison results were found to be satisfactory. National Advisory Committee for Aeronautics (NACA) airfoils [33][34][35] and David Taylor Model Basin (DTMB) propellers [2,[36][37][38] have been extensively adopted for the validation of the numerical schemes thanks to the availability of the experimental data for several different foil profiles and propeller models.

Surface Panel Method
Panel methods involve numerical models based on simplified assumptions regarding inviscid and incompressible flow problems. In principle, if a problem can be solved by distributing the unknown quantities on the boundary surface surrounding a foil under simplified potential flow, the characteristic coefficients of the foil can be obtained. Under these assumptions, the velocity vector that describes the flow field can be represented as the gradient of a scalar velocity potential. 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, the continuity equation can be expressed as follows [1]:  is the Laplacian and  is the velocity potential. The solution of Equation (1) is based on the Green formula, a well-known function method, which can be derived from Gauss's law and the divergence theorem. For a finite fluid domain V, the Green formula can be obtained using the divergence theorem [9]: can represent any function in the finite field that has a continuous first-order partial derivative. The velocity field function in the fluid domain is represented in the calculation of the surface element, where S represents the outer boundary of the volume V and n is the outer normal vector of the boundary. By Green's third formula for hydrofoil and propeller problems, the velocity potential at point ) , can be expressed as [32]: where Q represents the singularity point in the field and R is the distance between the collocation point P and the singularity point Q. Hence, the basic assumptions of the potential panel method for a lifting body can be illustrated as in Figure 1. In the nonpenetrating and perturbation velocity potential attenuation conditions, the integral equation for the surface of the object can be expressed as [28,29] where   represents the distribution of the doublet of the object's wake surface The doublet strength on the wake surface can be determined by the Kutta condition and the wake surface, and the perturbation velocity potential of the object surface can then be obtained. The perturbation velocity can be obtained by differentiating the perturbation velocity potential, and then the hydrodynamic performance can be obtained. So far, the main symbolic relationship of the above-mentioned potential flow theory is shown in Figure 1.

B-Spline Model Scheme
The motivation for using B-splines is to ensure that the B-spline curved panels have at least C 2 continuity and the B-splines can be provided with local control, i.e., changing a single point does not affect the entire curve. While natural cubic splines provide that level of continuity, they are also subject to global control, i.e., changing a single point affects the entire curve.
The B-spline curve method inherits the advantages of the Bezier curve method, preserving the basic modelling features such as defining the vertices of curves and surfaces. The B-spline curve uses a set of basic functions that is different from the one used in the Bezier curve.
The B-spline curve is expressed as follows [19]: is the de Boor point, which is the control point of the curve. The polygon defined by all of the control vertices is defined as the control polygon, for instance, a control polygon covering the NACA4410 hydrofoil geometry [33] can be calculated to generate a discreted cubic B-Spline curve or surface as shown in Figure 2. In Equation (5), is a base function of order k. Every base function is a polynomial of order k, consisting of increasing vector nodes u. The base function can be expressed using the de Boor-Cox recursive formula [19,24] as follows: In practice, the calculation of the entire curve is rarely performed in the application of the B-spline curve. To calculate a series of normalized B-spline basis function from Equations (6)-(8), a B-spline function of a specific order can be obtained; some of the curves are shown in Figure 3.
, and all the control vertices constitute the control surface grid. The order of u is k, and that of v is l. The corresponding node vectors are order tensor product B-spline surface can be expressed as [19]: For the surface problem, the inverse solution process involves calculating the inverse of the tensor product. Therefore, the inverse of the surface can be expressed by 2 inverse curves. The B-spline surface equation that needs to be inversed can be expressed as follows [19]: where an equation similar to the curve equation [19] is constructed as given below: The inverse of the applied base function is used to obtain the control vertices on the isoparametric lines. The B-spline curve is defined by the resulting control vertices. The control vertices of the required curves can be solved using the control vertices obtained as the initial points.

Numerical Representation of Surface Panel Method
The numerical calculation process for a surface panel method is shown in Figure 4. First, the hydrodynamic surface is defined to generate grid coordinates and points to discretize the hydrodynamic surface into N panels. Second, the number and type of singularities (solutions of the Laplace equation), such as the source, doublet, and vortex, are chosen. Third, the collocation points are placed on the curved panels to satisfy the surface-tangency condition by using the cubic B-spline interpolation scheme. Fourth, the Kutta condition is imposed and the influence of all of the singularities and the freestream at a fixed collocation point is summed. Fifth, solving the singularity strength from the system of equations generated by step 4, i.e., the zero normal flow boundary condition on each of the collocation points resulting in the set of algebraic equations as shown in Equations (12)- (14), is carried out. Finally, the estimated performance to be improved or satisfied is determined, including pressure, velocity, and load. The related numerical representation of the panel method can be referenced in Appendix.

Modelling of 3-D Hydrofoil
The Dirichlet boundary condition is used for a hydrofoil with a certain thickness. The governing equation can be expressed as follows (the internal disturbance potential of the object is equal to zero) [28]: Equation (12) is applied to every collocation point P. In the equation, l and k are the count numbers of the singularity elements. The geometry of the hydrofoil is shown in Figure 5. The influence coefficients of the doublets k c , the wake doublets l c , and the source k b can be calculated. The Kutta condition is expressed by the relationship between the strength of the doublets on the upper and lower surfaces and the doublet strength on the wake surface [28,32]: The matrix equation under the Kutta condition can be expressed as follows: In this case, μw is replaced with μN -μ1. Therefore, the order will be reduced to N. Only the first and the Nth columns will change because of the term ±ciw. Hence, the influence of the doublet can be expressed as follows: After the known matrix multiplication is moved to the right-hand side (RHS) of this equation, the RHS vector can be obtained. Based on the theory of the panel method shown above, a program based on MATLAB was developed to calculate the strength of the singularities. Then, the coefficient of pressure on the surface of the foil can be obtained. The external potential u  can be represented as the internal potential i  plus the doublet strength μ. Then, the local external tangential velocity component at each collocation point can be expressed by the following equations: where l is the distance between 2 adjacent collocation points.

Numerical Modelling of 3-D Propeller
For geometric boundary discretization and the motion of the propeller, 2 sets of coordinate systems are established to describe the propeller geometry accurately. The rectangular coordinate system and the cylindrical coordinate system are defined as shown in Figure 6. The conversion between the 2 sets of coordinate systems is as follows [32]: The above equation expresses the coordinates of any point on the camber surface of the propeller blade. If we use to replace f in the equation, we can obtain the numerical expression for the upper and lower surfaces of the propeller blade. Figure 7 shows the 3-D propeller generation process for the surface panel method. The surface of the propeller blade is discretized using hyperbolic quadrilateral panels. Each propeller blade is divided into 2 where max s δ and 1 s δ are the spanwise dimensions of the panel with the maximum spanwise dimension and the first panel, respectively. The cosine-clustering method was adopted to optimize the sparse density of the discrete surfaces of the propeller panels.

Propeller Wake Model
The propeller wake model has an influence on the perturbation velocity potential flow of the propeller. The helical spiral surface can be used to model the wake surface of a propeller. The propeller pitch angle can be expressed by a linear combination of the blade pitch angle and that based on lift-line theory. Greeley and Kerwin [36] proposed a model of radial shrinkage in the wake region and achieved good results. In this wake model, the radius of the wake surface decreases and eventually shrinks into 2 vortices, the hub vortex and the tip vortex. Kanemaru and Ando [16] improved the geometry of the wake surface by making it more refined, combining Gaschler's experimental results [40] on the wake surface into the modelling.
In this study, a wake shrinkage model was applied to improve the prediction of the performance of a propeller in which the pitch of the wake surface changes with shrinkage ( Figure 8). To obtain a more reasonable panel distribution on the wake surface, the discrete coordinates of the wake surface can be calculated using Equations (21)- (25). Specifically, when a downstream coordinate is given, the region from the blade to the downstream coordinate w  is the near-wake region, and the region farther than the downstream coordinate w  is the far-wake region. The distance from the angular coordinate TE  of the trailing edge to the downstream coordinate w  is defined by s, and the radius of the wake is contracted to w r at 1  s . It is assumed that the wake radius of the central part changes smoothly, and that the pitch angle of the wake surface changes smoothly in the direction of the blade radius and wake flow. The radius of the far-wake region at 1  s does not shrink, while the pitch angle does not change in the direction of the wake flow. The geometric parameters of the wake surface are defined as follows: (1) Tip vortex radius [ Figure 9 shows discrete panels on a propeller and wake surface for the surface panel method. The geometric model of the wake surface can be constructed accurately using the characteristic quantities of the wake surface. Because the influence of the panels at the trailing edge is greater, the cosine-clustering method was applied for the wake surface at a given radial profile. The corner point coordinates of the panels in the wake transition region were determined using the following equation [37]:

Validation of the Panel Method Using 2-D Hydrofoil Performance
To validate the two-dimensional hydrofoil panel method, the hydrofoil profiles NACA0012, NACA2410, NACA6412, and NACA4418 [33] were selected for calculation. The lift coefficient of NACA0012 was calculated at different angles of attack using the two-dimensional hydrofoil panel method and compared with the experimental data; the results are shown in Figure 10a-d.   [33], and the average error using the B-spline panel method is within 6.3% whereas the average error was 7.1% using the cubic spline panel method. The results for NACA0012 and NACA2410 have the highest degree of agreement. There are some differences between the calculated results and experimental data for NACA6412 and NACA4418, but the trend of the pressure coefficient curve obtained through calculation is in good agreement with the experimental data [33]. The higher the degree of the camber line, the greater the effect on the calculation accuracy; the number of panels also has some influence on the calculation results. For the two-dimensional hydrofoil, the number of panels should be controlled in the range of 40-80. A lower grid density would lead to less accuracy. Figure 11 shows a comparison of the calculated and experimental results of lift coefficient for NACA0012 at different angles of attack using the B-spline and cubic spline panel methods. The calculated results using the B-spline panel method show a higher degree of agreement with the experimental data.

Validation of the Panel Method Using 3-D Hydrofoil Performance
In this study, a three-dimensional rectangular hydrofoil was selected for calculation using the panel method. The hydrofoil type selected was NACA0012 [33]. The aspect ratio was set to 10, the angle of attack was set to 5°, and the number of panels was 180. The surface pressure distribution on the hydrofoil is shown in Figure 12. The surface pressure distribution of the hydrofoil can be directly compared with the experimental data reported by Lee and Jang [34]. In the present study, the number of chordwise panels per section is 18, and the number of spanwise panels per section is 10. With this discretization method, it is possible that some regions of the panel shape become relatively slender, resulting in a lower calculation accuracy. It is necessary to develop some algorithm to judge the size of the panels to ensure that they are not very slender. The results obtained using the hydrofoil panel method are compared with the experimental results for the sections defined by a span fraction (s) to the wing span (S) as s/S = 0.25, 0.5, 0.65, and 0.85 as shown in Figure 13a

Hydrodynamic Analysis of Sweptback and Sweptforward Hydrofoils
Compared with the straight rectangular hydrofoil, the sweptback hydrofoil can effectively reduce the impact force of the flow and improve the hydrodynamic performance. To test the adaptability and stability of the panel method, rectangular hydrofoils with different sweptback and sweptforward angles were selected for analysis. The hydrodynamic performance was calculated with different sweepback angles at a limited angle of attack using the panel method. Discretization of the hydrofoil was performed using the cosine-clustering method. The aspect ratio was set to 6 and the angle of attack was set to 5°. A comparison of the calculation results and the experimental data of the sweptback hydrofoil performance [41] is given in Figure 14. The experimental data considered the ground effect on the sweptback hydrofoils whereas the prediction of the sweptback hydrofoil performance in the present panel method did not consider the ground effect, hence, the value of the experimental lift coefficient at the zero sweptback angle is greater than those of the panel methods as shown in Figure 14a.

Hydrodynamic Analysis of Dihedral Hydrofoils
The stability of the hydrofoil can be improved by using a dihedral hydrofoil as shown in Figure  16. The calculation results of the hydrodynamic performance of the dihedral hydrofoil with different aspect ratios based on the cubic B-spline panel method are shown in Figure 17.  For low-speed hydrofoils, the hydrodynamic performance can be improved by increasing the lift coefficient and the lift-to-drag ratio; in this case, it is possible to increase the hydrofoil aspect ratio, reduce the sweptback angle, and keep the hydrofoil as straight as possible. At higher speeds, the hydrofoil must not only have a good hydrodynamic performance, but it may also be resistant to stress fatigue; hence, it is not appropriate to reduce the aspect ratio and increase the sweptback angle.

Results of Propeller Performance
In this study, the DTMB series propellers with different values of rake and skew were selected for calculation of hydrodynamic performance using the panel methods with the present method and the conventional cubic spline scheme. The selected propellers were DTMB4119, DTMB4497, DTMB4498, DTMB4381, DTMB4382, DTMB4383, and DTMB4384.

Validation of the Panel Method Using DTMB4119 Propeller Performance
DTMB4119 was approved by the 20 th International Towing Tank Conference as a standard propeller to verify the accuracy of the panel method. The geometric parameters of the propeller and the experimental data pertaining to DTMB4119 were taken from the conference reports [2,37], and the corresponding details for the other propellers were taken from previously published documents [29,36,43].
The propeller blade surface was discretized in the chordwise direction using the cosine-clustering method. The DTMB4119 propeller panel discretization model is shown below (Figure 18) and the relative characteristics of the DTMB4119 are shown in Table 1 [37]. The literature [2,37] offers the experimental data of the pressure coefficient distribution for DTMB4119 at different sections at a speed factor of J = 0.833. A comparison of the pressure coefficients obtained from the experimental data and the calculated results is shown in Figure 19a-d. In general, the calculated results of the surface pressure coefficient of the DTMB4119 propeller blade are in agreement with the experimental results [2,37], and the average error is within 7.3%. The selection of the wake model has some influence on the calculation results. The model selected in this study is the empirical wake model, which is different from the actual wake. Therefore, the wake model is not the main factor in the calculation results.

Results of Propeller Thrust and Toque
A comparison of the hydrodynamic coefficients for the propellers obtained from the calculation results and experimental results [2,36,38] is shown in Figure 20a-f. The number of panels is 280. Based on the cosine clustering method together with cubic B-spline interpolation, these panels were clustered more tightly and smoothly near leading and trailing edges of 3-D propeller hydrofoils placed on the propeller helix lines for the best results. Figure 20a-f show that the hydrodynamic performance of the propellers is in good agreement with the experimental data, and the average error is within 7.3% based on the cubic B-spline panel method whereas the average error 8.2% using the conventional cubic spline panel method. Especially for the DTMB4119 propeller, the accuracy of the hydrodynamic performance calculation is high, and the average error is within 6% based on the proposed method.

Conclusions
A computational technique using the B-spline panel method and the cubic spline panel method to predict the performance of a hydrofoil and propeller is presented in this paper. The surface pressure coefficient and lift performance of a two-dimensional hydrofoil are calculated and analyzed by the above numerical methods. The results are compared with the experimental data. In addition, the hydrodynamic performance of a three-dimensional hydrofoil and of a sweepback hydrofoil and a dihedral hydrofoil with different aspect ratios are analyzed, and the average error is within 4% in the present method, whereas it is 5.9% of the cubic spline method. The results confirm the rationality and feasibility of the method. The surface pressure distribution and the coefficients of thrust and torque of different propellers are calculated and compared with the experimental data, and the average error is within 7.3% in the present method whereas it is 8.2% using the conventional method. The comparison results are found satisfactory. In conclusion, the validity of all the calculated results from the cases studied in this paper was proved by comparing them to the available studies in the literature. Additionally, the calculated results using the proposed B-spline panel method show a higher degree of agreement with the experimental data than the cubic spline one. In the future, the panel method can be used for the inverse design based on the required load and modelling predictions of hydrofoil cavitation.
where ij  represents the Kronecker number and N represents the total number of panels on a single propeller blade.
The tangent vector of the panels can be expressed as The unit normal vector of the panels can be expressed as a a n a a (A6) The surface of the panels can be expressed as follows: The integral of the influence coefficient can be expressed as follows [24]: The integral function The following equation is obtained as follows [24]: The influence coefficients ij C and ij B can be expressed as follows [25]:  C  I  I  I  I  B  I  I  I where the influence coefficients can be expressed as follows [25]: R R a a R a n R a n R a a R a a R a R a R n R R a a R Q P , where C n is a normal vector on the collocation point and P is the coordinate of the collocation point. Thus, the governing equation can be solved.
A2. The Geometric Parameters of the Propeller for Surface Panel Code (see Table A1-A5)