張旭升，郭 亮，黃 勇，劉春龍

(中國(guó)科學(xué)院 長(zhǎng)春光學(xué)精密機(jī)械與物理研究所，吉林 長(zhǎng)春 130033)

0 引言

20 世紀(jì)20 年代，表面輻射換熱領(lǐng)域提出了用于表征幾何形狀、大小及其相對(duì)位置的物理概念—角系數(shù)［1］，并被廣泛應(yīng)用于高溫爐內(nèi)換熱［2－3］、航天器熱控系統(tǒng)設(shè)計(jì)［4］、輻射采暖及散熱［5－6］等諸多方面，在后續(xù)發(fā)展中衍生出輻射傳遞因子等類(lèi)似概念。作為表面輻射換熱領(lǐng)域的重要概念，國(guó)內(nèi)外學(xué)者已對(duì)角系數(shù)進(jìn)行了相關(guān)研究。張濤等［7］提出了計(jì)算角系數(shù)的能束均勻分步法，與有限元法結(jié)合后可用于求解復(fù)雜結(jié)構(gòu)的角系數(shù)。王天鵬等［8］建立了計(jì)算角系數(shù)的高精度有限體積模型，并將其應(yīng)用于計(jì)算地面建筑結(jié)構(gòu)的輻射換熱中。高亞峰等［9］采用等分圓弧射線(xiàn)掃描法實(shí)現(xiàn)了復(fù)雜建筑群內(nèi)地面對(duì)天空角系數(shù)的快速簡(jiǎn)化求解。張偉清等［4］提出了基于角系數(shù)計(jì)算漫反射表面輻射傳遞因子的新方法，且計(jì)算誤差主要源于角系數(shù)。M． R． Vuji i 等［10］采用蒙特卡洛法計(jì)算了平行矩形平面及各網(wǎng)格單元間的角系數(shù)，指出復(fù)雜結(jié)構(gòu)輻射換熱計(jì)算中，過(guò)分細(xì)化網(wǎng)格必將造成角系數(shù)計(jì)算誤差的增大。同樣，M．Mirhosseini 等［11］也采用蒙特卡洛法計(jì)算了圓柱面與相切矩形微元面間的角系數(shù)，文中著重分析了網(wǎng)格離散方案和單元光線(xiàn)數(shù)對(duì)計(jì)算精度和計(jì)算效率的影響。J．M． Cabeza －Lainez 等［12］應(yīng)用幾何光學(xué)原理推導(dǎo)了球冠面、圓錐面、橢圓面、拋物面等二次曲面間的角系數(shù)計(jì)算式。A． Bahadori 等［13］發(fā)展了一種計(jì)算矩形平面與微元面間角系數(shù)的簡(jiǎn)化方法，所得公式可用于快速解決工程領(lǐng)域的輻射換熱問(wèn)題。

綜上所述，各文獻(xiàn)的研究重點(diǎn)均集中在固定平面角系數(shù)計(jì)算或復(fù)雜結(jié)構(gòu)角系數(shù)計(jì)算簡(jiǎn)化等方面，并未對(duì)轉(zhuǎn)動(dòng)部件角系數(shù)的高精度計(jì)算進(jìn)行相關(guān)研究。鑒于此，本文將采用計(jì)算精度高、物理概念清晰、對(duì)復(fù)雜結(jié)構(gòu)適應(yīng)性強(qiáng)［6，14］的蒙特卡洛法建立單軸自旋平面對(duì)固定平面的旋轉(zhuǎn)角系數(shù)計(jì)算模型。通過(guò)與典型自旋平面位置下解析式計(jì)算值的比對(duì)，驗(yàn)證物理模型及算法的正確性。在此基礎(chǔ)上，分析旋轉(zhuǎn)角度、初始相對(duì)位置等因素對(duì)角系數(shù)大小的影響，探討發(fā)射光線(xiàn)數(shù)目和初始相對(duì)位置與統(tǒng)計(jì)標(biāo)準(zhǔn)差的關(guān)聯(lián)性。

1 計(jì)算模型及驗(yàn)證

1．1 計(jì)算模型

角系數(shù)是用于表征物體表面結(jié)構(gòu)特性對(duì)輻射換熱影響程度的無(wú)量綱因子，計(jì)算時(shí)需滿(mǎn)足以下限制條件:

(1)表面具有漫射特性且不透明;

(2)表面輻射熱流密度在空間均勻分布［11，14］。

在上述條件下，建立了旋轉(zhuǎn)角系數(shù)的物理模型及系統(tǒng)坐標(biāo)系，如圖1 所示。初始時(shí)刻，長(zhǎng)L ×寬W的兩矩形平面平行且中心正對(duì)，垂直距離為H。然后，單軸自旋平面沿著y 軸旋轉(zhuǎn)，旋轉(zhuǎn)方向與y 軸遵循左手螺旋定則，且90°為旋轉(zhuǎn)角度的極限位置，如圖1 中虛線(xiàn)框所示。

圖1 旋轉(zhuǎn)角系數(shù)物理模型及系統(tǒng)坐標(biāo)系

在系統(tǒng)坐標(biāo)系下，單軸自旋平面光線(xiàn)發(fā)射位置概率模型計(jì)算式如下

式中 α——單軸自旋平面旋轉(zhuǎn)角/°;

Rx、Ry——［0，1］區(qū)間內(nèi)獨(dú)立均布隨機(jī)數(shù)。

由于光線(xiàn)發(fā)射點(diǎn)位于單軸自旋平面上，故每次循環(huán)中滿(mǎn)足關(guān)系式Rx=Rz。

對(duì)于漫發(fā)射單軸自旋平面，光線(xiàn)發(fā)射方向概率模型計(jì)算式如下

式中 θ'、φ'——自旋平面局部坐標(biāo)系下發(fā)射光線(xiàn)

的天頂角和圓周角/°;

Rθ、Rφ——［0，1］區(qū)間內(nèi)獨(dú)立均布隨機(jī)數(shù)。根據(jù)坐標(biāo)變換關(guān)系可推導(dǎo)出如下角度變換關(guān)系式

再結(jié)合式(2)和反三角函數(shù)關(guān)系，即可求得系統(tǒng)坐標(biāo)系下發(fā)射光線(xiàn)的天頂角θ 和圓周角φ。至此，已具備蒙特卡洛法計(jì)算旋轉(zhuǎn)角系數(shù)的定解條件。為了便于對(duì)旋轉(zhuǎn)角系數(shù)的分析討論，規(guī)定L =W =H，均取無(wú)量綱單位長(zhǎng)度1．0。

1．2 模型驗(yàn)證

為了驗(yàn)證旋轉(zhuǎn)角系數(shù)物理模型及算法的正確性，將典型自旋平面位置下角系數(shù)與J．R．Howell 解析式計(jì)算值［15］進(jìn)行比對(duì)，如表1 所示。其中，本文蒙特卡洛法計(jì)算中采用單次循環(huán)、發(fā)射光線(xiàn)數(shù)目為1．0 ×108根。此時(shí)，角系數(shù)計(jì)算值與發(fā)射光線(xiàn)數(shù)目已呈現(xiàn)無(wú)關(guān)性。

表1 典型自旋平面位置下角系數(shù)計(jì)算值

分析上表可知，無(wú)論自旋平面處于平行或垂直的典型位置下，本文角系數(shù)計(jì)算值均與文獻(xiàn)［15］解析值大體一致，僅僅在小數(shù)點(diǎn)后5 位及以上存在微小差異，這也是基于概率統(tǒng)計(jì)的蒙特卡洛法所必然產(chǎn)生的隨機(jī)誤差［16］。

2 分析與討論

2．1 平移矢量

為了分析旋轉(zhuǎn)角度、初始相對(duì)位置等因素對(duì)角系數(shù)大小的影響，文中將固定平面沿矢量s 進(jìn)行平移。依據(jù)平面中心點(diǎn)及邊緣的相對(duì)位置關(guān)系，規(guī)劃出0．25、0．5、0．75 和1．0 等4 種平移距離和3 種平移方向，并將平移距離為0．0 的情況定義為基準(zhǔn)，據(jù)此所構(gòu)造算例的具體內(nèi)容見(jiàn)表2，角度β 為平移矢量s 的方位角。

表2 算例中固定平面的平移矢量

在各算例的模擬計(jì)算中，發(fā)射光線(xiàn)數(shù)目均取1.0 ×107根且為單次循環(huán)。角系數(shù)計(jì)算值如圖2、圖3、圖4 所示，圖中Δ(+x)=0．5 表示固定平面沿+x 方向平移0．5 個(gè)距離，其他以此類(lèi)推。

圖2 算例Ⅰ條件下的角系數(shù)計(jì)算值

圖3 算例Ⅱ條件下的角系數(shù)計(jì)算值

圖4 算例Ⅲ條件下的角系數(shù)計(jì)算值

分析基準(zhǔn)算例可知，隨著旋轉(zhuǎn)角的增大，角系數(shù)X (α，s)約從0．200 單調(diào)遞減至0．034，且在旋轉(zhuǎn)角大于45°后角系數(shù)呈線(xiàn)性下降趨勢(shì)。圖2 中給出固定平面沿+x 方向移動(dòng)時(shí)旋轉(zhuǎn)角度和初始位置對(duì)角系數(shù)的影響，即算例Ⅰ。隨著平移距離Δ(+x)的增大，角系數(shù)變化范圍逐漸減小，相同旋轉(zhuǎn)角度對(duì)應(yīng)的角系數(shù)明顯下降;當(dāng)Δ(+x)≥0．5 時(shí)，大角度角系數(shù)出現(xiàn)“零值現(xiàn)象”，該現(xiàn)象也從側(cè)面驗(yàn)證了本文物理模型及算法的正確性。圖3 中給出固定平面沿+y 方向移動(dòng)時(shí)旋轉(zhuǎn)角度和初始位置對(duì)角系數(shù)的影響，即算例Ⅱ。隨著平移距離Δ(+y)的增大，角系數(shù)變化范圍亦逐漸減小，但未出現(xiàn)“零值現(xiàn)象”;數(shù)據(jù)歸一化后發(fā)現(xiàn)，不同Δ(+y)的角系數(shù)變化趨勢(shì)與基準(zhǔn)算例大體一致、差異微小。圖4 中給出固定平面沿－x 方向移動(dòng)時(shí)旋轉(zhuǎn)角度和初始位置對(duì)角系數(shù)的影響，即算例Ⅲ。隨著平移距離Δ(－x)的增大，角系數(shù)變化范圍顯著減小，當(dāng)Δ(－x)=1．0 時(shí)，角系數(shù)僅分布在0．068 ～0．097 之間;而當(dāng)Δ(－x)≠0．0 時(shí)，存在非零旋轉(zhuǎn)角αmax使得角系數(shù)取得極大值，且Δ(－ x)與αmax呈正相關(guān)，如圖4 點(diǎn)劃線(xiàn)所示。

2．2 統(tǒng)計(jì)標(biāo)準(zhǔn)差

在蒙特卡洛法中，偽隨機(jī)數(shù)程序性能優(yōu)劣將直接關(guān)系到計(jì)算結(jié)果的精度［16］。為保證旋轉(zhuǎn)角系數(shù)在允許誤差范圍內(nèi)，文中將對(duì)不同旋轉(zhuǎn)角和初始相對(duì)位置下的統(tǒng)計(jì)標(biāo)準(zhǔn)差進(jìn)行分析。定義旋轉(zhuǎn)角為α?xí)r，旋轉(zhuǎn)角系數(shù)的統(tǒng)計(jì)標(biāo)準(zhǔn)差σ(α，k)計(jì)算式如下

式中

k——蒙特卡洛法中旋轉(zhuǎn)角系數(shù)的統(tǒng)計(jì)循環(huán)次數(shù);

Xn(α，s) ——第n 次循環(huán)中旋轉(zhuǎn)角為α 時(shí)的角系數(shù)。在統(tǒng)計(jì)標(biāo)準(zhǔn)差分析計(jì)算中，循環(huán)次數(shù)均為10 次。

圖5 基準(zhǔn)算例標(biāo)準(zhǔn)差與發(fā)射光線(xiàn)數(shù)目關(guān)系

圖5 給出發(fā)射光線(xiàn)數(shù)目分別為1．0 ×105、1．0 ×106和1．0 ×107時(shí)，基準(zhǔn)算例標(biāo)準(zhǔn)差隨發(fā)射光線(xiàn)數(shù)目和旋轉(zhuǎn)角度的變化關(guān)系。隨著發(fā)射光線(xiàn)數(shù)目的增大，相同旋轉(zhuǎn)角對(duì)應(yīng)的統(tǒng)計(jì)標(biāo)準(zhǔn)差明顯下降，但并非呈現(xiàn)線(xiàn)性關(guān)系;同時(shí)，在規(guī)定旋轉(zhuǎn)角范圍內(nèi)，標(biāo)準(zhǔn)差曲線(xiàn)變得更加均勻平滑，當(dāng)發(fā)射光線(xiàn)數(shù)目為1．0 ×107時(shí)，最大標(biāo)準(zhǔn)差也僅為1．957 ×10－4。因此，在分析初始相對(duì)位置對(duì)統(tǒng)計(jì)標(biāo)準(zhǔn)差的影響時(shí)，發(fā)射光線(xiàn)數(shù)目均取1．0 ×107根，如圖6 所示。研究表明，隨著平移距離Δ(－x)的增大，相同旋轉(zhuǎn)角的標(biāo)準(zhǔn)差未呈現(xiàn)出明顯變化規(guī)律，但整體波動(dòng)范圍略有減小，即角系數(shù)的角度均勻性有所提升。即:增大發(fā)射光線(xiàn)數(shù)目可同時(shí)提高角系數(shù)的計(jì)算精度和角度均勻性，而增大平移距離僅對(duì)角度均勻性略有改善。

圖6 算例Ⅲ初始相對(duì)位置對(duì)標(biāo)準(zhǔn)差的影響

3 結(jié) 論

通過(guò)上述計(jì)算分析，驗(yàn)證了蒙特卡洛法求解單軸自旋平面對(duì)固定平面旋轉(zhuǎn)角系數(shù)物理模型及算法的正確性。在此基礎(chǔ)上，分析了旋轉(zhuǎn)角度、初始相對(duì)位置、發(fā)射光線(xiàn)數(shù)目等因素對(duì)角系數(shù)大小和統(tǒng)計(jì)標(biāo)準(zhǔn)差的影響。研究表明

(1)當(dāng)固定平面平移矢量s 不存在－x 方向分量時(shí)，角系數(shù)隨旋轉(zhuǎn)角的增大而減小，且平移距離Δs 越大，角系數(shù)變化范圍越小，當(dāng)Δ(+ x)≥0． 5時(shí)，大角度角系數(shù)出現(xiàn)“零值現(xiàn)象”;

(2)當(dāng)固定平面平移矢量s 存在－x 方向分量時(shí)，對(duì)于任意非零平移距離Δs，存在非零旋轉(zhuǎn)角αmax使得角系數(shù)取得極大值，且Δs 與αmax呈正相關(guān);

(3)增大發(fā)射光線(xiàn)數(shù)目可同時(shí)提高角系數(shù)的計(jì)算精度和角度均勻性，而增大平移距離僅對(duì)角度均勻性略有改善。

［1］余其錚．輻射換熱原理［M］．哈爾濱:哈爾濱工業(yè)大學(xué)出版社，2000．

［2］岳雷．氧氣轉(zhuǎn)爐汽化冷卻煙道傳熱計(jì)算［J］．節(jié)能技術(shù)，2012，30(3):245 －248．

［3］陳朝松，張樹(shù)林，劉平元，等． 優(yōu)化壁溫計(jì)算模型及其在電站鍋爐壁溫在線(xiàn)監(jiān)測(cè)中的應(yīng)用［J］．動(dòng)力工程，2009，29(9):818 －822．

［4］張偉清，宣益民，韓玉閣．單元表面間輻射傳遞系數(shù)的新型計(jì)算方法［J］．宇航學(xué)報(bào)，2005，26(1):77 －80．

［5］王麗，劉志錚，王岳人．低溫地板輻射采暖埋管結(jié)構(gòu)尺寸設(shè)計(jì)實(shí)驗(yàn)研究［J］．節(jié)能技術(shù)，2010，28(1):11 －14．

［6］李景文．開(kāi)式直筒輻射管輻射換熱性能的數(shù)值計(jì)算［J］．節(jié)能技術(shù)，2013，31(6):512 －515．

［7］張濤，孫冰． 復(fù)雜結(jié)構(gòu)角系數(shù)計(jì)算方法［J］． 航空動(dòng)力學(xué)報(bào)，2009，24(4):753 －759．

［8］王天鵬，王良璧． 一種角系數(shù)的數(shù)值算法及在室內(nèi)壁面與窗口間輻射換熱中的應(yīng)用［J］． 蘭州交通大學(xué)學(xué)報(bào)，2012，31(6):126 －129．

［9］高亞鋒，湯小敏，許孟楠，等． 城市地面對(duì)天空熱輻射角系數(shù)的簡(jiǎn)化求解［J］． 太陽(yáng)能學(xué)報(bào)，2012，33(4):635 －639．

［10］M．R． Vuji i ，N．P． Lavery，S．G．R． Brown．View Factor Calculation Using the Monte Carlo Method and Numerical Sensitivity［J］． Communications in Numerical Methods in Engineering，2006，22(3):197 －203．

［11］M． Mirhosseini，A． Saboonchi． View Factor Calculation Using the Monte Carlo Method for a 3D Strip Element to Circular Cylinder［J］．International Communications in Heat and Mass Transfer，2011(38):821 －826．

［12］Jose M． Cabeza － Lainez，Jesus A． Pulido － Arcas．New Configuration Factors for Curved Surfaces［J］． Journal of Quantitative Spectroscopy ＆ Radiative Transfer，2013(117):70 －80．

［13］Alireza Bahadori，Sohrab Zendehboudi，Gholamreza Zahedi，et al．． Estimation of Configuration Factor for Radiation from a Rectangle to a Parallel Small Element of Surface Lying on the Perpendicular to a Corner of the Radiator［J］．Chemical Engineering Communications，2014，201(3):287 －299．

［14］談和平，夏新林，劉林華，阮立明．紅外輻射特性與傳輸?shù)臄?shù)值計(jì)算［M］． 哈爾濱:哈爾濱工業(yè)大學(xué)出版社，2006．

［15］J． R． Howell． A Catalog of Radiation Heat Transfer Configuration Factors (Version 4)［M］． New York:McGraw Hill，2010．

［16］Yong Shuai，Hao－Chun Zhang，He －Ping Tan．Radiation Symmetry Test and Uncertainty Analysis of Monte Carlo Method Based on Radiative Exchange Factor ［J］． Journal of Quantitative Spectroscopy ＆ Radiative Transfer，2008(109):1281 －1296．