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

方茜茜，方 偉 王 凱

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

1 引 言

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

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

2 黑體空腔的有效發射率

為了簡化運算，近似認為腔體內表面的光學性質均一，與位置和波長無關。對均勻的漫-鏡反射模型做了如下近似:

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

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

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

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

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

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

對于等溫灰體腔:

平均垂直有效發射率可以認為是當Hd→∞，Rd=Ra時的積分有效發射率。

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

3 熱輻射物理基礎

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

在計算中求得定向有效發射率具有很重要的意義。A.Ono 曾在文獻［10］中進行詳細的理論推導。

圖1 所示為面元dx1和dx3通過任意面元dx2進行輻射交換的情況。反射率r1，32表示從dx1方向入射，經面元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發出經過dx2反射到半球空間的反射率。

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

式中，e1

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

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

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

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

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

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

根據相互性原理，上式簡化為:

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

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

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

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

4 光線追跡算法

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

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

4.1 逆向光線追跡

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.2 正向光線追跡

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

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

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

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

半球有效發射率的計算方程:

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

5 有效發射率的計算

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

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

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

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

腔口徑面元dSa出射的輻射沿ψ 向入射到探測器dSd，L是P點沿PQ方向的輻亮度。如果繼續采用逆向光線追跡的方法，認為光線由均勻分布在探測器圓面上的Qi發出，經過均勻分布在腔口徑上的Pi點入射到腔體內部，之后的光線追跡過程和上面介紹的定向有效發射率完全相同。那么上面的積分運算可以變為求和運算:

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

式中，Fa-d為理想黑體對探測器圓面Sd的角系數。

此時黑體空腔的積分有效發射率:

6 不等溫腔的定向有效發射率

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

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

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

7 結 論

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

［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］ 張宏.黑體空腔發射率求解中的Monte-Carlo 法［J］.哈爾濱科學技術大學學報，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］ 黃東濤，陸家欽，段宇寧.黑體空腔發射率計算的蒙特卡羅模型［J］.清華大學學報，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 方法在黑體空腔有效發射率計算中的應用［J］.宇航計測技術，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.