Hydrodynamic Analysis of 3D Hydrofoil Using Nonuniform Rational B-Spline and Boundary Element Method

In this paper, two different 3Dhydrofoils with profiles NACA0012 are simulated in the potential flow. Boundary elementmethod (BEM) and nonuniform rational B-spline (NURBS) are coupled to reduce error and increase accuracy. )e computer code is developed in different submergence depths (d), flow velocities (U), and various angles of attack (AoA), and the pressure is obtained by NURBS formulation. )e pressure on a 3D hydrofoil with NACA412 profile iscompared with other existing methods. )e validity of result is revealed. )e accuracy of the results is acceptable. )e competition of the two models’ results indicates that the increasing chord length leads to increase in Cpmin, and the decrease in depth and angle of attack leads to the growing value of Cpmin. Moreover, when the flow velocity is changed, the changes of potential and pressure coefficient distribution do not follow the specific trend. NURBS is a basic equation in different CAD packages because it is able tomesh surfaces.)is study demonstrates that this algorithm doesmesh surface of high quality, so it can be developed to generate mesh on the submerged three-dimensional bodies .


Introduction
Lifting bodies are an integral part of analyzing ship performance. e effect of free surface and flow velocity investigation on these bodies is an interesting issue. 3D hydrofoils are one of the most common lifting bodies in the hydrodynamics field. Many factors such as profile, chord length, span size, and angle of attack affect hydrodynamic performance. erefore, solution methods have a crucial role in raising the accuracy of results and reducing computational time. ere are several numerical methods in marine hydrodynamics; one of the most common is the boundary element method (BEM), which is solved based on Green's theory. e boundary element method is a beneficial tool for meshing the objects with complex geometries, which increases the accuracy of computations.
In this paper, the boundary element method is combined with the nonuniform rational B-spline (NURBS) modeling to reduce computational time and error. NURBS modeling is based on basic object modeling and can model the objects at a high level of quality, which increases computational accuracy.
One of the earliest studies in this field is by Giesing and Smith, who studied a 2D hydrofoil. In this investigation, the Kelvin wave source was given out on the body surface, and the linearized free surface condition was satisfied [1]. e results of this study were an integral equation that was attained by using the Neumann condition. Comparison of the boundary element method (BEM) and the finite element for two-dimensional hydrofoil was carried out by Wu and Eatock Taylor [2]. e tandem hydrofoils in different submergence depths using Reynolds-averaged Navier-Stokes (RANS) equations were investigated by ArianMaram et al. [3]. Yeung and Bouger [4] introduced a hybrid method by applying Green's theorem, which overcomes the steady-state wave-resistance.
is method is a practical tool to study lifting and nonlifting 2D bodies. Also, the effect of the nonlinear free surface was carried out by Refs. [5][6][7].
Zhang et al. [8] divided the computational domain into subdomains of viscosity and potential; they used the hybrid surface method and Euler-Lagrange hybrid method for free surface tracking that the practical advantage of this method is to avoid the singularity caused by the collision of liquids.
In addition to the potential flow-based methods mentioned above, the finite volume method can achieve free surface adsorption.
Another research in this field was addressed by Kostas et al. [9]. ey studied two-dimensional hydrofoil by combining the boundary element methods and the NURBS modeling on the asymmetric NACA 4412 in the submerged state. In their works, the genetic algorithm approach has been used to control the methods and achieve optimal conditions. Zhang et al. carried out a numerical investigation on cavitation around a 3D hydrofoil [10].
is project was performed by improving the FBM method and mass transfer cavitation model considering the maximum liquid and vapor density ratio. e numerical results in this project show that the predicted patterns and developments for cavities are well-compared. Ghiasi and Islam computed the total potential, tangential velocity, and pressure coefficient around a sphere and Wigley hull using the NURBS and Gaussian points, which were attained from NURBS [11]. e results show that this method has good accuracy even when the number of Gaussian points is low and the order of primary function in the NURBS method is low. Zang et al. coupled the isogeometric boundary element method (IGA BEM) with the nonuniform rational B-splines (NURBS) and boundary element method (BEM) with NURBS for solving the axisymmetric tanks with porous baffles in inviscid, irrotational, and incompressible flow [12]. ey investigated using NURBS for defining the geometries of the domain with lower error as possible. Several parameters were investigated to analyze the accuracy of this coupling method, whose results were utterly satisfying. Also, Sun et al. [13] combined the isogeometric boundary element method (IGABEM) and Bézier extraction of nonuniform rational B-splines (NURBS) for studying a crack in two-dimension infinite isotropic. en, they implemented their investigation because of exactly described boundaries and the crack, especially. Finally, their result again demonstratedwhen IGBEM and NURBS were coupled,they modified the accuracy of results. Abbasnia and Ghiasi [14] used the NURBS formulation to simulate a fully nonlinear wave, and they also improved a high-order fully nonlinear 2D numerical wave tank (NWT).
e Laplace equation was solved by Green theory, and then NURBS was used to simulate free surface boundary and nodes modification with the lowest error possible. So, various B-spline basis functions were employed for computing the volume, momentum, and energy conservation. Also, two damping zones were added in NWT for declining the reflected wave energy. Results indicated that the increase in the order of NURBS improves the outcomes of the volume conservation. In general, results show that adding NURBS to the solutions process is satisfactory.
e review of previous works shows that submerged and three-dimensional objects did not draw much attention.
Moreover, the mesh surface on the leading and trailing edge of 2D hydrofoils was faced with errors. Consequently, in this paper, submerged 3D hydrofoils are chosen. Because this kind of an object has different applications in hydromechanics problems, then an algorithm, which is able to calculate pressure with a lower error, can be useful to predict body performance. Hence, the NURBS formulation is used to mesh and calculate pressure, especially on areas, in which points are much closer to each other, such as middle of upper on lower surface and the leading and trailing edge, because this formulation obtains an accurate result at an acceptable time.
is paper's primary focus is the coupled BEM and NURBS to investigate some factors on pressure distribution. NURBS formulations are used to increase the accuracy of results, especially at the leading and trailing edge. NACA0012 profile was used as a model for simulating two types of a 3D hydrofoil. e first model has a constant chord length spanwise, while the second model's chord length reduces from root to tip. Figure 1 demonstrates sketch of hydrofoil in two types. e angle of attack (AoA), submergence depths (d), and flow velocity (U) are factors studied in this paper. e remainder of the paper is organized as follows. Section 2 presents the mathematical formulation to describe the governing equations and boundary conditions, and then the governing equations are developed to calculate the pressure and hydrodynamic coefficients. Section 3 presents and discusses the results of two 3D hydrofoils with the constant and varying chord length in different conditions, while summary and conclusion are drawn in Section 4.

Mathematical Formulation
In this paper, the pressure distribution on 3D hydrofoils has been investigated by considering an incompressible, nonviscous, and irrotational flow. e total potential can be written as follows: where Φ is the total potential, U is the flow velocity, and ∅ is the perturbation velocity potential. Also, the governing equation is the Laplace equation which must be satisfied by the perturbation velocity potential and the flow velocity everywhere in the fluid domain.
(2) Figure 2 indicates submergence depths of hydrofoil (d), the direction of the axis (X, Z), flow velocity (U), normal vector ( n → ), angle of attack (α), and chord length (C). Also, it introduces free surface, body surface (S H ), and wake surface (S W ). Incoming uniform stream direction with velocity (U) is left-to-right.
To satisfy equation (2), kinematic and dynamic boundary conditions on the free surface and kinematic boundary conditions on the body have been used, and the general equation can be written as follows: 2 Mathematical Problems in Engineering where Pa, ρ, V, η, and g are the atmospheric pressure, density, total velocity, free surface deformation, and gravitational acceleration, respectively. Based on the potential velocity, the kinematic condition on the body surface is expressed as follows: where n → is a normal vector.
e Kutta condition is one of the most important conditions applied at the trailing edge of a hydrofoil. It puts equal pressure on the upper and lower surfaces at the trailing edge. Hence, the circulation with constant strength value has been assumed around the hydrofoil.
e main equation is based on Green's theory defined as follows: Adding the Kutta condition in equation (9), the following equation can be derived: where p and q are field point and source point; G and ∅ are Green function and perturbation velocity potential; z/zn is the first derivative by normal vector; and Δ∅ w is a potential difference on the surface of wake at the trailing edge.
ree-dimensional Green's function is defined as G ij � 1/4πr ij where r ij is the distance between the field point and source point on the hydrofoil. e typical method to obtain the pressure coefficient can be written as the following equations: where P is the pressure at field element, P ∞ is the pressure at infinity, V and U are total velocity vector and incoming flow velocity, u, v, and w are velocity on a surface in x, y, and z directions, respectively, and z∅/zx, z∅/zy, and z∅/zz are the first derivative of potential by the surface in x, y, and z directions, respectively. e typical prior methods to compute the first derivative have an error at two edges of the hydrofoil. Hence, in this present work, the NURBS equation has been used to solve this problem.
Initially, the NURBS equation has been used to define potential on the body surface, in which the general equation for the definition of the surface is as follows: where N i,p (u) and N j,q (v) are B-spline basis functions in u and v directions. P ij is the control point, and m and n are the numbers of control points in both u and v directions. Both N i,p (u) and N j,q (v) are defined by p-degree and q-degree. erefore, the first-order potential derivative by the surface has been carried out to calculate pressure coefficients.
is method helps to improve results and to reduce errors for computing velocity on the surface. In final, with the use of equations (12a) and (11b), the coefficient of pressure will be counted.

Numerical Results
To demonstrate the validity of the results, the pressure coefficient distribution in the middle of the 3D hydrofoil was compared with the 2D-hydrofoil result reported by Refs. [4,15], which is shown in Figure 3. Flow velocity is U � 3.13 m/s, AoA � 5°, and d � 1 m.
Comparing results from Figure 3 indicates that the accuracy of the present method for hydrofoils is acceptable, particularly at the leading and trailing edge.
For further investigation, the subsequent validation has been carried out on a 3D hydrofoil with the profile of NACA 0006; in this case study, "sections" chord length is constant from root to tip, and validation has been performed based on cavitation on the hydrofoil surface. e cavitation number is an important parameter to analyze cavitation. is is because it is used to calculate the start point and zone of cavitation.
It is directly related to the difference between local pressure (vapor pressure of water) and total static pressure. e cavitation number σ can be written as follows: where p 1 and p v are the total static pressure and the vapor pressure of water, respectively. Figure 4 shows a comparison of the hydrofoil surface's cavity between the present method and reported data by simulations.
In this numerical computation, NACA 0006 hydrofoil has been analyzed at AoA � 5°and d � 1 m, and flow velocity has been considered 2.35 m/s according to σ � 0.25.
As shown in Figure 4, the accuracy of the presented result of the hydrofoil surface cavity comparing to the Bal measurement result is acceptable. It can be seen from Figure 4 that the cavity has not occurred at the leading edge while the cavity has happened near the trailing edge and also has been grown. Additionally, it was found that the cavity has been made on 80% of the hydrofoil surface.

e Analysis on 3D Hydrofoil with the Constant Chord
Length.
e first model of study is the 3D hydrofoil with a profile of NACA 0012, and the chord length is constant from root to tip. is model has 20 sections, and each section has 67 nodes. Figure 5 shows a perspective view of the first 3D hydrofoil case study. e effect of flow velocity and angle of attack on the potential and the pressure coefficient distribution has been analyzed in this model. In this subsection, all results presented are from the middle section of the model. It can be seen from Figures 6 and 7 that the potential value and the pressure coefficient distribution have similar behavior versus AoA. In the case of AoA � 0°, the potential and the pressure coefficient distribution are equal on the upper and lower surface, while the increase in angle of attack leads to the reduce in minimum of pressure coefficient distribution. Figures 8 and 9 demonstrate the effects of flow velocity on potential and the pressure coefficient distribution value for NACA 0012 sectional geometry, which have been investigated at AoA of 5°and depth of 1 m. Figures 8 and 9 demonstrate that the velocity change does not especially influence potential value and the pressure coefficient distribution. Additionally, in these simulations, cavitation was controlled, but cavitation has not happened for the first model in any angles and speeds.   Mathematical Problems in Engineering

3.2.
e Analysis on 3D Hydrofoil with varying Chord Length. Figure 10 shows the perspective view of the second 3D hydrofoil, which has a different shape. e second model has six sections so that "sections" chord length changes from 0.5 m to 1 m, and each section has 100 nodes. Figure 11 shows the pressure coefficient distribution around six sections of the model at d � 1 m, U � 2 m/s, and AoA � 5°. Figure 11 indicates that the increase in the length of the chord leads to the increase in C p min .
To investigate the effect of velocity flow, depth, and AoA on the hydrofoil, the analysis was performed for different      Figure 12. Figure 12 shows that increasing AoA has a direct impact on the pressure coefficient distribution which means the value of minimum of pressure coefficient at AoA � 15°is lower than that of AoA � 10°and AoA � 5°, respectively. Figure 13 shows the effect of flow velocity on the pressure coefficient distribution in the fourth section in constant AoA and depth.
Results reported in Figures 9 and 13indicate that increase in flow velocity does not influence on increasing pressure coefficient distribution. Value of minimum of pressure coefficient at U = 0.7, 1, and 2 m/s is −2.09, −1.39, and −1.26, respectively. Figure 14 shows the effect of submergence depths on the pressure coefficient distribution in the fourth section.
It can be seen from Figure 14 that by increasing submergence depths, minimum of pressure coefficient distribution around the hydrofoil surface decreases from −1.26 (at d � 1 m) to −1.55 (at d � 2) and −1.84 (at d � 3 m).
It should be noted that the cavitation does not occur in any of the cases studied as same as the first model.

Summary and Conclusion
In this paper, the boundary element method (BEM) and nonuniform rational B-spline (NURBS) were used; both methods are widely used to solve hydrodynamic problems.
Two different types of hydrofoils with NACA0012 profiles were applied in this study. In the first case study, chord length was constant in the span's direction, whereas in the second case study, chord length decreased from 1 m to 0.5 m in a spanwise order toward the wing tip.   BEM and NURBS methods were employed to investigate the effect of submergence depths, angle of attack, and flow velocity on the potential and the pressure coefficient distribution. e first model was simulated in velocity 0.7, 1, and 2 m/s and AoA � 0°, 5°, 10°, and 15°to compute the potential value and pressure coefficient distribution. Furthermore, the second model was simulated in the submergence depths of 1, 2, and 3 m, AoA � 5°, 10°, and 15°, and U � 0.7, 1, and 2 (m/s) for computing pressure distribution. e conclusions are summarized as follows: (i) e potential value and pressure coefficient distribution have similar behavior in the change of flow velocity and angle of attack; in other words, both increase or decline (ii) e increasing angle of attack causes decreasing minimum of pressure coefficient distribution around the body (iii) Increasing the submergence depths leads to reduce minimum of pressure coefficient distribution on the body Besides, the second model's results show that changing each "sections" length causes the change of pressure coefficient distribution on each section; in other words, increasing the length chord causes reduction in C p min . erefore, for the first model, midsection results were presented, and for the second model, the fourth section was chosen as an example.
In addition, high quality of mesh has been generated on middle of upper on lower surface and the leading and trailing edge of hydrofoils to reduce error in calculated pressure in these areas by NURBS formulation. Results of this algorithm reveal that although NURBS is a useful tool to generate the high quality of mesh surface, it has a limitation.
is method is not a suitable tool to generate mesh for the body with complex geometry because it requires more time to analyze the complex geometries.

Symbols
AoA: Angle of attack C: Chord length D: Submergence depths g: Gravitational acceleration G: Green function L: Lower surface of hydrofoil m: Numbers of control points in u direction n: Numbers of control points in v direction where i and j are defined on the central point of field and source element, respectively, and k is defined on the wake surface's element, which is created at trailing edge.

C. Calculation of Cp
To compute pressure coefficient, the derivative of the NURBS surface equations can be used as follows:

D. Study of Cavitation
To study cavitation, equation (13) is rewritten as the following equation: e start point of cavitation occurs when σ � −C P min while σ > − C P min relates to a zone of cavitation.

Data Availability
Some data or models used during the study are available from the first author upon request.