CHEN Qing-ren,YE Heng-kui,GUAN Yan-min

(College of Naval Architecture and Ocean Engineering,Huazhong University of Science and Technology,Wuhan 430074,China)

1 Introduction

Hydrofoils have being widely used in high-speed boats.The hydrofoil performance is more related to the submergence depth,moving speed and aspect ratio.The effect of free surface is unnegligible as the submergence depth becomes small.

In early time researches,Giesing and Smith[1],Yeung and Bouger[2],Bal and Han[3]carried out the theoretical investigations based on potential theory for two-dimensional hydrofoil moving beneath free surface.Park and Chun[4]used a higher order boundary element method to investigate the flow around a three dimensional hydrofoil.Wang and Zou[5]applied a panel method to calculate steady three-dimensional lifting potential flows about hydrofoils moving beneath the free surface.Zhang and Wang[6]utilized a B-spline method to predict the performance for three-dimensional hydrofoils with emphasis on the tip region.Xie and Vassalos[7]used a potential-based panel method to deal with the steady moving three-dimensional hydrofoil by Rankine sources distribution and double-body model approach.

In this paper,a time domain panel method as in the Refs.[8-9]is developed for the hydrodynamic analysis of a three-dimensional hydrofoils under free surface.The Rankine sources and dipoles are distributed on the body surface and the free surface.The dipoles are also shed on the wake surface in each time step.Kutta condition is satisfied.In the numerical computation,an analytical method[10-11]is introduced to calculate all the singular integrals and nonsingular integrals of influence coefficient matrix for improving the computation accuracy and speed of the numerical integrations.The time-stepping[12-13]iterative computation is performed to convergence.The numerical results of pressure distribution,lifting force,wave resistance and wave patterns are calculated for different Froude number,submergence depth and aspect ratio to examine the effect of free surface on performance of the hydrofoil.

2 Theoretical formulation

Consider the flow of a hydrofoil moving at a constant velocity in an inviscid and incompressible fluid.The o-xyz Cartesian coordinate system is chosen,as shown in Fig.1.The oxy plane is on the undisturbed free surface and oz-axis is positive upward and through the centre of leading edge of the foil.If the field point P （ x,y, ）z is situated on the boundary surface S,Green’s theorem may be written as[14]

Fig.1 Coordinate system

where φPis the velocity potential on the boundary surface,r is the distance between the field point and the source point,and S consists of the hydrofoil body surface SB,the free surface SFand the wake surface SW.is the unit normal vector of the boundary surface directing outwards to the fluid domain.δPis the solid angle facing to the flow domain on the position of point P.δPis defined as[15]

The Eq.(1)may be rewritten as

where Δφ is the potential jump across the wake surface.

The kinematic boundary condition on foil surface is

The wave elevation on the free surface SFmay be defined by a function η=η （ x,y, ）t.Being to Beinoulli’s equation,the dynamic boundary condition on the free surface is that,

where p0represents the pressure on the free surface.Neglecting the second-order quantity in fluid velocity,Eq.(5)is rewritten as

A Kutta condition which requires the velocity at the trailing edge to be finite should be satisfied:

The initial conditions are given as:

3 Numerical implementation

To develop a numerical method for the integral equation Eq.(3)which determines the potential φ on the foil surface and the normal velocityon the free surface,the boundary surface is divided into N quadrilateral panels and the control point is on the center of the panel.

where N is the total number of the panels,NB,NFand NWrepresent the number of the panel on the foil body surface,the free surface and the wake surface separately.

Eq.(3)can be discrectized as

where influence coefficient matrix:

where φ+and φ-are the instantaneous velocity potential at the upper and lower trailing edge of the foil,respectively.The doublets are shed down from trailing edge of the foil with the strength Δφjat each time step to advance the solution.The doublets on the wake surface move downstream at the distance lwwhich is shown in Fig.1.

Generally,a Gaussian integration is used to calculate the influence coefficients in expression(11).In this paper,an analytical method is employed to perform all the singular integral and non-singular integral of influence coefficient computation to improve accuracy of integrals for r approaching to zero.These types of surface integrals often occur on the body surface near the trailing edge of the foil.As shown in Fig.2,the quadrilateral element P1P2P3P4is the integration element,point O is the field point and O′is the projected point of O on the plane P1P2P3P4.The directed area Sp1p2p3p4can be defined as

Fig.2 The subdivision of the integration element and the polar coordinate

where SO′P3P4is negative and the others are positive.Then the integral over the quadrilateral is equal to the sum of the integrals over four triangles respectively.

By introducing the polar coordinate system （ρ, ）θ with the origin at point O′and the polar axis O′X′perpendicular to line P1P2,the integral of the integrandover triangle O′P1P2becomes

where h is the perpendicular distance from field point O to the integration plane and d is the perpendicular distance from point O′to line P1P2.The integrals over the other triangles can be carried out in the same way.

Applying the body boundary condition(4)and initial conditions(8),the unknown potential on the foil surface and normal velocity on the free surface at initial time can be obtained by solving Eq.(10).Once the wave elevation has been determined by the normal velocity on the free surface,the potential on the new free surface at a subsequent time t+Δt can be obtained using iterative expression(6).Also the doublets on the wake surface are updated by expressions(12)and(13).The procedure is repeated to get a solution at the next time step until the free surface wave profile trends to steady.

Using the solution of Eq.(10),the potential at the nodes of elements on the foil surface can be computed as follows:

The velocity distribution on foil surface can be obtained by employing the interpolation of corner-point potentialwith the finite element method.The pressure distribution on foil body can be calculated from the Bernoulli’s equation.

The wave resistance RWand lifting force L of the hydrofoil can be obtained by integrating pressure over the body surface SB.The wave resistance and lifting force coefficient are defined as:

where SCis the center plane area of the hydrofoil.

4 Numerical results

In the numerical computation example,NACA4412 section is selected for the hydrofoil.In order to test this numerical method,firstly,the flow around a high aspect ratio foil is calculated,and the results at the middle longitudinal section could be regarded as a 2D result,which are compared with analytic results.Secondly,the numerical results for 3D foil are analyzed.Finally,the effects on lifting force and wave resistance caused by submergence depth,aspect radio and moving velocity of the foil are discussed.

In the numerical computation,using symmetry,the number of panels are chosen as NB=800,NF=2 200,NW=1 200,the panel spacing of wake surface in x-direction lw=V·Δt,time step Δt=0.01s,the total time t=6s.

4.1 Comparison with the results of a 2D hydrofoil

The time domain panel method computation is employed to analyze the hydrodynamic force and wave elevation for a rectangular wing with aspect ratio AR=16,chord length Froude number Fn=1,submergence depth h/c=1 and angle of attack α=5°.The calculated pressure distribution on the body surface of the middle section of hydrofoil and the comparisons with that in 2-D computation given by Yeung and Bouger[2]are shown in Fig.3.Fig.4 shows the nondimensional wave elevation on the free surface for Fn=0.9 at y=0(the longitudinal central plane of hydrofoil).The calculated results are in good agreement with those in 2D computation.

Fig.5 and Fig.6 show the comparisons of the lifting force and wave resistance coefficients at middle section of the hydrofoil for different Froude number.The aspect ratio in present com-putation is AR=16 and AR=6 in Xie and Vassalos[7].It can be seen that the present results agree well with 2D data and this agreement is much closer than the values calculated by Xie and Vassalos because of the higher aspect ratio.

Fig.3 Pressure distribution at middle section for NACA4412 hydrofoil with AR=16

Fig.4 Wave elevation at middle section for NACA4412 hydrofoil with AR=16

Fig.5 Comparison of lift forcing coefficient at middle section of the hydrofoil

Fig.6 Comparison of wave resistance coefficient at middle section of the hydrofoil

4.2 The analysis of the character of the numerical results for 3D foil

Figs.7 and 8 display the contour of potential and pressure distribution on the half body surface of the hydrofoil respectively.It can be seen that the potential and pressure distribution on the hydrofoil vary significantly on the suction surface.The relation between the coefficient CW,CLand the computation time is shown in Fig.9.The wave resistance and lifting force change fast at the beginning of the computation,then attain slowly the steady state with computation time,in which this stable value may be regarded as the steady result in time domain.

Fig.7 The potential distribution on the hydrofoil

Fig.8 Pressure distribution on the hydrofoil

Fig.9 CWand CL-computation time curve in time domain

Fig.10 Instantaneous wave elevations at y=0 for the hydrofoil with AR=6

In Fig.10,the instantaneous wave elevations in the longitudinal central section of the free surface are given to show the wave generation process,which is caused by a hydrofoil with h/c=1 and Fn=1.5.The wave crest and trough grow gradually as the hydrofoil moves,until they achieve the steady state.Fig.11 expresses the divergent wave pattern of the half free surface.

Fig.11 Divergent wave contour for the hydrofoil

4.3 Effects on lifting force and wave resistance caused by various parameters

In Fig.12,lifting force coefficients of foil are shown as a function of Froude number for two aspect ratios of AR=4,6,respectively.As the Froude number increases,the lifting forces decrease and trend to an invariable value.It can be seen that the margin of lifting force between two hydrofoils with different aspect ratios decreases at a higher Fn,which shows that the three dimensional effect of hydrofoil becomes stronger as the moving velocity of hydrofoil increases.

Figs.13 and 14 show the effect of submergence depth on lifting force and wave resistance coefficient at various Froude numbers.It can be seen that the free surface effect almost does not exist at a large submergence depth,while as submergence depth decreases,the effect becomes stronger and stronger.At a range of submergence depth,the free surface has a positive effect on the lifting force at fixed Froude numbers of Fn=0.6,0.8,1.0,which is similar to the ground effect,but the lifting force decreases at all submergence depths for a higher Froude number,i.e.,Fn=1.5.It is found that the negative effect of the free surface on wave resistance becomes stronger for smaller submergence ratios.

Fig.12 Effects of Froude numbers on lift coefficient for two aspect radios

Fig.13 Effects of submergence depths on lift coefficients at various Froude numbers

Fig.14 Effects of submergence depths on wave resistance coefficients at various Froude numbers

5 Conclusions

A time domain panel method based on Green’s theorem has been developed in the present paper.The computed results show that this numerical method can be applicable for a three-dimensional hydrofoil moving under the free surface.An analytical method used to compute the integrals of influence coefficients can improve the computation accuracy.By comparing with other results given in the literature,good agreements are obtained to show the availability of this method.

For a specified submergence depth,the lifting force coefficient with the Froude number increasing,decreases and trends to an invariable value.The three dimensional effect of hydrofoil becomes stronger at a higher speed.It is found that the effects of the free surface on lifting force and wave resistance become stronger for smaller submergence depth ratios.However,the free surface may increase the lifting force at some submergence depth ratios.

The time-stepping scheme enables this method having the flexibility to be extended to resolve the nonlinear problems.The 3D numerical wave tank based on this method has bright application prospect in study of 3D lifting body with free surface,and it can be extended to many other nonlinear fields of hydrodynamics research,such as flapping hydrofoil,combined hull-foil case and so on.

[1]Giesing J P,Smith A M O.Potential flow about two-dimensional hydrofoils[J].Journal of Fluid Mechanics,1967,28:113-129.

[2]Yeung R W,Bouger Y C.Hybrid integral-equation method for the steady ship-waves[J].International Journal for Numerical Methods in Engineering,1979,14:317-336.

[3]Bai K J,Han J H.A localized finite element method for the nonlinear steady waves due to a two-dimensional hydrofoil[J].Journal of Ship Research,1994,38:42-51.

[4]Park I K,Chun H H,et al.Free surface flow analysis around three dimensional hydrofoils by a higher order boundary element method[C]//Proceedings of China-Korea Marine Hydrodynamics Meeting.Shanghai,China,1997.

[5]Wang Chun,Zou Zaojian.A 3D panel method for hydrodynamic forces on hydrofoils[J].Journal of Ship Mechanics,2001,5(3):18-25.

[6]Zhang Lijun,Wang Yanying.The investigation of hydrodynamic performance for three-dimensional hydrofoil based on B-spline grid[J].Journal of Hydrodynamics,Ser.A,2006,21(3):381-387.(in Chinese)

[7]Xie N,Vassalos D.Performance analysis of 3D hydrofoil under free surface[J].Ocean Engineering,2007,34:1257-1264.

[8]Yang Xianghui,Ye Hengkui,et al.Time domain research on the wave making of moving submerged body[J].Journal of Ship Mechanics,2007,11(6):845-858.

[9]Kara F,Tang C Q,Vassalos D.Time domain three-dimensional fully nonlinear computations of steady body-wave interaction problem[J].Ocean Engineering,2007,34(5):776-789.

[10]Yuan Zhengqiang,Huang Jian,Zhu Jialin.Exact integration of the linear element of Laplace equation in boundary element method[J].Chinese Journal of Applied Mechanics,2004,21(3):117-120.

[11]Sladek V,Sladek J,Tanaka M.Evaluation of 1/r integrals in BEM formulations for 3-D problems using coordinate multitransformations[J].Engineering Analysis with Boundary Elements,1997,20:229-244.

[12]Singh S P,Sen D.A comparative linear and nonlinear ship motion study using 3-D time domain methods[J].Ocean Engineering,2007,34(13):1863-1881.

[13]Li Guoan,Ye Hengkui.Wave making computation in time domain for multi-hull ships[J].Journal of Hydrodynamics,Ser.B,2006,18(4):225-230.

[14]Newman J N.Marine hydrodynamics[M].Cambridge:The MIT Press,1977.

[15]Yang S A.An investigation into integral equation methods involving nearly singular kernels for acoustic scattering[J].Journal of Sound and Vibration,2000,234(2):225-239.