五、 垂向二維水流與熱流偶合數值模式模擬與驗證
5.1 熱流案例驗證與模擬
本節為非偶合案例之展示，將分為五個不同的案例進行數值驗 證與討論。案例 5.1-1 為採用定溫 Dirichlet 邊界均質案例；案例 5.1-2 為同時採用定熱通量 Neumann 邊界與定溫 Dirichlet 邊界均質案例；
案例 5.1-3 與 5.1-4 均為非均質案例模擬，分別為水平向材質分割與 垂直向材質分割；案例 5.1-5 則另外增加一個沈流點(sink point)，用 以驗證源流項(source term)或沈流項(sink term)之均質案例；案例 5.1-6 設定同案例 5.1-5，差異處在於其網格在抽取點附近進行網格加密；
案例 5.1-7 為案例 5.1-6 配置之暫態模擬。表 5.1-1 為非偶合案例列表。
表 5.1-1 非偶合案例列表
案例 5.1-1 Dirichlet B.C.
案例 5.1-2 單一材質
Dirichlet + Neumann B.C.
案例 5.1-3 材質水平向切割
案例 5.1-4 複合材質
材質垂直向切割 案例 5.1-5
規則網 格
單一材質 sink point 案例 5.1-6
穩態
sink point 案例 5.1-7 暫態
不規則
網格 單一材質
sink point
1. 案例 5.1-1 案例說明：
整體模擬區域如圖 5.1-1 所示，為 11 公尺見方之垂向二維方形 薄板。在網格切割上，由原點開始規則地每 1 公尺配置一個計算節 點，透過 Voronoi Diagram 可以產生 121 個(11×11)規則形狀之矩形網 格，因此每網格之長與寬均為 1 公尺。在材質方面，整體區塊設定為 相同材質，其中土壤孔隙率為 0.38，土壤熱傳導係數為 172,800 (J Kday m
⋅ )，水流熱傳導係數為 48038(
m Kday
J ⋅ )，土壤密度為 2,323(kg/m3)，土壤比熱為 710(
K kg
J )，水流密度為 1000( )。
此外，由於本節案例僅單獨模擬熱流部分，並不模擬水流部分，因此 令水流流速為 0，故模擬上僅有傳導項而無對流項。
/m3
kg
在邊界條件設定上，左方邊界溫度設定為 70(℃)，右方邊界溫 度設定為 25(℃)，上下方邊界則設定為無熱流量邊界(No Heat Flux Boundary)，因此本案例溫度將由左往右逐步遞減，亦即熱流方向應 為由左往右傳遞。本案例為穩態模擬，全體網格之初始猜值雖任意設 定下，最終仍可收斂至真解，但過於懸殊之設定值易影響模擬收斂性 與效率，在此設定為 50(℃)。透過模式求解，網格數值逐步依據邊界 條件之設定值，收斂於前述之熱流型態。
數值結果：
圖 5.1-2 為溫度分佈之模擬等值圖，從圖面看來溫度由左向右逐 步遞減，其等值線分佈均勻，顯示遞減趨勢趨於直線。由於案例邊界 條件設計極為簡單，因此流量估算解析解即為傅立葉公式，將左右兩 方之邊界溫度差值、左右兩端距離與熱傳導係數帶入公式中即可求 得，其數值為 5,642,576.64(
Jday)。
表 5.1-2 為當外迭代收斂標準為 時，解析解與數值解於左 右兩側邊界之穿越熱量比較表，數值解方面是將左右兩側網格所計算 之穿越流量累加所得，由表上所示其數值誤差已達 0.002%。
10-7
1×
表 5.1-2 左右邊界穿越流量比較表(往左為正) 左側邊界流量(
Jday) 右側邊界流量(
Jday) 數值解 5,642,472 (99.998%) -5,642,681 (100.002%) 解析解 5,642,576.64(100%) -5,642,576.64 (100%)
0 1 2 3 4 5 6 7 8 9 10
0 1 2 3 4 5 6 7 8 9 10
溫度 25℃
溫度 70℃
圖 5.1-1 案例 5.1-1 格網與模式配置圖
0 1 2 3 4 5 6 7 8 9 10 0
1 2 3 4 5 6 7 8 9 10
圖 5.1-2 溫度分佈等值圖 2. 案例 5.1-2
案例說明：
由於 Dirichlet 邊界與 Neumann 邊界為兩種最基本之邊界條件，
本案例用以驗證 Neumann 邊界正確性，因此大部分設定均沿用案例 5.1-1 設定，左方邊界之熱通量設定為 5,642,576.64 ( )，該數值 即為案例 5.1-1 數值解計算所得之熱通量，右方與上下兩邊界則沿用 案例 5.1-1 設定，右方邊界設定為 25(℃)，上下方邊界則設定為無熱 流量邊界(No Heat Flux Boundary)，整體模擬區域如圖 5.1-3 所示，本 案例溫度分布應與案例 5.1-1 依樣，由左往右逐步遞減，亦即熱流方 向應為由左往右流動。
day J /
數值結果：
本案例主要是為驗證 Neumann 邊界條件之模擬正確性，因此與 案例 5.1-1 設定即為相似，僅在邊界條件上取不同之邊界條件型式，
因此模擬結果兩者應該一致。圖 5.1-4 為案例 5.1-1 與案例 5.1-2 之溫 度分佈圖，由圖面顯示兩者溫度分佈等值線完全重合，顯示雖然設定 為不同邊界條件型式，但其模擬結果仍然一致，可以證明 Neumann 邊界條件型式之正確性。
0 1 2 3 4 5 6 7 8 9 10
0 1 2 3 4 5 6 7 8 9 10
溫度 25℃
熱通量
512961 J/day
圖 5.1-3 案例 5.1-2 格網與模式配置圖
00 1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
:為案例5.1-1之溫度等值圖 :為案例5.1-2之溫度等值圖
圖 5.1-4 溫度分佈等值圖(案例 5.1-1 與案例 5.1-2) 3. 案例 5.1-3
案例說明：
本案例探討主題為非均質與非均質含水層差異之探討，因此除 了材質以外的方面均與案例 5.1-1 一致，例如網格切割、邊界條件等，
材質方面則區分為兩個區塊，圖 5.1-5 為本案例之格網與材質配置 圖，其中材質 1 與案例 5.1-1 之性質相同，孔隙率為 0.38、密度為 2323(kg/m3)、熱傳導係數為 172800(
m Kday
J ⋅ )、比熱為 710( )；
材質 2 大部分的性質均與材質 1 相同，其熱傳導係數為 691200(
/m3
kg
m Kday
J ⋅ )，約為材質 1 之熱傳導係數之 4 倍，在水流熱傳導 係數部分，其數值為 48038(J Kday⋅m)。
數值結果：
圖 5.1-6 為案例 5.1-1 與案例 5.1-3 之溫度模擬等值圖，圖 5.1-7
為兩者之總熱容量模擬等值圖。由於右邊區塊之熱傳導係數左邊區塊 之 4 倍大，因此左邊需要透過更高之溫度梯度方可維持相同的熱通 量。因此左方溫度變化較右方溫度變化緩和。若觀察圖 5.1-7，由於 兩區塊材質比熱相等，因此總熱容量分佈與溫度分佈呈正比，故分佈 上趨勢相同。
在熱能守衡分析部份，左邊界之熱通量為 8,619,116 (J day)，右 邊界之熱通量為 8,619,490 (J day)，兩者之相對誤差為 0.0043%，符 合熱能守恆之條件。
0 1 2 3 4 5 6 7 8 9 10
0 1 2 3 4 5 6 7 8 9 10
材質一 材質二
溫度 25℃
溫度 70℃
圖 5.1-5 案例 5.1-3 格網與材質配置圖
:為案例5.1-3之溫度等值圖 :為案例5.1-1之溫度等值圖
0 1 2 3 4 5 6 7 8 9 10
0 1 2 3 4 5 6 7 8 9 10
圖 5.1-6 溫度分佈等值圖(案例 5.1-1 與案例 5.1-3)
00 1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
圖 5.1-7 總熱容量分布等值圖(案例 5.1-3)
4. 案例 5.1-4 案例說明：
本案例探討主題亦為非均質與非均質含水層差異，上下兩部分區 隔成不同材質，本案例之網格切割、邊界條件、初始條件等設定與案 例 5.1-1 相同。在材質設定方面，其中材質 1 與案例 5.1-1 之性質相 同，孔隙率為 0.38、密度為 2323( )、熱傳導係數為 172,800 (
/m3
kg
m Kday
J ⋅ )、比熱為 710( )；材質 2 大部分的性質均與材質 1 相同，其比熱為 860( )，約為材質 1 之比熱之 1.2 倍，在水流 熱傳導係數部分，其數值為 48038(
K kg J /
K kg J /
m Kday
J ⋅ )。圖 5.1-8 為本案例之 材質設定圖。
數值結果：
圖 5.1-9 為案例 5.1-4 之溫度模擬等值圖，由於溫度遞減與熱傳 導係數有關，由於材質一與材質二之熱傳導係數相同，該數值又與案 例 5.1-1 之熱傳導係數相同，因此兩案例之溫度坡降型式相等，從圖 5.1-9 所示，兩者之等值線重合。圖 5.1-10 為兩者之總熱容量，總熱 容量為溫度、控制體積與對應之等效比熱之乘積，由於兩區塊之材質 比熱不同，因此會有不同之數值，原則上其仍與溫度分佈呈現正比關 係。
而在熱能量守衡部份，本案左邊界之流入總量為 5,642,476 (J day)，右邊界之流出總量為 5,642,677(J day)，兩者之相對誤差為 0.003562266%，符合熱能守恆之條件。
0 1 2 3 4 5 6 7 8 9 10
0 1 2 3 4 5 6 7 8 9 10
材質二
材質一
溫度 25℃
溫度 70℃
圖 5.1-8 案例 5.1-4 格網與模式配置圖
:為案例5.1-4之溫度等值圖 :為案例5.1-1之溫度等值圖
0 1 2 3 4 5 6 7 8 9 10
0 1 2 3 4 5 6 7 8 9 10
圖 5.1-9 溫度分佈等值圖(案例 5.1-1 與案例 5.1-4)
00 1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
圖 5.1-10 熱容量分布等值圖(案例 5.1-1 與案例 5.1-4)
5. 案例 5.1-5 案例說明：
本案例探討熱量匯流項對於熱流模擬之影響。整體模擬區域如 圖 5.1-11 所示，模擬範圍、模擬網格、材質設定等與案例 5.1-1 一致。
在邊界條件設定上，左右邊界設定為 70(℃)，上下邊界則設定為無熱 通量邊界(No Heat Flux Boundary)，顯示本身並無背景熱流流動。於 點位(5, 5)處配置匯流點，並以 5,000,000(J/day)之熱通量進行熱量抽 取。因此本案例本身並無背景熱流流動，然因受到匯流項之影響，溫 度分佈將以匯流點為中心，以同心圓型式往內降溫。本案例為穩態模 擬，全體網格之初始猜值雖可任意設定，但過於懸殊之設定值易影響
式求解，各網格數值必逐步依據邊界條件之設定值，而收斂於前述之 熱流型態。
數值結果：
圖 5.1-12 為溫度之模擬等值圖，從圖面看來的確如預期一般，
形成一個橢圓形的溫降錐，其溫降錐頂的溫度改變量約為 22.77(℃)。
此外，本研究進一步作熱量平衡分析，由於系統中有個匯流點，
且上下邊界設定為無流量邊界，因此在穩態系統中，必須透過上下邊 界之來補充熱量，其補充量必須等於抽取之熱量。上下邊界之總熱通 量為 4,999,722 (J/day)，與匯流點之熱通量相比，其計算誤差約為 0.00556%，因此符合熱量平衡條件。
0 1 2 3 4 5 6 7 8 9 10
0 1 2 3 4 5 6 7 8 9 10
Sink & Source
溫度 70℃
溫度 70℃
圖 5.1-11 案例 5.1-5 格網與模式配置圖
00 1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
圖 5.1-12 溫度分佈等值圖
6. 案例 5.1-6 案例說明：
本案例修改自案例 5.1-5，為局部加密案例，因此材質區塊、邊 界條件、初始條件等之設定均一致，唯一差異處在於網格切割型式不 同，圖 5.1-13 為本案例之網格配置圖。
數值結果：
圖 5.1-14 為溫度之模擬等值圖，從圖面看來的確如預期一般，
形成一個橢圓形的溫降錐，其溫降錐頂的溫度改變量為 30.62(℃)。
與案例 5.1-5 比較，此案例在匯流點附近之溫度改變量較大，如圖 5.1-16 所示，兩案例之椎頂溫度值差異約 5(℃)。
一般而言，在數值變化較為劇烈之處，若可以較細之模擬網格 計算模擬，將會有較正確之模擬結果。然較細之模擬網格意味著需要
更多的計算節點方可佈滿整個模擬區域，因此其計算量亦隨之提高。
故透過局部加密可以在可接受的計算增加量下，提升計算精度。
以下為本案例之熱量平衡分析，上下邊界之總熱通量為
4,999,843(J/day)，與匯流點之熱通量相比，其計算誤差為 0.00314%，
因此符合熱量平衡條件。
0 1 2 3 4 5 6 7 8 9 10
0 1 2 3 4 5 6 7 8 9 10
Sink & Source
溫度 70℃
溫度 70℃
圖 5.1-13 案例 5.1-6 格網與模式配置圖
0 10 20 30 40 50 60 70
0 1 2 3 4 5 6 7 8 9
距離(公尺)
溫度(攝氏溫度)
10 案例5.1-5(規則網格)
案例5.1-6(在sink點周圍加密)
圖 5.1-16 CC'剖面溫度分布圖(案例 5.1-5 與案例 5.1-6) 7. 案例 5.1-7
案例說明：
本案例修改自案例 5.1-6，為其暫態案例。在邊界條件的設定上，
左右邊界之溫度設定為 70(℃)，上下邊界則設定為無流量邊界(No Heat Flux Boundary)。在初始條件上，所有位置之初始溫度為 70(℃)。
於點位(5, 5)處配置抽取熱量井，並以 5,000,000(J/day)之熱量抽取量 進行抽取。另外，模擬間距為 1 天，總模擬時刻數為 1500 個時刻。
數值結果：
由於此案例為案例 5.1-6 之暫態模擬，因此長期模擬後應趨於穩 態，由圖 5.1-17 可知，在趨於穩態時其模擬結果與案例 5.1-6 幾乎完 全一致，故此暫態模式通過驗證。
圖 5.1-18 為匯流點隨時刻之溫度變化圖，從圖面顯示，溫度隨 時間逐漸遞減，直到 50 天後則趨近於穩定。
:為案例5.1-7之溫度等值圖