行政院國家科學委員會專題研究計畫 期中進度報告
含渦流作用之重力水波研究(1/2)
計畫類別: 個別型計畫 計畫編號: NSC92-2611-E-002-030- 執行期間: 92 年 08 月 01 日至 93 年 07 月 31 日 執行單位: 國立臺灣大學水工試驗所 計畫主持人: 黃良雄 報告類型: 精簡報告 處理方式: 本計畫可公開查詢中 華 民 國 93 年 5 月 31 日
Å
Å
Åo
o
o»
»
»
I
I.
I
.
.
õ
õ
õº
º
ºÝ
Ý
ÝÞ
Þ
Þ@
@
@~
~
~
i
i
i
×
×
×
â
â
â:
:
:ø
ø
ø®
®
®à
à
à
¥
¥
¥æ
æ
æi
i
i®
®
®@
@
@~
~
~
(1/2)
The Study of Gravitational Water Wave
with Vortex Effect (1/2)
i_r
: NSC 92-2611-E-002-030Ƨ
: 92
O
8
`
1
^
93
O
7
`
31
^
x
x
x¹
¹
¹ß
ß
ß
:
?
?
?
.
.
.>
>
>0
0
0
»
»
Ȗ
ñ
ñ
È
È
È
.
.
.þ
þa
þ
a
a
.
.
.
Abstract
In this study, a grid-free numerical model for the simulation of gravitational water wave with vortex effect is proposed. The computational algorithm 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 vorticity 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.
Up to the present we have finished the establishment of mathematical formulations, the numerical models about potential flow and linear wave, and the numerical tests of vorticity generation from solid boundary. The unfinished parts are the numerical tests of vorticity generation from free surface, model integration, verifications and applica-tions.
Keywords: gravitational water wave, vortex, boundary integral equation method,
vor-tex method, numerical simulation
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 proposed. 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, and the propagating of wave groups.
2
Mathematical Formulations
2.1
General Formulations of Two-Dimensional Linear Wave
Prob-lems
In the two-dimensional fluid domain V (t) with boundary S(t), 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 + ν∇
2u in V , (1)
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)
and to consider the linearized vorticity equation
∂ω
∂t = ν∇
2ω (4)
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.
1. Kinematic free surface boundary condition:
∂ζ
∂t = v ; (5)
2. Dynamic free surface boundary conditions, normal stress:
−p + 2µ∂v
∂y = −p
s; (6)
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 config-uration 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 Theorem
The Helmholtz decomposition theorem states that, an arbitrary differentiable velocity field u can be decomposed into its irrotational part ue and rotational part uv asu = 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 equations, we get
p ρ = − ∂φ ∂t − gy , (14) and ∂Ψ ∂t = ν∇ 2Ψ . (15)
Using (13), (15) can be rewritten as
∂Ψ
∂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ω (22)
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)
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 FieldThe 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πln1r for two - dimensional problems
− 1
4πr for three - dimensional problems
, (31)
where r is the magnitude of the vector (x − x0).
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 rotational part of velocity vector from equation (32) as (Wu and Thompson, 1973)
uv = Z V ω0× ∇0G dV0− Z S [(uv0· n0)∇0G − (uv0× 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 ω0× (x − x0) |x − x0|d dV0− Z S (uv0· n0)(x − x0) |x − x0|d dS0 + Z S (uv0× 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 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
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 ω0× ∇0G dV0− 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).
3
Numerical 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 vorticity 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 problems is shown in Fig. 1, and the brief introductions of numerical methods for vortex field simulation are shown in the following subsections.
3.1
The Vortex Method
The vortex method is a particle method in which fluid particles carrying concentrations of vorticity are followed as their positions and concentrations evolve with the motion of the fluid. Let
ω(x, t) =X
j
fσ(x − xj(t))Γj (40)
be an approximation to the vorticity at some time t = k∆t. The jth term on the right hand side of (40) is referred to as the jth vortex element or jth vortex blob, xj(t) is its position at
time t and Γj its strength. fσ is the vorticity distribution function. Here we use the Gaussian
distribution function fσ(x − xj) = 1 πσ2 j exp µ −|x − xj| 2 σ2 j ¶ , (41)
where the quantity σj is clearly a measure of the core size of the vortex element.
3.2
The Vortex Sheet Method
Near solid and free-surface boundaries, vortex blobs do not provide good approximations of the vorticity because the flow has a large velocity gradient in the normal direction that is not well-captured by the radially symmetric blob representation. A better approximation of
the vorticity in these regions consists of a collection of smoothed one-dimensional segments parallel to the boundary called vortex sheets.
In the vortex sheet method the vorticity at time t = k∆t is approximated by a sum of linear concentrations of vorticity
ω(x, t) =X
j
ωjbl(x − xj(t))δ(yj(t) − y) , (42)
where (x, y) denote coordinates which are parallel and perpendicular to the boundary re-spectively. Each term of the sum in (42) is referred to as a vortex sheet. The jth sheet has center (xj, yj) and strength ωj. Here δ is the Dirac delta function, and bl = b(x/l) is what
we refer to as the smoothing or cutoff function in analogy with the vortex method, and l is the length of sheet. The most commonly used cutoff is the hat or tent function originally proposed by Chorin (1978)
b(x) =
(
1 − |x| |x| ≤ 1,
0 otherwise. (43)
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 the vorticity distribution at the free surface is determined from (26), the vorticity distribution at solid boundary can be solved from (39). Fig. 2 is the velocity field of vortex sheets, which are just shedding from the flat plate. Fig. 3 is the numerical result of a uniform flow passing a rectangular dike. Compared to the experimental result (Van Dyke, 1982), the flow pattern seems to be reasonable.
References
Batchelor, G. K. (1967), An Introduction to Fluid Dynamics, Cambridge University Press, Cambridge.
Chorin, A. J. (1978), “Vortex Sheet Approximation of Boundary Layers”, J. Comp. Phys., 27, 428-442.
Van Dyke, M. (1982), An Album of Fluid Motion, The Parabolic Press, California.
Wu, J. C. and Thompson, J. F. (1973), “Numerical Solutions of Time-Dependent Incom-pressible Navier-Stokes Equations Using An Integral-Differential Formulation”, Computers
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.
¾
x (m) y (m ) -1 -0.5 0 0.5 1 0 0.001 0.002 0.003 0.004
Figure 2: Vortex sheets created from the solid boundary
U = 1.0, B = 3.0, H = 1.4, Re = 1000000 Time = 40 sec.
Figure 3: Numerical simulation of a uniform flow passing a rectangular dike