蒙特卡洛法基本原理

∗
(2.10)
Φ i , j =σ ⋅ RDi*, j (Ti 4 − T j4 )
根据定义， RDi , j 可表示为：
(2.11)
RDi , j =Φ i , j / Φ i =N i , j / N i
(2.12)
式中， N i 是代表表面单元 Ai 发射辐射能的抽样能束数， N i , j 是统计获得的最终被表
= Φ i , j σ ( Aiε iTi 4 ⋅ RDi , j − A j ε jT j4 ⋅ RD j , i )
可以证明，表面间的辐射传递系数 RDi , j 与 RD j , i 存在相对关系：
(2.9)
Aiε i ⋅ RDi , j =A j ε j ⋅ RD j , i
引入归一化辐射传递系数 RDi= RDi , j ⋅ ε i ⋅ Ai ，式(2.10)可改写为： ,j
r= 0
2 2 2 rmin + Rr ⋅ (rmax − rmin )
ϕ 0 = ϕ min + Rϕ ⋅ (ϕ max − ϕ min )
(2.21)
式中， Rr 、 Rϕ 分别是发射点的径向和圆周方向分布随机数。相应的发射点直角坐标 为 (r0 cos ϕ 0 , r0 sin ϕ 0 , 0) 。 对圆柱表面、圆锥(台)表面、球形(球冠、球带) 表面等典型旋转曲面，在圆柱坐标系 下建立发射点的分布概率模型比较方便。 以球形表面为例， 如图 2.5 所示。 若球面半径为 r0 ， 极角与圆周角的值域分别为 [θ min , 概率模型：

n= ±( Fx i + Fy j + Fz k ) /
Fx2 + Fy2 + Fz2
(2.15)
式中， Fx 、 Fy 、 Fz 是函数 F ( x, y, z ) 的偏导数； i 、 j 、 k 分别是 x 、 y 、 z 三个
坐标轴的方向矢量。 另外，由于某一部件表面只是其方程所表示表面中的部分区域，其边界约束可以通 过将相关的约束表面方程改写成解析不等式来表示。即
y = f1 ( x) 、 y = f 2 ( x) 分 别 是 y 方 向 上 游 、 下 游 两 边 界 的 直 线 方 程 ， 对 矩 形 平 面 ， f1 ( x) = ymin 、 f 2 ( x) = ymax 。
对三角形、梯形等其它形状平面，可先按矩形平面进行发射点随机抽样，然后根据边 界约束条件剔除边界外的点。 这样处理比较简便， 缺点是需进行一系列的边界约束条件判断， 增加了计算量。 对大规模的计算， 由此引起的计算量增加是比较明显。 为了直接在指定形状、 大小的平面内进行发射点抽样，可结合具体的形状、尺寸构造发射点概率模型。 另外一类典型的发射平面是圆形、扇形、环型平面，可统一为环扇形平面。对这类平 面， 先按极坐标形式的发射点概率模型进行抽样， 然后再将发射点坐标转换为直角坐标[27-30]。 极坐标形式的发射点概率模型如下：
x0 = r0 ⋅ sin θ 0 ⋅ cos ϕ 0 y0 = r0 ⋅ sin θ 0 ⋅ sin ϕ 0 z0 = r0 ⋅ (cos θ 0 − cos θ max )
除圆柱面、球形表面外，还有一些表面是旋转抛物面、旋转椭球面或旋转双曲面等复 杂的旋转曲面。 对这些复杂曲面， 原则上也可采用与前面类似的方法建立发射点分布概率模 型。如某曲面在柱坐标系中的方程为 z = f ( r ) ，则表面上发射点的径向分布概率模型原则 上可表示为：
，25]
。
在辐射传递计算中应用的 MCM 主要有两种方法： (1) 抽样能束携带能量， 概率模拟和 能量平衡方程的求解没有分离；(2) 抽样能束不携带能量，概率模拟和温度场的迭代计算分 离。前一种方法中，温度场每变化一次就需要重新进行一次 MCM 模拟计算，计算量巨大； 后一种方法中，由 MCM 模拟计算的是各表面之间的辐射换热能量份额，通常以辐射传递系 数（或辐射传递因子、辐射交换因子、辐射网络系数）表示，只要各表面的辐射物性不变， 辐射传递系数不变。当温度场变化时，只需重新求解能量平衡方程即可，因此，目前普遍采 用后一种方法。 改进的 MCM 鉴于分解难点的思想，抽样能束本身不携带能量，利用概率模拟求辐射 传递系数 RDi , j ，然后将 RDi , j 代入能量方程。 表面间的辐射传递系数 RDi , j 的定义为：在一个换热系统中，由表面单元 Ai 发射的辐 射能经直接投射以及系统中任意表面的一次或多次散射（反射、折射和衍射）后，最终被表 面单元 A j 吸收的份额。显然，若已知 RDi , j ，两表面间的辐射换热量为：
r0 ∂f Rr = ⋅ rdr ∫rmin 1 + ∂r 2
(2.23)
∫r
rmax
min
∂f 1 + ⋅ rdr ∂r
2
(2.24)
ϕmax
x
θmax
ϕmin
r0
• ( x0 , y0 , z0 )
θmin
z
图 2.5 球形表面发射点坐标系 其中， Rr 是发射点的径向分布随机数， rmax 和 rmin 分别是曲面的最大和最小半径。在
2.8.1 能束光线的发射点分布概率模型
温度、发射率均匀的表面或表面单元的发射点分布均匀。选择合适的发射点坐标系、 并根据表面的具体形状构造发射点分布概率模型， 对保证模拟计算中发射点抽样的均匀、 随 机分布很重要。 对平面，以平面内某点为发射点坐标系原点、其表面法向为某一坐标轴比较合适，如 为 z 轴，则 x 、 y 两坐标轴在平面内。最简单的是矩形与平行四边形平面，可使其相对的两 条边与 x 轴平行，则发射点 ( x0 , y0 , z0 ) 的分布概率模型为：
x0 = xmin + Rx ⋅ ( xmax − xmin ) z0 = 0
y = f1 ( x0 ) + Ry ⋅ [ f 2 ( x0 ) − f1 ( x0 ) ] 0
(2.20)
式中， [ xmin , xmax ] 是 x 值域， Rx 、 Ry 分别是 x 、 y 方向的发射点分布随机数，
m = m1 i + m2 j + m3k 。
= m1 sin θ ⋅ cos ϕ = m2 sin θ ⋅ sin ϕ m3 = cos θ
(2.19)
2.8 表面辐射传递的主要概率模型
表面辐射传递的概率模型是描述表面发射、 反射、 吸收统计规律的数学模型， 根据热辐射的基本定律可构造辐射传递各子过程的概率模型。
Fk ( x, y, z ) ≥ 0
(2.16)
其中， Fk ( x, y, z ) = 0 是第 k 个约束表面方程。这样，由表面方程、正法向和边界约 束三部分构成了对一个实际部件表面的完整数学描述。 用几何射线来模拟辐射能束，以参数方程形式表示。设辐射能束的源点坐标为
( x0 , y0 , z0 ) 、方向矢量为 m = m1 i + m2 j + m3k ，则辐射能束的参数方程为：
对结构复杂的系统，特别是各物体表面之间存在复杂的遮挡、部分遮挡关系，应用上 述求解方法困难很大。 目前， 蒙特卡罗法已成为复杂几何系统内辐射换热计算的最主要方法。
2.6 蒙特卡罗法计算辐射换热的基本原理
蒙特卡罗方法（Monte－Carlo Method 或MCM）是一种概率模拟方法，它是通过随机 变量的统计试验来求解数学物理或工程技术问题的一种数值方法。早期的随机试验是用投 针、 掷骰子、 掷钱币等方法进行， 由于受模拟试验工具的限制， 能真正解决的实际问题很少。 20 世纪 40 年代中期， 由于电子计算机的发明，MCM首先在核武器的研制中得到应用。 1964 年，Howell将MCM引入到辐射换热计算领域[18-24]。 MCM模拟计算的基本思想是：将热辐射的传输过程分解为发射、反射、吸收、散射等 一系列独立的子过程， 并建立每个子过程的概率模型。 令每个单元发射一定量的能束， 跟踪、 统计每个能束的归宿（被哪些单元吸收，或反射） ，从而得到该单元辐射能量分配的统计结 果[24
F ( x, y, z ) = C1 x 2 + C2 y 2 + C3 z 2 + C4 xy + C5 xz + C6 yz + C7 x + C8 y + C9 z + C10 = 0
式中， C1 ~ C10 是方程的系数与常数项。
(2.13)
为了区别表面的朝向，还必须确定表面的正法向。对不透明表面，通常定义指向辐 射能束传递空间的法向为表面正法向。对式(2.14)描述的标准表面，表面正法向 n 为：
第2章 基本原理
2.5 表面辐射换热原理与计算方法简介
物体表面间的辐射换热是每个表面发射、吸收、反射辐射能的综合作用结果，取决于 每个表面的热辐射能发射能力、吸收能力、反射方式以及物体之间的相对空间几何关系[15]。 对辐射换热而言，物体的热辐射特性包括：热辐射发射特性、对投入辐射能的吸收、反射特 性三个方面[14,15]。 基尔霍夫定律表明：物体的光谱方向发射率等于其同一温度下的光谱方向吸收率，即：
面单元 A j 吸收的抽样能束数。 从上述介绍可看出，采用 MCM 进行辐射换热计算的关键在于建立系统内物体表面的 数学描述、各种表面的热辐射统计行为概率模型、能束抽样、跟踪与统计。
2.7 辐射能束与表面的数学描述
采用标准的二次方程来描述所有的表面，这种标准化处理对软件的通用性非常有 益
[21]
。在直角坐标系中，表面方程标准形式为：
x = x0 + m1 ⋅ t y = y0 + m2 ⋅ t z = z0 + m3 ⋅ t
其中，参数 t ≥ 0 ，代表辐射能束到达点与源点之间的距离。辐射能束的源点坐标 (2.17)
( x0 , y0 , z0 ) 由发射点概率模型或入射线与表面的交点确定，而方向矢量 m 的确定比较复 杂。对镜反射能束，根据 Fresnel 反射定律， m 由入射线的方向矢量 m 0 与表面的正法向矢
X d 1, d 2 =
X d 2, d 1 =
dA2 cosθ 1 cosθ 2
π r2 dA1 cosθ 1 cosθ 2
π r2
dA2
(2.7)
(2.8)
θ
n2
n1
dΩ
θ
r
dA1
图 2.4 两表面单元间的辐射换热角系数 对两个有限大小的表面，可在此基础上，通过积分求出二者的角系数，或利用角系数的 性质（相对性、完整性）通过代数分析求得。
量 n 通过矢量合成得到。经推导，镜反射辐射能束的方向矢量 m 为：
(2.18) m = m 0 − 2(n ⋅ m 0 )n 对漫发射、漫反射及各向异性发射和反射， m 由相应的概率模型确定。通常根据概率
再通过下式转换为方向矢量 模型得出光线传播方向在当地坐标系下的纬度角 θ 与经度角 ϕ ，
θ max ] 、 [ϕ min , ϕ max ] ，可推导出发射点的圆柱坐标系
= θ 0 arccos[cos θ min + Rθ ⋅ (cos θ max − cos θ min )]
ϕ 0 = ϕ min + Rϕ ⋅ (ϕ max − ϕ min )
(2来自百度文库22)
式中， Rθ 、 Rϕ 分别是球面上发射点沿极角与圆周方向的分布随机数。转换为直角坐 标形式为：
ε λ (T ,θ , ϕ ) = α λ (T ,θ , ϕ ) 。对漫灰体而言，则有： ε (T ) = α (T ) 。
在飞机内部物体表面之间的辐射换热中，将所涉及的物体表面均视为漫灰体表面。多 个漫灰体表面之间的辐射换热方法主要有： 净热量法、 Gebhart法、 网络法、 蒙特卡罗法[16]等。 网络法：利用电路来比拟辐射热流的传递路径，在投入辐射、有效辐射以及角系数概 念的基础上，将辐射换热系统模拟为由表面辐射热阻、空间辐射热阻构成的网络系统，通过 求表面的有效辐射，获得各表面的净辐射换热量[15-18]。 上述方法都要首先获得各表面之间的角系数。某一表面 Ai 发射的辐射能投射到另一表 对黑体或漫灰体表面， 角系数是纯几何因子， 面 A j 上的份额称为前者对后者的角系数 X i , j 。 仅反映两表面的几何特性及空间几何关系。 如图 2.1 所示， 根据兰贝特定律， 可推导出温度、 发射率分别均匀的两个表面微元 dA1 、 dA2 之间的角系数为：