基于改進(jìn)蒙特卡洛法求解輻射傳熱問(wèn)題的計(jì)算精度評(píng)價(jià)

李國(guó)軍，鐘佳琪，李定勇，王曉東

（1.東北大學(xué)冶金學(xué)院，遼寧沈陽(yáng) 110000；2.吉利汽車自動(dòng)化研究所，浙江寧波 315000）

能源動(dòng)力、航空航天等領(lǐng)域的傳熱過(guò)程通常以輻射傳熱為主，大多數(shù)情況該類問(wèn)題難于進(jìn)行解析求解而多借助于數(shù)值方法進(jìn)行求解，故而對(duì)該類問(wèn)題的數(shù)值計(jì)算方法研究具有重要意義［1-5］.蒙特卡洛方法（Monte Carlo method，MCM）可以精確地處理光譜特性、非均勻介質(zhì)、各向異性散射及復(fù)雜幾何形狀等復(fù)雜輻射計(jì)算問(wèn)題，已經(jīng)成為解決輻射傳熱的主要數(shù)值計(jì)算方法之一［6-9］.MCM 是一種統(tǒng)計(jì)模擬方法，其數(shù)值結(jié)果的精度隨著抽樣能束數(shù)目的增加而提高，而計(jì)算時(shí)間隨著追蹤能束數(shù)的增加而增加.為兼顧計(jì)算精度及計(jì)算時(shí)間，筆者提出了一種改進(jìn)MCM［10］，該方法可以在只進(jìn)行一次抽樣情況下完成輻射換熱求解，且其計(jì)算精度及計(jì)算效率高于傳統(tǒng)MCM.隨著MCM 在輻射傳熱數(shù)值模擬計(jì)算應(yīng)用范圍的不斷擴(kuò)展，如何定量分析和評(píng)估其計(jì)算結(jié)果的誤差及精度已成為關(guān)注的焦點(diǎn)，建立公認(rèn)的數(shù)值誤差分析和精度評(píng)估方法成為MCM 主要研究?jī)?nèi)容之一.Siegel 等［11］首次用不確定度來(lái)估計(jì)MCM 計(jì)算輻射傳熱問(wèn)題的統(tǒng)計(jì)誤差.Planas Almazan［12］運(yùn)用MCM 射線路徑軌跡量化了混合網(wǎng)交換公式固有的統(tǒng)計(jì)誤差.Plotnikov 和Shkarupa［13］應(yīng)用直接模擬MCM 求解稀薄氣體動(dòng)力學(xué)問(wèn)題的統(tǒng)計(jì)誤差.阮立明等［14］利用輻射交換因子的守恒性和互易性的檢驗(yàn)來(lái)評(píng)估MCM的計(jì)算誤差，得到一種間接評(píng)價(jià)MCM計(jì)算精度的方法.Yarahmadi等［15］提出一種改進(jìn)平均相對(duì)不確定評(píng)價(jià)方法，驗(yàn)證由于表面溫度和發(fā)射率的不確定性而引起的凈熱流密度不確定度的新表達(dá)式.Wang等人［16-17］提出了一種直接定量評(píng)價(jià)MCM 精確度的評(píng)價(jià)方法.

為評(píng)價(jià)改進(jìn)MCM 計(jì)算精度，本文擬采用改進(jìn)MCM 分別對(duì)漫灰表面、參與性各項(xiàng)同性介質(zhì)的方形及圓形封閉腔體內(nèi)的輻射傳熱問(wèn)題進(jìn)行研究，建立表面和體積微元輻射熱通量的誤差計(jì)算模型，得出其最小誤差與能束數(shù)的函數(shù)關(guān)系.采用直接定量評(píng)價(jià)方法開展基于輻射熱通量的改進(jìn)MCM計(jì)算誤差的分析和精度評(píng)價(jià)，研究網(wǎng)格密度及采樣能束數(shù)變化對(duì)改進(jìn)MCM求解輻射傳熱計(jì)算精度及效率的影響.

1 改進(jìn)蒙特卡洛法簡(jiǎn)介

應(yīng)用區(qū)域法求解輻射換熱時(shí)，當(dāng)表面段和體積段的參數(shù)確定后，各段之間的輻射傳遞因子也可通過(guò)計(jì)算獲得.若假定微元發(fā)射的全部能量到達(dá)其他各段能量的比例與該微元反射能量到達(dá)其他各段能量的比例相同，與發(fā)射能量份額及反射能量的份額無(wú)關(guān)，則采用MCM 求解輻射傳遞因子時(shí)，對(duì)全部微元發(fā)射的能束只需進(jìn)行一次采樣追蹤，以確定微元段發(fā)射能量到達(dá)其他微元段的比例.當(dāng)表面微元反射時(shí)，將反射能束按發(fā)射能束處理，且到達(dá)其他段的比例在之前已經(jīng)確定，無(wú)需再計(jì)算，散射情況也類似處理.該思路即為改進(jìn)MCM，具體求解方法如下.

1.1 微元發(fā)射能束直接到達(dá)其他微元段比例

為方便計(jì)算，將表面微元按順序依次命名為1，2，…，Ns，體積微元安排在表面微元之后為Ns+1，Ns+2，…，Ns+Ng，其中Ng和Ns分別為體積和表面微元總數(shù)，則總微元數(shù)為Ns+Ng.則熱交換場(chǎng)中微元i發(fā)射能束直接到達(dá)微元j比例為

式中：kUi,j（k=1，2，3，…）表示第i個(gè)微元發(fā)射的能量第k次循環(huán)到達(dá)第j個(gè)微元的能束數(shù)，Ni為第i段發(fā)射的總能束數(shù)，當(dāng)k=1時(shí)表示直接到達(dá).

1.2 微元發(fā)射能束到達(dá)其他微元段的總能束數(shù)

第i個(gè)微元發(fā)射的能量最終到達(dá)第j個(gè)微元的能量Ui,j分為直接到達(dá)被吸收的能量及經(jīng)k次反射后到達(dá)被吸收的能量的累加，則

式中：Ui,j為第i微元段發(fā)射能束到達(dá)第j微元段能束總數(shù)，其中包含直接到達(dá)與反射到達(dá)情況，ε為微元段黑度.

對(duì)體積微元j=Ns+1，Ns+2，…，Ns+Ng有

式中：ω為散射反照率，Ui,j為第i微元段發(fā)射能束到達(dá)第j微元段能束總數(shù)，其中包含直接到達(dá)與散射到達(dá)情況.

當(dāng)介質(zhì)為各項(xiàng)同性散射介質(zhì)時(shí)，

反照率為

式中：κe為含粒子介質(zhì)系的衰減系數(shù)，κsp為含粒子介質(zhì)系的粒子的衰減系數(shù)，ω表示散射介質(zhì)的散射反照率.

首先確定一個(gè)隨機(jī)數(shù)R，如果R＞ω，則能束被吸收，否則被散射.

確定散射方向是MCM 研究含粒子系輻射傳遞的關(guān)鍵.本文僅計(jì)算各向同性散射，已知散射相函數(shù)的歸一化條件

對(duì)于各項(xiàng)同性散射的散射相函數(shù)表達(dá)為

對(duì)于各項(xiàng)同性散射，即b=0時(shí)

1.3 計(jì)算終止條件

當(dāng)微元段i發(fā)射的能束經(jīng)過(guò)多次的反射或散射后剩余能量逐漸減少，當(dāng)剩余能量與發(fā)射能量的比值滿足式（10）時(shí)，式中ξ為無(wú)窮小量，則認(rèn)為計(jì)算精度已經(jīng)滿足要求.

2 直接定量評(píng)價(jià)法

將輻射體系劃分為M個(gè)面元和N個(gè)體積微元，則表面和體積微元凈輻射熱通量

因計(jì)算過(guò)程中存在誤差，則其真值可以描述為計(jì)算值與計(jì)算誤差之和，即

式中：Fm→i為表面微元m對(duì)表面微元i的輻射交換因子；Fn→j為體積微元n對(duì)體積微元j的輻射交換因子；ε是表面微元i的黑度；κa是氣體光學(xué)厚度；A、V分別為表面微元的面積和體積微元的體積；qai、qvj分別為表面微元i、體積微元j上的凈輻射熱通量；為計(jì)算值；Δqai、Δqvj為計(jì)算誤差.

在等溫輻射平衡狀態(tài)下，凈輻射通量理論上為0，則微元段吸收的輻射能量與發(fā)射的輻射能量理論上絕對(duì)相等，從而得到表面微元無(wú)因次方程為

由式（11）（13）（15），可得表面微元i的無(wú)因次凈輻射熱流與精確值的偏差

同理可得，體積微元j的無(wú)因次凈輻射熱流與精確值的偏差

則表面微元和體積微元的輻射熱通量相對(duì)誤差為

表面和體積微元輻射通量的相對(duì)均方根誤差表示為

3 結(jié)果與討論

3.1 網(wǎng)格密度對(duì)計(jì)算誤差的影響

以漫灰表面、參與性各向同性介質(zhì)的方形及圓形封閉腔體內(nèi)的輻射傳熱問(wèn)題為例，如圖1 所示，將邊長(zhǎng)為L(zhǎng)的方腔離散按長(zhǎng)度方向Nx和高度方向Ny均勻劃分網(wǎng)格；將半徑為R的圓形腔體沿徑向Nr及周向Nz劃分為均勻網(wǎng)格.本節(jié)中引入光學(xué)厚度τg=并分析網(wǎng)格劃分密度對(duì)計(jì)算誤差影響.計(jì)算設(shè)定條件為：設(shè)方腔離散網(wǎng)格數(shù)為Nx=Ny=5，10，20，40，圓形腔體離散網(wǎng)格數(shù)為Nz=4，14，28，40，Nr=2.微元發(fā)射能束N=105，表面發(fā)射率ε=0.5，介質(zhì)的散射反照率ω=0.5.為減小偽隨機(jī)數(shù)的隨機(jī)性對(duì)計(jì)算結(jié)果的影響，本文采用2 次計(jì)算結(jié)果及其平均值進(jìn)行研究.

圖1 算例圖Fig.1 The physical model of two-dimensional cavity

由于偽隨機(jī)數(shù)的隨機(jī)性，不同的離散網(wǎng)格密度的計(jì)算結(jié)果隨機(jī)誤差均值表現(xiàn)出不同的噪聲.以圖2（a）（b）（c）（d）和圖2（e）（f）（g）（h）所示，離散網(wǎng)格密度較小則誤差波動(dòng)明顯，離散網(wǎng)格密度越大，結(jié)果越穩(wěn)定，并且離散網(wǎng)格密度的增加使改進(jìn)MCM 計(jì)算結(jié)果的隨機(jī)性和任意性最小化到更小的程度.從圖2可以看出，方形網(wǎng)格數(shù)取Nx=Ny=40，圓形網(wǎng)格數(shù)取Nz=40，Nr=2即可滿足精度條件.

圖2 不同網(wǎng)格密度下體積微元輻射熱通量計(jì)算誤差Fig.2 Calculation error of radiant heat flux of gas elements with different grid density

3.2 輻射傳熱量的精確值和計(jì)算值

以漫灰表面、參與性各向同性介質(zhì)的方形及圓形封閉腔體內(nèi)的輻射傳熱問(wèn)題為例，將方形離散成均勻表面微元和體積微元，其中Nx=Ny=40.將圓形腔體沿徑向及周向分別劃分為均勻網(wǎng)格，即Nz=40，Nr=2.網(wǎng)格劃分如圖1所示.

當(dāng)微元發(fā)射能束數(shù)N=105，表面發(fā)射率ε方=0.8，ε圓=0.5，介質(zhì)的散射反照率ω=0.5，氣體光學(xué)厚度分別為τg=0.000 5，τg=0.5，τg=50 時(shí)，表面微元無(wú)因次熱通量計(jì)算結(jié)果分別如圖3所示.

由圖3 可知，光學(xué)厚度不同時(shí)，分別采用改進(jìn)MCM 及MCM 計(jì)算方腔及圓形腔內(nèi)表面微元無(wú)因次熱通量的計(jì)算值與真實(shí)值誤差較小.當(dāng)τg=0.00 5時(shí)，采用改進(jìn)MCM 求解得到的方腔及圓形腔的相對(duì)均方根誤差Ea分別為0.002 5 及0.002 3（MCM 求解得到的方腔及圓形腔的相對(duì)均方根誤差Ea分別為0.008 0 及0.003 7）；當(dāng)τg=0.5 時(shí)，其改進(jìn)MCM 求解的Ea分別為0.002 9 及0.0028（MCM 求解的Ea分別為0.014 1 及0.011 6）；當(dāng)τg=50 時(shí)，其改進(jìn)MCM 求解的Ea分別為0.003 2 及0.0031（MCM 求解的Ea分別為0.027 0及0.023 8）.由此可以看出，改進(jìn)MCM對(duì)方腔及圓形腔具有相同適用性，在相同條件下其求解精度高于MCM.

圖3 表面微元吸收的輻射熱通量的計(jì)算值和精度值Fig.3 Calculation and accuracy of radiative heat flux absorbed by surface elements

圖4給出了氣體光學(xué)厚度分別為τg=50，τg=0.5，τg=0.000 5 時(shí)方腔及圓形腔體積微元輻射熱通量計(jì)算值與真實(shí)值關(guān)系.由圖可知，當(dāng)τg=0.000 5 時(shí)，用改進(jìn)MCM 求解得到的方腔及圓形腔的相對(duì)均方根誤差Ev分別為0.0013 8 及0.026 8（MCM 求解的Ev分別為0.040 3 及0.097 8）；當(dāng)τg=0.5 時(shí)，方腔及圓形腔的相對(duì)均方根誤差Ev分別為0.014 2及0.051 2（MCM求解的Ev分別為0.040 3 及0.120 7）；當(dāng)τg=50 時(shí)，方腔及圓形腔的相對(duì)均方根誤差Ev分別為0.499 7 及0.146 2（MCM求解的Ev分別為1.001 8及0.478 2）.可見，當(dāng)其他條件不變時(shí)，隨著光學(xué)厚度增加，MCM 及改進(jìn)MCM計(jì)算結(jié)果誤差有增大趨勢(shì).

圖4 體積微元吸收的輻射熱通量的計(jì)算值和精度值Fig.4 Calculation and accuracy of radiative heat flux absorbed by gas elements

3.3 微元發(fā)射能束數(shù)與計(jì)算誤差關(guān)系

設(shè)方形及圓形的表面發(fā)射率ε=0.5，介質(zhì)的散射反照率ω=0.5，當(dāng)微元發(fā)射能束數(shù)N分別為N=1 000，3 000，5 000，10 000，100 000，500 000時(shí)，求解表面和體積微元無(wú)因次熱通量最小誤差值.

表1 給出了Ea，min、Ev，min與微元發(fā)射能束數(shù)之間的關(guān)系.可以看出，隨著微元發(fā)射能束數(shù)增加，最小誤差逐漸減小.對(duì)圓形封閉腔，無(wú)因次熱通量最小誤差值Ea，min與微元發(fā)射能束數(shù)N的關(guān)系可擬合為如式（22）所示.

表1 Ea，min和Ev，min隨N的變化Tab.1 The change of Ea，min and Ev，min with different N

由圖5 可知，計(jì)算精度隨著微元能束數(shù)的增加而提高，但是當(dāng)達(dá)到一定的微元能束數(shù)時(shí)，能束數(shù)繼續(xù)增加對(duì)計(jì)算精度提升逐漸不明顯.通常對(duì)于燃燒模擬，輻射傳熱計(jì)算中可接受的Ea，min誤差水平要求約為0.01，所以對(duì)于用改進(jìn)MCM 求解輻射傳熱的方形和圓形算例分別需要微元能束數(shù)N=5 000 和N=10 000 即為滿足精度要求.通過(guò)式（23）可以獲得不同最小誤差要求情況下需要的微元能束數(shù)目.

圖5 Ea，min和Ev，min與N的函數(shù)關(guān)系Fig.5 The functional relationship between Ea，min，Ev，min and N

4 結(jié)論

本文闡述了基于改進(jìn)MCM 進(jìn)行求解輻射換熱問(wèn)題的思路，并引入直接定量評(píng)價(jià)方法對(duì)不同光學(xué)厚度情況下的圓形腔及方腔內(nèi)參與性介質(zhì)的輻射熱通量計(jì)算精度進(jìn)行了評(píng)價(jià).研究結(jié)果如下：

1）改進(jìn)MCM 對(duì)方腔及圓形腔內(nèi)輻射問(wèn)題具有適用性，在相同計(jì)算條件下，其直接評(píng)價(jià)結(jié)果明顯高于采用改進(jìn)MCM 進(jìn)行計(jì)算的精度.如對(duì)于光學(xué)厚度為0.005 時(shí)，對(duì)方腔及圓形腔，采用改進(jìn)MCM 的Ea值分別0.002 5 及0.002 3 而采用MCM 時(shí)Ea分別為0.008 0及0.003 7.

2）研究了微元發(fā)射能束數(shù)與計(jì)算誤差關(guān)系，給出了能束數(shù)與Ea關(guān)系擬合式，為采用改進(jìn)MCM 進(jìn)行輻射換熱計(jì)算提供追蹤能束數(shù)選取的依據(jù).