基于蒙特卡洛方法的红外光源辐照特性模拟

宋伟浩,周贤文,张 薇

(上海海事大学航运技术与控制工程交通行业重点实验室,上海 201306)

1 引 言

产品在运输、贮存和使用过程中,温度对产品的使用寿命及可靠性有着重要影响[1-2]。为确保产品的可靠性和使用的安全性,有必要进行温度可靠性环境试验。以红外灯阵列作为辐射热源进行温度测试,通过优化灯阵可有效提高受热平面温度分布均匀程度,进而提高温度环境试验的可信度与准确性。灯阵优化工作的基础是建立红外灯的辐射照度分布模型。

Ziemke等采用有限元法计算了红外灯阵的辐射照度分布[3]。刘守文等利用蒙特卡洛方法对灯管半侧镀有反射膜的管状红外灯进行了建模,实验结果显示所建立模型具有较高精度[4]。杨晓宁等采用蒙特卡洛方法模拟了无反射涂层的管状红外灯辐照分布情况,通过优化灯阵改善了受热面上热流分布均匀性[5]。杨国巍等学者用蒙特卡洛模拟了管状红外灯热流分布,建立了热流分布数据库,最后采用遗传算法优化了红外灯阵[6]。

张坤、陈海清等针对带有抛物面型反射器的红外灯泡辐射特性进行了模拟仿真研究[7],未考虑对光线方向抽样半球空间的方位变化。

许多学者多以无反射涂层或在灯管半圆周涂有反射层的管状红外灯为研究对象,采用蒙特卡洛方法模拟红外单灯的热流分布,对于带有抛物面反射器红外灯的建模与研究工作较少。

2 点光源的蒙特卡洛模拟

利用蒙特卡洛方法对红外灯辐照特性进行模拟,要按照一定的概率模型对出射光线进行重复的随机抽样生成光线的相关参数。利用计算机生成服从特定概率分布的关于光线参数的随机数,是对红外灯辐照特性进行蒙特卡洛模拟的基础。

在点光源的理想模型中,其出射光线均匀散布在以该点光源为球心的球体空间内,对光线的方向向量随机抽样,可唯一确定某条随机光线[8]。设辐射线的天顶角为θ;圆周角为φ;与球面交点为P;则向量OP即为抽样光线的方向向量,示意图如图1所示。对点光源的出射光线进行随机抽样,可等同于对单位球面上随机点P抽样。

图1 点光源光线随机抽样示意图

设P点位置服从球面上的均匀分布,且圆周角φ和天顶角θ相互独立,取服从(0,1)均匀分布的随机数rφ与rθ,可表示出点光源出射光线方向向量nx,ny,nz,如式(1):

(1)

3 扩展光源的蒙特卡洛模拟

扩展光源可以视作点光源的集合体,但发光特点又与点光源有所区别,需对其光线出射位置与方向进行抽样。为简化分析,将钨丝红外灯灯丝简化为圆柱体扩展光源,仅由圆柱侧面向外出射光线。假设:红外灯圆柱发光体为漫灰体,各处温度均匀一致；每条出射红外光携带的能量相等,与出射位置、角度以及波长无关。

对圆柱面上发光点的位置抽样,假定圆柱发光体半径为r,高度为h,轴线与z轴重合,下表面平面方程为z=h0。对光线出射点随机抽样的示意图如图2所示。随机发光点M位置可由圆周角α和竖坐标z表示。

图2 圆柱面发光点坐标抽样示意图

M点服从圆柱面上均匀分布,其落在阴影部分面积的联合概率分布函数,如式(2)所示:

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

(2)

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

由式(2)可知,α和z相互独立,故可取服从(0,1)均匀分布的随机数r1和r2分离变量,α和z可以由式(3)表示:

(3)

进而点M(x0,y0,z0)坐标可由式(4)描述。

(4)

圆柱体扩展光源的出射光线散布于以圆柱在M点外向法线为旋转对称轴的半球空间内,且出射方向遵循兰贝特定律。半球空间随抽样点M位置变动而变动。为确定M点出射光线的方向向量,需以M点为原点O′,建立辅助三维空间直角坐标系O′-x′y′z′。辅助坐标系O′-x′y′z′的建立示意图如图3。

图3 辅助坐标系示意图

设在坐标系O-xyz下,圆柱的方程为F(x,y,z)=0,x′、y′、z′轴方向向量分别为a、b、c。以圆柱体在点M处的外法向量为z′轴方向向量c,可由式(5)表示:

c=k1,k2,k3

(5)

a=b×c

(6)

将a,b,c单位化,分别得到O-xyz坐标系下的x′轴、y′轴和z′轴的单位方向向量:l1,l2,l3,m1,m2,m3和n1,n2,n3。

在O′-x′y′z′坐标系内,对M点也即O′点出射光线进行蒙特卡洛抽样,示意图如图4所示。

图4 圆柱面出射光线随机抽样示意图

光线O′Q方向可由圆周角β与天顶角γ表示,且可按照兰贝特定律处理。取服从(0,1)均匀分布的随机数r3和r4,则β与γ可由式(7)表示:

(7)

进而可得在O′-x′y′z′坐标系内光线方向向量(n′x,n′y,n′z),如式(8)所示:

(8)

为进行后续的光线追迹,需将该方向向量按照式(9)转化为坐标系O-xyz下单位光线方向向量(i0,j0,k0):

(9)

已知随机点M坐标(x0,y0,z0)和光线在O-xyz坐标系下的单位方向向量(i0,j0,k0),则光线在O-xyz坐标系内的直线参数方程可由式(10)表示:

(10)

至此,运用蒙特卡洛方法完成了对圆柱发光体一条光线的随机抽样。为研究受热平面上的辐射照度分布,需对大量光线进行随机抽样,并对每条光线进行追迹,获取其在受热平面上的落点坐标。在受热平面上划分网格,统计各网格内光线落点数量,进一步根据单光线携带能量值和网格面积,计算得出红外灯在受热平面的辐射照度分布[7]。

对光线进行追迹,需要根据所选红外灯的光学结构特点进行解算。选用带有玻璃反射器的红外灯泡,其示意图如图5所示,其中图5(a)是红外灯的实物图,图5(b)是红外灯等效的数学模型图。

图5 红外灯实物图与数学模型简图

忽略玻璃灯罩对光路的影响与发光圆柱体对光线的遮挡作用。则圆柱发光体的出射光线会直接或者经过反射器反射后抵达受热平面。参考实际反射器形状可将其等效为一旋转抛物面,圆柱发光体的轴线与抛物面旋转对称轴重合,且圆柱几何中心与和抛物面焦点相重合。开口向下且旋转对称轴为z轴的旋转抛物面可由式(11)描述:

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

(11)

z≥h1

式中,F为旋转抛物面的焦距,m；h1为反射器下边缘高度,m；h2为红外灯的反射器的高度,m。

联立光线所在直线方程和抛物面方程,可解出直线与抛物面的一个或者两个交点。当存在一个交点时,该交点为光线所在射线的反向延长线与抛物面的交点,属无效交点,此时光线将直接投射到受热平面。当存在两个交点时,其中一个为无效交点,另一个为有效交点。由发光点指向有效交点的向量与光线的方向向量夹角为零度,据此可筛选出入射点位置坐标。

光线在反光罩上发生发射时,对光线进行追迹需要求解其反射光线。联立式(10)和式(11)求解并经过判别得入射点Q0坐标。根据反射定律可解得反射光线的方向向量。

(12)

随机光线可能在反射器内经过多次反射才能抵达受热平面。当光线经过多次反射时,可将求解的反射光线作为下次反射的入射光线。当判定光线可直达受热平面时,联立光线所在直线方程和接收平面方程,可解得光线落点坐标。对红外灯的大量出射光线追迹,可获得红外单灯在受热平面上的光线落点坐标数据库。对N条光线进行随机抽样并追迹获得光线落点数据库的程序流程图如图6所示。

图6 光线追迹的流程图

4 红外单灯辐照模拟

采用利用蒙特卡洛方法对光线进行随机抽样,抽样光线的数量会对结果的准确度产生很大影响[9],而误差与抽样光线总数的平方根成反比,故若将误差控制在一定范围内,需要追迹分析大量光线。当追迹光线数量增加到100万条以上时,可实现红外灯辐照特性的高精度仿真。利用计算机对红外灯出射光线进行蒙特卡洛模拟、追迹的时间随追迹光线的数量增多而成倍增长,使时效性变差。孟庆龙等为缩短灯阵优化时间,对单个镝灯光学模型追迹光线数量为104条[10]。为兼顾精确度与时效性,在实际的模拟计算中对单灯追迹光线数量为10万条。

选用带有玻璃反射器的红外灯,其额定电压为230 V,额定功率为175 W,可将90%的能耗转化为红外热。反射器最大处直径为120 mm,反射器的等效旋转抛物面高度为60 mm。假设灯丝为长10 mm,直径为5 mm的圆柱体。圆柱发光体的轴线中心与抛物面的焦点相重合,且两者的旋转对称轴同为z轴。对红外单灯最低处距受热平面0.8 m时辐照情形进行蒙特卡洛模拟,光线在受热平面上的落点分布情况如图7所示。

图7 单灯光线落点散点图

图7中(a)、(b)、(c)、(d)分别为受热平面上40 m×40 m、10 m×10 m、2 m×2 m以及0.6 m×0.6 m范围的光线落点分布图。由图可知,在以红外灯轴线与受热平面交点位置(原点O)为中心的区域光线落点密度最大,由该点向外扩展光线落点密度逐渐下降,也即红外灯辐射能量以环状形式由中心位置向四周逐渐递减。

为对受热平面受到的辐照分布进行量化分析,需在受热平面上划分网格,对大量的落点坐标进行统计分析,网格划分示意图如图8所示。图中实线网格用于统计统计投射在相应网格内光线落点的数量。以图中边线加粗的网格为例,判断某一光线落点是否在某网格内的准则如下:设落点坐标是(x,y),若满足:Ximin≤x

图8 受热平面网格划分示意图

按照上述网格划分与相应物理量的对应关系,统计各网格内的光线数量,仿真计算出红外灯在受热平面上0.6 m×0.6 m范围内的光线落点密度分布三维图如图9(a)所示，光线落点密度分布伪彩图如图9(b)所示。由图9可以看出:红外灯在目标平面上靠近中心区域接收到的光线数量多,且呈环状向四周逐渐递减。也即对于随机抽样的光线而言,落在中心区域的概率更大。

图9 光线分布密度分布图

辐射照度W/m2分布揭示了红外灯辐射能量在受热平面上的分布情况。单条光线携带能量可由红外灯辐射功率除以随机抽样光线总数得出,也即单位时间内红外灯向外辐射的能量由该时间内出射的各抽样光线均分。统计网格内光线落点数量,根据每条光线携带能量即可计算出该网格内接收的。网格的辐射能量除以网格面积,可得该受热平面在该网格处的辐射照度。解算出受热平面上各网格对应的辐射照度,可得出红外灯在受热平面上的辐射照度分布情况。红外灯在受热面上0.6 m×0.6 m范围内的辐射照度分布情况如图10,其中图10(a)为照度分布三维图,图10(b)为辐射照度分布伪彩图。

图10 模拟辐射照度分布图

为对建模的有效性进行验证,依照图11对红外单灯的辐射照度进行测量。选用的红外灯具有以灯轴线为对称轴的旋转对称特点,可以认为其对外辐射也具备旋转对称的特点。以红外灯轴线与受热平面的交点为圆心O,在受热平面上绘制一系列同心圆,半径步长为Δr。通过圆心绘制四条直线,相邻两直线夹角为60°。测量直线与各圆交点处的辐射照度,并将同圆上的测量值取平均值,得出距离圆心不同间距处的辐射照度。

图11 单灯辐射照度测量原理图

将测量结果进行曲线拟合,得出辐射照度分布与距离轴线距离的函数关系,进一步转化为辐射照度与二维空间位置的关系。根据该函数关系可得受热平面辐射照度分布图如图12,其中图12(a) 、图12(b)分布为照度分布三维图与伪彩图。

图12 测量辐射照度分布图

将模拟结果和实验测量结果进行对比,两者的辐射照度轮廓图如图13所示。其中,三角代表测得的辐射照度值,虚线为实测辐射照度与距灯轴线距离的拟合曲线,实线为利用计算机进行蒙特卡洛模拟与追迹的辐射照度模拟曲线。

图13 模拟结果与实验数据对比图

由图13可知,计算机模拟结果与实验测得结果具有一定差别,主要体现在灯轴线附近区域以及边缘区域。实验结果显示在红外灯正下中心区域辐射照度最高且为尖峰,而模拟结果显示该处呈现带有凹陷的平台。在边缘位置,实测辐射照度值下降较为平缓且未降至0,模拟结果显示辐射照度随距灯轴线距离的增加而迅速下降直至为0。造成差异原因在于实际灯丝形状为“C”形,且处于工作状态时温度分布不均匀,另外红外灯前罩等结构会对光线的传播路径造成影响,进而影响到受热平面上的辐射照度分布。

灯高0.8 m时,在距红外灯轴线0.2 m范围内,测量辐射照度平均值为832.9 W/m2,最大值为1769.6 W/m2；模拟辐射照度平均值为824.9 W/m2,轴线与受热平面交点处辐射照度值为1654.5 W/m2。模拟值的平均值误差小于1%,最大值误差为6.5%。由此可知,该模型可以较好地模拟红外灯在受热平面的辐射照度分布情况。

5 结 论

对红外灯大量出射光线进行了蒙特卡洛模拟与追迹,建立了红外单灯在受热平面上的光线落点坐标数据库。在受热平面上划分网格并统计各网格内光线落点数量,得出红外单灯在受热平面上的辐射照度分布。实验测量了红外单灯的辐射照度,并将其与计算机模拟结果进行了对比,结果显示利用蒙特卡洛方法建立的红外灯辐射照度模型具有较高精度。可将该模型结合温度解算模型,优化红外灯阵参数改善受热平面上的温度分布均匀程度,提高产品温度环境试验的准确性与可信度。

[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)

陈循,陶俊勇,张春华.可靠性强化试验与加速寿命试验综述[J].国防科技大学学报,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)

蒋培,陈循,张春华,等.可靠性强化试验技术综述[J].强度与环境,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].宇航学报,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)

杨晓宁,孙玉玮,余谦虚.提高红外灯阵热流模拟均匀性的优化设计方法[J].航天器环境工程,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)

杨国巍,苏新明,裴一飞,等.红外灯阵热流分布仿真优化研究[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)

黄德云,任飞,凌一品,等.光学照明系统设计中若干问题的研究[J].数学的实践与认识,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)

郑晓东,汪扬春,秦文红.非序列光线追迹程序照度分布计算的随机误差分析[J].光子学报,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)

孟庆龙,王元.全光谱日光模拟器的辐照特性仿真研究与优化[J].太阳能学报,2012,33(1):73-80.