鐵礦石燒結(jié)過(guò)程傳熱傳質(zhì)數(shù)值模擬

張小輝，張家元，張建智，蘇浩，周孑民

(中南大學(xué) 能源科學(xué)與工程學(xué)院，湖南 長(zhǎng)沙，410083)

將鐵礦粉與石灰石、燃料等混合后進(jìn)行燒結(jié)造塊的過(guò)程中，會(huì)發(fā)生一系列復(fù)雜的物理化學(xué)反應(yīng)，這些物理化學(xué)反應(yīng)均伴隨著傳熱和傳質(zhì)[1-3]，在這一過(guò)程中，燒結(jié)礦的溫度與這些反應(yīng)互相影響。目前，對(duì)于燒結(jié)傳熱過(guò)程的研究主要集中在對(duì)燒結(jié)料層溫度的準(zhǔn)確預(yù)測(cè)以及燒結(jié)終點(diǎn)的判斷[4-5]，而計(jì)算過(guò)程多是進(jìn)行相應(yīng)的簡(jiǎn)化。夏德宏等[6]對(duì)燒結(jié)過(guò)程碳燃燒與傳熱過(guò)程建立數(shù)學(xué)模型，進(jìn)行數(shù)值模擬，并開(kāi)發(fā)出可視化軟件。但是，該模型只考慮碳的燃燒。龍紅明等[7]對(duì)燒結(jié)過(guò)程中各層分別建立傳熱方程即分段計(jì)算模型，對(duì)燒結(jié)料層溫度進(jìn)行模擬計(jì)算，但實(shí)際的燒結(jié)過(guò)程是連續(xù)的，該模型也沒(méi)有考慮點(diǎn)火、保溫階段以及燃燒帶可能發(fā)生碳的不完全燃燒反應(yīng)，而這種反應(yīng)是煙氣中含有CO的主要原因。Yang等[8-9]建立一維非穩(wěn)態(tài)的計(jì)算模型，考慮氣固之間的對(duì)流、傳熱以及輻射，子模型包含氣固反應(yīng)、氣體反應(yīng)以及孔隙率的變化。作者利用這一模型對(duì)不同碳質(zhì)量分?jǐn)?shù)、不同抽風(fēng)量的工況進(jìn)行優(yōu)化模擬。Tan等[10]針對(duì)燒結(jié)過(guò)程建立包含有15個(gè)物理化學(xué)反應(yīng)的計(jì)算流體力學(xué)計(jì)算模型，用以對(duì)燒結(jié)過(guò)程二惡英的產(chǎn)生進(jìn)行研究，但是，該研究沒(méi)有考慮石灰石的分解對(duì)料層溫度及廢氣成分的影響。為此，本文作者以局部非熱力學(xué)平衡理論和組分傳輸理論為基礎(chǔ)，結(jié)合燒結(jié)過(guò)程中各反應(yīng)的動(dòng)力學(xué)方程，全面考慮燒結(jié)過(guò)程中的反應(yīng)以及這些反應(yīng)對(duì)于傳熱和傳質(zhì)的影響，建立三維非穩(wěn)態(tài)計(jì)算模型，對(duì)燒結(jié)傳熱和傳質(zhì)過(guò)程進(jìn)行研究。

1 模型的建立

對(duì)于氣體流過(guò)料層時(shí)所產(chǎn)生的壓力降，采用多孔介質(zhì)模型來(lái)計(jì)算[11-12]。對(duì)于多孔介質(zhì)條件下所產(chǎn)生的動(dòng)量方程附加源項(xiàng)采用Ergun公式進(jìn)行修正[12-13]。

1.1 傳熱模型

對(duì)于計(jì)算單元內(nèi)氣、固溫度的求解采用局部非熱力學(xué)平衡雙能量方程[13-14]模型，分別建立能量方程，通用方程組[15-16]如式(1)和(2)所示。

固相：

氣相：

其中：hp為氣、固傳熱系數(shù)，W/(m2·K)；dp為固體粒徑，m。

1.2 組分傳輸模型

組分輸運(yùn)方程[18]如下：

1.3 燒結(jié)過(guò)程反應(yīng)動(dòng)力學(xué)模型

燒結(jié)過(guò)程所發(fā)生的物理化學(xué)反應(yīng)主要有物料的干燥、碳燃燒、石灰石分解、鐵礦石的融化及凝固[11]等。

1.3.1 物料干燥

對(duì)于物料的干燥(水分蒸發(fā)，即(H2O)l=(H2O)g)過(guò)程，考慮2個(gè)階段：恒速干燥階段和減速干燥階段。這2個(gè)階段以物料的含水量是否達(dá)到臨界含水量為標(biāo)準(zhǔn)[17]。

當(dāng)料層含水量大于臨界值時(shí)，處于恒速干燥階段。此時(shí)由上層氣體傳遞的熱量完全用于水的蒸發(fā)，認(rèn)為T(mén)s為一定值。這一階段的干燥速率為：

其中：a為顆粒與床層的比表面積，m-1；ΔH1為水的汽化潛熱，J/kg；Mw為水相對(duì)分子質(zhì)量。

當(dāng)料層含水量小于臨界值大于平衡值時(shí)，處于降速干燥階段，此時(shí)由上層氣體傳遞的熱量用于固體顆粒溫度的升高和水的汽化潛熱。在燒結(jié)高溫條件下，達(dá)到平衡含水量時(shí)已非常接近干燥狀態(tài)，因此，將這一階段適當(dāng)簡(jiǎn)化，視平衡含水量為 0，這一階段的干燥速率為：

其中：w和 wc分別為含水質(zhì)量分?jǐn)?shù)、臨界含水質(zhì)量分?jǐn)?shù)。

1.3.2 碳反應(yīng)

在燒結(jié)過(guò)程中，碳的反應(yīng)[10]有：C+O2=CO2，C+CO2=CO，C+H2O=CO+H2。

對(duì)于C+O2=CO2，反應(yīng)速率[19]為：

其中：kf,1為傳質(zhì)系數(shù)，m/s；k1為反應(yīng)速率常數(shù)，m/s。對(duì)于C+CO2=CO，反應(yīng)速率[19]為：

其中：φc為燃料顆粒形狀系數(shù)；P為料層內(nèi)壓力；y為煙氣中CO2體積分?jǐn)?shù)；Ef,2為反應(yīng)有效系數(shù)；kf,2為界膜傳質(zhì)系數(shù)，m/s；k2為反應(yīng)速率常數(shù)，m3/(kg·s)。

對(duì)于C+H2O=CO+ H2，反應(yīng)速率[19]為

其中：u為煙氣中水蒸氣體積分?jǐn)?shù)；kf,3為界膜傳質(zhì)系數(shù)，m/s；Ef,3為反應(yīng)有效系數(shù)；k3為反應(yīng)速率常數(shù)，m3/(kg·s)。

1.3.3 石灰石分解

石灰石的分解反應(yīng)不受內(nèi)擴(kuò)散限制，燒結(jié)料層中氣流速度快，外擴(kuò)散阻力較小，可以認(rèn)為反應(yīng)是處于化學(xué)反應(yīng)的控制范圍之內(nèi)，反應(yīng)速度主要受溫度的影響。根據(jù)熱平衡，石灰石分解反應(yīng)速率[3]為：

其中：al為石灰石顆粒比表面積，m-1；Tl為分解開(kāi)始溫度，K；Ul為與分解度相關(guān)的系數(shù)。

1.3.4 鐵氧化物的反應(yīng)

燒結(jié)過(guò)程中主要發(fā)生鐵的氧化物的還原反應(yīng)，鐵礦石主要成分不同，其所發(fā)生的反應(yīng)也不同。這里以赤鐵礦為例。根據(jù)文獻(xiàn)[19]，主要發(fā)生的反應(yīng)為：Fe2O3+3CO=2Fe+3CO2，F(xiàn)e2O3+3H2=2Fe+3H2O。

對(duì)于Fe2O3+3CO=2Fe+3CO2，反應(yīng)速率為：

其中：kf,5為界膜傳質(zhì)系數(shù)，m/s；x為煙氣中CO體積分?jǐn)?shù)；fs為鐵礦石還原度；k5為反應(yīng)速率常數(shù)，m/s；K5為反應(yīng)平衡常數(shù)。

對(duì)于Fe2O3+3H2=2Fe+3H2O，反應(yīng)速率為：

其中：kf,6為界膜傳質(zhì)系數(shù)，m/s； DH2為H2分子擴(kuò)散系數(shù)，m2/s；k6為反應(yīng)速率常數(shù)，m/s；K6為反應(yīng)平衡常數(shù)。

1.3.5 液相形成與凝固

根據(jù)文獻(xiàn)[20]，燒結(jié)過(guò)程中液相形成速率為：

凝固速率為：

其中：Tmlt為開(kāi)始熔化(凝固)溫度，K。

1.3.6 其他氣相反應(yīng)

氣相中含有少量生成的H2，在高溫區(qū)與O2反應(yīng)：H2+O2=(H2O)g，反應(yīng)速率[20]為：

其中：φ為煙氣中H2體積分?jǐn)?shù)。

在高溫區(qū)還會(huì)發(fā)生水煤氣轉(zhuǎn)化反應(yīng) CO+H2O=CO2+H2，反應(yīng)速率[19]為：

其中：

其中：R為氣體常數(shù)，J/(mol·K)。

2 模型的求解

以上模型基于FLUNT計(jì)算平臺(tái)建立和求解。對(duì)氣體流過(guò)燒結(jié)礦時(shí)由黏性阻力和慣性阻力引起的附加動(dòng)量源項(xiàng)[21]，利用Ergun公式計(jì)算黏性阻力系數(shù)和慣性阻力系數(shù)[13]；對(duì)于基于局部非熱力學(xué)平衡雙能量方程所建立的氣、固兩相能量方程，首先利用自定義標(biāo)量(UDS)分別定義Ts和Tf，然后結(jié)合以上反應(yīng)的動(dòng)力學(xué)方程，利用自定義功能(UDF)功能對(duì)能量方程的源項(xiàng)分別在不同的溫度條件下進(jìn)行定義；對(duì)于燒結(jié)廢氣成分的模擬基于組分輸運(yùn)和有限速率化學(xué)反應(yīng)模型，基于 UDF功能同時(shí)結(jié)合以上反應(yīng)的反應(yīng)速率分別進(jìn)行定義。

用于該計(jì)算模型的參數(shù)及初始值如表1所示，部分參數(shù)的取值參考漣鋼360 m2燒結(jié)機(jī)實(shí)際操作參數(shù)。

表1 計(jì)算參數(shù)初始值Table 1 Initial value of simulation parameters

由于燒結(jié)的傳熱傳質(zhì)是非穩(wěn)態(tài)過(guò)程，因此，在求解過(guò)程中，根據(jù)所處階段的不同調(diào)整邊界條件，在點(diǎn)火、保溫、燒結(jié)階段，調(diào)整入口風(fēng)量、風(fēng)溫；對(duì)于傳熱計(jì)算，碳顆粒燃燒的溫度條件為923 K，石灰石分解的溫度條件為993 K。鐵礦石開(kāi)始熔化及凝固溫度為1 473 K。其他在本文中所用到的參數(shù)及計(jì)算式如表2所示。

表2 計(jì)算模型所用到的參數(shù)Table 2 Parameters in simulation model

3 計(jì)算結(jié)果及分析

利用以上所建立的模型進(jìn)行鐵礦石燒結(jié)傳熱傳質(zhì)的模擬計(jì)算，圖1所示為不同料層高度處燒結(jié)料層與氣體溫度變化情況。料層高度以料面為起點(diǎn)計(jì)算，由圖1可知：在氣體溫度達(dá)到最高點(diǎn)之前，同一點(diǎn)處的氣體溫度大于固體溫度，此時(shí)氣體向固體傳熱，加熱料層并使燃料燃燒；隨著料層中燃料的快速燃燒，固體溫度超過(guò)氣體溫度，此時(shí)固體向氣體傳熱，隨著燃燒的完成，由上而下的氣流對(duì)料層進(jìn)行冷卻；由于蓄熱作用，燒結(jié)速度逐漸加快，并且料層最高溫度逐漸增加。

圖2所示為燒結(jié)料層在整個(gè)燒結(jié)過(guò)程中的溫度分布情況。由圖2可見(jiàn)：料層最高溫度在前700 s偏低，在700 s之后由于料層的蓄熱作用，最高溫度逐漸升高，同時(shí)高溫帶的厚度不斷增加。

圖1 不同高度處氣、固溫度Fig.1 Gas and solid temperatures in different heights

圖2 燒結(jié)料層溫度分布Fig.2 Temperature distribution of sinter

圖3所示是從點(diǎn)火到燒結(jié)完成整個(gè)過(guò)程中料層最高溫度以及計(jì)算域出口處的氣體溫度的變化趨勢(shì)。固體溫度的測(cè)點(diǎn)沿臺(tái)車高度等距離分布， 即在0.6 m料層高度上平均分布6個(gè)測(cè)點(diǎn)，從點(diǎn)火開(kāi)始每60 s取1次值，測(cè)量3次，將每1點(diǎn)處測(cè)得的固體溫度值進(jìn)行平均，然后取6個(gè)點(diǎn)的最大值。氣體溫度的測(cè)點(diǎn)設(shè)在位于臺(tái)車下部的風(fēng)箱入口處，每個(gè)風(fēng)箱設(shè)1測(cè)點(diǎn)，取5次測(cè)量平均值作為每1風(fēng)箱的氣體平均溫度。

從圖3可見(jiàn)：保溫段結(jié)束之后固體料層最高溫度下降約80 K，之后又穩(wěn)定上升，時(shí)間從240 s持續(xù)到400 s左右，可稱為“過(guò)渡階段”。這一固體最高溫度降低的階段發(fā)生在離料面約0.04～0.1 m的高度范圍內(nèi)，這一階段最高溫度的突降，是由進(jìn)入料層氣流的溫度突降和流量增加所引起。在整個(gè)燒結(jié)過(guò)程中料層溫度最高點(diǎn)位于倒數(shù)第2個(gè)風(fēng)箱位置左右，此后料層溫度下降，因此，可以認(rèn)為此時(shí)整個(gè)料層的燒結(jié)反應(yīng)結(jié)束，料層溫度達(dá)到最高點(diǎn)時(shí)燒結(jié)臺(tái)車所處的位置可認(rèn)為是燒結(jié)終點(diǎn)位置。

圖3 燒結(jié)過(guò)程中氣、固溫度變化曲線Fig.3 Temperature curves of solid and gas in sintering process

由圖3可見(jiàn)：對(duì)于臺(tái)車底部廢氣溫度的變化，在1 200 s之前，出口處的氣體溫度上升緩慢，在燒結(jié)反應(yīng)接近料層底部的時(shí)候，氣體溫度迅速上升，氣體溫度達(dá)到最高點(diǎn)的位置與料層溫度達(dá)到最高點(diǎn)的位置接近，因此，燒結(jié)臺(tái)車廢氣溫度達(dá)到最高點(diǎn)的位置也可以作為燒結(jié)終點(diǎn)位置的判斷標(biāo)準(zhǔn)。由于風(fēng)箱中氣流溫度是幾個(gè)臺(tái)車底部氣流的混合，故測(cè)量值與達(dá)到燒結(jié)終點(diǎn)時(shí)臺(tái)車底部氣流溫度相差較大，但是，溫度發(fā)展趨勢(shì)與計(jì)算結(jié)果一致。

圖4所示為料層底部出口處氣體成分體積分?jǐn)?shù)在燒結(jié)過(guò)程中的變化以及與實(shí)測(cè)值的比較。測(cè)點(diǎn)分布于臺(tái)車下部的各風(fēng)箱入口處，與廢氣溫度的測(cè)量共用測(cè)點(diǎn)。由水蒸氣含量的變化可知在點(diǎn)火啟動(dòng)之后，出口處水蒸氣體積分?jǐn)?shù)快速上升，之后在保溫階段完成時(shí)達(dá)到平穩(wěn)狀態(tài)；隨著氣體及燒結(jié)料層溫度的不斷上升，料層含水量不斷減少，出口氣體中水蒸氣體積分?jǐn)?shù)也逐漸減小在1 000 s時(shí)開(kāi)始減小，在1 200 s左右料層水分蒸發(fā)完畢。由O2，CO和CO2體積分?jǐn)?shù)的變化可知，在點(diǎn)火和保溫階段，從料面到離料面0.04 m的高度內(nèi)溫度高，在短時(shí)間內(nèi)快速發(fā)生上述一系列物理化學(xué)反應(yīng)，而該階段內(nèi)抽風(fēng)負(fù)壓小，氣體流量小，故前240 s這3種氣體體積分?jǐn)?shù)變化較大，O2快速減少，而CO和CO2快速增加；在保溫階段后，氣體成分趨于穩(wěn)定，在燒結(jié)前沿到達(dá)離料面約0.55 m處時(shí)，O2，CO和CO2的體積分?jǐn)?shù)又逐漸恢復(fù)到與入口處的濃度達(dá)到一致的水平。

圖3和圖4中計(jì)算值與測(cè)試值的對(duì)比可知，計(jì)算值與測(cè)試值吻合較好，說(shuō)明該模型在預(yù)測(cè)燒結(jié)過(guò)程的氣、固溫度變化以及燒結(jié)廢氣成分上是可靠的。

圖4 出口處氣體組分體積分?jǐn)?shù)Fig.4 Volume concentrations of gas composition in outflow

4 結(jié)論

(1) 通過(guò)對(duì)鐵礦石燒結(jié)過(guò)程傳熱傳質(zhì)的研究，結(jié)合多孔介質(zhì)模型、k-ε湍流模型、基于局部非熱力學(xué)平衡的雙方程能量模型、組分輸運(yùn)模型以及化學(xué)反應(yīng)動(dòng)力學(xué)模型，以 FLUENT軟件為計(jì)算平臺(tái)，利用 UDF和UDS功能，建立了燒結(jié)過(guò)程傳熱傳質(zhì)的數(shù)學(xué)計(jì)算模型。

(2) 利用所建立的模型對(duì)燒結(jié)過(guò)程氣、固溫度以及廢氣成分進(jìn)行計(jì)算和分析，通過(guò)與測(cè)試值的比較，證明了該模型的可靠性。研究結(jié)果表明：在點(diǎn)火到正常抽風(fēng)燒結(jié)的“過(guò)渡階段”， 在離料面0.04～0.1 m的高度范圍內(nèi)，燒結(jié)礦最高溫度降低約80 K，燒結(jié)速度減慢，導(dǎo)致出口處氣體成分發(fā)生相應(yīng)突變。

(3) 利用該模型可以模擬計(jì)算在不同配炭比、不同點(diǎn)火及燒結(jié)條件下的傳熱傳質(zhì)過(guò)程，并對(duì)燒結(jié)過(guò)程的控制、廢氣成分的控制以及燒結(jié)終點(diǎn)的預(yù)測(cè)有指導(dǎo)意義。

[1] Yang W, Choi S, Choi E S, et al. Combustion characteristics in an iron ore sintering bed: Evaluation of fuel substitution[J].Combustion and flame, 2006, 145(3): 447-463.

[2] Sadrnezhaad S K, Ferdowsi A, Payab H. Mathematical model for a straight grate iron ore pellet induration process of industrial scale[J]. Computational Materials Science, 2008, 44(2):296-302.

[3] 張玉柱, 艾立群. 鋼鐵冶金過(guò)程的數(shù)學(xué)解析與模擬[M]. 北京:冶金工業(yè)出版社, 1997: 95-107 ZHANG Yuzhu, AI Liqun, Mathematical analysis and simulation of metallurgical process[M]. Beijing: Metallurgical Industry Press, 1997: 95-107.

[4] Mitterlehner J, Loeffler G, Winter F, et al. Modeling and simulation of heat front propagation in the iron ore sintering process[J]. ISIJ International, 2004, 44(1): 11-20.

[5] 張家元, 李黎, 張小輝. 基于數(shù)值模擬的燒結(jié)終點(diǎn)在線控制[J]. 燒結(jié)球團(tuán), 2011, 36(3): 5-9.ZHANG Jiayuan, LI Li, ZHANG Xiaohui. On-line control of burning through point based on numerical simulation[J].Sintering and Pelletizing, 2011, 36(3): 5-9.

[6] 夏德宏, 焦紅蕾, 張剛, 等. 燒結(jié)工序中燃燒與傳熱過(guò)程的數(shù)值模擬[J]. 工業(yè)爐, 2006, 28(3): 30-32.XIA Dehong, JIAO Honglei, ZHANG Gang, et al. Numerical simulation of combustion and heat transfer in sintering process[J]. Industrial Furnace, 2006, 28(3): 30-32.

[7] 龍紅明, 范曉慧, 毛曉明, 等. 基于傳熱的燒結(jié)料層溫度分布模型[J]. 中南大學(xué)學(xué)報(bào): 自然科學(xué)版, 2008, 39(3): 436-442.LONG Hongming, FAN Xiaohui, MAO Xiaoming, et al.Sintering bed temperature distribution model based on heat transfer[J]. Journal of Central South University: Science and Technology, 2008, 39(3): 436-442.

[8] Yang W, Ryu C, Choi S, et al. Mathematical model of thermal processes in an iron ore sintering bed[J]. Metals and Materials International, 2004, 10(5): 493-500.

[9] Yang W, Ryu C, Choi S, et al. Modeling of combustion and heat transfer in an iron ore sintering bed with considerations of multiple solid phases[J]. ISIJ International, 2004, 44(3):492-499.

[10] Tan P, Neuschütz D. CFD modeling of sintering phenomena during iron ore sintering[C]//Multiphase Phenomena and CFD Modeling and Simulation in Materials Processes. Charlotte:TMS, 2004: 451-459.

[11] Li H G, Zhang J L, Pei Y D, et,al. Melting characteristics of iron ore fine during sintering process[J]. Journal of Iron and Steel Research, 2011, 18(5): 11-15.

[12] 張家元, 張小輝, 周孑民, 等. 燒結(jié)礦冷卻過(guò)程正交模擬優(yōu)化試驗(yàn)研究[J]. 鋼鐵, 2011, 46(7): 86-89.ZHANG Jiayuan, ZHANG Xiaohui, ZHOU Jiemin, et al.Optimal orthogonal simulation research of sinter cooling[J]. Iron and Steel, 2011, 46(7): 86-89.

[13] 張小輝, 張家元, 戴傳德, 等. 燒結(jié)礦冷卻過(guò)程數(shù)值仿真與優(yōu)化[J]. 化工學(xué)報(bào), 2011, 62(11): 3081-3087.ZHANG Xiaohui, ZHANG Jiayuan, DAI Chuande, et al.Optimal and simulation of sinter cooling process[J]. CIESC Journal, 2011, 62(11): 3081-3087.

[14] 劉偉, 范愛(ài)武, 黃曉明. 多孔介質(zhì)傳熱傳質(zhì)理論與應(yīng)用[M].北京: 科學(xué)出版社, 2006: 28-32.LIU Wei, FAN Aiwu, HUANG Xiaoming. Porous media theory and application of heat and mass transfer[M]. Beijing: Science Press, 2006: 28-32.

[15] Alazmi B, Vafai K. Analysis of fluid flow and heat transfer interfacial conditionsbetween a porous medium and a fluid layer[J]. International Journal of Heat and Mass Transfer, 2001,44(9): 1735-1749.

[16] 李菊香, 涂善東. 考慮局部非熱平衡的流體層流橫掠多孔介質(zhì)中恒熱流平板的傳熱分析[J]. 化工學(xué)報(bào), 2010, 61(1): 10-14.LI Juxiang, TU Shandong. Heat transfer of laminar flow over a plate embedded in porous medium with a constant heat flux under local non-equilibrium condition[J]. CIESC Journal, 2010,61(1): 10-14.

[17] Ramos M V, Kasai E, Kano J, et al. Numerical simulation model of the iron ore sintering process directly describing the agglomeration phenomenon of granules in the packed bed[J].ISIJ International, 2000, 40(5): 448-454.

[18] 周萍, 周乃君, 蔣愛(ài)華, 等. 傳遞過(guò)程原理及其數(shù)值仿真[M].長(zhǎng)沙: 中南大學(xué)出版社, 2006: 300.ZHOU Ping, ZHOU Naijun, JIANG Aihua, et, al. Principle and numerical simulation of transport processes[M]. Changsha:Central South University Press, 2006: 300.

[19] 鞭巖, 森山昭. 冶金反應(yīng)工程學(xué)[M]. 蔡志鵬, 謝裕生, 譯. 北京: 科學(xué)出版社, 1981: 282-288.Muchi I, Moriyama A. Reaction engineering in metallurgical[M].CAI Zhipeng, XIE Yusheng, trans. Beijing: Science Press, 1981:282-288.

[20] Tan P, Neuschütz D. Study on polychlorinated dibenzo-p-dioxin/furan formation in iron ore sintering process[J].Metallurgical and Materials Transactions B, 2004, 35(5):983-991.

[21] Carvalho M G, Farias T, Fontes P. Predicting radiative heat transfer in absorbing, emitting, and scattering media using the discrete transfer method[C]//Heat Transfer Division, New York:American Society of Mechanical Engineers, 1991, 160, 17-26.