1
行政院國家科學委員會專題研究計畫成果報告
以有限元素法模擬波浪之變形
Simulations of Wave Tr ansfor mations by Finite Element Method
計畫編號：NSC 89-2611-E-006-071
執行期限：89年8月1日至90年7月31日
主持人：許泰文 國立成功大學水利及海洋工程系(所)
一、中文摘要 本計畫擬利用有限元素方法，建立近 岸波場模式，模擬波浪以大角度入射近岸 地區時，波浪因地形、人工結構影響產生 之折射、繞射、反射、淺化及碎波等波浪 之變形。模式將以包含波浪繞射及反射之 緩坡方程式建立數值模式，利用有限元素 法易於處理不規則邊界的優點，引用二階 近似輻射邊界條件，擴展波場模式處理大 角度波浪入射的能力，探討在未碎波的情 形下，二階輻射邊界條件對於任意角度波 浪入射之適用情形。 關鍵詞：有限元素法，大角度波浪入射， 緩坡方程式，輻射邊界條件 Abstr actIn this study, a finite element method is used to develop a numerical model for the simulation of combined wave refraction, diffraction, reflection, shoaling and breaking in the nearshore zone. The governing equation is the mild-slope equation. Finite element method has advantage to treat complicated boundary. The second order radiation boundary condition is expanded to improve the wave model for large incident wave angles. The validity of the second order radiation condition is examined through varying arbitrary wave incident angles.
Keywor ds: finite element method, large-angle incident waves, mild-slope equation, radiation condition. 二、緣由與目的 台灣四面環海，當波浪由外海向近岸 傳播時，由於地形水深之變化與近岸結構 物之影響，波浪會發生折射、繞射、反 射、淺化及碎波等變形，為了描述此一系 列之物理現象，以數值方法來模擬分析近 岸波場為重要的方法之一。一般而言，利 用數值方法模擬近岸波場之變化時，在不 考慮碎波情況下，經常遇到的問題為不規 則邊界之處理與波浪大角度入射問題；目 前常用的數值方法為有限差分法，但受限 其離散化過程之影響，其矩形正交網格不 易與不規則邊界契合，且對於結構物附近 局部波場模擬之網格較為粗糙，而造成解 析度相對降低，此為有限差分法在使用上 之限制，因此，為了避免上述之問題，以 有限元素法來模擬近岸波場則為可行方法 之一。有限元素法的優點在於處理不規則 邊界時，可將其網格做適當之調整，以符 合不規則邊界之特性，因此，其網格處理 比有限差分法更具靈活性。 往昔學者利用有限元素法處理波浪斜 向入射之問題時，對於外海開放邊界之處 理，大多採用之輻射邊界條件是假設波浪 以垂直邊界的方向通過邊界。但事實上， 數值模式計算時，由於波浪受到地形水深 的改變或結構物之影響，使得波浪通過邊 界時，並不一定會以垂直方向通過邊界， 因而導致計算結果與實際情形之差距。為 了解決此一問題，必須對外海開放邊界輔 以 解 析 解 或 其 他 數 值 方 法 (Berkhoff ， 1972 、 Chen 和 Mei ， 1974 、 Bettess 和 Zienkiewicz，1977、陳和蔡，1990)，亦或 設置特殊之邊界，如圓形或半圓形邊界 (蘇，1993)，但此種處理方式必須滿足設 定邊界處為等水深，或劃定之邊界需符合 等水深波浪反射之條件，如此一來，邊界 條件的設定更趨於複雜。基於上述的種種 因素，引發了本研究之動機，希望能以一 般化之輻射邊界條件，處理波浪大角度入 射之問題，而不需將外海開放邊界再經過 特別的處理，如此將可簡化模式的複雜
2 性，並提升模式效率。 本研究計畫共分三年，本年度為第一 年計畫，主要發展具有二階輻射邊界條件 之波浪場模式，探討在未碎波情形下，二 階輻射邊界條件對於任意角度波浪入射情 況之適用情形。 三、理論基礎與數值方法 本文中模式計算所採用之座標系統如 圖1所示，包含卡氏座標系統 (x,y)，以及 計算邊界上之局部座標系統 (n,s)，其中n 為垂直邊界之法線方向，s為平行邊界之 切線方向，(n,s) 座標系統在外部計算領域 邊界遵循逆時針法則，內部結構物邊界則 為順時針方向。 x y n s n s n n s s 結 構物邊界 計 算領域邊界 結構物 邊界 計算領 域邊界 圖1 計算領域之 (x,y) 與 (n,s) 座標關係圖 模式採用橢圓型態之緩坡方程式作為 控制方程式，其表示如下： (CC ) CC k2 0 g g∇η + η= ⋅ ∇ (1) 式中，η=η(x,y)為波浪流速勢，C為波
速 (wave celerity)，C 為群波速度 (group g
velocity)，k為週波數；k=2π/L，ω為角 頻率；ω=2π/T， L 為波長，T為週期， 而∇為水平方向梯度；∇=(∂/∂x,∂/∂y)， x 表垂直岸線座標；y 為平行岸線座標。 模式計算採用之邊界條件，依邊界條 件特性可分為全反射邊界、部分反射邊 界、輻射邊界條件以及外海波浪入射邊界 條件等。其中全反射邊界、部分反射邊界 則可以由輻射邊界條件表示之。本文利用 二項式展開法將輻射邊界條件展開至二 階，如下式表示： ) 0 s k 2 1 ( k i n 2 2 2 ∂ = η ∂ + η α − ∂ ∂η (2) 式(2)中當α=0及α=1時，分別代表全反 射邊界和透過邊界，而0<α<1則代表部 分反射邊界。 至於在外海開放邊界處，邊界上包含 了入射波浪ηi及散射波浪ηs，波浪流速勢 可表為η=ηi +ηs。其中散射波ηs需滿足 輻射邊界條件。根據式(2)，可得到外海入 射邊界條件二階近似表示式如下： ∂ η ∂ α − η α − ∂ η ∂ = ∂ η ∂ 2 i 2 i i s k 2 i k i n n ∂ η ∂ α + η α + 22 s k 2 i k i (3) 上述外海入射邊界條件式(3)中，若令 i η 項為零，則可簡化為式(2)。因此，本文 中將以外海入射邊界條件式(3)作為通用邊 界條件，當選用的邊界條件為全反射邊 界、部分反射邊界及輻射邊界條件時，只 需令ηi為零即可。 本研究以 Galerkin 有限元素法作為數 值計算的方法，利用加權殘差法則以及散 度定理將控制方程式改寫為加權殘差式。 模式採用三角形元素建構數值格網，由於 邊界條件式式(2)與式(3)中均存在二階偏 微 分 項 ， 因 此 利 用 二 次 形 狀 函 數 (the quadratic shape function) 描述元素中各物 理量之分佈。 四、結果與討論 (一) 波浪通過等水深直立圓柱 選取波浪通過等水深直立圓柱的問題 進 行 模 式 驗 證 。 計 算 條 件 為 波 高 m 053 . 0 H0 = ，週期T=1.0sec，圓柱直徑 m 1 D= ， 水 深 h=0.35m ， 並 選 擇 0o 、 o 30 、45o等不同的波浪入射角度分別計算 一階及二階輻射邊界條件之波高分佈圖， 其計算結果如圖2、圖3所示。 由圖2(a)∼2(c)可以看出波浪正向入射 時，理論解 (MacCamy 和 Fuchs, 1954) 與 二階輻射邊界條件所計算出之波高分佈相 當接近，結果頗令人滿意；以一階幅射邊
3 界條件計算的結果則與理論解有較大的差 異。另外，當波浪入射角度45o時，圖3(b) 可看出利用二階輻射邊界條件所計算出之 結果，皆沿著入射波方向呈左右兩側對 稱；然而以一階輻射邊界條件計算所得之 結果（圖3(a)），波高分佈並不完全對稱 於入射波方向，波高分佈在尾水端邊界時 會垂直邊界的切線方向，此乃所設定的一 階輻射邊界條件是假設波浪以垂直邊界方 向通過邊界所致。二階輻射邊界條件多了 一個允許波浪斜向傳遞出去的修正項，因 此改善了此一不合理之現象。 圖4、5為沿波浪入射方向中央線剖面 之波高分佈圖。綜合比較結果可看出在圓 柱前端反射效應明顯之處，以一階、二階 輻射邊界條件計算所得之結果均接近理論 解；至於在圓柱後端遮蔽區處，由一階輻 射邊界條件所計算之結果有些許的振盪現 象，而二階輻射邊界條件則改善了很多， 尤其是入射角度為 o 45 時（圖7），一階輻 射邊界所產生之振盪現象更是明顯，這是 由於一階輻射邊界條件是讓波浪以垂直邊 界方向通過邊界，使計算領域內產生多餘 的反射波，所造成的誤差。 根據以上之討論，本文模式採用二階 輻射邊界條件模擬波浪大角度斜向入射 時，對於結構物突出水面的繞射效應之模 擬，有著非常不錯之結果。 (二) 波浪通過半球形淺灘地形 選取波浪通過一半球形淺灘地形問題 來驗證本模式模擬波浪之折、繞射現象。 計算條件採用 Ito和Tanimoto (1972) 之試 驗條件，地形水深分佈如下式所示： R r for h h R r for h h 0 m ≥ = < + = 2 0r e (6) 2 m 2 m 2 (x x ) (y y ) r = − + − (7) 2 m 0 0 (h h )/R e = − (8) 其 中 淺 灘 中 心 頂 部 水 深 hm =0.05m ， m 15 . 0 h0 = ， 淺灘半徑R=0.8m ，x 、m m y 代表淺灘之圓心位置，而計算之週期 511 . 0 T= sec， 入 射 波高H0 =0.0208m， 波長L0 =0.4m。 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 X (m) -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 Y (m) (a) 理論解 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 X (m) -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 Y (m) -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 X (m) -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 Y (m) (b) (一階解) (c) (二階解) 圖2 正向入射之波高分佈 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 X (m) -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 Y (m) -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 X (m) -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 Y (m) (b) (一階解) (c) (二階解) 圖3 斜向入射 ( 45° ) 之波高分佈 0.0 0.4 0.8 1.2 1.6 2.0 H/H0 θ = 0i o analytical solution 1 order solution 2 order solution -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0
Axis-dir ectional distance (m)
圖4 軸向波高分佈，正向入射 0.0 0.4 0.8 1.2 1.6 2.0 H/H0 θ = 45i o analytical solution 1 order solution 2 order solution -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5
Axis-dir ectional distance (m)
4 0.0 1.0 2.0 3.0 4.0 5.0 6.0 X (m ) 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Y (m ) 0.0 1.0 2.0 3.0 4.0 5.0 6.0 X (m ) 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Y (m ) 0.0 1.0 2.0 3.0 4.0 5.0 6.0 X (m) 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Y (m) 0.0 1.0 2.0 3.0 4.0 5.0 6.0 X (m) 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Y (m) 0.0 1.0 2.0 3.0 4.0 5.0 6.0 X (m) 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Y (m) 0.0 1.0 2.0 3.0 4.0 5.0 6.0 X (m) 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Y (m) 圖6∼圖8則為模式利用一階與二階邊 界條件在不同入射角度時之波高分佈圖， 其中左圖為一階幅射邊界計算結果，右圖 則為二階輻射邊界條件計算結果。波浪正 向入射時（圖6），可以看到一階輻射邊 界條件所計算出之波高分佈會有抖動的現 象產生，而二階輻射邊界條件則沒有；當 波 浪 以 45o 與80o斜 向 入 射 時 （ 圖 7 、 圖 8），一階輻射邊界條件所計算出之波高 分佈已有明顯不對稱情形，但是以二階輻 射邊界條件所計算出之結果（圖7、圖8之 右圖）卻相當良好，波高分佈幾乎無抖動 的現象。 圖6 以一階（左圖）、二階（右圖）輻射 邊界條件計算之波高分佈（正向入射） 圖7 以一階（左圖）、二階（右圖）輻射 邊界條件計算之波高分佈（入射角45度） 圖8 以一階（左圖）、二階（右圖）輻射 邊界條件計算之波高分佈（入射角80度） 五、計畫成果自評 本研究以有限元素法解析橢圓型態的 緩坡方程式，計算波浪的折、繞、反射等 效應，並以二階輻射邊界條件改善大角度 波浪斜向入射問題，經驗證比較與討論 後，歸納出以下幾點結論： 1. 有關波浪正向入射之問題，採用一階與 二階輻射邊界條件所計算出來之結果， 遠離開放邊界處其差異性並不大，但在 靠近開放邊界處，二階輻射邊界條件比 一階輻射邊界條件更具精確性。 2. 對於大角度波浪斜向入射問題，以二階 輻射邊界條件來處理，確實能大幅度改 善一階輻射邊界條件之缺點。 3. 模式對於波浪的折、繞、反射效應，以 二階輻射邊界來處理，均能正確地模擬 出結果。 六、參考文獻 [1] Berkhoff, J.C.W., 1972: “Computation of Combined Refraction-Diffraction”, Proceedings of the 13th International Conference on Coastal Engineering, 471-490.
[2] Bettess, P. and O.C. Zienkiewice, 1977: “Diffraction and Refraction of Surface Waves Using Finite and Infinite Element”, International Journal for Numerical Methods in Engineering, 11, 1271-1290.
[3] Chen, H.S. and C.C. Mei, 1974: “Oscillation and Wave Force on an Offshore Harbor”, Ralph M. Parsons Laboratory, Massachussetts Institute of Technology, Report No. 190.
[4] Ito, Y. and K. Tanimoto, 1972: “A Method of Numerical Analysis of Wave Propagation Application to Wave Diffraction and Refraction”, Proceedings of the 13th International Conference on Coastal Engineering, 503-522.
[5] MacCamy, R.C. and R.A. Fuchs, 1954: “Wave Force on Piles：a Diffraction Theory”, Institute of Engineering Research, Waves Investigation Laboratory, Serier 3, Issue 334, Berkeley California. [6] 陳伯旭, 蔡丁貴, 1990: “局部輻射邊界條件在水 波數值模式上之應用”, 第十二屆海洋工程研討 會論文集, 1-9. [7] 蘇青和, 1993: “多孔消波體波能消散解析及其 應用於港池之研究”, 國立成功大學水利及海洋 工程研究所博士論文.
