WAVE FIELD GENERATED BY FINITE-SPAN HYDROFOILS OPERATING BENEATH A FREE SURFACE

Summary The present paper focuses on the numerical investigation of the flow around the fully submerged 2D and 3D hydrofoils operating close to a free surface. Iterative boundary element method is implemented to predict the flow field. This study aims to investigate the aspect ratio effect on the free surface interactions and hydrodynamic performance of the hydrofoils under a free surface by using potential flow theory. Three different submergence depths and aspect ratios are studied in the wide range of Froude Numbers. In 3D cases, spanwise width of the numerical wave tank model is selected both equal and wider to the foil span, to observe the sidewall effects. Wave field seems to be two dimensional at low Froude numbers. On the other hand, signs of three dimensionalities are observed on the free surface structure for higher Fn, even the predicted wave elevations are very close to 2D calculations in the midsection. Increment in the Fn give a rise to the amplitude of the generated waves first, however a further increase in Fn has a lowering effect with the beginning of waves spill in the spanwise direction in the form of Kelvin waves. Free surface proximity and resultant wave field are also seeming to be linked with the lift force on the hydrofoil. As aspect ratio of the foil increase, 3D lift values are getting closer to those of 2D calculations. However, it is seen that, 3D BEM predictions of a hydrofoil under free surface effect cannot be considered two-dimensional even the aspect ratio is equal to 8.


Introduction
Using hydrofoils to generate additional lift is a common practice in marine applications. Submerged hydrofoils attached to the moving body produce lift and reduce the displacement, as well the resistance of the body, which means reducing the required energy to drive the hull. However, the dynamic mechanism of shallow submergence is quite different than unbounded flow, due to the existence of the free surface. The effect of free surface and its proximity on the performance of the hydrofoil is significant [1][2][3] and should be taken into consideration in the design stage.
Shallow submergence alters the flow field around the hydrofoil in term of pressure balance. Hydrofoils generate lift with the lower net pressure on its suction side. While the foil approaches through the free surface, it creates a wave field. This wave field is the result of the lower pressure zone at the upper side of the foil: lower pressure between the foil and the surface causes the free surface to deform, reduces the suction and therefore, reduces the lift [4]. In addition, generating a wave field unavoidably produces a drag force. On the other hand, Filippas and Belibassakis [5] showed that, free surface proximity acts upon the foil performance alternatively, and in some conditions (operating velocity and submergence ratio) the lift force may increase.
Influence of the free surface on the foils moving close to it has been subject of many the researchers in hydrodynamics field. In his well-known experimental work, Duncan [6] evaluated the wave elevations and the drag resistance of a NACA0012 hydrofoil operating beneath the free surface for various depths and three different foil speeds. Tank width in the experiment is almost equal to the foil span. Thus, a periodic wave structure (no wave spill) is obtained, and the foil was approximated as two-dimensional. They used the free surface deformation measurements to make a relation with the breaking and non-breaking wave resistance act on the foil. Duncan's experimental results are numerically validated [7] and extended for different submergence depths and a range of Froude Numbers [8,9] using 2D potential flow theory. Both numerical works point out the growing influence of free surface on lift force with decreasing free surface proximity. Di Mascio et al. [10] applied the singlephase level set method to simulate the viscous flow field around 2D NACA0012 hydrofoil at a constant submergence depth. Chen [11] conducted the numerical simulation of 2-D NACA4412 hydrofoil submerged under the free surface with a new developed vortex based panel method. They concluded that energy dissipation consideration in their method provides better agreement for numerical results with experimental ones. Karim et al. [12] addressed the 2D finite volume method (FVM) to solve the RANS equations in the viscous flow around a hydrofoil moving near a free surface. They have tested their numerical approach with Duncan's case than using seven different submerge ratios at a single velocity to capture the effect of free surface on the submerged NACA0015 hydrofoil. Prasad et al. [13] presented a numerical simulation of 2D unsteady incompressible viscous flow around a hydrofoil under the free surface at different submergence depth ratios with using OpenFOAM. They have compared the performance of the hydrofoil for different submergence ratio-Froude Number pairs to investigate the lift and drag coefficients. Wu and Chen [14] analyzed viscous uniform flow past a 2D partially cavitating hydrofoil placed at a final depth near the free surface by using well-known  − turbulence model and FVM method. They have tested the accuracy of their procedure to capture the cavitation behaviour near the free surface. Also, a detailed literature survey about cavitation near free surface can be found in their study. Hoque et al. [15] solved 2D RANS equations to analyze the free surface effect on shallowly submerged NACA4412 hydrofoil. Their results revealed that, for submergence ratio equal to 5 or more, the effect of the impact of the free surface almost vanishes and the shallow the deep water condition can be considered, for constant Froude Number. Xu et al. [16] studied the effect of free surface on the performance of 2D hydrofoils in tandem using boundary element method and revealed the dynamics of hydrodynamic force variation on the hydrofoils. Some other researchers also applied boundary element method to investigate the flow around 2D hydrofoils near a free surface [17,18].
The other important parameter that affects the performance of the shallowly submerged hydrofoil is its three-dimensional geometry; in other words, aspect ratio. It is worthy of note that, the term "shallowly submerged" is used to indicate the foil is operating near a free surface, not in very small submergence. In very close proximity to free surface, different flow phenomenas like wave breaking and cavitation may occur. These conditions change the physics, thus the numerical modelling of the problem. For this reason, it should be kept in mind that, shallow submergence term is used to refer the free surface effect, and the paper deals with the moderate submergence values. The assumption of two-dimensionality has some advantages in terms of solution time and procedure and provides useful information to the general insight into the problem of submerged objects near a free surface. But in real practice, the foil is always three dimensional. As mentioned before, the proximity of the foil to the free surface creates a wave pattern on the fluid surface. Literature survey so far shows that this interaction changes the flow field around the foil and the hydrodynamic forces that act upon it. With the absence of free surface, it is known that three-dimensionality effects and aspect ratio has a strong influence on the performance of a foil. When the foil is moving close a free surface, three dimensionality of the geometry might also alter the interaction between the foil and the free surface. Besides, the aspect ratio of the foil will also affect the resulting wave structure and make the problem more complicated. Thus, investigation of the threedimensionality and aspect ratio effects has importance for reliable determination of the performance of a hydrofoil moving under a free surface.
Daskovsky [4] experimentally and theoretically investigated the problem of the hydrofoil in surface proximity. He used the horseshoe vortex model and biplane model, which are both based on potential flow theory. Predictions of both models showed that decrease of the aspect ratio of the foil gives rise to the lift reduction factor, for submergence depths up to h/c=5 (here, h is the free surface distance of the foil and c is the chord length). Bal et al. [19] applied panel method to study cavitating hydrofoils under a free surface. They analysed the flow around 2D and 3D hydrofoil geometries and compare the results with other methods in the literature. Zhu et al. [1] studied horizontal and vertical 3D oscillating hydrofoils near the free surface by a hybrid method. They found that, in horizontal foil case, kelvin waves at the free surface, which are resultant of forwarding motion of the foil are dominated by the unsteady waves generated by the oscillatory motion. Effect of oscillation for submerged hydrofoils near a free surface is also studied by Esmaeilifar et al. [20], too. In his potential flow based numerical study, Bal [3] investigated the effect of free surface for high speed submerged 3D hydrofoils. He also included the tandem case into the study. He reported that the free surface could have either increasing or decreasing effect on lift and drag values of the hydrofoils with respect to the longitudinal distance between the tandem hydrofoils. Xie and Vassalos [2] analysed the performance of a 3D NACA4412 hydrofoil under the free surface using the potential-based panel method. They have tested different Froude Numbers and submergence depths in the range of 1<h/c<5. Their results indicate that the lift force acting on the foil decreases with the decrease of the aspect ratio. Ali and Karim [21] carried out a RANS based numerical study to investigate the effect of free surface on the shallowly submerged NACA0012 hydrofoil. They have mounted their foil geometry from wall to wall in the solution domain, similar to the 2D experiment of Duncan. They have presented the wave elevations at the centerline of the foil for different submergence ratios. Ghassemi and Kohansal [22] numerically investigated the wave generated by a NACA 4412 hydrofoil in various forms and presented the free surface deformations in various hydrofoil forms. Different geometry or flow configurations such as hydrofoil with winglets [23], flapping-foil [5] and slotted hydrofoil [24] are investigated by several researchers for hydrofoils operating near the free surface. De silva and Yamaguchi [25] utilized a viscous solver to examine the effect of design parameters on the performance of a oscillating hydrofoil under a free surface. Filippas [26] performed the numerical solution of the 3D, unsteady and nonlinear problem of Yavuz Hakan Ozdemir, Taner Cosgun, Wave Field Generated by Finite-Span Hydrofoils Baris Barlas Operating Beneath a Free Surface flapping foil under a free surface and in waves using an efficient GPU accelerated boundary element methodology. Filippas et al. [27] compared BEM and RANS based numerical solvers to test the flapping foils beneath a free surface. Ship interactions of actively controlled flapping foils were investigated by Belibassakis and Filippas [28]. The effect of variable bathymetry and sheared currents presented in Filippas et al. [29]. Even some previous research presents the investigation of the performance parameters and/or some flow features like cavitating behaviour in hydrofoils under free surface, the structure of the wave field and free surface interactions of the shallowly submerged hydrofoils are not fully understood.
The paper deals with the numerical solution of the flow around 2-D and 3-D NACA0012 hydrofoils moving beneath a free surface. Potential flow theory is applied to analyze the flow field. The 2D hydrofoil is simulated first, to validate the present numerical methodology. The numerical results obtained by the present study are compared with those of Duncan [6]. 3D calculations are carried out for three different submergence depths (where 0.75≤h/c≤1.5) and aspect ratios, and a wide range of Froude Numbers for each of these parameters. 3D hydrofoils are analysed with a free surface width of larger than foil span, to allow the free surface to deform in third direction, and thus, make the problem truly treedimensional. Numerical results are reported in terms of lift coefficients, pressure distributions around the foil surface, wave elevations on the midsection of the foils and three-dimensional free surface deformations.

Geometry and Conditions
The flow around the fully submerged two and three-dimensional NACA0012 hydrofoils working in free surface proximity is numerically analyzed using iterative boundary element method. Schematical description of the problem can be seen in Figure.1. Hydrofoil has a chord length of c and span of s for three-dimensional cases. Submergence depth is measured as the vertical distance between the leading edge of the foil and the calm water-free surface level. The incoming flow passes the foil with a uniform speed of U. The solution domain is consisting of the flow field, which is bounded by the hydrofoil surface and the free surface. It is assumed that flow is incompressible, inviscid and irrotational, except the wake region of the foil. A Cartesian coordinate system with a positive -z axes pointing the upwards direction is located at the foil surface. Positive -x-axes have the same direction as the flow velocity. framework of the potential flow. Kinematic and dynamic boundary conditions are applied on the free surface to model the interactions with the foil in close vicinity and to ensure the pressure on the free surface is equal to the atmospheric pressure. Velocity on the surface is tangent to the wall owing to the kinematic boundary condition on the surface. Kutta condition is used at the trailing edge of the foil to satisfy the speed at where the flow leaves the trailing edge is finite. Mathematical details of the boundary conditions are described in detail in the following sections.

Solution Strategy
The numerical background of this study is based on the boundary element method. 2D and 3D NACA0012 hydrofoils in an unbounded fluid are analyzed first to determine the accuracy of the code. After that, hydrofoil beneath the free surface is solved with the similar conditions of Duncan's [6] experimental arrangement, both in two and three dimensions. Tank width in Duncan's experiment is nearly equal to the foil span. Therefore, free surface measurements have no signs of three-dimensionality effects and experimental results are assumed twodimensional. These results are also widely used in the literature for 2D validations [7,10,13].
In the present study, 2D and 3D simulations of Duncan's case was performed to compare the wave structure. Numerical results are validated with Duncan's experimental data, and similar periodic wave structure is obtained with 3D simulations. The periodic and 2D-like wave structure is an inevitable result when the foil extends from wall-to-wall in the solution domain.
Later on, free surface width is expanded in the spanwise direction, to allow three-dimensional free surface deformations. Different aspect ratios are tested for various submersion depths and Froude numbers using this configuration. The purpose here is to investigate the effect of different aspect ratios of hydrofoil while allowing the water surface to deform freely in three dimensions. Two-dimensional versions of final configuration are also simulated for comparisons. Details of studied cases are summarized in Table 1.

Formulation
It is accepted that the fluid is incompressible, inviscid, and the flow is irrotational. A Cartesian coordinate system is placed on the foil surface. The components of the free stream velocity U in the x-z frame of reference. The angle of attack α is defined as the angle between the free stream velocity and the x-axis.
The total velocity potential can be splitted into a sum of a free stream potential and a perturbation potential [30,31]; Where is the perturbation potential. The total velocity potential  should satisfy the Laplace's equation in the fluid domain "  ".
( ) The domain  is bounded by the hydrofoil surface  The problem can be constructed by specifying the boundary conditions as follows: I-Kinematic body boundary condition: Where is the unit normal vector, which points out of the fluid domain.  [22,31]: Therefore, the gradient of perturbation potential should go to zero. Trailing edge of the foil should also be treated with Kutta condition, implies that flow leaves the sharp trailing edge smoothly;

III.
Kutta condition at the trailing edge of a foil: To ensure the wake surface is force-free, following assumption of the wake surface should be made, which states that the pressure on the two sides of the wake is equal [30]: Therefore, streamwise velocity must be continuous across the surface [30,31]: subscript U and L are the upper and lower surface of the wake, and V and P are the velocity and pressure respectively. t is a unit vector in the direction of the free stream velocity. The wake surface has zero thickness and the pressure jump across SW is zero, while there is a jump in the potential: where the constant  is the circulation around the body. We consider the hydrofoil at a constant velocity U at the surface of a fluid of infinite depth. Additional four boundary conditions on the free surface are given as follows:

IV. Kinematic free surface condition:
Kinematic free surface boundary condition can be written: V. Dynamic free surface condition: From Bernoulli's equation, the following equation can be given as, 22 1 where g is the gravitational acceleration. Substituting Eq. (10) in Eq. (9) and omitting the second order terms, following linearized free surface equation can be found as, VI. Radiation boundary condition It is necessary to impose a condition to ensure that the free surface waves vanish upstream of the disturbance. The upstream radiation condition gives [32]:

Numerical implementation
The velocity potential  can be written an integral equation on the foil surface by applying Green's theorem: where The hydrofoil is discretized by 60 line elements in which the mesh is clustered near the leading and trailing edge of the foil to employ the 2-D boundary element method. In 3-D simulations, the body is discretized into 12 strips in spanwise and 30 strips in the chordwise direction.
These elements estimate the actual geometry by a straight line while the unknown quantity supposed to be constant. While ij G is evaluated numerically, using a standard Gaussian quadrature ij H is calculated analytically. A more detailed discussion of this point is provided in Katsikadelis [33].

Determination of the influence coefficient the perturbation potentials
The above set of algebraic equations is solved by using the Gaussian elimination method. Iterative boundary element method is used to solve the unknown perturbation potentials on the hydrofoil. According to Green's third identity the perturbation potential on the hydrofoil surface and on the free surface can be expressed as, where SFS boundary of the free surface. The iterative boundary element method comprised of two parts: (i) the hydrofoil part, which solves for the unknown perturbation potentials on the hydrofoil surface, and (ii) the free surface part, which solves for the unknown perturbation potentials on the free surface. The potential in the flow domain due to the influence of the hydrofoil can be written as follows: On the other hand, the potential in the flow domain due to the influence of the free surface can be written as follows: By substituting equation (20)  A more detailed discussion of this point is provided in Bal et al. [19]. In this study, 22   x term is calculated by using fourth-order backward finite difference method. Free surface panels was created same dimensions and the first and second derivatives of  with respect to x is given as, where Dx is the free surface panel length. 2 2   x can be calculated by using equation (26). In order to avoid upstream waves first and second derivatives of  with respect to x are forced to be equal to zero [19].

NUMERICAL VALIDATION
Numerical validation and accuracy of the present iterative boundary element methodology were determined step-by-step. Computed lift coefficients of 2D and 3D NACA0012 hydrofoils in an unbounded fluid were compared with Xfoil results, as the first step of the numerical validation. Then, wave elevation of a point vortex (source) travelling under a free surface was validated with analytical results, to see the reliability of the free surface calculations. After that, the wave profile of the two-dimensional hydrofoil moving close the free surface was numerically obtained and compared with the experimental data of Duncan [6]. Additionally, the three-dimensional flow field, which is similar to the experimental arrangement of the Duncan, was analyzed to investigate the free surface deformation.
Comparison of the lift coefficients from present BEM computations and numerical results from the XFoil in an unbounded fluid is given in Figure 3. Present BEM perfectly captures the linear increase of the L C with the angle of attack. The dependency of lift on the aspect ratio is visible. In an unbounded fluid, it is well known that lift coefficient and attack of angle is linear before the stall.  The 2D lift coefficient calculated by the iterative method is given with the iteration number in Figure 5. Through computations, maximum of five iterations was found to be sufficient to obtain reliable results. Maximum number of iteration was determined by continuously observing the difference between the iteratively computed values of the wave elevation at each point and requiring that the difference be less than 10 The wave profiles of present 2D calculations compared with that of Duncan [6] is shown in Figure 6. The agreement of the present numerical results and experimental ones are quite good. Although some distinctions are seen at the first wave crest, wavelength and amplitudes are in good agreement in general. The difference between the experimental data and the numerical results is due to the neglect of viscosity and the use of the linearized free water surface boundary condition. Furthermore, Karim MM et al. [12], Wu and Chen [14], Prasad et al. [13] modelled turbulent free surface flow in Duncan [6]'s using viscous solvers. As can be seen from these studies, viscous solutions also modelled the wave deformations with some discrepancies, as the potential solutions With the same flow conditions of 2D validation, iterative boundary element method with the linear free surface condition is applied to the submerged rectangular 3D NACA 0012 hydrofoil. The views of panels used on the free surface and rectangular hydrofoil are shown in Figure 7. The number of panels on the free surface is 960 (80 panels in the x-direction (NX=80) and 12 panels (NY=12) in the y-direction). In this study, a planar quadrilateral panel and constant-strength singularity similar with reference [19] was used which is the simplest and most basic three-dimensional element. As mentioned before, the hydrofoil is discretized into 12 strips in spanwise and 30 strips in the chordwise direction. Upstream distance has initially been taken as one λ at the front of the hydrofoil, three λ lengths behind the hydrofoil and S length along the beam respectively.  The wave deformation obtained by iterative method is in Figure 8. It can be seen that free surface deformation is very similar to 2D results. Hydrofoil generates a periodic wave field on the free surface. Wave deformations have no signs of three dimensionalities, as mentioned by Duncan [6]. 2D like wave structure is an expected result of the wall-to-wall extension of the hydrofoil.

RESULTS
The detailed investigation of the flow field around the 3D NACA 0012 hydrofoil moving beneath the free surface is presented in this section. All results are generated with a 5° angle of attack. Differently from the previous analysis, here, free surface width is expanded to 2.0 times of foil span, to observe the wave field created by a finite-span hydrofoil.  Figure 9. Numerical predictions indicate that performance of the hydrofoil significantly affected by the free surface proximity.  The lift coefficients seem to be considerably affected by the aspect ratio, especially for low and moderate Fn's. We can see that the lift coefficient for 3D NACA hydrofoil growths with its aspect ratio. Similar trend was observed by Xie and Vassalos [2] for NACA4412 hydrofoil in surface proximity. As a well-known behaviour in foil theory, a lift coefficient of a 3D hydrofoil approaches to the 2D value, as the aspect ratio increases.
One may look for the grounds of the different hydrodynamic lift with increasing Fn's at the pressure distribution around the hydrofoil. The velocity potential along the tangential direction is differentiated to compute the local external tangential velocity on each collocation point, as follows [36]: Then, pressure coefficient can be obtained by using the Bernoulli Equation: Wave Field Generated by Finite-Span Hydrofoils The evaluation of surface pressure coefficients on 2D and 3D hydrofoils at different Froude numbers for fixed submergence ratio h/c=0.75 is depicted in Figure 10. The reason to choose this single submersion ratio is having the strongest variations of L C with Fn. So the closest position to the free surface, h/c=0.75, will be used for investigations at the rest of the study. As shown in Figure 10, the pressure on the lower surface does not show the considerable difference and always higher than the upper surface of the hydrofoil. However, the pressure on the upper side (the free surface side) is strongly influenced by the Froude number and this effect advanced gradually from the leading edge to trailing edge of the hydrofoil. As the Froude numbers getting higher, the difference between the suction and pressure sides of the hydrofoil decreases. This trend is matching up with the L C values at Figure 9 explains the increase or decrease in lift values.  Figure 9.
The pressure distribution on the free surface side of 3D hydrofoil is shown in Figure 11. Flow leakage can be clearly seen from the figure. Free surface proximity creates lower values of P C at Fn=0.5 and the region of the low pressure outspreads towards the leading edge. Pressure distribution at Fn=0.9 has a more uniform view, also higher in magnitude. Operating Beneath a Free Surface Figure 12 presents the wave elevations in the midsection of the 3D hydrofoil, which is very similar to the 2D computations. A small difference can be seen in the wave amplitudes at Fn=0.7. Increment of Fn from 0.5 to 0.7 gives rise to the wave amplitudes. On the other hand, a further increase in operating velocity lowers the amplitude of the generated wave field. A similar tendency was observed for NACA 0006 hydrofoil in the numerical work of Bal and Kinnas [7]. Amplitudes in their wave contours rise with the increment of Fn from 0.5 to 1.0., and then reduces again at Fn=1.5.
Free surface deformations and wave shades for 3D hydrofoil are given in Figure 13 and Figure 14, respectively. Note that, Figure 14 presents the perspective view of the free surface, and generated waves are symmetric with respect to -xz plane, as shown in Figure 13. As can be seen in both figures, hydrofoil creates a periodic wave pattern at Fn=0.5. Wavefield in that speed is very similar to the case where the foil is located from wall-to-wall, at Figure 8. With the rise of Fn to 0.7, periodic wave formation turn into diverging Kelvin wave form, and evolve into far-field of the wake. Note that, Fn=0.7 is also in the velocity range where lift reduction begins at Figure 9. For Fn=0.9, Kelvin-type wave pattern becomes more apparent.
Wave structures become very different from 2D cases for Fn=0.7 and 0.9, even they are in very good agreement in the midsection of the hydrofoil (see Figure 12). Generally speaking, flow needs to exert its kinetic energy to generate waves. Thus, one can expect a rise in the amplitude of the waves created by a hydrofoil at high Fn, due to its greater kinetic energy. However, free surface results indicate that generated waves are also beginning to spill in the spanwise direction in the form of Kelvin waves, as the Fn increase. This mechanism might be linked with the change of wave amplitudes in Figure 12 and variable difference of lift values of 2D and 3D hydrofoils at Figure 9, especially in the range of 0.5<Fn<0.9.

CONCLUSIONS
The paper covers the numerical simulation of flow around shallowly submerged 2D and 3D NACA 0012 hydrofoils. Calculations are carried out using an iterative boundary element method. The effect of the operating velocity and aspect ratio on the foil performance and wavefield are investigated at three different free surface proximities. In the preliminary 3D analysis, numerical towing tank width is taken as equal to the foil span, for the purpose of validation with the available experimental data from the literature. Then, free surface width is extended to see the effect of finite-span hydrofoil (in other words aspect ratio) on the wave field. The outcomes of the study can be summarized by the following: • Fn has a significant effect on the lift of the hydrofoil. Different values of Fn acts upon the lift force of the hydrofoil alternatively, within the same h/c ratio. Lower values of Fn, lift gives rise to the hydrodynamic lift, but after a specific point, this effect is shifted and increasing Fn reduces the lift. The effect of Fn on the lift force diminishes with the increasing submersion depth. • Lift coefficient for 3D hydrofoil growths with its aspect ratio, as expected. Two dimensional and high aspect ratio results are close to each other for higher values of Fn. However, in the range of for low and moderate Fn, strong variations are observed with the potential flow based lift results of 2D and 3D hydrofoils, even in AR=8. Their variations become apparent with increasing free surface proximity. • Pressure distribution on the upper side of the 2D hydrofoil is strongly influenced by the free surface, as the lower side Cp values do not show the considerable difference. Change of the pressure coefficient with Fn is in accordance with the variation of lift As the Fn increases, the Kelvin wave pattern is formed and extends in the broader region in the far-field. • Wave spill due to the Kelvin pattern might be related to the change of wave amplitudes and lift forces with Fn.