主题047：燃烧数值模拟中的不确定性量化

目录

1. 引言
2. 不确定性来源与分类
3. 概率统计基础
4. 敏感性分析方法
5. 不确定性传播方法
6. Python仿真实现
7. 工程应用案例
8. 总结与展望

---

1. 引言

1.1 不确定性量化的重要性

燃烧数值模拟在工程设计和科学研究中发挥着越来越重要的作用，但模拟结果往往存在不确定性：

• 模型不确定性：湍流模型、燃烧模型、辐射模型的简化假设
• 参数不确定性：边界条件、物性参数、化学反应速率
• 数值不确定性：网格离散、时间步长、收敛判据

**不确定性量化(UQ)**的目标：

1. 识别和量化不确定性的来源
2. 评估不确定性对结果的影响
3. 提高模拟结果的可靠性和可信度
4. 为决策提供风险评估依据

1.2 燃烧模拟中的典型不确定性

物理模型不确定性：

• 湍流-燃烧相互作用模型
• 详细反应机理的简化
• 辐射传热模型

输入参数不确定性：

• 燃料组分和热值波动
• 入口流速和温度测量误差
• 几何尺寸的制造公差

数值方法不确定性：

• 网格密度和类型
• 数值格式和格式参数
• 迭代收敛标准

1.3 不确定性量化方法概述

主要方法分类：

1. 敏感性分析(SA)：

   • 局部敏感性分析
   • 全局敏感性分析
   • 筛选方法

2. 不确定性传播(UP)：

   • 蒙特卡洛方法
   • 多项式混沌展开
   • 随机配点法

3. 模型验证与确认(V&V)：

   • 代码验证
   • 模型确认
   • 不确定性量化

---

2. 不确定性来源与分类

2.1 认知不确定性 vs 偶然不确定性

2.1.1 认知不确定性(Epistemic)

定义：由于缺乏知识或信息不足导致的不确定性

特点：

• 可以通过更多研究减少
• 通常用区间或模糊数表示
• 主观性较强

燃烧模拟中的例子：

• 反应机理的缺失路径
• 湍流模型的适用性
• 边界条件的精确值

2.1.2 偶然不确定性(Aleatory)

定义：由于系统固有的随机性导致的不确定性

特点：

• 不可通过研究消除
• 通常用概率分布描述
• 客观存在

燃烧模拟中的例子：

• 燃料组分的自然波动
• 环境温度的随机变化
• 测量仪器的随机误差

2.2 不确定性分类框架

不确定性
├── 模型不确定性
│   ├── 模型形式不确定性
│   ├── 模型参数不确定性
│   └── 模型适用性不确定性
├── 输入不确定性
│   ├── 边界条件不确定性
│   ├── 初始条件不确定性
│   └── 物性参数不确定性
└── 数值不确定性
    ├── 离散化不确定性
    ├── 迭代不确定性
    └── 舍入不确定性

2.3 燃烧模拟中的关键不确定参数

2.3.1 化学反应动力学参数

反应速率常数：

$$k=A\cdot T^n\cdot \exp \left(-\frac{E_a}{RT}\right)$$

不确定性来源：

• 指前因子A：±20-50%
• 活化能Ea：±5-10%
• 温度指数n：±0.1-0.3

示例：甲烷燃烧关键反应

反应	A (不确定性)	Ea (不确定性)
CH₄ + OH → CH₃ + H₂O	±30%	±5%
CH₃ + O₂ → CH₂O + OH	±50%	±10%
CO + OH → CO₂ + H	±20%	±3%

2.3.2 湍流模型参数

k-ε模型参数：

参数	标准值	典型不确定性
C_μ	0.09	±10%
C_ε1	1.44	±5%
C_ε2	1.92	±5%
σ_k	1.0	±10%
σ_ε	1.3	±10%

影响：

• 射流扩散角度变化：±15%
• 回流区长度变化：±20%
• 混合效率变化：±10%

2.3.3 边界条件不确定性

典型边界条件不确定性：

参数	测量精度	不确定性范围
入口流速	±1% FS	±0.5-2%
入口温度	±1°C	±5-20K
燃料流量	±0.5%	±1-3%
压力	±0.1% FS	±100-500 Pa

---

3. 概率统计基础

3.1 概率分布类型

3.1.1 常用连续分布

正态分布(Gaussian)：

$$f(x)=\frac{1}{\sigma \sqrt{2\pi }}\exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)$$

适用场景：

• 测量误差
• 制造公差
• 自然变异

对数正态分布：

$$f(x)=\frac{1}{x\sigma \sqrt{2\pi }}\exp \left(-\frac{(\ln x-\mu {)}^2}{2{\sigma }^2}\right)$$

适用场景：

• 反应速率常数
• 颗粒尺寸分布
• 寿命数据

均匀分布：

$$f(x)=\{\begin{array}{ll}\frac{1}{b-a}& a\le x\le b\\ 0& \text{其他}\end{array}$$

适用场景：

• 缺乏数据时的保守估计
• 模型参数的先验分布
• 舍入误差

三角分布：

$$f(x)=\{\begin{array}{ll}\frac{2(x-a)}{(b-a)(c-a)}& a\le x\le c\\ \frac{2(b-x)}{(b-a)(b-c)}& c<x\le b\\ 0& \text{其他}\end{array}$$

适用场景：

• 专家判断
• 有限数据
• 保守估计

3.1.2 分布选择指南

数据特征	推荐分布	理由
对称、无界	正态分布	中心极限定理
正、右偏	对数正态	物理约束
有明确范围	均匀/三角	最大熵原理
极端值重要	威布尔/极值	尾部特性

3.2 统计量与置信区间

3.2.1 描述性统计量

集中趋势：

• 均值(Mean)：$\mu =\frac{1}{N}{\sum }_{i=1}^N{x}_i$
• 中位数(Median)：排序后中间值
• 众数(Mode)：出现频率最高值

离散程度：

• 标准差：$\sigma =\sqrt{\frac{1}{N}{\sum }_{i=1}^N(x_i-\mu {)}^2}$
• 变异系数：$CV=\frac{\sigma }{\mu }\times 100\%$
• 四分位距：$IQR=Q_3-Q_1$

形状特征：

• 偏度：${\gamma }_1=\frac{1}{N}{\sum }_{i=1}^N{\left(\frac{x_i-\mu }{\sigma }\right)}^3$
• 峰度：${\gamma }_2=\frac{1}{N}{\sum }_{i=1}^N{\left(\frac{x_i-\mu }{\sigma }\right)}^4-3$

3.2.2 置信区间

正态分布的置信区间：

$$CI=\bar{x}\pm z_{\alpha /2}\cdot \frac{s}{\sqrt{n}}$$

常用置信水平：

置信水平	z值
90%	1.645
95%	1.96
99%	2.576

Bootstrap置信区间：

1. 从原始样本中有放回地抽取n个样本
2. 计算统计量（如均值）
3. 重复B次（通常B=1000-10000）
4. 取统计量的分位数作为置信区间

3.3 假设检验

3.3.1 参数检验

t检验（均值比较）：

$$t=\frac{{\bar{x}}_1-{\bar{x}}_2}{s_p\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}}$$

F检验（方差比较）：

$$F=\frac{s_1^2}{s_2^2}$$

适用条件：

• 数据近似正态分布
• 样本独立
• 方差齐性（t检验）

3.3.2 非参数检验

Kolmogorov-Smirnov检验：

$$D=\underset{x}{\max }|F_n(x)-F(x)|$$

Mann-Whitney U检验：

• 不假设分布形式
• 比较两个独立样本
• 基于秩次

---

4. 敏感性分析方法

4.1 局部敏感性分析

4.1.1 单因素扰动法

基本原理：

$$S_i=\frac{\partial y}{\partial x_i}\approx \frac{y(x_i+\Delta x_i)-y(x_i)}{\Delta x_i}$$

无量纲化敏感性系数：

$$S_i^{\ast }=\frac{\partial y/y}{\partial x_i/x_i}=\frac{x_i}{y}\cdot \frac{\partial y}{\partial x_i}$$

实施步骤：

1. 选择基准点（通常是 nominal 值）
2. 对每个参数进行小扰动（±1-5%）
3. 计算输出变化
4. 计算敏感性系数

优缺点：

• 优点：计算量小，易于实现
• 缺点：仅适用于线性或弱非线性系统，忽略参数间交互作用

4.1.2 有限差分法

前向差分：

$$\frac{\partial y}{\partial x_i}\approx \frac{y(x_i+h)-y(x_i)}{h}$$

中心差分：

$$\frac{\partial y}{\partial x_i}\approx \frac{y(x_i+h)-y(x_i-h)}{2h}$$

步长选择：

• 太大：截断误差增大
• 太小：舍入误差增大
• 经验法则：$h\approx \sqrt{{ϵ}_{mach}}\cdot x_i$

4.2 全局敏感性分析

4.2.1 Morris筛选法

基本思想：

• 在参数空间内随机采样
• 计算基本效应(Elementary Effects)
• 通过统计量识别重要参数

基本效应计算：

$$E{E}_i=\frac{y(x_1,...,x_i+\Delta ,...,x_k)-y(x_1,...,x_i,...,x_k)}{\Delta }$$

统计量：

• 均值$\mu$：参数的总体影响
• 标准差$\sigma$：参数的非线性程度和交互作用

参数分类：

类型	特征	处理策略
不重要	低$\mu$，低$\sigma$	固定为名义值
线性重要	高$\mu$，低$\sigma$	精确建模
非线性重要	高$\mu$，高$\sigma$	重点研究
交互作用	低$\mu$，高$\sigma$	考虑交互项

4.2.2 Sobol敏感性分析

方差分解：

$$V(Y)=\sum _i{V}_i+\sum _{i<j}{V}_{ij}+...+V_{\mathrm{12...}k}$$

一阶敏感性指数：

$$S_i=\frac{V_i}{V(Y)}=\frac{V_{X_i}(E_{X_{\sim i}}(Y|X_i))}{V(Y)}$$

总敏感性指数：

$$S_{Ti}=\frac{E_{X_{\sim i}}(V_{X_i}(Y|X_{\sim i}))}{V(Y)}=1-\frac{V_{X_{\sim i}}(E_{X_i}(Y|X_{\sim i}))}{V(Y)}$$

解释：

• $S_i$：参数$X_i$的单独贡献
• $S_{Ti}$：参数$X_i$的总贡献（包括交互作用）
• $S_{Ti}-S_i$：参数$X_i$与其他参数的交互作用

计算方法：

• 蒙特卡洛积分
• 准蒙特卡洛（Sobol序列、Halton序列）
• 多项式混沌展开

4.2.3 FAST方法

Fourier Amplitude Sensitivity Test

原理：

• 将参数映射为频率
• 通过傅里叶分析识别各参数的贡献
• 每个参数分配不同特征频率

优点：

• 计算效率高于Sobol方法
• 可以同时计算一阶和总敏感性指数

4.3 敏感性分析在燃烧模拟中的应用

4.3.1 反应机理简化

目标：识别对目标量影响最大的基元反应

步骤：

1. 定义目标量（如点火延迟时间、火焰速度）
2. 对所有反应速率常数进行敏感性分析
3. 识别高敏感性反应
4. 保留关键反应，简化次要反应

示例结果：

反应	敏感性系数	重要性
H + O₂ → OH + O	0.85	关键
CH₄ + OH → CH₃ + H₂O	0.62	重要
CH₃ + O → CH₂O + H	0.15	次要
…	<0.1	可简化

4.3.2 模型参数标定

目标：确定模型参数的最优值和不确定性范围

方法：

1. 选择待标定参数
2. 定义目标函数（模拟值与实验值的偏差）
3. 进行敏感性分析，识别关键参数
4. 使用优化算法调整关键参数
5. 量化参数不确定性

---

5. 不确定性传播方法

5.1 蒙特卡洛方法

5.1.1 基本蒙特卡洛

算法流程：

1. 定义输入参数的概率分布
2. For i = 1 to N:
   a. 从输入分布中随机采样
   b. 运行模拟计算
   c. 记录输出结果
3. 统计分析输出结果

样本量确定：

对于95%置信水平、5%相对误差：

$$N={\left(\frac{z_{\alpha /2}\cdot CV}{ϵ}\right)}^2\approx {\left(\frac{1.96\cdot CV}{0.05}\right)}^2$$

典型样本量：

• 低维问题(k<10)：10³-10⁴
• 中维问题(10<k<50)：10⁴-10⁵
• 高维问题(k>50)：10⁵-10⁶

5.1.2 方差缩减技术

拉丁超立方采样(LHS)：

• 将每个参数的范围分成N个等概率区间
• 在每个区间中随机采样一个点
• 组合成N个样本点

优点：

• 比纯随机采样更均匀地覆盖参数空间
• 可以用更少的样本达到相同的精度
• 通常可以减少10-50%的样本量

重要性采样：

$$E[g(X)]=\int g(x)f(x)dx=\int g(x)\frac{f(x)}{h(x)}h(x)dx$$

适用场景：

• 小概率事件估计
• 尾部特性分析
• 失效概率计算

5.2 谱方法

5.2.1 多项式混沌展开(PCE)

基本思想：
将随机输出表示为随机输入的正交多项式展开：

$$Y(\xi )=\sum _{\alpha }{c}_{\alpha }{\Psi }_{\alpha }(\xi )$$

其中：

• $\xi$：标准随机变量
• ${\Psi }_{\alpha }$：正交多项式基函数
• $c_{\alpha }$：展开系数

多项式选择：

输入分布	正交多项式	权重函数
正态分布	Hermite	$e^{-x^2/2}$
均匀分布	Legendre	1
伽马分布	Laguerre	$x^{\alpha }{e}^{-x}$
贝塔分布	Jacobi	$(1-x{)}^{\alpha }(1+x{)}^{\beta }$

系数计算：

投影法：

$$c_{\alpha }=\frac{⟨Y,{\Psi }_{\alpha }⟩}{⟨{\Psi }_{\alpha },{\Psi }_{\alpha }⟩}=\frac{E[Y{\Psi }_{\alpha }]}{E[{\Psi }_{\alpha }^2]}$$

回归法：

• 在配置点上计算输出
• 使用最小二乘法拟合系数

5.2.2 随机配点法

基本思想：

• 在精心选择的配点上求解确定性问题
• 使用插值或回归构造响应面
• 基于响应面进行统计分析

配点选择：

• 全因子设计（低维）
• 稀疏网格（中维）
• 自适应配点（高维）

优点：

• 可以使用现有的确定性求解器
• 非侵入式方法
• 适用于复杂模型

5.3 矩方法

5.3.1 泰勒展开法

一阶近似：

$$Y\approx Y({\mu }_X)+\sum _{i=1}^k\frac{\partial Y}{\partial X_i}(X_i-{\mu }_{X_i})$$

均值：

$${\mu }_Y\approx Y({\mu }_X)$$

方差：

$${\sigma }_Y^2\approx \sum _{i=1}^k{\left(\frac{\partial Y}{\partial X_i}\right)}^2{\sigma }_{X_i}^2$$

适用条件：

• 输入不确定性小
• 系统接近线性
• 快速估计

5.3.2 点估计法

Rosenblueth方法：

对于k个参数，需要计算2^k个点：

$$Y_{\pm }=Y({\mu }_{X_1}\pm {\sigma }_{X_1},...,{\mu }_{X_k}\pm {\sigma }_{X_k})$$

均值和方差：

$${\mu }_Y=\frac{1}{2^k}\sum Y_{\pm }$$

$${\sigma }_Y^2=\frac{1}{2^k}\sum (Y_{\pm }-{\mu }_Y{)}^2$$

改进方法：

• Harr方法
• 改进Rosenblueth方法
• 仅适用于低维问题

---

6. Python仿真实现

6.1 概率分布与采样仿真

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy import stats

print("=" * 70)
print("主题047：不确定性量化方法 - Python仿真")
print("=" * 70)

# ============================================================================
# 实验1: 概率分布与采样方法
# ============================================================================
print("\n" + "=" * 70)
print("实验1: 概率分布与采样方法")
print("=" * 70)

# 设置随机种子以保证可重复性
np.random.seed(42)

# 定义燃烧模拟中常见的概率分布
print("\n常见概率分布定义:")
print("-" * 50)

# 1. 正态分布 - 测量误差
temp_mean = 300  # K
temp_std = 10    # K
temp_dist = stats.norm(temp_mean, temp_std)
print(f"入口温度: 正态分布 N({temp_mean}, {temp_std}²)")

# 2. 对数正态分布 - 反应速率常数
A_mean = 1e10
A_sigma = 0.5  # 对数标准差
A_dist = stats.lognorm(s=A_sigma, scale=A_mean)
print(f"指前因子: 对数正态分布 LN({A_mean}, {A_sigma}²)")

# 3. 均匀分布 - 缺乏精确数据时
phi_low = 0.8
phi_high = 1.2
phi_dist = stats.uniform(phi_low, phi_high - phi_low)
print(f"当量比: 均匀分布 U({phi_low}, {phi_high})")

# 4. 三角分布 - 专家判断
pressure_min = 100000  # Pa
pressure_mode = 101325  # Pa
pressure_max = 102000  # Pa
pressure_dist = stats.triang(
    c=(pressure_mode - pressure_min)/(pressure_max - pressure_min),
    loc=pressure_min,
    scale=pressure_max - pressure_min
)
print(f"环境压力: 三角分布 T({pressure_min}, {pressure_mode}, {pressure_max})")

# 生成样本
n_samples = 10000
temp_samples = temp_dist.rvs(n_samples)
A_samples = A_dist.rvs(n_samples)
phi_samples = phi_dist.rvs(n_samples)
pressure_samples = pressure_dist.rvs(n_samples)

print(f"\n生成样本统计 (n={n_samples}):")
print("-" * 50)
print(f"温度样本: 均值={np.mean(temp_samples):.2f}, 标准差={np.std(temp_samples):.2f}")
print(f"A因子样本: 均值={np.mean(A_samples):.2e}, 标准差={np.std(A_samples):.2e}")
print(f"当量比样本: 均值={np.mean(phi_samples):.3f}, 标准差={np.std(phi_samples):.3f}")
print(f"压力样本: 均值={np.mean(pressure_samples):.0f}, 标准差={np.std(pressure_samples):.0f}")

# 绘制分布
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 温度分布
axes[0, 0].hist(temp_samples, bins=50, density=True, alpha=0.7, color='blue', edgecolor='black')
x_temp = np.linspace(temp_mean - 4*temp_std, temp_mean + 4*temp_std, 100)
axes[0, 0].plot(x_temp, temp_dist.pdf(x_temp), 'r-', linewidth=2, label='PDF')
axes[0, 0].axvline(temp_mean, color='g', linestyle='--', label=f'Mean={temp_mean}')
axes[0, 0].set_xlabel('Temperature (K)', fontsize=11)
axes[0, 0].set_ylabel('Probability Density', fontsize=11)
axes[0, 0].set_title('Normal Distribution: Inlet Temperature', fontsize=12)
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 对数正态分布
axes[0, 1].hist(A_samples, bins=50, density=True, alpha=0.7, color='green', edgecolor='black')
x_A = np.linspace(A_mean*0.1, A_mean*3, 100)
axes[0, 1].plot(x_A, A_dist.pdf(x_A), 'r-', linewidth=2, label='PDF')
axes[0, 1].axvline(A_mean, color='g', linestyle='--', label=f'Mean={A_mean:.0e}')
axes[0, 1].set_xlabel('Pre-exponential Factor', fontsize=11)
axes[0, 1].set_ylabel('Probability Density', fontsize=11)
axes[0, 1].set_title('Log-normal Distribution: Rate Constant', fontsize=12)
axes[0, 1].set_xscale('log')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# 均匀分布
axes[1, 0].hist(phi_samples, bins=50, density=True, alpha=0.7, color='orange', edgecolor='black')
x_phi = np.linspace(phi_low-0.1, phi_high+0.1, 100)
axes[1, 0].plot(x_phi, phi_dist.pdf(x_phi), 'r-', linewidth=2, label='PDF')
axes[1, 0].axvline((phi_low+phi_high)/2, color='g', linestyle='--', label=f'Mean={(phi_low+phi_high)/2}')
axes[1, 0].set_xlabel('Equivalence Ratio', fontsize=11)
axes[1, 0].set_ylabel('Probability Density', fontsize=11)
axes[1, 0].set_title('Uniform Distribution: Equivalence Ratio', fontsize=12)
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# 三角分布
axes[1, 1].hist(pressure_samples, bins=50, density=True, alpha=0.7, color='purple', edgecolor='black')
x_pressure = np.linspace(pressure_min-500, pressure_max+500, 100)
axes[1, 1].plot(x_pressure, pressure_dist.pdf(x_pressure), 'r-', linewidth=2, label='PDF')
axes[1, 1].axvline(pressure_mode, color='g', linestyle='--', label=f'Mode={pressure_mode}')
axes[1, 1].set_xlabel('Pressure (Pa)', fontsize=11)
axes[1, 1].set_ylabel('Probability Density', fontsize=11)
axes[1, 1].set_title('Triangular Distribution: Pressure', fontsize=12)
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('exp1_probability_distributions.png', dpi=150, bbox_inches='tight')
print("\n结果已保存至: exp1_probability_distributions.png")

plt.close()

6.2 敏感性分析仿真

# ============================================================================
# 实验2: 敏感性分析
# ============================================================================
print("\n" + "=" * 70)
print("实验2: 敏感性分析")
print("=" * 70)

def ignition_delay(T, P, phi, Ea=50000, A=1e10):
    """
    简化的点火延迟时间模型
    
    参数:
    T: 温度 (K)
    P: 压力 (Pa)
    phi: 当量比
    Ea: 活化能 (J/mol)
    A: 指前因子 (1/s)
    
    返回:
    点火延迟时间 (s)
    """
    R = 8.314  # J/(mol·K)
    # 简化Arrhenius模型
    tau = (1 / A) * np.exp(Ea / (R * T)) * (1e5 / P)**0.5 * (1 + abs(phi - 1))
    return tau

def flame_speed(T, P, phi, alpha=1.5, beta=0.3):
    """
    简化的层流火焰速度模型
    
    参数:
    T: 温度 (K)
    P: 压力 (Pa)
    phi: 当量比
    alpha, beta: 模型参数
    
    返回:
    火焰速度 (m/s)
    """
    # 简化模型
    S_l = 0.4 * (T / 300)**alpha * (P / 1e5)**(-beta) * np.exp(-2*(phi - 1.05)**2)
    return S_l

# 局部敏感性分析
print("\n局部敏感性分析:")
print("-" * 50)

# 基准值
T0, P0, phi0 = 300, 101325, 1.0
h = 0.01  # 扰动幅度

# 点火延迟时间的敏感性
tau0 = ignition_delay(T0, P0, phi0)

# 温度敏感性
tau_T_plus = ignition_delay(T0*(1+h), P0, phi0)
tau_T_minus = ignition_delay(T0*(1-h), P0, phi0)
S_T = (tau_T_plus - tau_T_minus) / (2*h*T0) * T0 / tau0

# 压力敏感性
tau_P_plus = ignition_delay(T0, P0*(1+h), phi0)
tau_P_minus = ignition_delay(T0, P0*(1-h), phi0)
S_P = (tau_P_plus - tau_P_minus) / (2*h*P0) * P0 / tau0

# 当量比敏感性
tau_phi_plus = ignition_delay(T0, P0, phi0*(1+h))
tau_phi_minus = ignition_delay(T0, P0, phi0*(1-h))
S_phi = (tau_phi_plus - tau_phi_minus) / (2*h*phi0) * phi0 / tau0

print(f"基准点火延迟时间: {tau0:.6f} s")
print(f"温度敏感性系数: {S_T:.3f}")
print(f"压力敏感性系数: {S_P:.3f}")
print(f"当量比敏感性系数: {S_phi:.3f}")

# 全局敏感性分析 - Morris方法简化版
print("\n全局敏感性分析 (Morris方法):")
print("-" * 50)

n_trajectories = 10
n_params = 3
delta = 0.1

param_names = ['Temperature', 'Pressure', 'Equivalence Ratio']
param_ranges = [(250, 350), (80000, 120000), (0.6, 1.4)]

EE_matrix = np.zeros((n_trajectories, n_params))

for t in range(n_trajectories):
    # 随机起点
    x_base = np.array([
        np.random.uniform(param_ranges[i][0], param_ranges[i][1])
        for i in range(n_params)
    ])
    
    for i in range(n_params):
        # 计算基本效应
        x_plus = x_base.copy()
        x_minus = x_base.copy()
        
        param_range = param_ranges[i][1] - param_ranges[i][0]
        x_plus[i] = min(x_base[i] + delta * param_range, param_ranges[i][1])
        x_minus[i] = max(x_base[i] - delta * param_range, param_ranges[i][0])
        
        y_plus = ignition_delay(x_plus[0], x_plus[1], x_plus[2])
        y_minus = ignition_delay(x_minus[0], x_minus[1], x_minus[2])
        
        EE_matrix[t, i] = (y_plus - y_minus) / (x_plus[i] - x_minus[i])
        EE_matrix[t, i] *= x_base[i] / ignition_delay(x_base[0], x_base[1], x_base[2])

# 统计量
mu = np.mean(EE_matrix, axis=0)
sigma = np.std(EE_matrix, axis=0)

print("参数\t\t\t 均值(μ)\t 标准差(σ)")
print("-" * 50)
for i, name in enumerate(param_names):
    print(f"{name:20s}\t{mu[i]:8.3f}\t{sigma[i]:8.3f}")

# 绘制敏感性分析结果
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 局部敏感性
params = ['Temperature', 'Pressure', 'Equivalence Ratio']
local_sens = [abs(S_T), abs(S_P), abs(S_phi)]

axes[0].bar(params, local_sens, color=['red', 'blue', 'green'], edgecolor='black')
axes[0].set_ylabel('Sensitivity Coefficient |S|', fontsize=11)
axes[0].set_title('Local Sensitivity Analysis', fontsize=12)
axes[0].grid(True, alpha=0.3, axis='y')
axes[0].tick_params(axis='x', rotation=15)

# Morris散点图
axes[1].scatter(mu, sigma, s=200, c=['red', 'blue', 'green'], alpha=0.7, edgecolors='black')
for i, name in enumerate(param_names):
    axes[1].annotate(name[:4], (mu[i], sigma[i]), xytext=(5, 5), 
                     textcoords='offset points', fontsize=10)
axes[1].axhline(y=0.1, color='gray', linestyle='--', alpha=0.5)
axes[1].axvline(x=0, color='gray', linestyle='--', alpha=0.5)
axes[1].set_xlabel('Mean (μ)', fontsize=11)
axes[1].set_ylabel('Standard Deviation (σ)', fontsize=11)
axes[1].set_title('Morris Sensitivity Analysis', fontsize=12)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('exp2_sensitivity_analysis.png', dpi=150, bbox_inches='tight')
print("\n结果已保存至: exp2_sensitivity_analysis.png")

plt.close()

6.3 蒙特卡洛不确定性传播仿真

# ============================================================================
# 实验3: 蒙特卡洛不确定性传播
# ============================================================================
print("\n" + "=" * 70)
print("实验3: 蒙特卡洛不确定性传播")
print("=" * 70)

def combustion_model(T, P, phi, X_O2=0.21):
    """
    简化的燃烧性能模型
    
    返回多个输出量
    """
    # 点火延迟时间
    tau = ignition_delay(T, P, phi)
    
    # 火焰速度
    S_l = flame_speed(T, P, phi)
    
    # 绝热火焰温度（简化计算）
    T_ad = T + 2000 * phi * np.exp(-0.5*(phi-1)**2) * (P/1e5)**0.1
    
    # 燃烧效率
    eta = 0.95 * np.exp(-0.1*abs(phi-1)) * (T/300)**0.05
    
    return tau, S_l, T_ad, eta

# 定义输入不确定性
print("\n输入参数不确定性:")
print("-" * 50)

# 温度: 正态分布 N(300, 15²)
T_dist = stats.norm(300, 15)

# 压力: 正态分布 N(101325, 2000²)
P_dist = stats.norm(101325, 2000)

# 当量比: 三角分布 T(0.8, 1.0, 1.2)
phi_dist = stats.triang(c=0.5, loc=0.8, scale=0.4)

# 氧气浓度: 均匀分布 U(0.19, 0.23)
XO2_dist = stats.uniform(0.19, 0.04)

print("温度: N(300, 15²) K")
print("压力: N(101325, 2000²) Pa")
print("当量比: Triangular(0.8, 1.0, 1.2)")
print("氧浓度: U(0.19, 0.23)")

# 蒙特卡洛模拟
n_mc = 5000
print(f"\n蒙特卡洛模拟 (n={n_mc}):")
print("-" * 50)

# 生成输入样本
T_samples = T_dist.rvs(n_mc)
P_samples = P_dist.rvs(n_mc)
phi_samples = phi_dist.rvs(n_mc)
XO2_samples = XO2_dist.rvs(n_mc)

# 计算输出
tau_results = []
Sl_results = []
Tad_results = []
eta_results = []

for i in range(n_mc):
    tau, S_l, T_ad, eta = combustion_model(T_samples[i], P_samples[i], phi_samples[i], XO2_samples[i])
    tau_results.append(tau)
    Sl_results.append(S_l)
    Tad_results.append(T_ad)
    eta_results.append(eta)

tau_results = np.array(tau_results)
Sl_results = np.array(Sl_results)
Tad_results = np.array(Tad_results)
eta_results = np.array(eta_results)

# 统计分析
print("\n输出统计量:")
print("-" * 50)
print(f"{'输出量':<20s} {'均值':<12s} {'标准差':<12s} {'CV(%)':<8s}")
print("-" * 50)
print(f"{'点火延迟时间(s)':<20s} {np.mean(tau_results):<12.6f} {np.std(tau_results):<12.6f} {np.std(tau_results)/np.mean(tau_results)*100:<8.2f}")
print(f"{'火焰速度(m/s)':<20s} {np.mean(Sl_results):<12.4f} {np.std(Sl_results):<12.4f} {np.std(Sl_results)/np.mean(Sl_results)*100:<8.2f}")
print(f"{'绝热火焰温度(K)':<20s} {np.mean(Tad_results):<12.2f} {np.std(Tad_results):<12.2f} {np.std(Tad_results)/np.mean(Tad_results)*100:<8.2f}")
print(f"{'燃烧效率':<20s} {np.mean(eta_results):<12.4f} {np.std(eta_results):<12.4f} {np.std(eta_results)/np.mean(eta_results)*100:<8.2f}")

# 置信区间
print("\n95%置信区间:")
print("-" * 50)
print(f"点火延迟时间: [{np.percentile(tau_results, 2.5):.6f}, {np.percentile(tau_results, 97.5):.6f}] s")
print(f"火焰速度: [{np.percentile(Sl_results, 2.5):.4f}, {np.percentile(Sl_results, 97.5):.4f}] m/s")
print(f"绝热火焰温度: [{np.percentile(Tad_results, 2.5):.2f}, {np.percentile(Tad_results, 97.5):.2f}] K")
print(f"燃烧效率: [{np.percentile(eta_results, 2.5):.4f}, {np.percentile(eta_results, 97.5):.4f}]")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 点火延迟时间
axes[0, 0].hist(tau_results, bins=50, density=True, alpha=0.7, color='blue', edgecolor='black')
axes[0, 0].axvline(np.mean(tau_results), color='r', linestyle='--', linewidth=2, label=f'Mean={np.mean(tau_results):.6f}')
axes[0, 0].axvline(np.percentile(tau_results, 2.5), color='g', linestyle=':', linewidth=2, label='95% CI')
axes[0, 0].axvline(np.percentile(tau_results, 97.5), color='g', linestyle=':', linewidth=2)
axes[0, 0].set_xlabel('Ignition Delay Time (s)', fontsize=11)
axes[0, 0].set_ylabel('Probability Density', fontsize=11)
axes[0, 0].set_title('Ignition Delay Time Distribution', fontsize=12)
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 火焰速度
axes[0, 1].hist(Sl_results, bins=50, density=True, alpha=0.7, color='green', edgecolor='black')
axes[0, 1].axvline(np.mean(Sl_results), color='r', linestyle='--', linewidth=2, label=f'Mean={np.mean(Sl_results):.4f}')
axes[0, 1].axvline(np.percentile(Sl_results, 2.5), color='g', linestyle=':', linewidth=2, label='95% CI')
axes[0, 1].axvline(np.percentile(Sl_results, 97.5), color='g', linestyle=':', linewidth=2)
axes[0, 1].set_xlabel('Flame Speed (m/s)', fontsize=11)
axes[0, 1].set_ylabel('Probability Density', fontsize=11)
axes[0, 1].set_title('Flame Speed Distribution', fontsize=12)
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# 绝热火焰温度
axes[1, 0].hist(Tad_results, bins=50, density=True, alpha=0.7, color='red', edgecolor='black')
axes[1, 0].axvline(np.mean(Tad_results), color='b', linestyle='--', linewidth=2, label=f'Mean={np.mean(Tad_results):.2f}')
axes[1, 0].axvline(np.percentile(Tad_results, 2.5), color='g', linestyle=':', linewidth=2, label='95% CI')
axes[1, 0].axvline(np.percentile(Tad_results, 97.5), color='g', linestyle=':', linewidth=2)
axes[1, 0].set_xlabel('Adiabatic Flame Temperature (K)', fontsize=11)
axes[1, 0].set_ylabel('Probability Density', fontsize=11)
axes[1, 0].set_title('Adiabatic Flame Temperature Distribution', fontsize=12)
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# 燃烧效率
axes[1, 1].hist(eta_results, bins=50, density=True, alpha=0.7, color='purple', edgecolor='black')
axes[1, 1].axvline(np.mean(eta_results), color='r', linestyle='--', linewidth=2, label=f'Mean={np.mean(eta_results):.4f}')
axes[1, 1].axvline(np.percentile(eta_results, 2.5), color='g', linestyle=':', linewidth=2, label='95% CI')
axes[1, 1].axvline(np.percentile(eta_results, 97.5), color='g', linestyle=':', linewidth=2)
axes[1, 1].set_xlabel('Combustion Efficiency', fontsize=11)
axes[1, 1].set_ylabel('Probability Density', fontsize=11)
axes[1, 1].set_title('Combustion Efficiency Distribution', fontsize=12)
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('exp3_monte_carlo.png', dpi=150, bbox_inches='tight')
print("\n结果已保存至: exp3_monte_carlo.png")

plt.close()

6.4 响应面与代理模型仿真

# ============================================================================
# 实验4: 响应面方法
# ============================================================================
print("\n" + "=" * 70)
print("实验4: 响应面方法")
print("=" * 70)

from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score

# 定义火焰速度的响应面
print("\n构建火焰速度响应面:")
print("-" * 50)

# 生成训练数据（使用LHS设计）
from pyDOE2 import lhs

n_train = 50
# 拉丁超立方采样
lhs_samples = lhs(2, samples=n_train)

# 将样本映射到参数空间
T_train = 250 + lhs_samples[:, 0] * 100  # 250-350 K
phi_train = 0.6 + lhs_samples[:, 1] * 0.8  # 0.6-1.4

# 计算训练输出
Sl_train = np.array([flame_speed(T, 101325, phi) for T, phi in zip(T_train, phi_train)])

# 添加噪声模拟实验误差
noise = np.random.normal(0, 0.02, len(Sl_train))
Sl_train_noisy = Sl_train + noise

# 构建多项式响应面
poly = PolynomialFeatures(degree=2)
X_train = np.column_stack([T_train, phi_train])
X_train_poly = poly.fit_transform(X_train)

# 拟合模型
model = LinearRegression()
model.fit(X_train_poly, Sl_train_noisy)

# 预测
Sl_pred = model.predict(X_train_poly)

print(f"响应面模型: 二阶多项式")
print(f"R² = {r2_score(Sl_train_noisy, Sl_pred):.4f}")

# 生成预测网格
T_grid = np.linspace(250, 350, 50)
phi_grid = np.linspace(0.6, 1.4, 50)
T_mesh, phi_mesh = np.meshgrid(T_grid, phi_grid)

# 响应面预测
X_grid = np.column_stack([T_mesh.ravel(), phi_mesh.ravel()])
X_grid_poly = poly.transform(X_grid)
Sl_mesh = model.predict(X_grid_poly).reshape(T_mesh.shape)

# 真实值
Sl_true = np.array([flame_speed(T, 101325, phi) for T, phi in zip(T_mesh.ravel(), phi_mesh.ravel())])
Sl_true_mesh = Sl_true.reshape(T_mesh.shape)

# 不确定性传播
print("\n使用响应面进行不确定性传播:")
print("-" * 50)

# 定义输入分布
T_dist_rsp = stats.norm(300, 15)
phi_dist_rsp = stats.triang(c=0.5, loc=0.8, scale=0.4)

# 蒙特卡洛（使用响应面）
n_rsp = 10000
T_samples_rsp = T_dist_rsp.rvs(n_rsp)
phi_samples_rsp = phi_dist_rsp.rvs(n_rsp)

X_samples = np.column_stack([T_samples_rsp, phi_samples_rsp])
X_samples_poly = poly.transform(X_samples)
Sl_samples_rsp = model.predict(X_samples_poly)

print(f"响应面MC结果:")
print(f"  均值: {np.mean(Sl_samples_rsp):.4f} m/s")
print(f"  标准差: {np.std(Sl_samples_rsp):.4f} m/s")
print(f"  95% CI: [{np.percentile(Sl_samples_rsp, 2.5):.4f}, {np.percentile(Sl_samples_rsp, 97.5):.4f}]")

# 与直接MC比较（少量样本）
n_direct = 500
T_direct = T_dist_rsp.rvs(n_direct)
phi_direct = phi_dist_rsp.rvs(n_direct)
Sl_direct = np.array([flame_speed(T, 101325, phi) for T, phi in zip(T_direct, phi_direct)])

print(f"\n直接MC结果 (n={n_direct}):")
print(f"  均值: {np.mean(Sl_direct):.4f} m/s")
print(f"  标准差: {np.std(Sl_direct):.4f} m/s")
print(f"  95% CI: [{np.percentile(Sl_direct, 2.5):.4f}, {np.percentile(Sl_direct, 97.5):.4f}]")

print(f"\n计算效率:")
print(f"  响应面MC: {n_rsp}次预测（几乎瞬时）")
print(f"  直接MC: {n_direct}次模拟（耗时）")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 响应面3D图
ax1 = axes[0, 0]
contour1 = ax1.contourf(T_mesh, phi_mesh, Sl_mesh, levels=20, cmap='viridis')
ax1.scatter(T_train, phi_train, c='red', s=30, edgecolors='black', label='Training points')
ax1.set_xlabel('Temperature (K)', fontsize=11)
ax1.set_ylabel('Equivalence Ratio', fontsize=11)
ax1.set_title('Response Surface: Flame Speed', fontsize=12)
plt.colorbar(contour1, ax=ax1, label='Flame Speed (m/s)')
ax1.legend()

# 真实值对比
ax2 = axes[0, 1]
contour2 = ax2.contourf(T_mesh, phi_mesh, Sl_true_mesh, levels=20, cmap='viridis')
ax2.set_xlabel('Temperature (K)', fontsize=11)
ax2.set_ylabel('Equivalence Ratio', fontsize=11)
ax2.set_title('True Model: Flame Speed', fontsize=12)
plt.colorbar(contour2, ax=ax2, label='Flame Speed (m/s)')

# 残差图
ax3 = axes[1, 0]
residuals = Sl_train_noisy - Sl_pred
ax3.scatter(Sl_pred, residuals, alpha=0.6, edgecolors='black')
ax3.axhline(y=0, color='r', linestyle='--')
ax3.set_xlabel('Predicted Values', fontsize=11)
ax3.set_ylabel('Residuals', fontsize=11)
ax3.set_title('Residual Plot', fontsize=12)
ax3.grid(True, alpha=0.3)

# 不确定性传播结果对比
ax4 = axes[1, 1]
ax4.hist(Sl_samples_rsp, bins=50, density=True, alpha=0.5, color='blue', 
         label=f'RSM (n={n_rsp})', edgecolor='black')
ax4.hist(Sl_direct, bins=30, density=True, alpha=0.5, color='red', 
         label=f'Direct (n={n_direct})', edgecolor='black')
ax4.set_xlabel('Flame Speed (m/s)', fontsize=11)
ax4.set_ylabel('Probability Density', fontsize=11)
ax4.set_title('Uncertainty Propagation Comparison', fontsize=12)
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('exp4_response_surface.png', dpi=150, bbox_inches='tight')
print("\n结果已保存至: exp4_response_surface.png")

plt.close()

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