非規則形狀介質內輻射-導熱耦合傳熱的間斷有限元求解*

王存海 鄭樹 張欣欣

1) (北京科技大學能源與環境工程學院, 北京 100083)

2) (冶金工業節能減排北京市重點實驗室, 北京 100083)

3) (華北電力大學能源動力與機械工程學院, 北京 102206)

采用間斷有限元法(discontinuous finite element method, DFEM)求解非規則形狀介質內的輻射導熱耦合傳熱問題, 得到了典型非規則形狀介質內輻射導熱耦合傳熱問題的高精度數值結果.和傳統連續型有限元方法不同, DFEM將計算區域劃分成相互獨立的離散單元, 形函數的構造、未知量的加權近似以及控制方程的求解均在每一個離散單元上進行.通過在單元之間施加迎風格式的數值通量, DFEM保證了整個計算區域的連續性, 因此這種方法兼具良好的幾何靈活性和局部守恒性.推導了輻射傳輸方程和能量擴散方程的DFEM離散格式, 驗證了DFEM求解輻射導熱耦合傳熱問題的正確性; 同時研究了不同幾何形狀介質內輻射導熱耦合傳熱問題, 得到了典型非規則形狀介質內輻射導熱耦合傳熱的高精度數值結果.

1 引 言

參與性介質內輻射導熱耦合傳熱過程廣泛存在于眾多工程實際應用中[1?3], 如高溫熔巖融化、材料加工等, 因此輻射導熱耦合問題的研究對相關傳熱分析和工程設計具有重要的意義.輻射傳輸方程是積分?微分型方程, 這類方程多通過數值解法進行求解, 如蒙特卡羅法[4?8]和修正的歐拉算法[9,10].由于輻射傳輸問題常常涉及復雜邊界條件和非規則計算區域[11,12], 國內外學者發展了多種數值方法[13?27]用以研究輻射傳熱問題.Mishra等[28,29]采用不同數值方法分別求解輻射傳輸方程和能量擴散方程, 發展了求解輻射導熱耦合問題的混合方法.

雖然目前多種數值方法被成功應用于輻射導熱耦合問題的求解, 然而非規則形狀介質內輻射導熱耦合傳熱的高精度數值求解仍舊面臨較大的挑戰[30].對于含有彎曲壁面的非規則形狀介質而言,輻射強度在界面附近等局部區域的變化率較大, 而傳統連續型數值方法則是基于全局形函數對未知量進行加權近似, 沒有考慮局部區域的大梯度輻射強度分布, 計算結果往往存在較大誤差.如文獻[31]中即使采用較為細密的網格劃分, 采用有限體積法所得輻射熱流和溫度分布曲線仍舊存在非常明顯的振蕩, 因此, 擴展高精度數值算法求解輻射導熱耦合傳熱方程, 獲得非規則形狀介質內輻射導熱耦合問題的高精度數值解具有重要意義.

間斷有限元法(discontinuous finite element method, DFEM)[32]是在傳統連續型有限元法(finite element method, FEM)基礎上發展起來的一種高精度數值方法.利用傳統FEM的空間網格劃分, DFEM利用相鄰單元邊界的數值通量從而釋放了FEM強加的內部單元連續性.由于DFEM形函數構造和控制方程求解均在每個離散單元上進行, 因此該方法兼具較好的幾何靈活性和局部守恒性, 并被逐漸應用于輻射傳輸問題的求解[33?38].本文采用DFEM求解輻射傳輸方程和能量擴散方程, 將其擴展應用于輻射導熱耦合傳熱問題的求解, 驗證了數值結果的精確性并給出了典型非規則形狀介質內輻射導熱耦合傳熱問題的高精度數值解.

2 控制方程及DFEM離散

輻射導熱耦合傳熱過程的控制方程包括輻射傳輸方程和能量擴散方程, 本文的求解思路是先求解輻射傳輸方程得到輻射強度信息, 然后再把輻射源項代入到能量擴散方程以得到溫度場分布.

首先分析輻射傳輸方程.參與性介質內輻射傳輸方程的離散坐標形式可寫為

式中, I表示輻射強度, r表示坐標, Ω表示離散方向, β為衰減系數, 上標m表示離散方向.為了書寫簡潔, I (r,?m) 在以下表達中縮寫成Im.考慮介質發射和散射, (1)式中的源項 S (r,?m) 表示為

式中, Ib表示黑體輻射強度, κa表示發射系數, κs表示散射系數, M = Nθ× Nφ表示角度離散個數,νm′表示 Ωm′方向的權重, Φm′,m表示從方向 Ωm′到方向Ωm的散射相函數.

在DFEM應用過程中, 計算區域首先被劃分為緊密相鄰的離散單元, 如圖1(a)所示.待求解的目標量(如輻射強度)在相鄰單元邊界上被認為是非連續的, 每個單元上的數值解是相互獨立的, 如圖1(b)所示.DFEM在每個離散單元上求解控制方程, 相鄰單元之間通過數值通量相互連接, 從而保證了計算區域的連續性.

圖1 (a)空間網格; (b)相鄰單元間數值通量示意圖Fig.1.(a) Spatial mesh; (b) sketch of numerical flux across the adjacent elements.

以圖1(b)中的單元e作為研究對象進行分析,在該單元上對(1)式使用權函數 w (r) 進行加權積分可得

考慮輻射強度在相鄰邊界上的非連續性, 將高斯散度定理應用于(3)式可得

式中, nΓe表示圖1(b)所示的離散單元邊界外法向單位向量, Γe表示離散單元邊界,表示相鄰邊界上的數值通量, 定義為

本文數值通量的格式選取為迎風格式

因此(4)式中的 -ImnΓe·?m可以寫為

在每一個單元上對未知量采用形函數近似, 單元e上的輻射強度分布可表示為

式中, ?i(r) 表示定義在單元e上的形函數,表示第i個節點在第m個離散方向上的輻射強度,Nn= 3表示單元e的3個節點.將(9)式代入(4)式中, 并采用形函數 ?i(r) 作為權重函數 w (r) , 可得單元e上輻射傳輸方程的DFEM離散格式為

式中Im= [I1, I2, I3]T表示離散節點上的輻射強度;矩陣Km和Hm的元素表示為

其中Nb= 3表示單元e的3個邊界.

然后分析能量擴散方程.能量擴散方程表示為

式中, k表示熱導率; ? ·q(r) 為熱輻射源項, 表示為

與輻射強度類似, 單元e上溫度的形函數近似表示為

溫度的數值通量表示為

最后可得能量擴散方程的DFEM離散格式為

式中, T = [T1, T2, T3]T包含離散節點上的溫度;矩陣M和N的元素表示為

至此, 輻射導熱耦合傳熱控制方程的DFEM離散已經完成.依次求解(10)式所示的輻射傳輸方程和(16)式所示的能量擴散方程即可得到離散節點的溫度值.

3 模型驗證

基于上述DFEM離散格式, 發展了MATLAB數值計算程序用以求解(11)和(17)式中的各項元素, 進而通過求解方程(10)和(16)得到輻射強度和溫度分布.本節將DFEM應用于求解二維方形介質內的輻射導熱耦合傳熱以驗證數值模型的正確性.方形介質邊長為L, 衰減系數為β, 光學厚度為 τ = βL = 1.0.所有壁面均為發射率 ε = 1.0的黑壁面, 底邊溫度為Tb= 1000 K, 其他壁面溫度為 500 K, 普朗克數為 Npl= kβ/(4σTb3), 其中k表示熱導率, σ 表示 Stefan?Boltzmann 常數, σ =5.67 × 10—8.

圖2 方形介質對稱線x/L = 0.5上無量綱溫度T/Tb分布 (a)不同空間網格劃分; (b)不同角度劃分Fig.2.Dimensionless temperature T/Tb along the sym?metry line x/L = 0.5 for the cases：(a) Different spatial dis?cretization schemes; (b) different angular discretization schemes.

首先進行網格無關性檢驗.分別將方形區域的邊界均勻劃分為5, 10和15個單元, 二維計算區域離散為Ne= 62, 226和504個三角形單元.圖2(a)所示為角度劃分為 Nθ× Nφ= 10 × 20, 不同離散單元條件下介質對稱線x/L = 0.5上的無量綱溫度T/Tb的分布規律.圖2(b)所示為離散單元數Ne= 226, 不同角度劃分 Nθ× Nφ= 3 × 6, 5 ×10, 10 × 20和15 × 30條件下對稱線上的無量綱溫度分布.由圖2所示的結果可知, 離散單元Ne=226, 離散角度 Nθ× Nφ= 10 × 20 條件下所得的結果滿足網格無關性要求, 在該網格劃分條件下采用DFEM求解本算例所需時間為25.38 s.

針對不同普朗克數Npl, DFEM所得無量綱溫度分布如圖3(a)所示, 并和文獻[39]采用離散傳遞法(discrete transfer method, DTM)所得結果進行了比較.針對不同的介質散射反照率ω = κa/β,DFEM和DTM所得無量綱溫度的對比結果如圖3(b)所示.由圖3(b)結果可知DFEM結果和文獻DTM結果符合得很好, 驗證了DFEM求解輻射導熱耦合問題的正確性.

圖4所示為在Npl= 0.01, β = 1.0以及ω = 0.0條件下, 不同數值方法所得對稱線x/L = 0.5上無量綱溫度分布的對比結果.DTM是一種基于光線跟蹤的數值方法[39], 避免了控制方程近似處理引起的誤差, 因此采用該方法所得結果可視為該問題的基準解.以DTM結果作為基準, 采用有限體積法(finite volume method, FVM)[29]所得結果最大誤差為2.21%, 而本文DFEM結果的最大誤差僅為0.68%, 說明DFEM結果更加精確.

圖3 方形介質對稱線x/L = 0.5上無量綱溫度的DFEM結果和DTM結果對比 (a)不同普朗克數; (b)不同散射反照率Fig.3.Comparison of dimensionless temperature along the square medium symmetry line x/L = 0.5 obtained by DFEM and DTM for the cases：(a) Different Planck num?bers; (b) different scattering albedos.

圖4 不同數值方法得到的方形介質對稱線上的無量綱溫度對比Fig.4.Comparison of dimensionless temperature along the square medium symmetry line x/L = 0.5 obtained by differ?ent numerical methods.

4 結果和討論

首先考慮如圖5(a)所示的內含圓形壁面(半徑Rc= 0.2 m)的半圓形(半徑Rs= 1.0 m)介質內的輻射導熱耦合問題.介質衰減系數β =1.0 m—1, 內部圓形壁面溫度Tc= 400 K, 外部半圓形壁面和底邊溫度Ts= 300 K, 壁面均為黑壁面且發射率ε = 1.0.采用DFEM方法求解沿著半徑方向的無量綱溫度分布.方向離散為Nθ× Nφ=20 × 40, 空間離散為如圖5(b)所示的1768個三角形單元, 所得結果滿足網格無關性要求, 采用DFEM求解本算例所需時間為182.69 s.

在普朗克數Npl= kβ/(4σTs3) = 0.1條件下,圖6所示為采用不同方法(Fluent[31], 不同網格處理方式的FVM[31], 混合格子Boltzmann?有限體積法(LBM?FVM)[40]以及本文DFEM)所得的半圓形介質對稱線x = 0上的溫度分布.由圖6可知采用Fluent和LBM?FVM以及本文DFEM所得結果符合得較好; 而文獻[31]采用FVM所得結果則具有較大的偏差, 這是文獻[21]采用四邊形網格處理彎曲壁面引起的計算誤差所致.圖6所示結果表明, 和FVM相比, 本文發展的DFEM求解非規則形狀介質內的輻射導熱耦合問題更加精確.

圖7所示為普朗克數Npl= 0.1和1.0條件下半圓形介質底邊y = 0上的總熱流密度(負號表示方向)的分布規律.由圖7可知底邊中點位置處的熱流密度最強, 這是由于內部圓環溫度較高,底邊中點距離高溫圓環最近, 在輻射和導熱共同作用條件下, 該位置處的換熱強度達到最大.由圖7(a)結果可知, 當Npl= 0.1時, DFEM結果和Fluent結果吻合良好, 而LBM?FVM結果存在明顯誤差.當普朗克數增大到Npl= 1.0時, 采用不同方法所得結果符合得比較好, 如圖7(b)所示.圖7所示結果表明, 對于普朗克數較小即輻射占優的輻射導熱耦合傳熱問題, DFEM所得結果更加精確.

圖5 內含圓形熱壁面的半圓形介質 (a)結構示意圖; (b)網格劃分Fig.5.Semicircle medium with an inner circle hot boundary：(a) Geometry sketch; (b) spatial discretization.

圖6 不同數值方法得到的半圓介質對稱線上溫度分布 (a) y = [0.0, 0.2]; (b) y = [0.6, 1.0]Fig.6.Temperature distributions along the symmetric line of the semicircle medium obtained by different numerical algorithms：(a) y = [0.0, 0.2]; (b) y = [0.6, 1.0].

圖7 不同普朗克數條件下半圓介質底邊上總熱流密度分布 (a) Npl = 0.1, (b) Npl = 1.0Fig.7.Total flux distributions along the bottom boundary of the semicircle medium under the situations with different Plank num?bers：(a) Npl = 0.1; (b) Npl = 1.0.

圖8 內含圓形熱邊界的非規則形狀介質 (a)結構示意圖; (b)網格劃分Fig.8.Irregular medium with an inner hot boundary：(a) Geometry sketch; (b) spatial discretization.

接下來考慮如圖8(a)所示的內含圓形熱邊界的非規則形狀介質, 該模型可用來研究圓形熱管和環境之間的輻射導熱耦合換熱.介質尺寸(單位為m)在圖8(a)中給出, 介質衰減系數為β = 1.0 m—1,內部圓形邊界溫度為400 K, 其他邊界溫度為自然環境溫度300 K, 所有壁面均為黑壁面且發射率為ε = 1.0.計算區域離散為圖8(b)所示的1142個三角形單元、方向離散為 Nθ× Nφ= 20 × 40條件下所得結果達到網格無關性要求, 采用DFEM求解本算例所需時間為118.38 s.

圖9(a)所示為不同普朗克數Npl= 0.1和 1.0條件下非規則介質中線x = 0.5上的溫度分布.由圖9(a)可知, 對于Npl= 0.1的輻射占優問題, 從壁面到圓環之間的溫度梯度先減小后增加, 說明靠近圓環和彎曲壁面處的介質溫度變化較為劇烈, 而中間區域介質的溫度變化較為平緩; 對于Npl= 1.0的導熱占優問題, 該區域內的溫度梯度逐漸增加,說明彎曲壁面附近處的介質溫度變化較為平緩, 高溫圓環附近處的介質溫度變化劇烈; 圓環上方的介質溫度變化劇烈程度也有所區別, 但由于該區間垂直高度較小, 因此該區間范圍的介質溫度沿y方向接近于線性減小.

在中心線上的同一位置處, 普朗克數Npl=0.1對應的介質溫度總是高于Npl= 1.0對應的介質溫度.這是由于當Npl= 0.1時, 輻射傳熱作用顯著, 相同位置的輻射源項較強, 因此介質溫度更高.圖9(b)和圖9(c)所示分別為同普朗克數Npl=0.1和 1.0條件下的介質溫度分布的等值線圖, 圖9中相鄰兩條等值線之間的溫度差為10 K.由圖9(b)和圖9(c)的對比結果可知, 當普朗克數Npl= 0.1時, 強烈的輻射效應能夠使高溫圓環能量傳播得更遠, 因此介質內的能量分布更加均勻; 當普朗克數Npl= 1.0時, 顯著的導熱效應導致溫度從高溫圓環向外迅速降低, 因此外部邊界附近較大范圍內的介質溫度處于較低水平.

圖9 (a)普朗克數Npl = 0.1和1.0時內含圓形熱邊界的非規則形狀介質中線上溫度分布; (b) Npl = 0.1時介質溫度分布; (c) Npl =1.0時介質溫度分布Fig.9.(a) Temperature distributions along the centerline of the irregular medium with an inner hot boundary; (b) temperature dis?tribution within the computation domain for the case of Npl = 0.1; (c) temperature distribution within the computation domain for the case of Npl = 1.0.

圖10 內含兩個圓形熱邊界的矩形介質 (a)結構示意圖; (b)網格劃分Fig.10.The medium the square medium with two circular hot boundaries：(a) Geometry sketch; (b) spatial discr etization.

進一步考慮如圖10(a)所示的內含兩個高溫圓環的矩形介質, 該模型可用來研究兩個高溫管道及其所處環境之間的輻射導熱耦合換熱.矩形介質邊長為 Lx× Ly= 1.0 m × 1.0 m, 內部圓環半徑為R = 0.1 m, 圓心位置分別為(x, y) = (0.3 m,0.3 m)和(x, y) = (0.7 m, 0.7 m).矩形介質邊界和圓環邊界之間充滿衰減系數為β = 1.0 m—1的參與性介質, 內部兩個圓形邊界溫度均為400 K, 其他邊界溫度為300 K, 所有壁面均為黑壁面且發射率為ε = 1.0.計算區域離散為圖10(b)所示的1262 個三角形單元、方向離散為 Nθ× Nφ= 20 ×40條件下所得到的結果達到了網格無關性的要求,采用DFEM求解本算例所需時間為131.75 s.

在不同普朗克數Npl= 0.1和 1.0條件下, 中線x = 0.5上的溫度分布如圖11(a)所示, 可知介質中線溫度在y = [0, 0.3]逐漸升高且溫度變化劇烈程度逐漸降低, 再y = [0.3, 0.7]的介質保持較高的溫度, 在y = [0.7, 1.0]逐漸降低.在同一位置處, Npl= 0.1對應的介質溫度高于Npl= 1.0對應的介質溫度, 說明輻射效應作用對該模型中線上的介質溫度具有顯著的提升作用.圖11(b)和圖11(c)所示分別為同普朗克數Npl= 0.1和 1.0條件下計算區域內溫度分布的等值線圖, 相鄰兩條等值線之間的溫度差為5 K.由圖11(b)和圖11(c)所示結果可知, 介質中心點位置(x, y) = (0.5 m, 0.5 m)周圍存在一個溫度梯度較小的區域, 且普朗克數越小, 該區域的面積越大.以上結果表明輻射效應會弱化中心點位置處的溫度梯度, 而導熱效應則會強化該位置處的溫度梯度.

圖11 (a)普朗克數Npl = 0.1和1.0時內含兩個圓形熱邊界的矩形介質中線上溫度分布; (b) Npl = 0.1時介質溫度分布; (c) Npl =1.0時介質溫度分布Fig.11.(a) Temperature distributions along the centerline of the square medium with two circular hot boundaries; (b) temperature distribution within the computation domain for the case of Npl = 0.1; (c) temperature distribution within the computation domain for the case of Npl = 1.0.

5 結 論

本文將DFEM擴展應用于求解輻射導熱耦合傳熱問題, 采用非結構化三角形單元劃分計算區域, 給出了輻射傳輸方程和能量擴散方程的DFEM離散格式, 驗證了數值算法的正確性并研究了非規則形狀介質內的輻射導熱耦合傳熱問題.研究結果表明：本文發展的DFEM能夠得到非規則形狀介質內輻射導熱耦合傳熱問題的高精度數值結果, 尤其是對于普朗克數較小的輻射占優問題; 非規則形狀介質內的溫度分布規律表明輻射效應會弱化熱邊界附近的溫度梯度, 而導熱效應則會強化熱邊界附近的溫度梯度.