蒙特-卡羅法計算黑體空腔有效發(fā)射率

方茜茜，方 偉 王 凱

(1.中國科學(xué)院 長春光學(xué)精密機械與物理研究所，吉林 長春130033;2.中國科學(xué)院 研究生院，北京100039)

1 引 言

黑體空腔的有效發(fā)射率在光度學(xué)、輻射度學(xué)等領(lǐng)域具有重要的應(yīng)用。有效發(fā)射率與腔的幾何尺寸、內(nèi)表面的光學(xué)性質(zhì)、腔壁溫度分布以及觀察條件密切相關(guān)。從20 世紀(jì)60 年代開始，科學(xué)家便通過實驗測量的方式確定黑體空腔的有效發(fā)射率［1-5］，由于輻射收集不完全、腔體溫度等重要參數(shù)測量不準(zhǔn)確、應(yīng)用不方便等因素，實驗測量并不是得到黑體空腔有效發(fā)射率的有效手段，而同步發(fā)展起來的理論計算方法克服了上述困難［6-8］。20 世紀(jì)80 年代，隨著計算機的廣泛應(yīng)用，蒙特-卡羅方法開始應(yīng)用于輻射熱傳遞等領(lǐng)域。該方法的優(yōu)點是易于處理復(fù)雜的空腔形狀和壁面輻射特性，物理解釋清晰，便于計算機處理，且不存在收斂性問題，所以是計算黑體空腔發(fā)射率的有效方法。對于任意形狀、內(nèi)表面光學(xué)性質(zhì)均一( 內(nèi)表面可能是均勻漫反射、鏡面反射或者均勻漫-鏡反射) 的等溫和不等溫腔體，該方法均取得了較好的結(jié)果［9-13］。

1973 年，Heinisch［14］首次將蒙特-卡羅法應(yīng)用于圓錐形腔的半球有效發(fā)射率的計算，其研究重點是腔內(nèi)光線追跡算法、不等溫腔修正項的計算以及如何提高計算速度等。目前，蒙特-卡羅法在計算黑體空腔發(fā)射率上已經(jīng)擁有扎實的理論基礎(chǔ)和相對固定的算法，并受到國內(nèi)外相當(dāng)多研究人員的關(guān)注，如清華大學(xué)、中國計量科學(xué)研究院等單位的研究人員［15-17］在20 世紀(jì)90 年代即采用蒙特-卡羅法對圓錐腔、圓柱腔、帶錐底的圓柱腔等不同形狀的黑體空腔進(jìn)行了有效發(fā)射率計算，并通過優(yōu)化隨機數(shù)產(chǎn)生算法，添加環(huán)境輻射修正、不等溫修正等得到了較好的計算結(jié)果，但在蒙特-卡羅法的物理思想和光線追跡算法上基本沿用國外的理論。目前，隨著計算機硬件的發(fā)展，蒙特-卡羅法的計算速度與以前相比得到很大提高，算法本身也得到進(jìn)一步的完善和優(yōu)化，但基本的光線追跡思想仍保持不變。隨著計算機軟件技術(shù)的迅速發(fā)展和廣泛使用，不僅能通過編程實現(xiàn)光線追跡算法，還可以使用optiCAD，Tracepro 等軟件實現(xiàn)，使得蒙特-卡羅方法的實施過程變得更加便捷。本文具體闡述了運用蒙特-卡羅法計算黑體空腔有效發(fā)射率的物理基礎(chǔ)和兩種光線追跡算法，討論了兩種光線追跡思想的利與弊，對前人的工作進(jìn)行了總結(jié)，為今后的研究提供參考。

2 黑體空腔的有效發(fā)射率

為了簡化運算，近似認(rèn)為腔體內(nèi)表面的光學(xué)性質(zhì)均一，與位置和波長無關(guān)。對均勻的漫-鏡反射模型做了如下近似:

(1) 半球反射率是指由某一方向入射的光在半球空間內(nèi)的反射率，它由鏡面反射ρs和漫反射ρd兩部分組成( ρ=ρs+ρd=1 -ε) ，近似認(rèn)為半球反射率ρ 不依賴于入射角。

(2) 表面漫反射率定義為D=ρd/ρ，也不依賴于入射角。

對于任意形狀的腔體，從腔口徑出射的光譜輻亮度隨方向變化而變化。溫度為T0的等溫腔的光譜定向有效發(fā)射率定義為:

式中:Lλ( λ，T0，ξ，ω) 為腔壁ξ 點沿ω 方向在波長λ 處的光譜輻亮度，Lλ，bb( λ，T0) 為理想黑體在相同溫度和波長處的光譜輻亮度。對于等溫灰體腔，εe( ξ，ω) 只與腔體的幾何形狀、光學(xué)性質(zhì)以及出射條件有關(guān)，與波長無關(guān)。

對于軸對稱黑體空腔，觀察方向平行于腔體的軸線并與腔體口徑相交于坐標(biāo)(x，y) 時的定向有效發(fā)射率被稱作垂直有效發(fā)射率εe，n(x，y) ，它是定向有效發(fā)射率的特例。

另一個重要的物理量是積分有效發(fā)射率。在距離腔體口徑Hd處放置一個半徑為Rd的共軸探測器，積分有效發(fā)射率是由溫度為T0的腔體發(fā)出，入射到探測器上的光譜輻射通量Φλ( 或總輻射通量Φ) 與和腔口徑面積相同、溫度相同的理想黑體輻射到相同探測器上的光譜輻射通量Φλbb( 或總輻射通量Φbb) 的比值。

對于等溫灰體腔:

平均垂直有效發(fā)射率可以認(rèn)為是當(dāng)Hd→∞，Rd=Ra時的積分有效發(fā)射率。

Hd=0，Rd=Ra時的積分有效發(fā)射率定義為腔的半球有效發(fā)射率。

3 熱輻射物理基礎(chǔ)

黑體空腔的各個有效發(fā)射率都直接或間接地與定向有效發(fā)射率相關(guān)，垂直有效發(fā)射率是入射光線平行于腔體軸線的定向有效發(fā)射率，對垂直有效發(fā)射率在整個腔體孔徑上積分就可得到平均垂直有效發(fā)射率，從后面的討論中還可以看到，對定向有效發(fā)射率積分還可以得到半球有效發(fā)射率。

在計算中求得定向有效發(fā)射率具有很重要的意義。A.Ono 曾在文獻(xiàn)［10］中進(jìn)行詳細(xì)的理論推導(dǎo)。

圖1 所示為面元dx1和dx3通過任意面元dx2進(jìn)行輻射交換的情況。反射率r1，32表示從dx1方向入射，經(jīng)面元dx2反射到dx3方向單位立體角的反射率，θ12 表示面元dx2的法線與dx1和dx2中心連線之間的夾角，θ32 表示面元dx2的法線與dx2和dx3中心連線之間的夾角，由相互性原理可得:

圖1 面元x1和x3通過中間表面x2的輻射熱交換Fig.1 Radiation heat exchange between surface element x1 and x3 through x2

定向半球反射率r1

2 定義為由面元dx1發(fā)出經(jīng)過dx2反射到半球空間的反射率。

對于不透明體，由基爾霍夫輻射定律得到:

式中，e1

2 是由dx2向dx1方向的輻射通量與相同面積、相同溫度的黑體在相同方向的輻射通量的比值，稱為定向有效發(fā)射率。

任意形狀，開口面積為A的等溫腔如圖2 所示，觀察方向為D，則由x0處發(fā)出到達(dá)D方向的定向光輻射通量為:

圖2 任意形狀的等溫腔體Fig.2 Isothermal cavity with arbitrary shape

式中，Lb是理想黑體的輻亮度。dΦD0 由面元dx0直接發(fā)射至D方向以及經(jīng)x0反射到立體角兩部分輻射通量組成。因此:

式中，W是x0可視的所有面積集合。

由式(8) 、(9) 可得定向有效發(fā)射率的簡化積分方程:

根據(jù)相互性原理，上式簡化為:

式中，ρD0表示由方向D入射到腔體x0處，并最終反射到半球空間的定向半球反射率。將式( 6) 、(7) 、(11) 帶入式(12) 中，可得:

式(13) 的解可以表示為:

式中，fi表示從方向D入射到x0處，經(jīng)過i次反射射出腔體的輻射量。

在過去的30 多年間，研究人員已經(jīng)成功地對球體、圓錐、圓柱及組合形狀等腔的定向有效發(fā)射率進(jìn)行蒙特-卡羅計算，雖然在具體計算過程中采用的算法不同，但通過對每束光線每次反射出腔口徑的輻射量進(jìn)行求和，用于計算腔的半球反射率，然后由基爾霍夫輻射定律得到腔的定向有效發(fā)射率的物理思想是基本相同的。

4 光線追跡算法

光線追跡算法是計算腔體有效發(fā)射率的核心。在國際上曾經(jīng)有正向追跡和逆向追跡兩種方法，二者的區(qū)別是起始發(fā)出光線的位置不同。正向追跡是光線從腔內(nèi)壁一點發(fā)出，然后對光線進(jìn)行追跡，逆向追跡是光線由腔外觀察點發(fā)出，當(dāng)光線入射到腔內(nèi)部后對光線追跡。由于逆向追跡的思想在計算定向有效發(fā)射率時很方便，所以到后期一般都采用逆向追跡的方法。

不等溫腔的有效發(fā)射率只是在等溫情況下做了一些修正，這部分討論均是針對等溫腔，不等溫情況將在后面討論。假設(shè)腔內(nèi)表面是光學(xué)性質(zhì)均勻的漫-鏡反射體，反射率ρ=ρd+ρs( ρd為漫反射率，ρs為鏡面發(fā)射率) 。

4.1 逆向光線追跡

逆向光線追跡是光線由腔外觀察點發(fā)出，入射到腔內(nèi)后對光線追跡。該方法可以方便地計算定向有效發(fā)射率，經(jīng)過積分運算后可以得到平均垂直有效發(fā)射率和半球有效發(fā)射率。

逆向光線追跡的算法一般分為以下4 個步驟:

(1) 假設(shè)一束光線從D方向入射至腔的ξ0點，在ξ0處一部分光線被吸收，一部分光線被反射。由偽隨機數(shù)產(chǎn)生器生成( 0，1) 之間的隨機數(shù)η，如果η ＜D，則腔體表面發(fā)生漫反射，否則發(fā)生鏡面反射。

如果光線第一次和腔壁相互作用發(fā)生鏡面反射，則反射光線的方向可以通過下式來確定:

式中，ωi，ωr，n 分別表示入射光矢量、反射光矢量以及腔壁ξ0的法線方向。

漫反射情況稍微復(fù)雜一些，近似認(rèn)為純漫反射的光能量均勻分布在半球空間中，所以在半球空間隨機均勻地產(chǎn)生一個方向代表該束光的漫射方向，如果反射光直接射出腔體，則結(jié)束該束光的追跡，否則繼續(xù)對該束光線追跡。

在半球空間中同等概率地產(chǎn)生漫射方向可以歸結(jié)為在半球上均勻地產(chǎn)生一點，該點的坐標(biāo)就可以作為漫射的方向矢量。Marsaglia［13］曾全面系統(tǒng)地總結(jié)了關(guān)于在球面上均勻地產(chǎn)生一點的問題，他不僅對算法的合理性、可行性予以考慮，還從計算速度、收斂特性上進(jìn)行了討論。通常采用的算法也是最容易理解的算法，即是在( -1，1)之間隨機產(chǎn)生3 個數(shù)時，球面上隨機點的坐標(biāo)為:

Marsaglia 提出了一種新的方法，基本思想是在球體z軸上任意選取一點，過該點作平行于xy面的平面，該平面與半徑為1 的球的交線是半徑為的圓周。在圓周上隨機產(chǎn)生一點。具體算法是產(chǎn)生偽隨機數(shù)V1、V2∈( -1，1) ，當(dāng)S=時，球面隨機點的坐標(biāo)為( 2V1( 1 -由于該算法與以往相比減少了偽隨機數(shù)的個數(shù)和平方根的運算次數(shù)，計算速度提高了約2 倍。

(2) 追跡第一次反射的光線與腔壁作用點的坐標(biāo)。目前研究人員基本采用相同的算法追跡光線與腔壁作用點的坐標(biāo)。

式中:Φ(x，y，z) =0 為腔表面方程，ξ0為發(fā)出光線的位置坐標(biāo)，ξ 為作用點位置坐標(biāo)，ω 為反射方向矢量，t為系數(shù)。

(3) 光線第二次和腔壁發(fā)生相互作用，重復(fù)(1) ～(2) 的分析過程，直到光線射出腔口徑或者經(jīng)多次反射腔內(nèi)的光能量減小到可以忽略為止。( 假設(shè)入射光束總能量是1，多次反射能量衰減為10-5或10-6可結(jié)束追跡) 。

(4) 結(jié)束本條光束的追跡，進(jìn)行下一束光的追跡。為了提高蒙特-卡羅法的計算精度，對同一條件入射的情況，一般需要追跡105～107條光束。

根據(jù)上述逆向光束追跡過程，可以推導(dǎo)出黑體空腔定向有效發(fā)射率的計算方程。Sapritsky［9］給出的方程如下:

Sapritsky 總共對N條光束進(jìn)行追跡，其中第i條光束經(jīng)過M次反射后射出腔體。F是漫反射角因子，具體表達(dá)式如下:

式中:Ω 表示ξ 與腔口徑所成的立體角; θξ是腔壁ξ 處的法線與立體角dΩ 軸線之間的夾角。F是漫反射光通過腔口徑射出腔體的部分。

由此可以看出，Sapritsky 認(rèn)為每一條光束只要經(jīng)歷一次漫反射，能量就減小1 -F( ξ) 倍。這對鏡面反射也是成立的，如果光束經(jīng)鏡面反射剛好射出腔體，則F( ξ) =1，否則F( ξ) =0。

Prokhorov［11］采用更簡單的算法:

上式各物理量的含義同式( 18) ，只是Prokhorov 將漫反射和鏡面反射同等對待，只是在確定兩者的反射方向上有所區(qū)別。

4.2 正向光線追跡

正向光線追跡可以方便地計算具有漫-鏡反射表面腔的半球有效發(fā)射率。Heinisch［14］曾采用正向追跡方法較精確計算漫反射等溫圓錐形腔的半球有效發(fā)射率。這里將計算推廣到任意形狀、具有漫-鏡反射表面的等溫腔:

(1) 腔壁發(fā)出光束的起始點坐標(biāo)ξ 由偽隨機數(shù)確定。理想情況下，在趨于無窮次的光束追跡過程中，起始發(fā)出光線點的位置應(yīng)均勻地分布在整個腔壁上。

(2) 從ξ 發(fā)出的光能量分為兩部分，一部分能量E·Fi直接從腔口射出，另一部分能量E(1 -Fi) 在腔內(nèi)沿隨機方向(Rθ、Rφ) 傳播。傳至下一點時，光束能量將被隨機地吸收或反射，這取決于0 ～1 的隨機數(shù)Rα，當(dāng)Rα≤ε 時，能量被吸收，否則光束被反射。反射的能量E(1 -Fi) 再次分為兩部分，分別為E(1 -Fi)Fi1和E(1 -Fi) ·(1 -Fi1) 。重復(fù)上面的過程直至能量在某一點處被吸收。

(3) 選另一發(fā)光點，重復(fù)上面的過程。

半球有效發(fā)射率的計算方程:

式中:Gi=(1 -Fi)Fi1+(1 -Fi) (1 -Fi1)Fi2+…+(1 -Fi) ( 1 -Fi1)Fi2…( 1 -Fim-1)Fim;m=1，2…表示能量在第m次反射之后被吸收;Fi為第i個發(fā)光點所在微元對腔口的角系數(shù);Fim為第i個發(fā)光點發(fā)出的能量在第m次反射時所在微元對腔口的角系數(shù)。

5 有效發(fā)射率的計算

上面講述了采用逆向追跡的方法計算黑體空腔定向有效發(fā)射率的具體算法。垂直有效發(fā)射率是光束垂直于孔徑的特殊定向有效發(fā)射率，它的算法和定向有效發(fā)射率相同。

積分有效發(fā)射率εe(Rd，Hd) 是半徑為Rd的圓形探測器共軸放置在與腔口徑相距Hd的位置時，由探測器測量得到的腔體發(fā)射率。它是計算平均垂直有效發(fā)射率和半球有效發(fā)射率的基礎(chǔ)。

圖3 計算積分有效發(fā)射率簡圖Fig.3 Diagram of calculating integrated effective emissivity

Prokhorov［11］給出了積分有效發(fā)射率的計算公式，如圖3 所示，由探測器接收到的總輻射通量表達(dá)式如下:

腔口徑面元dSa出射的輻射沿ψ 向入射到探測器dSd，L是P點沿PQ方向的輻亮度。如果繼續(xù)采用逆向光線追跡的方法，認(rèn)為光線由均勻分布在探測器圓面上的Qi發(fā)出，經(jīng)過均勻分布在腔口徑上的Pi點入射到腔體內(nèi)部，之后的光線追跡過程和上面介紹的定向有效發(fā)射率完全相同。那么上面的積分運算可以變?yōu)榍蠛瓦\算:

將黑體空腔換為面積為Sa的理想黑體，在相同條件下，探測器接收的輻射通量:

式中，F(xiàn)a-d為理想黑體對探測器圓面Sd的角系數(shù)。

此時黑體空腔的積分有效發(fā)射率:

6 不等溫腔的定向有效發(fā)射率

Sapritsky［9］提出了不等溫腔的定向有效發(fā)射率，它是在等溫情況下加上修正因子，如下所示:

式中:Tξ為腔壁ξ 點的實際溫度，T0為參考恒定溫度。在逆向光線追跡的過程中，不等溫修正項如下:

式中:Tk為光線第k次在腔壁上反射時腔壁作用點處的實際溫度;Le( λ，Tk) 為溫度為Tk，波長為λ 的光譜輻亮度。

7 結(jié) 論

本文綜述了運用蒙特-卡羅方法計算黑體空腔有效發(fā)射率的理論基礎(chǔ)和相應(yīng)的光線追跡算法。不難看出，該方法能夠方便、精確地計算黑體空腔的有效發(fā)射率。在計算具有漫反射表面的圓柱形腔的積分有效發(fā)射率時( 包括半球有效發(fā)射率和平均垂直有效發(fā)射率) ，得到的不確定度為0.000 1 ～0.000 2，與一些研究人員采用相關(guān)理論計算方法得到的不確定度大致相同［18-19］。大部分采用理論計算或者實驗測試的方法得到的腔體發(fā)射率結(jié)果與蒙特-卡羅法得到的結(jié)果均符合得很好。目前，蒙特-卡羅方法基本上已經(jīng)完全取代理論計算的方法，成為獲得黑體空腔有效發(fā)射率的最主要途徑。

［1］ KELLY F J，MOORE D G. A test of analytical expressions for the thermal emissivity of shallow cylindrical cavities［J］.Appl． Opt．，1965，4(1) ：31-40.

［2］ HEINISCH R P，SCHMIDT R N. Development and application of an instrument for the measurement of directional emittance of blackbody cavities［J］.Appl． Opt．，1970，9(8) ：1920-1925.

［3］ BAUER G，BISCHOFF K. Evaluation of the emissivity of a cavity source by reflection measurements［J］，Appl． Opt．，1971，10(12) ：2639-2643.

［4］ BALLICO M. Limitations of the Welch-Satterthwaite approximation for measurement uncertainty calculations［J］.Metrologia，2000，37(1) ：295-300.

［5］ GALAL Y S，SPERFELD P，METZDORF J. Measurement and calculation of the emissivity of a high-temperature black body［J］.Metrologia，2000，37(5) ：365-368.

［6］ BEDFORD R E .Temperature：Its Measurement and Control in Science and Industry［M］. New York：Reinhold Publishing Corp.，1962.

［7］ BARTELL F O，WOLFE W L. Cavity radiators：an ecumenical theory［J］.Appl． Opt．，1976，15(1) ：84-88.

［8］ GEIST J. Theoretical analysis of laboratory blackbodies 1：a generalized integral equation［J］.Appl． Opt．，1973，12(6) ：1325-1330.

［9］ SAPRITSKY V I，PROKHOROV A V. Calculation of the effective emissivities of specular-diffuse cavities by the Monte-Carlo method［J］.Metrologia，1992，29(1) ：9-14.

［10］ ONO A. Calculation of the directional emissivities of cavities by the Monte-Carlo method［J］.J． Opt． Soc． Am．，1980，70(5) ：547-554.

［11］ PROKHOROV A V，HANSSEN L M. Effective emissivity of a cylindrical cavity with an inclined bottom： I. Isothermal cavity［J］.Metrologia，2004，41(6) ：421-431.

［12］ PROKHOROV A V，HANSSEN L M. Effective emissivity of a cylindrical cavity with an inclined bottom：II.non-isothermal cavity［J］.Metrologia，2010，47(1) ：33-46.

［13］ MARSAGLIA G. Choosing a point from the surface of a sphere［J］.The Annals Mathematical Statistics，1972，43(2) ：645-646.

［14］ HEINISCH R P. Radiant emission from baffled conical cavities［J］.J． Opt． Soc． Am．，1973，63(2) ：152-158.

［15］ 張宏.黑體空腔發(fā)射率求解中的Monte-Carlo 法［J］.哈爾濱科學(xué)技術(shù)大學(xué)學(xué)報，1996，20(3) :72-76.ZHANG H. Monte-Carlo method in calculating the radiant emission of blackbody cavity［J］.J． Harbin University Sci．Technol．，1996，20(3) :72-76.( in Chinese)

［16］ 黃東濤，陸家欽，段宇寧.黑體空腔發(fā)射率計算的蒙特卡羅模型［J］.清華大學(xué)學(xué)報，1997，37(2) :28-31.HUANG D T，LU J Q，DUAN Y N. Monte-Carlo model for calculation of cavity effective emissivities［J］.J． Tsinghua University，1997，37(2) :28-31.( in Chinese)

［17］ 孫富韜.Monte-Carlo 方法在黑體空腔有效發(fā)射率計算中的應(yīng)用［J］.宇航計測技術(shù)，2010，30(2) :26-29.SUN F T. Calculation of the effective emissivities of blackbody cavity by the Monte-Carlo method［J］.J． Astronautic Metrology and Measurement，2010，30(2) :26-29.( in Chinese)

［18］ CHEN S，CHU Z，CHEN H. Precise calculation of the integrated emissivity of baffled blackbody cavities［J］.Metrologia，1980，16(2) :69-72.

［19］ CHANDOS R J，CHANDOS R E. Radiometric properties of isothermal，diffuse wall cavity sources［J］.Appl． Opt．，1974，13(9) :2142-2152.