传热学仿真-主题080-参与性介质辐射

kkchenjj

165人浏览 · 2026-03-21 04:41:42

kkchenjj · 2026-03-21 04:41:42 发布

第八十篇：参与性介质辐射

摘要

参与性介质辐射是指辐射在能够吸收、发射和散射的介质中传输的过程，是高温气体、燃烧产物、大气等介质的辐射传热基础。本教程系统介绍参与性介质辐射的基本理论，讨论辐射传输方程、吸收系数模型和数值求解方法。通过Python实现参与性介质的辐射传输计算，包括一维平板介质、球对称介质和复杂几何的辐射传热分析。本教程涵盖辐射传输方程、光学厚度、发射率计算、离散坐标法等核心内容。通过本教程的学习，读者将掌握参与性介质辐射的分析方法，能够进行高温气体和燃烧系统的辐射传热计算。

关键词

参与性介质，辐射传输方程，吸收系数，散射，光学厚度，离散坐标法，辐射传热

在这里插入图片描述

1. 参与性介质辐射基础

1.1 辐射传输方程

辐射传输方程（RTE）描述辐射在参与性介质中的传输：

$dIλds=−κλIλ−σsλIλ+κλIbλ+σsλ4π∫4πIλ(s^i)Φ(s^i,s^)dΩi\frac{dI_{\lambda}}{ds} = -\kappa_{\lambda} I_{\lambda} - \sigma_{s\lambda} I_{\lambda} + \kappa_{\lambda} I_{b\lambda} + \frac{\sigma_{s\lambda}}{4\pi} \int_{4\pi} I_{\lambda}(\hat{s}_i) \Phi(\hat{s}_i, \hat{s}) d\Omega_i$

其中：

$κλ\kappa_{\lambda}$ ：吸收系数
$σsλ\sigma_{s\lambda}$ ：散射系数
$βλ=κλ+σsλ\beta_{\lambda} = \kappa_{\lambda} + \sigma_{s\lambda}$ ：消光系数
$Φ\Phi$ ：散射相函数

1.2 光学厚度

光学厚度定义：
$τλ=∫0Lβλds\tau_{\lambda} = \int_0^L \beta_{\lambda} ds$

光学厚度分类：

光学薄： $τ≪1\tau \ll 1$
光学厚： $τ≫1\tau \gg 1$

1.3 发射率与吸收率

介质发射率：
$ε=1−e−τ\varepsilon = 1 - e^{-\tau}$

介质吸收率（基尔霍夫定律）：
$α=ε\alpha = \varepsilon$

2. 气体辐射特性

2.1 谱带模型

窄谱带模型：
$εˉ=1−exp⁡(−SdX)\bar{\varepsilon} = 1 - \exp\left(-\frac{S}{d} X\right)$

其中 $S / d$ 为平均吸收系数， $X$ 为光程。

宽谱带模型：
$Aˉ=A0ln⁡(1+uA0)\bar{A} = A_0 \ln\left(1 + \frac{u}{A_0}\right)$

2.2 全谱模型

加权求和灰气体模型（WSGGM）：
$ε=∑iai(T)(1−e−κipL)\varepsilon = \sum_{i} a_i(T) \left(1 - e^{-\kappa_i p L}\right)$

2.3 颗粒辐射

颗粒消光系数：
$β=Np⋅Cext=Np⋅πrp2Qext\beta = N_p \cdot C_{ext} = N_p \cdot \pi r_p^2 Q_{ext}$

其中 $Q_{ext}$ 为消光效率因子。

3. Python仿真实现

3.1 一维辐射传输

"""
主题080：参与性介质辐射
参与性介质辐射传输计算
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_bvp
import os

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

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

sigma = 5.670374419e-8

def planck_blackbody(T):
    """黑体辐射强度"""
    return sigma * T**4 / np.pi

def solve_1d_radiation_transfer(L, N, kappa, sigma_s, T_wall1, T_wall2, T_gas):
    """
    一维辐射传输方程求解（简化模型）
    L: 介质厚度
    N: 网格数
    kappa: 吸收系数
    sigma_s: 散射系数
    T_wall1, T_wall2: 壁面温度
    T_gas: 气体温度（假设均匀）
    """
    dx = L / (N - 1)
    x = np.linspace(0, L, N)
    beta = kappa + sigma_s
    omega = sigma_s / beta if beta > 0 else 0
    
    # 简化：使用扩散近似
    # d²G/dx² - 3κ(1-ω)G = -3κ(1-ω)4σT⁴
    
    # 构建矩阵
    A = np.zeros((N, N))
    b = np.zeros(N)
    
    # 内部节点
    for i in range(1, N-1):
        A[i, i-1] = 1 / dx**2
        A[i, i] = -2 / dx**2 - 3 * kappa * (1 - omega)
        A[i, i+1] = 1 / dx**2
        b[i] = -3 * kappa * (1 - omega) * 4 * sigma * T_gas**4
    
    # 边界条件（Marshak边界）
    eps1, eps2 = 1.0, 1.0  # 壁面发射率
    
    # 左边界
    A[0, 0] = 1
    b[0] = sigma * T_wall1**4
    
    # 右边界
    A[-1, -1] = 1
    b[-1] = sigma * T_wall2**4
    
    # 求解
    G = np.linalg.solve(A, b)
    
    # 计算辐射热流
    q = np.zeros(N)
    for i in range(1, N-1):
        q[i] = -1/3 * (G[i+1] - G[i-1]) / (2*dx)
    
    return x, G, q

def gas_emissivity(kappa, L, T):
    """
    气体发射率计算
    """
    tau = kappa * L
    eps = 1 - np.exp(-tau)
    return eps

def wsggm_emissivity(T, p_co2, p_h2o, L):
    """
    加权求和灰气体模型
    """
    # 简化模型：使用4个灰气体组分
    a = [0.3, 0.3, 0.2, 0.2]  # 权重
    k = [0.1, 0.5, 1.0, 2.0]  # 吸收系数 (1/m·atm)
    
    p_total = p_co2 + p_h2o
    eps = 0
    for i in range(len(a)):
        eps += a[i] * (1 - np.exp(-k[i] * p_total * L))
    
    return eps

# 1. 一维辐射传输
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

L = 1.0  # m
N = 50
kappa = 1.0  # 1/m
sigma_s = 0.5  # 1/m
T_wall1 = 1000  # K
T_wall2 = 500  # K
T_gas = 800  # K

x, G, q = solve_1d_radiation_transfer(L, N, kappa, sigma_s, T_wall1, T_wall2, T_gas)

# 计算温度
T_radiation = (G / sigma)**0.25

axes[0, 0].plot(x, T_radiation, 'b-', linewidth=2, label='Radiation Temperature')
axes[0, 0].axhline(y=T_gas, color='r', linestyle='--', linewidth=2, label=f'Gas T = {T_gas}K')
axes[0, 0].axhline(y=T_wall1, color='g', linestyle=':', linewidth=2, label=f'Wall1 T = {T_wall1}K')
axes[0, 0].axhline(y=T_wall2, color='orange', linestyle=':', linewidth=2, label=f'Wall2 T = {T_wall2}K')
axes[0, 0].set_xlabel('Position (m)', fontsize=11)
axes[0, 0].set_ylabel('Temperature (K)', fontsize=11)
axes[0, 0].set_title('1D Radiation Transfer in Participating Medium', fontsize=12)
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 2. 辐射热流分布
axes[0, 1].plot(x, q/1000, 'r-', linewidth=2)
axes[0, 1].set_xlabel('Position (m)', fontsize=11)
axes[0, 1].set_ylabel('Radiative Heat Flux (kW/m²)', fontsize=11)
axes[0, 1].set_title('Radiative Heat Flux Distribution', fontsize=12)
axes[0, 1].grid(True, alpha=0.3)

# 3. 光学厚度对发射率的影响
kappa_range = np.linspace(0.01, 5.0, 100)
L_values = [0.1, 0.5, 1.0, 2.0]
colors = ['blue', 'green', 'orange', 'red']

for L_val, color in zip(L_values, colors):
    eps_values = [gas_emissivity(k, L_val, 1000) for k in kappa_range]
    axes[1, 0].plot(kappa_range, eps_values, color=color, linewidth=2, 
                    label=f'L = {L_val}m')

axes[1, 0].set_xlabel('Absorption Coefficient κ (1/m)', fontsize=11)
axes[1, 0].set_ylabel('Gas Emissivity', fontsize=11)
axes[1, 0].set_title('Effect of Optical Thickness on Emissivity', fontsize=12)
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# 4. WSGGM模型
T_range = np.linspace(500, 2000, 100)
p_co2 = 0.1  # atm
p_h2o = 0.2  # atm
L = 1.0  # m

eps_wsggm = [wsggm_emissivity(T, p_co2, p_h2o, L) for T in T_range]

axes[1, 1].plot(T_range, eps_wsggm, 'b-', linewidth=2)
axes[1, 1].set_xlabel('Temperature (K)', fontsize=11)
axes[1, 1].set_ylabel('Gas Emissivity', fontsize=11)
axes[1, 1].set_title(f'WSGGM Emissivity (CO₂={p_co2}atm, H₂O={p_h2o}atm)', fontsize=12)
axes[1, 1].grid(True, alpha=0.3)

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

3.2 散射效应分析

"""
散射效应分析
"""

def scattering_effect(kappa, sigma_s, L, T_hot, T_cold):
    """
    分析散射对辐射传热的影响
    """
    beta = kappa + sigma_s
    omega = sigma_s / beta if beta > 0 else 0
    tau = beta * L
    
    # 纯吸收情况（无散射）
    q_absorption_only = sigma * (T_hot**4 - T_cold**4) * (1 - np.exp(-kappa * L))
    
    # 有散射情况（简化模型）
    # 散射使有效吸收减少
    kappa_eff = kappa * (1 - omega)
    q_with_scattering = sigma * (T_hot**4 - T_cold**4) * (1 - np.exp(-kappa_eff * L))
    
    return q_absorption_only, q_with_scattering, omega

# 参数分析
kappa = 1.0  # 1/m
L = 1.0  # m
T_hot = 1500  # K
T_cold = 500  # K

sigma_s_range = np.linspace(0, 5.0, 100)

q_abs = []
q_scat = []
omega_values = []

for sigma_s in sigma_s_range:
    q_a, q_s, omega = scattering_effect(kappa, sigma_s, L, T_hot, T_cold)
    q_abs.append(q_a/1000)
    q_scat.append(q_s/1000)
    omega_values.append(omega)

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

# 1. 散射对热流的影响
axes[0].plot(omega_values, q_scat, 'b-', linewidth=2, label='With Scattering')
axes[0].axhline(y=q_abs[0], color='r', linestyle='--', linewidth=2, label='Absorption Only')
axes[0].set_xlabel('Scattering Albedo ω', fontsize=11)
axes[0].set_ylabel('Radiative Heat Flux (kW/m²)', fontsize=11)
axes[0].set_title('Effect of Scattering on Heat Transfer', fontsize=12)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 2. 不同光学厚度下的发射率
tau_values = np.linspace(0.1, 10, 100)
eps_values = [1 - np.exp(-tau) for tau in tau_values]

axes[1].plot(tau_values, eps_values, 'g-', linewidth=2)
axes[1].axhline(y=0.632, color='r', linestyle='--', linewidth=2, label='τ=1 (ε=0.632)')
axes[1].axvline(x=1, color='r', linestyle='--', linewidth=2)
axes[1].set_xlabel('Optical Thickness τ', fontsize=11)
axes[1].set_ylabel('Emissivity ε', fontsize=11)
axes[1].set_title('Emissivity vs Optical Thickness', fontsize=12)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/scattering_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("散射效应分析图已保存")

print("\n" + "="*60)
print("仿真3：一维辐射传递方程求解")
print("="*60)

# 求解一维辐射传递方程（RTE）

# 介质参数
L_medium = 1.0  # m，介质厚度
kappa = 1.0  # 1/m，吸收系数
sigma_s = 0.5  # 1/m，散射系数
beta = kappa + sigma_s  # 消光系数
omega = sigma_s / beta  # 散射反照率

# 边界条件
T_wall1 = 1000  # K，左壁温度
T_wall2 = 500  # K，右壁温度
eps_wall = 0.8  # 壁面发射率

# 网格划分
N_z = 50
z = np.linspace(0, L_medium, N_z)
dz = z[1] - z[0]

# 离散方向（2流近似）
mu = [0.5, -0.5]  # 方向余弦
w = [0.5, 0.5]  # 权重

# 源项（黑体辐射）
T_gas = 800  # K，气体温度
sigma = 5.67e-8
I_b = sigma * T_gas**4 / np.pi  # 黑体辐射强度

# 求解辐射强度
I_plus = np.zeros(N_z)  # 正向
I_minus = np.zeros(N_z)  # 反向

# 边界条件
I_plus[0] = eps_wall * sigma * T_wall1**4 / np.pi + (1 - eps_wall) * I_minus[0]
I_minus[-1] = eps_wall * sigma * T_wall2**4 / np.pi + (1 - eps_wall) * I_plus[-1]

# 迭代求解
max_iter = 1000
tolerance = 1e-6

for iteration in range(max_iter):
    I_plus_old = I_plus.copy()
    I_minus_old = I_minus.copy()
    
    # 正向扫描
    for i in range(1, N_z):
        # 源项
        S = (1 - omega) * I_b + omega * (I_plus[i-1] + I_minus[i-1]) / 2
        # 离散方程
        I_plus[i] = (I_plus[i-1] * mu[0] + dz * beta * S) / (mu[0] + dz * beta)
    
    # 反向扫描
    for i in range(N_z-2, -1, -1):
        S = (1 - omega) * I_b + omega * (I_plus[i] + I_minus[i+1]) / 2
        I_minus[i] = (-I_minus[i+1] * mu[1] + dz * beta * S) / (-mu[1] + dz * beta)
    
    # 检查收敛
    error = np.max(np.abs(I_plus - I_plus_old)) + np.max(np.abs(I_minus - I_minus_old))
    if error < tolerance:
        print(f"RTE求解收敛，迭代次数: {iteration+1}")
        break

# 计算辐射热流
q_rad = np.pi * (I_plus - I_minus)

# 计算辐射能量密度
G = 2 * np.pi * (I_plus + I_minus)

print("一维辐射传递方程求解结果:")
print(f"介质厚度: {L_medium} m")
print(f"吸收系数: {kappa} m⁻¹")
print(f"散射系数: {sigma_s} m⁻¹")
print(f"散射反照率: {omega:.3f}")
print(f"\n边界条件:")
print(f"  左壁温度: {T_wall1} K")
print(f"  右壁温度: {T_wall2} K")
print(f"  气体温度: {T_gas} K")

print(f"\n辐射热流:")
print(f"  左壁: {q_rad[0]:.2f} W/m²")
print(f"  右壁: {q_rad[-1]:.2f} W/m²")
print(f"  中间: {q_rad[N_z//2]:.2f} W/m²")

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

# 辐射强度分布
axes[0].plot(z, I_plus, 'b-', linewidth=2, label='I⁺ (Forward)')
axes[0].plot(z, I_minus, 'r--', linewidth=2, label='I⁻ (Backward)')
axes[0].fill_between(z, I_plus, alpha=0.3, color='blue')
axes[0].fill_between(z, I_minus, alpha=0.3, color='red')
axes[0].set_xlabel('Position z (m)', fontsize=11)
axes[0].set_ylabel('Radiation Intensity (W/(m²·sr))', fontsize=11)
axes[0].set_title('Radiation Intensity Distribution', fontsize=12, fontweight='bold')
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)

# 辐射热流
axes[1].plot(z, q_rad, 'g-', linewidth=2)
axes[1].axhline(y=0, color='k', linestyle='--', linewidth=1)
axes[1].fill_between(z, q_rad, alpha=0.3, color='green')
axes[1].set_xlabel('Position z (m)', fontsize=11)
axes[1].set_ylabel('Radiative Heat Flux (W/m²)', fontsize=11)
axes[1].set_title('Radiative Heat Flux Distribution', fontsize=12, fontweight='bold')
axes[1].grid(True, alpha=0.3)

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

print("\n图3：一维辐射传递方程求解已保存")

print("\n" + "="*60)
print("仿真4：锅炉炉膛辐射换热")
print("="*60)

# 模拟锅炉炉膛内的辐射换热

# 炉膛几何参数
W_furnace = 10  # m，炉膛宽度
H_furnace = 15  # m，炉膛高度
D_furnace = 8  # m，炉膛深度

# 燃烧参数
Q_combustion = 50e6  # W，燃烧热功率
T_flame = 1800  # K，火焰温度
T_wall_furnace = 600  # K，水冷壁温度

# 烟气辐射特性（简化模型）
def gas_emissivity(T, L, P_CO2=0.1, P_H2O=0.1):
    """计算烟气发射率"""
    # Leckner简化模型
    P_total = P_CO2 + P_H2O
    C = 1.0  # 修正系数
    
    # 发射率计算
    eps_gas = C * (1 - np.exp(-0.5 * P_total * L * (1000/T)**0.5))
    return min(eps_gas, 1.0)

# 平均射线长度
L_mean = 0.9 * (4 * W_furnace * H_furnace * D_furnace / 
                (2*(W_furnace*H_furnace + W_furnace*D_furnace + H_furnace*D_furnace)))

print(f"平均射线长度: {L_mean:.2f} m")

# 烟气发射率
eps_gas = gas_emissivity(T_flame, L_mean)
print(f"烟气发射率: {eps_gas:.3f}")

# 水冷壁有效发射率（考虑积灰）
eps_wall_eff = 0.8  # 有效发射率

# 炉膛辐射换热计算
A_wall_total = 2 * (W_furnace*H_furnace + W_furnace*D_furnace + H_furnace*D_furnace)
sigma = 5.67e-8

# 火焰辐射
Q_rad_flame = eps_gas * sigma * A_wall_total * (T_flame**4 - T_wall_furnace**4)

# 水冷壁吸热
Q_wall_abs = eps_wall_eff * sigma * A_wall_total * (T_flame**4 - T_wall_furnace**4)

# 辐射热效率
eta_rad = Q_wall_abs / Q_combustion * 100

print("\n锅炉炉膛辐射换热分析:")
print(f"炉膛尺寸: {W_furnace}×{H_furnace}×{D_furnace} m")
print(f"燃烧功率: {Q_combustion/1e6:.1f} MW")
print(f"火焰温度: {T_flame} K")
print(f"水冷壁温度: {T_wall_furnace} K")
print(f"\n辐射换热:")
print(f"  火焰辐射: {Q_rad_flame/1e6:.2f} MW")
print(f"  水冷壁吸热: {Q_wall_abs/1e6:.2f} MW")
print(f"  辐射热效率: {eta_rad:.1f}%")

# 不同火焰温度的影响
T_flame_range = np.linspace(1500, 2200, 50)
Q_abs_vs_T = []
eta_vs_T = []

for T_f in T_flame_range:
    eps_g = gas_emissivity(T_f, L_mean)
    Q_abs = eps_wall_eff * sigma * A_wall_total * (T_f**4 - T_wall_furnace**4)
    Q_abs_vs_T.append(Q_abs / 1e6)  # MW
    eta_vs_T.append(Q_abs / Q_combustion * 100)

# 不同过量空气系数的影响（影响烟气成分）
excess_air_furnace = np.linspace(1.0, 1.5, 20)
Q_abs_vs_air = []

for lambda_air in excess_air_furnace:
    # 过量空气增加，CO2和H2O浓度降低
    P_CO2_eff = 0.12 / lambda_air
    P_H2O_eff = 0.11 / lambda_air
    eps_g = gas_emissivity(T_flame, L_mean, P_CO2_eff, P_H2O_eff)
    Q_abs = eps_wall_eff * sigma * A_wall_total * eps_g * (T_flame**4 - T_wall_furnace**4)
    Q_abs_vs_air.append(Q_abs / 1e6)

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

# 吸热量vs火焰温度
ax1 = axes[0]
ax1.plot(T_flame_range, Q_abs_vs_T, 'b-', linewidth=2, label='Heat Absorption')
ax1.fill_between(T_flame_range, Q_abs_vs_T, alpha=0.3, color='blue')
ax1.axvline(x=T_flame, color='r', linestyle='--', linewidth=2, label=f'Design: {T_flame}K')
ax1.set_xlabel('Flame Temperature (K)', fontsize=11)
ax1.set_ylabel('Heat Absorption (MW)', fontsize=11, color='b')
ax1.tick_params(axis='y', labelcolor='b')
ax1.legend(loc='upper left', fontsize=10)
ax1.grid(True, alpha=0.3)

ax2 = ax1.twinx()
ax2.plot(T_flame_range, eta_vs_T, 'r--', linewidth=2, label='Efficiency')
ax2.set_ylabel('Radiation Efficiency (%)', fontsize=11, color='r')
ax2.tick_params(axis='y', labelcolor='r')
ax2.legend(loc='lower right', fontsize=10)

axes[0].set_title('Furnace Heat Absorption vs Flame Temperature', fontsize=12, fontweight='bold')

# 吸热量vs过量空气
axes[1].plot(excess_air_furnace, Q_abs_vs_air, 'g-', linewidth=2)
axes[1].fill_between(excess_air_furnace, Q_abs_vs_air, alpha=0.3, color='green')
axes[1].set_xlabel('Excess Air Coefficient', fontsize=11)
axes[1].set_ylabel('Heat Absorption (MW)', fontsize=11)
axes[1].set_title('Effect of Excess Air on Heat Absorption', fontsize=12, fontweight='bold')
axes[1].grid(True, alpha=0.3)

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

print("\n图4：锅炉炉膛辐射换热已保存")
print("\n所有仿真完成！")