方茜茜,方 偉 王 凱
(1.中國(guó)科學(xué)院 長(zhǎng)春光學(xué)精密機(jī)械與物理研究所,吉林 長(zhǎng)春130033;2.中國(guó)科學(xué)院 研究生院,北京100039)
黑體空腔的有效發(fā)射率在光度學(xué)、輻射度學(xué)等領(lǐng)域具有重要的應(yīng)用。有效發(fā)射率與腔的幾何尺寸、內(nèi)表面的光學(xué)性質(zhì)、腔壁溫度分布以及觀察條件密切相關(guān)。從20 世紀(jì)60 年代開(kāi)始,科學(xué)家便通過(guò)實(shí)驗(yàn)測(cè)量的方式確定黑體空腔的有效發(fā)射率[1-5],由于輻射收集不完全、腔體溫度等重要參數(shù)測(cè)量不準(zhǔn)確、應(yīng)用不方便等因素,實(shí)驗(yàn)測(cè)量并不是得到黑體空腔有效發(fā)射率的有效手段,而同步發(fā)展起來(lái)的理論計(jì)算方法克服了上述困難[6-8]。20 世紀(jì)80 年代,隨著計(jì)算機(jī)的廣泛應(yīng)用,蒙特-卡羅方法開(kāi)始應(yīng)用于輻射熱傳遞等領(lǐng)域。該方法的優(yōu)點(diǎn)是易于處理復(fù)雜的空腔形狀和壁面輻射特性,物理解釋清晰,便于計(jì)算機(jī)處理,且不存在收斂性問(wèn)題,所以是計(jì)算黑體空腔發(fā)射率的有效方法。對(duì)于任意形狀、內(nèi)表面光學(xué)性質(zhì)均一( 內(nèi)表面可能是均勻漫反射、鏡面反射或者均勻漫-鏡反射) 的等溫和不等溫腔體,該方法均取得了較好的結(jié)果[9-13]。
1973 年,Heinisch[14]首次將蒙特-卡羅法應(yīng)用于圓錐形腔的半球有效發(fā)射率的計(jì)算,其研究重點(diǎn)是腔內(nèi)光線追跡算法、不等溫腔修正項(xiàng)的計(jì)算以及如何提高計(jì)算速度等。目前,蒙特-卡羅法在計(jì)算黑體空腔發(fā)射率上已經(jīng)擁有扎實(shí)的理論基礎(chǔ)和相對(duì)固定的算法,并受到國(guó)內(nèi)外相當(dāng)多研究人員的關(guān)注,如清華大學(xué)、中國(guó)計(jì)量科學(xué)研究院等單位的研究人員[15-17]在20 世紀(jì)90 年代即采用蒙特-卡羅法對(duì)圓錐腔、圓柱腔、帶錐底的圓柱腔等不同形狀的黑體空腔進(jìn)行了有效發(fā)射率計(jì)算,并通過(guò)優(yōu)化隨機(jī)數(shù)產(chǎn)生算法,添加環(huán)境輻射修正、不等溫修正等得到了較好的計(jì)算結(jié)果,但在蒙特-卡羅法的物理思想和光線追跡算法上基本沿用國(guó)外的理論。目前,隨著計(jì)算機(jī)硬件的發(fā)展,蒙特-卡羅法的計(jì)算速度與以前相比得到很大提高,算法本身也得到進(jìn)一步的完善和優(yōu)化,但基本的光線追跡思想仍保持不變。隨著計(jì)算機(jī)軟件技術(shù)的迅速發(fā)展和廣泛使用,不僅能通過(guò)編程實(shí)現(xiàn)光線追跡算法,還可以使用optiCAD,Tracepro 等軟件實(shí)現(xiàn),使得蒙特-卡羅方法的實(shí)施過(guò)程變得更加便捷。本文具體闡述了運(yùn)用蒙特-卡羅法計(jì)算黑體空腔有效發(fā)射率的物理基礎(chǔ)和兩種光線追跡算法,討論了兩種光線追跡思想的利與弊,對(duì)前人的工作進(jìn)行了總結(jié),為今后的研究提供參考。
為了簡(jiǎn)化運(yùn)算,近似認(rèn)為腔體內(nèi)表面的光學(xué)性質(zhì)均一,與位置和波長(zhǎng)無(wú)關(guān)。對(duì)均勻的漫-鏡反射模型做了如下近似:
(1) 半球反射率是指由某一方向入射的光在半球空間內(nèi)的反射率,它由鏡面反射ρs和漫反射ρd兩部分組成( ρ=ρs+ρd=1 -ε) ,近似認(rèn)為半球反射率ρ 不依賴于入射角。
(2) 表面漫反射率定義為D=ρd/ρ,也不依賴于入射角。
對(duì)于任意形狀的腔體,從腔口徑出射的光譜輻亮度隨方向變化而變化。溫度為T(mén)0的等溫腔的光譜定向有效發(fā)射率定義為:
式中:Lλ( λ,T0,ξ,ω) 為腔壁ξ 點(diǎn)沿ω 方向在波長(zhǎng)λ 處的光譜輻亮度,Lλ,bb( λ,T0) 為理想黑體在相同溫度和波長(zhǎng)處的光譜輻亮度。對(duì)于等溫灰體腔,εe( ξ,ω) 只與腔體的幾何形狀、光學(xué)性質(zhì)以及出射條件有關(guān),與波長(zhǎng)無(wú)關(guān)。
對(duì)于軸對(duì)稱黑體空腔,觀察方向平行于腔體的軸線并與腔體口徑相交于坐標(biāo)(x,y) 時(shí)的定向有效發(fā)射率被稱作垂直有效發(fā)射率εe,n(x,y) ,它是定向有效發(fā)射率的特例。
另一個(gè)重要的物理量是積分有效發(fā)射率。在距離腔體口徑Hd處放置一個(gè)半徑為Rd的共軸探測(cè)器,積分有效發(fā)射率是由溫度為T(mén)0的腔體發(fā)出,入射到探測(cè)器上的光譜輻射通量Φλ( 或總輻射通量Φ) 與和腔口徑面積相同、溫度相同的理想黑體輻射到相同探測(cè)器上的光譜輻射通量Φλbb( 或總輻射通量Φbb) 的比值。
對(duì)于等溫灰體腔:
平均垂直有效發(fā)射率可以認(rèn)為是當(dāng)Hd→∞,Rd=Ra時(shí)的積分有效發(fā)射率。
Hd=0,Rd=Ra時(shí)的積分有效發(fā)射率定義為腔的半球有效發(fā)射率。
黑體空腔的各個(gè)有效發(fā)射率都直接或間接地與定向有效發(fā)射率相關(guān),垂直有效發(fā)射率是入射光線平行于腔體軸線的定向有效發(fā)射率,對(duì)垂直有效發(fā)射率在整個(gè)腔體孔徑上積分就可得到平均垂直有效發(fā)射率,從后面的討論中還可以看到,對(duì)定向有效發(fā)射率積分還可以得到半球有效發(fā)射率。
在計(jì)算中求得定向有效發(fā)射率具有很重要的意義。A.Ono 曾在文獻(xiàn)[10]中進(jìn)行詳細(xì)的理論推導(dǎo)。
圖1 所示為面元dx1和dx3通過(guò)任意面元dx2進(jìn)行輻射交換的情況。反射率r1,32表示從dx1方向入射,經(jīng)面元dx2反射到dx3方向單位立體角的反射率,θ12 表示面元dx2的法線與dx1和dx2中心連線之間的夾角,θ32 表示面元dx2的法線與dx2和dx3中心連線之間的夾角,由相互性原理可得:
圖1 面元x1和x3通過(guò)中間表面x2的輻射熱交換Fig.1 Radiation heat exchange between surface element x1 and x3 through x2
定向半球反射率r1
2 定義為由面元dx1發(fā)出經(jīng)過(guò)dx2反射到半球空間的反射率。
對(duì)于不透明體,由基爾霍夫輻射定律得到:
式中,e1
2 是由dx2向dx1方向的輻射通量與相同面積、相同溫度的黑體在相同方向的輻射通量的比值,稱為定向有效發(fā)射率。
任意形狀,開(kāi)口面積為A的等溫腔如圖2 所示,觀察方向?yàn)镈,則由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ā)射率的簡(jiǎn)化積分方程:
根據(jù)相互性原理,上式簡(jiǎn)化為:
式中,ρD0表示由方向D入射到腔體x0處,并最終反射到半球空間的定向半球反射率。將式( 6) 、(7) 、(11) 帶入式(12) 中,可得:
式(13) 的解可以表示為:
式中,fi表示從方向D入射到x0處,經(jīng)過(guò)i次反射射出腔體的輻射量。
在過(guò)去的30 多年間,研究人員已經(jīng)成功地對(duì)球體、圓錐、圓柱及組合形狀等腔的定向有效發(fā)射率進(jìn)行蒙特-卡羅計(jì)算,雖然在具體計(jì)算過(guò)程中采用的算法不同,但通過(guò)對(duì)每束光線每次反射出腔口徑的輻射量進(jìn)行求和,用于計(jì)算腔的半球反射率,然后由基爾霍夫輻射定律得到腔的定向有效發(fā)射率的物理思想是基本相同的。
光線追跡算法是計(jì)算腔體有效發(fā)射率的核心。在國(guó)際上曾經(jīng)有正向追跡和逆向追跡兩種方法,二者的區(qū)別是起始發(fā)出光線的位置不同。正向追跡是光線從腔內(nèi)壁一點(diǎn)發(fā)出,然后對(duì)光線進(jìn)行追跡,逆向追跡是光線由腔外觀察點(diǎn)發(fā)出,當(dāng)光線入射到腔內(nèi)部后對(duì)光線追跡。由于逆向追跡的思想在計(jì)算定向有效發(fā)射率時(shí)很方便,所以到后期一般都采用逆向追跡的方法。
不等溫腔的有效發(fā)射率只是在等溫情況下做了一些修正,這部分討論均是針對(duì)等溫腔,不等溫情況將在后面討論。假設(shè)腔內(nèi)表面是光學(xué)性質(zhì)均勻的漫-鏡反射體,反射率ρ=ρd+ρs( ρd為漫反射率,ρs為鏡面發(fā)射率) 。
逆向光線追跡是光線由腔外觀察點(diǎn)發(fā)出,入射到腔內(nèi)后對(duì)光線追跡。該方法可以方便地計(jì)算定向有效發(fā)射率,經(jīng)過(guò)積分運(yùn)算后可以得到平均垂直有效發(fā)射率和半球有效發(fā)射率。
逆向光線追跡的算法一般分為以下4 個(gè)步驟:
(1) 假設(shè)一束光線從D方向入射至腔的ξ0點(diǎn),在ξ0處一部分光線被吸收,一部分光線被反射。由偽隨機(jī)數(shù)產(chǎn)生器生成( 0,1) 之間的隨機(jī)數(shù)η,如果η <D,則腔體表面發(fā)生漫反射,否則發(fā)生鏡面反射。
如果光線第一次和腔壁相互作用發(fā)生鏡面反射,則反射光線的方向可以通過(guò)下式來(lái)確定:
式中,ωi,ωr,n 分別表示入射光矢量、反射光矢量以及腔壁ξ0的法線方向。
漫反射情況稍微復(fù)雜一些,近似認(rèn)為純漫反射的光能量均勻分布在半球空間中,所以在半球空間隨機(jī)均勻地產(chǎn)生一個(gè)方向代表該束光的漫射方向,如果反射光直接射出腔體,則結(jié)束該束光的追跡,否則繼續(xù)對(duì)該束光線追跡。
在半球空間中同等概率地產(chǎn)生漫射方向可以歸結(jié)為在半球上均勻地產(chǎn)生一點(diǎn),該點(diǎn)的坐標(biāo)就可以作為漫射的方向矢量。Marsaglia[13]曾全面系統(tǒng)地總結(jié)了關(guān)于在球面上均勻地產(chǎn)生一點(diǎn)的問(wèn)題,他不僅對(duì)算法的合理性、可行性予以考慮,還從計(jì)算速度、收斂特性上進(jìn)行了討論。通常采用的算法也是最容易理解的算法,即是在( -1,1)之間隨機(jī)產(chǎn)生3 個(gè)數(shù)時(shí),球面上隨機(jī)點(diǎn)的坐標(biāo)為:
Marsaglia 提出了一種新的方法,基本思想是在球體z軸上任意選取一點(diǎn),過(guò)該點(diǎn)作平行于xy面的平面,該平面與半徑為1 的球的交線是半徑為的圓周。在圓周上隨機(jī)產(chǎn)生一點(diǎn)。具體算法是產(chǎn)生偽隨機(jī)數(shù)V1、V2∈( -1,1) ,當(dāng)S=時(shí),球面隨機(jī)點(diǎn)的坐標(biāo)為( 2V1( 1 -由于該算法與以往相比減少了偽隨機(jī)數(shù)的個(gè)數(shù)和平方根的運(yùn)算次數(shù),計(jì)算速度提高了約2 倍。
(2) 追跡第一次反射的光線與腔壁作用點(diǎn)的坐標(biāo)。目前研究人員基本采用相同的算法追跡光線與腔壁作用點(diǎn)的坐標(biāo)。
式中:Φ(x,y,z) =0 為腔表面方程,ξ0為發(fā)出光線的位置坐標(biāo),ξ 為作用點(diǎn)位置坐標(biāo),ω 為反射方向矢量,t為系數(shù)。
(3) 光線第二次和腔壁發(fā)生相互作用,重復(fù)(1) ~(2) 的分析過(guò)程,直到光線射出腔口徑或者經(jīng)多次反射腔內(nèi)的光能量減小到可以忽略為止。( 假設(shè)入射光束總能量是1,多次反射能量衰減為10-5或10-6可結(jié)束追跡) 。
(4) 結(jié)束本條光束的追跡,進(jìn)行下一束光的追跡。為了提高蒙特-卡羅法的計(jì)算精度,對(duì)同一條件入射的情況,一般需要追跡105~107條光束。
根據(jù)上述逆向光束追跡過(guò)程,可以推導(dǎo)出黑體空腔定向有效發(fā)射率的計(jì)算方程。Sapritsky[9]給出的方程如下:
Sapritsky 總共對(duì)N條光束進(jìn)行追跡,其中第i條光束經(jīng)過(guò)M次反射后射出腔體。F是漫反射角因子,具體表達(dá)式如下:
式中:Ω 表示ξ 與腔口徑所成的立體角; θξ是腔壁ξ 處的法線與立體角dΩ 軸線之間的夾角。F是漫反射光通過(guò)腔口徑射出腔體的部分。
由此可以看出,Sapritsky 認(rèn)為每一條光束只要經(jīng)歷一次漫反射,能量就減小1 -F( ξ) 倍。這對(duì)鏡面反射也是成立的,如果光束經(jīng)鏡面反射剛好射出腔體,則F( ξ) =1,否則F( ξ) =0。
Prokhorov[11]采用更簡(jiǎn)單的算法:
上式各物理量的含義同式( 18) ,只是Prokhorov 將漫反射和鏡面反射同等對(duì)待,只是在確定兩者的反射方向上有所區(qū)別。
正向光線追跡可以方便地計(jì)算具有漫-鏡反射表面腔的半球有效發(fā)射率。Heinisch[14]曾采用正向追跡方法較精確計(jì)算漫反射等溫圓錐形腔的半球有效發(fā)射率。這里將計(jì)算推廣到任意形狀、具有漫-鏡反射表面的等溫腔:
(1) 腔壁發(fā)出光束的起始點(diǎn)坐標(biāo)ξ 由偽隨機(jī)數(shù)確定。理想情況下,在趨于無(wú)窮次的光束追跡過(guò)程中,起始發(fā)出光線點(diǎn)的位置應(yīng)均勻地分布在整個(gè)腔壁上。
(2) 從ξ 發(fā)出的光能量分為兩部分,一部分能量E·Fi直接從腔口射出,另一部分能量E(1 -Fi) 在腔內(nèi)沿隨機(jī)方向(Rθ、Rφ) 傳播。傳至下一點(diǎn)時(shí),光束能量將被隨機(jī)地吸收或反射,這取決于0 ~1 的隨機(jī)數(shù)Rα,當(dāng)Rα≤ε 時(shí),能量被吸收,否則光束被反射。反射的能量E(1 -Fi) 再次分為兩部分,分別為E(1 -Fi)Fi1和E(1 -Fi) ·(1 -Fi1) 。重復(fù)上面的過(guò)程直至能量在某一點(diǎn)處被吸收。
(3) 選另一發(fā)光點(diǎn),重復(fù)上面的過(guò)程。
半球有效發(fā)射率的計(jì)算方程:
式中:Gi=(1 -Fi)Fi1+(1 -Fi) (1 -Fi1)Fi2+…+(1 -Fi) ( 1 -Fi1)Fi2…( 1 -Fim-1)Fim;m=1,2…表示能量在第m次反射之后被吸收;Fi為第i個(gè)發(fā)光點(diǎn)所在微元對(duì)腔口的角系數(shù);Fim為第i個(gè)發(fā)光點(diǎn)發(fā)出的能量在第m次反射時(shí)所在微元對(duì)腔口的角系數(shù)。
上面講述了采用逆向追跡的方法計(jì)算黑體空腔定向有效發(fā)射率的具體算法。垂直有效發(fā)射率是光束垂直于孔徑的特殊定向有效發(fā)射率,它的算法和定向有效發(fā)射率相同。
積分有效發(fā)射率εe(Rd,Hd) 是半徑為Rd的圓形探測(cè)器共軸放置在與腔口徑相距Hd的位置時(shí),由探測(cè)器測(cè)量得到的腔體發(fā)射率。它是計(jì)算平均垂直有效發(fā)射率和半球有效發(fā)射率的基礎(chǔ)。
圖3 計(jì)算積分有效發(fā)射率簡(jiǎn)圖Fig.3 Diagram of calculating integrated effective emissivity
Prokhorov[11]給出了積分有效發(fā)射率的計(jì)算公式,如圖3 所示,由探測(cè)器接收到的總輻射通量表達(dá)式如下:
腔口徑面元dSa出射的輻射沿ψ 向入射到探測(cè)器dSd,L是P點(diǎn)沿PQ方向的輻亮度。如果繼續(xù)采用逆向光線追跡的方法,認(rèn)為光線由均勻分布在探測(cè)器圓面上的Qi發(fā)出,經(jīng)過(guò)均勻分布在腔口徑上的Pi點(diǎn)入射到腔體內(nèi)部,之后的光線追跡過(guò)程和上面介紹的定向有效發(fā)射率完全相同。那么上面的積分運(yùn)算可以變?yōu)榍蠛瓦\(yùn)算:
將黑體空腔換為面積為Sa的理想黑體,在相同條件下,探測(cè)器接收的輻射通量:
式中,F(xiàn)a-d為理想黑體對(duì)探測(cè)器圓面Sd的角系數(shù)。
此時(shí)黑體空腔的積分有效發(fā)射率:
Sapritsky[9]提出了不等溫腔的定向有效發(fā)射率,它是在等溫情況下加上修正因子,如下所示:
式中:Tξ為腔壁ξ 點(diǎn)的實(shí)際溫度,T0為參考恒定溫度。在逆向光線追跡的過(guò)程中,不等溫修正項(xiàng)如下:
式中:Tk為光線第k次在腔壁上反射時(shí)腔壁作用點(diǎn)處的實(shí)際溫度;Le( λ,Tk) 為溫度為T(mén)k,波長(zhǎng)為λ 的光譜輻亮度。
本文綜述了運(yùn)用蒙特-卡羅方法計(jì)算黑體空腔有效發(fā)射率的理論基礎(chǔ)和相應(yīng)的光線追跡算法。不難看出,該方法能夠方便、精確地計(jì)算黑體空腔的有效發(fā)射率。在計(jì)算具有漫反射表面的圓柱形腔的積分有效發(fā)射率時(shí)( 包括半球有效發(fā)射率和平均垂直有效發(fā)射率) ,得到的不確定度為0.000 1 ~0.000 2,與一些研究人員采用相關(guān)理論計(jì)算方法得到的不確定度大致相同[18-19]。大部分采用理論計(jì)算或者實(shí)驗(yàn)測(cè)試的方法得到的腔體發(fā)射率結(jié)果與蒙特-卡羅法得到的結(jié)果均符合得很好。目前,蒙特-卡羅方法基本上已經(jīng)完全取代理論計(jì)算的方法,成為獲得黑體空腔有效發(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é)報(bào),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ì)算的蒙特卡羅模型[J].清華大學(xué)學(xué)報(bào),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ā)射率計(jì)算中的應(yīng)用[J].宇航計(jì)測(cè)技術(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.