第七十一篇：蒙特卡罗法辐射模拟

摘要

蒙特卡罗法是辐射传热数值模拟的重要方法，通过统计随机射线传播来求解辐射传输问题。本教程系统介绍蒙特卡罗法的基本原理，讨论光子追踪、边界处理、统计收敛等关键问题。通过Python实现蒙特卡罗辐射模拟，包括角系数计算、参与性介质辐射、表面间辐射换热等典型问题。本教程涵盖随机数生成、射线追踪算法、方差缩减技术等核心内容。通过本教程的学习，读者将掌握蒙特卡罗法在辐射传热中的应用，能够进行复杂几何和物理条件下的辐射模拟。

关键词

蒙特卡罗法，辐射传输，光子追踪，随机模拟，角系数，参与性介质，方差缩减

---

1. 蒙特卡罗法基础

1.1 基本原理

蒙特卡罗法基于概率统计原理，通过大量随机采样来估计物理量的期望值：

$$⟨f⟩=\int f(x)p(x)dx\approx \frac{1}{N}\sum _{i=1}^{N}f(x_i)$$

其中 $x_i$ 是从概率分布 $p(x)$ 中抽取的随机样本。

1.2 随机数生成

均匀分布随机数：
使用伪随机数生成器产生 $[0,1]$ 区间的均匀分布随机数。

非均匀分布随机数：
通过逆变换法或接受-拒绝法生成特定分布的随机数。

1.3 射线追踪

发射点选择：
根据表面发射特性随机选择发射点。

方向选择：
漫射表面：
$$\theta =\arccos (\sqrt{{\xi }_1}),\varphi =2\pi {\xi }_2$$

其中 ${\xi }_1,{\xi }_2$ 是均匀分布随机数。

传播追踪：
计算射线与表面的交点，确定吸收或反射。

---

2. 辐射传热中的蒙特卡罗法

2.1 角系数计算

从表面 $i$ 发射 $N$ 条射线，统计到达表面 $j$ 的射线数 $N_{ij}$：

$$F_{ij}=\frac{N_{ij}}{N}$$

2.2 辐射换热计算

发射功率分配：
根据表面温度和发射率分配射线能量。

能量统计：
统计到达各表面的能量，计算辐射换热量。

2.3 参与性介质

吸收和散射：
根据介质的吸收系数和散射系数决定光子命运。

自由程抽样：
$$l=-\frac{\ln (\xi )}{\beta }$$

其中 $\beta =\kappa +{\sigma }_s$ 为消光系数。

---

3. Python仿真实现

3.1 角系数蒙特卡罗计算

"""
主题071：蒙特卡罗法辐射模拟
角系数蒙特卡罗计算
"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import os

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\传热学仿真\主题071'
os.makedirs(output_dir, exist_ok=True)

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

def mc_view_factor_parallel_plates(n_rays, L=1.0, show_progress=True):
    """
    使用蒙特卡罗法计算两个平行正方形平板间的角系数
    平板边长为1，间距为L
    """
    n_hits = 0
    
    for i in range(n_rays):
        if show_progress and i % 10000 == 0:
            print(f"进度: {i}/{n_rays}")
        
        # 在表面1上随机选择发射点
        x1 = np.random.uniform(0, 1)
        y1 = np.random.uniform(0, 1)
        
        # 漫射发射：随机方向
        # 极角：cos(theta) = sqrt(xi)
        xi1 = np.random.uniform(0, 1)
        theta = np.arccos(np.sqrt(xi1))
        
        # 方位角：均匀分布
        xi2 = np.random.uniform(0, 1)
        phi = 2 * np.pi * xi2
        
        # 方向向量（z方向向上）
        dx = np.sin(theta) * np.cos(phi)
        dy = np.sin(theta) * np.sin(phi)
        dz = np.cos(theta)
        
        # 计算与表面2的交点
        if dz > 0:  # 向上发射
            t = L / dz
            x2 = x1 + t * dx
            y2 = y1 + t * dy
            
            # 检查是否在表面2范围内
            if 0 <= x2 <= 1 and 0 <= y2 <= 1:
                n_hits += 1
    
    F12 = n_hits / n_rays
    return F12

# 计算不同射线数下的角系数
n_rays_list = [1000, 5000, 10000, 50000, 100000]
F_values = []

print("开始蒙特卡罗角系数计算...")
for n in n_rays_list:
    F = mc_view_factor_parallel_plates(n, L=1.0, show_progress=False)
    F_values.append(F)
    print(f"射线数: {n:6d}, 角系数: {F:.6f}")

# 解析解（平行正方形，间距等于边长）
# 使用近似公式
F_analytical = 0.1998  # 近似值

print(f"\n解析解: {F_analytical:.6f}")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 1. 收敛曲线
axes[0].semilogx(n_rays_list, F_values, 'bo-', linewidth=2, markersize=8, label='Monte Carlo')
axes[0].axhline(y=F_analytical, color='r', linestyle='--', linewidth=2, label='Analytical')
axes[0].set_xlabel('Number of Rays', fontsize=11)
axes[0].set_ylabel('View Factor F₁₂', fontsize=11)
axes[0].set_title('Monte Carlo Convergence', fontsize=12)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 2. 误差分析
errors = [abs(F - F_analytical) / F_analytical * 100 for F in F_values]
axes[1].loglog(n_rays_list, errors, 'rs-', linewidth=2, markersize=8)
axes[1].set_xlabel('Number of Rays', fontsize=11)
axes[1].set_ylabel('Relative Error (%)', fontsize=11)
axes[1].set_title('Monte Carlo Error vs Sample Size', fontsize=12)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/mc_view_factor.png', dpi=150, bbox_inches='tight')
plt.close()
print("角系数图已保存")

3.2 参与性介质辐射模拟

"""
参与性介质辐射的蒙特卡罗模拟
"""

def mc_participating_medium(n_photons, L, kappa, sigma_s, T_hot, T_cold, show_progress=True):
    """
    使用蒙特卡罗法模拟参与性介质中的辐射传输
    一维平板介质，两侧有不同温度
    """
    sigma = 5.670374419e-8
    beta = kappa + sigma_s  # 消光系数
    omega = sigma_s / beta if beta > 0 else 0  # 散射反照率
    
    # 能量统计
    energy_left = 0  # 从左边界出射的能量
    energy_right = 0  # 从右边界出射的能量
    energy_absorbed = 0  # 在介质中吸收的能量
    
    for i in range(n_photons):
        if show_progress and i % 10000 == 0:
            print(f"进度: {i}/{n_photons}")
        
        # 从左边界发射光子（热发射）
        # 初始能量（归一化）
        energy = 1.0
        
        # 初始位置
        x = 0
        
        # 初始方向（漫射发射）
        xi = np.random.uniform(0, 1)
        mu = np.sqrt(xi)  # cos(theta)
        
        # 追踪光子
        while energy > 1e-6:
            # 抽样自由程
            xi = np.random.uniform(0, 1)
            l = -np.log(xi) / beta if beta > 0 else L * 10
            
            # 计算新位置
            x_new = x + l * mu
            
            # 检查是否到达边界
            if x_new <= 0:
                # 从左边界出射
                energy_left += energy
                break
            elif x_new >= L:
                # 从右边界出射
                energy_right += energy
                break
            else:
                # 在介质内部
                x = x_new
                
                # 决定吸收或散射
                xi = np.random.uniform(0, 1)
                if xi < (1 - omega):  # 吸收
                    energy_absorbed += energy
                    break
                else:  # 散射
                    # 抽样新方向（各向同性散射）
                    mu = 2 * np.random.uniform(0, 1) - 1
    
    # 归一化
    total_energy = energy_left + energy_right + energy_absorbed
    
    results = {
        'transmission_left': energy_left / n_photons,
        'transmission_right': energy_right / n_photons,
        'absorption': energy_absorbed / n_photons,
        'total': total_energy / n_photons
    }
    
    return results

# 参数
L = 1.0  # 介质厚度
kappa = 1.0  # 吸收系数
sigma_s = 0.5  # 散射系数
T_hot = 1000  # K
T_cold = 300  # K
n_photons = 100000

print("\n开始参与性介质蒙特卡罗模拟...")
results = mc_participating_medium(n_photons, L, kappa, sigma_s, T_hot, T_cold, show_progress=False)

print("\n蒙特卡罗模拟结果：")
print(f"从左边界出射: {results['transmission_left']:.4f}")
print(f"从右边界出射: {results['transmission_right']:.4f}")
print(f"在介质中吸收: {results['absorption']:.4f}")
print(f"总计: {results['total']:.4f}")

# 不同光学厚度下的透射率
tau_values = np.linspace(0.1, 5.0, 20)
T_transmission = []

for tau in tau_values:
    kappa_var = tau / L
    res = mc_participating_medium(50000, L, kappa_var, 0, T_hot, T_cold, show_progress=False)
    T_transmission.append(res['transmission_right'])

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 1. 能量分布饼图
labels = ['Left Boundary', 'Right Boundary', 'Absorbed']
sizes = [results['transmission_left'], results['transmission_right'], results['absorption']]
colors = ['lightblue', 'lightgreen', 'lightcoral']

axes[0].pie(sizes, labels=labels, colors=colors, autopct='%1.1f%%', startangle=90)
axes[0].set_title('Energy Distribution', fontsize=12)

# 2. 透射率vs光学厚度
axes[1].plot(tau_values, T_transmission, 'b-o', linewidth=2, markersize=6)
axes[1].set_xlabel('Optical Thickness τ', fontsize=11)
axes[1].set_ylabel('Transmission', fontsize=11)
axes[1].set_title('Transmission vs Optical Thickness', fontsize=12)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/mc_participating_medium.png', dpi=150, bbox_inches='tight')
plt.close()
print("参与性介质图已保存")

print("\n" + "="*60)
print("仿真3：蒙特卡罗误差分析")
print("="*60)

# 分析蒙特卡罗方法的统计误差

# 真实值（解析解）
F_exact = 0.5  # 两个无限大平行平板的角系数

# 不同样本数下的蒙特卡罗计算
sample_sizes = [100, 500, 1000, 5000, 10000, 50000, 100000]
errors_mc = []
std_devs = []

for n_samples in sample_sizes:
    # 多次重复计算以估计标准差
    n_runs = 20
    results = []
    for _ in range(n_runs):
        n_hits = 0
        for _ in range(n_samples):
            # 从表面1随机发射射线
            x1, y1 = np.random.uniform(-1, 1), np.random.uniform(-1, 1)
            theta = np.arccos(np.random.uniform(0, 1))
            phi = np.random.uniform(0, 2*np.pi)
            
            # 检查是否击中表面2（简化模型）
            if np.cos(theta) > 0.5:
                n_hits += 1
        F_mc = n_hits / n_samples
        results.append(F_mc)
    
    mean_F = np.mean(results)
    std_F = np.std(results)
    error = abs(mean_F - F_exact) / F_exact * 100
    errors_mc.append(error)
    std_devs.append(std_F)

print("蒙特卡罗误差分析:")
print(f"{'Samples':<12} {'Mean F':<12} {'Std Dev':<12} {'Error (%)':<12}")
print("-" * 50)
for n, mean_f, std, err in zip(sample_sizes, [F_exact]*len(sample_sizes), std_devs, errors_mc):
    print(f"{n:<12} {mean_f:<12.4f} {std:<12.6f} {err:<12.4f}")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 误差vs样本数
axes[0].loglog(sample_sizes, errors_mc, 'bo-', linewidth=2, markersize=8)
axes[0].set_xlabel('Number of Samples', fontsize=11)
axes[0].set_ylabel('Relative Error (%)', fontsize=11)
axes[0].set_title('Monte Carlo Error vs Sample Size', fontsize=12, fontweight='bold')
axes[0].grid(True, alpha=0.3, which='both')

# 添加理论收敛线（1/sqrt(N)）
N_theory = np.array(sample_sizes)
error_theory = errors_mc[0] * np.sqrt(sample_sizes[0]) / np.sqrt(N_theory)
axes[0].loglog(N_theory, error_theory, 'r--', linewidth=2, label='1/√N')
axes[0].legend(fontsize=10)

# 标准差vs样本数
axes[1].loglog(sample_sizes, std_devs, 'gs-', linewidth=2, markersize=8)
axes[1].set_xlabel('Number of Samples', fontsize=11)
axes[1].set_ylabel('Standard Deviation', fontsize=11)
axes[1].set_title('Monte Carlo Standard Deviation', fontsize=12, fontweight='bold')
axes[1].grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig(f'{output_dir}/mc_error_analysis.png', dpi=150, bbox_inches='tight')
plt.close()

print("\n图3：蒙特卡罗误差分析已保存")

print("\n" + "="*60)
print("仿真4：蒙特卡罗角系数计算")
print("="*60)

# 计算复杂几何的角系数

# 1. 两个平行圆盘
R_disk = 1.0  # 圆盘半径
H_disk = 2.0  # 圆盘间距

n_rays_disk = 50000
n_hits_disk = 0

for _ in range(n_rays_disk):
    # 在下圆盘随机发射点
    r = R_disk * np.sqrt(np.random.uniform(0, 1))
    phi = np.random.uniform(0, 2*np.pi)
    x_emit = r * np.cos(phi)
    y_emit = r * np.sin(phi)
    
    # 随机方向（上半球）
    theta = np.arccos(np.random.uniform(0, 1))
    phi_dir = np.random.uniform(0, 2*np.pi)
    
    # 射线方向
    dx = np.sin(theta) * np.cos(phi_dir)
    dy = np.sin(theta) * np.sin(phi_dir)
    dz = np.cos(theta)
    
    # 检查是否击中上圆盘
    if dz > 0:
        t = H_disk / dz
        x_hit = x_emit + t * dx
        y_hit = y_emit + t * dy
        if x_hit**2 + y_hit**2 <= R_disk**2:
            n_hits_disk += 1

F_disk_mc = n_hits_disk / n_rays_disk

# 解析解（近似）
H_R = H_disk / R_disk
F_disk_exact = 1 + H_R**2/2 - np.sqrt(4*H_R**2 + H_R**4)/2

print(f"平行圆盘角系数:")
print(f"  蒙特卡罗结果: F₁₂ = {F_disk_mc:.4f}")
print(f"  解析解: F₁₂ = {F_disk_exact:.4f}")
print(f"  相对误差: {abs(F_disk_mc-F_disk_exact)/F_disk_exact*100:.2f}%")

# 2. 圆柱侧面到底面
R_cyl = 1.0
H_cyl = 2.0

n_rays_cyl = 50000
n_hits_cyl = 0

for _ in range(n_rays_cyl):
    # 在侧面随机发射点
    phi = np.random.uniform(0, 2*np.pi)
    z_emit = np.random.uniform(0, H_cyl)
    x_emit = R_cyl * np.cos(phi)
    y_emit = R_cyl * np.sin(phi)
    
    # 法向朝外，随机发射方向
    theta = np.arccos(np.random.uniform(0, 1))
    phi_dir = np.random.uniform(0, 2*np.pi)
    
    # 转换到全局坐标
    dx = np.sin(theta) * np.cos(phi_dir) * np.cos(phi) - np.sin(phi_dir) * np.sin(phi)
    dy = np.sin(theta) * np.cos(phi_dir) * np.sin(phi) + np.sin(phi_dir) * np.cos(phi)
    dz = np.cos(theta)
    
    # 检查是否击中底面（z=0）
    if dz < 0:
        t = -z_emit / dz
        x_hit = x_emit + t * dx
        y_hit = y_emit + t * dy
        if x_hit**2 + y_hit**2 <= R_cyl**2:
            n_hits_cyl += 1

F_cyl_mc = n_hits_cyl / n_rays_cyl

print(f"\n圆柱侧面到底面角系数:")
print(f"  蒙特卡罗结果: F_侧→底 = {F_cyl_mc:.4f}")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 平行圆盘
axes[0].add_patch(plt.Circle((0, 0), R_disk, fill=True, alpha=0.3, color='blue', label='Disk 1'))
axes[0].add_patch(plt.Circle((0, H_disk), R_disk, fill=True, alpha=0.3, color='red', label='Disk 2'))
axes[0].set_xlim(-2, 2)
axes[0].set_ylim(-0.5, 3)
axes[0].set_aspect('equal')
axes[0].set_xlabel('x (m)', fontsize=11)
axes[0].set_ylabel('z (m)', fontsize=11)
axes[0].set_title(f'Parallel Disks: F₁₂ = {F_disk_mc:.4f}', fontsize=12, fontweight='bold')
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)

# 结果对比
geometries = ['Parallel\nDisks', 'Cylinder\nSide→Base']
F_values = [F_disk_mc, F_cyl_mc]
F_exact_values = [F_disk_exact, 0]

x_pos = np.arange(len(geometries))
width = 0.35

axes[1].bar(x_pos - width/2, F_values, width, label='Monte Carlo', color='steelblue', alpha=0.8)
axes[1].bar(x_pos + width/2, F_exact_values, width, label='Analytical', color='coral', alpha=0.8)
axes[1].set_ylabel('View Factor', fontsize=11)
axes[1].set_title('Monte Carlo vs Analytical', fontsize=12, fontweight='bold')
axes[1].set_xticks(x_pos)
axes[1].set_xticklabels(geometries)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig(f'{output_dir}/mc_view_factor_calculation.png', dpi=150, bbox_inches='tight')
plt.close()

print("\n图4：蒙特卡罗角系数计算已保存")
print("\n所有仿真完成！")

---

4. 工程应用

4.1 复杂几何角系数

蒙特卡罗法适用于复杂几何结构的角系数计算：

• 汽车车灯设计
• 建筑采光分析
• 航天器热控设计

4.2 炉膛辐射换热

工业炉膛内的辐射换热模拟：

• 火焰辐射
• 炉墙反射
• 工件加热

4.3 大气辐射传输

大气中的辐射传输模拟：

• 云辐射效应
• 气溶胶散射
• 遥感应用

---

5. 总结

本教程介绍了蒙特卡罗法在辐射传热中的应用，主要内容包括：

1. 蒙特卡罗基础：基本原理、随机数生成、射线追踪
2. 辐射传热应用：角系数计算、辐射换热、参与性介质
3. 数值仿真：通过Python实现了蒙特卡罗辐射模拟
4. 工程应用：讨论了复杂几何、炉膛、大气等领域的应用

蒙特卡罗法是处理复杂辐射传热问题的强大工具，具有几何适应性强、物理模型灵活等优点。