國 立 交 通 大 學
土木工程學系
博 士 論 文
非線性自由液面流場的數值模擬
以及其在探討微尺度上層海洋動力之應用
Numerical Simulation of a Fully-Nonlinear Free-Surface Flow
and its Applications on the Studies of
Microscale Dynamics in the Upper Oceans
研究生：洪立萍
指導教授：蔡武廷
非線性自由液面流場的數值模擬
以及其在探討微尺度上層海洋動力之應用
Numerical Simulation of a Fully-Nonlinear Free-Surface Flow and its
Application on the Studies of Microscale Dynamics in the Upper Oceans
研 究 生：洪立萍 Student：Li-ping Hung 指導教授：蔡武廷 Advisor：Wu-ting Tsai 國 立 交 通 大 學 土 木 工 程 學 系 博 士 論 文 A Thesis
Submitted to Department of Civil Engineering College of Engineering
National Chiao Tung University in Partial Fulfillment of the Requirements
for the Degree of Doctor of Philosophy
in
Civil Engineering December 2006
Hsinchu, Taiwan, Republic of China
非線性自由液面流場的數值模擬
以及其在探討微尺度上層海洋動力之應用
學生：洪立萍 指導教授：蔡武廷
國立交通大學
土木工程學系博士班
摘 要
本論文建立一完整描述三維非線性自由液面流場運動的數值模式。為 於未知的自由液面邊界上滿足非線性邊界條件，本模式利用座標轉換，將 非穩態且不規則物理區間的流場控制方程式與邊界條件，映射至規則的計 算區間，以進行數值解析。模式的主要數值方法包括：分別以擬頻譜法以 及二階有限差分法近似控制方程與邊界條件中的水平與垂直導數；以二階 Runge-Kutta 法進行非穩態的自由液面位移與流場速度的時間積分；利用改 良式割線疊代法求解計算區間中不可分離的壓力Poisson 方程式，以滿足不 可壓縮流條件。 本論文首先應用所發展的數值模式，探討二維重力–寄生表面張力波之 發展與演化。模擬結果顯示，當寄生表面張力波完全發展時，強烈的渦漩 自各個表面張力波之波谷剝離，形成高旋度的邊界層，進而大幅提升流場 能量的黏滯消散。為進一步顯示本模式解析非線性波浪運動之能力，我們 進行三維「短峰型重力波」與「新月型重力波」的計算；其結果完整模擬 波浪非線性交互作用之過程，所形成的波型結構與水槽實驗之觀測相同。 最後，我們應用本數值模式，模擬正向風壓與風剪應力作用下之三維自由 液面邊界層紊流場。模擬結果顯示：寄生表面張力波列自重力波波峯之前高低流速交錯的條痕結構。此流場結構與實驗觀測相同。以上三例有關三 維自由液面流場之應用，驗證本論文所發展之模式，可以正確解析自由液 面流場不同尺度的物理過程；包括液面重力波、表面張力波、黏滯邊界層 流場以及紊流結構。
Numerical Simulation of a Fully-Nonlinear Free-Surface Flow
and its Applications on the Studies of Microscale Dynamics
in the Upper Oceans
Student：Li-ping Hung Advisor：Professor Wu-ting Tsai
Department of Civil Engineering National Chiao Tung University
Abstract
A numerical model for simulating a three-dimensional flow bounded by a fully-nonlinear, free-surface boundary is developed. The free-surface boundary conditions are satisfied on the free-moving surface exactly without any approximation or linearization. The governing equations and the boundary conditions on the time-dependent, physical domain are transformed onto a regular domain for numerical computations. The major implementations of the developed model include: Pseudo-spectral and finite-difference schemes are employed to approximate the spatial-differential operators in the horizontal and vertical directions, respectively. A second-order Runge-Kutta method is adopted to integrate in time for the velocity field and the free-surface elevation. The solenoidal condition of the velocity field is satisfied by solving the transformed pressure Poisson equation. A modified secant method is used to accelerate the iterative solution of the non-separable Poisson equation.
We first apply the numerical model to study the development and evolution of parasitic capillary ripples on a two-dimensional gravity-capillary wave. The simulation results reveal that strong vortices are shed from the troughs of the capillary ripples, convected backward, and form a strong boundary layer. As a result, energy dissipation is enhanced with the generation of parasitic capillaries. The model capabilities for resolving nonlinear interactions among
the surface waves are then demonstrated by simulating the three-dimensional evolutions of the short-crested waves and the crescent waves. Finally, the developed model is applied to simulate a wind-driven turbulent boundary layer bounded by a dynamic water surface. The simulation results reveal fine surface structures, including parasitic capillary waves developed from the front face near the crest of the dominant gravity wave and streamwise velocity streaks on the back gravity-wave surface. These surface structures have been widely observed in the laboratory and field experiments. The three simulation examples also indicate the capabilities of the present model in resolving flow processes of various length scales, including gravity and capillary waves, viscous sublayer and coherent eddies of a turbulent boundary layer.
致
謝
本論文得以完成首要感謝指導老師蔡武廷教授六年半來的悉心指導。 受教於蔡老師的求學過程中，不僅提升學術成果品質更增廣對於科學研究 的視野，在此對蔡老師致上最深的謝意。感謝黃維信教授、傅武雄教授、 黃清勇教授、賴明治教授與葉克家教授於口試期間的指正與建議，您的指 教使本論文更加嚴謹並啟發我省思未來之研究方向。感謝於交通大學時期 與中央大學時期的實驗室夥伴，有了你們的陪伴讓我渡過美好的求學生 涯。感謝我的家人，你們給予我無限的包容與支持，使我能致力於完成學 業。目
次
中文摘要 ……… 英文摘要 ……… 致謝 ……… 目次……… 表目次……… 圖目次……… 符號表……… i iii v vi x xi xv 第一章 前言……… 1.1 研究背景……… 1.2 三維自由液面紊流邊界層數值模式的發展回顧……… 1.3 模式發展與應用……… 1 1 2 5 第二章 模式之數學理論 ……… 7 2.1 控制方程式與邊界條件……… 2.1.1 流場之控制方程式 ……… 2.1.2 自由液面的運動邊界條件 ……… 2.1.3 作用於自由液面上的表面張力 ……… 2.1.4 作用於自由液面上的風應力 ……… 2.1.5 自由液面上的應力平衡 ……… 2.1.6 自由滑移之流場底部邊界 ……… 2.2 座標轉換……… 7 8 8 9 10 11 12 122.2.1 座標轉換關係式……… 2.2.2 計算區間之流場控制方程式與邊界條件 ………… 2.3 能量守恆……… 12 14 16 第三章 模式之數值方法……… 17 3.1 網格點配置與空間導數的離散……… 3.2 Runge-Kutta 時間積分與壓力 Poisson 方程式………… 3.3 改良割線疊代法求解 Poisson 方程式……… 3.3.1 疊代形式之壓力 Poisson 方程式……… 3.3.2 直接疊代法與割線疊代法 ……… 3.3.3 疊代法的收斂測試 ……… 3.4 自由液面邊界條件之離散……… 3.4.1 自由液面之垂直速度 (wζ=1) ……… 3.4.2 自由液面速度在垂直方向上之一次與二次導數… 3.5 數值方法之驗證……… 17 21 27 27 30 31 32 33 34 36 第四章 寄生表面張力波……… 41 4.1 波形與流場結構……… 4.2 生成機制之探討……… 4.2.1 表面張力的影響……… 4.2.2 風壓的作用……… 4.2.3 初始波形斜率的影響……… 4.3 旋度場之生成機制 ……… 4.3.1 自由液面旋度場 ……… 41 46 46 51 55 73 73
4.3.2 通過自由液面之旋度通量 ……… 4.3.3 渦旋之傳輸 ……… 4.4 能量消散……… 4.4.1 寄生表面張力波引起之能量變化 ……… 4.4.2 平均液面斜率與能量消散的關係 ……… 76 78 80 80 80 第五章 三維非線性自由液面流場之模擬……… 83 5.1 短峰型重力波 ……… 5.1.1 模擬之初始條件 ……… 5.1.2 模擬結果 ……… 5.2 新月型重力波 ……… 5.2.1 模擬之初始條件 ……… 5.2.2 模擬結果……… 5.3 風波下之剪紊流場 ……… 5.3.1 模擬之初始條件……… 5.3.2 模擬結果……… 83 83 84 87 87 87 91 91 92 第六章 結論與未來展望 ……… 99 參考文獻 ……… 101 附錄一 流場能量守恆方程……… 附錄二 擬頻譜法之基礎理論－離散傅立葉級數……… 附錄三 低儲存量之二階 Runge-Kutta 數值積分法……… 附錄四 非線性重力波近似解之三階 Stokes 波……… 附錄五 全域座標系統與局部座標系統之轉換……… 108 113 125 127 129
附錄六 流場邊界之旋度通量……… 附錄七 短峰型重力波之三階近似解……… 附錄八 新月型重力波之模擬初始流場……… 131 134 137
表目次
表4.1：不同波長與初始波形斜率之重力波，其所產生之寄生
表面張力波波形特性。………
56
圖目次
圖2.1：物理區間與計算區間之網格分布圖。……… 13 圖3.1：計算區間垂直方向非等間距網格分布情形。………… 18 圖3.2：垂直方向非等間距之正規網格於物理區間與計算區間 分布圖。……… 18 圖3.3：正規網格系統與交錯網格系統配置圖。……… 19 圖3.4：於垂直方向使用正規網格系統以及共同正規網格與交 錯網格系統之配置示意圖。……… 20 圖3.5：以疊代法求解 Poisson 壓力方程式的測試結果。…… 32 圖3.6：模擬過程之平均水面高度、流場之平均速度散度以及 最大速度散度。……… 37 圖3.7：模擬過程中總能量之演化以及數值能量損失與黏滯消 散能量之比。……… 39 圖3.8：前進波之波形模擬結果。……… 40 圖4.1：波長λ =5cm，初始波形陡度ak=0.25之前進波在考慮 表面張力影響但不受任何風壓之作用下，其波形( )
η 、 波形斜率( )
ηx 與波形二次導數( )
ηxx 的演化過程。……… 42 圖4.2：波長λ =5cm，初始波形陡度ak=0.25之前進波在考慮 表面張力影響但不受任何風壓之作用下，於時間 0 0 0 0, 0.5 , , 2 25 . 0 T T T T t= 之瞬時速度場與旋度場。……… 44 圖4.3：寄生表面張力波之速度場向量圖與旋度分布。……… 45 圖4.4：波長λ =5cm，初始波形陡度ak=0.25之前進波在不考 慮表面張力影響且不受任何風壓之作用下，其波形( )
η 、波形斜率( )
ηx 與波形二次導數( )
ηxx 的演化過程。 47 圖4.5：波長λ =5cm，初始波形陡度ak=0.25之前進波在不考 慮表面張力影響且不受任何風壓之作用下，於時間 2 , , 5 . 0 , 25 . 0 T T T T t= 之瞬時速度場與旋度場。……… 48圖4.6：表面張力作用下之正向應力分布情形。……… 50 圖4.7：波長λ =5cm，初始波形陡度ak=0.25之前進波在考慮 表面張力影響並受正向風壓之作用下，其波形
( )
η 、波 形斜率( )
ηx 與波形二次導數( )
ηxx 的演化過程。……… 53 圖4.8：波長λ =5cm，初始波形陡度ak=0.25之前進波在考慮 表面張力影響且受正向風壓之作用下，於時間 0 0 0 0, 0.5 , , 2 25 . 0 T T T T t= 之瞬時速度場與旋度場。……… 54 圖4.9：波長λ =5cm，初始波形陡度ak=0.12之前進波在考慮 表面張力影響但不受任何風壓之作用下，其波形( )
η 與 波形斜率( )
ηx 的演化過程。……… 59 圖4.10：波長λ =5cm，初始波形陡度ak =0.16之前進波在考慮 表面張力影響但不受任何風壓之作用下，其波形( )
η 與 波形斜率( )
ηx 的演化過程。……… 60 圖4.11：波長λ =5cm，初始波形陡度ak =0.20之前進波在考慮 表面張力影響但不受任何風壓之作用下，其波形( )
η 與 波形斜率( )
ηx 的演化過程。……… 61 圖4.12：波長λ=5cm，初始波形陡度ak =0.25之前進波在考慮 表面張力影響但不受任何風壓之作用下，其波形( )
η 與 波形斜率( )
ηx 的演化過程。……… 62 圖4.13：波長λ =5cm，初始波形陡度ak =0.28之前進波在考慮 表面張力影響但不受任何風壓之作用下，其波形( )
η 與 波形斜率( )
ηx 的演化過程。……… 63 圖4.14：波長λ =7.5cm，初始波形陡度ak =0.20之前進波在考 慮表面張力影響但不受任何風壓之作用下，其波形( )
η 與波形斜率( )
ηx 的演化過程。……… 64 圖4.15：波長λ =7.5cm，初始波形陡度ak=0.23之前進波在考 慮表面張力影響但不受任何風壓之作用下，其波形( )
η 與波形斜率( )
ηx 的演化過程。……… 65圖4.16：波長λ =7.5cm，初始波形陡度ak=0.25之前進波在考 慮表面張力影響但不受任何風壓之作用下，其波形
( )
η 與波形斜率( )
ηx 的演化過程。……… 66 圖4.17：波長λ =7.5cm，初始波形陡度ak =0.28之前進波在考 慮表面張力影響但不受任何風壓之作用下，其波形( )
η 與波形斜率( )
ηx 的演化過程。……… 67 圖4.18：波長λ =7.5cm，初始波形陡度ak =0.30之前進波在考 慮表面張力影響但不受任何風壓之作用下，其波形( )
η 與波形斜率( )
ηx 的演化過程。……… 68 圖4.19：波長λ =10cm，初始波形陡度ak=0.25之前進波在考 慮表面張力影響但不受任何風壓之作用下，其波形( )
η 與波形斜率( )
ηx 的演化過程。……… 69 圖4.20：波長λ =10cm，初始波形陡度ak=0.28之前進波在考 慮表面張力影響但不受任何風壓之作用下，其波形( )
η 與波形斜率( )
ηx 的演化過程。……… 70 圖4.21：波長λ =10cm，初始波形陡度ak=0.30之前進波在考 慮表面張力影響但不受任何風壓之作用下，其波形( )
η 與波形斜率( )
ηx 的演化過程。……… 71 圖4.22：波長λ =10cm，初始波形陡度ak=0.31之前進波在考 慮表面張力影響但不受任何風壓之作用下，其波形( )
η 與波形斜率( )
ηx 的演化過程。……… 72 圖4.23：自由液面旋度分布圖。 ……… 74 圖4.24：波長λ =5cm，初始波形陡度ak =0.25之前進波，其自 由液面旋度通量之各項分量延著邊界之分布圖。…… 77 圖4.25：自由液面之瞬時旋度變化率。……… 79 圖4.26：波長λ =5cm，初始波形陡度ak =0.25之重力波總能量 (Ew)隨時間之演化。……… 81 圖4.27：平均液面斜率與能量消散率的瞬時關係。………… 82 圖5.1：短峰型重力波之三維菱形結構波形。……… 85圖5.2：在考慮表面張力作用下，短峰型重力波之三維菱形結 構波形。……… 86 圖5.3：二維前進波列受三維不安定擾動而形成三維新月型波 之演化過程。……… 89 圖5.4 : 在考慮表面張力作用下，二維前進波列受三維不安定 擾動而形成三維新月型波之演化過程。……… 90 圖5.5：三維風波之液面高度(η)、沿流向之液面斜率(ηx)與橫 向之液面斜率(ηy)。……… 94 圖5.6：三維風波之液面紊流速度場。沿流紊流速度(u′)、橫 向紊流速度(v′)與垂直紊流速度(w′)。……… 95 圖5.7：由波背相與波前相之風波三維液面變化檢視圖。……… 96 圖5.8：三維風波之橫向紊流旋度場(ω′y) 在位於y=0、0.24λy、 y λ 47 . 0 與0.82λy垂直沿流方向之流場剖面。……… 98
符號表
) , , (x y z 流場位置座標系統 ) , , (ξ ψ ζ 計算區間之位置座標系統 ) , , (u v w 流場速度 ) , , (n t1 t2 自由液面上之局部座標 } ˆ {n 自由液面之單位法線方向向量 } ˆ {t1 投影至 x-z 平面之自由液面單位切線方向向量 } ˆ {t2 投影至 y-z 平面之自由液面單位切線方向向量 η 自由液面高度 h 流場水深 (平均水平面與流場底部之垂直距離) p 流場壓力 (動力壓力 Hydrodynamic pressure)P 總壓力 ( Total pressure 或稱為 Thermodynamic pressure , P= p−ρgz)
ρ 流體密度 g 重力加速度 ν 流場運動黏滯係數 µ 流場動力黏滯係數 (µ=ρν ) σ 表面張力係數 κ 自由液面之曲率 s n τ 自由液面所受之風應力於液面法線方向{nˆ}之分量 s t1 τ 自由液面所受之風剪應力於液面切線方向{tˆ1}之分量 s t2 τ 自由液面所受之風剪應力於液面切線方向{tˆ2}之分量 0 T 波形之線性理論週期 Θ 風壓分布與液面波形之相位差
0 p 風壓參數 0 ω 自由液面旋度值，其指向頁內為 s u 自由液面切線方向速度分量 r u 自由液面法線方向速度分量 sˆ 與自由液面邊界相切之切線方向 rˆ 正向於自由液面之法線方向，其指向流場外為正 ϕ 前進波入射方向與風吹拂方向間的夾角 φ 自由液面之切線方向與全域( global )座標x方向之夾角 Re Reynolds 數 Fr Froude 數 We Weber 數
第一章
前言
1.1 研究背景
在非常接近海-氣交界面的紊流邊界層內，其流場動力特性為通過界面 間之動量、氣體、熱量與水氣等傳輸過程的主要機制。了解這些由接近液 面紊流結構所主控的交換過程為海-氣交互作用研究之主要課題。然而直到 近期發展之實驗技術才能合理地直接量測接近邊界處之紊流數值，因此以 往的相關研究皆以環境變數參數化界面之傳輸特性；例如：建立風速與通 過界面之傳輸率的關係（Liss and Merlivat，1986 與 Wanninkhof，1992）， 或將全球風場模式化以推算全球海-氣間二氧化碳通量（Takahashi et al.， 1997）。此外，由已知的流場物理背景建立流場結構之概念圖像以推測傳 輸過程的機制；稱之為概念式的理論，亦普遍地使用於海-氣交互作用研究。 其中包括：Coantic (1986) 之小渦旋模式（small eddy model）；以局部的紊 流消散係數參數化邊界層之紊流通量。Higbie (1935) 與 Danckwerts (1951) 之表面更新模式（surface renewal model）；以更新渦旋描述邊界層中之水 體交換。然而，風吹拂過海面所引起的傳輸過程是相當複雜的，除了風引起的 剪紊流外，海面上亦會出現短風波，而風波於液面上傳遞不僅影響邊界層 內紊流結構，其兩者間的交互作用亦可能改變液面狀態因而對界面傳輸有
相當大的衝擊。為深入地了解液面運動，Okuda (1982)、Ebuchi et al.(1987)、
Banner and Peirson(1998)、Caulliez and Collard (1999) 以及 Veron and Melville (2001) 等在實驗水槽中，對於由風場趨動之海面運動進行模擬與量 測，其對於液面演化過程與流場結構特性提供定性與定量的研究成果。但 有關海-氣傳輸的實驗量測與模式皆僅能建立部份重要環境參數；如風引起 之液面阻力速度、紊流擴散參數與液面粗糙度等和傳輸速度的相關性，並 無法透視邊界層內紊流結構或說明控制海氣-界面間之動量、質量與能量傳 輸的物理機制。
隨著近期遙測技術的發展與普及，液面紅外線影像的實驗與分析技術 被應用於可視化流場以提供控制海-氣交換機制的定量資訊（Jessup et al.， 1997、Jahne and Haussecker，1998、Veron and Melville，2001、Zappa et al.， 2001、Asher et al.，2004、Garbe et al.，2004 與 Schimpf et al.，2004）。 以紅外線影像量測流場的實驗方法之兩大優點為：其能同時取得所有傳輸 過程與接近液面紊流間之定量關係，以及所有傳輸速度的量測對於時間與 空間皆可提供足夠的解析度而可對外在環境參數於時間與空間之變異性進 行比對。Schimpf et al. (2004) 即指出液面溫度擾動分布與溫度次層的溫差 以及紊流邊界層中之表面更新結構有直接的關係，而其結論亦可由概念式 理論之推論值得到應證。 無論是應用液面紅外線影像實驗技術或是概念式理論模式，海洋表面 上的流場結構特徵為液面運動與紊流邊界層間交互作用的重要指標。而以 偵測海面流場結構特徵來推斷其下流場之動力過程的研究方法，其可信度 與準確性必須仰賴對於自由液面運動與紊流邊界層流場間相互作用的物理 特性以及相關的定量數據等具備相當深入且精確的知識背景。然而，流場 內複雜的非線性作用與傳輸過程之物理機制尚有許多無法由量測數據得以 解釋。欲解析流場內的細部過程，數值模擬即是最佳的研究工具，其可提 供流場演化過程之完整數據並可詳盡地檢視其動力特性。
1.2 三維自由液面紊流邊界層數值模式的發展回顧
早期相關於自由液面紊流邊界層的數值模擬主要包括：探討初始均質 等向性紊流與自由液面交互作用；Perot and Moin (1995) 與 Walker et al. (1996)、由底部無滑移邊界演生之紊流與自由液面交互作用之明渠紊流模 式；Leighton et al. (1991)、Lam and Banerjee (1992)、Handler et al. (1993)、 Komori et al. (1993)、Borue et al. (1995) 與 Pan and Banerjee (1995)，以及考慮平均剪力流中紊流與自由液面交互作用；Tsai (1998)、Shen et al. (1999) 與
Banerjee (1992)、Handler et al. (1993)、Perot and Moin (1995)、Pan and Banerjee (1995) 與 Walker et al. (1996)皆理想化自由液面邊界，將其視為一 自由滑移的水平邊界。而Borue et al. (1995)、Tsai (1998)、Shen et al. (1999)
以及Tsai et al. (2005) 等則線性化自由液面條件，將自由液面邊界條件滿足 於平均水面上而非真實的運動液面。Sullivan et al. (2004) 亦視平均液面高 度為流場邊界，將自由液面運動參數化為一流場邊界上的外力。由於對自 由液面運動的簡化，使得模擬結果只能反應自由液面的存在對於紊流運動 的影響，而關於自由液面運動對紊流場的作用以及紊流對於液面邊界變動 的貢獻等非線性物理過程皆無法描述。然而建立能真實地呈現自由液面瞬 時運動的數值模式之困難處在於：模擬的流場區間為不規則幾何區域；其 將增加計算過程中繁複的代數運算，以及流場區間之時變性；此則需要付 出大量的計算成本，並須考慮模式滿足不可壓縮流條件的穩定性。 隨著計算能力的進步，欲解析液面波動對於紊流結構的影響之數值模 式持續發展：Maass and Schumann (1994)、 De Angelis et al. (1997)、 Cherukat et al. (1998) 與 Calhoun and Street (2001)模擬無滑移靜止波形固 體邊界與自由滑移水平固體邊界間之紊流場，以了解邊界形變對於紊流場
的的影響。而Sullivan et al. (2000) 與 Shen et al. (2003) 進一步考慮等速移
動的正旋函數波形固體邊界以加入運動邊界對流場的影響機制。這些波形 邊界流場模擬對於計算過程中區域轉換運算技術已相當成熟，模擬結果亦 指出於接近邊界處的局部流場特性與水平邊界並不相同，但對於紊流場貢
獻於液面邊界變動之物理過程仍無法描述。Komori et al. (1994) 首次求解非
穩態自由液面運動邊界的模擬，但只限於微小的液面變動且其結果與線性
化自由液面條件之模擬並無太大不同。Hodges and Street (1999) 為近期成功
求解完全非線性自由液面邊界之數值模擬，但受限於龐大計算量，此模式
僅能使用網格數32×32×64，並不足以解析紊流邊界層。Fulgosi et al. (2003)
可由自由液面運動方程式追蹤液面變化，但其模式之液面最大波形斜率ak
受限於ak ≤0.01。
有兩種：固定網格（fixed grid）與移動網格（moving grid）。於固定網格模 式中，控制方程離散於物理區間中之固定網格點，並於網格中描述自由液 面邊界。於移動網格模式中，物理區間的邊界跟隨著自由液面運動，因此 必須根據每個離散時間點的液面高度變化更新對應邊界之物理區間計算網 格。固定網格模式有利於兩相流之自由液面；如氣泡或水滴等，或是可能 產生巨大形變或複雜的幾何形狀之自由液面；如波浪破碎等。其主要用以 捕捉自由液面高度變化的方法有：（1）marker-and-cell 方法；追蹤流場內 標記之質點位置（Harlow and Welch，1965；Raad and Bidoae，2005 等）。 （2）the level-set 方法；決定網格內自由液面以下的距離（Osher and Sethian，1988；Sussman，1994；Sussman et al.，2006 等）。（3）volume of fluid 方法；以流體體積守恆決定網格點內流場體積的比例（Hirt and Nichols，1981；Rhee et al.，2005 等）。上述固定網格模式所解析之自由液 面，由於其並無滿足質量守恆與動量守恆，且不具流體運動或動力特性， 因此無法應用於探討自由液面邊界層或是紊流場與自由液面交互作用之物 理過程。而移動網格模式，由於其在每個離散時間點都需重新製作網格， 因此需要耗費較多的計算成本，但自由液面變化完全決定於其運動與動力
邊界條件，因此較能合理地解析其邊界層流場。Komori et al. (1994)、Sullivan
et al. (2000)、Shen et al. (2003)、Hodges and Street (1999) 與 Fulgosi et al.
(2003) 等探討自由液面紊流邊界層的數值模擬皆使用移動網格模式描述自 由液面運動。
在移動網格模式之計算網格配置方法有兩大類；分別是使用 restricted
curvilinear coordinate grid 系統（簡稱為 RCCG）以及 generalized curvilinear coordinate grid 系統（簡稱為 GCCG）。此兩種網格系統皆可將不規則的物 理區間轉換至規則且各網格線間相互線性正交的計算區間中。RCCG 系統 主要是將物理區間之網格點在計算區間垂直方向上重新分佈，它的優點是 格點的產生較容易，轉換後的形式也較簡單，然因自由液面上之格點僅限 於垂直之移動，因此在大波幅之自由液面運動時，格點的扭曲亦變大。如 Dommermuth (1993)、De Angelis et al. (1997)、Sullivan et al. (2000)、 Lin and
Li (2002)、Shen et al. (2003) 與 Fulgosi et al. (2003) 皆使用 RCCG 來處理變 動的邊界。GCCG 系統中物理區間之網格點則是皆正交於邊界，所以格點
的扭曲較小，但格點的生成與轉換後的形式皆較複雜。如Komori et al. (1993)
與Hodges and Street (1999) 使用 GCCG 來處理變動的邊界。
為滿足不可壓縮流條件，需求解在計算區間中經座標轉換的壓力 Poisson 方程式，此方程式為不可分離之橢圓形式偏微分方程式，其求解方
法包括：直接求解全係數矩陣方程式（如Komoriet al.，1993、Lin and Li，
2002）、使用疊代法（如 Dommermuth，1993、Sullivan et al.，2000），或
應用multi-grid 法（如 Hodges and Street，1999）等。但在實際進行數值模
擬時，以上方法在計算時間與適用性都有某些限制。如若使用直接求解全 矩陣的方式，則空間離散格點數將受限制；若使用疊代法，則只適用於微 小變形的自由液面邊界。為了以疊代法有效求解大波幅自由液面運動的壓
力 Poisson 方程式，Zhang (1996) 以改良牛頓疊代法（modified Newton’s
iteration method）求解轉換壓力 Poisson 方程式。其不僅能加速原直接疊代 法的收斂速度，亦提高了可求解之邊界形變的最大波幅。但此改良牛頓疊 代法極易因疊代條件而導致疊代無法收斂，對此本模式將利用交替使用割 線疊代法與直接疊代法，以確定求解過程之收斂並能達到有效的加速效果。
1.3 模式發展與應用
本論文建立一能完整描述自由液面非線性運動之三維紊流模式，改進 Tsai (1998) 非靜水壓-不可壓縮流場模式中所運用之數值方法，將流場運動 之控制方程與非線性邊界條件滿足於一伴隨著自由液面運動的非穩態流場 區間。模式所解析的流場區間包含垂直方向之完全非線性自由液面邊界與 底部自由滑移邊界，以及水平方向之週期邊界。在數值計算的過程，ζ -座標轉換（或稱σ -座標轉換；Lin and Li，2002）應用於將非穩態不規則流場
物理區間中之流場控制方程與邊界條件轉換至穩態之規則計算區間中。在 經過轉換的計算區間中，利用改良割線疊代法求解不可分離的流場壓力
Poisson 方程式，其精確的壓力值將使流場得以滿足不可壓縮流條件。 Runge-Kutta 時間積分與自由液面非線性運動邊界條件將應用於追蹤非穩態 的自由液面變動。模式中水平方向導數以擬頻譜法近似離散而垂直方向導 數以二階有限差分法近似離散。垂直方向之非均勻網格分布提升接近液面 處的網格密度以增加紊流邊界層空間解析度。本模式的目的為解析自由液 面運動以及其與紊流邊界層間之交互作用，因此，本模式必須解析空間中 不同尺度的物理過程；包括液面重力波、表面張力波、黏滯邊界層流場以 及紊流結構。此結合擬頻譜法於週期區間之水平方向導數以及非均勻網格 分布之有限差分法於垂直方向導數的數值近似法具有最佳之空間解析能 力，且應用於離散三維的壓力Poisson 方程式可將其轉換為頻譜空間中的一 維全微分方程以提升模式的計算能力。此求解紊流邊界層流場的數值離散 方法廣範地使用於明渠紊流場（Moin and Kim，1982）與大氣紊流邊界層 （Moeng，1984）等。 非線性自由液面模式首先將應用於模擬一海面上常見之無包覆空氣微 小破碎波（micro-breaking）流場。在此二維模擬中完整地描述寄生表面張 力波之波形演化與流場結構特性，對於流場中強烈渦旋的生成機制與傳輸 亦有詳盡地探討。模擬結果證實寄生表面張力波的生成將引起流場中強烈 的能量黏滯消散，同時建立流場平均表面粗糙度與能量消散的線性關係。 對於三維自由液面運動的模擬；包括短峰型重力波（short-crested wave）與 新月型重力波（crescent wave），其液面波形特徵的演化將驗證本模式解析 波浪非線性交互作用過程的能力。而在風波與自由液面紊流邊界層的三維 數值模擬中，解析在風應力作用下液面的非線性運動以及邊界層紊流場結 構。模擬過程中，發展於波前相之寄生表面張力波列與波背相之沿流向條 痕結構與 Ebuchi et al. (1987)、Caulliez and Collard (1999) 以及 Veron and Melville (2001) 等風波實驗結果相符，而本論文對於此流場表面特徵的演化 與生成機制以及其內部流場結構皆有更詳盡地分析。
第二章
模式之數學理論
本模式模擬一包含非線性自由液面運動的不可壓縮黏性牛頓流之非靜 水壓-非穩態流場，其流場運動之控制方程式為質量守恆與動量守恆方程 式，邊界條件為自由液面運動條件與應力平衡條件、底部自由滑動條件以 及水平方向之週期條件。其中作用於自由液面上的應力包括壓力、黏滯應 力、表面張力以及由風場引起之正向應力與剪應力。非線性之自由液面運 動邊界條件與動力邊界條件皆真實地滿足於隨著「時間」與「空間」變動 之液面邊界上。由於自由液面的運動將使流場為一不規則形狀的區域空 間，本模式將應用空間座標轉換，將流場之控制方程式與邊界條件由包含 自由液面變動的不規則物理區間轉換至一立方體形的規則計算區間，並在 計算區間中進行流場速度與壓力場之求解。2.1 控制方程式與邊界條件
考慮一三維、非穩態之不可壓縮牛頓流體之流場，其底部為自由滑移 邊界，頂部為可運動之自由液面，側邊為週期邊界。在卡式座標系統{ }
xr =(x,y,z)中，x與y表示水平方向分量，z表示向上為正之垂直方向分量， 自由液面邊界表示為 z=η(x,y,t)；其平均水面位於 z =0，底部邊界為 h z =− 。流場中之主要變數為速度場{ }
ur ={u,v,w}；u(x,y,z)與v(x,y,z)為水平 速度分量，w(x,y,z)為垂直速度分量，以及壓力場 p(x,y,z)。流場之特徵長 度 L 與特徵速度 U 使用於無因次化控制方程式與邊界條件。無因次流場參數包括：Reynolds 數；Re，定義為慣性力與黏滯力之比 Re=UL/ν，Froude
數；Fr，定義為慣性力與重力之比 Fr =U/(gL)1/2，與 Weber 數；We，定義
為慣性力與表面張力之比 We=ρU2L/σ 。其中ν 為運動黏滯係數，g為重力
2.1.1 流場控制方程式 根據不可壓縮流之質量守恆，流場滿足無因次連續方程式如下： 0 = ∂ ∂ + ∂ ∂ + ∂ ∂ z w y v x u 。 (2.1) 流場之無因次動量方程式為Navier-Stokes 方程式，以保守形式表示如下： u Re x p z uw y uv x uu t u ( ) ( ) ( ) + 1 ∇2 ∂ ∂ − = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ ， (2.2) v Re y p z vw y vv x uv t v ( ) ( ) ( ) + 1 ∇2 ∂ ∂ − = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ ， (2.3) w Re z p z ww y vw x uw t w ( ) ( ) ( ) + 1 ∇2 ∂ ∂ − = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ 。 (2.4) 其中∇2為 Laplacian 運算元：∇2 =∂2/∂x2 +∂2/∂y2+∂2/∂z2， p 為動力壓力： 2 /Fr z P p= + ，其中P為總壓力，z/Fr2為靜水壓力。 2.1.2 自由液面的運動邊界條件 流場的頂部邊界為一隨時間變動的自由液面 η(x,y,t)，空間中描述此液 面為： 0 ) , , ( ) , , , (x y z t =z− x y t = F η 。 (2.5) 而自由液面之法線方向單位向量為： k j i F F n y x y x y y x x ˆ ) 1 ( 1 ˆ ) 1 ( ˆ ) 1 ( } ˆ { 2 1 2 1 2 1 2 2 2 2 2 2 + + + + + − + + + − = ∇ ∇ ≡ η η η η η η η η 。 (2.6) 自由液面之切面由兩個切線方向單位向量{tˆ1}與{tˆ2}組成，此兩向量之必要 條件為與自由液面之法線方向正交，如下所示： 0 } ˆ { } ˆ {t1 ⋅ n = ， (2.7) 0 } ˆ { } ˆ {t2 ⋅ n = 。 (2.8)
此兩切線向量並非唯一，因此我們定義{tˆ1}為與{nˆ}正交並投影於x-z 平面上 之自由液面切線單位向量，其表示為： k j i t x x x ˆ ) 1 ( ˆ 0 ˆ ) 1 ( 1 } ˆ { 2 1 2 1 2 2 1 ≡ + + + η + η η ， (2.9) 而{tˆ2}為與{nˆ}正交並投影於y-z 平面上之自由液面切線單位向量，其表示為： k j i t y y y ˆ ) 1 ( ˆ ) 1 ( 1 ˆ 0 } ˆ { 2 1 2 1 2 2 2 ≡ + + + η + η η 。 (2.10) 空間中自由液面隨時間產生形變但其始終保持為一個面，因此液面函數對 於時間的全微分必定為零： 0 ) , , , ( = Dt t z y x DF ，或表示為 0 = ∇ ⋅ + ∂ ∂ F u t F r 。 (2.11) 將自由液面表示式(2.5)式代入上式，並運用(2.6)式即可得： F n u t = ⋅ × ∇ ∂ ∂ } ˆ { } {r η 。 (2.12) 其物理意義為自由液面位置隨時間的變化率等於流場速度在自由液面法線 方向上之速度分量。將式(2.12)展開即為自由液面運動邊界條件： 0 = − ∂ ∂ + ∂ ∂ + ∂ ∂ w y v x u t η η η 。 (2.13) 2.1.3 作用於自由液面上的表面張力 當兩不互溶之液體接觸時會產生一密度不連續之交界面；例如空氣與 水之交界面，此兩種流體間之分子鍵結特性將使交界面處產生一張力作 用；亦即於界面間產生一壓力差，而張力將驅使界面形變為一曲面；如雨 滴或是水中的氣泡，此現象稱之為表面張力作用。此張力定義為表面張力 係數σ ；為單位長度所受的力，單位為N/cm或是N/m，其量值決定於交界 面 間 之 液 體 種 類 與 溫 度 。 在 20o C 下 ， 水 與 空 氣 界 面 之 表 面 張 力 N/m 073 . 0 = σ 。於自由液面流場中，當液面產生曲度時，流場界面出現由表
面張力引起之壓力差，此正向壓力差與表面張力係數以及液面曲率κ關係如 下： σκ = ∆p ， (2.14) 其中，自由液面之曲率κ為：
(
xx y yy x x y xy)
y x yy y xy x y x y xy y xx x y x x y x yy y x xx y x y x y y x x z y x n η η η η η η η η η η η η η η η η η η η η η η η η η η η η η η η η η η η η η κ 2 ) 1 ( ) 1 ( ) 1 ( 1 ) 2 2 ( ) 1 ( 2 1 ) 2 2 ( ) 1 ( 2 1 ) 1 ( ) 1 ( ) 1 ( 1 ) 1 ( ) 1 ( 2 2 2 2 2 2 2 2 2 2 2 2 2 / 1 2 2 2 / 1 2 2 2 / 1 2 2 2 3 2 3 2 3 2 1 2 1 + + − + − + + = + ⋅ + + ⋅ + + ⋅ + + ⋅ + + + − + + − = + + ∂ ∂ + + + − ∂ ∂ + + + − ∂ ∂ = ⋅ ∇ ≡ − − − − r 。 (2.15) 2.1.4 作用於自由液面上的風應力 風作用於流場之應力張量[τw]於自由液面上之應力分量，可表示為向量 型式{τrs}如下所示： T T [ ] {ˆ} } {τrs = τw ⋅ n ， (2.16) 其{nˆ}為應力作用面（亦即自由液面）之法線方向。因此，作用於自由液面 上之風應力又可分解為作用於nˆ方向之正向應力分量 s n τ ，以及分別作用於切 線方向ˆt1與ˆt2之剪應力分量 s t1 τ 與 s t2 τ ，表示如下： 2 1 ˆ ˆ ˆ } { 2 1t t n s t s t s n s τ τ τ τr = + + 。 (2.17) 其各分量分別為：{ }
[ ]
{ }
s n w n nˆ τ ˆ Τ =τ ， (2.18){ }
tˆ1[ ]
{ }
nˆ 1 2tˆ1 tˆ2 s t s t w = + ⋅ Τ τ τ τ ，{ }
[ ]
{ }
s t s t w n t t tˆ2 τ ˆ Τ =τ1ˆ1⋅ˆ2 +τ 2。 (2.19) (2.20) 2.1.5 自由液面上的應力平衡 自由液面的動力邊界條件即為應力平衡條件。而自由液面運動過程中 所受的應力包括：流體之總壓力P、流體黏滯應力、表面張力以及風所引起 之應力。流場運動過程中，流體所受之應力張量[T]為：[ ]
∂ ∂ + − ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ + − ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ + − = z w P y w z v x w z u y w z v y v P y u x v x w z u y u x v x u P T µ µ µ µ µ µ µ µ µ 2 2 2 ， (2.21) 其中µ為流體動力黏滯係數。而作用於自由液面上之應力在分別於沿著自由 液面法線方向{ }
nr 與切線方向{ }
tr1 和{ }
tr2 之分量須與表面張力以及風應力達到 平衡，如下所示：{ }
[ ]
{ }
s n n T nˆ ˆ T =−σκ+τ ，{ }
tˆ1[ ]
T{ }
nˆ 1 2tˆ1 tˆ2 s t s t + ⋅ = Τ τ τ ，{ }
[ ]
{ }
s t s tt t n T tˆ2 ˆ Τ =τ1ˆ1⋅ˆ2 +τ2。 (2.22) (2.23) (2.24) 將式(2.22)至(2.24)展開，並無因次化即可得自由液面動力邊界條件： s n 2 3 2 2 2 2 2 2 2 2 2 ) 1 ( 2 ) 1 ( ) 1 ( 1 ) 1 ( 1 2 τ η η η η η η η η η η η η η η η η η η + + + − + + + = ∂ ∂ + ∂ ∂ + ∂ ∂ − ∂ ∂ + ∂ ∂ − ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + + + − y x xy y x x yy y xx y x y y x x y x We z w y w z v x w z u y v y u x v x u Re p Fr ， (2.25) + + + + ⋅ ⋅ + + + = ∂ ∂ + ∂ ∂ − ∂ ∂ + ∂ ∂ − + ∂ ∂ + ∂ ∂ − ∂ ∂ − ∂ ∂ s t y x y x s t y x x y x x y x Re y w z v x w z u y u x v x u z w 2 1 2 1 2 1 2 1 2 2 1 2 2 2 2 2 ) 1 ( ) 1 ( ) 1 ( ) 1 ( ) 1 ( 2 τ η η η η τ η η η η η η η η ， (2.26) + + + + ⋅ ⋅ + + + = ∂ ∂ + ∂ ∂ − ∂ ∂ + ∂ ∂ − + ∂ ∂ + ∂ ∂ − ∂ ∂ − ∂ ∂ s t s t y x y x y x y y x y x y Re x w z u y w z v y u x v y v z w 2 1 2 1 2 1 2 1 2 2 1 2 2 2 2 2 ) 1 ( ) 1 ( ) 1 ( ) 1 ( ) 1 ( 2 τ τ η η η η η η η η η η η η 。 (2.27) 2.1.6 流場底部自由滑移邊界條件 流場底部之自由滑移條件如下： 0 = ∂ ∂ z u ， 0 = ∂ ∂ z v ， 0 = w ， 0 = ∂ ∂ z p 。 (2.28) (2.29) (2.30) (2.31)
2.2 座標轉換
由於我們所探討的流場是一個包含自由液面運動的流場，其物理區間 並不是一個規則幾何形狀區間，因此，我們使用座標轉換將原本包含自由 液面邊界的不規則「物理區間」轉換至一規則形狀的「計算區間」。 2.2.1 座標轉換關係式 此座標轉換定義物理區間(x,y,z,t)與計算區間(ξ,ψ,ζ,τ)之格點位置關 係，如圖2.1 所示：其關係式為： t = τ ， x = ξ ， y = ψ ， h t y x h z + + = ) , , ( η ζ 。 (2.32) (2.33) (2.34) (2.35) 由式(2.35)可知，此座標轉換與時間演化相依，而物理區間在垂直方向上之 邊界；z=η(x,y,t) 與 z=−h，即對應至計算區間之 ζ =1 與 ζ =0。在物理 區間中所有導數運算皆將運用連鎖法則轉換至計算區間，其一次導數之轉 換分別如下所示： ζ η ζ η τ ζ ζ τ τ τ ∂ ∂ + − ∂ ∂ = ∂ ∂ + ∂ ∂ = ∂ ∂ h t t t ， ζ η ζ η ξ ζ ζ ξ ξ ξ ∂ ∂ + − ∂ ∂ = ∂ ∂ + ∂ ∂ = ∂ ∂ h x x x ， ζ η ζ η ψ ζ ζ ψ ψ ψ ∂ ∂ + − ∂ ∂ = ∂ ∂ + ∂ ∂ = ∂ ∂ h y y y ， ζ η ζ ζ ∂ ∂ + = ∂ ∂ = ∂ ∂ h z z 1 。 (2.36) (2.37) (2.38) (2.39) 圖2.1：物理區間(a)經座標轉換至計算區間(b)之二維示意圖。 (a) (b)
2.2.2 計算區間之流場控制方程式與邊界條件 依據式(2.37)至(2.39)之轉換法則，物理區間之無因次連續方程式(2.1)經 座標轉換後於計算區間中之形式為： 0 1 = ∂ ∂ + + ∂ ∂ + − ∂ ∂ + ∂ ∂ + − ∂ ∂ ζ η ζ η ζ η ψ ζ η ζ η ξ ξ ψ w h v h v u h u 。 (2.40) 而Navier-Stokes 方程式為： ∂ ∂ + ∂ ∂ + ∂ ∂ + + ∂ ∂ + − ∂ ∂ − ∂ ∂ − ∇ = ∂ ∂ + − ∂ ∂ + ∂ ∂ ζ η ζ η ζ η η ζ ζ η ψ ξ ζ η ζ η ξ τ ψ ξ τ ξ ) ( ) ( ) ( 1 ) ( ) ( 1 2 uv uu u h uw h uv uu u Re p h p u ， ∂ ∂ + ∂ ∂ + ∂ ∂ + + ∂ ∂ + − ∂ ∂ − ∂ ∂ − ∇ = ∂ ∂ + − ∂ ∂ + ∂ ∂ ζ η ζ η ζ η η ζ ζ η ψ ξ ζ η ζ η ψ τ ψ ξ τ ψ ) ( ) ( ) ( 1 ) ( ) ( 1 2 vv uv v h vw h vv uv v Re p h p v ， ∂ ∂ + ∂ ∂ + ∂ ∂ + + ∂ ∂ + − ∂ ∂ − ∂ ∂ − ∇ = ∂ ∂ + + ∂ ∂ ζ η ζ η ζ η η ζ ζ η ψ ξ ζ η τ ψ ξ τ ) ( ) ( ) ( 1 ) ( ) ( 1 1 2 vw uw w h ww h vw uw w Re p h w 。 (2.41) (2.42) (2.43) 其轉換後之Laplacian 運算元為：
[
]
[
η η η η η]
η ζ ζ ζ ψ η ζ η ζ ξ η ζ η ζ η ζ η η ψ ξ ψψ ξξ ψ ξ ψ ξ ψ ξ ∂ ∂ + + + − + + ∂ ∂ ∂ + − ∂ ∂ ∂ + − ∂ ∂ + + + + ∂ ∂ + ∂ ∂ = ∇ 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ) ( ) )( ( ) ( 2 2 2 ) ( ) ( 1 h h h h h 。 (2.44) 自由液面運動邊界條件式(2.13)經過座標轉換後為： ψ ξ η η τ η v u w− − = ∂ ∂ ， (2.45)自由液面動力邊界條件式(2.25)至式(2.27)經過座標轉換後為： s n We w w v u v u Re Fr p τ η η η η η η η η η ψ η ξ η ξ ψ η η ψ η ξ η η η η ψ ξ ξψ ψ ξ ξ ψψ ψ ξξ ψ ξ ψ ξ ξ ψ ψ ξ − + + − + + + − ∂ ∂ − ∂ ∂ − ∂ ∂ + ∂ ∂ + ∂ ∂ + − ∂ ∂ + − + + + = 2 / 3 2 2 2 2 2 2 2 2 2 ) 1 ( 2 ) 1 ( ) 1 ( 1 ) ( ) 1 ( ) 1 ( ) 1 ( 2 ， (2.46) + + + + ⋅ ⋅ + + + = ∂ ∂ + + + + ∂ ∂ − ∂ ∂ − + ∂ ∂ − ∂ ∂ + + + + ∂ ∂ − ∂ ∂ − s t s t Re w h w w v u h u u 2 1 2 1 2 1 2 1 2 2 1 2 2 2 2 2 2 2 2 2 ) 1 ( ) 1 ( ) 1 ( ) 1 ( ) 1 ( ) 1 ( 1 2 τ η η η η τ η η η ζ η η η η ψ η η ξ η ξ η ζ η η η ψ η ξ η ψ ξ ψ ξ ψ ξ ξ ψ ξ ξ ψ ξ ξ ψ ψ ξ ψ ξ ， (2.47) + + + + ⋅ ⋅ + + + = ∂ ∂ + + + + ∂ ∂ − + ∂ ∂ − ∂ ∂ + + + + ∂ ∂ − ∂ ∂ − ∂ ∂ − s t s t Re w h w w v h v v u 2 1 2 1 2 1 2 1 2 2 1 2 2 2 2 2 2 2 2 2 ) 1 ( ) 1 ( ) 1 ( ) 1 ( ) 1 ( ) 1 ( 1 2 τ τ η η η η η η η ζ η η η η ψ η ξ η η ζ η η η ψ η ξ η ψ η ψ ξ ψ ξ ψ ξ ψ ψ ξ ψ ψ ψ ξ ψ ξ ψ ξ ξ 。 (2.48) 自由滑移之底部邊界條件式(2.28)至式(2.31)經過座標轉換後為： 0 = ∂ ∂ ζ u ， 0 = ∂ ∂ ζ v ， 0 = w ， 0 = ∂ ∂ ζ p 。 (2.49) (2.50) (2.51) (2.52)
2.3 能量守恆
流場能量守恆即為流場中各形態能量；包括動能、勢能、表面張力所 做的功以及黏滯能量消散之變化率總合必須等於外力對流場做功所輸入之 能量增加率。欲了解流場運動過程中能量分布情形，我們將速度向量與動 量方程式進行內積，即可得流場能量方程式，其推導過程詳述於附錄一。 所得之無因次能量方程式表示如下：(
2)
(
2)
(
2 2)
0 0 2 2 2 2 2 2 2 2 2 0 0 2 2 2 2 2 2 1 0 0 2 2 1 1 2 1 0 0 1 1 1 2 1 ) ( 1 2 1 ) ( 2 1 dS w v w u dS t dV x w z u y w z v x v y u z w y w x w z v y v x v z u y u x u Re dS t We dS t Fr dV w v u Dt D y x S S s t y y s t x x s n V S S V + + ⋅ + + + ⋅ + + + ∂ ∂ = ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ − ∂ ∂ + ∂ ∂ + + +∫∫
∫∫
∫∫∫
∫∫
∫∫
∫∫∫
η η τ η η τ η η τ η η κ η 。 (2.53) 其中S0為自由液面投射至x-y 平面之面積。等號左邊第一項為流體運動 引起之動能變化率（≡dEk dt），第二項為自由液面變動所產生之勢能變化 率（≡dEη dt），第三項為表面張力所做功引起之能量變化率（≡dEσ dt）， 第四項為由流體黏滯性所引起之能量消散變化率（≡dEν dt）。而等號右邊 第一項為自由液面上正向風應力做功產生之能量變化率（ dE s dt n τ ≡ ），第二 項為風剪應力做功產生之能量變化率（ dE s dt t τ ≡ ）。若流場運動中無任何外 力作用，則等號右邊即為零。第三章
模式之數值方法
3.1 網格點配置與空間導數的離散
本模式計算區間的邊界為水平方向（ξ 與ψ ）之四個週期邊界面以及 垂直方向之自由液面映射邊界：ζ =1，與底部邊界：ζ =0。模式中同時使 用兩種網格系統：正規網格系統與交錯網格系統。其中，正規網格系統在 水平方向為均勻分布，而垂直方向則為非均勻分布；亦即網格間為非等間 距，其分布型式如下( Gavrilakis，1992)： 8417 . 1 tanh 1 8417 . 1 tanh ) ( − × = z N k k ζ 。 (3.1) 其中Nz為垂直方向之總格點數，k為垂直位置指標。根據式(3.1)，在垂直方 向上越接近液面，空間解析度愈高，其網格點位置與間距關係顯示於圖3.1。 垂直方向非等間距的正規網格於物理區間與計算區間之分布，如圖 3.2 所 示，而此網格分布將更有利於解析自由液面邊界層。交錯網格系統之水平 位置與正規網格重疊，但其垂直方向則位移至相鄰正規網格之間距中心， 如圖3.3 所示。兩種網格系統之使用方式為：水平速度 u、v 與壓力 p 置於 正規網格點，垂直速度 w 則置於交錯網格點。模式中之連續方程式、壓力 Poisson 方程式，以及水平方向之動量方程式皆滿足於正規網格點上，而垂 直方向之動量方程式則滿足於交錯網格點上。 水平方向的導數：如∂u ∂ξ、∂ u2 ∂ξ2 、∂v ∂ψ 、∂ v2 ∂ψ2 與∂(uv) ∂ξ 等， 使用擬頻譜法進行離散，模式中將利用大量的快速傅立葉轉換來進行水平 導數的運算。有關於擬頻譜法之離散型式傅立葉級數理論將於附錄二中說 明。而垂直方向導數計算：如∂u ∂ζ 與∂ u2 ∂ζ2等，則使用二階中間差分法進 行離散。 模式中同時使用兩種網格系統的原因是為了維持連續方程式（同時亦為壓力方程式）之空間解析度。由於映射法中連續方程式的滿足需要壓力 在垂直方向的二次導數，因此必須將垂直方向之動量方程式取垂直梯度（其 詳細步驟將於下節中說明）。 16 32 48 64 80 96 112 128 0.003 0.006 0.009 0.012 0.015 0 0.2 0.4 0.6 0.8 1 圖3.1：計算區間之垂直方向非等間距網格分佈情形。ζk為第k個網格點位置(實 線)，∆ζk為第k個網格間距(虛線)。越接近上邊界ζ =1處之網格分佈越密。以垂直 方向網格點數NZ =129為例，最大網格距約為0.015；最小網格距約為 0.0015。 圖3.2：垂直方向非等間距之正規網格於物理區間(a)與計算區間(b)之二維示意圖。 k k ζ k ζ ∆ k ζ k ζ ∆ (a) (b)
若不使用交錯網格，則壓力方程式之空間解析度就無法維持在一個垂直網 格距。若將垂直動量方程式(2.4)經時間積分後簡化表示為： B z p A w + ∂ ∂ = ， (3.2) 其中A與B為速度場相關係數。而假設在垂直方向為正規網格分布的情形下 (參考圖 3.4(a) )，在垂直位置指標k之網格點上，其連續方程式之垂直速度 在垂直方向上的導數 (∂w ∂z) 使用中間差分法離散過程之示意情形如下： 圖 3.3：正規網格系統與交錯網格系統配置圖。模式中水平速度 u, v 與壓力 p 置 於正規網格點，垂直速度w 則置於交錯網格點。而連續方程式、壓力 Poisson 方程 式，以及水平方向之動量方程式皆滿足於正規網格點上，垂直方向之動量方程式 則滿足於交錯網格點上。 ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ k i u, ξ ζ ○ 正規網格 ▲交錯網格 k i u+1, k i u−1, ui+2,k 1 ,k+ i u ui+ k1, +1 ui+ k2,+1 1 , 1 + − k i u 2 1 ,k+ i w 2 1 , 1 + + k i w 2 1 , 2 + + k i w 2 1 , 1 + − k i w 2 3 ,k+ i w 2 3 , 1 + + k i w 2 3 , 2 + + k i w 2 3 , 1 + − k i w
[
]
C p p A p p A B z p A B z p A w w z w z k k k z k k k z k k k k z k k z k + ∆ − − ∆ − ∆ ≈ − ∂ ∂ − + ∂ ∂ ∆ ≈ − ∆ ≈ ∂ ∂ − − + + + − + + − + 2 2 2 1 2 1 2 1 2 1 2 1 1 1 1 1 1 1 。 (3.3) 其中C為速度場相關係數。由式(3.3)可知，其離散後所需位於 k+2、k 與 k-2 之壓力值已跨過了鄰近於k之網格點。因此，連續方程式之垂直方向空間解 析度將由∆z擴大至2∆z。而在垂直方向為交錯網格分布的情形下(參考圖 3.4(b) )，離散過程之示意情形如下： 圖3.4：於垂直方向使用(a)正規網格系統以及(b)共同正規網格與交錯網格系統之配 置示意圖。文中說明使用(b)網格配置法才能使壓力 Poisson 方程式之空間解析度維 持為∆z。 ○ ○ ○ ○ ▲ 正規網格 交錯網格 ○ ○ ○ ○ ○ ▲ ▲ ▲ k p 1 + k p 2 1 + k w 2 3 + k w ○ pk−1 ▲ ○ pk−2 2 2 , + + k k w p pk+2 2 1 − k w 2 3 − k w (a) (b) z ∆ z ∆ 2 z ∆ 1 1 , + + k k w p k k w p , 1 1 , − − k k w p 2 2 , − − k k w p[
]
C p p A p p A B z p A B z p A w w z w z k k k z k k k z k k k k z k k z k + ∆ − − ∆ − ∆ ≈ − ∂ ∂ − + ∂ ∂ ∆ ≈ − ∆ ≈ ∂ ∂ − − + + − − + + − + 1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 1 1 1 。 (3.4) 其離散後所需位於k+1、k 與 k-1 之壓力值，則恰好與二階中間有限差分法 之空間解析度相符。因此，垂直速度 w 與壓力 p 必須分別配置於交錯網格 系統與正規網格系統才可維持一致的空間解析度。3.2 Runge-Kutta 時間積分與壓力 Poisson 方程式
本節中我們將利用二階 Runge-Kutta 對 Navier-Stokes 方程式(2.41)至 (2.43)以及自由液面運動方程式(2.45)進行時間積分，並運用連續方程式之映 射法（Projection）得到壓力 Poisson 方程式。為達到最佳演算法，在此將 Runge-Kutta 積分計算流程分解為兩個步驟，表示如下：[
( , ) ( , )]
2 ) , ( ) , ( * 1 i i i i i i i u f u f u u u f u u u f u τ τ τ τ τ τ τ − ∆ + ∆ + = ⋅ ∆ + = = ′ ∗ + ∗ 。 (3.5) 上式之 ui為時間τi時之速度，ui+1為下一個時間離散點τi +∆τ 時之速度。在 積分第一個步驟中，先以ui與 f(τi,ui)求得一中繼速度u∗，第二個步驟再以u∗ 與 f(τi +∆τ,u∗)求得ui+1。有關於Runge-Kutta 積分式之演化與其低儲存量的 運用將詳細說明於附錄三。 為了便利於說明，我們將Navier-Stokes 方程式(2.41)至(2.43)表示為： FU p h p u + ∂ ∂ + + ∂ ∂ − = ∂ ∂ ζ η ζ η ξ τ ξ ， (3.6)FV p h p v + ∂ ∂ + + ∂ ∂ − = ∂ ∂ ζ η ζ η ψ τ ψ ， FW p h w + ∂ ∂ + − = ∂ ∂ ζ η τ 1 。 (3.7) (3.8) 其中 FU、FV 與 FW 為速度平流項與黏滯項，其皆為速度(u,v,w)以及自由 液面高度η的函數，如下所示： ∂ ∂ + ∂ ∂ + ∂ ∂ + + ∂ ∂ + − ∂ ∂ − ∂ ∂ − ∇ = ζ η ζ η ζ η η ζ ζ η ψ ξ ψ ξ τ ) ( ) ( ) ( 1 ) ( ) ( 1 2 uv uu u h uw h uv uu u Re FU ， ∂ ∂ + ∂ ∂ + ∂ ∂ + + ∂ ∂ + − ∂ ∂ − ∂ ∂ − ∇ = ζ η ζ η ζ η η ζ ζ η ψ ξ ψ ξ τ ) ( ) ( ) ( 1 ) ( ) ( 1 2 vv uv v h vw h vv uv v Re FV ， ∂ ∂ + ∂ ∂ + ∂ ∂ + + ∂ ∂ + − ∂ ∂ − ∂ ∂ − ∇ = ζ η ζ η ζ η η ζ ζ η ψ ξ ψ ξ τ ) ( ) ( ) ( 1 ) ( ) ( 1 2 vw uw w h ww h vw uw w Re FW 。 而自由液面運動方程式為： ψ ξ η η τ η v u w− − = ∂ ∂ 。 (3.9) 在Runge-Kutta 第一步驟中： + ∂ ∂ + + ∂ ∂ − ⋅ ∆ + = ∗ i i i i i i p FU h p u u ζ η ζ η ξ τ ξ ， (3.10) + ∂ ∂ + + ∂ ∂ − ⋅ ∆ + = ∗ i i i i i i p FV h p v v ζ η ζ η ψ τ ψ ， (3.11)
+ ∂ ∂ + − ⋅ ∆ + = ∗ i i i i p FW h w w ζ η τ 1 ， (3.12) ] [ i i i i i i w u v ψ ξ η η τ η η∗ = +∆ ⋅ − − 。 (3.13) 其式(3.13)中之η∗可由已知速度場
(
ui,vi,wi)
與自由液面高度ηi順利求得。但 i p 尚為未知數，因此欲求速度場(
u∗,v∗,w∗)
必須先求解壓力。而流場壓力可由 強制欲求解之速度場滿足連續方程式，即可演化成為壓力方程式，如下所 示： 0 1 = ∂ ∂ + + ∂ ∂ + − ∂ ∂ + ∂ ∂ + − ∂ ∂ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ζ η ζ η ζ η ψ ζ η ζ η ξ ψ ξ w h v h v u h u 。 (3.14) 此即為連續方程式之映射法。將式(3.10)、(3.11)與(3.12)代入上式如下： 0 1 1 = + ∂ ∂ + − ⋅ ∆ + ∂ ∂ + + + ∂ ∂ + + ∂ ∂ − ⋅ ∆ + ∂ ∂ + − + ∂ ∂ + + ∂ ∂ − ⋅ ∆ + ∂ ∂ + + ∂ ∂ + + ∂ ∂ − ⋅ ∆ + ∂ ∂ + − + ∂ ∂ + + ∂ ∂ − ⋅ ∆ + ∂ ∂ ∗ ∗ ∗ ∗ ∗ i i i i i i i i i i i i i i i i i i i i i i i i i i i i FW p h w h FV p h p v h FV p h p v FU p h p u h FU p h p u ζ η τ ζ η ζ η ζ η ψ τ ζ η ζ η ζ η ζ η ψ τ ψ ζ η ζ η ξ τ ζ η ζ η ζ η ζ η ξ τ ξ ψ ψ ψ ξ ξ ξ 。 (3.15) 將上式進行代數運算並整理其型式可表示為：i i i i i i i i i i i i i i i i i S p p h p h p h p h h p h h p p = ∂ ∂ ∂ + ∂ ∂ ∂ + − ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ − ∂ ∂ ∂ ∂ + + + + ∂ ∂ + + + ∂ ∂ + ∂ ∂ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ζ ψ η ζ ξ η η ζ ζ η η ψ ζ η η ξ ζ ζ ζ ζ η η ζ η η η η ζ η η ψ ξ ψ ξ ψ ξ ψ ψ ξ ξ 2 2 2 2 2 2 2 2 ) )( ( ) ( ) )( ( 1 。 (3.16) 其中
