方程式離散化

為傳統之完全隱示法(fully implicit method)。

dt

最簡單的推估方式包含算數平均(式 3.2.3-4a)、幾何平均(式 3.2.3-4b) 與調和平均(3.2.3-4c)等，因此上述三式之變數(x)若於達西公式中，即

代表水力傳導係數(K)。此外，傳統上兩串連之非均質區塊，其等效

( ) (

^α

) ( )

^τ

cond Area N

L

解。

圖 3.2.4-1 為內迭代執行流程圖，圖中之變數(V )若於地下水流問 題中，其代表壓力水頭；若於熱流問題中，其代表單位熱容量。首先，

給予變數初始解作為微分形式解法之第一步驟，以初始解作為初始搜 尋點。其次代入方程式集合中，最後透過各自之守恆方程式，意即水 流問題之質量守恆方程式與熱流問題之能量守恆方程式，可以算出變 數數值對應之守恆誤差。接著再透過差分方式，進一步計算該守恆誤 差之一階微分值。最後，判斷是否已經滿足收斂標準，若尚未收斂則 依據微分值修改變數數值，反覆執行此一動作，直至完全滿足內迭代 收斂標準為止。

在收斂條件方面，其收斂條件共有三種，滿足其一即停止內迭 代演算。首先是迭代次數之最大限制量，其次是守恆誤差之一階微分 值趨近於零，最後則是守恆誤差本身小於設定之內迭代收斂標準。

內迭代起點

取得節點[i]之初始搜尋值

計算守衡誤差

差分，計算守衡誤差微分值

是否符合收斂標準

內迭代收歛 n V⁽ⁿ⁾,

)

V(

' _n

∂

= ∂ε ε

1

' ,

) ( ) ( ) 1

( ⁺ = + n=n+

V V

V ^{n n} α ε_n

)

V(

' _n

∂

= ∂ε ε

圖 3.2.4-1 模式內迭代流程圖 3.2.5 資訊同步流程

在資訊同步流程部分，於前述之概念模式建立中，水流熱流偶 合模式是由地下水流模式及熱流模式所組成，從概念模式中即可瞭 解，兩模式彼此交互影響。其中水流模式可以計算水流之質量穿越 量，其數值則與熱流模式中之對流項估算影響極大；而熱流模式所計

算之溫度變化，則會影響水流之膨脹與壓縮。圖 3.2.5-1 為兩模式之 變數關係圖，當每次完成內迭代過程後均會將互相同步彼此所計算之 數值。

圖 3.2.5-1 地下水流模式與熱流模式之變數關係圖

3.2.6 外迭代處理方法

內迭代流程是透過最佳化方法求得各節點數值，其概念是假設 周遭節點之數值為正確，然而實際上僅有緊鄰邊界之計算節點，其周 遭邊界因為透過邊界條件之設定（Dirichlet B.C.），其數值方為正確。

意即內迭代過程僅訂定了單一節點與相鄰節點之關係方程式，如何求 得一組解，可以同時滿足聯立的關係方程組，在傳統上多採用矩陣解 法，來求得各節點對應之數值。因此，本研究採取外迭代流程求得各 計算節點之數值，其外迭代流程如圖 3.2.6-1 所示，反覆執行各節點 之資訊同步流程及內迭代流程，每次外迭代流程均可將整體數值往真 解逼近，當節點數值每次之迭代改變量漸漸趨緩，即表示已經趨近真 解。因此若其改變量小於外迭代收斂標準，則認定外迭代流程已經收

斂。

是 否

外迭代起始

同步取得資訊 [i]

水流節點[i]

內迭代

檢查最大結點數值改變幅度 小於外迭代收歛標準

外迭代收歛 同步取得資訊

[ii]

同步取得資訊 [n]

熱流節點[i]

內迭代

水流節點[ii]

內迭代

熱流節點[ii]

內迭代

水流節點[n]

內迭代

熱流節點[n]

內迭代

圖 3.2.6-1 外迭代流程圖

3.2.7 整體數值模擬流程

整體數值模擬流程如下圖 3.2.7-1，首先依據模式設定檔，讀入 空間切割相關資訊、水文地質參數、邊界條件與方程式集合等資訊，

接著則依據模擬之模擬型態、起始時刻、結束時刻與模擬間距，開始 進行模擬，流程中則當外迭代流程收斂後，則進行是否進行下一時刻 之模擬，若否則完成模擬。若穩態形式之模擬，則僅執行外迭代流程 一次；若非穩態形式之模擬，則依據起始時刻、結束時刻與模擬間距 等資訊進行判斷。

圖 3.2.7-1 整體數值模擬流程圖

四、垂向二維地下水流與熱流偶合數值模式開發

1. Func_Porosity_Affeected_by_Pressure (斜體字代表方程式實做函式 名稱)

2. Func_Character_PS

式 4.1-2 與式 4.1-3 是以 van Genuchten(1980)所提出之壓力水頭 與含水量的特性曲線關係式為基礎，修改為可同時應用於飽和與未飽 和之地下水流問題。未飽和問題中，透過壓力水頭與土壤參數計算對 應之土壤飽和度；飽和問題中，其數值等於土體孔隙率。本函式為節

(

r

)

3. Func_Calc_Water_Density

式 4.1-4 為引用 Rana A.與 Frank J. Millero 在 1973 年經實驗迴歸

4. Func_Set_Temperature_Const

本函式用以設定地下水流溫度，本函式適用於地下水流問題非 偶合運算。對於偶合運算時，則不使用本函式，其數值則透過熱流模 式運算，並透過資訊同步過程取得。本函式為節點形式之方程式。

5. Func_Calc_TotalHead_DensityDependent

本函式用於計算總水頭，本方程式可應用於變密度條件下，其

0

0

γ

∫

γ +

=

z z dz p

h ... (4.1-5)

6. Func_Calc_TotalHead_DensityIndependentl

本函式用於計算總水頭，本方程式僅可應用於定密度條件下，

其位置水頭之估算可直接以位置高程訂定之，再與壓力水頭之總和即 為總水頭。本函式為節點類型之方程式。其中 為總水頭、 為位置 水頭、 為壓力水頭。

h z

p p z

h= + ... (4.1-6)

7. Func_WaterCapacity

式 4.1-7 則用以將土壤含水量轉換為控制體積內的蓄水質量。本 函式為節點類型之方程式。本函式為節點類型之方程式。其中M& 為控 制體積內之蓄水質量、ρ_f 為流體密度、θ為含水量、Vol則為控制體 積。

Vol

M& =(ρ_fθ) ... (4.1-7)

8. Func_InterporlateDensity_UpWind

本函式透過連結之兩端點水流密度值，以上風法之概念推估連 結中央之水流密度，其數值可代表控制表面上的水流密度。本函式為 連結類型之方程式。

9. Func_InterporlateK

本函式透過連結之兩端點實際水力傳導係數，以調和平均推估 連結中央之實際水力傳導係數，其數值可代表控制表面上的水力傳導 係數。本函式為連結類型之方程式。

10. Func_DarcyLaw_Flux

11. Func_Continuity_Universal

本函式為地下水模式中之守恆方程式，依據模擬問題可分為穩 態與暫態兩種模擬方式，其中暫態模擬又可依據顯示法(explicit method)、完全隱示法(fully implicit method)與混合隱示法

(Crank-Nicholson method)，建立不同時間項之處理方式。

在穩態模擬上，連續方程式如式 4.1-10 所示，總穿越量與源項

(

α

)

τ

Func_Interpolate_Density 與 Func_InterporlateDensity_UpWind 都屬於 水流密度之空間推估函式，為相同功能但內部處理方法不同之函式，

因此可以依據需求擇一使用。此外，方程式 Func_Character_PS 為 van Genuchten 所提出之特性曲線模型，因此未來亦可因應不同之需求，

選用其他特性曲線模型，例如 Brook 等人所提出之方法。

4.2 熱流數值模式實做

1. Func_Calculate_Flow_SpecificHeat

本函式用以設定地下水流流體之比熱，由於本問題僅探討純水 之問題，故在此將地下水流流體之比熱訂定為常數值 4179(J kgK)。

未來若問題延伸至污染傳輸議題，流體比熱則應隨著污染值濃度而有 所變化，屆時可再行擴充。本函式為節點類型之方程式。

2. Func_Set_Porosity_const

本函式用以設定材質之孔隙率，本函式適用於熱流問題非偶合 運算。對於偶合運算時，則不使用本函式，其數值則透過水流模式運 算，並透過資訊同步過程取得。

3. Func_Set_Flow_Density_const

本函式用以設定地下水流流體之密度，其數值固定為

3

則不使用本函式，其數值則透過水流模式運算，並透過資訊同步過程 取得。本函式為節點類型之方程式。

4. Func_Set_FlowMassFlux_const

本函式以設定地下水流流體之質量穿越量，其數值固定為 0( )，本函式適用於熱流問題非偶合運算。對於偶合運算時，則 不使用本函式，其數值則透過水流模式運算，並透過資訊同步過程取 得。本函式為連結類型之方程式。

day kg /

5. Func_Calculate_Flow_HeatConductivity_const

本函式以設定地下水流流體之熱傳導係數，其數值固定為 48,038.4(^{J K}_day⋅m)。未來若搭配污染傳輸模式，地下水流流體之熱 傳導係數設定則可隨污染傳輸模式結果而變化。本函式為節點類型之 方程式。

6. Func_InterporlatePorosity

本函式透過連結兩端點之孔隙率數值，以算數平均推估連結中 央之孔隙率，其數值可代表控制表面上的孔隙率。本函式為連結類型 之方程式。

7. Func_Interporlate_Flow_Density

本函式透過連結之兩端點水流密度值，以算數平均推估連結中 央之水流密度，其數值可代表控制表面上的水流密度。本函式為連結 類型之方程式。

8. Func_Heatcapacity_Transfer_To_Temperature

式 4.2-1 是單位水體熱容量與溫度的轉換關係式，透過單位體積 之水體所含熱容量與水體比熱做轉換可得此單位水體所相對應之溫 度。本函式為節點類型之方程式。

t i i f H t

i s T

uh = _{, ,} × ... (4.2-1) 其中uh為單位水體之熱容量、s_H,f 為流體比熱、T為水體之熱容量之 溫度。

9. Func_CalcMaterial_Constant

式 4.2-2 是計算控制體積內之多孔介質之質量，因本模式有考慮 到多孔介質之壓密性，然而此壓縮性是假定反映在孔隙率上，即在控 制斷面上無多孔介質之質量流量，固在此計算在控制體積內之多孔介 質之質量僅需考慮受壓後之孔隙率與多孔介質密度之乘積即可。本函 式為節點類型之方程式。

Vol n

M_s =(1− )*ρ_s* ... (4.2-1) 其中M_s為多孔介質質量、n為孔隙率、ρ_s為多孔介質密度、Vol為控制 體積。

10. Func_ Temperature y_Transfer_To_Total Heatcapacit

式 4.2-1 是總熱容量與溫度的轉換關係式，首先透過流體與多孔 介質之質量計算土體之等效比熱，再乘上控制體積數值與溫度求得對 應的總熱容量。本函式為節點類型之方程式。

( ) [ ( ) ]

{

^{n s n s}

}

^{Vol T}

H = ×ρ_f × _H,f + 1− ×ρ_s× _H,s × × ... (4.2-1) 其中H為總熱容量、 為孔隙率、Vol為控制體積、T 為多孔介 質與流體之平均溫度、

n

ρf 為流體密度、s_H,f 為流體比熱、ρ_s為材質密 度、s_H,s為多孔介質比熱。

11. Func_Calculate_EfficientHeatK

式 4.2-2 是以 Bear(1988)所提出轉換方程式，以土體孔隙率、水 流熱傳導係數與多孔介質熱傳導係數，計算土體之等效熱傳導係數。

其中 為等效熱傳導係數、n為孔隙率、 為流體熱傳導係數、 為多孔介質熱傳導係數。

Ke K_H,f K_H,s

(

H f

) ( ( )

Hs

)

eq

H n K n K

K _, = * _, + 1− * _, ... (4.2-2)

12. Func_Interporlate_Flow_SpecificHeat

本函式透過連結兩端點之水流比熱數值，以算數平均推估連結 中央之水流比熱，其數值可代表控制表面上的水流比熱。本函式為連 結類型之方程式。

本函式透過連結兩端點之水流比熱數值，以算數平均推估連結 中央之水流比熱，其數值可代表控制表面上的水流比熱。本函式為連 結類型之方程式。

在文檔中 應用細胞自動機於垂向二維地下水流與熱流偶合數值模式之建立 (頁 31-0)