李國軍 鐘佳琪 李定勇 王曉東

摘要：針對傳統(tǒng)蒙特卡洛法計算輻射傳輸耗時問題，提出了一種改進蒙特卡洛方法，通過比例迭代累加法來求解反射及散射能量，從而大幅減少了計算時間.引入直接評價方法，以包含參與性介質的密閉系統(tǒng)（方形和圓形為例）為例，分析了網格密度、發(fā)射能束數(shù)及物性參數(shù)對改進蒙特卡洛法計算精度的影響.當光學厚度為0.005時，采用改進蒙特卡洛方法求得方腔及圓形腔的表面微元輻射通量的相對均方根誤差值分別為0.0025及0.0023，而采用傳統(tǒng)蒙特卡洛方法時對應誤差分別為0.0080及0.0037.可見，相同計算條件下，改進蒙特卡洛方法對輻射換熱問題具有更高的精度.進一步研究了追蹤能束數(shù)對計算誤差的影響規(guī)律，給出了計算誤差與追蹤能束數(shù)擬合關系式，為計算該類問題的能束數(shù)選取提供了支撐.

關鍵詞：蒙特卡洛法;評價;計算精度;計算效率

中圖分類號：TK124

文獻標志碼：A

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

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

1改進蒙特卡洛法簡介

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

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

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

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

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

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

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

對體積微∑元j=Ns+1，Ns+2，...，Ns+Ng有

式中：ω為散射反照率，U達第j微元段能束總數(shù)，其中包含直接到達與散射到達情況.

當介質為各項同性散射介質時，

式中：K。為含粒子介質系的衰減系數(shù)，Km為含粒子介質系的粒子的衰減系數(shù)，ω表示散射介質的散射反照率.

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

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

對于各項同性散射的散射相函數(shù)表達為

對于各項同性散射，即b=0時

1.3計算終止條件

當微元段i發(fā)射的能束經過多次的反射或散射后剩余能量逐漸減少，當剩余能量與發(fā)射能量的比值滿足式（10）時，式中ξ為無窮小量，則認為計算精度已經∑滿足要求.

2直接定量評價法

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

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

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

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

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

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

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

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

式中：Q0ai、Q0vj分別為表面微元i、體積微元j上吸收無因次熱流的計算值;ΔQai、ΔQvj分別為表面微元i、體積微元j上吸收無因次熱流的計算誤差.

3結果與討論

3.1網格密度對計算誤差的影響以漫灰表面、參與性各向同性介質的方形及圓

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

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

3.2輻射傳熱量的精確值和計算值

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

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

由圖3可知，光學厚度不同時，分別采用改進MCM及MCM計算方腔及圓形腔內表面微元無因次熱通量的計算值與真實值誤差較小.當τg=0.005時，采用改進MCM求解得到的方腔及圓形腔的相對均方根誤差Ea分別為0.0025及0.0023（MCM求解得到的方腔及圓形腔的相對均方根誤差Ea分別為0.0080及0.0037）;當τg=0.5時，其改進MCM求解的Ea分別為0.0029及0.0028（MCM求解的Ea分別為0.0141及0.0116）;當τg=50時，其改進MCM求解的Ea分別為0.0032及0.0031（MCM求解的Ea分別為0.0270及0.0238）.由此可以看出，改進MCM對方腔及圓形腔具有相同適用性，在相同條件下其求解精度高于MCM.

圖4給出了氣體光學厚度分別為τg=50，τg=0.5，τg=0.0005時方腔及圓形腔體積微元輻射熱通量計算值與真實值關系.由圖可知，當τg=0.0005時，用改進MCM求解得到的方腔及圓形腔的相對均方根誤差Ev分別為0.00138及0.0268（MCM求解的Ev分別為0.0403及0.0978）;當τg=0.5時，方腔及圓形腔的相對均方根誤差Ev分別為0.0142及0.0512（MCM求解的Ev分別為0.0403及0.1207）;當τg=50時，方腔及圓形腔的相對均方根誤差Ev分別為0.4997及0.1462（MCM求解的Ev分別為1.0018及0.4782）.可見，當其他條件不變時，隨著光學厚度增加，MCM及改進MCM計算結果誤差有增大趨勢.

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

設方形及圓形的表面發(fā)射率ε=0.5，介質的散射反照率ω=0.5，當微元發(fā)射能束數(shù)N分別為N=1000，3000，5000，10000，100000，500000時，求解表面和體積微元無因次熱通量最小誤差值.

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

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

4結論

本文闡述了基于改進MCM進行求解輻射換熱問題的思路，并引入直接定量評價方法對不同光學厚度情況下的圓形腔及方腔內參與性介質的輻射熱通量計算精度進行了評價.研究結果如下：

1）改進MCM對方腔及圓形腔內輻射問題具有適用性，在相同計算條件下，其直接評價結果明顯高于采用改進MCM進行計算的精度.如對于光學厚度為0.005時，對方腔及圓形腔，采用改進MCM的Ea值分別0.0025及0.0023而采用MCM時Ea分別為0.0080及0.0037.

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

參考文獻

[1]李國軍，陳海耿，伊智.積分降重法求解輻射直接交換面積[J].東北大學學報（自然科學版），2010，31（1）：80-83.

[2] MODEST M F. Radiative heat transfer[M]. New York：Academic Press，1993：639-666.

[3] RUAN L M，TAN H P，YAN Y Y. A Monte Carlo（mc）method applied to the medium with nongray absorbing-emitting- anisotropic scattering particles and gray approximation[J].Nu?merical Heat Transfer，Part A：Applications，2002，42（3）： 253-268.

[4] TAN H P，SHUAI Y，XIA X L. Reliability of stray light calcula?tion code by the Monte Carlo method[J]. Optical Engineering， 2005，44（2）：1-2.

[5] KOVTANYUK A E，BOTKIN N D，HOFFMANN K H.Numerical simulations of a coupled radiative-conductive heat transfer model using a modified Monte Carlo method[J].International Journal of Heat and Mass Transfer，2012，55（4）：649-654.

[6] TSENG H Y，HONG Y B，WU C Y，et al.Monte Carlo simulation of radiative transfer in a medium with varying refractive index specified at discrete points[J].Applied Mathematical Modelling， 2 0 1 7 ，4 8 ：8 7 0 - 8 8 4 .

[7] WANG C H，ZHANG Y，YI H L，et al.Transient radiative trans?fer in two-dimensional graded index medium by Monte Carlo method combined with the time shift and superposition principle [J]. Numerical Heat Transfer，Part A：Applications，2016，69 （ 6 ）：5 7 4 - 5 8 8 .

[8]龔光彩，劉佳.空氣載能輻射空調混合通風協(xié)同運行研究[J].湖南大學學報（自然科學版），2019，46（5）：148-156.

[9] KHANNA S，YADAV U. Compution of total exchange ares using Monte Carlo method[J]. Journal of Engineering Thermophysics， 2013，8（3）：1-7.

[10] LI G J，ZHONG J Q，WANG X D. An improved Monte Carlo method for radiative heat transfer in participating media[J].Jour?nal of Quantitative Spectroscopy and Radiative Transfer，2020， 2 5 1 ：1 0 7 0 8 1 .

[11] SIEGEL R，HOWELL J R. Thermal radiation heat transfer[M]. NewYork：Taylor&Francis，2002：23-566.

[12] PLANAS ALMAZAN P. Accuracy of Monte Carlo ray-tracing ther?mal radiation calculations：A practical discussion[J]. Noordwijk： European Space Agency，1997：579-591.

[13] PLOTNIKOV M Y，SHKARUPA E V.Estimation of the statisti?cal error of the direct simulation Monte Carlo method[J].Compu?tational Mathematics and Mathematical Physics，2010，50（2）： 335-344.

[14]阮立明，譚建宇，董士奎，等.輻射傳遞蒙特卡洛法精度分析及數(shù)值試驗[J].工程熱物理學報，2003，24（5）：813-816.

[15] YARAHMADI M，MAHAN J R，PRIESTLEY J. Estimation and use of the radiation distribution factor median for predicting uncer?tainty in the Monte Carlo ray-trace method[J]. Journal of Heat Transfer，2019，141（6）：1- 7.

[16] WANG D D，ZHOU H C.Quantitative evaluation of the computa?tional accuracy for the Monte Carlo calculation of radiative heat transfer[J].Journal of Quantitative Spectroscopy and Radiative T r a n s f e r ，2 0 1 9 ，2 2 6 ：1 0 0 - 1 1 4 .

[17] WANG D D，CHENG Q，ZHUO X S，et al. Effects of surface emissivity and medium scattering albedo on the computational ac?curacy of radiative heat transfer by MCM[J].Journal of Quantita?tive Spectroscopy and Radiative Transfer，2020，240：106712.