間步驟，延散步驟(advection and diffusion step)及傳播步驟(propagation-step)，分別 求得 n + 1/2 與 n +1 時刻之流場，並利用隱式數值方法求解。延散步驟求解移流 項(advection terms)和擴散項(diffusion terms)，傳播步驟求解壓力項、底床剪應力 項和連續方程式。水理控制方程式先就時間部分離散如下：
變數；∆t = tn+1 −tn；n 表示 n∆t 時刻之已知變數；n + 1/2 表示在(n +1)∆t 與 n∆t 間
針對 n+1 時刻的水深值(D n+1)做線性化處理，且僅保留一階項，(3.8)式可改寫成
項則可將式(2.15)與式(2.16)改寫為：
Horizontal Diffusion in f
Dh Dh D
1 1 1
Horizontal Diffusion in f
Dh Dh D
其中 Horizontal Diffusion in ξ 與 Horizontal Diffusion in η 為水平方向之剪應力，
如第二章式(2.22)、式(2.23)所述，在此不加贅述。
式(3.10)與式(3.11)針對垂直方向格網，也就是圖 3.3 的 T、P、B 三個格點進
本模式採用控制體積(control volume)法的觀念來離散控制方程式之空間項，
控制體積法的基本概念如圖 3.1 所示，其中(a)圖為實際區域，(b)圖為計算區域，
E、W、N、S 表相鄰格點，e、w、n、s 表控制面。模式計算之變數則放置在交 錯網格(staggered grid)上，如圖 3.2。在控制方程式中，除了移流項採用一階精度
混合型上風法(hybrid upwind scheme) (Spalding 1972)差分外，所有空間差分均採 中央差分法。至於移流效應重要性的判斷，則採用格網雷諾數(mesh Reynolds number)Rx、Ry作為判斷的因子，當|Rx|或|Ry|大於 2 時，代表移流效應重要，差分
0 2 所示，E、W、N、S、T、B 為相鄰格點，e、w、n、s、t、b 為控制面。模式計 算之變數亦放置在交錯網格(staggered grid)上，如前述圖 3.2。(3.10)及(3.11)等式 左邊採用 Crank-Nicolson method，而等式右邊Mu和Mv的空間差分則採用中央差
下，並無法適用在過於劇烈變化的底床，否則誤差累積後，數值穩定性將受到影
時間間距的計算與流速差異量相同。 表粒徑之比例在作用層中的和為 100%，以結合演算法(Spasojevic and Holley, 1990)概念求解如下：
式(3.26)與式(3.27)式分別對應式(2.40)與式(2.41)，均為一非線性代數式，加 以線性化後，利用 Newton-Raphson 法疊代求解：
1 1
式(3.28)與式(3.29)中，[∂F/∂s]為 Jacobian 係數矩陣中之列向量；lsn+1為前
[ ]
採用中央差分法離散式(3.28)與式(3.29)的空間項後，可分別表示為：
( ) ( ) ( ) ( )
3.3 模式演算流程模式演算流程模式演算流程 模式演算流程
本研究數值模式演算流程如圖 3.4 所示。
1. 首先藉由水平二維模式計算水深(式 3.9)與水深平均流速，此部分包含輸入 邊界條件，考慮流變關係對底床剪力之影響以及計算水體中因含砂濃度產 生的密度差異量。水深平均流速先以式(3.4)、式(3.5)計算延散步驟之流速，
再由式(3.6)、式(3.7)計算傳播步驟之流速。
2. 以流速差異量模式，求解式(3.10)、式(3.11)得到三維流場之流速差異量。三 維流場之流速係以u = +uɶ u 、v = +vɶ v 得到。
3. 以式(2.1)之連續方程式計算w。
4. 以式(2.12)、式(2.13)與式(2.14)計算有效剪應力Tij，以反饋深度平均方程式。
5. 若水理部分尚未收斂，或是在變量流過程尚未達到下一個計算階段，則回 到步驟 1 繼續計算。
6. 以式(3.34)求解質量傳輸方程式。
7. 藉由式(3.31)與式(3.32)之概念，聯立求解式(3.37)與式(3.38)，以解析作用層 中各個代表粒徑的百分比與底床沖淤量。
8. 若尚未達到計算時間則回到步驟 1，否則結束。
(a)
(洪聖翔, 2011)
圖 3.2 交錯格網示意圖
i、j分別代表水平格網上任一點之縱向及橫向位置；k代表在垂向格網上的位置(洪聖翔, 2011)
圖 3.3 三維流場與濃度離散之控制體積法示意圖(計算域)
&
u uɶ
&
v vɶ
& s w z
uɶ w
圖 3.4 數值模式演算流程圖