二維水深積分模式應用於三維亞臨界分流分析

(1)

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

二維水深積分模式應用於三維亞臨界分流分析

計畫類別：個別型計畫

計畫編號： NSC92-2211-E-006-030-

執行期間： 92 年 08 月 01 日至 93 年 07 月 31 日

執行單位：國立成功大學水利及海洋工程學系（所）

計畫主持人：唐啟釗

報告類型：精簡報告

處理方式：本計畫可公開查詢

中華民國 93 年 11 月 9 日

(2)

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

二維水深積分模式應用於三維亞臨界分流分析

Application of two-dimensional depth-integral model on

three-dimensional sub-critical dividing flow analysis

計畫編號：NSC92-2211-E006-030

執行期限： 92 年 8 月 1 日至 93 年 7 月 31 日

主持人：唐啟釗國立成功大學水利及海洋工程系

電子信箱:

cjtang@mail.ncku.edu.tw

一、中文摘要

明渠分流現象在實際工程應用非常普遍且

重要。實際應用於流況分析時，常應用一維斷面

平均模式，或二維水深積分模式來模擬；自然

的，後者之計算結果優於前者。然而傳統水深積

分模式假設了「靜水壓分佈」與「忽略垂直向加

速度效應」，限制了分析較短波長流速變化與垂

直（三維）流動的流況。本計畫擬針對這些假設

之缺失，提出修正方式，擴充二維水深積分模式

之應用限制，以分析三維亞臨界分流現象。數值

模式之建立將先以垂直座標之羃次級數展開三

維流場之速度與壓力，並由代入動量與連續方程

式以獲得各級數（二維）係數之回歸關係偏維分

方程組。此關係方程組之維度將比欲解析之流場

減少一維（垂直向），故可將三維流場簡化成二

維模式組，其方程式個數與裁切之級數項數相

等。在使用數值方法求解各方程組後，即可求得

三維亞臨界分流流場。本計畫亦進行實驗工作，

實驗量測工作將在淡江大學水環系進行，（因協

同主持人許中杰教授已建立良好之分流實驗設

備，可節省經費與設置時間）。斷面速度與壓力

的細化量測，將與目前發展之數值模式相互驗證

比較，以確認模式之正確性。

關鍵詞：亞臨界流、分流、明渠流、水深積分式、

多項式級數解。

Abstract

The use of diving flows in open channel

system is common and important for practical

hydraulic engineering application. In practice, the

popular models to simulate such flow motion are

the one-dimensional section-averaged equation and

the two-dimensional depth-integral equations,

between which the latter undoubtedly obtains the

better result than the former. It is also known that

the assumptions of the hydrostatic pressure

distribution with negligible vertical acceleration are

made in these traditional integral-equation models.

The assumptions primarily restrict to the utilization

of model for flows with rapid variation in velocity

within the length scale shorter than the grid size

(which is equivalent to a half of the wavelength to

be resolved for all variables) and for the

three-dimensional flows with significant vertical

momentum transport. The present project is

proposed to release those assumptions through

modification of the 2D depth-integral model

equations and extending it to apply for the 3D

sub-critical diving flow in open channel. We will

expand the 3D velocity and pressure by using

truncated power series in vertical coordinates with

which the recursive formulations for those (2D)

coefficients can be obtained from the substitution

into the governing equations and boundary

conditions. The reduction of one-dimensionality in

the modified model due to using series expansion

method would simplify the 3D flow modeling to a

two-dimension equation set. The number of this

equation set is right equal to the truncated order of

the power series. Using numerical method to find

out the flow solution, one can then solve the 3D

sub-critical diving flow in the open channel system

as desired.

Keywords: sub-critical flow, dividing flow, open

(3)

expansion.

二、緣由與目的

本研究延續主持人之前國科會專題計畫，如

波統合理論之發展與應用」等之研究成果，延伸

於明渠流動現象之分析。本年度計畫主要以原有

之波場致流（以波動為主）之數值模式，擴充至

明渠亞臨界分流（以流動為主）之模式發展與應

用。雖然兩模式有部份類同，但因考慮之現象與

其衍生之流動機制極不相同，故仍具甚大之挑戰

性。本研究係在目前已有之研究成果，進一步架

構更精微之分析研究方法，並依理論模式的需

求，比對更精密的量測工作以供比較，期能對此

複雜現象，提供方便學理應用的分析工具。

明渠分流現象在實際渠流流動亦是非常普

遍且重要，工程應用分流之目的，主要用於控制

輸送水體與分配水量等之應用，如灌溉、導洪、

截流、越域引水、工業冷卻配水、污水處理、發

電、給排水、航行等等，無論在都市或郊區都可

見其應用，十分廣泛。分流(dividing flow)與匯流

(junction flow 或 combining flow)兩者均屬歧流

(branch flow)之渠流現象，但分流與匯流流動特

性冏異：分流常在主支渠交會處，引致強烈流離

(flow separation)與產生側向二次流動(secondary

flow)現象，不但常使交會點附近產生極為複雜的

三維流況，且因渦動流況局部影響使水面急劇變

化，亦阻塞通水斷面而形成局部束縮流動，使主

支渠分配流量之機制亦受其影響而變化。尤其當

較大交角（

θ

）之分流，上述之流離與二次流況

之影響甚為明顯。亞臨界流之分流產生之流離常

出現於分流支渠進口偏上游處；而匯流產生之流

離現象則出現於與支渠交會點上游之主渠對側

壁端，及交會點下游主渠邊壁處。若渠道中出現

超臨界流，則發生在歧流點的流況將更複雜，分

離與二次流流動極可能伴隨水躍現象而形成交

替水深變化的複雜流況情形，或亦可能仍然維持

超臨界流況，端視 Froude 數與主支渠寬，分流

交角配置等而定。由亞臨界流之分流實驗可觀察

到：此兩渠交會處下游端，將形成不連續之水位

變化，較強烈的流離常使水面產生急遽塌陷，由

近壁水面處強烈的渦流，逐漸向支渠下游與垂直

向渠底延伸，或擴張或束縮。因垂直向（水深方

向）之加速度效應與水平速度梯度動量的交換，

可能形成非穩態的分離區；且因非均勻之非線性

慣性動量分佈，而產生局部動壓力作用，使流離

區之封閉泡狀範圍不但非為垂直向下柱狀之水

平渦流結構，且亦因垂直向渦狀動量交換，而構

成二次流動出現在渦流區內即外流場區。後者肇

因於流入支渠的水體被束縮加速，但因流離泡區

曲面形狀在垂直與水平向都不均勻，故沿此曲面

流動之水體受離心力與本身慣性交互作用而呈

現不均勻加速分佈。

早期針歧流之研究多為物理模型實驗，一維

理論或數值分析為主要方法。1944 年 Taylor 可能

是第一個應用解析模式分析歧流況之學者，他依

此推得此模式獲得交匯點上游之水深；他亦使用

4 吋寬之渠道進行實驗，最終結論認為：直接使

用（一維）解析方法來處理動量方程式，以求解

分流流況的方式為不可行（Hsu 等，2002）

。Grace

與 Priest（1958）使用 5 吋寬之渠道來觀察不同

交角之分流及渠道寬度比，並利用圖形表示了無

因次化的結果。他們將分流區分為：存在與不存

在局部駐波兩類流況。當無駐波時，上下游水深

比與支渠下游水深比均近於 1，且這些比率將隨

支主渠流量比之增加而呈緩和減少。Modi 等

（1981）使用非黏性流體與無旋流的假設，利用

複變函數保角變換方法計算流離區大小，但他們

所得結果較實驗觀測值為大。Best 與 Reid（1984）

兩人考量不同交角與上下游流量比之效應進行

分析，並定義流離區形狀指數（separation zone

shape index）為流離區水平長度與渠寬之比。

Ramamurthy 等（1988, 1990）量測兩渠交會處之

壓差來說明側邊動量傳遞情形。Hsu 等（1998a，

1998b，2002）之系列實驗亦獲得不同交角之交

匯流最大束縮流況之束縮係數等關係，與其他學

者實驗觀測結果符合。這些研究都是亞臨界流

況，可能近於非線性拋物線形或橢圓型態偏微分

方程組的特性，而有別於超臨界流況的雙曲線型

態偏微分方程組特性，或更複雜的穿臨界流流

況。這顯然與波傳現象有很大不同，而亦有別於

主持人近年之淺水長波研究。

計算和分析分流或匯流流況的數值模擬研

(4)

究，大抵以一維流動模式（如 Hanger, 1989 等），

或二維水深模式為主（如張正熙，1995 等），較

特別值得一提的是，愛荷華大學水利研究所(IIHR)

的研究群如 Huang 等（2002）應用 Wilcox 發展

之紊動動能（k）與比耗率（

ω ）模式進行三維

匯流數值模擬，這是計算流體力學實際應用水利

工程問題的一例。他們宣稱可獲得三個量階誤差

以內的精度，但在他們應用 SGI Indigo 工作站

（R8000，單一 75MHZ 時脈的 CPU 平台）於計算

主支渠同寬，直角匯流的計算例中，用了 18 萬

格點與瞬時移動格網技巧計算三維流速分量與

壓力，共計使用約 230 CPU 小時才獲得結果，(雖

然文章中描述非穩態的計算模式但並未見討論

流場之瞬變效應，極可能僅是穩態的計算結

果？)，顯然並非很有效率。以個人淺見認為：

使用上述之一維斷面平均模式或者是二維傳統

水深模式，去分析三維結構的流離現象，當然有

很大精度的問題；而另一方面，使用三維 K-ω 紊

流模式去分析，雖能獲得某程度精度的改善卻又

太過複雜，而不能符合實際工程應用的需要。到

底是選用過度簡化而不夠精確的模式，或是需用

太通用但很繁複模式，來分析這種不具有十分複

雜幾何形狀的實際工程問題，這種兩難狀態，其

實也是目前計算流體力學（CFD）學者共同面臨

的困境：既要能獲得計算結果的精度與權衡計算

工作效率。顯然兼顧兩者的合理性，方能符合工

程實用的目標。

三、理論推導

定義水平向卡氏座標為

(

x

,

y

)

，垂直向上為

z ，時間為

t

。設未擾動之自由水面表為

z

=

0 ，

擾動後水面表為

z

=

Z

f

( , , )

x y t

；底床高程表為

( , )

b

z

=

Z x y

，因此，任何時間之水深為

( , , )

f

( , , )

b

( , )

h x y t

≡

Z

x y t

−

Z x y

。若考慮水體為密

度均勻分佈，水平速度向量表為

u

=

( ,

u

v

)

，垂

直速度分量為

w

，相關水平向之梯度為

1

(

x

,

y

)

∇ = ∂ ∂ ∂ ∂ ， p 為流體之壓力。為了分析

方便，將所有變數以下列參考尺度無因次化：參

考長度為特性水深

h （如可取為無擾動之定水

₀

深，

Z

b

= − =

h

₀

const

）；參考時間為

h

0

/

g

；參

考速度為

gh ；參考壓力為

₀

ρ

gh

0

；其餘變量之

參考量可由上述參考尺度導出。因此將連續方程

式與水平動量方程式對水深積分（即沿 z 向由

Z

_b

積至

Z ）

_f

，帶入相關液面與底床邊界條件後，可

寫成

1

:

t

(

)

0 C h

+ ∇ ⋅

h

u

=

, (1)

( )

[

]

[

]

1 1 1 f f f

:

(

)

(

)

(

)

(

)

(

)

t z b z f

h

h p

Z

p

+ ∇ ⋅

= −∇ ⋅

−

+

−

− ∇

⋅

−

x xx x x xx

u

uu

I σ

σ

I

σ

M

(2)

式中

z z

σ

zz

⎡

⎤

= ⎢

⎥

⎣

⎦

xx x x

σ

,

xx xy yx yy

σ

⎡

⎤

= ⎢

⎥

⎣

⎦

xx

σ

,

z xz yz

σ

⎡

⎤

= ⎢ ⎥

⎣

⎦

x

σ

各

為

3 3×

黏性應力張量，

2 2×

水平黏性應力張量，

與水平面上黏性剪應力張量(

2 1×

向量)。

(

T

)

zx

=

xz

σ

為垂直面上水平黏性剪應力張量(

1 2×

向量)，

σ

_zz

為水平面上黏性正向應力，

p ,

_f

p 為

₀

液面與床底壓力。式(1)與式(2)所有含上標線者

為水深平均量，如

u 為水平向水深平均流速，p 為

水深平均壓力等；下標 b 和 f 各代表在液面與底

床之邊界值。至此所得的關係式(1)(2)並未介入

任何假設，因此應可適用於任何符合動量與質量

守恆定律之不可壓縮流動，當然也包括二維或三

維的黏性流場。惟因速度與壓力在垂直向的分佈

尚未未知，因此無法直接應用於三維流場的計算

模擬。式(1)為各模式之通式，普遍適用；式(2)

則各有不同。式(2)等號右方含壓力與各種黏性

應力的水深平均值，和液面與底床邊界值等項，

若忽略所有正向黏性應力

σ

_xx

, (

σ

_xx

) , (

_b

σ

_xx

)

_f

，及

液面的黏性切向應力 (

σ

_x_z

)

_f

效應，再使用靜水壓

分佈

/ 2

p

=

ρ

gh

，

p

₀

=

ρ

gh

和

p

_f

=

0 ，

(3)

並使用

uu u u (垂直向均勻速度分布或忽略其

=

高階差項)與帶入式(1)後，式(2)即為傳統的水深

平均動量方程式

1 1

1 (

)

t

Z

f z b

h

+ ⋅∇ + ∇

=

_x

u

σ

(4)

但更一般化的壓力分佈須由垂直向動量方程式

積分得到，表示為

f f 2 2 f f f f z

( , , )

(

)

(

) (

)

(5)

zz Z Z zz _z z

p

z t

p

g Z

z

w

w dz

dz

t

ρ

=

+

− +

−

∂

⎛

⎞

+

_⎜

+ ∇ ⋅

_⎟

+

∇ ⋅

∂

⎝

⎠

∫

x

σ

u

σ

(5)

為了簡便應用，上式常用

p

_f

=

(

σ

_zz

)

_f

=

σ

_zz

= 以

0 忽略表面張力，不均勻氣壓分佈液面空氣正向應

力、與風剪力等效應。由此式可得底床與水深平

均壓力為

(

)

(

)

b z b z b

hw

p

gh

h

w

Z

t

ρ

∂

⎡

ρ

⎤

=

+

+ ∇ ⋅

_⎣

−

_⎦

− ∇ ⋅

∂

u

σ

x

σ

x

(6)

( )

(

)

2 2

2

z

gh

hp

hzw

h

z w

z

h w

t

∂

ρ

∂

⎡

⎤

=

+

+ ∇ ⋅

_⎣

u

−

σ

_x

_⎦

−

(7)

式中

w

可由連續方程式積分求得：

( , , )

b z h

w

z t

w

dz

−

=

−

∫

∇ ⋅

x

u

(8)

可獲得較傳統完整之水深積分式。假設級數展開

0

, , }( , , )

,

}( , )

N n n n n n

w p

η

t

w p

t

η

=

∑

{u

x

{u

x

, (9)

式中

η

= −

(

z

Z

_b

) /

h

，因此 0

≤ ≤ 。將壓力級數

η

1 再分解為靜壓

ρ

gh

(1

−

η

)

與動壓級數

n n

d

ρ

∑

η

之

和。帶入 Euler 方程組(黏性項可於獲得級數解後

再帶入應力應變率關係式中由數值疊代求收斂

解)可得各級數之係數如下：

1 1

1 ,

n n n

n

h

w

h

n

−

n

−

=

u

⋅∇ − ∇ ⋅

u

w

_o

= , (10)

0

1 1 1 1 1 ( ) n n n n n n h n h d a w w h n t − − n − t − ∂ − ∂ ⎛ ⎞ ⎛ ⎞ = + _⎜ + ∇ ⋅ _⎟− _⎜ + ⋅∇ _⎟ ∂ ∂ ⎝ c ⎠ ⎝ c ⎠

,

1 ≤ ≤

n

2 N

+

2 ,

(11)

2 2 0 0 1 N n n

p

gh

d

ρ

+ =

=

+

∑

= −

(12)

式中

a

_n

∑

> < = − < >

∆

−

=

/2 0 2 2 /

2

n k n n k n k

w

, (13)

/ 2 / 2 / 2 0

(

)

n n k n k n k k n n n k

w

< > − − < > < > =

=

∑

+

− ∆

c

u

(14)

/ 2

Int( / 2),

n

<

>=

∆ = −

n

1 Mod( , 2)

n

；Int 與

Mod 為 Fortran 整數函數。至於 u 級數須滿足

0

(

0

)

n n n n

n

h

t

h

t

∂

₊

_{⋅∇ −}

⎛

∂

_{∇ ⋅}

⎞

₌

⎜

⎟

∂

⎝

∂

⎠

u

+

u

S

(15)

1

)

(

1)

n k n k n k k k n n

n k

h

h k

n

d

h

− − =

⎡

−

⎤

= −

_⎢

⋅∇

−

∇ ⋅

_⎥

+

⎣

⎦

+ ∇ −

∇

∑

n

S

u

( u

(16)

至於剪應力可依據牛頓定律可由剪應變率級數

2 0 1 0 ( 1) N n n n n n n n h Z h h +

η

=

+

⎧

⎫

∇ = ∇ ⋅

_⎨

∇ − ∇

−

∇

_⎬

⎩

⎭

∑

u u u u

, (17)

1 1 0 1 1 ( 1) Re N n z n n n z h

η

− + =

∂

=

+

∂

∑

x u σ u

數值計算時，在每一時階中，依序先計算式(12a),

(10), (13), (14), (11), (12)與(16), (15)，離散式採用

二階中央差分法來表示時空微分項，非線性係數

使用兩時階平均法，使離散關係完全滿足二階

2

(

,

)

O

∆ ∆ ∆

x

y

t

之精度。因存在非線性耦合項，故

每時階均需疊代計算已獲得時間精確解。因此本

模式可應用穩流與非穩流況計算。為了加速收

斂，使用斑馬線是交互逐線疊代，併用超鬆弛技

巧。本模式除了使用單機 PC-Windows XP 平台

外，尚使用本系 8 節點 PC-Cluster on Linux 與國

網中心 PC-Cluster on Linux 平台系統。

四、結果與討論

圖一為測試一初始孤立波(波峰值為 0.5 倍靜水

深

h 由

₀

x

=

10 h

₀

開始傳遞)經過不同主支渠寬鄰

近分流點

x

=

35 h

₀

處瞬時速度場與徑線分佈圖。

此圖(a)顯示主渠寬與水深比

B h

/

₀

= 與支渠寬與

1 水深比

b h

/

₀

= 之流況；圖 (b) 為

1 B h

/

₀

= 但

1

0

/

2 b h

= ；圖(c)為

B h

/

₀

= 和

2 b h

/

₀

= 。這些計

1 算用以測試格網與時階的尺度，由孤立波在直渠

道的傳遞特性(本結果前段本報告未顯示)以大

致了解計算參數的合適性。因為級數展開的第零

階模式簡化可獲得緩變水深如 Peregrine(1967)的

長波方程組，或定水深之 Boussinesq 方程組，如

前所述，這可由傳統的水深模式回復非零的暫態

(transient)動壓力或垂直加速度部份，因此也可以

作為發展程式檢驗之用。(本文忽略底床變化效

應。)

壓力場之分佈取至二次冪次(即拋物線分佈)，其

結果如圖二所示。剪力分佈如圖三所示。本文計

算成果仍需再仔細評估。

此法與三維模式之計算效益評估：

設水平平面格網需 100 x 100=1 萬格點以計算變

數，三維模式需再於三維各個格點處計算 4 個變

(6)

數(3 個速度分量 , ,

u v w ，壓力 p 計算於流場內

部，水深方向解析度取 10 格點)；另外平面二維

格點計算液面高（水深

h

），因此需解 41 萬方程

式組內外迴圈疊代至收斂。反之，因垂直向速度

w

與壓力 p 直接由回歸代數式計算，無須疊代，

只需先計算 ,

u v 即可。若水平面每格點速度 ,

u v 級

數取十項（應該不需要，但保守估計使與 3D 流

況變數量相當，乍視需算 20 萬方程組疊代至收

斂，但實情是每階級數之係數方程式 ,

u v 各自偶

合於非線性傳遞項必須內外迴圈疊代（2 萬點），

而其他階級數之 ,

u v 為源項(18 萬點)，無須內迴

圈疊代只需外迴圈疊代，水深

h

亦是只需外迴圈

疊代，因此僅為 3D 程式計算時間約十分之一以

內，可知其效率。附帶一提，傳統水深方程則相

當為本模式級數只取一項之情形，當然節省另外

十分之一計算時間，但精度可能太差而難以接

受。本計畫工作將作一比較以具體說明此觀點之

正確性。

X Y 34 35 36 37 38 39 -2 -1 0 1 2

(a)

B h

/

₀

= 與

1 b h

/

₀

= 之流況

1

X Y 34 35 36 37 38 39 -2 -1 0 1 2

(b)

B h

/

₀

= 與

1 b h

/

₀

= 之流況

2

X Y 34 35 36 37 38 39 -2 -1 0 1 2

(c)

B h

/

₀

= 與

2 b h

/

₀

= 之流況

1 圖一二維水深模式分流場瞬時速度與徑線圖

（主渠寬與水深比

B h ；支渠寬與水深比

/

₀

b h ）

/

₀

五、參考文獻：

2 Best, JL, and Reid, I (1984). ‘‘Separation zone

at open-channel junctions,’’ J Hydrau Eng,,

110(11), 1588–1594.

3 Grace, JL, and Priest, MS (1958) ‘‘Division of

flow in open channel junctions,’’ Bulletin No.

31, Engineering Experimental Station, Alabama

Polytechnic Institute.

4 Hager, WH. (1983). ‘‘An approximate treatment

of flow in branches and bends,’’ Proc Inst Mech

Engrs London, UK, 198C(4), pp63-69.

5 Hager, WH. (1989) ‘‘Transitional flow in

channel junctions,’’ J Hydrau Eng, 115(2),

pp243–259.

(7)

6 Hager, WH. (1992) ‘‘Discussion of ‘Dividing

flow in open channels,’ By A. S. Ramamurthy,

D. M. Tran, and L. B. Carballada,’’ J Hydrau.

Eng, 118~4, 634–637.

7 Hsu, CC., Wu, FS., and Lee, WJ. (1998a)

‘‘Flow at 90° equal-width open-channel

junction,’’ J Hydraul. Eng., 124(2), 186–191.

8 Hsu, CC, Lee, WJ., and Chang, CH. (1998b)

‘‘Subcritical open-channel junction flow,’’ J

Hydrau Eng, 124(8), 847–855.

9 Huang, J. (2000) ‘‘Development and validation

of a three-dimensional numerical model for

application to river flow,’’ PhD thesis, Univ. of

Iowa, Iowa.

10 Modi, PN., Ariel, PD., and Dandekar, MM.

(1981) ‘‘Conformal Mapping for channel

junction flow,’’ J Hydrau Div, ASCE, 107(12),

1713–1733.

11 Ramamurthy, AS., Carballada, LB., and Tran,

DM (1988) ‘‘Combining open-channel flow at

right-angled junctions,’’ J Hydrau Eng, 114(12),

1449–1460.

12 Ramamurthy, A. S., and Satish, M. G., (1988)

‘‘Division of flow in short open channel

branches.’’ J. Hydrau Eng, 114, 4, 428–438.

13 Ramamurthy, AS., Tran, DM., and Carballada,

LB (1990) ‘‘Dividing flow in open channels.’’ J.

Hydrau Eng, 116, 3, 449–455.

14 Taylor

, E. H., (1944) ‘‘Flow characteristics at

rectangular open-channel junctions,’’ Trans. Am.

Soc. Civ. Eng., 109, 893–912.

15 Wilcox, D. C. (1993). “Turbulence modeling for

(8)

附錄: 本計劃部分成果

(9)

1.國立成功大學水利及海洋工程學系副教授。 2 國立成功大學水工所技術員。

孤立波遇斜坡棚台崩解現象之數值研究

唐啟釗

1

鄭志維

2

摘要

本文以數值模擬孤立波經斜坡(坡度 s=0.1)到較淺等深水域時，因斜坡底床作用引致孤立波崩裂而分解成數個不等高度之孤立波列與震盪尾波現象。此現象可由 Pregrine(1967) 發展之通用 Boussinesq 方程組長波模式描述之，數值分析採用二階時空精度之中央差分與兩時階平均方法離散後，疊代求得時間精確解。因孤立波沿底床斜坡而增加垂直向加速度非線性效應，產生額外的頻散而改變波形，本文以數值觀點，量化相關波傳動力性質，以觀察物理現象過程細部結構。當孤立波在較淺水域崩解分裂後，重組非線性與頻散平衡關係，由透射波列、震盪尾波群及反射波列等之波能與質量之暫態分配，了解非線性淺化互制作用。

Numerical analysis on solitary waves scattering

by a sloping shelf

Chii-Jau Tang

1

Chih-Wai Chang

2

ABSTRACT

This study attemps to analyze the soliton disintegration process during a solitary wave traveling over a sloping shelf of smaller constant water depth by numerical calculation. Using a numerical method of second-order central difference in space and time with time-averaging scheme, we can discretize the generalized Boussinesq equations (gB, Peregrine in 1967) to seek the numerical solution through time-accurate iteration. The action of the plane slope introduce not only the nonolinear convection in the vertical direction but dispersion to deform the waveform in the course of interaction. The detailed structure in dynamic wave motion can be depicted quantitatively along with time in the present study. The evolved distribution of wave energy and mass transport by the transmitted, oscillatory trailing and reflected wave pattern over slope does reasonally respond to the inherent reconstruction of nonlinearity and dispersion of these fission solitons.

一、緒論

波浪自深海中向岸傳遞時，隨著水深變淺，其非線性與頻散效應逐漸增加，波浪變形引發之現象對海岸地區影響甚鉅；另外，因颱風產生之暴潮、地震發生而產生的海嘯和受月球引力影響的潮汐等長波現象亦是威脅海岸結構物安全性重要因素。因此非線性長波淺化的研究不但具有學理價值，且對工程中港灣、河口的設計以致於海灘的防護開發有其必要性。長波在等水深水域傳遞時，若其弱非線性 ( 用微擾參數 0/ 0 1 A h

ε =

表示，式中 0 A 與 0 h 為波幅與靜水深)與頻散性(用微擾參數 2 0 (kh) 1

µ =

表示，

_k

為波數)相互平衡，則常能保持恆久波形。典型的例子如：單一

_sech

2_{波形之孤立波與具有空間}

周期性的 cn 波，兩者均為 KdV (Korteweq & de Vries, 1895)方程式 的正合解。此兩種具恆久波形特性的長波均只能單向運行；波速大小隨其波幅高低而定；兩波追撞互制後仍能各自保持原有波形 而僅有相位移遞。除了 KdV 模式描述的水波外，尚有更一般性的 Boussinesq (1872)方程組長波模式，可允許雙向波動而仍保有上述長波特性。Boussinesq 方程組長波模式可由非常廣泛形式的族群構成，以不同微分項與二階微擾階次形式出現於頻散效應項中。這些形式通常以水平流速(或速度勢)及水波高兩變數構成聯立系統，其中水平流速可為斷面平均量、底床量、水面量、甚或特定水深量(Wu, 2001)。不同的微擾展開使其符合微小振幅線性波的完整或近似頻散關係，亦可獲得包含更高微分項的 Boussinesq 模式。雖然這些模式都不存在封閉型正合解析解，但可用數值計算傳統二階 Boussinesq 方程的積分式 (Tseng, 1997)或以近似解析法(Yih & Wu, 1985)求得互制過程演變的波形。較有名的 Boussinesq 方程組除了 Boussinesq 本人發展的原始形式外，尚可因考慮底床變化對長波變形或更高階效應的研究，如 Green(1837)、Airy(1845)、Lin & Clark (1959)、Long (1964)、Mei & Mehaute (1966)、Pregrine (1967)、Benjamin(1972)、 Wu (1981)、Nwogu (1993)、Wei (1995)等學者不勝枚舉。本文採用最經典的 Pregrine(1967)推導的 Boussinesq 模式求解。Peregrineg 使用水平向水深平均流速與瞬時水深為變數，頻散效應為二階精度，可包含底床深度、斜率與曲率變化之影響。因本文僅考慮單一斜面底床連接兩不同等水深之水域，因此頻散的效應，除了一般平面底床的頻散外，尚有底床坡度的影響(但無底床曲率效應)，這些都是垂直加速度引致水平動壓力的作用，將與水面坡度 (靜水壓作用力)及水體水平慣性動量(包含非穩態與非線性流傳) 抗衡過程演變的波動現象。孤立波入射斜坡前之波形，可由初始 KdV 正合解(

_sech

2_{孤立波)由上游遠端演進至 Boussinesq 的近似} 解。尾隨主波後方之震盪尾波因波幅甚小而逐漸落後，對該主孤立波與斜坡互制預期將無影響。主波到達斜坡受其作用，擾動安定的入射孤立波，使其垂直加速度效應變得明顯，因而加強非線性水平動量，使波前變陡而波峰後拉長，使波形空間對稱性破壞。Lamb(1932, §185)討論 Green(1837)定律：描述直角渠道緩變斜坡上線性長波波幅值與該處靜水深之

−

_{1/ 4}

冪次成比率，

(10)

Pregrine(1967)證實此種結果尚為合理，但亦討論 Ippen & Kulin (1954)實驗對不同斜坡可回歸獲得不同冪次值結果。此種淺化過程因同時考慮頻散作用而緩和趨向碎波之波形。若斜坡較長，則主波峰後方將產生另一較低波峰，以較低波速跟隨主峰前進，此種崩解機制將持續發生，而有可能發生系列次波峰。當波列到達較淺棚台時，局部頻散與非線性重組平衡機構，並使波峰高者有較快傳遞速度，而陸續脫離後方次波峰，形成分裂的孤立波(子) 列現象。孤立子經斜坡崩解成特定數量個孤立子在 Mei(1989, p563)一書亦有討論，並有引述其他學者的著作，在此不另作討論。孤立波經斜坡亦會產生反射波，根據 Pregrine(1967)以線性化 Airy 方程分析，反射波成一近乎平台狀水位台高，長約 2 倍入射波長，其高約 0/ 3 / 2 s A ，式中

s

為斜坡坡度。這些分析都限制斜坡的影響為高階微擾量，對實際應用各種海岸坡度較為不便，最終分析仍需回到 Boussinesq 方程組模式之數值分析。雖然 Boussinesq 方程組有許多形式，基本數值方法大同小異，本文將採用 Peregrine 推導的模式分析。數值求解 Boussinesq 方程組時，必須適當處理時間微分與非線性項，為了減少衍生數值頻散與耗散誤差，本文使用時間與空間均為二階精度的驢散方法：時空微分項均以中央差分離散(時間離散點為兩時階之半，而線性化項則用兩時階平均法進似。這個離散法之優點將於後文討論之。本文將探討單一孤立波傳遞至斜面斜坡連接到較淺陸棚的波形演變結果，雖然計算的例子非常簡單，但因包含非線性與頻散交替演進而變化波形，除上段討論之斜坡地形與非線性長波互制作用所產生的物理現象外，數值計算衍生的誤差波與上述交互作用尚可構成複雜的影響，本文試圖以各種量化觀點來分析討論這些演變影響，或許可提供實際數值應用時之參考。

二、數值模式

應用非擾動水深h 來正常化所有長度尺度，₀ h₀/g 來表示時 間尺度(壓力使用靜水壓ρgh0尺度已經被隱藏合併在Boussinesq 模式中，式中ρ 為水密度與 g 為重力加速度)。因此Boussinesq方 程組可表為下列無因次形式(Peregrine, 1967)： 1: t ( )x 0 L h + hu = (1) 2 2: 2 3 d d t x x xx t x xt xxt L u +uu +η = d u +dd u + u (2) 上兩式中座標下標表對其微分，其中 h≡ + 表(總)水深， d 表d η 靜水深， u 表 x 方向之深度平均速度。KdV方程式在單位靜水深 渠道之非週期性恆久波(permanent wave)正合解為單一孤立波，若其波峰位於x=x₀，峰高為A₀，其波形以水深表示時，可寫為： 2 0 0 1 sech [ ( )] h= +A k x− +ct x (3)

η

圖 1 本文探討之孤立波淺化示意圖式中波數k= 3A₀/ 2，波速c= d+A₀ ，均為波幅A 之函數，₀ 此為非線性水波的特性。本文將以此孤立波形為Boussinesq方程組模式之初始條件。須特別注意：此解僅為Boussinesq方程組在 定靜水深(d為常數)時之近似解；因此，即使在水平渠道，此波形 亦會演化到穩定的特定波形(且包含微小尾波現象)，此時的波形方為Boussinesq方程組離散式之數值解。如圖1，本文考慮演化後之孤立波由較深 (d= ) 水平底床經過 1:10 斜面坡坡到較淺1 (d=0.5)水平底床，探討孤立波因與斜面斜坡交互作用 (非線性淺化)而產生崩解成特定數量之孤立子列現象。數值計算式(1)(2)模式必須先將此連續的方程組對離散空間與時間格點進行離散化處理，而獲得代數方程組聯立系統。本文採用二階精度離散法，將所有變數在不同時階與格點位置，以泰 勒展開式近似之，如h以下式展開： 1 / 2 1 / 2 1 / 2 1 / 2 1 / 2 1 1₍ ₎2 2 2 2 1 ₃ 1 ₄ ₅ ( ) ( ) 6 2 24 2 ( ) n i n t n t n i t i tt i n i t _n t _n ttt i tttt i h h h h h h h O t + ∆ ∆ + + + ∆ ₊ ∆ ₊ ⎧ ⎫ ⎪ _{⎪ =} _± ₊ ⎨ ⎬ ⎪ ⎪ ⎩ ⎭ ± + + ∆ (4) 2 1 1 1 1 1 2 3 2 1 1 5 6 2 ( ) x n n n n i i x i xx i x _n x _n xxx i xxxx i h h x h h h h O x ∆ + + + + ± ∆ + ∆ + = ± ∆ + ± + ∆ (5) u 之展開亦同。因此模式之時空微分項皆以中央差分對x=i∆x， t n t=( +1/2)∆ 格點處離散，至於所有非線性向中各微分式之係數均以二時階平均值表之。如此將使離散精度達到二階 ) , ( 2 2 t x O∆ ∆ ，因此Boussinesq方程組離散為：連續方程： 0 8 ) )( ( 8 ) )( ( 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 = ∆ − − + + + ∆ − − + + + ∆ − − + − + + + + − + − + + + + + x h h h h u u x u u u u h h t h h n i n i n i n i n i n i n i n i n i n i n i n i n i n i (6) 動量方程： 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 2 ( )( ) 8 ( ) 4 2 ( 2 ) ( ) 2 ( ) ( ) [ ] 2 2 [ 3 n n n n n n n n i i i i i i i i n n n n i i i i i i n n i i i i i i n n n n i i i i i i i i u u u u u u u u t x h h h h d d x x d d d d u u x t d d u u u u d x x t d u + + + + + + − − + + + + − − + − + + − + + + − + − + − − ₊ + + − − ∆ ∆ + − − − + − ∆ ∆ − + − = ⋅ ⋅ ∆ ∆ − − − − + ⋅ ⋅ ∆ ∆ ∆ + 11 11 1 1 1 2 2 ( 2 ) ] 3 n n n n n n i ui ui ui ui ui x t + + + + + − − − + + − − ∆ ∆ (7) 模式(6)(7)高於二階以上之微分項已截切(truncate)。然此模式因此而被改變，將比原(1)(2)式多了數值頻散與耗散項(水平底床之例參見唐啟釗等，2000)。僅簡述計算流程如下：在某一時階，先從 式(7)疊代計算 u 至收斂(因含非線性項故需疊代)，再由此 u 值代 入式(6)計算 h (直接用三角元素矩陣消去法求得而無須疊代)，再 將此 h 解代至式(7)求 u 值，反覆上列步驟直到 u 與 h 兩值之改變 絕對值和(L 模數1 )小於設定誤差界限(取10−5以下)，再進行下時階計算。因時階每個時階都需使求解變數收斂到某精度範圍，而終能在該時階同時滿足二階離散聯立代數式，故可稱為二階精度之時間精確(time accurate)解。

三、數值結果與討論

本文考慮初始條件為峰高 0 A =0.5 之入射孤立波由左方 (

x

=

0

)進入，其波形為 KdV 方程解析解之 sech2函數，如式(3) 所示，僅為 Boussinesq 方程組近似解，因此為了允許演進到 Boussinesq 離散模式解，必須設定足夠遠之位置x 。一般數值經₀ 驗約需 4~5 倍波長，以A =0.5 之孤立波可定義波長為高度到₀ 0 A /10 之寬度(Lamb, 1932)，約為水深 5.6 倍，故以x0=10(初始 d=1 s=1/10 d=0 5

(11)

孤立波峰至斜坡趾部距離為 40 倍水深，斜坡趾部位於 x=50 處)，行進到斜坡趾部前約有 30 單位時間左右足供發展到 Boussinesq 解(參見圖二計算結果)。此後波形開始受斜坡淺化影響而與之交 互作用，並由主峰後端開始分裂出次波峰的情形。由於格網 x∆ 或 時階 t∆ 的大小均會直接影響求解變數的精度與計算結果的品質，如格點過大將產生過大的數值頻散與耗散。頻散為相位誤差，特指不同波長成分的波以不同波速傳遞；耗散為振幅誤差，指其傳遞隨時間逐漸變小的現象。本文作者(2000)曾對不同的時間與空間尺寸計算結果，探討數值頻散與耗散等影響，並用孤立波的特性來確認模式之可靠性。經過系列研究發現：使用上述的離散法可以適切的計算高振幅的孤立波。有趣的是，這種固定的格網在波傳過程，只有很偶然的情況使波峰點恰好落在格點位置，因此應該無法瞬時捕捉到波峰的位置，而容易形成忽高忽低的震盪形式。然而經過格網尺寸的評估，適當選擇格網，可以減少震盪的情形。本文因此使用空間尺寸值_{∆ =}_x _0.05，時間尺寸值 0.1 t ∆ = 。 30 40 50 60 70 -0.8 -0.4 0 0.4 0.8 h t=0 slope=1:10 t=10 t=69 t=40 t=35 condition: dt=0.1 dx=0.05 0 10 20 30 40 -0.8 -0.4 0 0.4 0.8 h 60 70 80 90 100 domain x -0.8 -0.4 0 0.4 0.8 h t=20 t=30 t=45 t=50 t=60 (a) (b) (c) 圖 2 初始孤立波經斜坡陸棚各區之波形演變圖 (a)斜坡前等深段；(b)斜坡段；(c)斜坡後等深段圖 2 為孤立波經過坡度 0.1 波形不同時間波形演變情形，由計算結果發現：當孤立波於斜坡前等水深段時，僅有些微誤差所產生之尾波(誤差波)追隨其後，而不影響本文探討斜坡互制作用與孤立波分裂現象的分析。由圖 2 亦可觀察到孤立波傳播至斜坡時，波形漸成不對稱性(_t_≈₃₅_s)，且主波峰後端開始分裂出較低的次波峰，注意此分裂點並非位於峰點而是後方近反曲點處，即二次空間微分項開始變號處，顯然的是頻散作用構成的波形演變。當波峰抵達斜坡頂端(_t_≈₄₀_s)時，因非線性使較高波高的主波峰以較快速度移動，致其脫離次波峰，波形逐漸分裂(fission) 成兩獨立孤立波，相同的，落後的次波峰又尾隨著更低的次波峰，最後於堤頂等深段再度分裂，形成三個孤立子各自行進的現象(_t_≈₆₉_s)。可以想像的是，每個孤立子成形後各自維持恆久波形，意味每個子系統之頻散與非線性均已達到局部平衡。自然界之神奇實在妙不可言。另外，將式(2)動量方程中等號右方之頻散項可分成斜坡曲率(α)、斜坡斜率(β)與局部水深(γ)項等三種主要頻散影響及非線性(λ)，分別定義為： t xxu d d 2 = α (8) xt xu dd = β (9) xxt u d 3 2 = γ

(10)

x uu λ= (11) 式(10)就是傳統水平底床(取d= )的頻散項。孤立波互制過程，1 斜坡曲率(α)、斜坡斜率(β)與局部水深(γ)三頻散項及非線性 (λ)效應，各以時距 t=5 倍數之時間演變分佈。因本文僅採用平面斜坡，只有斜坡之首尾端斜率不連續，因此以有限差分法計算 時，僅在此兩點有不為零之曲率值，而且此值將與 x∆ 有關。目前的斜坡條件計算所得之曲率的頻散項α 值約介於 [-4.e-10, 1.2e-9]，幾乎看不出其效應，是可以理解的。至於斜率的頻散項β 值只在於斜坡段不為零，並由計算結果發現：β值僅於接近斜坡頂端時產生最大影響量，約介於[-4.e-9, 1.2e-8]，亦是微不足道。而局部水深之頻散影響γ 值，在三者中影響最大，其值約為[-3.e-6, 3.e-6]。為進一步了解發展成熟(符合Boussinesq解)之單一孤立波接近斜坡時波形演變情形，我們選取t=30與40二時間來觀察各頻散項與非線性之空間變化，由於這些項彼此階次不同，因此將不同空間對應結果取絕對值並以對數表示，繪如圖3。整體而言，其曲率與斜率之頻散影響量α β/ =O( )ε ，約差一階次左右，而且皆可忽略，此與Mei(1982)論述相同。另外，局部水深之頻散與斜率之頻散影響量約差二階次左右，β γ_/ =_O₍ε2₎_。γ _{對頻散的影} 響不但在斜坡上很重要，且由較深水域延伸至較淺水域均都很重要。這個結論與Mei (1989，11.8小節)使用的推導假設(dx=O( ),ε 2 ( ) xx d =Oε )獲得_u=_{O d}₍ −1/ 2η₎_{的結論十分吻合。另外由本文計算} 結果可知：單一孤立波在較深水平底床(d= )尚未上斜坡時，即1 已經存在與非線性在低階相抗衡的頻散量γ|_d₌₁，其高階差量亦是造成_{sech 波形孤立波因須滿足Boussinesq模式而產生尾波的重}2 要機制。此頻散量γ與非線性量會持續變化於斜坡上，這兩量的空間分佈顯然不是局部平衡的。在斜坡上，因淺化使非線性λ值增加的影響而使孤立波之波前變陡，且波幅持續增加。相對地，主要頻散項γ值的影響，如同Mei之分析 ₍ 1/ 2 ₎ xxt xxt u =O d− η ，加速度的二次空間變率 ( )u_{t xx}與液面曲率η_xx時變率同階次，因而在較弱非線性的背波面反曲點上γ 值，次波峰因此產生。故當傳遞至較淺水平底床時，非線性與頻散λ|_d₌_0.5的再度重建平衡結構，促成斜坡後孤立波分裂的主因。 35 40 45 50 55 60 65 70 x domain -40 -20 0 -40 -20 0 -40 -20 0 -40 -20 0 圖 3 孤立波曲率、坡度、斜率等頻散與非線性項演變圖示(實線為t=30；虛線t=40) 上述結果僅討論波形在空間的演變情形。現若欲模擬「實驗時在某些特定位置量測波高資料的情況」，來研判分析時間演變對這種互制作用現象。這種分析觀點的討論將有助於對非線性波波形演變資料收集分析與相關現象了解。今假設波高計分別置放 ' α ' β ' γ ' λ 6 10− × 6 10− × 6 10− × 6 10− ×

(12)

於 x=15、40、50、55、65、75、85 等七種位置求取其 0~70 單位之時間函數之水位值，圖 4 為波高時序列圖。 0 10 20 30 40 50 60 70 time domain 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 -0.2 0 0.2 0.4 0.6 0.8 1 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 圖 4 孤立波淺化過程不同位置點水位演變圖 由上圖可知於當 x=15 時 KdV 孤立波已完整發展為 Boussinesq 孤立波，雖然尾隨其後尚有震盪尾波，但其量甚微，對主孤立波與斜坡(坡趾位於x=50)互制應無影響。且每一位置隨時間演進波形之數值誤差尾波雖都可察覺其存在，但與主波形相較都是微小量而可忽略。由這些時間變化可發覺峰值發生時間前後的對稱性較空間波形分佈(參考圖 2)的對稱性高。較明顯不對稱的位置在_x₌₅₅處，對應單位時間發生在 t=35 左右。圖 2 與圖 4 交叉比對可以提供孤立波非線性淺化互制現象很好的參考，對於了較淺陸棚使孤立波分裂的現象也很容易由此兩圖看到。有趣的是孤立波的分裂必須登上較淺陸棚才會發生，斜坡的作用似乎只在產生第二個次波峰，經過許多的數值實驗與不同斜度的斜坡，亦可得到相同的結論。顯然孤立波重建頻散與非線性的平衡機構必須有足夠長度與時間，而且必須在水平底床的條件下這個平衡機構才能被維持住。 0 40 80 Nf 0.001 0.01 0.1 1 10 100 A(f ) x=55 x=65 x=75 x=65, 75 x=15, 40 x=50 圖 5 孤立波在相對圖四不同位置之頻域波譜演變圖 (轉換資料 1024 筆，最高頻位於 Nf=256) 頻散現象主要由非線性頻散關係_{ω ω}₌ _{( )}_k 決定，且預期孤立波之分裂將使頻域資料改變，因此頻域的分析有其必要性。本文以 x=40、50、55、65、75、85 不同位置之波形時序列資料各選取 1024 筆，利用快速富立葉轉換(FFT)將其轉成頻域波幅譜，繪如圖 5(橫軸 Nf 為對應之頻率點編號，因 Nyquist 頻率故最高頻編號位於 max(Nf)=256；縱軸為富立葉係數，與真正波譜密度成比率)。總體而言，孤立波為長波現象，尤其是單一孤立波之頻率集中於低頻帶(常數項與波形到靜水位間之總體積成比率)，至高頻將呈指數衰減，為了明瞭較寬頻域範圍受地形淺化互制作用之影響，仍舊使用對數表示各譜分量之波幅量，由圖 5 可看到 x=15 至 x=40 為 d=1 較深之水平底床區，尾波的振幅僅改變高頻區的分佈，且都是較小量階的變化。X=50-55 為斜坡帶，這是淺化過程波形改變產生次波峰的效應，高頻處有較大幅度的變動。尤其在 x=55 處 Nf>30 以上看出明顯改變(事實上因目前使用半對數座標，高頻處的明顯變化僅變化振幅最大量約 10%左右而已)，這應是非線性淺化的作用，因為線性波與結構物的交互作用都是低階效應，而高階的交互作用必須由低頻的非線性組合或高次時空微分頻散項構成，來一起改變頻域的資訊。當 x=65 與 x=75 兩位置均已在較淺的陸棚上(d=0.5)，孤立波的分裂使低頻分佈產生明顯變化(與最大振幅同量階之變化)，已呈現在波譜圖上。 10 20 30 40 50 60 70 80 90 100 x domain 0.4 0.5 0.6 0.7 0.8 0.9 1 ev ol u ti on s of w ave p eak 0 10 20 30 40 50 60 70 transmission of A0=0.5 圖 6 不同時間下孤立波尖峰波高與傳遞位置演變圖最後在圖 6 繪出主孤立波峰隨時間傳遞到不同位置，其峰高演變圖。在斜坡區 X=50-55 間波峰急劇增加。然而必須注意恆久孤立波最大波高約為靜水深之 0.82 附近，本圖顯示為波峰高程，陸棚(較淺之水平底床)高程在 0.5 處，因此分裂的第一個孤立波峰發展到穩定時(t=50-70)，已超過恆久孤立波最大波高值。這是 Boussinesq 模式的缺點，顯示第二階的微擾模式有高估波峰的情形，預期更高階的頻散效應將可抑制波峰成長的極限，而能符合恆久孤立波最大波高值之理論。另外，由其他學者實驗亦證實在目前初始波高條件下，孤立波將會碎波。然而因計算方法的安定性，在如此強烈的非線性作用，頻散項卻能有效的抑制碎波的發生，模式本身與計算方法的太過安定，雖然不能反應適切的物理現象，但無疑地，佐證現行計算方法的「可用性」。更進一步選擇使用更好且通用的波場模式，顯然是改善數值模擬「實用性」很重要的步驟，而不僅是數值方法本身的安定與準確性而已。

四、結論

1 ch η 2 ch η 7 ch η 6 ch η 3 ch η 4 ch η 5 ch η x =15 x =40 x =50 x =55 x =75 x =65 x =85

(13)

本文以數值方法分析孤立波經斜坡(坡度 s=0.1)到較淺等深水域的淺化過程，因斜坡底床作用引致波形變化並使孤立波在陸棚上崩裂而分解成數個不等高度之孤立波列與震盪尾波現象。本文採用 Pregrine(1967) 發展之通用 Boussinesq 方程組長波模式，以二階時空精度之中央差分與兩時階平均方法求得時間精確解。本文以數值觀點，量化相關波傳動力性質，以觀察孤立波在較淺水域崩解分裂現象過程細部結構，了解非線性淺化與頻散的互制其實是一種微弱的作用。經過一些波形與頻域資料的分析，本文亦簡單地討論使用 Boussinesq模式的缺點，雖然目前的數值方法能夠安定有效地計算數值解，但亦建議選擇適切模式的重要性，或許可提供後續研究進一步探討非線性波淺化模擬的建議。

五、謝誌

本研究部分成果由國科會計畫NSC 92-2211- E006-030經費補助，特此致謝。

參考文獻

1.Airy, Sir G.B. (1845) “Tides and waves,” Encyc. Metrop., Art. 260. 2.Benjamin, T.B., Bona, J.L. & Mahony, J.J. (1972) “Model

equations for long waves in nonlinear dispersive systems,” Philos.

Trans. Roy. Soc. A, 272, pp47.

3.Boussinesq, J (1872). "Théorie des ondes et des Ramous qui se Propageant le Long d’un Canal Rectangulaire Horizontal, en Communiquant au Liquide Contenu dans ce Canal des Vitesses Sensiblement Pareilles de la Surface au Fond," J Math Pures Appl, (2) Vol 17, pp 55-108.

4.Green, G. (1837) “On the motion of waves in a variable canal of small depth and width,” Camb. Trans. vi, Papers, p225.

5.Korteweg, D.J. & De Vries, G. (1895). "On the Change of Form of Long Waves Advancing in a Rectangular Canal and on a New Type of Long Stationary Waves," Phil Mag, (5), Vol 39, pp 422-443. 6.Lamb, H (1932). "Hydrodynamics," Cambridge University Press. 7.Long, R.R. (1964) “The initial-value problem for long waves of

finite amplitude,” J. Fluid Mech., vol. 20, pp 161-170.

8.Madsen, O. S. and C. C. Mei (1969) “The trans- formation of a Solitary Wave over an uneven bottom,” J. Fluid Mech., Vol. 39, pp.781-791.3.

9.Mei, C.C. & B., Mehaute (1966) “Note on the equations of long waves over an uneven bottom,” J. Geophys. Res., vol. 71, pp509-517.

10.Mei, CC (1983). "The Applied Dynamics of Ocean Surface Waves," Wiley, N.Y.

11.Nwogu, O (1993). "Alternative Form of Boussinesq Equations for Nearshore Wave Propagation,” J. Waterway Harbor Coast Eng, ASCE, Vol. 119, No 6, pp 618-38.

12.Peregrine, D.H. (1967) “Long wave on a beach,” J. Fluid Mech. Vol. 27, pp815-827.

13.Tang, CJ., Hsu,CC, and Chang, DL (1999) “Error analyses on the solitary-wave solution of the Boussinesq equations,” The 21

conference on Ocean Engineering, pp175-182 (in Chinese).

14.Tang, CJ, Hsu, CS, and Lo, SY (2000) “Interaction of physical and error waves in numerical analysis on waves,” The 2000 symposium

on Cross-Strait Port and Coastal Development, pp74-81 (in Chinese).

15.Tseng, M.H. (1997) “Solitary wave solution to Boussinesq equations,” J. Wtr. Port, Coastal Ocean Engrg., ASCE, Vol. 123, No.3, pp138-141.

16.Yih, C.S. & Wu, T.Y. (1985) “General solution for interaction of solitary waves including head-on collisions,” Acta Mechanica

Sinica, Vol. 11, No.3, pp193-199.

17.Lin, C.C. & Clark, A. (1959) “TsingHwa, J. Chinese Studies, Spec. No. 1, Nat. Sci., 54.

18.Wei, G., Kirby, J.T., & Grilli, S.T. (1995) “A fully nonlinear Boussinesq model for surface waves, Part I,” J Fluid Mech., vol 294, pp 71-92.

19.Wu, T.Y. (1981) “Long waves in ocean and coastal waters,” J.

Engrg. Mech., ASCE, vol. 107, pp501-522.

20.Wu, TY (2001) "A unified theory for modeling water waves,"

Advances in Appl. Mech., 37, pp 1-88.

21.唐啟釗、許正昇、羅聖源(2000) ” 數值波浪分析中誤差波與物理波交互作用之現象，” 2000 兩岸港口及海岸開發研討會論文集，74-82 頁。

二維水深積分模式應用於三維亞臨界分流分析

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