第二十八篇：污染物生成机理

摘要

燃烧污染物生成机理是理解和控制燃烧过程中有害物质排放的关键科学基础。本教程系统介绍燃烧过程中主要污染物（NOx、SOx、CO、颗粒物等）的生成机理、影响因素和控制策略。内容涵盖热力型NOx、快速型NOx、燃料型NOx的形成机制，SOx的氧化过程，CO的生成与氧化，碳烟的形成与演化等核心理论。通过Python仿真实验，读者将深入理解污染物生成的化学动力学过程，掌握NOx生成模型、碳烟形成模拟、排放预测等关键技术，为进行低污染燃烧系统的设计和优化奠定坚实的理论基础。

关键词

污染物生成，NOx，SOx，CO，碳烟，热力型NO，快速型NO，燃料型NO，Zeldovich机理

---

1. 燃烧污染物概述

1.1 主要污染物类型

氮氧化物（NOx）：

• NO（一氧化氮）：主要组分，占NOx的90-95%
• NO₂（二氧化氮）：次要组分，占5-10%
• N₂O（一氧化二氮）：温室气体，痕量

硫氧化物（SOx）：

• SO₂（二氧化硫）：主要组分
• SO₃（三氧化硫）：次要组分，形成硫酸雾

一氧化碳（CO）：

• 不完全燃烧产物
• 有毒气体

颗粒物（PM）：

• 碳烟（Soot）
• 灰分（Ash）
• 有机颗粒物

未燃碳氢化合物（UHC）：

• 燃料不完全燃烧
• 壁面淬熄效应

1.2 污染物形成的影响因素

温度影响：

• NOx：指数增长（Zeldovich机理）
• CO：高温有利于氧化
• 碳烟：高温促进生成但加速氧化

当量比影响：

• 贫燃：NOx高，CO低
• 富燃：NOx低，CO高，碳烟高
• 化学计量比附近：复杂竞争

停留时间影响：

• 长停留时间：有利于NOx生成
• 短停留时间：有利于降低NOx

混合程度影响：

• 均匀混合：NOx高
• 分层混合：可降低NOx

---

2. NOx生成机理

2.1 热力型NO（Thermal NO）

Zeldovich机理（扩展）：

$$N_2+O\rightleftharpoons NO+N(R1)$$
$$N+O_2\rightleftharpoons NO+O(R2)$$
$$N+OH\rightleftharpoons NO+H(R3)$$

反应速率常数（cm³/mol/s）：

$$k_{f1}=1.8\times 10^{14}\exp (-38370/T)$$
$$k_{f2}=1.8\times 10^{10}T\exp (-4680/T)$$
$$k_{f3}=7.1\times 10^{13}\exp (-450/T)$$

NO生成速率：

$$\frac{d[NO]}{dt}=2k_{f1}[N_2][O]$$

温度依赖性：

• 活化能高（约75 kcal/mol）
• 温度>1800K时显著生成
• 温度每升高100K，NOx增加约3-5倍

2.2 快速型NO（Prompt NO）

Fenimore机理：

$$CH+N_2\rightleftharpoons HCN+N(R4)$$
$$C+N_2\rightleftharpoons CN+N(R5)$$

HCN氧化路径：

$$HCN+O\to NCO+H$$
$$NCO+H\to NH+CO$$
$$NH+H\to N+H_2$$
$$N+OH\to NO+H$$

特征：

• 在火焰前锋快速生成
• 与碳氢自由基有关
• 低温区（<1800K）也可生成
• 占总NOx的10-30%

2.3 燃料型NO（Fuel NO）

燃料氮转化路径：

1. 挥发分氮释放：
   $$Fuel-N\to HCN+N{H}_3+N_2$$

2. HCN氧化：
   $$HCN+O\to NCO+H$$
   $$NCO+O\to NO+CO$$

3. NH₃氧化：
   $$N{H}_3+OH\to N{H}_2+H_{2}O$$
   $$N{H}_2+O\to NH+OH$$
   $$NH+O\to N+OH$$
   $$N+OH\to NO+H$$

影响因素：

• 燃料氮含量
• 挥发分/焦炭比例
• 氧浓度
• 温度

2.4 NOx还原机理

再燃（Reburning）：

$$NO+C{H}_i\to HCN+\dots$$
$$NO+CO\to \frac{1}{2}{N}_2+C{O}_2$$

SNCR（选择性非催化还原）：

$$NO+N{H}_3+\frac{1}{4}{O}_2\to N_2+\frac{3}{2}{H}_{2}O$$

SCR（选择性催化还原）：

$$NO+N{H}_3\stackrel{catalyst}{\to }{N}_2+H_{2}O$$

---

3. SOx生成机理

3.1 燃料硫氧化

主要反应：

$$S+O_2\to S{O}_2$$
$$H_{2}S+\frac{3}{2}{O}_2\to S{O}_2+H_{2}O$$
$$COS+\frac{3}{2}{O}_2\to S{O}_2+C{O}_2$$

SO₂向SO₃转化：

$$S{O}_2+O+M\to S{O}_3+M$$
$$S{O}_2+OH+M\to HOS{O}_2+M$$
$$HOS{O}_2+O_2\to H{O}_2+S{O}_3$$

3.2 SOx排放特性

燃料硫含量与SO₂排放：

$$[S{O}_2{]}_{ppm}=\frac{S_{fuel}\times 10^6}{32\times (1+\frac{A}{F})\times M_{exhaust}}$$

其中：

• $S_{fuel}$：燃料硫质量分数
• $A/F$：空燃比
• $M_{exhaust}$：排气分子量

---

4. CO生成与氧化

4.1 CO生成机理

不完全燃烧来源：

• 富燃区缺氧
• 低温区反应冻结
• 混合不均匀

主要生成路径：

$$CH+O\to CO+H$$
$$C_{2}H+O\to 2CO+H$$
$$HCO+M\to CO+H+M$$

4.2 CO氧化机理

主要反应：

$$CO+OH\to C{O}_2+H(R6)$$
$$CO+O_2\to C{O}_2+O(R7)$$
$$CO+H{O}_2\to C{O}_2+OH(R8)$$

速率常数：

$$k_6=2.76\times 10^4{T}^{1.5}\exp (116/T)$$

CO平衡浓度：

$$[CO{]}_{eq}=K_p\frac{[C][O]}{[C{O}_2]}$$

---

5. 碳烟生成机理

5.1 碳烟形成过程

四个阶段：

1. 成核（Nucleation）：

   • 多环芳烃（PAH）聚合
   • 形成初始颗粒（~1-2 nm）

2. 表面生长（Surface Growth）：

   • 乙炔加成
   • 芳香族化合物沉积

3. 凝聚（Coagulation）：

   • 颗粒碰撞合并
   • 数量减少，尺寸增大

4. 氧化（Oxidation）：

   • OH、O₂氧化
   • 颗粒烧蚀

5.2 碳烟生成模型

经验模型 - Hiroyasu模型：

$$\frac{d{m}_{soot}}{dt}=A_{sf}{m}_{fuel}^{0.5}{P}^{1.8}\exp (-\frac{E_a}{RT})-A_{sc}{m}_{soot}(\frac{P_{O_2}}{P}{)}^{1.8}\exp (-\frac{E_{sc}}{RT})$$

半经验模型 - Moss模型：

$$\frac{d{Y}_{soot}}{dt}=C_f{Y}_{fuel}^n\exp (-\frac{T_a}{T})-C_o{Y}_{soot}(\frac{Y_{O_2}}{Y_{O_2,ref}}{)}^m$$

详细机理模型：

• 基于PAH化学
• 包含数百个反应
• 计算成本高

5.3 碳烟颗粒动力学

颗粒数密度方程：

$$\frac{\partial N}{\partial t}+\nabla \cdot (\mathbf{u}N)={\dot{\omega }}_{nucleation}+{\dot{\omega }}_{coagulation}$$

颗粒质量方程：

$$\frac{\partial M}{\partial t}+\nabla \cdot (\mathbf{u}M)={\dot{\omega }}_{growth}+{\dot{\omega }}_{oxidation}$$

方法 of Moments：

$$M_r={\int }_0^{\infty }{d}^{r}n(d)dd$$

---

6. 污染物控制策略

6.1 燃烧中控制

降低热力型NOx：

• 降低峰值温度
• 分级燃烧
• 烟气再循环
• 富氧燃烧

降低快速型NOx：

• 快速混合
• 降低碳氢自由基浓度

降低燃料型NOx：

• 燃料预处理
• 再燃技术

6.2 燃烧后控制

SNCR：

• 温度窗口：850-1100°C
• 还原剂：氨水或尿素
• 脱硝效率：30-50%

SCR：

• 温度窗口：300-400°C
• 催化剂：V₂O₅-WO₃/TiO₂
• 脱硝效率：80-95%

湿法脱硫：

• 石灰石-石膏法
• 海水脱硫
• 氨法脱硫

6.3 碳烟控制

燃烧优化：

• 改善混合
• 提高燃烧温度
• 延长停留时间

后处理：

• 柴油颗粒过滤器（DPF）
• 再生技术
• 催化氧化

---

7. Python仿真实验

7.1 实验1：热力型NOx生成模拟

实验目的：基于Zeldovich机理模拟热力型NOx的生成过程。

Python代码：

import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import os

output_dir = r'd:\文档\燃烧仿真\主题028'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("="*60)
print("主题028：污染物生成机理")
print("="*60)

# ============================================================
# 实验1：热力型NOx生成模拟
# ============================================================
print("\n实验1：热力型NOx生成模拟")
print("-" * 40)

# Zeldovich机理参数
# 反应1: N2 + O <=> NO + N
k_f1 = lambda T: 1.8e14 * np.exp(-38370/T)  # cm^3/mol/s
k_r1 = lambda T: 3.8e13  # cm^3/mol/s

# 反应2: N + O2 <=> NO + O
k_f2 = lambda T: 1.8e10 * T * np.exp(-4680/T)
k_r2 = lambda T: 3.8e7 * T * np.exp(-20820/T)

# 反应3: N + OH <=> NO + H
k_f3 = lambda T: 7.1e13 * np.exp(-450/T)
k_r3 = lambda T: 4.1e13 * np.exp(-20460/T)

# 氧原子平衡浓度（基于O2解离平衡）
def O_equilibrium(T, P_atm=1.0):
    """计算氧原子平衡浓度"""
    # O2 <=> 2O 的平衡常数
    Kp_O2 = 3.0e21 * T**(-0.5) * np.exp(-59360/T)  # atm
    X_O2 = 0.21 * (1 - 0.1)  # 假设10%燃料消耗
    X_O = np.sqrt(Kp_O2 * X_O2)
    
    # 转换为mol/cm^3
    R_atm = 82.057  # cm^3 atm / mol / K
    c_total = P_atm / (R_atm * T)
    return X_O * c_total

# OH平衡浓度
def OH_equilibrium(T, P_atm=1.0, phi=1.0):
    """计算OH平衡浓度"""
    # 简化模型
    if phi < 1.0:
        X_OH = 1e-3 * (T/2000)**2 * np.exp(-5000/T)
    else:
        X_OH = 5e-4 * (T/2000)**2 * np.exp(-3000/T)
    
    R_atm = 82.057
    c_total = P_atm / (R_atm * T)
    return X_OH * c_total

# NO生成ODE
def NO_formation(y, t, T, P_atm):
    NO, N = y
    
    # 平衡浓度
    O = O_equilibrium(T, P_atm)
    OH = OH_equilibrium(T, P_atm)
    
    # 假设N2和O2浓度基本不变
    R_atm = 82.057
    c_total = P_atm / (R_atm * T)
    N2 = 0.79 * c_total
    O2 = 0.21 * c_total * 0.9
    
    # 反应速率
    r1_f = k_f1(T) * N2 * O
    r1_r = k_r1(T) * NO * N
    r2_f = k_f2(T) * N * O2
    r2_r = k_r2(T) * NO * O
    r3_f = k_f3(T) * N * OH
    r3_r = k_r3(T) * NO * 1e-10  # H浓度假设很小
    
    # NO和N的变化率
    dNOdt = r1_f - r1_r + r2_f - r2_r + r3_f - r3_r
    dNdt = r1_f - r1_r - r2_f + r2_r - r3_f + r3_r
    
    return [dNOdt, dNdt]

# 模拟不同温度下的NO生成
T_range = [1600, 1800, 2000, 2200, 2400]  # K
t_max = 0.1  # s
t = np.linspace(0, t_max, 1000)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. NO生成曲线
ax1 = axes[0, 0]
colors = plt.cm.viridis(np.linspace(0, 1, len(T_range)))

NO_final = []
for i, T in enumerate(T_range):
    y0 = [1e-12, 1e-15]  # 初始NO和N浓度
    solution = odeint(NO_formation, y0, t, args=(T, 1.0))
    NO_conc = solution[:, 0]
    
    # 转换为ppm
    R_atm = 82.057
    c_total = 1.0 / (R_atm * T)
    NO_ppm = NO_conc / c_total * 1e6
    
    ax1.semilogy(t*1000, NO_ppm, linewidth=2.5, color=colors[i], label=f'T={T}K')
    NO_final.append(NO_ppm[-1])

ax1.set_xlabel('Time (ms)', fontsize=11)
ax1.set_ylabel('NO Concentration (ppm)', fontsize=11)
ax1.set_title('Thermal NO Formation', fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 2. 温度对最终NO浓度的影响
ax2 = axes[0, 1]
ax2.semilogy(T_range, NO_final, 'ro-', linewidth=2.5, markersize=10)
ax2.set_xlabel('Temperature (K)', fontsize=11)
ax2.set_ylabel('Final NO Concentration (ppm)', fontsize=11)
ax2.set_title('Temperature Effect on NO', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 3. 活化能分析
ax3 = axes[1, 0]
inv_T = 1000 / np.array(T_range)
ln_NO = np.log(NO_final)
coeffs = np.polyfit(inv_T, ln_NO, 1)
Ea_eff = -coeffs[0] * 8.314  # J/mol

ax3.plot(inv_T, ln_NO, 'bo', markersize=10, label='Simulation')
ax3.plot(inv_T, np.polyval(coeffs, inv_T), 'r--', linewidth=2, 
         label=f'Ea={Ea_eff/1000:.1f} kJ/mol')
ax3.set_xlabel('1000/T (K^-1)', fontsize=11)
ax3.set_ylabel('ln(NO)', fontsize=11)
ax3.set_title('Arrhenius Plot for NO Formation', fontsize=12, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 4. 氧原子浓度
ax4 = axes[1, 1]
T_plot = np.linspace(1500, 2500, 100)
O_conc = [O_equilibrium(T) for T in T_plot]

ax4.semilogy(1000/T_plot, O_conc, 'b-', linewidth=2.5)
ax4.set_xlabel('1000/T (K^-1)', fontsize=11)
ax4.set_ylabel('O Atom Concentration (mol/cm^3)', fontsize=11)
ax4.set_title('O Atom Equilibrium Concentration', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)
ax4.invert_xaxis()

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

print("✓ 实验1完成：热力型NOx生成模拟")
print(f"  有效活化能: {Ea_eff/1000:.1f} kJ/mol")

# ============================================================
# 实验2：NOx生成路径分析
# ============================================================
print("\n实验2：NOx生成路径分析")
print("-" * 40)

# 三种NOx生成机理的贡献
T_analysis = np.linspace(1200, 2400, 50)

# 热力型NO生成速率
def thermal_NO_rate(T, P_atm=1.0):
    O = O_equilibrium(T, P_atm)
    R_atm = 82.057
    c_total = P_atm / (R_atm * T)
    N2 = 0.79 * c_total
    return 2 * k_f1(T) * N2 * O  # mol/cm^3/s

# 快速型NO生成速率（简化模型）
def prompt_NO_rate(T, phi=1.0):
    # 与CH自由基和N2反应相关
    if phi >= 1.0:
        # 富燃区快速型NO显著
        return 1e-12 * np.exp(-3000/T) * (phi - 0.8)**2
    else:
        return 1e-15

# 燃料型NO生成速率（简化模型）
def fuel_NO_rate(T, N_fuel=0.01):
    # 假设燃料氮含量1%
    return 1e-10 * N_fuel * np.exp(-5000/T)

thermal_rate = [thermal_NO_rate(T) for T in T_analysis]
prompt_rate = [prompt_NO_rate(T, phi=1.2) for T in T_analysis]
fuel_rate = [fuel_NO_rate(T, N_fuel=0.005) for T in T_analysis]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. 三种机理的生成速率
ax1 = axes[0, 0]
ax1.semilogy(T_analysis, thermal_rate, 'r-', linewidth=2.5, label='Thermal NO')
ax1.semilogy(T_analysis, prompt_rate, 'b-', linewidth=2.5, label='Prompt NO')
ax1.semilogy(T_analysis, fuel_rate, 'g-', linewidth=2.5, label='Fuel NO')
ax1.set_xlabel('Temperature (K)', fontsize=11)
ax1.set_ylabel('Formation Rate (mol/cm^3/s)', fontsize=11)
ax1.set_title('NOx Formation Mechanisms', fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 2. 各机理贡献比例（归一化）
ax2 = axes[0, 1]
total_rate = np.array(thermal_rate) + np.array(prompt_rate) + np.array(fuel_rate)
thermal_frac = np.array(thermal_rate) / total_rate
prompt_frac = np.array(prompt_rate) / total_rate
fuel_frac = np.array(fuel_rate) / total_rate

ax2.stackplot(T_analysis, thermal_frac, prompt_frac, fuel_frac,
              labels=['Thermal', 'Prompt', 'Fuel'],
              colors=['red', 'blue', 'green'], alpha=0.7)
ax2.set_xlabel('Temperature (K)', fontsize=11)
ax2.set_ylabel('Fraction', fontsize=11)
ax2.set_title('NOx Formation Mechanism Contributions', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10, loc='center right')
ax2.set_xlim([1200, 2400])
ax2.set_ylim([0, 1])

# 3. 当量比对快速型NO的影响
ax3 = axes[1, 0]
phi_range = np.linspace(0.8, 1.5, 50)
prompt_vs_phi = [prompt_NO_rate(1800, phi) for phi in phi_range]

ax3.plot(phi_range, prompt_vs_phi, 'b-', linewidth=2.5)
ax3.set_xlabel('Equivalence Ratio phi', fontsize=11)
ax3.set_ylabel('Prompt NO Rate (mol/cm^3/s)', fontsize=11)
ax3.set_title('Prompt NO vs Equivalence Ratio', fontsize=12, fontweight='bold')
ax3.axvline(x=1.0, color='r', linestyle='--', linewidth=1.5, label='Stoichiometric')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 4. 燃料氮含量对燃料型NO的影响
ax4 = axes[1, 1]
N_content = np.linspace(0, 0.02, 50)  # 0-2%
fuel_vs_N = [fuel_NO_rate(1600, N) for N in N_content]

ax4.plot(N_content*100, fuel_vs_N, 'g-', linewidth=2.5)
ax4.set_xlabel('Fuel Nitrogen Content (%)', fontsize=11)
ax4.set_ylabel('Fuel NO Rate (mol/cm^3/s)', fontsize=11)
ax4.set_title('Fuel NO vs Fuel Nitrogen Content', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

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

print("✓ 实验2完成：NOx生成路径分析")

# ============================================================
# 实验3：CO生成与氧化
# ============================================================
print("\n实验3：CO生成与氧化")
print("-" * 40)

# CO氧化机理
# CO + OH -> CO2 + H
k_CO_OH = lambda T: 2.76e4 * T**1.5 * np.exp(116/T)  # cm^3/mol/s

# 模拟CO氧化过程
def CO_oxidation(y, t, T, phi):
    CO, OH = y
    
    # 反应速率
    r_ox = k_CO_OH(T) * CO * OH
    
    # CO消耗
    dCOdt = -r_ox
    
    # OH再生（简化）
    dOHdt = r_ox - 0.1*OH  # 假设OH有一定消耗
    
    return [dCOdt, dOHdt]

# 不同温度下的CO氧化
T_CO = [1200, 1400, 1600, 1800, 2000]
t_CO = np.linspace(0, 0.01, 1000)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. CO氧化曲线
ax1 = axes[0, 0]
colors_CO = plt.cm.plasma(np.linspace(0, 1, len(T_CO)))

for i, T in enumerate(T_CO):
    y0 = [1e-6, 1e-8]  # 初始CO和OH浓度 (mol/cm^3)
    solution = odeint(CO_oxidation, y0, t_CO, args=(T, 1.0))
    CO_conc = solution[:, 0]
    
    # 归一化
    CO_norm = CO_conc / CO_conc[0]
    
    ax1.semilogy(t_CO*1000, CO_norm, linewidth=2.5, color=colors_CO[i], label=f'T={T}K')

ax1.set_xlabel('Time (ms)', fontsize=11)
ax1.set_ylabel('Normalized CO Concentration', fontsize=11)
ax1.set_title('CO Oxidation at Different Temperatures', fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 2. CO氧化速率常数
ax2 = axes[0, 1]
T_k = np.linspace(1000, 2200, 100)
k_values = [k_CO_OH(T) for T in T_k]

ax2.semilogy(1000/T_k, k_values, 'b-', linewidth=2.5)
ax2.set_xlabel('1000/T (K^-1)', fontsize=11)
ax2.set_ylabel('Rate Constant (cm^3/mol/s)', fontsize=11)
ax2.set_title('CO + OH Rate Constant', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.invert_xaxis()

# 3. 当量比对CO排放的影响
ax3 = axes[1, 0]
phi_CO = np.linspace(0.6, 1.4, 50)

# 简化模型：CO排放与当量比的关系
CO_emission = []
for phi in phi_CO:
    if phi < 0.9:
        # 贫燃区CO低
        CO_emission.append(100 * phi)
    elif phi < 1.05:
        # 接近化学计量比
        CO_emission.append(500 * abs(phi - 1.0) + 100)
    else:
        # 富燃区CO高
        CO_emission.append(2000 * (phi - 1.0)**2 + 500)

ax3.plot(phi_CO, CO_emission, 'r-', linewidth=2.5)
ax3.set_xlabel('Equivalence Ratio phi', fontsize=11)
ax3.set_ylabel('CO Emission (ppm)', fontsize=11)
ax3.set_title('CO Emission vs Equivalence Ratio', fontsize=12, fontweight='bold')
ax3.axvline(x=1.0, color='k', linestyle='--', linewidth=1.5, label='Stoichiometric')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 4. CO-NOx权衡曲线
ax4 = axes[1, 1]
# 模拟不同当量比下的CO和NOx排放
phi_trade = np.linspace(0.7, 1.3, 50)

NOx_trade = []
CO_trade = []

for phi in phi_trade:
    # NOx：贫燃高，富燃低
    T_ad = 2200 * np.exp(-(phi-0.95)**2/0.1)
    NOx = 1000 * (T_ad/2000)**4 * np.exp(-phi) if phi < 1.2 else 100
    
    # CO：富燃高，贫燃低
    if phi < 0.9:
        CO = 50 * phi
    elif phi < 1.05:
        CO = 200 + 1000 * abs(phi - 1.0)
    else:
        CO = 500 + 3000 * (phi - 1.0)**2
    
    NOx_trade.append(NOx)
    CO_trade.append(CO)

ax4.plot(phi_trade, NOx_trade, 'b-', linewidth=2.5, label='NOx')
ax4_twin = ax4.twinx()
ax4_twin.plot(phi_trade, CO_trade, 'r-', linewidth=2.5, label='CO')

ax4.set_xlabel('Equivalence Ratio phi', fontsize=11)
ax4.set_ylabel('NOx (ppm)', fontsize=11, color='b')
ax4_twin.set_ylabel('CO (ppm)', fontsize=11, color='r')
ax4.set_title('CO-NOx Trade-off', fontsize=12, fontweight='bold')
ax4.tick_params(axis='y', labelcolor='b')
ax4_twin.tick_params(axis='y', labelcolor='r')
ax4.grid(True, alpha=0.3)

# 合并图例
lines1, labels1 = ax4.get_legend_handles_labels()
lines2, labels2 = ax4_twin.get_legend_handles_labels()
ax4.legend(lines1 + lines2, labels1 + labels2, loc='upper left', fontsize=10)

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

print("✓ 实验3完成：CO生成与氧化")

# ============================================================
# 实验4：碳烟生成模拟
# ============================================================
print("\n实验4：碳烟生成模拟")
print("-" * 40)

# Hiroyasu碳烟模型参数
A_sf = 100.0  # 表面生长系数
A_sc = 10.0   # 氧化系数
Ea_sf = 8000  # 表面生长活化能 (K)
Ea_sc = 12000 # 氧化活化能 (K)

# 碳烟生成ODE
def soot_formation(y, t, T, phi, P_atm=1.0):
    m_soot = y[0]
    
    # 燃料质量（简化假设随时间减少）
    m_fuel = 1e-3 * np.exp(-10*t)
    
    # 表面生长项
    growth = A_sf * np.sqrt(m_fuel) * P_atm**1.8 * np.exp(-Ea_sf/T)
    
    # 氧化项（与O2浓度相关）
    if phi < 1.0:
        Y_O2 = 0.23 * (1 - phi)
    else:
        Y_O2 = 0.01  # 富燃区O2很低
    
    oxidation = A_sc * m_soot * (Y_O2/0.23)**1.8 * np.exp(-Ea_sc/T)
    
    dmdt = growth - oxidation
    return [dmdt]

# 模拟不同条件下的碳烟生成
conditions = [
    (1800, 1.0, 'T=1800K, phi=1.0'),
    (2000, 1.0, 'T=2000K, phi=1.0'),
    (1800, 1.2, 'T=1800K, phi=1.2'),
    (1800, 0.9, 'T=1800K, phi=0.9'),
]

t_soot = np.linspace(0, 0.05, 500)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. 碳烟质量演化
ax1 = axes[0, 0]
colors_soot = plt.cm.coolwarm(np.linspace(0, 1, len(conditions)))

for i, (T, phi, label) in enumerate(conditions):
    y0 = [1e-9]  # 初始碳烟质量
    solution = odeint(soot_formation, y0, t_soot, args=(T, phi))
    m_soot = solution[:, 0]
    
    ax1.plot(t_soot*1000, m_soot*1e6, linewidth=2.5, color=colors_soot[i], label=label)

ax1.set_xlabel('Time (ms)', fontsize=11)
ax1.set_ylabel('Soot Mass (mg)', fontsize=11)
ax1.set_title('Soot Formation and Oxidation', fontsize=12, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 2. 温度对碳烟的影响
ax2 = axes[0, 1]
T_soot = np.linspace(1400, 2200, 50)
soot_vs_T = []

for T in T_soot:
    y0 = [1e-9]
    solution = odeint(soot_formation, y0, t_soot, args=(T, 1.2))
    soot_vs_T.append(solution[-1, 0] * 1e6)

ax2.plot(T_soot, soot_vs_T, 'k-', linewidth=2.5)
ax2.set_xlabel('Temperature (K)', fontsize=11)
ax2.set_ylabel('Final Soot Mass (mg)', fontsize=11)
ax2.set_title('Soot vs Temperature (phi=1.2)', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 3. 当量比对碳烟的影响
ax3 = axes[1, 0]
phi_soot = np.linspace(0.7, 1.5, 50)
soot_vs_phi = []

for phi in phi_soot:
    y0 = [1e-9]
    solution = odeint(soot_formation, y0, t_soot, args=(1800, phi))
    soot_vs_phi.append(solution[-1, 0] * 1e6)

ax3.plot(phi_soot, soot_vs_phi, 'b-', linewidth=2.5)
ax3.set_xlabel('Equivalence Ratio phi', fontsize=11)
ax3.set_ylabel('Final Soot Mass (mg)', fontsize=11)
ax3.set_title('Soot vs Equivalence Ratio (T=1800K)', fontsize=12, fontweight='bold')
ax3.axvline(x=1.0, color='r', linestyle='--', linewidth=1.5, label='Stoichiometric')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 4. 碳烟生成-氧化平衡
ax4 = axes[1, 1]
# 表面生长和氧化速率对比
T_balance = np.linspace(1400, 2200, 50)

# 固定条件
m_fuel_fix = 1e-4
m_soot_fix = 1e-6
phi_fix = 1.2
P_fix = 1.0

growth_rate = []
oxid_rate = []

for T in T_balance:
    gr = A_sf * np.sqrt(m_fuel_fix) * P_fix**1.8 * np.exp(-Ea_sf/T)
    Y_O2 = 0.01 if phi_fix > 1.0 else 0.23 * (1 - phi_fix)
    ox = A_sc * m_soot_fix * (Y_O2/0.23)**1.8 * np.exp(-Ea_sc/T)
    growth_rate.append(gr)
    oxid_rate.append(ox)

ax4.semilogy(1000/T_balance, growth_rate, 'g-', linewidth=2.5, label='Growth')
ax4.semilogy(1000/T_balance, oxid_rate, 'r-', linewidth=2.5, label='Oxidation')
ax4.set_xlabel('1000/T (K^-1)', fontsize=11)
ax4.set_ylabel('Rate (mg/s)', fontsize=11)
ax4.set_title('Soot Growth vs Oxidation Rates', fontsize=12, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
ax4.invert_xaxis()

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

print("✓ 实验4完成：碳烟生成模拟")

# ============================================================
# 实验5：污染物控制策略
# ============================================================
print("\n实验5：污染物控制策略")
print("-" * 40)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. EGR对NOx的影响
ax1 = axes[0, 0]
EGR_rate = np.linspace(0, 0.4, 50)  # 0-40% EGR

# 简化模型：EGR降低峰值温度
NOx_EGR = []
for egr in EGR_rate:
    T_peak = 2200 * (1 - 0.3*egr)  # EGR降低温度
    NOx = 2000 * (T_peak/2200)**4 * np.exp(-egr)  # EGR也稀释反应物
    NOx_EGR.append(NOx)

NOx_baseline = NOx_EGR[0]
NOx_reduction = [(1 - n/NOx_baseline)*100 for n in NOx_EGR]

ax1.plot(EGR_rate*100, NOx_reduction, 'b-', linewidth=2.5)
ax1.set_xlabel('EGR Rate (%)', fontsize=11)
ax1.set_ylabel('NOx Reduction (%)', fontsize=11)
ax1.set_title('EGR Effect on NOx Reduction', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 2. 分级燃烧效果
ax2 = axes[0, 1]
# 主燃区当量比
phi_primary = np.linspace(0.8, 1.2, 50)

NOx_staged = []
for phi_p in phi_primary:
    # 贫燃主燃区降低NOx
    if phi_p < 1.0:
        NOx = 500 * phi_p**2
    else:
        NOx = 500 + 2000*(phi_p - 1.0)**2
    NOx_staged.append(NOx)

# 对比非分级燃烧
NOx_unstaged = [1500] * len(phi_primary)

ax2.plot(phi_primary, NOx_unstaged, 'r--', linewidth=2.5, label='Unstaged')
ax2.plot(phi_primary, NOx_staged, 'b-', linewidth=2.5, label='Staged')
ax2.set_xlabel('Primary Zone Equivalence Ratio', fontsize=11)
ax2.set_ylabel('NOx Emission (ppm)', fontsize=11)
ax2.set_title('Staged Combustion Effect', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# 3. SNCR效率
ax3 = axes[1, 0]
T_SNCR = np.linspace(700, 1300, 100)  # °C

# SNCR效率曲线（温度窗口）
SNCR_eff = []
for T in T_SNCR:
    if 850 <= T <= 1100:
        # 最佳温度窗口
        eff = 40 + 30 * np.exp(-((T-950)/150)**2)
    elif T < 850:
        # 温度太低
        eff = 40 * np.exp(-((T-850)/100)**2)
    else:
        # 温度太高（NH3氧化）
        eff = 70 * np.exp(-((T-1100)/100)**2)
    SNCR_eff.append(eff)

ax3.plot(T_SNCR, SNCR_eff, 'g-', linewidth=2.5)
ax3.axvline(x=850, color='r', linestyle='--', linewidth=1.5, alpha=0.7)
ax3.axvline(x=1100, color='r', linestyle='--', linewidth=1.5, alpha=0.7)
ax3.fill_between([850, 1100], [0, 0], [100, 100], alpha=0.2, color='green')
ax3.set_xlabel('Temperature (°C)', fontsize=11)
ax3.set_ylabel('DeNOx Efficiency (%)', fontsize=11)
ax3.set_title('SNCR Temperature Window', fontsize=12, fontweight='bold')
ax3.set_ylim([0, 100])
ax3.grid(True, alpha=0.3)

# 4. 综合控制策略对比
ax4 = axes[1, 1]
strategies = ['Baseline', 'EGR 20%', 'Staged', 'SNCR', 'Combined']
NOx_values = [1500, 800, 500, 300, 100]
CO_values = [200, 250, 300, 200, 220]

x = np.arange(len(strategies))
width = 0.35

bars1 = ax4.bar(x - width/2, NOx_values, width, label='NOx', color='red', alpha=0.7)
bars2 = ax4.bar(x + width/2, CO_values, width, label='CO', color='blue', alpha=0.7)

ax4.set_ylabel('Emission (ppm)', fontsize=11)
ax4.set_title('Emission Control Strategies Comparison', fontsize=12, fontweight='bold')
ax4.set_xticks(x)
ax4.set_xticklabels(strategies, fontsize=10)
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3, axis='y')

# 添加数值标签
for bar in bars1:
    height = bar.get_height()
    ax4.text(bar.get_x() + bar.get_width()/2., height,
             f'{int(height)}', ha='center', va='bottom', fontsize=9)

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

print("✓ 实验5完成：污染物控制策略")

# ============================================================
# 实验6：综合排放预测
# ============================================================
print("\n实验6：综合排放预测")
print("-" * 40)

# 建立简化的排放预测模型
def emission_model(T, phi, residence_time, EGR=0):
    """
    简化的排放预测模型
    输入: T (K), phi (当量比), residence_time (s), EGR (0-1)
    输出: NOx (ppm), CO (ppm), Soot (mg/m^3)
    """
    # 调整温度（EGR效应）
    T_eff = T * (1 - 0.3*EGR)
    
    # NOx预测（热力型主导）
    if phi < 1.2:
        NOx = 100 * (T_eff/1800)**4 * np.exp(-phi) * (1 - EGR)
    else:
        NOx = 50 * (T_eff/1800)**2  # 富燃区NOx降低
    
    # CO预测
    if phi < 0.9:
        CO = 50 + 100*(0.9-phi)
    elif phi > 1.05:
        CO = 200 + 1000*(phi-1.0)**2
    else:
        CO = 50
    
    # 碳烟预测
    if phi > 1.0:
        Soot = 10 * (phi-1.0)**2 * (T_eff/1800)**2 * 1000  # mg/m^3
    else:
        Soot = 1
    
    # 停留时间效应（不完全燃烧）
    if residence_time < 0.001:
        CO *= 2
        Soot *= 1.5
    
    return NOx, CO, Soot

# 参数空间分析
T_range = np.linspace(1600, 2200, 30)
phi_range = np.linspace(0.7, 1.4, 30)

T_grid, phi_grid = np.meshgrid(T_range, phi_range)
NOx_grid = np.zeros_like(T_grid)
CO_grid = np.zeros_like(T_grid)
Soot_grid = np.zeros_like(T_grid)

for i in range(len(phi_range)):
    for j in range(len(T_range)):
        NOx, CO, Soot = emission_model(T_grid[i,j], phi_grid[i,j], 0.01)
        NOx_grid[i,j] = NOx
        CO_grid[i,j] = CO
        Soot_grid[i,j] = Soot

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. NOx分布
ax1 = axes[0, 0]
contour1 = ax1.contourf(T_grid, phi_grid, NOx_grid, levels=20, cmap='Reds')
ax1.set_xlabel('Temperature (K)', fontsize=11)
ax1.set_ylabel('Equivalence Ratio phi', fontsize=11)
ax1.set_title('NOx Emission Map', fontsize=12, fontweight='bold')
cbar1 = plt.colorbar(contour1, ax=ax1)
cbar1.set_label('NOx (ppm)', fontsize=10)
ax1.axhline(y=1.0, color='k', linestyle='--', linewidth=1)

# 2. CO分布
ax2 = axes[0, 1]
contour2 = ax2.contourf(T_grid, phi_grid, CO_grid, levels=20, cmap='Blues')
ax2.set_xlabel('Temperature (K)', fontsize=11)
ax2.set_ylabel('Equivalence Ratio phi', fontsize=11)
ax2.set_title('CO Emission Map', fontsize=12, fontweight='bold')
cbar2 = plt.colorbar(contour2, ax=ax2)
cbar2.set_label('CO (ppm)', fontsize=10)
ax2.axhline(y=1.0, color='k', linestyle='--', linewidth=1)

# 3. 碳烟分布
ax3 = axes[1, 0]
contour3 = ax3.contourf(T_grid, phi_grid, Soot_grid, levels=20, cmap='Greys')
ax3.set_xlabel('Temperature (K)', fontsize=11)
ax3.set_ylabel('Equivalence Ratio phi', fontsize=11)
ax3.set_title('Soot Emission Map', fontsize=12, fontweight='bold')
cbar3 = plt.colorbar(contour3, ax=ax3)
cbar3.set_label('Soot (mg/m^3)', fontsize=10)
ax3.axhline(y=1.0, color='r', linestyle='--', linewidth=1)

# 4. 综合排放指数
ax4 = axes[1, 1]
# 综合排放指数（加权）
Emission_Index = NOx_grid/1000 + CO_grid/500 + Soot_grid/50
contour4 = ax4.contourf(T_grid, phi_grid, Emission_Index, levels=20, cmap='viridis')
ax4.set_xlabel('Temperature (K)', fontsize=11)
ax4.set_ylabel('Equivalence Ratio phi', fontsize=11)
ax4.set_title('Combined Emission Index', fontsize=12, fontweight='bold')
cbar4 = plt.colorbar(contour4, ax=ax4)
cbar4.set_label('Emission Index', fontsize=10)
ax4.axhline(y=1.0, color='r', linestyle='--', linewidth=1)

# 标记最佳区域（低排放区）
min_idx = np.unravel_index(np.argmin(Emission_Index), Emission_Index.shape)
ax4.plot(T_grid[min_idx], phi_grid[min_idx], 'r*', markersize=15, 
         label=f'Optimal: T={T_grid[min_idx]:.0f}K, phi={phi_grid[min_idx]:.2f}')
ax4.legend(fontsize=9)

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

print("✓ 实验6完成：综合排放预测")
print(f"  最优工况: T={T_grid[min_idx]:.0f}K, phi={phi_grid[min_idx]:.2f}")

print("\n" + "="*60)
print("所有实验完成！")
print("="*60)