第三章 模式之數值方法
3.5 數值方法之驗證
本數值模式的建立是基於求解三維流場的 Navier-Stokes 方程式與連續 方程式，並伴隨著非穩態的自由液面運動與動力邊界條件。由於我們欲模 擬的物理問題為海洋表面下邊界層流場，其空間上具有在水平方向尺度遠 大於垂直方向尺度的特性，因此設定水平方向為週期邊界。而數值方法中，
空間之水平方向導數運算以擬頻譜法離散，而垂直方向導數則以二階中間 差分法離散。二階Runge-Kutta 應用於控制方程式與自由液面運動邊界條件 之時間積分。自由液面動力邊界條件分別運用於壓力方程式之邊界條件與 自由液面水平方向動量方程式中。壓力Poisson 方程式得以使流場達到質量 守恆。本節將以模擬一二維重力波之傳遞驗證數值離散方法的正確性與模 式的可用性。此模擬之各項參數說明如下：前進波波長λ =10cm，初始波形 陡度ak =0.25；a為波幅，k =2π λ為波數。模擬物理區間為水平方向長度
cm
10 ，而平均水深8cm 之流場。流場之特徵長度定義為L≡1/k，特徵速度 定義為波形傳遞之相速度U ≡ g k 。流場之無因次參數包括：Froude 數
=1
=U gL
Fr ，Reynolds 數Re=ULν ≈6286。自由液面初始高度為一三階 Stokes 波波形，而初始流場速度為 Stokes 波之速度場，並忽略表面張力的 影響。
模擬過程之平均水面高度，流場之平均速度散度以及最大速度散度分 別顯示於圖3.6(a) 至圖 3.6(c)。圖中顯示水平方向網格數與垂直方向網格數 分別為16×48、32×64以及128×128之不同解析度的模擬結果。以解析度為
128
128× 為例，當模擬結束時，其平均水面高度已達到數值計算之機器精確 度：k3
∫∫
ηdxdy ≈10−15。流場之平均速度散度：2
2
) 1
(gk − ∇⋅vr L ，在整個計算過 程中皆保持在小於10−11，而最大速度散度：
⋅ ∞
− ∇ v L
gk) 12 r
( ，亦小於10−8。其
平均水面高度的守恆顯示模式中處理自由液面邊界條件之離散方式符合流 場之物理特性。模擬過程中流場之質量守恆是取決於壓力Poisson 方程式之 收斂條件，而本模式皆將疊代之收斂條件設定為10−6。因此，流場之速度散 度確認質量守恆，亦驗證了求解壓力Poisson 方程式之正確性。
t/T0
Maximumvelocitydivergence
0 1 2 3 4 5 6
10-13 10-12 10-11 10-10 10-9 10-8 10-7 10-6
(c)
t/T0
Errornormofvelocitydivergence
0 1 2 3 4 5 6
10-13 10-12 10-11 10-10 10-9 10-8 10-7 10-6
(b)
t/T0
Meanfreesurface
0 1 2 3 4 5 6
10-20 10-19 10-18 10-17 10-16 10-15 10-14 10-13
(a)
圖 3.6：模擬過程之平均水面高度(a)、流場之平均速度散度(b)以及最大速度散度
(c)。圖中顯示不同解析度的結果：其中(Nx,Nz) =(16,48)為點線、(32,64)為折線、
) 128 , 128
( 為實線，N 為水平方向網格數；x Nz為垂直方向網格數。T 為前進波之線0
性週期。
圖3.7 (a) 顯示模擬過程中總能量之變化，其中包括流場動能、勢能以 及因黏滯消散所損失的能量。圖 3.7(b) 顯示模擬過程中數值能量損失與黏 滯消散能量之比。圖中顯示經過模擬初期調整，總能量之演化隨時間呈現 穩定，並隨著空間解析度增加而減少。以解析度為128×128為例，當模擬結 束時，由數值計算引起之能量損失為流場初始能量之0.03%。流場之能量守 恆驗證了求解Navier -Stokes 方程式之數值方法的正確性。而圖 3.7(b) 顯示 模擬過程中數值能量損失與黏滯消散能量之比約為0.01，其顯示本模式應用 的數值方法在足夠的解析度下具有解析流場邊界層的能力。
圖3.8 為自由液面上之波形演化模擬結果並顯示其與 Fenton(1988) 高 階近似理論解之比較。Fenton(1988) 之穩態高階 Stokes 波近似解是由解析 非黏性勢流場理論之Laplace 方程式所得，其自由液面邊界條件為 Bernoulli 方程式。而本節中所模擬之前進波，其動力機制為重力與其慣性力，黏滯 效應非常微弱。因此，由圖3.8 可知，模擬波形與 Fenton(1988)近似解非常 相近，表示本模式以求解Navier-Stokes 方程式並以自由液面運動方程式與 應力平衡為邊界條件亦能精準地反應自由液面之非線性運動特性。
t/T0
[E(t)-E(0)]/E(0)
0 1 2 3 4 5 6
10-8 10-7 10-6 10-5 10-4 10-3 10-2 10-1
(a)
t/T0
[E(t)-E(0)]/E(t)
0 1 2 3 4 5 6
10-5 10-4 10-3 10-2 10-1 100
(b)
ν
圖3.7：模擬過程中總能量之演化(a)以及數值能量損失與黏滯消散能量之比(b)。其
中總能量包括流場動能、勢能以及因黏滯消散所損失的能量。圖中顯示不同解析 度的模擬結果：其中(Nx,Nz) =(16,48)為點線、(32,64)為折線、(128,128)為實線，
N 為水平方向網格數；x Nz為垂直方向網格數。T 為前進波之線性週期。 0
xk
ηk
-3 -2 -1 0 1 2 3
-0.4 -0.2 0 0.2 0.4
128x128 32x64 16x48 Fenton
圖3.8：前進波之波形模擬結果。圖中顯示不同解析度於時間 t T =5.0 的模擬結 果：(Nx,Nz)=(128,128)為點線；(32,64)為折-點線；(16,48)為折線、Fenton(1988) 非
線性理論穩態解為實線。其N 為水平方向網格數；x Nz為垂直方向網格數，T 為前
進波之非線性週期。