單軸自旋平面對固定平面的旋轉(zhuǎn)角系數(shù)研究

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

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

0 引言

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

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

1 計算模型及驗證

1．1 計算模型

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

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

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

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

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

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

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

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

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

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

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

的天頂角和圓周角/°;

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

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

1．2 模型驗證

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

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

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

2 分析與討論

2．1 平移矢量

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

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

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

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

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

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

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

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

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

式中

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

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

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

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

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

3 結(jié) 論

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

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

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

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

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

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

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

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

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

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

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

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

［9］高亞鋒，湯小敏，許孟楠，等． 城市地面對天空熱輻射角系數(shù)的簡化求解［J］． 太陽能學(xué)報，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ù)值計算［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．