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 , and the decrease in depth and angle of attack leads to the growing value of . 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 .
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. Figure 1 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.
The 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.
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: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.
Figure 2 indicates submergence depths of hydrofoil (d), the direction of the axis (X, Z), flow velocity (U), normal vector (), 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.
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:where , , , , and 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 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.
The 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; is the first derivative by normal vector; and is a potential difference on the surface of wake at the trailing edge.
Three-dimensional Green’s function is defined as where 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:where P is the pressure at field element, is the pressure at infinity, V and U are total velocity vector and incoming flow velocity, , and are velocity on a surface in x, y, and z directions, respectively, and 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:where and are B-spline basis functions in and directions. is the control point, and m and n are the numbers of control points in both and directions. Both and 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, , 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.
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:where and 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 and d = 1 m, and flow velocity has been considered 2.35 m/s according to .
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.
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. Figure 5 shows a perspective view of the first 3D hydrofoil case study.
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.
Figures 6 and 7 show the potential and the pressure coefficient distribution for different , , , and , U = 1 m/s, and d = 1 m.
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 , 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 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.
3.2. The Analysis on 3D Hydrofoil with varying Chord Length
Figure 10 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 11 shows the pressure coefficient distribution around six sections of the model at d = 1 m, U = 2 m/s, and .
Figure 11 indicates that the increase in the length of the chord leads to the increase in .
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 (, , and ). The results are presented in Figures 12–14.
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 at U = 1 m/s and d = 1 m. The obtained results are shown in 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 is lower than that of and , 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.
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 and , , , and 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, , , and , 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 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 . 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.
Appendix
A. Application of the Free Surface
To investigate the effect of the depth, an imaging method has been used. So, one term was added on Gij, and it was rewritten to where is the distance between the control element on the hydrofoil and the singular element on an image of the hydrofoil.where x, y, and z are coordinates of the central point of the field element, x0, y0, and z0 are coordinates of the central point of the source element on the body, and ξ, , and ζ are coordinates of central point of field element on the body image.
B. Calculation of the Velocity Potential
To calculate the velocity potential, equation (10) is discretizated, and it is rewritten as the following equation: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:
The start point of cavitation occurs when while relates to a zone of cavitation.
Symbols
| AoA: | Angle of attack |
| C: | Chord length |
| D: | Submergence depths |
| : | Gravitational acceleration |
| G: | Green function |
| L: | Lower surface of hydrofoil |
| m: | Numbers of control points in direction |
| n: | Numbers of control points in direction |
| : | Normal vector |
| : | B-spline basis function in direction |
| : | B-spline basis function in direction |
| p: | Field point |
| P: | Pressure at field element |
| : | Atmospheric pressure |
| P ij: | Control point |
| p 1: | Total static pressure |
| : | Vapor pressure of water |
| : | Pressure at infinity |
| q: | Source point |
| : | Distance between the field point and source point on the hydrofoil |
| : | Distance between the control element on the hydrofoil and the singular element on an image of the hydrofoil |
| S H: | Body surface |
| S W: | Wake surface |
| U: | Flow velocity |
| Up: | The upper surface of hydrofoil |
| : | Total velocity |
| : | Angle of attack |
| : | Circulation value |
| : | Density of water |
| : | The cavitation number |
| : | Potential difference on the surface of wake at the trailing edge |
| : | The total potential |
| ∅: | The perturbation velocity potential |
| : | The free surface deformation. |
Data Availability
Some data or models used during the study are available from the first author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
Moloud ArianMaram conceptualized the study, developed methodology, validated the study, performed formal analysis, and wrote the original draft. Mahmoud Ghiasi conceptualized the study, developed methodology, validated the study, performed formal analysis, and reviewed and edited the original draft. Hassan Ghassemi conceptualized the study, developed methodology, validated the study, and reviewed and edited the original draft. Hamid Reza Ghafari conceptualized the study, validated the study, performed formal analysis, and wrote the original draft.