王存海 鄭樹 張欣欣

1) (北京科技大學(xué)能源與環(huán)境工程學(xué)院, 北京 100083)

2) (冶金工業(yè)節(jié)能減排北京市重點(diǎn)實(shí)驗(yàn)室, 北京 100083)

3) (華北電力大學(xué)能源動(dòng)力與機(jī)械工程學(xué)院, 北京 102206)

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

1 引 言

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

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

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

2 控制方程及DFEM離散

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

首先分析輻射傳輸方程.參與性介質(zhì)內(nèi)輻射傳輸方程的離散坐標(biāo)形式可寫為

式中, I表示輻射強(qiáng)度, r表示坐標(biāo), Ω表示離散方向, β為衰減系數(shù), 上標(biāo)m表示離散方向.為了書寫簡潔, I (r,?m) 在以下表達(dá)中縮寫成Im.考慮介質(zhì)發(fā)射和散射, (1)式中的源項(xiàng) S (r,?m) 表示為

式中, Ib表示黑體輻射強(qiáng)度, κa表示發(fā)射系數(shù), κs表示散射系數(shù), M = Nθ× Nφ表示角度離散個(gè)數(shù),νm′表示 Ωm′方向的權(quán)重, Φm′,m表示從方向 Ωm′到方向Ωm的散射相函數(shù).

在DFEM應(yīng)用過程中, 計(jì)算區(qū)域首先被劃分為緊密相鄰的離散單元, 如圖1(a)所示.待求解的目標(biāo)量(如輻射強(qiáng)度)在相鄰單元邊界上被認(rèn)為是非連續(xù)的, 每個(gè)單元上的數(shù)值解是相互獨(dú)立的, 如圖1(b)所示.DFEM在每個(gè)離散單元上求解控制方程, 相鄰單元之間通過數(shù)值通量相互連接, 從而保證了計(jì)算區(qū)域的連續(xù)性.

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

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

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

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

本文數(shù)值通量的格式選取為迎風(fēng)格式

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

在每一個(gè)單元上對(duì)未知量采用形函數(shù)近似, 單元e上的輻射強(qiáng)度分布可表示為

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

式中Im= [I1, I2, I3]T表示離散節(jié)點(diǎn)上的輻射強(qiáng)度;矩陣Km和Hm的元素表示為

其中Nb= 3表示單元e的3個(gè)邊界.

然后分析能量擴(kuò)散方程.能量擴(kuò)散方程表示為

式中, k表示熱導(dǎo)率; ? ·q(r) 為熱輻射源項(xiàng), 表示為

與輻射強(qiáng)度類似, 單元e上溫度的形函數(shù)近似表示為

溫度的數(shù)值通量表示為

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

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

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

3 模型驗(yàn)證

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

圖2 方形介質(zhì)對(duì)稱線x/L = 0.5上無量綱溫度T/Tb分布 (a)不同空間網(wǎng)格劃分; (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.

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

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

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

圖3 方形介質(zhì)對(duì)稱線x/L = 0.5上無量綱溫度的DFEM結(jié)果和DTM結(jié)果對(duì)比 (a)不同普朗克數(shù); (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 不同數(shù)值方法得到的方形介質(zhì)對(duì)稱線上的無量綱溫度對(duì)比Fig.4.Comparison of dimensionless temperature along the square medium symmetry line x/L = 0.5 obtained by differ?ent numerical methods.

4 結(jié)果和討論

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

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

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

圖5 內(nèi)含圓形熱壁面的半圓形介質(zhì) (a)結(jié)構(gòu)示意圖; (b)網(wǎng)格劃分Fig.5.Semicircle medium with an inner circle hot boundary：(a) Geometry sketch; (b) spatial discretization.

圖6 不同數(shù)值方法得到的半圓介質(zhì)對(duì)稱線上溫度分布 (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 不同普朗克數(shù)條件下半圓介質(zhì)底邊上總熱流密度分布 (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 內(nèi)含圓形熱邊界的非規(guī)則形狀介質(zhì) (a)結(jié)構(gòu)示意圖; (b)網(wǎng)格劃分Fig.8.Irregular medium with an inner hot boundary：(a) Geometry sketch; (b) spatial discretization.

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

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

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

圖9 (a)普朗克數(shù)Npl = 0.1和1.0時(shí)內(nèi)含圓形熱邊界的非規(guī)則形狀介質(zhì)中線上溫度分布; (b) Npl = 0.1時(shí)介質(zhì)溫度分布; (c) Npl =1.0時(shí)介質(zhì)溫度分布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 內(nèi)含兩個(gè)圓形熱邊界的矩形介質(zhì) (a)結(jié)構(gòu)示意圖; (b)網(wǎng)格劃分Fig.10.The medium the square medium with two circular hot boundaries：(a) Geometry sketch; (b) spatial discr etization.

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

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

圖11 (a)普朗克數(shù)Npl = 0.1和1.0時(shí)內(nèi)含兩個(gè)圓形熱邊界的矩形介質(zhì)中線上溫度分布; (b) Npl = 0.1時(shí)介質(zhì)溫度分布; (c) Npl =1.0時(shí)介質(zhì)溫度分布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 結(jié) 論

本文將DFEM擴(kuò)展應(yīng)用于求解輻射導(dǎo)熱耦合傳熱問題, 采用非結(jié)構(gòu)化三角形單元?jiǎng)澐钟?jì)算區(qū)域, 給出了輻射傳輸方程和能量擴(kuò)散方程的DFEM離散格式, 驗(yàn)證了數(shù)值算法的正確性并研究了非規(guī)則形狀介質(zhì)內(nèi)的輻射導(dǎo)熱耦合傳熱問題.研究結(jié)果表明：本文發(fā)展的DFEM能夠得到非規(guī)則形狀介質(zhì)內(nèi)輻射導(dǎo)熱耦合傳熱問題的高精度數(shù)值結(jié)果, 尤其是對(duì)于普朗克數(shù)較小的輻射占優(yōu)問題; 非規(guī)則形狀介質(zhì)內(nèi)的溫度分布規(guī)律表明輻射效應(yīng)會(huì)弱化熱邊界附近的溫度梯度, 而導(dǎo)熱效應(yīng)則會(huì)強(qiáng)化熱邊界附近的溫度梯度.