基于蒙特卡洛方法的紅外光源輻照特性模擬

宋偉浩,周賢文,張 薇

(上海海事大學(xué)航運(yùn)技術(shù)與控制工程交通行業(yè)重點(diǎn)實(shí)驗(yàn)室,上海 201306)

1 引 言

產(chǎn)品在運(yùn)輸、貯存和使用過(guò)程中,溫度對(duì)產(chǎn)品的使用壽命及可靠性有著重要影響[1-2]。為確保產(chǎn)品的可靠性和使用的安全性,有必要進(jìn)行溫度可靠性環(huán)境試驗(yàn)。以紅外燈陣列作為輻射熱源進(jìn)行溫度測(cè)試,通過(guò)優(yōu)化燈陣可有效提高受熱平面溫度分布均勻程度,進(jìn)而提高溫度環(huán)境試驗(yàn)的可信度與準(zhǔn)確性。燈陣優(yōu)化工作的基礎(chǔ)是建立紅外燈的輻射照度分布模型。

Ziemke等采用有限元法計(jì)算了紅外燈陣的輻射照度分布[3]。劉守文等利用蒙特卡洛方法對(duì)燈管半側(cè)鍍有反射膜的管狀紅外燈進(jìn)行了建模,實(shí)驗(yàn)結(jié)果顯示所建立模型具有較高精度[4]。楊曉寧等采用蒙特卡洛方法模擬了無(wú)反射涂層的管狀紅外燈輻照分布情況,通過(guò)優(yōu)化燈陣改善了受熱面上熱流分布均勻性[5]。楊國(guó)巍等學(xué)者用蒙特卡洛模擬了管狀紅外燈熱流分布,建立了熱流分布數(shù)據(jù)庫(kù),最后采用遺傳算法優(yōu)化了紅外燈陣[6]。

張坤、陳海清等針對(duì)帶有拋物面型反射器的紅外燈泡輻射特性進(jìn)行了模擬仿真研究[7],未考慮對(duì)光線方向抽樣半球空間的方位變化。

許多學(xué)者多以無(wú)反射涂層或在燈管半圓周涂有反射層的管狀紅外燈為研究對(duì)象,采用蒙特卡洛方法模擬紅外單燈的熱流分布,對(duì)于帶有拋物面反射器紅外燈的建模與研究工作較少。

2 點(diǎn)光源的蒙特卡洛模擬

利用蒙特卡洛方法對(duì)紅外燈輻照特性進(jìn)行模擬,要按照一定的概率模型對(duì)出射光線進(jìn)行重復(fù)的隨機(jī)抽樣生成光線的相關(guān)參數(shù)。利用計(jì)算機(jī)生成服從特定概率分布的關(guān)于光線參數(shù)的隨機(jī)數(shù),是對(duì)紅外燈輻照特性進(jìn)行蒙特卡洛模擬的基礎(chǔ)。

在點(diǎn)光源的理想模型中,其出射光線均勻散布在以該點(diǎn)光源為球心的球體空間內(nèi),對(duì)光線的方向向量隨機(jī)抽樣,可唯一確定某條隨機(jī)光線[8]。設(shè)輻射線的天頂角為θ;圓周角為φ;與球面交點(diǎn)為P;則向量OP即為抽樣光線的方向向量,示意圖如圖1所示。對(duì)點(diǎn)光源的出射光線進(jìn)行隨機(jī)抽樣,可等同于對(duì)單位球面上隨機(jī)點(diǎn)P抽樣。

圖1 點(diǎn)光源光線隨機(jī)抽樣示意圖

設(shè)P點(diǎn)位置服從球面上的均勻分布,且圓周角φ和天頂角θ相互獨(dú)立,取服從(0,1)均勻分布的隨機(jī)數(shù)rφ與rθ,可表示出點(diǎn)光源出射光線方向向量nx,ny,nz,如式(1):

(1)

3 擴(kuò)展光源的蒙特卡洛模擬

擴(kuò)展光源可以視作點(diǎn)光源的集合體,但發(fā)光特點(diǎn)又與點(diǎn)光源有所區(qū)別,需對(duì)其光線出射位置與方向進(jìn)行抽樣。為簡(jiǎn)化分析,將鎢絲紅外燈燈絲簡(jiǎn)化為圓柱體擴(kuò)展光源,僅由圓柱側(cè)面向外出射光線。假設(shè):紅外燈圓柱發(fā)光體為漫灰體,各處溫度均勻一致；每條出射紅外光攜帶的能量相等,與出射位置、角度以及波長(zhǎng)無(wú)關(guān)。

對(duì)圓柱面上發(fā)光點(diǎn)的位置抽樣,假定圓柱發(fā)光體半徑為r,高度為h,軸線與z軸重合,下表面平面方程為z=h0。對(duì)光線出射點(diǎn)隨機(jī)抽樣的示意圖如圖2所示。隨機(jī)發(fā)光點(diǎn)M位置可由圓周角α和豎坐標(biāo)z表示。

圖2 圓柱面發(fā)光點(diǎn)坐標(biāo)抽樣示意圖

M點(diǎn)服從圓柱面上均勻分布,其落在陰影部分面積的聯(lián)合概率分布函數(shù),如式(2)所示:

Fα,z=α/2π·z-h0/h

(2)

0≤α≤2π,h0≤θ≤h0+h

由式(2)可知,α和z相互獨(dú)立,故可取服從(0,1)均勻分布的隨機(jī)數(shù)r1和r2分離變量,α和z可以由式(3)表示:

(3)

進(jìn)而點(diǎn)M(x0,y0,z0)坐標(biāo)可由式(4)描述。

(4)

圓柱體擴(kuò)展光源的出射光線散布于以圓柱在M點(diǎn)外向法線為旋轉(zhuǎn)對(duì)稱軸的半球空間內(nèi),且出射方向遵循蘭貝特定律。半球空間隨抽樣點(diǎn)M位置變動(dòng)而變動(dòng)。為確定M點(diǎn)出射光線的方向向量,需以M點(diǎn)為原點(diǎn)O′,建立輔助三維空間直角坐標(biāo)系O′-x′y′z′。輔助坐標(biāo)系O′-x′y′z′的建立示意圖如圖3。

圖3 輔助坐標(biāo)系示意圖

設(shè)在坐標(biāo)系O-xyz下,圓柱的方程為F(x,y,z)=0,x′、y′、z′軸方向向量分別為a、b、c。以圓柱體在點(diǎn)M處的外法向量為z′軸方向向量c,可由式(5)表示:

c=k1,k2,k3

(5)

a=b×c

(6)

將a,b,c單位化,分別得到O-xyz坐標(biāo)系下的x′軸、y′軸和z′軸的單位方向向量:l1,l2,l3,m1,m2,m3和n1,n2,n3。

在O′-x′y′z′坐標(biāo)系內(nèi),對(duì)M點(diǎn)也即O′點(diǎn)出射光線進(jìn)行蒙特卡洛抽樣,示意圖如圖4所示。

圖4 圓柱面出射光線隨機(jī)抽樣示意圖

光線O′Q方向可由圓周角β與天頂角γ表示,且可按照蘭貝特定律處理。取服從(0,1)均勻分布的隨機(jī)數(shù)r3和r4,則β與γ可由式(7)表示:

(7)

進(jìn)而可得在O′-x′y′z′坐標(biāo)系內(nèi)光線方向向量(n′x,n′y,n′z),如式(8)所示:

(8)

為進(jìn)行后續(xù)的光線追跡,需將該方向向量按照式(9)轉(zhuǎn)化為坐標(biāo)系O-xyz下單位光線方向向量(i0,j0,k0):

(9)

已知隨機(jī)點(diǎn)M坐標(biāo)(x0,y0,z0)和光線在O-xyz坐標(biāo)系下的單位方向向量(i0,j0,k0),則光線在O-xyz坐標(biāo)系內(nèi)的直線參數(shù)方程可由式(10)表示:

(10)

至此,運(yùn)用蒙特卡洛方法完成了對(duì)圓柱發(fā)光體一條光線的隨機(jī)抽樣。為研究受熱平面上的輻射照度分布,需對(duì)大量光線進(jìn)行隨機(jī)抽樣,并對(duì)每條光線進(jìn)行追跡,獲取其在受熱平面上的落點(diǎn)坐標(biāo)。在受熱平面上劃分網(wǎng)格,統(tǒng)計(jì)各網(wǎng)格內(nèi)光線落點(diǎn)數(shù)量,進(jìn)一步根據(jù)單光線攜帶能量值和網(wǎng)格面積,計(jì)算得出紅外燈在受熱平面的輻射照度分布[7]。

對(duì)光線進(jìn)行追跡,需要根據(jù)所選紅外燈的光學(xué)結(jié)構(gòu)特點(diǎn)進(jìn)行解算。選用帶有玻璃反射器的紅外燈泡,其示意圖如圖5所示,其中圖5(a)是紅外燈的實(shí)物圖,圖5(b)是紅外燈等效的數(shù)學(xué)模型圖。

圖5 紅外燈實(shí)物圖與數(shù)學(xué)模型簡(jiǎn)圖

忽略玻璃燈罩對(duì)光路的影響與發(fā)光圓柱體對(duì)光線的遮擋作用。則圓柱發(fā)光體的出射光線會(huì)直接或者經(jīng)過(guò)反射器反射后抵達(dá)受熱平面。參考實(shí)際反射器形狀可將其等效為一旋轉(zhuǎn)拋物面,圓柱發(fā)光體的軸線與拋物面旋轉(zhuǎn)對(duì)稱軸重合,且圓柱幾何中心與和拋物面焦點(diǎn)相重合。開口向下且旋轉(zhuǎn)對(duì)稱軸為z軸的旋轉(zhuǎn)拋物面可由式(11)描述:

x2+y2+4Fz-h1+h2=0

(11)

z≥h1

式中,F為旋轉(zhuǎn)拋物面的焦距,m；h1為反射器下邊緣高度,m；h2為紅外燈的反射器的高度,m。

聯(lián)立光線所在直線方程和拋物面方程,可解出直線與拋物面的一個(gè)或者兩個(gè)交點(diǎn)。當(dāng)存在一個(gè)交點(diǎn)時(shí),該交點(diǎn)為光線所在射線的反向延長(zhǎng)線與拋物面的交點(diǎn),屬無(wú)效交點(diǎn),此時(shí)光線將直接投射到受熱平面。當(dāng)存在兩個(gè)交點(diǎn)時(shí),其中一個(gè)為無(wú)效交點(diǎn),另一個(gè)為有效交點(diǎn)。由發(fā)光點(diǎn)指向有效交點(diǎn)的向量與光線的方向向量夾角為零度,據(jù)此可篩選出入射點(diǎn)位置坐標(biāo)。

光線在反光罩上發(fā)生發(fā)射時(shí),對(duì)光線進(jìn)行追跡需要求解其反射光線。聯(lián)立式(10)和式(11)求解并經(jīng)過(guò)判別得入射點(diǎn)Q0坐標(biāo)。根據(jù)反射定律可解得反射光線的方向向量。

(12)

隨機(jī)光線可能在反射器內(nèi)經(jīng)過(guò)多次反射才能抵達(dá)受熱平面。當(dāng)光線經(jīng)過(guò)多次反射時(shí),可將求解的反射光線作為下次反射的入射光線。當(dāng)判定光線可直達(dá)受熱平面時(shí),聯(lián)立光線所在直線方程和接收平面方程,可解得光線落點(diǎn)坐標(biāo)。對(duì)紅外燈的大量出射光線追跡,可獲得紅外單燈在受熱平面上的光線落點(diǎn)坐標(biāo)數(shù)據(jù)庫(kù)。對(duì)N條光線進(jìn)行隨機(jī)抽樣并追跡獲得光線落點(diǎn)數(shù)據(jù)庫(kù)的程序流程圖如圖6所示。

圖6 光線追跡的流程圖

4 紅外單燈輻照模擬

采用利用蒙特卡洛方法對(duì)光線進(jìn)行隨機(jī)抽樣,抽樣光線的數(shù)量會(huì)對(duì)結(jié)果的準(zhǔn)確度產(chǎn)生很大影響[9],而誤差與抽樣光線總數(shù)的平方根成反比,故若將誤差控制在一定范圍內(nèi),需要追跡分析大量光線。當(dāng)追跡光線數(shù)量增加到100萬(wàn)條以上時(shí),可實(shí)現(xiàn)紅外燈輻照特性的高精度仿真。利用計(jì)算機(jī)對(duì)紅外燈出射光線進(jìn)行蒙特卡洛模擬、追跡的時(shí)間隨追跡光線的數(shù)量增多而成倍增長(zhǎng),使時(shí)效性變差。孟慶龍等為縮短燈陣優(yōu)化時(shí)間,對(duì)單個(gè)鏑燈光學(xué)模型追跡光線數(shù)量為104條[10]。為兼顧精確度與時(shí)效性,在實(shí)際的模擬計(jì)算中對(duì)單燈追跡光線數(shù)量為10萬(wàn)條。

選用帶有玻璃反射器的紅外燈,其額定電壓為230 V,額定功率為175 W,可將90%的能耗轉(zhuǎn)化為紅外熱。反射器最大處直徑為120 mm,反射器的等效旋轉(zhuǎn)拋物面高度為60 mm。假設(shè)燈絲為長(zhǎng)10 mm,直徑為5 mm的圓柱體。圓柱發(fā)光體的軸線中心與拋物面的焦點(diǎn)相重合,且兩者的旋轉(zhuǎn)對(duì)稱軸同為z軸。對(duì)紅外單燈最低處距受熱平面0.8 m時(shí)輻照情形進(jìn)行蒙特卡洛模擬,光線在受熱平面上的落點(diǎn)分布情況如圖7所示。

圖7 單燈光線落點(diǎn)散點(diǎn)圖

圖7中(a)、(b)、(c)、(d)分別為受熱平面上40 m×40 m、10 m×10 m、2 m×2 m以及0.6 m×0.6 m范圍的光線落點(diǎn)分布圖。由圖可知,在以紅外燈軸線與受熱平面交點(diǎn)位置(原點(diǎn)O)為中心的區(qū)域光線落點(diǎn)密度最大,由該點(diǎn)向外擴(kuò)展光線落點(diǎn)密度逐漸下降,也即紅外燈輻射能量以環(huán)狀形式由中心位置向四周逐漸遞減。

為對(duì)受熱平面受到的輻照分布進(jìn)行量化分析,需在受熱平面上劃分網(wǎng)格,對(duì)大量的落點(diǎn)坐標(biāo)進(jìn)行統(tǒng)計(jì)分析,網(wǎng)格劃分示意圖如圖8所示。圖中實(shí)線網(wǎng)格用于統(tǒng)計(jì)統(tǒng)計(jì)投射在相應(yīng)網(wǎng)格內(nèi)光線落點(diǎn)的數(shù)量。以圖中邊線加粗的網(wǎng)格為例,判斷某一光線落點(diǎn)是否在某網(wǎng)格內(nèi)的準(zhǔn)則如下:設(shè)落點(diǎn)坐標(biāo)是(x,y),若滿足:Ximin≤x

圖8 受熱平面網(wǎng)格劃分示意圖

按照上述網(wǎng)格劃分與相應(yīng)物理量的對(duì)應(yīng)關(guān)系,統(tǒng)計(jì)各網(wǎng)格內(nèi)的光線數(shù)量,仿真計(jì)算出紅外燈在受熱平面上0.6 m×0.6 m范圍內(nèi)的光線落點(diǎn)密度分布三維圖如圖9(a)所示，光線落點(diǎn)密度分布偽彩圖如圖9(b)所示。由圖9可以看出:紅外燈在目標(biāo)平面上靠近中心區(qū)域接收到的光線數(shù)量多,且呈環(huán)狀向四周逐漸遞減。也即對(duì)于隨機(jī)抽樣的光線而言,落在中心區(qū)域的概率更大。

圖9 光線分布密度分布圖

輻射照度W/m2分布揭示了紅外燈輻射能量在受熱平面上的分布情況。單條光線攜帶能量可由紅外燈輻射功率除以隨機(jī)抽樣光線總數(shù)得出,也即單位時(shí)間內(nèi)紅外燈向外輻射的能量由該時(shí)間內(nèi)出射的各抽樣光線均分。統(tǒng)計(jì)網(wǎng)格內(nèi)光線落點(diǎn)數(shù)量,根據(jù)每條光線攜帶能量即可計(jì)算出該網(wǎng)格內(nèi)接收的。網(wǎng)格的輻射能量除以網(wǎng)格面積,可得該受熱平面在該網(wǎng)格處的輻射照度。解算出受熱平面上各網(wǎng)格對(duì)應(yīng)的輻射照度,可得出紅外燈在受熱平面上的輻射照度分布情況。紅外燈在受熱面上0.6 m×0.6 m范圍內(nèi)的輻射照度分布情況如圖10,其中圖10(a)為照度分布三維圖,圖10(b)為輻射照度分布偽彩圖。

圖10 模擬輻射照度分布圖

為對(duì)建模的有效性進(jìn)行驗(yàn)證,依照?qǐng)D11對(duì)紅外單燈的輻射照度進(jìn)行測(cè)量。選用的紅外燈具有以燈軸線為對(duì)稱軸的旋轉(zhuǎn)對(duì)稱特點(diǎn),可以認(rèn)為其對(duì)外輻射也具備旋轉(zhuǎn)對(duì)稱的特點(diǎn)。以紅外燈軸線與受熱平面的交點(diǎn)為圓心O,在受熱平面上繪制一系列同心圓,半徑步長(zhǎng)為Δr。通過(guò)圓心繪制四條直線,相鄰兩直線夾角為60°。測(cè)量直線與各圓交點(diǎn)處的輻射照度,并將同圓上的測(cè)量值取平均值,得出距離圓心不同間距處的輻射照度。

圖11 單燈輻射照度測(cè)量原理圖

將測(cè)量結(jié)果進(jìn)行曲線擬合,得出輻射照度分布與距離軸線距離的函數(shù)關(guān)系,進(jìn)一步轉(zhuǎn)化為輻射照度與二維空間位置的關(guān)系。根據(jù)該函數(shù)關(guān)系可得受熱平面輻射照度分布圖如圖12,其中圖12(a) 、圖12(b)分布為照度分布三維圖與偽彩圖。

圖12 測(cè)量輻射照度分布圖

將模擬結(jié)果和實(shí)驗(yàn)測(cè)量結(jié)果進(jìn)行對(duì)比,兩者的輻射照度輪廓圖如圖13所示。其中,三角代表測(cè)得的輻射照度值,虛線為實(shí)測(cè)輻射照度與距燈軸線距離的擬合曲線,實(shí)線為利用計(jì)算機(jī)進(jìn)行蒙特卡洛模擬與追跡的輻射照度模擬曲線。

圖13 模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對(duì)比圖

由圖13可知,計(jì)算機(jī)模擬結(jié)果與實(shí)驗(yàn)測(cè)得結(jié)果具有一定差別,主要體現(xiàn)在燈軸線附近區(qū)域以及邊緣區(qū)域。實(shí)驗(yàn)結(jié)果顯示在紅外燈正下中心區(qū)域輻射照度最高且為尖峰,而模擬結(jié)果顯示該處呈現(xiàn)帶有凹陷的平臺(tái)。在邊緣位置,實(shí)測(cè)輻射照度值下降較為平緩且未降至0,模擬結(jié)果顯示輻射照度隨距燈軸線距離的增加而迅速下降直至為0。造成差異原因在于實(shí)際燈絲形狀為“C”形,且處于工作狀態(tài)時(shí)溫度分布不均勻,另外紅外燈前罩等結(jié)構(gòu)會(huì)對(duì)光線的傳播路徑造成影響,進(jìn)而影響到受熱平面上的輻射照度分布。

燈高0.8 m時(shí),在距紅外燈軸線0.2 m范圍內(nèi),測(cè)量輻射照度平均值為832.9 W/m2,最大值為1769.6 W/m2；模擬輻射照度平均值為824.9 W/m2,軸線與受熱平面交點(diǎn)處輻射照度值為1654.5 W/m2。模擬值的平均值誤差小于1%,最大值誤差為6.5%。由此可知,該模型可以較好地模擬紅外燈在受熱平面的輻射照度分布情況。

5 結(jié) 論

對(duì)紅外燈大量出射光線進(jìn)行了蒙特卡洛模擬與追跡,建立了紅外單燈在受熱平面上的光線落點(diǎn)坐標(biāo)數(shù)據(jù)庫(kù)。在受熱平面上劃分網(wǎng)格并統(tǒng)計(jì)各網(wǎng)格內(nèi)光線落點(diǎn)數(shù)量,得出紅外單燈在受熱平面上的輻射照度分布。實(shí)驗(yàn)測(cè)量了紅外單燈的輻射照度,并將其與計(jì)算機(jī)模擬結(jié)果進(jìn)行了對(duì)比,結(jié)果顯示利用蒙特卡洛方法建立的紅外燈輻射照度模型具有較高精度?？蓪⒃撃Ｐ徒Y(jié)合溫度解算模型,優(yōu)化紅外燈陣參數(shù)改善受熱平面上的溫度分布均勻程度,提高產(chǎn)品溫度環(huán)境試驗(yàn)的準(zhǔn)確性與可信度。

[1] CHEN Xun,TAO Junyong,ZHANG Chunhua.Reliability enhancement testing and accelerated life testing:an introductory review[J].Journal of National University of Defense Technology,2002,24(4):29-32.(in Chinese)

陳循,陶俊勇,張春華.可靠性強(qiáng)化試驗(yàn)與加速壽命試驗(yàn)綜述[J].國(guó)防科技大學(xué)學(xué)報(bào),2002,24(4):29-32.

[2] JIANG Pei,CHEN Xun,ZHANG Chunhua,et al.Reliability enhancement testing[J].Structure & Environment Engineering,2003,30(1):58-64.(in Chinese)

蔣培,陳循,張春華,等.可靠性強(qiáng)化試驗(yàn)技術(shù)綜述[J].強(qiáng)度與環(huán)境,2003,30(1):58-64.

[3] Ziemke R A.Infrar ed heater used in qualification testing of international space station radiators[R].NASA,TM-2004-212332,2004.

[4] LIU Shouwen,YIN Xiaofang,PEI Yifei,et al.The study of heat flux distribution for infrared lamp based on monte carlo method[J].Journal of Astronautics,2010,31(2):608-614.(in Chinese)

劉守文,尹曉芳,裴一飛,等.基于蒙特卡羅方法的紅外燈熱流分布研究[J].宇航學(xué)報(bào),2010,31(2):608-614.

[5] YANG Xiaoning,SUN Yuwei,YU Qianxu.The optimized design for improving flux uniformity of infrared lamp array[J].Spacecraft Environment Engineering,2012,29(1):27-31.(in Chinese)

楊曉寧,孫玉瑋,余謙虛.提高紅外燈陣熱流模擬均勻性的優(yōu)化設(shè)計(jì)方法[J].航天器環(huán)境工程,2012,29(1):27-31.

[6] YANG Guowei,SUN Xinming,PEI Yifei,et al.Study of simulation and optimization for heat flux distribution of infrared lamp array[J].Spacecraft Engineering,2011,20(1):134-141.(in Chinese)

楊國(guó)巍,蘇新明,裴一飛,等.紅外燈陣熱流分布仿真優(yōu)化研究[J].航天器工程,2011,20(1):134-141.

[7] ZHANG Kun,CHEN Haiqing,LIAO Zhaoshu,et al.Monte Carlo simulation for the radiation characteristics of infrared light source[J].Laser & Infrared,2010,40(5):491-495.(in Chinese)

張坤,陳海清,廖兆曙,等.紅外光源輻射特性的蒙特卡洛模擬[J].激光與紅外,2010,40(5):491-495.

[8] HUANG Deyun,REN Fei,LING Yipin,et al.A study on some aspects in designing of illuminating system[J].Mathematics in Practice & Theory,2003,33(7):102-107.(in Chinese)

黃德云,任飛,凌一品,等.光學(xué)照明系統(tǒng)設(shè)計(jì)中若干問題的研究[J].數(shù)學(xué)的實(shí)踐與認(rèn)識(shí),2003,33(7):102-107.

[9] ZHENG Xiaodong,WANG Yangchun,QIN Wenhong.Random error analysis of illumination distribution calculated by non-sequential ray tracing programs[J].Acta Photonica Sinica,2008,37(10):1970-1974.(in Chinese)

鄭曉東,汪揚(yáng)春,秦文紅.非序列光線追跡程序照度分布計(jì)算的隨機(jī)誤差分析[J].光子學(xué)報(bào),2008,37(10):1970-1974.

[10] MENG Qinglong,WANG Yuan.Irradiance characteristic and design optimization of a new full-spectrum solar simulator[J].Acta Energiae Solaris Sinica,2012,33(1):73-80.(in Chinese)

孟慶龍,王元.全光譜日光模擬器的輻照特性仿真研究與優(yōu)化[J].太陽(yáng)能學(xué)報(bào),2012,33(1):73-80.