诱饵球镜面反射内壁的辐射传递系数计算*

厉夫兵，冷俊敏，李红莲

(北京信息科技大学 信息与通信工程学院，北京 100101)

诱饵球镜面反射内壁的辐射传递系数计算*

厉夫兵，冷俊敏，李红莲

(北京信息科技大学 信息与通信工程学院，北京 100101)

针对内表面具备镜面反射特性的空间诱饵球，基于蒙特卡罗方法建立漫发射模型，推导光线的最大反射次数，利用光线在球腔的传播特性求解光线与球壁的多次反射交点，计算球面微元吸收的光束能量，获取不同反射率情况下单一球面微元对全部球面微元的辐射传递系数。计算结果表明，单一微元对球面各微元的辐射传递系数中存在2个峰值，分别为微元对其自身的辐射传递系数，以及微元对其球对称微元的辐射传递系数。随着反射率的增大，微元对其自身的辐射传递系数逐渐接近于微元对其球对称微元的辐射传递系数；当反射率接近于1时，两者的值近似相等。

辐射传递系数；辐射换热；镜面反射；蒙特卡罗；角系数；球

0 引言

在弹道目标飞行中段释放气球诱饵的方法可以降低目标的探测概率[1-2]，增加天基红外监视系统对目标的识别与跟踪难度[3-4]。公开文献表明，空间气球多采用铝膜-聚酯薄膜-铝膜的超薄层结构[5]。在日光照射下，气球向阳面与背阴面的温差会加剧气球内部的辐射换热，而辐射换热又反过来影响着目标的温度分布。因此，提高辐射换热的计算精度对提高目标温度的计算准确度具有重要意义。

近年来，国内外许多学者对影响弹头、诱饵等空间目标的温度与红外辐射特性的主要因素进行了研究[6-8]，获得了目标在日照区及地球阴影区的整体温度或温度分布随时间的变化[9-10]。传统方法计算目标内部辐射换热时，通常将目标内壁视为灰体面元[11-12]，利用辐射角系数建立传热方程，通过多次迭代求解灰体面元间的辐射换热，整个求解过程较费时间。张伟清等先利用辐射角系数求解灰体面元间的辐射传递系数[13]，再直接利用辐射传递系数求解灰体面元间的辐射换热，避免了直接利用角系数求解辐射换热时的迭代计算过程，加速了计算进程。值得注意的是，当气球内壁铝膜足够光滑时，气球内表面将具备良好的镜面反射特性而非漫反射特性，此时利用基于面元灰体假设的计算方法将会给辐射换热计算引入新的误差，因此需要计算镜反射面元间的辐射传递系数。

由于镜面反射的存在，面元的辐射热量不可避免地在球腔内多次反射，因此优先采用蒙特卡罗方法计算面元间的辐射传递系数。经典的蒙特卡罗方法利用面元的发射率随机吸射光线，通过统计面元吸收光线数目的方法计算辐射传递系数[14-15]，这种随机吸收光线的方法在面元发射率ε较大时会导致反射光线数量的急剧减少。例如，当系统发射率ε=0.9、面元发射光线数量N=10 000条时，经过系统的1，2，3次反射后，剩余光线数量分别约为1 000，100，10。多次反射发生后，由于剩余光线的样本数量太少，难以对某面元的反射光线在半空间的分布情况进行分析。通过大量增加发射光线数量的方法虽然能够缓解上述问题，但同时也引入了大量的碰撞与求交计算，限制了蒙特卡罗方法的计算速度。

本文通过建立镜反射内壁面元的蒙特卡罗漫发射模型，根据光线的最大反射次数计算每次反射事件发生时光线损失的能量份额，利用光线在球腔内的传播特性直接求解光线与内壁的多次反射交点，通过统计各次反射事件中面元吸收的光线能量份额计算面元间的辐射传递系数，在此基础上分析镜反射内壁的辐射传递系数分布特性。

1 蒙特卡罗漫发射模型

1.1 目标建模与微元剖分

气球目标被建模为半径为R0的球面，并在天顶角和方位角方向进行等间隔剖分，获得一系列球面微元，并利用微元在天顶角、方位角方向的位序确定微元的位置。由球的性质可知，内壁微元表面任意位置点的法矢量为该点指向球心的单位矢量。

1.2 光线的随机发射位置、方向与能量

在蒙特卡罗方法中，面元向半空间的辐射能量被视为从面元随机位置、以随机方向发射的光线的集合，每束光线都具备一定的能量[14-15]。假定面元向半空间发射的光线总数为N，则每束光线携带的初始能量为

(1)

式中：T为面元的瞬态温度；ρ为面元的反射率；σ为Steven-Boltzmann常数。

在球坐标系中，光线在球面微元内表面的随机发射位置可表示为

(2)

式中：θa，θb分别是球面微元天顶角的上、下界；φa，φb分别是微元方位角的上、下界；Rθ，Rφ是在[0，1]区间满足均匀分布的随机变量。

光线在随机发射位置点上的随机发射方向的求解过程可参考文献[14]。通过分析面元向半空间范围内投射光束的二维概率密度，可推导光束的二维分布函数并得到光束的发射方向与随机数的关系为

(3)

式中：θf，φf分别为发射光束的天顶角、方位角；Rθ，Rφ为[0，1]区间内均匀分布的随机数。

2 辐射传递系数的计算

2.1 光线在球腔中的传播

图1所示为参考坐标系OXYZ中球面微元A上光束的某次随机发射位置p(θp，φp)。为计算方便，建立与参考坐标系OXYZ重合的本体坐标系Oxyz，并将本体坐标系Oxyz先沿z轴逆时针旋转角度φp，再沿y轴顺时针旋转角度π-θp，使该坐标系的z轴负半轴恰好通过光线的发射位置p。因此，从本体坐标系Oxyz向参考坐标系OXYZ的坐标转换矩阵为

(4)

式中：Ly(β)，Lz(γ)分别为坐标系沿y轴、z轴逆时针旋转β，γ角度对应的旋转矩阵，即

图1 面元随机发射位置p在坐标系OXYZ与Oxyz中的位置Fig.1 Random position ‘p’ on a spherical facet A in OXYZ coordinates and Oxyz coordinates

如图2所示，在坐标系Oxyz中，光线从点 p[0，0，-R0] 以角度 (θf，φf) 发射。光线与球壁第1次相交于点d (1次反射点)，并遵循镜面反射定律继续向前传播。由于球内表面上任意点的法矢量都指向球心，故过点d的反射光线必定在由球心、点p和点d所确定的大圆切平面内。容易推得，发射光线的多次反射光线也全部落在该圆内。从图2中可以看出，发射光线与球壁的1次反射点d在大圆中的圆心角α1=π-2θf，光线与球壁的k次反射点对应的圆心角为

(5)

式中：%表示求余操作。

容易求出，k次反射点在坐标系Oxyz中的天顶角θk和方位角φk为

(6)

将坐标系Oxyz中的k次反射点坐标 (R0，θk，φk) 先转换为直角坐标，再利用式(4)的坐标转换矩阵L求反射点在坐标系OXYZ中的直角坐标，并将该直角坐标转换为球坐标(R0，θr，φr)，最后利用反射点在坐标系OXYZ中的天顶角θr和方位角φr求出反射点所在的球面微元。

图2 光线在球腔中的多次反射 (坐标系Oxyz)Fig.2 Multi-reflections of a light ray inside the sphere in Oxyz coordinates

2.2 辐射传递系数的计算

辐射传递系数是离开一个表面到达另一个表面的能量份额，其中包含了从任意面元反射后到达目标表面的能量。因此，球面微元的辐射传递系数需要考虑光线的多次反射贡献。

本文中，单束光线携带的初始能量e0是被面元依吸收率逐渐吸收而衰减的。当光线发生第1次反射时，反射光束的能量变为ρe0，被交点所属微元吸收的能量为(1-ρ)e0。当光线发生第2次反射时，反射光束的能量变为ρ2e0，被交点所属微元吸收的能量为(1-ρ)ρe0。以此类推，第k次反射时，反射光束的能量变为ρke0，被交点所属微元吸收的能量为

(7)

给定阈值η，当第k次反射发生时，若反射光束的能量ρke0与初始能量e0之比低于阈值η时视为光束能量被完全吸收，即ρk≤η，则最小反射次数nr的取值可表示为

(8)

式中：「⎤表示上取整。

综上所述，对于球面微元的任意一束初始发射光线，利用式(5)～(7)可一次求得其nr个反射点的所属微元及微元吸收能量。对微元的全部N束发射光线进行遍历，统计每束光线的全部nr个反射点的所属微元，并计算微元吸收能量，可求得该微元对其它微元的辐射传递系数为

(9)

式中：Ej为微元i的辐射能量中最终被微元j吸收的部分，包括微元i对微元j的直接入射能量和被其它面元反射后到达微元j的能量。

3 数值计算结果与分析

本例的参数设置如下。气球半径R0=1 m，内壁镜面铝膜的反射率为ρ，球面沿天顶角剖分18份 (dθ=10°)，沿方位角方向剖分36份 (dφ=10°)。对特定的球面微元，随机发射光线的数目N=100万条。光束被多次吸收、反射后，当反射光束的能量与初始能量之比低于阈值η=0.01时，辐射能量被视为全部吸收，光线的多次反射计算过程结束。

3.1 辐射角系数

当球内壁为黑体时 (ε=1，ρ=0)，辐射传递系数等于辐射角系数。球内壁任意两球面微元m，n之间的角系数理论值为[15]

(10)

式中：Sn为微元n的面积；θa，θb分别为微元n的天顶角上、下界；φa，φb分别为微元的方位角上、下界。

由式(10)可知，球面微元间的角系数仅与微元n的面积有关，与发射微元m的位置和面积无关。

图3a)为球面发射微元A(天顶角、方位角中心为 (135°，115°)) 对球内壁全部微元的角系数计算值二维展开图，垂直方向为球内壁微元的天顶角，水平方向为方位角。图3b)为3个特定的方位角下 (对应于图3a)方框) 的辐射角系数计算值与理论值对比，图3c)为计算值与理论值的相对误差。

图3 球面微元的角系数理论值与计算值Fig.3 Theoretical and calculated shape factors for sphere facets

从图3a)可以看出，该微元对球面各微元的辐射角系数在方位角方向保持相对稳定，仅在天顶角方向变化。从图3b)～c)可以看出，辐射角系数计算值与理论值很接近。在本例的参数设置下，相对误差大部分处于±5%的范围内，最大相对误差不超过10%。由理论式(10)可知，当球面沿方位角方向均匀剖分时，φb，φa为常数，角系数与方位角无关，由此证明了本文方法计算辐射角系数的正确性。

3.2 辐射传递系数

图4a)～c)分别为镜面反射率ρ=0.40，0.70，0.99时 (反射次数nr=6，13，459)，发射微元A(135°，115°) 对全体微元的辐射传递系数矩阵。从图中可以看出，与图3a)的黑体面元不同，镜面反射时，辐射传递系数沿方位角方向也发生变化，并呈现出几个特性。①随着反射率的增大，辐射传递系数存在两个较为明显的峰值区域，一个恰为微元A(135°，115°)自身，另一个是其球心对称微元B(45°，295°) (注：球面两点a，b关于球心对称时，其天顶角满足θa+θb=180°，方位角满足 |φa-φb|=180°)。②当反射率较低时，微元A的辐射传递系数小于微元B，随着反射率的增加，其辐射传递系数逐渐接近微元B。当反射率接近于1时，两者的值近似相等。

图4 发射微元A对球面微元的辐射传递系数Fig.4 Radiative heat transfer coefficient between a manua-lly selected facet A and the whole spherical facet

图5a)～b)为反射率ρ=0.99时 (nr=459)，2个不同位置的发射微元A(5°，55°)，(85°，85°) 对全体球面微元的辐射传递系数矩阵。从图5中可以看出，与图4c)类似，辐射传递系数矩阵中的两个峰值所在的微元在空间上关于球心对称。由此可以推知，当反射率较大时，任意位置的发射微元对全体球面微元的辐射传递系数均存在2个峰值，且分别位于微元自身及其球心对称微元。

接下来首先分析图5中辐射传递系数存在2个关于球心对称的峰值的原因。在图2所示的坐标系Oxyz中，假设微元A，B分别为发射微元及其球心对称微元。当反射率ρ增大时，单束光线的反射次数将增加；当ρ接近于1时，反射次数nr将非常大，各次反射点将散布于大圆圆周上。虽然单一光线的多次反射点在大圆圆周上的散布位置是确定的；但对于发射天顶角θf∈[0，π/2]的大量随机发射光线而言，其多次反射点落在大圆圆周各处的概率将近似相等。由于大圆切割球心对称微元A，B形成的弧长相等，因此，在方位角φ=φf上，反射点落在微元A上弧线的概率等于其落在微元B上弧线的概率。特别地，由图2可知，由于任意方位角φ∈[0，2π)的大圆都经过微元A，B，因此反射点落在微元A，B上的概率远高于其他微元。即反射率接近于1时，微元A，B的辐射传递系数近似相等且远高于其他球面微元。

图5 ρ=0.99时不同位置的发射微元A对应的辐射传递系数矩阵Fig.5 Radiative heat transfer coefficient matrix corresponding to two different manually selected facets A when reflectivity ρ is 0.99

对图4a)～c)中微元A的辐射传递系数随反射率ρ的增大而逐渐逼近微元B的原因进行分析时，需要对微元A，B吸收的光线能量成份进行研究。通过对微元A的每条发射光线进行追踪，记录微元上的各个反射点都源于光线的第几次反射行为，然后统计微元上各次反射行为发生的次数。图6所示为在图4参数设置下，反射率ρ=0.40，0.70，0.83时，微元A，B上光线的第k次反射行为 (k=1，2，…，nr) 发生的次数统计。

图6 不同反射率时面元A，B上各次反射发生的次数Fig.6 Statistics of multi-reflection of rays on facet A and its spherically symmetric facet B

当η=0.01，ρ=0.4时，由式(8)可知光线发生6次反射即被完全吸收。如图6a)所示，当k=1时微元A，B上光线的第1次反射行为发生的次数近似相等，即微元A发出的光线中，与微元A，B直接发生碰撞的光线数量近似相等。由于微元A，B是关于球心对称的微元，面积相等，因此光线与两面元发生碰撞的概率相等，这与节3.1中辐射角系数的计算结论是一致的。

k=2表示光线的第2次反射行为，即发射微元A发出的光线经过球面微元的1次镜面反射后再与微元A或B相交。由图2可知，光线由微元A发出后，只有那些能够与微元B(或其附近微元) 发生直接碰撞的光线才有可能经1次反射后与微元A相交，此部分光线数量较少。但是，光线由微元A发出后，所有能够与“赤道”微元 (与Oxy平面相交的球面微元) 发生直接碰撞的光线经过1次镜面反射后，都有可能与微元B相交。因此，2次反射点落在微元B上的光线数量要远大于落在微元A上的光线数量，即微元B吸收的光线的2次碰撞能量远大于微元A，直接导致此时微元B的辐射传递系数较微元A大。

由图6a)～c)可以看出，当反射率ρ增大时，反射次数增多，且光线的第k次反射行为 (k≥3) 发生在微元A，B上的数目趋于相等。此时，虽然微元B吸收的光线的2次碰撞能量仍较微元A高，但相对于总能量而言，其所占的能量份额随反射率的增大而逐渐减小。因此，如图4所示，随着反射率的增大，微元A的辐射传递系数逐渐接近于微元B。

在本例中，光线被多次反射后，其携带的能量会逐渐降低，若反射光射的能量与初始能量之比低于阈值η=0.01，则光线能量被视为全部吸收，即光束被吸收的能量份额为0.99。当该阈值减小时，计算精度将会提高。图7所示为η=0.000 1、其他设置与图5相同时，得到的辐射传递系数矩阵。

图7 η=0.000 1，ρ=0.99时不同位置的微元A对应的辐射传递系数矩阵 Fig.7 Radiative heat transfer coefficient matrix corresponding to two different manually selected facets when η=0.000 1，ρ is 0.99

从图7与图5可以看出计算结果非常相近。这是因为当阈值η=0.000 1时，光束被吸收的能量份额为0.999 9；而当η=0.01时，光束被吸收的能量份额为0.99。由此可见，两者的能量变化不足1%，因此图7得到的结果与图5类似，但数据精度更高。

4 结束语

本文利用蒙特卡罗方法计算了球内壁的辐射角系数及不同镜面反射率情况下的辐射传递系数。辐射角系数的计算值与理论值一致。辐射传递系数的结果表明，首先，微元对全部球面微元的辐射传递系数矩阵中存在2个峰值，即微元对其自身、以及微元对其球心对称微元的辐射传递系数。其次，随着反射率的增大，微元对其自身的辐射传递系数逐渐接近于微元对其球心对称微元的辐射传递系数。当反射率接近于1时，两者的辐射传递系数值近似相等。由此可知，对于诱饵球镜反射内壁微元而言，当某微元的辐射出射能量一定时，反射率越大，辐射能量越倾向于向该微元及其球对称微元集中。本文的结论有助于分析空间诱饵球向阳面高温微元与背阴面低温微元之间的辐射换热与温度变化过程。

[1] 郦苏丹，任萱.大气层外诱饵释放研究[J].宇航学报，2001，22(2)：100-104. LI Su-dan，REN Xuan.Research of Releasing Decoy Outside Atmosphere[J].Journal of Astronautics，2001，22(2)：100-104.

[2] 林两魁，谢恺，徐晖，等.中段弹道目标群的红外成像仿真研究[J].红外与毫米波学报，2009，28(3)：218-223. LIN Liang-kui，XIE Kai，XU Hui，et al.Research on Infrared Imaging Simulation of Midcourse Ballistic Object Target Complex[J].Journal of Infrared and Millimeter Waves，2009，28(3)：218-223.

[3] 王硕，张奕群，孙冰岩.红外点目标跟踪方法综述[J].现代防御技术，2016，44(2)：124-134. WANG Shuo，ZHANG Yi-qun，SUN Bing-yan.Review of Point Target Tracking Methods Based on Infrared Sensor[J].Modern Defence Technology，2016，44(2)：124-134.

[4] 浦甲伦，崔乃刚，郭继峰.天基红外预警卫星系统及其探测能力分析[J].现代防御技术，2008，36(4)：68-72. PU Jia-lun，CUI Nai-gang，GUO Ji-feng.Space-Based Infrared System and the Analysis of Its Detecting Capability[J].Modern Defence Technology，2008，36(4)：68-72.

[5] SESSLER A M，CORNWALL J M，DIETZ B，et al.Countermeasures：a Technical Evaluation of the Operational Effectiveness of the Planned US National Missile Defense System[R].Cambridge，MA：Union of Concerned Scientists and MIT Security Studies Program，April 2000.

[6] 姚连兴，侯秋萍，罗继强.弹道导弹中段目标表面温度与红外突防研究[J].航天电子对抗，2005，21(2)：5-16. YAO Lian-xing，HOU Qiu-ping，LUO Ji-qiang.Research on Surface Temperature and Infrared Signature of the Ballistic Warhead in Midcourse[J].Aerospace Electronic Warefare，2005，21(2)：5-16.

[7] 汪洪源，陈赟.天基空间目标红外动态辐射特性建模与仿真[J].红外与激光工程，2016，45(5)：0504002. WANG Hong-yuan，CHEN Yun.Modeling and Simulation of Infrared Dynamic Characteristics of Space-Based Targets[J].Infrared and Laser Engineering，2016，45(5)：0504002.

[8] 敬韩博，王英瑞.临近空间高超声速目标天基红外探测技术研究[J].现代防御技术，2016，44(6)：80-84. JING Han-bo，WANG Ying-rui.Space-Based Infrared Detection for Near Space Hyptrsonic Targets[J].Modern Defence Technology，2016，44(6)：80-84.

[9] ZHU D Q，SHEN W T，CAI G B，et al.Numerical Simulation and Experimental Study of Factors Influencing the Optical Characteristics of a Spatial Target[J].Applied Thermal Engineering，2013(50)：749-762.

[10] 张骏，杨华，凌永顺，等.弹道导弹中段弹头表面温度场分布理论分析[J].红外与激光工程，2005，34(5)：583-586. ZHANG Jun，YANG Hua，LING Yong-shun，et al.Theoretical Analysis of Temperature Field on the Surface of Ballistic Missile Warhead in Midcourse[J].Infrared and Laser Engineering，2005，34(5)：583-586.

[11] LI F B，XU X J.Modeling Time-Evolving Infrared Characteristics for Space Object with Mircromotions[J].IEEE Transactions on Aerospace and Electronic Systems，2012，48(4)：3567-3577.

[12] 厉夫兵，许小剑.计算超薄多层介质诱饵温度的快速算法[J].红外与激光工程，2011，40(6)：986-991.

LI Fu-bing，XU Xiao-jian.Fast Algorithm for Temperature Calculation of Multi-Layered Ultra-Thin Decoy[J].Infrared and Laser Engineering，2011，40(6)：986-991.

[13] 张伟清，宣益民，韩玉阁.单元表面间辐射传递系数的新型计算方法[J].宇航学报，2005，26(1)：77-81. ZHANG Wei-qing，XUAN Yi-min，HAN Yu-ge.A New Method of Radiative Transfer Coefficient between Unit Surfaces[J].Journal of Astronautics，2005，26(1)：77-81.

[14] HOWELL J R，PERLMUTTER M.Monte Carlo Solution of Thermal Transfer Through Radiant Media between Gray Walls[J].Journal of Heat Transfer，1964，86(1)：116-122.

[15] 卞伯绘.辐射换热的分析与计算[M].北京：清华大学出版社，1988：81-83，144-145. BIAN Be-hui.Analysis and Calculation of Radiative Heat Transfer[M].Beijing:Tsinghua Universtiy Press,1988:81-83,144-145.

Radiative Heat Transfer Coefficient Calculation of Inflatable Sphere with Specular Inner Surface

LI Fu-bing，LENG Jun-min，LI Hong-lian

(Beijing Information Science & Technology University，School of Information & Communication Engineering，Beijing 100101，China)

Monte Carlo technique is used to calculate the radiative heat transfer coefficient for an inflatable sphere with specular reflection inner surface. Light rays are emitted diffusely from a spherical facet and being traced. The maximum number of times a ray can be reflected is derived, and the multi-intersections of a ray with the inner side of spherical surface are obtained. Radiative heat transfer coefficients between a manually selected facet and the whole spherical facet are calculated at different reflectivities. The results indicate that there are two peaks in the calculated coefficients, the smaller one is the coefficient between the selected facet and itself (i.e., the self-to-self coefficient), and the larger one is between the selected facet and its spherically symmetric facet (i.e., the self-to-spherical-symmetry coefficient). The self-to-self radiative heat transfer coefficient gradually approaches and finally equals the self-to-spherical-symmetry coefficient as the reflectivity increases to 1.

radiation transfer coefficient；radiation heat transfer；specular reflection；Monte Carlo；shape factor；sphere

2016-12-01；

2017-02-11

北京市自然科学基金资助项目(3164043)

厉夫兵(1982-)，男，山东日照人。讲师，博士，主要从事目标红外辐射特性建模研究。

10.3969/j.issn.1009-086x.2017.03.030

TK124；TP391.9

A

1009-086X(2017)-03-0193-07

通信地址：100101 北京市朝阳区北四环中路35号北京信息科技大学信息与通信工程学院

E-mail：lifubing@bistu.edu.cn