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

Authors: Moloud ArianMaram; Mahmoud Ghiasi; Hassan Ghassemi; Hamid Reza Ghafari
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, two different 3D hydrofoils with profiles NACA0012 are simulated in the potential flow. Boundary element method (BEM) and nonuniform rational B-spline (NURBS) are coupled to reduce error and increase accuracy. The 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. The pressure on a 3D hydrofoil with NACA412 profile iscompared with other existing methods. The validity of result is revealed. The accuracy of the results is acceptable. The 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 to mesh surfaces. This study demonstrates that this algorithm does mesh surface of high quality, so it can be developed to generate mesh on the submerged three-dimensional bodies .


---

## Body

## 1. Introduction

Lifting bodies are an integral part of analyzing ship performance. The 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. Therefore, solution methods have a crucial role in raising the accuracy of results and reducing computational time. There 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. The 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]. The 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].The 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. This 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–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]. They 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]. This project was performed by improving the FBM method and mass transfer cavitation model considering the maximum liquid and vapor density ratio. The 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]. The 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]. They 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. Then, 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). The 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.The 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.This 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. The first model has a constant chord length spanwise, while the second model’s chord length reduces from root to tip. Figure1 demonstrates sketch of hydrofoil in two types. The angle of attack (AoA), submergence depths (d), and flow velocity (U) are factors studied in this paper.Figure 1
Sketch of hydrofoil: (a) constant chord length; (b) varying chord length.The remainder of the paper is organized as follows. Section2 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.


## 2. Mathematical Formulation

In this paper, the pressure distribution on 3D hydrofoils has been investigated by considering an incompressible, nonviscous, and irrotational flow. The total potential can be written as follows:(1)Φx,z=Ux+∅x,z,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)∇2∅=0.Figure2 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 (SH), and wake surface (SW). Incoming uniform stream direction with velocity (U) is left-to-right.Figure 2
2D geometrical definitions for hydrofoil and physical domain.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:(3)DηDt=0,(4)Paρ+U22=Paρ+V22+gη,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:(5)∂∅∂n=−U⟶.n⟶,where n⟶ is a normal vector.The 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.(6)PTE+=PTE−,(7)ΔP=P+−P−=0,(8)Δ∅=∅++∅−=Γ.The main equation is based on Green’s theory defined as follows:(9)2π∅pH=∫H∅q∂G∂nq−Gpq∂∅∂nqds.Adding the Kutta condition in equation (9), the following equation can be derived:(10)2π∅pH=∫H∅q∂G∂nq−Gpq∂∅∂nqds+∫SWΔ∅w∂G∂nds,where p and q are field point and source point; G and ∅ are Green function and perturbation velocity potential; ∂/∂n is the first derivative by normal vector; and Δ∅w is a potential difference on the surface of wake at the trailing edge.Three-dimensional Green’s function is defined asGij=1/4πrij where rij is the distance between the field point and source point on the hydrofoil.The typical method to obtain the pressure coefficient can be written as the following equations:(11a)Cp=P−P∞ρU2/2=1−VU2,(11b)V=u2+v2+W2=∂∅∂x2+∂∅∂y2+∂∅∂z2,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 ∂∅/∂x,∂∅/∂y,and∂∅/∂z are the first derivative of potential by the surface in x, y, and z directions, respectively.The 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:(12a)Su,v=∑i=0n∑j=0mNi,puNj,qvPi,j,0≤u,v≤1,(12b)Ni,0u=1,ui≤u≤ui+1,0,otherwise,(12c)Ni,pu=u−uiui+p−uiNi,p−1u+ui+p+1−uui+p+1−ui+1Ni+1,p−1u,where Ni,pu and Nj,qv are B-spline basis functions in u and v directions. Pij is the control point, and m and n are the numbers of control points in both u and v directions. Both Ni,pu and Nj,qv are defined by p-degree and q-degree. Therefore, the first-order potential derivative by the surface has been carried out to calculate pressure coefficients.This 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.


## 3. 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.Figure 3
Comparison of the pressure coefficient distribution between results of the present method with Bal and Yeung and Bouger.Comparing results from Figure3 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.The cavitation number is an important parameter to analyze cavitation. This 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. The cavitation numberσ can be written as follows:(13)σ=p1−pv0.5∗ρ∗U2,where p1 and pv are the total static pressure and the vapor pressure of water, respectively.Figure4 shows a comparison of the hydrofoil surface’s cavity between the present method and reported data by simulations.Figure 4
Comparison of the cavity on the hydrofoil surface between the presented method and Bal numerical result.In this numerical computation, NACA 0006 hydrofoil has been analyzed atAoA=5° and d = 1 m, and flow velocity has been considered 2.35 m/s according to σ=0.25.As shown in Figure4, 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.


### 3.1. The Analysis on 3D Hydrofoil with the Constant Chord Length

The first model of study is the 3D hydrofoil with a profile of NACA 0012, and the chord length is constant from root to tip. This model has 20 sections, and each section has 67 nodes. Figure5 shows a perspective view of the first 3D hydrofoil case study.Figure 5
Perspective view of the first 3D hydrofoil with profile NACA0012.The 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.Figures6 and 7 show the potential and the pressure coefficient distribution for different AoA=0°, 5°, 10°, and 15°, U = 1 m/s, and d = 1 m.Figure 6
Angle of attack on potential (U = 1 m/s and d = 1 m).Figure 7
The pressure coefficient distribution at different AoA (U = 1 m/s and d = 1 m).It can be seen from Figures6 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.Figures8 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 of5° and depth of 1 m.Figure 8
The effect of flow velocity on potential value (d = 1 m and AoA=5°).Figure 9
The effect of flow velocity on the pressure coefficient distribution (d = 1 m and AoA=5°).Figures8 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.


### 3.2. The Analysis on 3D Hydrofoil with varying Chord Length

Figure10 shows the perspective view of the second 3D hydrofoil, which has a different shape. The 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 10
Perspective view of the second hydrofoil.Figure11 shows the pressure coefficient distribution around six sections of the model at d = 1 m, U = 2 m/s, and AoA=5°.Figure 11
The pressure coefficient distribution around different hydrofoil sections (AoA=5°, d = 1 m, and U = 2 m/s).Figure11 indicates that the increase in the length of the chord leads to the increase in Cpmin.To investigate the effect of velocity flow, depth, and AoA on the hydrofoil, the analysis was performed for different depths (d = 1, 2, and 3 m), different flow velocities (U = 0.7, 1, and 2 m/s), and other AoA (AoA=5°, 10°, and 15°). The results are presented in Figures 12–14.Figure 12
The effect of AoA on the pressure coefficient distribution in the fourth section (U = 1 m/s and d = 1 m).Figure 13
The effect of flow velocity on pressure coefficient distribution in the fourth section (AoA=5° and d = 1 m).Figure 14
The effect of submergence depths on pressure coefficient distribution in the fourth section (AoA=5° and V = 2 m/s).The relation between AoA and the pressure coefficient distribution and velocity and the pressure coefficient distribution in this model is the same as that in the first model.Presented results in this part belong to the fourth section (chord length of 0.7 m). The effect of AoA on the pressure coefficient of the fourth section has been studied atU = 1 m/s and d = 1 m. The obtained results are shown in Figure 12.Figure12 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.Figure13 shows the effect of flow velocity on the pressure coefficient distribution in the fourth section in constant AoA and depth.Results reported in Figures9 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.Figure14 shows the effect of submergence depths on the pressure coefficient distribution in the fourth section.It can be seen from Figure14 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.


## 3.1. The Analysis on 3D Hydrofoil with the Constant Chord Length

The first model of study is the 3D hydrofoil with a profile of NACA 0012, and the chord length is constant from root to tip. This model has 20 sections, and each section has 67 nodes. Figure5 shows a perspective view of the first 3D hydrofoil case study.Figure 5
Perspective view of the first 3D hydrofoil with profile NACA0012.The 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.Figures6 and 7 show the potential and the pressure coefficient distribution for different AoA=0°, 5°, 10°, and 15°, U = 1 m/s, and d = 1 m.Figure 6
Angle of attack on potential (U = 1 m/s and d = 1 m).Figure 7
The pressure coefficient distribution at different AoA (U = 1 m/s and d = 1 m).It can be seen from Figures6 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.Figures8 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 of5° and depth of 1 m.Figure 8
The effect of flow velocity on potential value (d = 1 m and AoA=5°).Figure 9
The effect of flow velocity on the pressure coefficient distribution (d = 1 m and AoA=5°).Figures8 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.


## 3.2. The Analysis on 3D Hydrofoil with varying Chord Length

Figure10 shows the perspective view of the second 3D hydrofoil, which has a different shape. The 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 10
Perspective view of the second hydrofoil.Figure11 shows the pressure coefficient distribution around six sections of the model at d = 1 m, U = 2 m/s, and AoA=5°.Figure 11
The pressure coefficient distribution around different hydrofoil sections (AoA=5°, d = 1 m, and U = 2 m/s).Figure11 indicates that the increase in the length of the chord leads to the increase in Cpmin.To investigate the effect of velocity flow, depth, and AoA on the hydrofoil, the analysis was performed for different depths (d = 1, 2, and 3 m), different flow velocities (U = 0.7, 1, and 2 m/s), and other AoA (AoA=5°, 10°, and 15°). The results are presented in Figures 12–14.Figure 12
The effect of AoA on the pressure coefficient distribution in the fourth section (U = 1 m/s and d = 1 m).Figure 13
The effect of flow velocity on pressure coefficient distribution in the fourth section (AoA=5° and d = 1 m).Figure 14
The effect of submergence depths on pressure coefficient distribution in the fourth section (AoA=5° and V = 2 m/s).The relation between AoA and the pressure coefficient distribution and velocity and the pressure coefficient distribution in this model is the same as that in the first model.Presented results in this part belong to the fourth section (chord length of 0.7 m). The effect of AoA on the pressure coefficient of the fourth section has been studied atU = 1 m/s and d = 1 m. The obtained results are shown in Figure 12.Figure12 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.Figure13 shows the effect of flow velocity on the pressure coefficient distribution in the fourth section in constant AoA and depth.Results reported in Figures9 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.Figure14 shows the effect of submergence depths on the pressure coefficient distribution in the fourth section.It can be seen from Figure14 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.


## 4. 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. The first model was simulated in velocity 0.7, 1, and 2 m/s andAoA=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.The conclusions are summarized as follows:(i)
The 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)
The 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 bodyBesides, 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 inCpmin. Therefore, 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. This 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.



---

*Source: 7252949-2021-10-25.xml*