Analytical modelling for a three-dimensional hydrofoil with winglets operating beneath a free surface

Applied Mathematical Modelling 37 (2013) 2679–2701

Contents lists available at SciVerse ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Analytical modelling for a three-dimensional hydrofoil with winglets operating beneath a free surface H. Liang a, L. Sun a, Z. Zong a,⇑, L. Zhou a, L. Zou b,c a School of Naval Architecture Engineering, State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, PR China b School of Aeronautics and Astronautics, State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, PR China c Department of Mathematics, Imperial College London, London SW7 2AZ, UK

a r t i c l e

i n f o

Article history: Received 29 September 2011 Received in revised form 25 May 2012 Accepted 6 June 2012 Available online 15 June 2012 Keywords: Lifting surface theory Green’s function Hydrofoil Winglet Free surface

a b s t r a c t The current study focuses on establishing a theoretical lifting surface model for predicting the hydrodynamic loads acting on the three-dimensional hydrofoil with winglets, which is considerably influenced by the proximity to the free surface through finding the threedimensional Green’s function for the planar and vertical horseshoe vortices operating below a free surface. The hydrofoil surface is decomposed into a finite number of elements along the span direction and the chord directions, each of which can then be represented by a horseshoe vortex. The linearized free surface boundary condition is applied to analyze the influence of the free surface on the hydrofoil as well as the winglets. The thickness problem is considered using the source distribution among the hydrofoil and winglets surfaces and the analytical Green’s function that satisfies the linearized free surface boundary condition is used. As a sample application, numerical examples were conducted to show the performance of the hydrodynamic characteristics for the hydrofoil with winglets as a function of the Froude number. It was concluded that there are significant efficiency benefits from using winglets inside the free surface proximity effect. These results are substantiated by the comparison with the available published data. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction Determining the hydrodynamic loads acting on the hydrofoil with winglets in the fluid in the presence of a free surface can be quite challenging. However, this topic is of practical importance. Recent needs for the marine vehicles capable of transporting heavier loads and operating at higher speed have prompted the development of hydrofoil vessels [1] featuring high-speed and good sea-keeping performance. There have been studies on hydrofoils in the presence of a free surface [2–17] and winglets [18–27]. There has been a review on hydrofoils and hydrofoil vessels [2], which states the development process of hydrofoils as well as hydrofoil vessels since 1950s. In recent years, research on the hydrodynamic characteristics of the hydrofoil operating in the vicinity of a free surface has captured the attention of fluid dynamicists, applied mathematicians, computational scientists and engineers [28]. Despite of that, the experimental investigations on the hydrodynamic performance of the hydrofoil still draw attention to researchers especially with the emergency of the particle image vecilometry (PIV) technique [3–8], for the data obtained ⇑ Corresponding author. E-mail addresses: [email protected] (H. Liang), [email protected] (Z. Zong). 0307-904X/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.apm.2012.06.020

2680

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

from the experiments can be available as a benchmark for the validation of the numerical and theoretical analysis. Among them, Dular et al. [3] conducted an experimental investigation on the cavitating hydrofoil with the assistant of PIV and laser induced fluorescence (LIF) techniques to evaluate the capability of a commercial CFD code (Fluent). Daskovsky [4] carried out an experiment to study the free surface proximity effect on the lift performance of the hydrofoil. von Ellenrieder and Pothos [5] measured the wake field behind a two-dimensional heaving hydrofoil for different Strouhal numbers using PIV. The deflected wake was observed. Ducoin et al. [6] conducted an experimental and numerical investigation of the time–space distribution of the wall-pressure field on a hydrofoil undergoing a transient up-and-down pitching motion at four pitching velocities and a fixed Reynolds number. Foeth et al. [7] investigated the fully developed sheet cavitation on a hydrofoil with a spanwise varying angle of attack by the aid of PIV. Its recording revealed the growth of the attached cavity and the cloud shedding. Sarraf et al. [8] experimentally studied the hydrodynamic behavior of a two-dimensional hydrofoil, and the study was focused on the hysteretic behavior at a stall angle. In recent years, different numerical methods have been developed to treat hydrodynamic parameters generated by the hydrofoil. Especially the boundary element method (BEM) has been widely used to analyze the non-cavitating as well as cavitating hydrofoils. For instance, Xie and Vassalos [9] employed BEM to calculate the hydrodynamic behaviors of a three-dimensional hydrofoil, including lift coefficient, wave-making resistance and free surface elevation. And their work was extended to investigate the wave pattern generated by surface piercing strut and Wigley hull by Ghassemi et al. [10], and the results were examined with comparison of experimental data. For cavitating hydrofoils, Bal and Kinnas [11] investigated the performance of the two- and three-dimensional cavitating hydrofoils operating under a free surface steadily with varying Froude numbers using BEM. The influence of the number of iterations on the results was discussed. Later, Bal generalized his work to the cases of the surface piercing hydrofoil [12] and the hydrofoil in water of finite depth [13]. Moreover, Fridman [14] obtained the analytical solution to the two-dimensional nonlinear problem of the flow past a planning flat with stagnation zone in the spoiler vicinity. Saito et al. [15] utilized the homogeneous model to study the three-dimensional unsteady cavitating flow around a hydrofoil. Kelly and Xiong [16] presented a model for a free deforming hydrofoil featured with self propulsion. The qualitative agreement was found with the observed dynamics of fishlike locomotion. Wang and Ostoja-Starzewski [17] applied the large eddy simulation (LES) scheme to simulate the sheet/cloud cavitation. However, the former researches focused on the hydrofoil without winglets. It has been known that winglets can enhance the lifting characteristics of a wing while diminishing the induced drag effects [21]. There is therefore an increased requirement to investigate the hydrofoil with winglets operating beneath a free surface for the engineering purpose. In addition, the winglets are usually used on the hydrofoil [1]. The research into the winglet was pioneered by Whitcomb [18] in the mid 1970s. It was found that properly designed winglets can go a long way towards reducing the induced drag [19]. Theoretical predictions for the wing both with and without winglets were primarily carried out using the vortex lattice method developed by Ishimitsu [20]. Gall and Smith [21] generalized this approach to the biplane configuration with winglets, and their work was validated against the experimental data. Gerontakos and Lee [22] investigated the effect of winglet dihedral on the near field vortex flow experimentally, and the results were compared with a baseline wing. Gatto et al. [23,24] carried out experimental investigations to study the influence of bistable winglets and articulated winglets on enhancing the lift characteristics of a wing. Jansen et al. [25] and Weierman and Jacob [26] conducted the optimization for the nonplanar wing and the wing on the unmanned aerial vehicle (UAV) to improve the lift over drag ratio. Ogurek and Ashworth [27] conducted the experiment on a wing with winglets in and out of ground effect. Three different types of winglets were tested. Lifting surface theory is originally proposed by Glauert [29], and developed by Weissinger [30]. The key techniques of the lifting surface theory lie in dividing the thin hydrofoil surface into a finite number of elements, and each of them can be represented by a horseshoe vortex [31]. Similarly, the winglets attached to the tips of the hydrofoil can also be divided into elements. The bound vortex is expediently placed on the quarter-point line of the element, and the free vortices consist of two semi-finite vortex filaments stretching after forwards to infinity. Expediently, the collocation point locates at the three-quarter chord station measured from the leading edge. The kinematic flow condition that zero normal flow passes across the wing’s solid surface should be satisfied at the collocation point [32]. The work described in the present article investigates the combined effects of the winglets and the free surface proximity on the hydrodynamic properties of a three-dimensional hydrofoil though finding the three-dimensional free surface Green’s function for the horseshoe vortex element. The differentia between the present work and the previous study is that the Green’s functions satisfying the linearized free surface condition, rather than the Rankine source method [9,10], are derived and evaluated to deal with the lifting problem of the hydrofoil with winglets. In addition, we solve the lifting problem from the vortex point of view, rather than the dipole. Though the singularity may lead to the difficulties of the numerical calculation, this does not prove to be a problem. A great deal of work [33–38] has been performed and several alternative approximations have been given to calculate the singular integral effectively. Compared with the Rankine source method, we can just distribute the singularities and horseshoe vortex on the hydrofoil and winglets surfaces rather than the whole flow domain using the present approach, which can greatly reduce computing times. For the outline of this paper, in Section 2 we introduce the mathematical model for the hydrofoil with winglets operating beneath a free surface. In Section 3, we seek for the three-dimensional free surface Green’s function for the isolate horseshoe vortex on the hydrofoil as well as the winglet. In Sections 4, we deduce the formulae of the lifting force, wave-making resistance, induced drag and the free surface deformation. In addition, numerical results presented show the influences of three

2681

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

(a)

(b)

(c)

Fig. 1. Three types of winglet configuration: (a) upper winglet, (b) lower winglet and (c) full winglet.

types of winglets (shown in Fig. 1) and the free surface on the hydrodynamic properties of the hydrofoil compared with a baseline hydrofoil. In Section 5, a short conclusion is presented. 2. Governing equations and boundary conditions Throughout the paper the fluid is assumed to be incompressible and inviscid. The free surface disturbance is assumed sufficiently small, and thus the nonlinear effect can be neglected. Non-cavitating hydrofoil is accepted herein. It is supposed that the height of the upper winglet is less than the submergence, which means that the winglet cannot pierce the free surface. A right hand coordinate system o_xyz is employed, and let the o_xy plane be at the undisturbed free surface. In addition, oz axis is positive upward as shown in Fig. 2. Taking the hydrofoil as the frame of reference, the water can be regarded to flow past the hydrofoil at a constant velocity of U steadily. The velocity potential U(x, y, z) in the flow field satisfies following conditions [9]:

r2 U ¼ 0; jrUj < 1 at trailing edge of the hydrofoil; Uðx; y; zÞ ! Ux far away upstream;

ð1Þ ð2Þ ð3Þ

where Eq. (1) is the Laplace equation and Eq. (2) satisfies the Kutta condition. We suppose that a linear body boundary condition and the linear free surface condition are satisfied on the mean free surface (z = 0) [39]:

o2 U oU oU þk l ¼ 0; ox2 oz ox

ð4Þ

where l = l0 /U, and l0 is the Rayleigh ‘‘artificial’’ friction coefficient. In addition, k = g/U2, in which g is gravity acceleration, denotes wave number. In Eq. (4), U can be decomposed into three parts: the velocity potential of the isolate horseshoe vortex on the hydrofoil surface or the winglet u, the velocity potential of the image horseshoe vortex ui and the free surface velocity

y z

o Element

h

x Hydrofoil surface

Wingle Fig. 2. Schematic of a three-dimensional hydrofoil with winglets under a free surface.

2682

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

potential G. Here we replace the velocity potential by its induced velocity components, and the induced velocities of these components are computed from the Biot–Savart law [40]. Based on the potential flow theory, Eq. (4) can be rewritten in the form

o2 G oG oG o2 o o o þk l ¼  2 ðu þ ui Þ  k ðu þ ui Þ þ l ðu þ ui Þ ¼  ðu þ ui Þ  kðw þ wi Þ þ lðu þ ui Þ; 2 ox oz ox oz ox ox ox

ð5Þ

where u and w are the perturbation velocities induced by u in directions of ox and oz, and ui and wi denote the perturbation velocities induced by ui in directions of ox and oz. 3. Three-dimensional free surface Green’s function The perturbation velocity induced by a horseshoe vortex on the hydrofoil or the winglet can be calculated from the law of Biot–Savart. By substituting the perturbation velocities into Eq. (5), we can obtain the three-dimensional free surface Green’s function. For the purpose of simplifying the process of dealing with the boundary condition, we resort to the Fourier-type integration. The Fourier-type integral expression in the form of the inverse of the distance can be expressed as [41]:

1 1 ¼ R 2p

Z

Z p

1

emjzfjþim½ðxnÞ cos hþðygÞ sin h dh dm:

ð6Þ

p

0

3.1. Green’s function for the horseshoe vortex on the hydrofoil surface The flow past the horseshoe vortex on the hydrofoil surface is shown in Fig. 3. Based on Eq. (6), the perturbation velocity induced by the horseshoe vortex on the panel of the hydrofoil surface in the form of Fourier-type integration can be obtained (see Appendix A in details):

uH ¼

C 8p 2

Z

y2

y1

Z

1

Z p

memðzþhÞþim½x cos hþðygÞ sin h dh dm dg

ð7Þ

p

0

and

wH ¼

C

Z

8p2

y2

Z

y1

1

Z p

im cos hemðzþhÞþim½x cos hþðygÞ sin h dh dm dg þ

p

0

 sin hemðzþhÞþim½ðxnÞ cos hþðyy2 Þ sin h dh dm dn 

Z

C 8p

2

0

1

Z

1

0

C

Z

1

Z

1

Z p

im 8p2 0 p 0 Z p im sin hemðzþhÞþim½ðxnÞ cos hþðyy1 Þ sin h dh dm dn;

ð8Þ

p

where y1 and y2 are two end points coordinate in the direction of oy for the bound line vortex, and they satisfies y1 < y2. Then, C denotes the circulation of the horseshoe vortex on that element. The superscript H denotes the hydrofoil. Accordingly, the perturbation velocity induced by the image horseshoe vortex of the panel above the free surface in the form of Fourier-type integration can be obtained:

uiH ¼

C 8p2

Z

y2

y1

Z

1

Z p

memðzhÞþim½x cos hþðygÞ sin h dh dm dg;

ð9Þ

p

0

Horseshoe vortex on the hydrofoil surface

(x, y, z) B(x0, y2, –h)

R

dl

z A(x0, y1, –h) y o

x

Fig. 3. Schematic of the horseshoe vortex on the hydrofoil surface under a free surface.

2683

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

Horseshoe vortex on the winglet element

(x, y, z) R

D(x0, y0, z2) z

dl y

o

C(x0, y0, z1) x

Fig. 4. Schematic of the horseshoe vortex on the winglet moving beneath a free surface.

and

wiH ¼ 

Z

C 8p2

Z

y2

1

Z p p

0

y1

im cos hemðzhÞþim½x cos hþðygÞ sin h dh dm dg 

 sin hemðzhÞþim½ðxnÞ cos hþðyy2 Þ sin h dh dm dn þ

Z

C 8p

2

1

0

Z 0

1

Z p

Z

C 8p2

Z

1

0

Z p

1

im

p

0

im sin hemðzhÞþim½ðxnÞ cos hþðyy1 Þ sin h dh dm dn:

ð10Þ

p

The expression of GH can be obtained by substituting Eqs. (7)–(10) into boundary condition Eq. (5):

GH ¼

Z

C 4p2

y2

Z

y1

1

Z p

im cos h þ l emðzhÞþim½x cos hþðygÞ sin h dh dm dg: m cos2 h  k þ il cos h

p

0

ð11Þ

3.2. Green’s function for the horseshoe vortex on the winglets The flow past the horseshoe vortex on the winglet is shown in Fig. 4. Similarly, the perturbation velocity induced by the horseshoe vortex on the panel of the winglet in the form of Fourier-type integration can also be calculated:

uW ¼

Z

C 8p 2

z2

Z

z1

1

Z p

im sin hemðzfÞþ

im½ðxx0 Þ cos hþðyy0 Þ sin h

dh dm df

ð12Þ

p

0

and

C 8p2

wW ¼

Z

1

Z

x0

1

Z p

im sin hemðzz2 Þþim½ðxnÞ cos hþðyy0 Þ sin h dh dm dn þ

p

0

C

Z

8p 2

1

Z

x0

1

0

Z p

im

p

 sin hemðzz1 Þþim½ðxnÞ cos hþðyy0 Þ sin h dh dm dn;

ð13Þ

where z1 and z2 are two end points coordinate in the direction of oz for the bound line vortex, and they satisfy z1 < z2. The superscript W denotes the winglet. Simultaneously, the perturbation velocity induced by the image horseshoe vortex of the panel on the winglet above the free surface can be obtained:

C

uiW ¼

Z

8p2

z2

Z

1

im sin hemðzþfÞþim½ðxx0 Þ cos hþðyy0 Þ sin h dh dm df;

ð14Þ

p

0

z1

Z p

and

wiW ¼

C

Z

8p2

1

Z

x0

1

0

Z p

im sin hemðzþz2 Þþim½ðxnÞ cos hþðyy0 Þ sin h dh dm dn 

p

C 8p2

Z

1

x0

Z

1 0

 sin hemðzþz1 Þþim½ðxnÞ cos hþðyy0 Þ sin h dh dm dn:

Z p

im

p

ð15Þ

Thus, the expression of GW can be obtained by substituting Eqs. (7)–(15) into Eq. (5):

GW ¼

Z

C 4p

2

z2 z1

Z 0

1

Z p p

i sin hðim cos h þ lÞ emðzþfÞþim½ðxx0 Þ cos hþðyy0 Þ sin h dh dm df: m cos2 h  k þ il cos h

ð16Þ

2684

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

(a)

L1

L2

(b)

Fig. 5. Two integral paths around the singularity.

3.3. Simplification of the Green’s function There exist singularities when m = k sec2h  il sec h where value of h covers the interval [p,p]. While jhj < p/2, sec h is greater than 0, and the singularity locates under the real axis. Thus, the path L1 should be passed as shown in Fig. 5. With l ? 0, the singularity approaches the real axis from below. Similarly, when jhj > p/2, the path L2 should be passed to round the singularity. Thus, Eq. (9) can be rewritten according to the integral path:

GH ¼

Z

C 4p2

y2

Z p=2 Z p=2

y1

1

0ðL1 Þ

im cos h C emðzhÞþim½ðxx0 Þ cos hþðygÞ sin h dm dh dg þ 2 m cos2 h  k 4p

Z

y2

Z p=2 Z p=2

y1

1

0ðL2 Þ

im cos h  emðzhÞim½ðxx0 Þ cos hþðygÞ sin h dm dh dg: m cos2 h  k

ð17Þ

The integral paths L1 and L2 can be decomposed into the principle value integrals and their residues

GH ¼

Z

y2

Z p=2

Z

1

im sec h iC emðzhÞþim½ðxx0 Þ cos hþðygÞ sin h dm dh dg  m  k sec2 h 4p y1 p=2 4p 0 Z 1 Z y2 Z p=2 C 2 2 PV  sec3 hek sec hðzhÞþik sec h½ðxx0 Þ cos hþðygÞ sin h dh dg  2 4p y1 p=2 0 Z y2 Z p=2 im sec h i C  ik emðzhÞim½ðxx0 Þ cos hþðygÞ sin h dm dh dg  m  k sec2 h 4p y1 p=2

C

2

 sec3 hek sec

2

PV

hðzhÞik sec2 h½ðxx0 Þ cos hþðygÞ sin h

Z

y2

Z p=2

ik

p=2

y1

dh dg:

ð18Þ

Our purpose is to evaluate the real component of GH, and thus Eq. (18) can be simplified as

GH ¼

Z

C 2p2

y2

Z p=2

y1

PV

Z

p=2

 sec3 hek sec

2

0

1

im sec h C emðzhÞþim½ðxx0 Þ cos hþðygÞ sin h dm dh dg þ m  k sec2 h 2p

hðzhÞþik sec2 h½ðxx0 Þ cos hþðygÞ sin h

Z

y2

y1

Z p=2

k

p=2

dh dg:

ð19Þ

It has been known that the principle value integral in the form of

I ¼ PV

Z

1

0

1 emðzhÞþim½ðxx0 Þ cos hþðygÞ sin h dm m  k sec2 h

ð20Þ

can be rewritten as [42] k

I ¼ ek E1 ðkÞ  sgnðxÞpie ;

ð21Þ

where

k ¼ k sec2 hðz  hÞ þ ik sec2 h½ðx  x0 Þ cos h þ ðy  gÞ sin h and

x ¼ ðx  x0 Þ cos h þ ðy  gÞ sin h: The expression ekE1(k) can be evaluated using the Pade approximation [43] which can accurately and efficiently convert the improper integral to a polynomial reducing the complexity and computing time [44] (see Appendix B). Thus, Eq. (21) can be expressed using the complex exponential integral E1

GH ¼

Z p=2

C 2p

2

csc h sec h½ek1 E1 ðk1 Þ  ek2 E1 ðk2 Þ dh 

p=2

C 2p

Z p=2 p=2

where

ki ¼ k sec2 hðz  hÞ þ ik sec2 h½ðx  x0 Þ cos h þ ðy  yi Þ sin h and

xi ¼ ðx  x0 Þ cos h þ ðy  yi Þ sin h:

  i csc h sec h ½1 þ Hðx1 Þek1  ½1 þ Hðx2 Þ ek2 dh;

ð22Þ

2685

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

Similarly, Eq. (12) can be expressed in the form in accordance with the integral path

G

W

¼

Z

C 4p

2

z2

z1

Z p=2 Z p=2

1

0ðL1 Þ

m tan h  emðzþfÞ m  k sec2 h

m tan h emðzþfÞþ m  k sec2 h im½ðxx0 Þ cos hþðyy0 Þ sin h

im½ðxx0 Þ cos hþðyy0 Þ sin h

dm dhdf 

Z

C 4p

2

z2

z1

Z p=2 Z

1

p=2

0ðL2 Þ

dm dhdf

ð23Þ

and then, Eq. (21) is equivalent to the combination of the principle integral and its residue

Z

z2

Z p=2

Z

1

m tan h iC emðzþfÞþim½ðxx0 Þ cos hþðyy0 Þ sin h dm dh df þ 2h m  k sec 4p2 z1 p=2 4 p 0 Z 1 Z z2 Z p=2 C 2 2  tan h sec2 hek sec hðzþfÞþ ik sec h½ðxx0 Þ cos hþðyy0 Þ sin h dh df  2 PV 4p z1 p=2 0 Z z2 Z p=2 m tan h iC  k emðzþfÞim½ðxx0 Þ cos hþðyy0 Þ sin h dm dh df  m  k sec2 h 4p z1 p=2

GW ¼ 

C

PV

 tan h sec2 hek sec

2

hðzþfÞ ik sec2 h½ðxx0 Þ cos hþðyy0 Þ sin h

Z

z2

z1

Z p=2

k

p=2

dh df:

ð24Þ

Based on Eqs. (20) and (21), Eq. (24) can be rewritten as

GW ¼

C 2p 2

Z p=2

tan h½ ek1 E1 ðk1 Þ  ek2 E1 ðk2 Þ dh 

p=2

iC 2p

Z p=2

½1 þ Hðx0 Þ tan hðek1  ek2 Þ dh;

ð25Þ

p=2

where

ki ¼ k sec2 hðz þ zi Þ þ ik sec2 h½ðx  x0 Þ cos h þ ðy  y0 Þ sin h and

x0 ¼ ðx  x0 Þ cos h þ ðy  y0 Þ sin h: 3.4. Gradient of the Green’s function For the purpose of evaluating the induced velocity components, the gradient of the Green’s function should be considered. Taking the derivative of Eq. (22), we can obtain the velocity components induced by the horseshoe vortex on the hydrofoil surface

8 9 > < i cos h > = 1 1 3 k1 k2 v k csc h sec h½e E ðk Þ   e E ðk Þ þ   dh ¼ i sin h H 1 1 1 2 2 > > > k1 k2 > : ; : ; 2p p=2 wH 1 8 9 > Z p=2 < i cos h > =   C ik csc h sec3 h ½1 þ Hðx1 Þ ek1  ½1 þ Hðx2 Þek2  i sin h dh;  > > 2p p=2 : ; 1 8 9 > > < uH =

C

Z p=2

ð26Þ

where the relationship [45]

d Z 1 e E1 ðZÞ ¼ eZ E1 ðZÞ  dZ Z

ð27Þ

has been used. In the same way, the velocity components induced by the horseshoe vortex on the winglet can also be obtained

8 9 > i cos h >   < = 1 1 2 k1 k2 dh  v k sec h tan h e E ðk Þ   e E ðk Þ þ ¼ i sin h W 1 1 1 2 2 > > > k1 k2 > : ; : ; 2p p=2 wW 1 8 9 i cos h > > Z p=2 < = iC  ½1 þ Hðx0 Þk sec2 h tan hðek1  ek2 Þ  i sin h dh: > > 2p p=2 : ; 1 8 9 > < uW > =

C

Z p=2

ð28Þ

4. Results and discussion In this section, we investigate the hydrodynamic characteristics of the rectangular hydrofoils with and without winglets. The hydrofoil section we select is NACA-4412, and the winglet section we use is NACA-0012.

2686

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

4.1. The circulation distribution Generally, the hydrofoil surface and two winglets attached to the tips of the hydrofoil are symmetrical about the central plane. Suppose that the half hydrofoil surface and one winglet can be divided into N and Q elements in the spanwise direction, and M elements in the chord direction. Therefore there are 2N  M and 2Q  M elements placing on the hydrofoil and two winglets surface, respectively. The velocity of downwash on the control point includes two parts: the velocity induced by horseshoe vortices on all elements of the hydrofoil and winglets surfaces, and the velocity induced by the free surface disturbance velocity potential. Thus the velocity of downwash at the control point on the m, nth element of the hydrofoil surface can be expressed as:

2 wm;n s

13

0

2

13

0

6 6 C7 C7 B B C7 C7 M X 2N 6 M 2Q 6 B B X 6 m;n  m;n BoG H C7 XX6 m;n  m;n BoG W C7 ¼ 6wi;j þ wi;j þ ReB oz C7 þ 6wp;q þ wp;q þ ReB oz C7; 6 6 C7 C7 B B i¼1 j¼1 4 @z¼zm;n A5 p¼1 q¼1 4 @ z¼zm;n A5 y¼yn x¼xm

ð29Þ

y¼yn x¼xm

where wm;n represents the total velocity of downwash at the control point on the m, nth element of the hydrofoil surface. s  m;n wm;n i;j and wi;j denote the velocity of downwash at the control point on the m, nth panel induced by the horseshoe vortex on the i, jth element of the hydrofoil surface and its image horseshoe vortex.  m;n wm;n p;q and wp;q stand for the velocity of downwash at the control point on the m, nth panel induced by the horseshoe vortex on the p, qth element of the winglet and its image horseshoe vortex. In addition, xm, yn and zm,n represent the coordinate of the control point of the m, nth element. The derivative of GH and GW with respect to z can evaluated via Eqs. (26) and (28). Due to the presence of a free surface, the perturbation velocity in the direction of incoming flow also exists. Similarly, the velocity of incoming flow at each panel can be obtained:

2 13 13 0 6 6 C7 C7 B B 2Q 6  C7  C7 M X 2N 6 M X B B X 7 X 7 6 m;n 6 m;n   C B B oG oG m;n H W 7þ 6u 6u þ u  p;q þ ReB ox  C ¼Uþ C7 ox  C 7 7; 6  i;j þ ReB 6 p;q C C B B 6 i¼1 j¼1 6 @z¼z0m;n A7 @z¼z0m;n A7 5 p¼1 q¼1 4 5 4 2

um;n

0

y¼yn x¼x0m

ð30Þ

y¼yn x¼x0m

where um,n represents the velocity of the incoming flow at the centre point of the bound line vortex on the m, nth element of the hydrofoil surface.  m;n u i;j denotes the velocity at the centre point of the bound line vortex on the m, nth panel induced by the image horseshoe vortex on the i, jth element in the direction of ox.  m;n um;n p;q and up;q stand for the velocity at the centre point of the bound line vortex on the m, nth panel induced by the horseshoe vortex on the p, qth element of the winglet and its image horseshoe vortex in the direction of ox. In addition, x0m and z0m,n represent the coordinate of the centre point of the bound line vortex on the m, nth element. Based on the kinematic boundary condition requiring no flow across the hydrofoil surface, we can obtain a system of equations:

wm;n þ s

oGs o ¼ um;n f ðx; yÞ; ox oz

ð31Þ

where f(x, y) denotes camber function of the hydrofoil, and Gs, which is a known function, represents the velocity potential generated by the thickness effect (see Appendix C). However, Eq. (31) is an overdetermined system of linear equations with 2(N + Q)  M unknown numbers and 2N  M equations. Thus, the boundary conditions on the winglet should be required. The perturbation velocity in the direction of oy should be calculated due to the winglet plane vertical to oy axis which reads

2 13 13 0 6 6 C7 C7 B B 2Q 6 M X 2N 6  C7 M X  C7 B B X 7 X 7 6 g;h 6 g;h B B g;h g;h oG H  C7 oG W  C7 6 6 ¼ oy  C oy  C 7þ 7; 6v i;j þ v i;j þ ReB 6v p;q þ v p;q þ ReB C C B B 6 i¼1 j¼1 6 @ z¼zh A7 @ z¼zh A7 5 p¼1 q¼1 4 5 4 2

v g;h s

0

y¼r=2 x¼xg

ð32Þ

y¼r=2 x¼xg

v g;h denotes the total velocity at the control point on the g, hth element of the winglet in the direction of oy. s v and v g;h i;j denote the velocity in the direction of oy at the control point on the g, hth panel induced by the horseshoe

where

g;h i;j

vortex on the i, jth element of the hydrofoil surface and its image horseshoe vortex.

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

2687

 g;h v g;h p;q and v p;q stand for the velocity in the direction of oy at the control point on the g, hth panel induced by the horseshoe vortex on the p, qth element of the winglet and its image horseshoe vortex. In addition, xh, y = ±r/2 and zg represent the coordinate of the control point of the g, hth element, for the winglet locates at the tip of the hydrofoil, where r represents the span of the hydrofoil. The plane of the winglet is parallel to the incoming flow, therefore the boundary condition stating no normal flow across the solid surface can be expressed as:

Fig. 6. Circulation distribution for the hydrofoils with and without winglets along the half span at h/c = 1.0, Fr = 1.0, AR = 6.0 and a = 5.0°.

2688

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

Fig. 6 (continued)

v g;h s þ

oGs ¼ 0: oy

ð33Þ

Combining Eq. (31) with Eq. (33) yields a determined system of linear equation. We can obtain the circulation distribution of each element on the hydrofoil and the winglets. The circulation distribution across the span of the hydrofoil is investigated as shown in Fig. 6. In this case of calculation, we compare the circulation distribution for the hydrofoil with winglets with the case of no winglets at h/c = 1.0 and Fr = 1.0. The aspect ratio of the hydrofoil we set is 6.0, and the angle of attack is equal to 5.0 degrees. Fig. 6 illustrates the circulation distribution for the hydrofoils with full winglets and the case of no winglets. It can be clearly seen that the circulation is close to zero at the tips of hydrofoil without winglets. However, for the hydrofoil with winglets, the circulation at the tips is far away from zero due to the presence of winglets. Comparing with (a)-(j), we can find that the circulation is descending with farther distance from the leading edge. In addition, it is of interest here to note that the circulation loading is very close to elliptical for the hydrofoil without winglets. Nevertheless, the circulation distribution across the hydrofoil with winglets is in no way elliptical, and the circulation distribution for the presence of winglets is more uniform along the span than the case of no winglets, especially for the hydrofoil with larger winglets. Thus, the winglets force the flow to be more two-dimensional at the hydrofoils, which means larger lifting force. 4.2. Lift performance Given the circulation distribution, the lift coefficient for the hydrofoil with winglets can be calculated. The bound line vortex on the winglet is in the direction of oz, so the winglet cannot produce lifting force. Thus, the lift coefficient can be expressed as:

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

2689

M X 2N X

CL ¼

qU Ci;j Dri;j

L

qU 2 S=2

¼

i¼1 j¼1

qU 2 S=2

;

ð34Þ

where S denotes the area of the hydrofoil surface, and Dri,j represents the length of the bound line vortex on the i, jth panel. Consider a rectangular hydrofoil without winglet for different aspect ratios of 4.0, 5.0 and 6.0. The submergence-to-chord ratio h/c is 1.0. Based on Eq. (34) we can calculate the lift coefficient. Fig. 7 clearly shows the lift coefficient for a hydrofoil without winglets as a function of Froude number. The validity of the present work (solid line in Fig. 7) is examined against the published data which is obtained by Xie and Vassalos [9] (solid dots in Fig. 7). The angle of attack is 5 degrees, and submergence-to-chord ratio equals 1.0. By seeing Fig. 7(a)–(c), we observe that the lift coefficient for a three-dimensional hydrofoil without winglets increases with its aspect ratio as expected. The value of the lift coefficient varies dramatically as the Froude number is greater than 0.4 but less than 1.0, registered as increasing considerably before decreasing. In addition, there exists a peak at Froude number equaling 0.6 approximately in this interval. Thus, we can say that the free surface effect is significant while the Froude number is less than 1.0. Nice agreements have been found via comparison especially for a large Froude number (Fr P 0.8). Fig. 8 demonstrates the lift coefficient of the hydrofoil with AR = 6.0 and a = 5° varying with Froude number ranging from 1.0 to 5.0 for different Froude number. As can be easily seen from Fig. 8, the lift coefficient decrease with the immersion depth at Fr = 0.7, and is on the increase with the submergence at Fr = 1.5. As the Froude number is equal to 1.0, the lift coefficient is descending after ascending, and the peak value emerges at h/c = 2.0. Agreement between the present work and the results in [9] is reasonably good. As the submergence depth h goes to infinity, the hydrofoil degenerate to an airfoil which is not influenced by the free surface. Fig. 9 shows the lift performance of a biplane with winglets as a function of attack angle. The foil section is selected as NACA-0012. The attack angles of the upper and the lower foil are equal to 5°, and the height of the winglet is equal to the chord length of the foil. Good agreement can be observed from Fig. 9, and thus the present method can predict the flow characteristics of the hydrofoil with winglets. Fig. 10 prescribes the lift coefficient for the hydrofoil with and without winglets versus Froude number, and three types of winglet configurations are investigated. It is supposed that the height of the winglet is 0.5c. It can be clearly observed that the lift coefficient of the hydrofoil with winglets is greater than the case of no winglet at all times, and the hydrofoil with full winglet configuration is always larger than the other cases. However, the lift coefficient of the hydrofoil with upper winglets is larger than the case of the lower winglets while the Froude number is less than 1.2. As the Froude number is greater than 1.2, the reverse is happened (see Fig. 11). 4.3. Resistance performance As a non-lifting body operates beneath a free surface, the free surface is disturbed and energy is needed to push the water outward reflected as the wave-making resistance. Thus, wave-making resistance exists as long as the body is in proximity to a free surface. For a three-dimensional lifting body, the wake vortices induce velocity component downward, and thus the force is in the direction of the incoming flow due to the velocity of downwash which is registered as the induced drag. Even the threedimensional lifting body is in the free space, the induced drag still acts on the lifting body. Therefore, both of the induced drag and the wave-making resistance exert on the hydrofoil. 4.3.1. Wave-making resistance The wave-making resistance is caused by the pressure difference due to the presence of surface waves. Then, the wavemaking resistance can be expressed as [9]

W ¼

Z Z

qU

S

oUf cosðn; xÞ dS; ox

ð35Þ

where Uf denotes the velocity potential due to the free surface disturbance. Thus, the wave-making resistance coefficient can be obtained

CW ¼

W 2

qU S=2

¼

2 US

Z Z S

oUf fx ðx; yÞ dS: ox

ð36Þ

The wave-making resistance exists only if the hydrofoil locates close to the free surface. As the hydrofoil is far away from the free surface, the hydrofoil amounts to an airfoil whose wave-making resistance equals zero. From Eq. (36), we can obtain the wave-making resistance of a hydrofoil. Fig. 10 is the comparison of wave-making resistance at the central plane of the hydrofoil without winglets at h/c = 1.0, a = 5 degrees and AR = 6.0. The calculation results investigated by Xie and Vassalos [9] and two-dimensional calculations carried out by Yeung and Bouger [45] are used to examine the validity of the present approach. When the Froude number is larger than 0.4, the wave-making resistance is on the rise sharply, and it reaches the peak at

2690

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

Fig. 7. Lift coefficient versus Froude number for different aspect ratios compared with the results by Xie and Vassalos [9]: (a) AR = 4.0, (b) AR = 5.0, (c) AR = 6.0.

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

2691

Fig. 8. Lift coefficient varying with h/c for different Froude numbers compared with the results by Xie and Vassalos [9]: (a) Fr = 0.7, (b) Fr = 1.0, (c) Fr = 1.5.

2692

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

Fig. 9. Comparison of the lift coefficient of a biplane with winglets with the data obtained from Gall and Smith [46].

Fig. 10. Lift performance of the hydrofoil with and without winglets varying with Froude number at AR = 6.0, a = 5°.

Fig. 11. Comparison of wave-making resistance at the central plane of a hydrofoil without winglets at h/c = 1.0, a = 5 degrees and AR = 6.0.

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

2693

Fig. 12. Wave-making resistance acting on the hydrofoil with and without winglets at h/c = 1.0, a = 5 degrees and AR = 6.0.

Fr = 0.8. Then, the wave-making resistance decreases with Froude number as Froude number is greater than 0.8. The agreement is satisfactory. Fig. 12 presents the wave-making resistance for the hydrofoil with and without winglets, and three types of winglet configurations are investigated. Taking the hydrofoil without winglets as the benchmark, we can find that the wave-making resistance for cases of upper winglets and full winglets is greater than the case of no winglets. However, as the Froude number is larger than 1.0, the wave-making resistance for the hydrofoil with lower winglets is less than other cases observably.

4.3.2. Induced drag Apart from the wave-making resistance, the invisid induced drag for a three-dimensional hydrofoil should be taken into consideration on account of the existence of the velocity of downwash which is induced by the trailing vortex sheet and the velocity potential increment due to the presence of a free surface.

CD ¼

D 2

qU S=2

PM P2N ¼

i¼1

j¼1

qwi;js Ci;j Dri;j

qU 2 S=2

:

ð37Þ

Fig. 13(a) displays an example of the induced drag coefficient for a hydrofoil without winglets as a function of Froude number at h/c = 1.0, a = 5 degrees and AR = 6.0. As the Froude number is larger than 0.4 but less than 1.3, the induced drag coefficient goes down dramatically, and it reaches the trough at Fr = 1.3. Then, the induced drag coefficient is on the rise from quickly to slowly, and it will remain constant with increasing the Froude number. It should be noted that the induced drag coefficient is less than zero as the Froude number covers the interval of (0.9, 6.1) as shown in Fig. 13. At that time, the velocity of downwash is less than zero and becomes upwash under the influence of the free surface, registered as propulsion. Thus, we call this interval as the ‘‘propulsion zone’’. For a foil in free surface space, the velocity of downwash induced by the wake vortices sheet always points downward. However, the main difference between the hydrofoil and the foil is the free surface, and the free surface effect can also has an influence on the velocity of downwash. In a certain Froude number interval, the velocity of downwash induced by the free surface disturbance points upward, and the magnitude is larger than the velocity induced by the wake vortices sheet. At this time, the force in the direction of ox is less than zero which is the thrust rather than drag. The reason leading to the propulsion region is the fluctuation of the free surface. Comparing with Fig. 12, we can find that the induced drag coefficient and the wave-making resistance coefficient have the same order of the magnitude. However, the absolute value of the induced drag is much less than the wave-making resistance. Taking Fig. 13(a) as the benchmark, Fig. 13(b) exhibits the numerical examples of the induced drag for the hydrofoils with three types of winglets. For the hydrofoil with winglets, the induced drag is a little smaller (the thrust is larger) than the case without winglets, but it is not very obvious. The upper winglet has the better induced drag performance than the lower winglet, and the full winglet is the best. In addition, the presence of winglets can broaden the propulsion region. For the purpose of testing the validity of the code which is used to evaluate the induced drag for the hydrofoil with winglets, we investigate the relationship between the induced drag and the lift coefficient for a biplane with and without winglets as shown in Fig. 14. The satisfactory agreement with the results carried by Gall and Smith [46] indicates the validity of the present method.

2694

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

Fig. 13. Induced drag for the hydrofoil with and without winglets at h/c = 1.0, a = 5 degrees and AR = 6.0. (a) Induced drag for a hydrofoil without winglets; (b) Comparison of the induced drag between the case of winglets (hu = 0.5c; hl = 0.5c; hu = hl = 0.5c) and no winglets.

Fig. 14. Induced drag for the biplane with and without winglets versus lift coefficient at a = 5 degrees and AR = 5.0.

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

2695

4.4. Pressure coefficient The pressure distribution on the hydrofoil surface can be expressed as

CP ¼ 1 



2

rU U

¼

2u ; U

ð38Þ

where the second order terms can be neglected. Eq. (38) allows us to calculate the pressure coefficient on the hydrofoil. Fig. 15 presents the pressure coefficient at the central plan of the hydrofoil. The results performed by Xie and Vassalos [9] are given to examine the present approach, and Good agreement is achieved. Fig. 16 shows the pressure coefficients at the central section of the hydrofoil as a function of Froude number. It can be seen that the free surface has a more significant impact on the upper surface of the hydrofoil than the lower surface. Fig. 17 depicts the pressure coefficient at the tip section of the hydrofoil with and without winglets. A notable feature which can be easily observed is the pressure difference between the lower and upper surface for the hydrofoil without winglets is less than the case with winglets, and the largest difference emerges at the case of full winglets. This feature indicates that the lifting force at the tip section for the hydrofoil with winglets is much larger than the case of no winglets. By comparing with the cases of upper and lower winglets, we can observe that the upper winglets make the pressure coefficient of the hydrofoil’s upper surface higher, and conversely the pressure coefficient of the lower surface becomes larger due to the presence of the lower winglets.

Fig. 15. Comparison of the pressure coefficients at the central section for NACA4412 foil at a = 5°, AR = 10, h/c = 1.0 and Fr = 1.0.

Fig. 16. Pressure coefficients at the central section for a hydrofoil without winglets at a = 5°,AR = 10 for different Froude numbers.

2696

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

Fig. 17. Pressure coefficients at the tip section for the hydrofoil with and without winglets at a = 5°, AR = 6.

Fig. 18. Comparison of the free surface deformation at the central plane of the hydrofoil with Xie and Vassalos’ results [9] at Fr = 1.0, a = 5°, h/c = 1 and AR = 6.

Fig. 19. The free surface deformation at the central plane of the hydrofoil with and without winglets at Fr = 1.0, a = 5°, h/c = 1 and AR = 6.

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

2697

4.5. Free surface deformation Given the velocity potential in water, the free surface elevation can be obtained

f¼

U ou : g ox

ð39Þ

Fig. 18 is the comparison of the free surface elevation at the central plane of the hydrofoil, and the hydrofoil covers the interval of [0,1]. It can be clearly seen that there is no wave upstream. Compared with the results of Xie and Vassalos [9], the agreement is good. Fig. 19 presents the free surface deformations for the hydrofoil with full winglets and no winglets. We can find that the hydrofoil with winglets can induce a larger free surface deformation. 5. Conclusions In the present article, a theoretical lifting surface model is developed to analyze the hydrodynamic loads acting on the hydrofoil with three types of winglets. The difficulty but key task lies in seeking the three-dimensional free surface Green’s function for the horseshoe vortex on the hydrofoil surface and the winglet. The improper integration with infinite infamous singularities is tackled by converting it to a complex polynomial. The results of the analysis of the hydrofoil without winglets agree nicely with those of Xie and Vassalos [9] as well as Yeung and Bouger [45], and the results with winglets are examined with Gall and Smith [46]. These comparisons show the applicability of the present approach. From the results presented in this article, we can reach the following conclusions: (1) The winglets attached to the hydrofoil can produce increases in lifting force, and can also raise wave-making resistance. However, the winglets can reduce the induced drag. (2) The type of the winglet configuration has a significant influence on the hydrodynamic loads acting on the hydrofoil. In general, the full winglet configuration is better than other cases. (3) The circulation distribution across the hydrofoil without winglets is in proximity to elliptical. But for the hydrofoil with winglets, the circulation distribution is more uniform. The flow past the hydrofoil with winglets is more likely to the two-dimensional case, which can increase the lifting force and lessen the induced drag.

Acknowledgment The first author would like to thank Professor Francis Noblesse for discussing the evaluation of the Green’s function. The present work was supported by the National Natural Science Foundation of China (Grant No. 50921001 and 50909017), National Key Basic Research Special Foundation of China (Grant No. 2010CB832704) and Scientific Project for High-tech Ships: Key Technical Research on the Semi-planning Hybrid Fore-body Trimaran. Appendix A The horseshoe vortex lying on the hydrofoil is parallel to the o_xy plane, and the bound vortex segment is parallel to ox axis (as shown in Fig. 3). It is suppose that two end points of the bound line vortex are A and B, and their coordinate are (x0, y1, h) and (x0, y2, h). In accordance with Biot–Savart law, a three-dimensional vortex segment of strength C, lying on a closed curve C, can induce a velocity field:

V¼

Z

C 4p

R  dl R3

c

ðA:1Þ

:

For the bound vortex segment AB, the vectors R and dl are

R ¼ ðx  x0 ; y  g; z þ hÞ and

dl ¼ ð0; dg; 0Þ: Thus, the velocity components induced by the bound vortex segment AB can be expressed as

VAB ¼ 

C 4p

Z

y2

y1

ðz þ hÞi þ ðx  x0 Þk ½ðx  x0 Þ2 þ ðy  gÞ2 þ ðz þ hÞ2 3=2

dg:

ðA:2Þ

2698

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

For the wake vortex segment B1, the vectors R and dl are

R2 ¼ ðx  n; y  y2 ; z þ hÞ and

dl2 ¼ ðdn; 0; 0Þ: As a result, the velocity components induced by the wake line vortex B1 can be obtained

VAB ¼ 

Z

C 4p

y2

ðz þ hÞi þ ðx  x0 Þk ½ðx  x0 Þ2 þ ðy  gÞ2 þ ðz þ hÞ2 3=2

y1

dg:

ðA:3Þ

In the same way, the velocity components induced by the wake line vortex A1 can also be obtained

VA1 ¼

Z

C 4p

1

ðz þ hÞj  ðy  y1 Þk ½ðx  nÞ2 þ ðy  y1 Þ2 þ ðz þ hÞ2 3=2

x0

dn:

ðA:4Þ

By adding Eqs. (A.2)–(A.4) together, we can obtain the velocity components induced by the horseshoe vortex on the hydrofoil surface in directions of ox and oz

uH ¼

Z

C 4p

y2

ðz þ hÞ ½ðx  x0 Þ2 þ ðy  gÞ2 þ ðz þ hÞ2 3=2

y1

dg

ðA:5Þ

and

wH ¼

Z

C 4p

y2

ðx  x0 Þ

Z

C

dg 

4p ½ðx  x0 Þ þ ðy  gÞ þ ðz þ hÞ  Z 1 C ðy  y1 Þ þ dn: 4p x0 ½ðx  nÞ2 þ ðy  y1 Þ2 þ ðz þ hÞ2 3=2 2

y1

2

2 3=2

1

ðy  y2 Þ ½ðx  nÞ2 þ ðy  y2 Þ2 þ ðz þ hÞ2 3=2

x0

dn ðA:6Þ

Rewriting Eqs. (A.5) and (A.6) in the form of the reciprocal of the distance yields

uH ¼ 

Z

C 4p

y2

y1

o 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dg oz 2 ðx  x0 Þ þ ðy  gÞ2 þ ðz þ hÞ2

ðA:7Þ

and

wH ¼

Z

C 4p 

y2

y1

Z

C 4p

Z

o 1 C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dg þ ox 4p 2 2 2 ðx  x0 Þ þ ðy  gÞ þ ðz þ hÞ y2

y2 y1

o 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dn oy 2 ðx  nÞ þ ðy  y2 Þ2 þ ðz þ hÞ2

o 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dn: oy 2 ðx  nÞ þ ðy  y1 Þ2 þ ðz þ hÞ2

y1

ðA:8Þ

By substituting Eq. (6) into Eqs. (A.7) and (A.8), we can obtain the velocity components in the form of Fourier-type integral

uH ¼

C 8p 2

Z

y2

y1

Z

1

Z p

memðzþhÞþim½ðxx0 Þ cos hþðygÞ sin h dh dm dg

ðA:9Þ

p

0

and

wH ¼

C 8p2

Z

y2

y1

Z 0

1

Z p

im cos hemðzþhÞþim½ðxx0 Þ cos hþðygÞ sin h dh dm dg þ

p

 sin hemðzþhÞþim½ðxnÞ cos hþðyy2 Þ sin h dh dm dn   sin hemðzþhÞþim½ðxnÞ cos hþðyy1 Þ sin h dh dm dn:

Z

C 8p

2

y2

y1

Z

1 0

Z p

C 8p2

Z

y2

y1

Z 0

1

Z p

im

p

im

p

ðA:10Þ

It is supposed that the horseshoe vortex lying on the winglet is parallel to the o_yz plane, and the bound vortex segment is parallel to oz axis (as shown in Fig. 4). Two end points of the bound line vortex are C and D, and their coordinate are (x0,y0, z1) and (x0,y0,z2). For the bound vortex segment CD on the winglet element, the vectors R and dl are

R ¼ ðx  x0 ; y  y0 ; z  fÞ and

dl ¼ ð0; 0; dfÞ:

2699

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

Then, the velocity components induced by the bound vortex segment CD can be derived

V¼

Z

C 4p

z2

ðy  y0 Þi  ðx  x0 Þj ½ðx  x0 Þ2 þ ðy  y0 Þ2 þ ðz  fÞ2 3=2

z1

df:

ðA:11Þ

For the wake line vortex D1, the vectors R and dl are

R2 ¼ ðx  n; y  y0 ; z  z2 Þ and

dl2 ¼ ðdn; 0; 0Þ: Thus, the velocity components induced by the wake vortex segment D1 can be obtained

V2 ¼ 

Z

C 4p

1

ðz  z2 j  ðy  y0 ÞkÞ ½ðx  nÞ2 þ ðy  y0 Þ2 þ ðz  z2 Þ2 3=2

0

ðA:12Þ

:

From Eqs. (A.11) and (A.12), the velocity components in directions of ox and oz induced by the horseshoe vortex on the winglet can be expressed as:

uW ¼ 

Z

C 4p

z2

ðy  y0 Þ df

ðA:13Þ

½ðx  x0 Þ2 þ ðy  y0 Þ2 þ ðz  fÞ2 3=2

z1

and

wW ¼ 

Z

C 4p

1

ðy  y0 Þ 2

2

2 3=2

½ðx  nÞ þ ðy  y0 Þ þ ðz  z2 Þ 

x0

dn þ

C

Z

4p

1

x0

ðy  y0 Þ ½ðx  nÞ2 þ ðy  y0 Þ2 þ ðz  z1 Þ2 3=2

dn:

ðA:14Þ

Based on Eqs. (6), (A.13) and (A.14) can be expressed in the form of Fourier-type integral

uW ¼

Z

C 8p 2

Z

z2

Z p

1

im sin hemðzfÞþim½ðxx0 Þ cos hþðyy0 Þ sin h dh dm df

ðA:15Þ

p

0

z1

and

uW ¼

Z

C 8p2 þ

1

x0

C

Z

8p2

Z

1

Z p

im sin hemðzz2 Þþim½ðxnÞ cos hþðyy0 Þ sin h dh dm dn

0 p 1 Z 1 Z

x0

p

im sin hemðzz1 Þþim½ðxnÞ cos hþðyy0 Þ sin h dh dm dn:

ðA:16Þ

p

0

Appendix B The difficulty lies in seeking the solution to the complex function E1(k). There are several approaches to calculating the complex function E1(k). In the present article, we use the Pade approximation proposed by Smith et al. [43] which converts the expression to a polynomial. This approximation has higher precision and faster convergence compared with other approaches. The polynomial can be expressed as:

expðkÞE1 ðkÞ ¼

MþN þ e; D

ðB:1Þ

where

! 4 X i M ¼ ð1 þ m1 k þ m2 k þ m3 k þ m4 k Þ lnðkÞ ¼  mi k lnðkÞ; 2

3

4

ðB:2Þ

i¼0

! 5 X i ni k ; N ¼ cð1 þ n1 k þ n2 k þ n3 k þ n4 k þ n5 k Þ ¼ c 2

3

4

5

ðB:3Þ

i¼0

2

3

4

5

6

D ¼ ð1 þ d1 k þ d2 k þ d3 k þ d4 k þ d5 k þ d6 k Þ ¼

! 6 X i di k :

ðB:4Þ

i¼0

The parameters in Eqs. (B.2), (B.3) and (B.4) are shown in Table 1. In Eq. (B.3), c is known as Euler’s constant, c = 0.5772156649, and e is the error which satisfies jej 6 7  106. When the real part or image part of k is greater than 10 (Re (k)> 10 or Im (k) > 10), Eq. (B.1) can also be simplified as:

expðkÞE1 ðkÞ ¼

0:711093 0:278518 0:010389 þ þ þ e: k þ 0:415775 k þ 2:29428 k þ 6:2900

ðB:5Þ

2700

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

Table 1 Parameters for Pade approximation. m1

m2

m3

0.23721365

0.02065430

0.000763297

m4 0.000009768701

n1 1.4954589

n2 0.04180643

n3 0.03000059

n4 0.0019387339

n5 0.0005180156

d1 0.762736

d2 0.2838836

d3 0.066786

d4 0.012982

d5 0.0008701

d6 0.000299

Appendix C For a three-dimensional hydrofoil, the thickness problem can be addressed using the distribution of source. The strength of the source distributed over the hydrofoil surface can be expressed as [32]:

qðx; yÞ ¼ 2U

ogt ðx; yÞ; ox

ðC:1Þ

where gt denotes the thickness function. The influence of the free surface should be taken into consideration, for the hydrofoil operates near the free surface. The velocity potential for a point source with unit strength in proximity to the free surface is [47]

Gs ¼

Z 1 Z p=2 kðzhÞ 1 1 4m e cos½kðx  x0 Þ cos h cos½kðy  y0 Þ sin h   PV dh dk r r1 p k cos2 h  m 0 0 Z p=2 2  4m em sec hðzhÞ sin½m sec hðx  x0 Þ cos½mðy  y0 Þ sin h sec2 h sec2 h dh;

ðC:2Þ

0

where r and r1 are defined as

r¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx  x0 Þ2 þ ðy  y0 Þ2 þ ðz þ hÞ2

and

r1 ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx  x0 Þ2 þ ðy  y0 Þ2 þ ðz  hÞ2 :

A great deal of work on the numerical solution to Eq. (C.2) has been performed by dividing the flow field into the near-field, moderate-field and far-field according to the different distance [33–38]. Here, we use the analytical moving source Green’s function derived by Noblesse et al. [48] which is sufficiently accurate and can reduce the computing time significantly. The Green’s function for a source steadily advancing beneath the free surface in the analytical form is

4pGs ¼ Hðx  x0 Þ

4 F

Im 2

Z

t1

dtKeð1þt

2 ÞðzhÞ=F 2 þi

pffiffiffiffiffiffiffiffi

1þt 2 ½ðxx0 Þþtðyy0 Þ=F 2

t 1



1 þ GL ; r

ðC:3Þ

where 1/r is the fundamental free space singularity component, and GL is the local flow component which can be expressed analytically. Considering the derivative of the Green’s function in Eq. (C.3). The gradient of G can be expressed as

4prGs ¼ GW þ Gr  G0 r þ RL ;

ðC:4Þ G0r

where GW is the gradient of the wave component, and Gr and are the gradient of the fundamental free space singularity components, and RL is the gradient of the near-field component (see [48] in details). References [1] O.M. Faltinsen, Hydrodynamics of High-Speed Marine Vehicles, Cambridge University Press, Cambridge, 2005. [2] A.J. Acosta, Hydrofoils and hydrofoil craft, Annual Review of Fluid Mechanics 5 (1973) 161–184. [3] M. Dular, R. Bachert, B. Stoffel, B. Sirok, Experimental evaluation of numerical simulation of cavitating flow around hydrofoil, European Journal of Mechanics B/Fluids 24 (2005) 522–538. [4] M. Daskovsky, The hydrofoil in surface proximity, theory and experiment, Ocean Engineering 27 (2000) 1129–1159. [5] K.D. von Ellenrieder, S. Pothos, PIV measurements of asymmetric wake of a two dimensional heaving hydrofoil, Experiments in Fluids 44 (2008) 733– 745. [6] A. Ducoin, J.A. Astolfi, F. Deniset, J.F. Sigrist, Computational and experimental investigation of flow over a transient pitching hydrofoil, European Journal of Mechanics B/Fluids 28 (2009) 728–743.. [7] E.J. Foeth, C.W.H. van Doorne, T. van Terwisga, B. Wieneke, Time resolved PIV and flow visualization of 3D sheet cavitation, Experiments in Fluids 40 (2006) 503–513. [8] C. Sarraf, H. Djeridi, S. Prothin, J.Y. Billard, Thickness effect of NACA foils on hydrodynamic global parameters, boundary layer states and stall establishment, Journal of Fluids and Structure 26 (2010) 559–578. [9] N. Xie, D. Vassalos, Performance analysis of 3D hydrofoil under free surface, Ocean Engineering 34 (2007) 1257–1264.

H. Liang et al. / Applied Mathematical Modelling 37 (2013) 2679–2701

2701

[10] H. Ghassemi, M. Iranmanesh, A. Ardeshir, Simulation of free surface wave pattern due to the moving bodies, Iranian Journal of Science & Technology, Transaction B: Engineering 34 (2010). 177-134. [11] S. Bal, S.A. Kinnas, A BEM for the prediction of free surface effects on cavitating hydrofoils, Computational Mechanics 28 (2002) 260–274. [12] S. Bal, High-speed submerged and surface piercing cavitating hydrofoils, including tandem case, Ocean Engineering 34 (2007) 1935–1946. [13] S. Bal, The effect of finite depth on 2D and 3D cavitating hydrofoils, Journal of Marine Science and Technology 16 (2011) 129–142. [14] G. Fridman, Planing plate with stagnation zone in the spoiler vicinity, Journal of Engineering Mathematics 70 (2011) 225–237. [15] Y. Saito, R. Takami, I. Nakamori, T. Ikohagi, Numerical analysis of unsteady behavior of cloud cavitation around a NACA0015 foil, Computational Mechanics 40 (2007) 85–96. [16] S.D. Kelly, H. Xiong, Self-propulsion of a free hydrofoil with localized discrete vortex shedding: analytical modeling and simulation, Theoretical and Computational Fluid Dynamics 24 (2010) 45–50. [17] G. Wang, M. Ostoja-Starzewski, Large eddy simulation of a sheet/cloud cavitation on a NACA0015 hydrofoil, Applied Mathematical modeling 31 (2007) 417–447. [18] R.T. Whitcomb, A design approach and selected wind-tunnel results at high subsonic speeds for wing-tip mounted winglets, NASA TN D-8260 (1976). [19] H.H. Heyson, G.D. Riebe, C.L. Fulton, Theoretical parametric study of the relative advantages of winglets and wing-tip extensions, NASA TP-1020 (1977). [20] K.K. Ishimitsu, Aerodynamic design and analysis of winglets, AIAA Aircraft Systems and Technology Meeting, 27–29 September 1976, Dallas, Texas, AIAA-P-76-940. [21] P.D. Gall, H.C. Smith, Aerodynamic characteristics of biplanes with winglets, Journal of Aircraft, AIAA 24 (1987) 518–522. [22] P. Gerontakos, T. Lee, Effects of winglet dihedral on a tip vortex, Journal of Aircraft, AIAA 43 (2006) 117–124. [23] A. Gatto, F. Mattioni, M.I. Friswell, Experimental investigation of bistable winglets to enhance wing lift takeoff capability, Journal of Aircraft, AIAA 46 (2009) 647–655. [24] A. Gatto, P. Bourdin, M.I. Friswell, Experimental investigation into articulated winglet effects on flying wing surface pressure aerodynamics, Journal of Aircraft, AIAA 47 (2010) 1811–1815. [25] P.W. Jansen, R.E. Perez, J.R.R.A. Martins, Aerostructural optimization of nonplanar lifting surfaces, Journal of Aircraft, AIAA 47 (2010) 1490–1503. [26] J. Weierman, J.D. Jacob, Winglet design and Optimization for UAVs, 28th AIAA Applied Aerodynamics Conference, 28 June–1 July 2010, Chicago, Illinois, AIAA, 2010, pp. 2010–4224. [27] D.J. Ogurek, J. Ashworth, Experimental investigation of various winglet design for a wing in ground effect, in: 22nd Applied Aerodynamics Conference and Exhibit, 16–19 August 2004, Providence, Rhode Island, AIAA 2004-4720, 2004. [28] C. Pozrikidis, Fluid Dynamics: Theory, Computation, and Numerical Simulation, second ed., Springer, New York, 2009. pp. 680. [29] H. Glauert, The elements of airwing and airscrew theory, Cambridge University Press, Cambrige, 1947. [30] J. Weissinger, Uber eine Erweiterung der Prandtlschen Theorie der trahenden Linie, Math Nachr 2 (1949) 45–106. [31] H. Schlichting, E. Truckenbrodt, Aerodynamics of the Airplane, McGraw-Hill International Book Company, New York, 1979. pp. 121-125. [32] J. Katz, A. Plotkin, Low Speed Aerodynamics: From Wing Theory to Panel Methods, McGraw-Hill International Book Company, New York, 1991. pp. 386–390. [33] F. Noblesse, Alternative integral representations for the Green function of the theory of ship wave resistance, Journal of Engineering Mathematics 15 (1981) 241–265. [34] F. Noblesse, The Green function in the theory of radiation and diffraction of regular water waves by a body, Journal of Engineering Mathematics 16 (1982) 137–169. [35] J.N. Newman, Evaluation of the wave-resistance Green function: Part 1 – the double integral, Journal of Ship Research 31 (1987) 79–90. [36] F. Nobless, The fundamental solution in the theory of steady motion of a ship, Journal of Ship Research 21 (1977) 82–88. [37] F. Noblesse, On the fundamental function in the theory of steady motion of ships, Journal of Ship Research 22 (1978) 212–215. [38] F. Ursell, Mathematical note on the fundamental solution (Kelvin source) in ship hydrodynamics, Journal of Applied Mathematics, IMA 32 (1984) 335– 351. [39] I. Sahin, M.C. Hyman, Simulation of three-dimensional finite-depth wave phenomenon for moving pressure distributions, Ocean Engineering 28 (2001) 1621–1630. [40] J.N. Newman, Marine H ydrodynamics, Cambridge University Press, Cambridge, 1977. [41] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, Dover Publications, New York, 1972. [42] J.P. Giesing, A.M.O. Smith, Potential flow about two-dimensional hydrofoils, Journal of Fluid Mechanics 28 (1967) 113–129. [43] A.M.O. Smith, J.P. Giesing, J.L. Hess, Calculation of waves and wave resistance for bodies moving on or beneath the surface of the sea. Douglas Aircraft Co., Long Beach, California, Report No. 31488a, 1963. [44] D.C. Scullen, E.O. Tuck, Free-surface elevation due to moving pressure distributions in three dimensions, Journal of Engineering Mathematics 70 (2011) 29–42. [45] R.W. Yeung, Y.C. Bouger, Hybrid integral-equation method for the steady ship problem, International Journal for Numerical Methods in Engineering 14 (1979) 317–336. [46] P.D. Gall, H.C. Smith, Aerodynamic characteristics of biplanes with winglets, Journal of Aircraft, AIAA 24 (1987) 518–522. [47] J.V. Wehausen, E.V. Laitone, Surface waves, Surface waves online (1960). [48] F. Noblesse, G. Delhommeau, F. Huang, C. Yang, Practical mathematical representation of the flow due to a distribution of sources on a steadily advancing ship hull, Journal of Engineering Mathematics 71 (2011) 367–392.