含渦流作用之重力水波研究(2/2)

(1)

行政院國家科學委員會專題研究計畫成果報告

含渦流作用之重力水波研究(2/2)

計畫類別：個別型計畫計畫編號： NSC93-2611-E-002-012- 執行期間： 93 年 08 月 01 日至 94 年 07 月 31 日執行單位：國立臺灣大學水工試驗所計畫主持人：黃良雄報告類型：完整報告處理方式：本計畫可公開查詢

中華民國 94 年 10 月 31 日

(2)

Å

Åo

o

o»

»

I

I. .

.

õ

õº

º

ºÝ

Ý

_ÝÞ

_Þ

_Þ@

@

@~

~

i

Ï

_×

â

â:

:

:ø

ø

ø®

®

®à

à

¥

¥æ

¥

æ

i

i®

®

®@

@

@~

~

_(2/2)

The Study of Gravitational Water Wave

with Vortex Effect (2/2)

i_r

_{: NSC 93-2611-E-002-012}

Æ§

_{: 93}

O

₈

`

₁

^

₉₄

O

₇

`

₃₁

^

(3)

ê

1 Introduction 2

2 Mathematical Formulations 3

2.1 General Formulations of Two-Dimensional Linear Wave Problems . . 3

2.2 The Application of the Helmholtz Decomposition Theorem . . . 5

2.3 Integral Representations . . . 7

2.3.1 Irrotational Flow Field . . . 7

2.3.2 Rotational Flow Field . . . 8

3 The Vortex Methods 10 3.1 Vorticity Layer Approximation . . . 10

3.2 Vortex Blob Representation . . . 11

3.3 Vorticity Creation and Shedding . . . 12

4 Numerical Methods for Linear Wave Problems 13 5 Results and Discussions 14 5.1 Uniform flow over a rectangular obstacle . . . 14

5.2 Generation of surface waves by point vortex . . . 14

5.2.1 Analytical Solution . . . 14

5.2.2 Numerical results . . . 17

5.3 Linear wave over a rectangular obstacle . . . 17

(4)

Å

Åo

o

o»

»

I

I. .

.

õ

õº

º

ºÝ

Ý

_ÝÞ

_Þ

_Þ@

@

@~

~

i

Ï

_×

â

â:

:

:ø

ø

ø®

®

®à

à

¥

¥æ

¥

æ

i

i®

®

®@

@

@~

~

_(2/2)

The Study of Gravitational Water Wave

with Vortex Effect (2/2)

i_r

_{: NSC 93-2611-E-002-012}

Æ§

_{: 93}

O

₈

`

₁

^

₉₄

O

₇

`

₃₁

^

x

x¹

¹

¹ß

ß

_:

?

.

.>

>

>0

0

0 »

»

»ñ

ñ

È

.

.þ

þ

þa

a

.

Abstract

In this study, a grid-free numerical model for the simulation of gravita-tional water wave with vortex effect is proposed. The computagravita-tional algo-rithm is established by combining a potential flow field with a vortex flow field. The potential at the free surface and the free-surface displacement are solved by the time integrations. The vorticity field is the solution of the vor-ticity equation, which is solved using the vortex method and the vortex sheet method. The potential and velocity fields are solved by the boundary integral method.

Two cases are applied to validate the proposed model, and the vortic-ity shedding due to the rectangular obstacle in the linear wave is simulated. These results show the ability of the proposed model for the simulations of free-surface flows with vortex effect.

Keywords: gravitational water wave, vortex, boundary integral equation method, vortex method, numerical simulation

`

Í_{i2àTø:ø)ÕÝ]Pñ×myÕ½/xç}} Õ:øÅ(¥æi®óÂÿP Í@~2à Helmholtz DecompositionÞ >5_{irrotational part rotational part ËI 5½|TøC:ø} Ê óÂ]°_{O ¬ÎEyÌbã\&®Þ — A!Í@~"D} i®®ÞJ.mOãi«ºæ\&f.h®ÞÞÿ´

(5)

ÓÍ@~êÝÇ Xã\&fÝTø:øÝ5]Pñ ×_{ÊàyÞîaPi®®ÞTø/:ø)¿ÕÝóÂÿP} Í@_{~TàËÍÕ»?ÿPÝÑ@P¬ÕåaPi®Å(ìÝÎÎ} _{®ßÝ:ø¨éã9°Õ»|îÍÿP|b[2ÿaå:øÅ(Ý} Þ_îaPi®®Þ

1 Introduction

Because of the demand for security of constructions and the complex phenomenon in fluid dynamics, the research in wave interaction with structures, such as bridge piers and submerged dikes, is an important topic for numerical simulation.

Water wave passing a structure causes wave reflection and diffraction as well as the generation and shedding of vortices. Vortex shedding will lead to severe form drag force, and the wave height will be influenced by the energy dissipation and the vorticity near free surface. Owing to the complication of wave-vortex interaction and the consideration for the computational efficiency, it is difficult to simulate the water wave with vortex effect.

In this study, a grid-free numerical model using boundary integral equation method for the simulation of gravitational water wave with vortex effect is pro-posed. The computational algorithm is established by combining a potential flow field with a vortex flow field. In the first year, we develop the vortex separating mechanisms from rigid and free surface boundaries, and set up a two-dimensional vortex model. In the second year, we will improve the boundary conditions; deduce the boundary integral equations and calculating procedures to combine the two-dimensional potential model with vortex model. This hybrid model will be applied to the simulations of viscous water wave passing over a submerged dike.

The modern vortex method was born in the 1970s and the prominent investiga-tors involved in its early development are A. Chorin, A. Leonard, and C. Rehbach in France. Much interest in vortex methods during the early 1980s focused on math-ematical aspects such as the convergence properties (Hald et al., 1978; Hald, 1979; Beale and Majda, 1981, 1982a,b, 1985). Comprehensive reviews of the development of vortex methods and their applications can be found in Leonard (1980), Spalart (1988), Sarpkaya (1989) and Puckett (1993).

(6)

Recently, the difficulties of “N-body problem” for the particle vortex methods has been successfully addressed by the application of the Fast Multipole Method (Greengard and Rokhlin, 1987) for the calculation of the particle velocities, whereas some workers have bypassed the problem with mixed Eulerian-Lagrangian formula-tions such as the vortex-in-cell method (Couet et al., 1981; Rottman et al., 1987; Smith and Stansby, 1988; Brecht and Ferrante, 1990; Ebiana and Bartholomew, 1996; Liu and Doorly, 2000; Liu, 2001), although at the cost of adding interpolation errors.

2 Mathematical Formulations

2.1 General Formulations of Two-Dimensional Linear Wave

Problems

In the two-dimensional fluid domain V (t) with boundary S(t) (see Fig. 2 for definitions), the linearized Navier-Stokes equations for incompressible fluid with constant density ρ and kinematic viscosity ν are expressible in terms of the velocity u and the pressure p as

∂u ∂t = −

1

ρ∇p − gˆj + ν∇

2_u _{in V ,} ₍₁₎

where ˆj is the unit vector which points vertically upwards, and x = (x, y), u = (u, v). Together with the continuity equation

∇ · u = 0 in V , (2) the Navier-Stokes equations form a set of three scalar differential equations (for two-dimensional problems) sufficient for the determination of u and p, provided that adequate initial and boundary conditions are known.

It is, however, advantageous to introduce the concept of the vorticity ω = (0, 0, ω) defined by

ω= ∇ × u , (3)

(7)

∂ω ∂t = ν∇

2_ω ₍₄₎

obtained from equation (1) by taking the curl of each term in that equation and using equation (3). The set of equations (2) to (4), with u and ω as dependent variables, replaces the set of equations (1) and (2) in which u and p are dependent variables. There are several reasons for using ω in the formulation of the problem. The principal one is that the set of equations in terms of u and ω separate conveniently into two aspects: a dynamic aspect which deals with the change of the vorticity field ω with time described by equation (4), and a kinematic aspect which deals with the continuity of velocity field u governed by (2) and relates u at any instant to the vorticity field ω at that instant.

The linearized free surface boundary conditions are

1. Kinematic free surface boundary condition:

∂ζ

∂t = v ; (5)

2. Dynamic free surface boundary conditions, normal stress:

−p + 2µ∂v ∂y = −p

s_; ₍₆₎

3. Dynamic free surface boundary conditions, tangential stress:

∂u ∂y + ∂v ∂x = τs µ , (7)

where ζ(x, y, t) is the vertical displacement of the free surface from its unperturbed configuration consisting of the xy-plane, ps _{is the atmospheric pressure, τ}s _{is the}

wind stress, and µ is the dynamic viscosity. Because the vorticity ω = ∂v/∂x − ∂u/∂y, (7) can be rearranged to give the vorticity on the surface,

ω = 2∂v ∂x −

τs

(8)

2.2 The Application of the Helmholtz Decomposition

The-orem

The Helmholtz decomposition theorem states that, an arbitrary differentiable velocity field u can be decomposed into its irrotational part ue and rotational part

uv as

u = ue+ uv, (9)

and there exists a scalar function φ, called the scalar potential, and a vector function, Bv = (0, 0, Ψ) (for two-dimensional problems), called the vector potential, which

satisfies

ue = ∇φ , (10)

and

uv = ∇ × Bv. (11)

From the continuity equation in the incompressible flow, the scalar potential φ satisfies the Laplace equation

∇2φ = 0 . (12)

From the definition of vorticity, we get that Ψ satisfies the Poisson’s equation

∇2Ψ = −ω . (13)

Substituting the decomposition (9) to (11) into the linearized Navier-Stokes equa-tions, we get p ρ = − ∂φ ∂t − gy , (14) and ∂Ψ ∂t = ν∇ 2_{Ψ .} ₍₁₅₎

(9)

∂Ψ

∂t = −νω . (16)

Substituting for the surface vorticity from the shear stress boundary condition (8), at y = 0 the above equation becomes

∂Ψ ∂t = −2ν ∂v ∂x + τs ρ at y = 0 . (17) Also, we can substitute (5) into (17) and have

∂Ψ ∂t = −2ν ∂2_ζ ∂t∂x + τs ρ at y = 0 . (18) These are equivalent to the stress boundary condition, and can be regarded as the boundary condition of Ψ at the free surface.

Applying the decomposition (9) to (11), the KSFBC (5) can be expressed as

∂ζ ∂t = ∂φ ∂y − ∂Ψ ∂x at y = 0 . (19) and the normal stress boundary condition becomes

∂φ ∂t + gζ − 2ν ∂2_φ ∂x2 = − ps ρ + 2ν ∂2_Ψ ∂x∂y at y = 0 . (20) The results are summarized below. The problem requires to solve

∇2φ = 0 (21)

and the linearized vorticity equation

∂ω ∂t = ν∇

2_ω ₍₂₂₎

in the flow field, and at y = 0 the associated boundary conditions are

∂Ψ ∂t = −2ν ∂2_ζ ∂t∂x + τs ρ , (23) ∂ζ ∂t = ∂φ ∂y − ∂Ψ ∂x , (24) ∂φ ∂t + gζ − 2ν ∂2_φ ∂x2 = − ps ρ + 2ν ∂2_Ψ ∂x∂y. (25)

(10)

Now we focus on the special case where the wind stress τs _{is zero. From (8) the}

vorticity at the free surface is

ω = 2 ∂

2_φ

∂x∂y − 2 ∂2_Ψ

∂x2 , (26)

and (23) can be integrated with respect to time

Ψ = −2ν∂ζ

∂x + Ψ0 at y = 0 . (27) If the fluid was initially at rest, the constant of integration Ψ0 can be taken to be

zero; therefore

Ψ = −2ν∂ζ

∂x at y = 0 . (28) When these are substituted into (24), it becomes

∂ζ ∂t = ∂φ ∂y + 2ν ∂2_ζ ∂x2 at y = 0 . (29)

2.3 Integral Representations

2.3.1 Irrotational Flow Field

The governing equation of irrotational velocity field ue is the Laplace equation

(12). Applying the Green’s second identity, the integral representation of Laplace equation (12) can be obtained as

αφ(x) = Z

S

(φ0∇0G − G∇0φ0) n0 dS0, (30)

where the subscript “0” indicates a variable, or a differentiation, or an integration evaluated in the x0 space. α is the solid angle at x, n is the unit outward vector

normal to S, φ and G are finite and continuous and possess continuous first and second partial derivative in V . G is the fundamental solution of Laplace equation,

G =    − 1 2πln 1

r for two - dimensional problems

− 1

4πr for three - dimensional problems

, (31)

(11)

2.3.2 Rotational Flow Field

The second Green’s identity for vector states that, if vectors P and Q are single-valued, finite, and have continuous second derivatives, then

Z V [P · (∇ × ∇ × Q) − Q · (∇ × ∇ × P)] dV = Z S n · (Q × ∇ × P − P × ∇ × Q) dS , (32)

where S is the boundary of the fluid region V , n is the unit normal vector on S directed outward from V . If we let

Q = Bv and P = ∇ ×

a r

,

where a is a constant unit vector, we can get the integral presentation of the rota-tional part of velocity vector from equation (32) as (Wu and Thompson, 1973)

uv = Z V ω₀× ∇₀G dV₀− Z S [(uv 0· n0)∇0G − (uv 0× n0) × ∇0G] dS0, (33)

here the subscript “0” indicates that a variable or an integration is in the r0 space,

G is the fundamental solution of Laplace equation. The gradient of G is

∇0G = −

(x − x0)

2(d − 1)πrd, (34)

where d = 3 for three-dimensional problems and d = 2 for two-dimensional problems. Thus, the general expression for the rotational part of velocity is of the form

uv(x) = 1 2(d − 1)π Z V ω₀× (x − x₀) |x − x0|d dV0 − Z S (uv 0· n0)(x − x0) |x − x0|d dS0 + Z S (uv 0× n0) × (x − x0) |x − x0|d dS0 . (35)

The category of flows which were studied in Wu and Thompson (1973) is fluid in infinite domain bounded internally by a moving solid surface. In other words, besides the free stream velocity, the flow field is only driven by vorticity, and the boundary condition is simply the no-slip condition on the solid surface. In the

(12)

viscous wave problem, there exists an important potential effect – pressure gradient due to gravity waves, and the boundary conditions are more complicated. Owing to this, we have to formulate the relationship between the integral representation of potential flow and that of rotational flow.

In the present research, we find that the irrotational part of velocity can be expressed as ue(x) = − 1 4π∇ × Z S µ(x0)∇0 1 r × n0dS0, (36)

where µ(x0) is the strength of dipole located at x0. In other words,

ue(x) ≡ ∇ × Be, (37) where Be = − 1 4π Z S µ(x0)∇0 1 r × n0dS0. (38)

The above corollary indicates the result: the irrotational velocity field can also be expressed as a solenoidal field

ue= ∇ × Be,

i.e., the velocity field can be expressed as

u = ∇ × (Be+ Bv) = ∇ × B

without loss of generality. Therefore, we can let Q = B instead of Q = Bv in (32),

and get the integral representation of velocity field u as

u = Z V ω₀× ∇₀G dV₀− Z S [(u0· n0)∇0G − (u0× n0) × ∇0G] dS0. (39)

This indicates that if the vorticity field in the domain and the velocity on the boundaries are known, the velocity field in the domain can be evaluated from the integral equation (39).

(13)

3 The Vortex Methods

The potential φ at the free surface and the free-surface displacement ζ are solved by the time integrations of (25) and (29) respectively (ignoring the effects of wind stress and atmospheric pressure). The vorticity field ω is the solution of the vor-ticity equation (4), which is solved using the vortex method and the vortex sheet method. The potential field φ and the velocity field u are the solution of the integral representations (30) and (39) respectively, solved by the boundary integral method. The flow chart of computational scheme for two-dimensional linear wave prob-lems is shown in Fig. 1, and the brief introductions of numerical methods for vortex field simulation are shown in the following subsections.

3.1 Vorticity Layer Approximation

For convenience, in this subsection we only consider the vorticity-induced flowfield, i.e.

u = uv, ue = 0 .

For problems of external flow past finite bodies, we can consider the region V in Eq. (35) to be the entire (infinite) region occupied by the fluid. Then the boundary S is divided into two parts: the fluid–solid interface on which the no-slip condition (u = 0) applies and a surface infinitely remote from (and enclosing) the body on which the free-stream velocity boundary condition (u = u∞) applies. The surface

integrals in Eq. (35) can then be evaluated, giving

α(x) u(x) = 1 2(d − 1)π Z V ω₀× (x − x₀) |x − x0|d dV0+ u∞. (40)

The above equation can be recognized as the Biot-Savart law of induced velocities. Now we divide V into two regions: an interior (or exterior) region Vi located

away from S and a vorticity layer Vb located adjacent to S. In Vb, the vorticity

ω = ∂un ∂s −

∂us

∂n , (41)

where s and n is the tangential and normal coordinates in Vb. From the boundary

(14)

∂un ∂s ≪ ∂us ∂n in Vb, therefore, we have ω ≃ −∂us ∂n in Vb. (42)

From Eq. (72), we can get the tangential velocity on the interface of Vb and Vi is

us(s, n) = −

Z δ

0

ω(s, n) dn , (43) where δ is the thickness of Vb. Besides, the strength of vortex layer γ in Vb is

γ(s) = Z δ

0

ω(s, n) dn . (44) From Eqs. (43) and (44), we know that

us(s, δ) = −γ(s) . (45)

Otherwise, from Eq. (40) and the vorticity layer approximation, the velocity com-ponents are α ux(x, y) = − 1 2π Z Vi ω0(y − y0) (x − x0)2+ (y − y0)2 dV0 − 1 2π Z S+ γ0(y − y0) (x − x0)2+ (y − y0)2 dS0+ u ∞ x , (46) α uy(x, y) = 1 2π Z Vi ω0(x − x0) (x − x0)2+ (y − y0)2 dV0 + 1 2π Z S+ γ0(x − x0) (x − x0)2+ (y − y0)2 dS0+ u ∞ y , (47)

where S+ _{is the interface of V}

i and Vb.

3.2 Vortex Blob Representation

The vortex method is a particle method in which fluid particles carrying concentra-tions of vorticity are followed as their posiconcentra-tions and concentraconcentra-tions evolve with the motion of the fluid.

(15)

The vorticity field is represented by ω(x, t) = N X j=1 Γjγj[x − xj(t)] (48)

where γj is the vorticity distribution within the vortex located at xj with the

nor-malization

Z

V

γj(x)dV = 1. (49)

We assume that the differences in the distributions between the vortices depend only on a parameter σj, i.e., γj is given by

γj(x − xj) = 1 σ2 j f |x − xj| σj (50) where the shape or distribution function f is common to all vortex elements. The quantity σj is clearly a measure of the spread or core size of the vortex. The velocity

induced by the vorticity field is given by the Biot-Savart integral:

u(x, t) = − 1 2π N X j=1 (x − xj) × ˆezΓjg(|x − xj|/σj) |x − xj|2 (51) where g is given by g(y) = 2π Z y 0 f (z)z dz (52)

A number of core distribution functions been employed. For example the Gaus-sian distribution

γ(x) = 1

πσ2 exp(−|x|

2_/σ2_). ₍₅₃₎

3.3 Vorticity Creation and Shedding

Besides the diffusion of vorticity in the flow field, another important viscous effect is the vorticity creation at the boundary. At the free surface, amount of vorticity creation is determined by (26). At solid boundary, the creation of vorticity is to maintain the no-slip condition at the surface. While the vorticity field in the domain V (t) (except the boundary) is solved from the step 2 shown in Fig. 1, and

(16)

the vorticity distribution at the free surface is determined from (26), the vorticity distribution at solid boundary can be solved from (39).

4 Numerical Methods for Linear Wave Problems

The steps of computational scheme for two-dimensional linear wave problems (ignore the effects of wind stress and atmospheric pressure) is illustrated as a flow chart shown in Fig. 1, and is briefly described in the following subsections.

Free surface boundary conditions (25) and (29) are integrated at time t, to establish both the new wave height and the relevant boundary conditions on the free surface, at a subsequent time t + ∆t (with ∆t being a small time step). From the Taylor series expansion, we have

φ(t + ∆t) = φ(t) + ∂φ(t) ∂t ∆t + ∂2_φ(t) ∂t2 (∆t)2 2 + O(∆t) 3 ₍₅₄₎

for the potential, and

ζ(t + ∆t) = ζ(t) + ∂ζ(t) ∂t ∆t + ∂2_ζ(t) ∂t2 (∆t)2 2 + O(∆t) 3 ₍₅₅₎

for the wave height. The last terms in (54) and (55) represent truncation errors. The first Lagrangian time derivatives of φ and ζ can be computed from (25) and (29). The second time derivatives are

∂2_φ ∂t2 = −gζ,t+ 2ν ∂2_u ∂x∂t, (56) ∂2_ζ ∂t2 = ∂φ,t ∂y + 2ν ∂2_ζ ,t ∂x2 . (57) Here ∂2_u ∂x∂t = ∂2_φ ,t ∂x2 + ∂2_Ψ ,t ∂x∂y , φ,t = −gζ + 2ν ∂2_φ ∂x2 + ∂2_Ψ ∂x∂y , Ψ,t = −2ν ∂ζ,t ∂x .

(17)

5 Results and Discussions

5.1 Uniform flow over a rectangular obstacle

Here we show flow simulation results over a rectangular obstacle, to validate the rotational part of the present model. The rectangular obstacle of height 1.4 m and width 3 m was placed on the bottom. The inflow velocity is 1.0 m/s in the x direction.

Figs. 3 to 10 show the developed flow velocity field. The flow separates at the windward corner and clockwise recirculation eddies form and grow on the roof of the obstacle. Compared to the experimental result (Van Dyke, 1982), the flow pattern is reasonable.

5.2 Generation of surface waves by point vortex

A point vortex with constant strength Γ is moving at constant velocity U in the direction perpendicular to x axis and at a fixed distance below the undisturbed free surface, as shown in Fig. 12. The problem is to determine the surface waves generated by the point vortex. Utilizing this problem, we can validate the surface wave computation of the proposed model, by compared with the analytical solution.

5.2.1 Analytical Solution

A two-dimensional coordinate system moving with the vortex is used, with the y-axis pointing upward and the x-y-axis located in the undisturbed free surface, with x positive in the direction of the vortex motion. It is assumed that the fluid is inviscid. Formulating this problem in terms of the stream function ψ(x, y), one has that ψ satisfy the Laplace equation

ψxx+ ψyy = 0 for y ≤ ζ (x) (58)

everywhere in the fluid domain except at the point (0, −b) where ψ is required to be logarithmicly singular. The boundary conditions are:

1. Kinematic free surface boundary condition

(18)

2. Dynamic free surface boundary condition 1 2|∇ψ| 2 + gζ = 1 2U 2 _{on y = ζ (x)} ₍₆₀₎

3. Bottom boundary condition

~k × ∇ψ = −U~i as y → −∞ (61) where ~k = ~i × ~j.

4. Lateral boundary condition

~k × ∇ψ = −U~i as x → ∞ (62)

ζ (x) = 0 as x → ∞ (63)

and let

ψ (x, 0) = 0 as x → ∞ (64)

Far downstream there is a train of periodic waves with unknown wave length λ

ψ (x, y) = ψ (x + λ, y) as x → −∞ (65)

It is assumed that the point vortex is a weak disturbance so that in the neigh-borhood of the free surface the perturbations about the uniform flow are everywhere small, the stream function can be expanded by the regular perturbation method in terms of a perturbation parameter ǫ

ψ (x, y) = −U y + ψ(1)+ ψ(2)+ ψ(3)+ . . . (66) where ψ(n) _{is of O (ǫ}n_{). The free-surface elevation ζ (x) has the expansion}

(19)

and that the uniform-stream velocity U is an unknown of the problem with the expansion (Wehausen and Laitone, 1960).

U = u + u(1)+ u(2)+ . . . (68) where ζ(n) _{and u}(n) _{are both of O(ǫ}n₎

The first-order solution of the stream function is (Salvesen, 1969)

ψ(1) ₌ −Γ

2πRe

ln (z − ib) − ln (z + ib) + 2e−ikz

Z z ∞ eiku u + ibdu (69) where k = g/u2 _{is the wave number. If we define}

I (ζ) ≡ e−ζ E1(−ζ) ≡ e −ζ Z ∞ −ζ e−u u du (70)

where E1 is the exponential integral, (69) can be expressed as

ψ(1) = −Γ

2πRe [ln (z − ib) − ln (z + ib) + 2I {ik (z + ib)}] (71) Since

Re [ln (z − ib) − ln (z + ib) + 2I {ik (z + ib)}]_y=0 = Re [ln (x − ib) − ln (x + ib) + 2I {ik (x + ib)}] = 2Re [I {ik (x + ib)}]

the first-order wave elevation is

η(1) = −Γ

πuRe [I {ik (x + ib)}] (72) The uniform-stream velocity is

U = u 1 + 1 2k 2_α2 (73) where α is the first-order wave amplitude far downstream

α = 2Γ u e −kb (74) and lim x→−∞ζ (1)_{(x) = −α sin(kx).} ₍₇₅₎

(20)

5.2.2 Numerical results

The field length of numerical simulation was selected to be 17 m, and the water depth is 5 m (see Fig. 13). The strength of circulation Γ = π, and the velocity of point vortex is 3 m/s. The point vortex moved from (−6, −1.5) to (0, −1.5). The initial condition is given analytically. The downstream boundary conditions are assumed to be no mass flux. The boundary conditions at upstream boundary are given from the periodic wave (75).

Figs. 14 to 16 demonstrate time histories of water elevation at various points, where Point A is at x = −2.1, Point B at x = −0.1 and Point C at x = 1.9. It has been observed that the numerical results obtained by the proposed model agree with the the analytic solution very well. Figs. 17 and 18 are the velocity vector plots at time = 1.0 sec. and 1.7 sec., respectively.

5.3 Linear wave over a rectangular obstacle

The rectangular obstacle of height 0.125 m and width 0.25 m was placed on the bottom. The water depth is 0.25 m. Linear wave of amplutide 3.2 cm and period 1.5 sec. were generated. 180 panels with width 0.0125 m were placed on the boundary. Time interval of computation is 0.001 sec.

Figs. 19 to 28 demonstrate the instantaneous vortex particle paths and the corre-sponding velocity vectors. These figures show that when a wave crest arrives in front of the obstacle, vortices starts to form above the front corner of the obstacle. These vortices propagate to the downstream and form a larger eddy behind the obstacle, and interacting with the vortices formed from the back corner of the obstacle due to the change of velocity direction.

6 Conclusions

In this study, a grid-free numerical model combining potential and vortex flow fields for the simulating of gravitational water wave with vortex effect is proposed.

Two cases are applied to validate the proposed model: uniform flow over a rectangular obstacle, and the generation of surface waves by point vortex. Good agreements are found in the water elevation and vortex patterns between numerical and experimental/analytical results. The flowfield of vorticity shedding due to the

(21)

rectangular obstacle in the linear wave is simulated. These results show the ability of the proposed model for the simulations of free-surface flows with vortex effect.

References

Batchelor, G. K. An Introduction to Fluid Dynamics, Cambridge University Press, Cambridge, 1967.

Beale J. T. and A. Majda. Rates of convergence for viscous splitting of the Navier-Stokes equations. Math. Comp., 37, 243-259, 1981.

Beale J. T. and A. Majda. Vortex methods I: Convergence in three dimensions. Math. Comp., 39, 1-27, 1982.

Beale J. T. and A. Majda. Vortex methods II: High-order accuracy in two and three dimensions. Math. Comp., 39, 29-52, 1982.

Beale J. T. and A. Majda. High-order accurate vortex methods with explicit velocity kernels. J. Comp. Phys., 58, 188-208, 1985.

Brecht S. H. and J. R. Ferrante. Vortex-in-cell calculations in three dimensions. Comp. Phys. Comm, 58, 25-54, 1990.

Chorin, A. J. Vortex Sheet Approximation of Boundary Layers. J. Comp. Phys., 27, 428-442, 1978.

Cottet G.H. and P. Poncet. Advances in direct numerical simulations of 3D wall-bounded flows by vortex-in-cell methods. J. Comp. Phys., 193(1), 136-158, 2004.

Couet B., O. Buneman, and A. Leonard. Simulation of three-dimensional incom-pressible flows with a vortex-in-cell method. J. Comp. Phys., 39(2), 305-328, 1981.

Ebiana A. B. and R. W. Bartholomew. Design considerations for numerical filters used in vortex-in-cell algorithms. Comp. & Fluids, 25(1), 61-75, 1996.

Greengard L. and V. Rokhlin. A fast algorithm for particle simulations. J. Comp. Phys., 73, 325-348, 1987.

(22)

Hald O. and V. Mauceri Del Prete. Convergence of vortex methods for Euler’s equations. Math. Comp., 32, 791-801, 1978.

Hald O. Convergence of vortex methods for Euler’s equations II. SIAM J. Num. Anal., 16, 726-755, 1979.

Leonard A. Vortex methods for flow simulation. J. Comp. Phys., 37, 289-335, 1980.

Leonard A. Computing three-dimensional incompressible flows with vortex elements. Ann. Rev. Fluid Mech., 17, 523-559, 1985.

Liu, Chung Ho and Denis J. Doorly. Vortex particle-in-cell method for three-dimensional viscous unbounded flow computations. Int. J. Num. Meth. in Fluids, 32, 29-50, 2000.

Liu C. H. A three-dimensional vortex particle-in-cell method for vortex motions in the vicinity of a wall. Int. J. Num. Meth. in Fluids, 37, 501-523, 2001.

Puckett E. G. Vortex methods: An introduction and survey of selected research topics. In M. D. Gunzburger and R. A. Nicolaides, editors, Incompressible Com-putational Fluid Dynamics: Trends and Advances, pages 335-408. Cambridge Uni-versity Press, 1993.

Rottman J. W. , J. E. Simpson, and P. K. Stansby. The motion of a cylinder of fluid released from rest in a cross flow. J. Fluid Mech., 177, 307-337, 1987.

Spalart P. R. Vortex methods for separated flows. Technical Memorandum 100068, NASA, 1988.

Salvesen, N., On higher-order wave theory for submerged two-dimensional bodies. J, Fluid Mech., 38, 415-432, 1969.

Sarpkaya T. Computational methods with vortices. J. Fluids Eng., 11, 5-52, 1989.

P. A. Smith and P. K. Stansby. Impulsively started flow around a circular-cylinder by the vortex method. J. Fluid Mech., 194, 45-77, 1988.

(23)

Wehausen, J. V. and Laitone, E. V., Surface waves, Handbuch der Physik, V. 9, Berlin: Springer, 1960.

Wu, J. C. and Thompson, J. F. Numerical Solutions of Time-Dependent Incompress-ible Navier-Stokes Equations Using An Integral-Differential Formulation. Comput-ers and Fluids, 1, 197-215, 1973.

(24)

Input data and data preparation.

?

1. Time Integrations of φ,t and ζ,t.

Solve φ and ζ on the free surface.

?

2. Vorticity Transportation. Perform Euler convection and viscous diffusion of core radius

of all discrete vortices.

?

3. Vorticity Creation and Shedding. Evaluate the vorticity on boundaries. Shed the discrete vortices into flow field.

?

4. Potential Flow Computation. Solve the Laplace equations:

∇2_{φ = 0 and ∇}2_φ ,t = 0.

?

Advance time by ∆t.

(25)

S

_f

(t)

V(t)

S

_r1

S

b

S

_r2 t n x y z

Figure 2: Definition sketch for water wave problems.

U = 1.0, B = 3.0, H = 1.4, Re = 1000000 Time = 1.0 sec.

Figure 3: Numerical simulation of a uniform flow passing a rectangular dike (time = 1.0 sec)

(26)

U = 1.0, B = 3.0, H = 1.4, Re = 1000000 Time = 5.0 sec.

U = 1.0, B = 3.0, H = 1.4, Re = 1000000 Time = 10.0 sec.

(27)

U = 1.0, B = 3.0, H = 1.4, Re = 1000000 Time = 15.0 sec.

U = 1.0, B = 3.0, H = 1.4, Re = 1000000 Time = 20.0 sec.

(28)

U = 1.0, B = 3.0, H = 1.4, Re = 1000000 Time = 25.0 sec.

U = 1.0, B = 3.0, H = 1.4, Re = 1000000 Time = 30.0 sec.

(29)

U = 1.0, B = 3.0, H = 1.4, Re = 1000000 Time = 40 sec.

Figure 11: Experimental result of a uniform flow passing a rectangular dike (Van Dyke, 1982) x y U b ζ(x) Γ

(30)

x (m) ζ (m ) -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 -5 -4 -3 -2 -1 0 1 2 3 3.0

Figure 13: Setup of the simulation of surface wave generated by a point vortex

Time (sec) W a v e H e ig h t (m ) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 Numerical Analytic

(31)

Figure 15: Time histories of water elevation at point B

(32)

x (m) y (m ) -12 -10 -8 -6 -4 -2 0 2 4 -4 -2 0 2 _5.0

Figure 17: Velocity vector plot at time = 1.0 sec.

x (m) y (m ) -12 -10 -8 -6 -4 -2 0 2 4 -4 -2 0 2 5.0

(33)

x (m) y (m ) -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 3.0 t/T = 0

Figure 19: Wave over a rectangular dike: instantaneous particle paths and velocity vectors (t/T = 0) x (m) y (m ) -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 3.0 t/T = 1/10

Figure 20: Wave over a rectangular dike: instantaneous particle paths and velocity vectors (t/T = 1/10)

(34)

x (m) y (m ) -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 3.0 t/T = 2/10

Figure 21: Wave over a rectangular dike: instantaneous particle paths and velocity vectors (t/T = 2/10) x (m) y (m ) -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 3.0 t/T = 3/10

(35)

x (m) y (m ) -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 3.0 t/T = 4/10

(36)

x (m) y (m ) -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 3.0 t/T = 6/10

(37)

x (m) y (m ) -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 3.0 t/T = 8/10

含渦流作用之重力水波研究(2/2)

行政院國家科學委員會專題研究計畫 成果報告