第033篇：燃烧模型验证与确认

摘要

燃烧模型的验证与确认（Verification and Validation, V&V）是确保数值模拟结果可靠性和工程应用价值的关键环节。验证（Verification）关注"是否正确地求解了数学模型"，确认（Validation）关注"数学模型是否正确描述了物理现象"。本教程系统阐述燃烧仿真中V&V的理论基础、方法论体系和实践流程，涵盖误差分析、网格收敛性研究、不确定性量化、基准测试案例设计等核心技术。通过Python仿真实现，演示如何开展系统的模型验证与确认工作，包括数值误差估计、实验数据对比、敏感性分析和可信度评估。本教程为燃烧仿真工程师提供完整的V&V实施指南，帮助建立可靠的数值模拟能力。

关键词

燃烧模型验证与确认, 数值误差分析, 网格收敛性, 不确定性量化, 基准测试, 模型可信度, 实验对比, 敏感性分析

---

1. 引言

1.1 验证与确认的重要性

在燃烧数值模拟领域，模型验证与确认是连接理论研究与工程应用的桥梁。随着计算流体力学（CFD）技术的快速发展，燃烧仿真已成为发动机设计、工业燃烧器优化和火灾安全评估的重要工具。然而，数值模拟结果的可靠性始终是工程应用中的核心关切。一个未经充分验证的燃烧模型可能导致错误的设计决策，造成巨大的经济损失甚至安全事故。

验证与确认（V&V）的概念最早由美国计算机仿真学会（SCS）于1979年提出，经过四十余年的发展，已形成完整的理论体系和方法论框架。美国机械工程师学会（ASME）发布的V&V 20标准和美国航空航天学会（AIAA）的G-077指南为燃烧仿真V&V提供了规范性指导。这些标准强调：验证回答"我们是否正确地求解了方程"，确认回答"我们是否求解了正确的方程"。

燃烧过程的复杂性使得V&V工作面临独特挑战：

多物理场耦合：燃烧涉及流体流动、传热传质、化学反应和辐射等多个物理过程的强耦合，每个子模型都需要独立验证后再进行整体确认。

多尺度特性：从分子尺度的化学反应到设备尺度的流动结构，时间尺度跨越十几个数量级，难以在单一模拟中完整捕捉。

强非线性：燃烧过程中的温度指数依赖性和化学反应链式机制导致系统对初始条件和边界条件极为敏感。

实验测量困难：高温、高压和瞬态特性使得燃烧实验的精确测量极具挑战性，限制了确认数据的获取。

1.2 V&V的基本概念

1.2.1 验证（Verification）

验证是评估数值求解精度的过程，确保计算代码正确地实现了数学模型。验证分为两个层次：

代码验证（Code Verification）：检查数值算法是否正确实现，包括：

• 空间离散精度验证
• 时间积分精度验证
• 边界条件实现验证
• 并行计算正确性验证

计算验证（Calculation Verification）：评估特定算例的数值误差，包括：

• 离散化误差估计
• 迭代误差评估
• 舍入误差分析
• 网格收敛性研究

1.2.2 确认（Validation）

确认是评估数学模型预测物理现象准确性的过程，通过将数值结果与实验数据对比来实现。确认的核心问题包括：

模型适用性：数学模型是否包含了所有重要的物理机制？

参数校准：模型中的经验参数是否经过合理标定？

预测能力：模型能否准确预测未参与校准的工况？

不确定性量化：模型预测的不确定性范围是多少？

1.3 V&V方法论框架

完整的V&V流程遵循以下逻辑结构：

概念模型 → 数学模型 → 离散模型 → 计算机代码 → 数值解
    ↓           ↓           ↓            ↓           ↓
  概念确认    模型确认    代码验证    计算验证    实验对比

概念确认（Conceptual Validation）：评估概念模型是否完整描述了物理系统的关键特征。

模型确认（Model Validation）：评估数学模型是否准确描述了概念模型。

代码验证（Code Verification）：证明计算机代码正确求解了离散方程。

计算验证（Calculation Verification）：估计特定数值解的误差。

实验确认（Experimental Validation）：将数值结果与实验数据对比，评估模型的预测能力。

1.4 燃烧仿真V&V的特殊性

燃烧仿真V&V相比一般CFD具有独特挑战：

化学反应机理的不确定性：详细的化学反应机理包含数百个组分和数千个反应，其速率常数往往存在较大不确定性。简化机理的构建和验证是燃烧V&V的重要环节。

湍流-燃烧相互作用：湍流与化学反应的强耦合使得亚格子尺度模型难以验证。直接数值模拟（DNS）为大涡模拟（LES）和雷诺平均（RANS）模型提供了验证基准。

火焰结构的多尺度性：层流火焰厚度通常在毫米量级，而工业燃烧器尺度可达米级，要求自适应网格或火焰面追踪技术。

数值刚性与稳定性：化学反应的特征时间尺度远小于流动时间尺度，导致方程组刚性，需要特殊的数值处理方法。

---

2. 验证理论基础

2.1 数值误差分析

数值模拟中的总误差可分解为：

$${ϵ}_{total}={ϵ}_{model}+{ϵ}_{numerical}+{ϵ}_{input}$$

其中：

• ${ϵ}_{model}$：模型误差，源于数学模型对物理现象的简化
• ${ϵ}_{numerical}$：数值误差，包括离散误差、迭代误差和舍入误差
• ${ϵ}_{input}$：输入误差，来自边界条件、物性参数等的不确定性

2.1.1 离散误差

离散误差是数值解与精确解之间的差异，由空间和时间离散化引起。对于具有$p$阶精度的数值格式，离散误差可表示为：

$${ϵ}_h={\varphi }_h-{\varphi }_{exact}=C\cdot h^p+O(h^{p+1})$$

其中$h$为特征网格尺寸，$C$为误差常数，$p$为收敛阶数。

Richardson外推法：利用不同网格上的数值解估计精确解和离散误差。

$${\varphi }_{exact}\approx {\varphi }_{h_2}+\frac{{\varphi }_{h_2}-{\varphi }_{h_1}}{r^p-1}$$

其中$r=h_1/h_2$为网格细化比，通常取$r=2$。

网格收敛指数（GCI）：基于Richardson外推的误差估计方法，提供保守的误差界限。

$$GCI=\frac{F_s\cdot |ϵ|}{r^p-1}$$

其中$F_s$为安全因子，通常取$F_s=1.25$。

2.1.2 迭代误差

迭代误差是迭代求解过程中当前解与收敛解之间的差异：

$${ϵ}_{iter}={\varphi }^{(n)}-{\varphi }^{(\infty )}$$

对于线性迭代方法，误差随迭代次数呈指数衰减：

$$\Vert {ϵ}^{(n)}\Vert \le {\rho }^n\Vert {ϵ}^{(0)}\Vert$$

其中$\rho$为迭代矩阵的谱半径。

收敛判据：通常采用残差下降若干个数量级或解的变化量小于阈值作为收敛标准。

$$\frac{\Vert {\varphi }^{(n+1)}-{\varphi }^{(n)}\Vert }{\Vert {\varphi }^{(n)}\Vert }<{ϵ}_{tol}$$

2.1.3 舍入误差

舍入误差由计算机有限精度算术引起，对于双精度浮点数（约16位有效数字），相对舍入误差约为$10^{-16}$。在燃烧仿真中，舍入误差通常远小于离散误差，但在以下情况需要关注：

• 强刚性方程组的长时间积分
• 大梯度区域的差分计算
• 病态矩阵的求解

2.2 收敛阶数验证

收敛阶数是评估数值格式精度的关键指标。通过在不同网格上求解同一问题，可以验证数值格式的理论收敛阶数。

空间收敛阶数：

$$p=\frac{\log \left(\frac{\Vert {ϵ}_{h_1}\Vert }{\Vert {ϵ}_{h_2}\Vert }\right)}{\log (r)}$$

时间收敛阶数：

$$q=\frac{\log \left(\frac{\Vert {ϵ}_{\Delta t_1}\Vert }{\Vert {ϵ}_{\Delta t_2}\Vert }\right)}{\log (r_t)}$$

其中$r_t=\Delta t_1/\Delta t_2$为时间步长比。

2.3 制造解方法（Method of Manufactured Solutions, MMS）

MMS是代码验证的强大工具，通过构造满足控制方程的解析解来验证数值实现。

基本步骤：

1. 选择试验函数：选择足够光滑的解析函数作为"制造解"
   $${\varphi }_{MMS}(x,y,z,t)=f(x,y,z,t)$$

2. 计算源项：将制造解代入控制方程，计算所需的源项
   $$S=\frac{\partial {\varphi }_{MMS}}{\partial t}+\nabla \cdot (\mathbf{u}{\varphi }_{MMS})-\nabla \cdot (D\nabla {\varphi }_{MMS})$$

3. 数值求解：在计算域上求解带源项的控制方程

4. 误差分析：比较数值解与制造解，验证收敛阶数

MMS在燃烧仿真中的应用：

• 验证对流-扩散-反应方程的离散格式
• 检验复杂边界条件的实现
• 验证组分输运方程的求解器
• 检验化学反应源项的数值处理

---

3. 确认方法论

3.1 确认层级结构

燃烧模型确认采用分层递进策略：

单元问题（Unit Problems）：验证单个物理模型或子程序，如：

• 化学反应机理验证
• 物性计算子程序验证
• 湍流模型常数标定

基准案例（Benchmark Cases）：标准测试问题，具有解析解或高精度参考解：

• 一维层流火焰
• 二维涡-火焰相互作用
• 激波管点火延迟

子系统案例（Subsystem Cases）：简化几何的复杂物理问题：

• 对冲扩散火焰
• 旋流燃烧器冷态流动
• 射流火焰

完整系统（Complete Systems）：实际工程设备：

• 航空发动机燃烧室
• 工业燃气轮机
• 锅炉燃烧器

3.2 实验确认流程

3.2.1 实验设计原则

目标明确性：每个实验应有明确的确认目标，针对特定的模型假设或参数。

边界条件可控：实验边界条件应精确测量并可用于数值模拟。

测量完整性：应测量足够的物理量以全面评估模型性能。

不确定性量化：实验数据应包含测量不确定度估计。

3.2.2 确认度量指标

误差范数：

$$L_2\text{范数}:E_{L2}=\sqrt{\frac{1}{N}\sum _{i=1}^N({\varphi }_{sim,i}-{\varphi }_{exp,i}{)}^2}$$

$$L_{\infty }\text{范数}:E_{L\infty }=\underset{i}{\max }|{\varphi }_{sim,i}-{\varphi }_{exp,i}|$$

相对误差：

$$E_{rel}=\frac{\Vert {\varphi }_{sim}-{\varphi }_{exp}\Vert }{\Vert {\varphi }_{exp}\Vert }$$

确认度量（Validation Metric）：

$$d=\frac{|{\bar{\varphi }}_{sim}-{\bar{\varphi }}_{exp}|}{\sqrt{u_{sim}^2+u_{exp}^2}}$$

其中$\bar{\varphi }$为均值，$u$为不确定度。当$d<1$时，认为模型通过确认。

3.3 不确定性量化

3.3.1 不确定性来源

偶然不确定性（Aleatory Uncertainty）：固有的随机性，如：

• 湍流脉动
• 燃料组分波动
• 来流条件变化

认知不确定性（Epistemic Uncertainty）：知识缺乏导致，如：

• 模型参数不确定
• 边界条件测量误差
• 数值误差

3.3.2 不确定性传播方法

蒙特卡洛方法：通过大量随机采样传播不确定性。

$${\varphi }_{output}=f({\xi }_1,{\xi }_2,...,{\xi }_n)$$

其中${\xi }_i$为输入随机变量。

多项式混沌展开：利用正交多项式基函数表示随机输出。

$$\varphi (\xi )=\sum _{k=0}^P{\alpha }_k{\Psi }_k(\xi )$$

敏感性分析：评估输入参数对输出的影响程度。

Sobol指数：

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

其中$V$为方差，$E$为期望。

---

4. 燃烧模型验证实践

4.1 化学反应机理验证

4.1.1 机理简化与验证

详细化学反应机理包含数百个组分和数千个反应，直接用于CFD计算成本过高。机理简化是燃烧仿真的重要预处理步骤。

简化方法：

1. 敏感性分析：识别对目标量影响最大的反应
2. 准稳态近似（QSSA）：假设某些中间组分处于准稳态
3. 部分平衡假设：快速反应达到局部平衡
4. 直接关系图（DRG）：基于组分耦合度进行删减

验证指标：

• 点火延迟时间
• 层流火焰速度
• 火焰温度
• 主要组分浓度分布

4.1.2 机理比较与选择

常用燃烧机理数据库：

• GRI-Mech：天然气燃烧，约50种组分
• USC Mech：碳氢燃料，约100种组分
• LLNL机理：多组分燃料，详细机理
• Aramco Mech：宽范围燃料适用

4.2 湍流燃烧模型验证

4.2.1 RANS模型验证

标准验证案例：

1. Sandia火焰D/E/F：湍流扩散射流火焰
2. Sydney swirl flames：旋流燃烧
3. TNF Workshop案例：标准化测试火焰

验证指标：

• 平均温度和速度分布
• 湍流脉动强度
• 组分平均浓度和均方根值
• 混合分数分布

4.2.2 LES模型验证

DNS数据库对比：

• Sandia/TUD DNS：湍流射流火焰
• Cambridge DNS：层流-湍流转变
• Sandia DNS：喷射点火

验证指标：

• 瞬时火焰结构
• 条件平均统计量
• 火焰面密度分布
• 点火概率

4.3 辐射模型验证

辐射传热在高温燃烧中至关重要，特别是对于富燃料火焰和炉膛计算。

验证方法：

1. 解析解对比：简单几何的精确解
2. 蒙特卡洛参考解：高精度统计解
3. 实验测量：辐射热流测量

验证指标：

• 辐射源项分布
• 壁面辐射热流
• 介质温度分布

---

5. 基准测试案例设计

5.1 一维层流火焰

一维层流预混火焰是燃烧模型验证的基础案例。

控制方程：

连续性：
$$\frac{d(\rho u)}{dx}=0$$

能量：
$$\rho u{c}_p\frac{dT}{dx}=\frac{d}{dx}\left(\lambda \frac{dT}{dx}\right)-\sum _k{h}_k{\dot{\omega }}_k{W}_k-{\dot{q}}_{rad}$$

组分：
$$\rho u\frac{d{Y}_k}{dx}=-\frac{d}{dx}\left(\rho Y_k{V}_k\right)+{\dot{\omega }}_k{W}_k$$

验证参数：

• 层流火焰速度$S_L$
• 火焰厚度${\delta }_L$
• 绝热火焰温度$T_{ad}$
• 组分浓度分布

5.2 对冲扩散火焰

对冲扩散火焰用于验证非预混燃烧模型。

特征参数：

• 应变率：$a=\frac{2|U_{jet}|}{L}$
• 标量耗散率：$\chi =2D|\nabla Z{|}^2$
• 熄火应变率：$a_{ext}$

验证指标：

• 温度峰值
• 组分分布
• 熄火特性

5.3 激波管点火延迟

激波管实验提供高温高压下的点火延迟数据，用于验证化学反应机理。

点火延迟定义：

$${\tau }_{ign}=t(\max (dT/dt))-t_{shockpassage}$$

验证范围：

• 温度：1000-2500 K
• 压力：1-50 atm
• 当量比：0.5-2.0

---

6. Python仿真实现

6.1 数值误差分析仿真

以下代码演示如何分析数值格式的收敛阶数和估计离散误差。

# -*- coding: utf-8 -*-
"""
实验1: 数值误差分析与收敛阶数验证
演示Richardson外推法和GCI计算
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import warnings
warnings.filterwarnings('ignore')
import os

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

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

print("=" * 60)
print("实验1: 数值误差分析与收敛阶数验证")
print("=" * 60)

# 定义精确解（制造解）
def exact_solution(x, t):
    """对流-扩散方程的解析解"""
    return np.exp(-((x - 0.5 - 0.1*t)**2) / (0.01 + 0.004*t))

# 数值求解函数
def solve_convection_diffusion(N, dt, t_final, scheme='upwind'):
    """
    求解一维对流-扩散方程
    d(phi)/dt + u*d(phi)/dx = D*d2(phi)/dx2
    """
    dx = 1.0 / N
    x = np.linspace(0, 1, N+1)
    u = 0.1  # 对流速度
    D = 0.001  # 扩散系数
    
    # 初始条件
    phi = np.exp(-((x - 0.3)**2) / 0.01)
    
    # 时间推进
    n_steps = int(t_final / dt)
    
    for n in range(n_steps):
        phi_new = phi.copy()
        
        for i in range(1, N):
            if scheme == 'upwind':
                # 一阶迎风格式
                conv = -u * (phi[i] - phi[i-1]) / dx if u > 0 else -u * (phi[i+1] - phi[i]) / dx
            elif scheme == 'central':
                # 二阶中心格式
                conv = -u * (phi[i+1] - phi[i-1]) / (2*dx)
            
            # 二阶中心差分扩散项
            diff = D * (phi[i+1] - 2*phi[i] + phi[i-1]) / dx**2
            
            phi_new[i] = phi[i] + dt * (conv + diff)
        
        # 边界条件
        phi_new[0] = phi_new[1]
        phi_new[N] = phi_new[N-1]
        
        phi = phi_new
    
    return x, phi

# 网格收敛性研究
print("\n[1] 网格收敛性研究")
print("-" * 50)

grids = [20, 40, 80, 160, 320]
errors_l2 = []
errors_linf = []

# 参考解（最细网格）
x_ref, phi_ref = solve_convection_diffusion(640, 0.0001, 1.0, 'upwind')

for N in grids:
    x, phi = solve_convection_diffusion(N, 0.0001, 1.0, 'upwind')
    # 插值到参考网格
    phi_interp = np.interp(x_ref, x, phi)
    
    error_l2 = np.sqrt(np.mean((phi_interp - phi_ref)**2))
    error_linf = np.max(np.abs(phi_interp - phi_ref))
    
    errors_l2.append(error_l2)
    errors_linf.append(error_linf)
    
    print(f"  网格数 N={N:3d}: L2误差={error_l2:.6e}, Linf误差={error_linf:.6e}")

# 计算收敛阶数
p_l2 = np.polyfit(np.log(1/np.array(grids)), np.log(errors_l2), 1)[0]
p_linf = np.polyfit(np.log(1/np.array(grids)), np.log(errors_linf), 1)[0]

print(f"\n  收敛阶数: L2范数 p={p_l2:.2f}, Linf范数 p={p_linf:.2f}")
print(f"  (理论值: 一阶迎风格式 p=1.0)")

# 绘制收敛曲线
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 误差-网格关系
ax1 = axes[0]
ax1.loglog(1/np.array(grids), errors_l2, 'bo-', linewidth=2, markersize=8, label='L2 Error')
ax1.loglog(1/np.array(grids), errors_linf, 'rs-', linewidth=2, markersize=8, label='Linf Error')

# 参考线
h_ref = np.array([1/320, 1/20])
ax1.loglog(h_ref, h_ref * errors_l2[-1] / (1/grids[-1]), 'k--', linewidth=1.5, label='1st Order Ref')

ax1.set_xlabel('Grid Spacing h', fontsize=11)
ax1.set_ylabel('Error', fontsize=11)
ax1.set_title('Grid Convergence Study', fontsize=12, fontweight='bold')
ax1.legend(loc='upper left')
ax1.grid(True, which='both', linestyle='--', alpha=0.7)

# Richardson外推
ax2 = axes[1]
N_coarse, N_fine = 80, 160
x_c, phi_c = solve_convection_diffusion(N_coarse, 0.0001, 1.0, 'upwind')
x_f, phi_f = solve_convection_diffusion(N_fine, 0.0001, 1.0, 'upwind')

# Richardson外推解
r = 2  # 细化比
p = 1  # 一阶格式
phi_extrap = (r**p * phi_f[::2] - phi_c) / (r**p - 1)

ax2.plot(x_c, phi_c, 'b--', linewidth=2, label=f'Coarse (N={N_coarse})')
ax2.plot(x_f, phi_f, 'g-', linewidth=2, label=f'Fine (N={N_fine})')
ax2.plot(x_c, phi_extrap, 'r:', linewidth=2.5, label='Richardson Extrapolation')
ax2.plot(x_ref, phi_ref, 'k-', linewidth=1, alpha=0.5, label='Reference (N=640)')

ax2.set_xlabel('x', fontsize=11)
ax2.set_ylabel('phi', fontsize=11)
ax2.set_title('Richardson Extrapolation', fontsize=12, fontweight='bold')
ax2.legend(loc='upper right')
ax2.grid(True, linestyle='--', alpha=0.7)

plt.tight_layout()
plt.savefig(f'{output_dir}/exp1_error_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n  图1已保存: exp1_error_analysis.png")

print("\n" + "=" * 60)
print("实验1完成")
print("=" * 60)

6.2 制造解方法验证

以下代码演示如何使用制造解方法验证数值代码的正确性。

# -*- coding: utf-8 -*-
"""
实验2: 制造解方法（MMS）验证
验证数值格式的实现正确性
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
import os

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

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

print("\n" + "=" * 60)
print("实验2: 制造解方法（MMS）验证")
print("=" * 60)

# 制造解
def manufactured_solution(x, t):
    """制造解: 光滑的正弦函数"""
    return np.sin(np.pi * x) * np.exp(-t)

def manufactured_source(x, t, u, D):
    """计算制造解对应的源项"""
    phi = manufactured_solution(x, t)
    dphi_dt = -phi
    dphi_dx = np.pi * np.cos(np.pi * x) * np.exp(-t)
    d2phi_dx2 = -(np.pi**2) * phi
    
    # S = d(phi)/dt + u*d(phi)/dx - D*d2(phi)/dx2
    S = dphi_dt + u * dphi_dx - D * d2phi_dx2
    return S

# 带源项的数值求解器
def solve_with_source(N, dt, t_final, u, D):
    """求解带源项的对流-扩散方程"""
    dx = 1.0 / N
    x = np.linspace(0, 1, N+1)
    
    # 初始条件（与制造解一致）
    phi = manufactured_solution(x, 0)
    
    n_steps = int(t_final / dt)
    
    for n in range(n_steps):
        t = (n + 1) * dt
        phi_new = phi.copy()
        
        # 计算源项
        S = manufactured_source(x, t, u, D)
        
        for i in range(1, N):
            # 对流项（一阶迎风格式）
            conv = -u * (phi[i] - phi[i-1]) / dx if u > 0 else -u * (phi[i+1] - phi[i]) / dx
            
            # 扩散项（二阶中心）
            diff = D * (phi[i+1] - 2*phi[i] + phi[i-1]) / dx**2
            
            phi_new[i] = phi[i] + dt * (conv + diff + S[i])
        
        # Dirichlet边界条件
        phi_new[0] = manufactured_solution(0, t)
        phi_new[N] = manufactured_solution(1, t)
        
        phi = phi_new
    
    return x, phi

print("\n[2] 制造解方法验证")
print("-" * 50)

# 参数设置
u = 0.1
D = 0.01
t_final = 0.5

# 不同网格的误差分析
grids = [10, 20, 40, 80, 160]
errors_mms = []

for N in grids:
    x, phi_num = solve_with_source(N, 0.0001, t_final, u, D)
    phi_exact = manufactured_solution(x, t_final)
    
    error = np.sqrt(np.mean((phi_num - phi_exact)**2))
    errors_mms.append(error)
    
    print(f"  网格数 N={N:3d}: L2误差={error:.6e}")

# 计算观察到的收敛阶数
observed_order = np.polyfit(np.log(1/np.array(grids)), np.log(errors_mms), 1)[0]
print(f"\n  观察到的收敛阶数: p={observed_order:.2f}")
print(f"  理论收敛阶数: p=1.0 (一阶迎风格式)")

# 绘制结果
fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))

# 制造解形状
ax1 = axes[0]
x_fine = np.linspace(0, 1, 200)
t_values = [0, 0.25, 0.5]
colors = ['blue', 'green', 'red']

for t, color in zip(t_values, colors):
    phi_mms = manufactured_solution(x_fine, t)
    ax1.plot(x_fine, phi_mms, color=color, linewidth=2, label=f't={t}')

ax1.set_xlabel('x', fontsize=11)
ax1.set_ylabel('phi', fontsize=11)
ax1.set_title('Manufactured Solution', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, linestyle='--', alpha=0.7)

# 数值解与制造解对比
ax2 = axes[1]
N_plot = 40
x, phi_num = solve_with_source(N_plot, 0.0001, t_final, u, D)
phi_exact = manufactured_solution(x, t_final)

ax2.plot(x, phi_exact, 'k-', linewidth=2, label='Manufactured Solution')
ax2.plot(x, phi_num, 'ro--', markersize=6, linewidth=1.5, label=f'Numerical (N={N_plot})')
ax2.set_xlabel('x', fontsize=11)
ax2.set_ylabel('phi', fontsize=11)
ax2.set_title('Numerical vs Manufactured Solution', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, linestyle='--', alpha=0.7)

# 收敛阶数验证
ax3 = axes[2]
ax3.loglog(1/np.array(grids), errors_mms, 'bs-', linewidth=2, markersize=8, label='Observed Error')

# 一阶参考线
h_ref = np.array([1/160, 1/10])
ax3.loglog(h_ref, h_ref * errors_mms[-1] / (1/grids[-1]), 'k--', linewidth=1.5, label='1st Order (p=1.0)')

ax3.set_xlabel('Grid Spacing h', fontsize=11)
ax3.set_ylabel('L2 Error', fontsize=11)
ax3.set_title(f'Order of Accuracy (p={observed_order:.2f})', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, which='both', linestyle='--', alpha=0.7)

plt.tight_layout()
plt.savefig(f'{output_dir}/exp2_manufactured_solution.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n  图2已保存: exp2_manufactured_solution.png")

print("\n" + "=" * 60)
print("实验2完成")
print("=" * 60)

6.3 模型确认与实验对比

以下代码演示如何进行模型确认，将数值结果与"实验数据"进行对比分析。

# -*- coding: utf-8 -*-
"""
实验3: 模型确认与实验对比
模拟燃烧模型的验证与确认流程
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import warnings
warnings.filterwarnings('ignore')
import os

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

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

print("\n" + "=" * 60)
print("实验3: 模型确认与实验对比")
print("=" * 60)

# 模拟实验数据（带噪声和不确定度）
def generate_experimental_data(x, true_model, noise_level=0.05, systematic_error=0.0):
    """生成模拟实验数据"""
    y_true = true_model(x)
    noise = np.random.normal(0, noise_level * np.max(np.abs(y_true)), len(x))
    y_exp = y_true + noise + systematic_error * y_true
    uncertainty = noise_level * np.abs(y_exp) + 0.01 * np.max(np.abs(y_exp))
    return y_exp, uncertainty

# 真实物理模型（用于生成"实验数据"）
def true_combustion_model(x):
    """真实的燃烧模型（复杂反应网络）"""
    T = 300 + 1500 * np.exp(-((x - 0.5)**2) / 0.05) * (1 + 0.1 * np.sin(10*x))
    return T

# 简化数值模型（用于模拟）
def simplified_model(x, model_param=1.0):
    """简化的数值模型（用于CFD）"""
    T = 300 + 1500 * model_param * np.exp(-((x - 0.5)**2) / (0.05 * model_param))
    return T

print("\n[3] 模型确认流程演示")
print("-" * 50)

# 生成实验数据
np.random.seed(42)
x_exp = np.linspace(0, 1, 20)
T_exp, u_exp = generate_experimental_data(x_exp, true_combustion_model, noise_level=0.03)

# 数值模拟结果（不同模型参数）
x_sim = np.linspace(0, 1, 100)
param_values = [0.85, 0.95, 1.0, 1.05, 1.15]

print("\n  模型参数标定:")
print("  " + "-" * 45)

# 计算各参数下的确认度量
validation_metrics = []
for param in param_values:
    T_sim = simplified_model(x_sim, param)
    # 插值到实验点
    T_sim_interp = np.interp(x_exp, x_sim, T_sim)
    
    # 计算误差
    error = np.abs(T_sim_interp - T_exp)
    
    # 确认度量（考虑实验不确定度）
    metric = np.mean(error / u_exp)
    validation_metrics.append(metric)
    
    print(f"  参数={param:.2f}: 平均确认度量={metric:.3f}")

# 最优参数
optimal_param = param_values[np.argmin(validation_metrics)]
print(f"\n  最优模型参数: {optimal_param:.2f}")

# 统计检验
T_sim_optimal = simplified_model(x_sim, optimal_param)
T_sim_interp_opt = np.interp(x_exp, x_sim, T_sim_optimal)

# 计算各种误差指标
mse = np.mean((T_sim_interp_opt - T_exp)**2)
rmse = np.sqrt(mse)
mae = np.mean(np.abs(T_sim_interp_opt - T_exp))
max_error = np.max(np.abs(T_sim_interp_opt - T_exp))

print(f"\n  误差统计 (最优参数):")
print(f"    MSE  = {mse:.2f}")
print(f"    RMSE = {rmse:.2f}")
print(f"    MAE  = {mae:.2f}")
print(f"    Max Error = {max_error:.2f}")

# 决定系数 R^2
ss_res = np.sum((T_exp - T_sim_interp_opt)**2)
ss_tot = np.sum((T_exp - np.mean(T_exp))**2)
r_squared = 1 - ss_res / ss_tot
print(f"    R^2  = {r_squared:.4f}")

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

# 1. 模型与实验对比
ax1 = axes[0, 0]
ax1.errorbar(x_exp, T_exp, yerr=u_exp, fmt='ro', markersize=6, capsize=4, 
             label='Experimental Data', alpha=0.8)

for param in [0.85, 1.0, 1.15]:
    T_sim = simplified_model(x_sim, param)
    linestyle = '-' if param == 1.0 else '--'
    alpha = 1.0 if param == 1.0 else 0.5
    ax1.plot(x_sim, T_sim, linestyle, linewidth=2, alpha=alpha, 
             label=f'Model (param={param})')

ax1.set_xlabel('Position x', fontsize=11)
ax1.set_ylabel('Temperature (K)', fontsize=11)
ax1.set_title('Model Validation: Simulation vs Experiment', fontsize=12, fontweight='bold')
ax1.legend(loc='upper right')
ax1.grid(True, linestyle='--', alpha=0.7)

# 2. 参数敏感性分析
ax2 = axes[0, 1]
ax2.plot(param_values, validation_metrics, 'bo-', linewidth=2, markersize=8)
ax2.axhline(y=1.0, color='r', linestyle='--', linewidth=1.5, label='Validation Threshold')
ax2.axvline(x=optimal_param, color='g', linestyle=':', linewidth=1.5, label=f'Optimal={optimal_param}')
ax2.set_xlabel('Model Parameter', fontsize=11)
ax2.set_ylabel('Validation Metric', fontsize=11)
ax2.set_title('Parameter Sensitivity Analysis', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, linestyle='--', alpha=0.7)

# 3. 残差分析
ax3 = axes[1, 0]
residuals = T_exp - T_sim_interp_opt
ax3.scatter(T_sim_interp_opt, residuals, c='blue', s=60, alpha=0.6)
ax3.axhline(y=0, color='r', linestyle='-', linewidth=1.5)
ax3.axhline(y=np.mean(residuals), color='g', linestyle='--', linewidth=1.5, 
            label=f'Mean={np.mean(residuals):.2f}')
ax3.set_xlabel('Predicted Temperature (K)', fontsize=11)
ax3.set_ylabel('Residuals (K)', fontsize=11)
ax3.set_title('Residual Analysis', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, linestyle='--', alpha=0.7)

# 4. 误差分布
ax4 = axes[1, 1]
n_bins, bins, patches = ax4.hist(residuals, bins=10, density=True, alpha=0.7, color='steelblue', edgecolor='black')

# 拟合正态分布
mu, std = stats.norm.fit(residuals)
x_norm = np.linspace(residuals.min(), residuals.max(), 100)
ax4.plot(x_norm, stats.norm.pdf(x_norm, mu, std), 'r-', linewidth=2, 
         label=f'Normal Fit: μ={mu:.2f}, σ={std:.2f}')

ax4.set_xlabel('Residuals (K)', fontsize=11)
ax4.set_ylabel('Probability Density', fontsize=11)
ax4.set_title('Error Distribution', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, linestyle='--', alpha=0.7)

plt.tight_layout()
plt.savefig(f'{output_dir}/exp3_validation_comparison.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n  图3已保存: exp3_validation_comparison.png")

print("\n" + "=" * 60)
print("实验3完成")
print("=" * 60)

6.4 不确定性量化分析

以下代码演示如何进行不确定性量化和敏感性分析。

# -*- coding: utf-8 -*-
"""
实验4: 不确定性量化与敏感性分析
评估输入参数不确定性对燃烧模拟结果的影响
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import warnings
warnings.filterwarnings('ignore')
import os

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

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

print("\n" + "=" * 60)
print("实验4: 不确定性量化与敏感性分析")
print("=" * 60)

# 燃烧模型（简化的一维火焰温度预测）
def flame_temperature(A, Ea, T0, phi, alpha):
    """
    简化火焰温度模型
    A: 指前因子 (uncertain)
    Ea: 活化能 (uncertain)
    T0: 初始温度
    phi: 当量比
    alpha: 热损失系数 (uncertain)
    """
    R = 8.314
    # 简化的温度模型
    T_ad = T0 + 2000 * np.exp(-Ea / (R * 2000)) * A / 1e8 * (1 - 0.5 * (phi - 1)**2)
    T_flame = T_ad * (1 - alpha)
    return T_flame

print("\n[4] 不确定性量化分析")
print("-" * 50)

# 定义输入参数及其不确定性
# 格式: (均值, 标准差, 分布类型)
params = {
    'A': (1e8, 0.1e8, 'normal'),      # 指前因子
    'Ea': (50000, 2000, 'normal'),    # 活化能 (J/mol)
    'T0': (300, 10, 'normal'),        # 初始温度 (K)
    'phi': (1.0, 0.1, 'normal'),      # 当量比
    'alpha': (0.15, 0.03, 'normal')   # 热损失系数
}

# 蒙特卡洛模拟
n_samples = 10000
np.random.seed(123)

# 生成随机样本
samples = {}
for name, (mean, std, dist) in params.items():
    if dist == 'normal':
        samples[name] = np.random.normal(mean, std, n_samples)
    elif dist == 'uniform':
        samples[name] = np.random.uniform(mean - std, mean + std, n_samples)

# 计算输出分布
T_flame_samples = flame_temperature(
    samples['A'], samples['Ea'], samples['T0'], 
    samples['phi'], samples['alpha']
)

# 统计输出
T_mean = np.mean(T_flame_samples)
T_std = np.std(T_flame_samples)
T_p5 = np.percentile(T_flame_samples, 5)
T_p95 = np.percentile(T_flame_samples, 95)

print(f"\n  火焰温度统计:")
print(f"    均值: {T_mean:.2f} K")
print(f"    标准差: {T_std:.2f} K")
print(f"    95%置信区间: [{T_p5:.2f}, {T_p95:.2f}] K")

# 敏感性分析（相关系数）
print(f"\n  输入参数敏感性 (Pearson相关系数):")
print("  " + "-" * 45)
sensitivities = {}
for name in params.keys():
    corr = np.corrcoef(samples[name], T_flame_samples)[0, 1]
    sensitivities[name] = abs(corr)
    print(f"    {name:8s}: {corr:+.4f}")

# 排序敏感性
sorted_sens = sorted(sensitivities.items(), key=lambda x: x[1], reverse=True)
print(f"\n  参数重要性排序:")
for i, (name, sens) in enumerate(sorted_sens, 1):
    print(f"    {i}. {name}: {sens:.4f}")

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

# 1. 输出分布
ax1 = axes[0, 0]
n_bins, bins, patches = ax1.hist(T_flame_samples, bins=50, density=True, 
                                  alpha=0.7, color='steelblue', edgecolor='black')

# 拟合正态分布
mu, std = stats.norm.fit(T_flame_samples)
x_norm = np.linspace(T_flame_samples.min(), T_flame_samples.max(), 100)
ax1.plot(x_norm, stats.norm.pdf(x_norm, mu, std), 'r-', linewidth=2, 
         label=f'Normal Fit')

# 标记统计量
ax1.axvline(T_mean, color='g', linestyle='-', linewidth=2, label=f'Mean={T_mean:.1f}K')
ax1.axvline(T_p5, color='orange', linestyle='--', linewidth=1.5, label=f'5%={T_p5:.1f}K')
ax1.axvline(T_p95, color='orange', linestyle='--', linewidth=1.5, label=f'95%={T_p95:.1f}K')

ax1.set_xlabel('Flame Temperature (K)', fontsize=11)
ax1.set_ylabel('Probability Density', fontsize=11)
ax1.set_title('Output Uncertainty Distribution', fontsize=12, fontweight='bold')
ax1.legend(loc='upper left')
ax1.grid(True, linestyle='--', alpha=0.7)

# 2. 敏感性条形图
ax2 = axes[0, 1]
names = [s[0] for s in sorted_sens]
values = [s[1] for s in sorted_sens]
colors = plt.cm.viridis(np.linspace(0.2, 0.8, len(names)))

bars = ax2.barh(names, values, color=colors, edgecolor='black')
ax2.set_xlabel('Sensitivity (|Correlation|)', fontsize=11)
ax2.set_ylabel('Parameter', fontsize=11)
ax2.set_title('Parameter Sensitivity Ranking', fontsize=12, fontweight='bold')
ax2.grid(True, axis='x', linestyle='--', alpha=0.7)

# 添加数值标签
for bar, val in zip(bars, values):
    ax2.text(val + 0.01, bar.get_y() + bar.get_height()/2, f'{val:.3f}', 
             va='center', fontsize=9)

# 3. 散点图矩阵（关键参数）
ax3 = axes[1, 0]
key_params = ['Ea', 'alpha', 'phi']
for i, param in enumerate(key_params):
    ax3.scatter(samples[param], T_flame_samples, alpha=0.1, s=1, 
                label=f'{param} (r={sensitivities[param]:.3f})')

ax3.set_xlabel('Input Parameter Value', fontsize=11)
ax3.set_ylabel('Flame Temperature (K)', fontsize=11)
ax3.set_title('Input-Output Relationships', fontsize=12, fontweight='bold')
ax3.legend(markerscale=10)
ax3.grid(True, linestyle='--', alpha=0.7)

# 4. 累积分布函数
ax4 = axes[1, 1]
sorted_T = np.sort(T_flame_samples)
cdf = np.arange(1, len(sorted_T) + 1) / len(sorted_T)

ax4.plot(sorted_T, cdf, 'b-', linewidth=2)
ax4.axhline(0.05, color='r', linestyle='--', alpha=0.7)
ax4.axhline(0.95, color='r', linestyle='--', alpha=0.7)
ax4.axvline(T_p5, color='orange', linestyle=':', alpha=0.7)
ax4.axvline(T_p95, color='orange', linestyle=':', alpha=0.7)

ax4.set_xlabel('Flame Temperature (K)', fontsize=11)
ax4.set_ylabel('Cumulative Probability', fontsize=11)
ax4.set_title('Cumulative Distribution Function', fontsize=12, fontweight='bold')
ax4.grid(True, linestyle='--', alpha=0.7)

# 添加置信区间标注
ax4.fill_betweenx([0, 0.05], T_p5, T_p95, alpha=0.2, color='orange', label='90% CI')
ax4.legend()

plt.tight_layout()
plt.savefig(f'{output_dir}/exp4_uncertainty_quantification.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n  图4已保存: exp4_uncertainty_quantification.png")

print("\n" + "=" * 60)
print("实验4完成")
print("=" * 60)

6.5 网格收敛性指数（GCI）计算

以下代码演示如何计算网格收敛性指数（GCI），这是ASME V&V 20标准推荐的误差估计方法。

# -*- coding: utf-8 -*-
"""
实验5: 网格收敛性指数（GCI）计算
基于ASME V&V 20标准的误差估计
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
import os

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

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

print("\n" + "=" * 60)
print("实验5: 网格收敛性指数（GCI）计算")
print("=" * 60)

def calculate_gci(phi1, phi2, phi3, h1, h2, h3, p_theory=2.0, Fs=1.25):
    """
    计算网格收敛性指数（GCI）
    
    参数:
    phi1, phi2, phi3: 粗、中、细网格上的数值解
    h1, h2, h3: 网格尺寸（h1 > h2 > h3）
    p_theory: 理论收敛阶数
    Fs: 安全因子
    
    返回:
    p_obs: 观察到的收敛阶数
    phi_extrap: Richardson外推解
    GCI23: 细网格的GCI值
    GCI32: 中网格的GCI值
    """
    r21 = h2 / h1
    r32 = h3 / h2
    
    # 计算观察到的收敛阶数
    epsilon32 = phi3 - phi2
    epsilon21 = phi2 - phi1
    
    if abs(epsilon21) < 1e-14:
        p_obs = p_theory
    else:
        p_obs = np.log(abs(epsilon32 / epsilon21)) / np.log(r21)
    
    # Richardson外推
    phi_extrap = (r21**p_obs * phi2 - phi1) / (r21**p_obs - 1)
    
    # GCI计算
    GCI21 = Fs * abs(epsilon21) / (r21**p_obs - 1)
    GCI32 = Fs * abs(epsilon32) / (r32**p_obs - 1)
    
    return p_obs, phi_extrap, GCI21, GCI32

print("\n[5] 网格收敛性指数分析")
print("-" * 50)

# 模拟燃烧仿真数据（火焰温度峰值）
# 使用制造解模拟不同网格上的结果

def manufactured_flame_temp(h, true_value=1800):
    """制造解：模拟网格相关的数值误差"""
    # 模拟二阶精度的数值误差
    error = 50 * h**2 + np.random.normal(0, 5)
    return true_value + error

# 三套网格
grids = {
    'coarse': {'N': 40, 'h': 1/40},
    'medium': {'N': 80, 'h': 1/80},
    'fine': {'N': 160, 'h': 1/160}
}

# 模拟数值结果（多个工况）
n_cases = 5
case_names = ['Case A', 'Case B', 'Case C', 'Case D', 'Case E']

results = {name: {'coarse': [], 'medium': [], 'fine': []} for name in case_names}

np.random.seed(456)
for i, case in enumerate(case_names):
    true_val = 1700 + i * 50  # 不同工况的真实值
    for grid_name in ['coarse', 'medium', 'fine']:
        h = grids[grid_name]['h']
        results[case][grid_name] = manufactured_flame_temp(h, true_val)

# 计算GCI
print("\n  GCI计算结果:")
print("  " + "=" * 70)
print(f"  {'Case':<10} {'p_obs':<8} {'phi_fine':<12} {'phi_extrap':<12} {'GCI32 (%)':<12} {'Status':<10}")
print("  " + "-" * 70)

gci_results = {}
for case in case_names:
    phi1 = results[case]['coarse']
    phi2 = results[case]['medium']
    phi3 = results[case]['fine']
    h1 = grids['coarse']['h']
    h2 = grids['medium']['h']
    h3 = grids['fine']['h']
    
    p_obs, phi_extrap, GCI21, GCI32 = calculate_gci(phi1, phi2, phi3, h1, h2, h3, p_theory=2.0)
    
    # 检查收敛性
    ratio = GCI21 / (GCI32 * (h1/h2)**p_obs)
    status = "Converged" if 0.9 < ratio < 1.1 else "Check"
    
    gci_results[case] = {
        'p_obs': p_obs,
        'phi_fine': phi3,
        'phi_extrap': phi_extrap,
        'GCI32': GCI32,
        'status': status
    }
    
    print(f"  {case:<10} {p_obs:<8.3f} {phi3:<12.2f} {phi_extrap:<12.2f} {GCI32/phi3*100:<12.3f} {status:<10}")

# 绘制GCI分析结果
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. 网格收敛曲线
ax1 = axes[0, 0]
x_pos = np.arange(len(case_names))
width = 0.25

for i, grid_name in enumerate(['coarse', 'medium', 'fine']):
    values = [results[case][grid_name] for case in case_names]
    ax1.bar(x_pos + i*width, values, width, label=f'{grid_name.capitalize()} (h={grids[grid_name]["h"]:.4f})')

ax1.set_xlabel('Test Case', fontsize=11)
ax1.set_ylabel('Temperature (K)', fontsize=11)
ax1.set_title('Grid Convergence for Multiple Cases', fontsize=12, fontweight='bold')
ax1.set_xticks(x_pos + width)
ax1.set_xticklabels(case_names)
ax1.legend()
ax1.grid(True, axis='y', linestyle='--', alpha=0.7)

# 2. 观察到的收敛阶数
ax2 = axes[0, 1]
p_values = [gci_results[case]['p_obs'] for case in case_names]
colors = ['green' if 1.8 < p < 2.2 else 'orange' for p in p_values]

bars = ax2.bar(case_names, p_values, color=colors, edgecolor='black', alpha=0.7)
ax2.axhline(y=2.0, color='r', linestyle='--', linewidth=2, label='Theoretical p=2.0')
ax2.set_xlabel('Test Case', fontsize=11)
ax2.set_ylabel('Observed Order of Accuracy', fontsize=11)
ax2.set_title('Observed Order of Convergence', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, axis='y', linestyle='--', alpha=0.7)

# 添加数值标签
for bar, p in zip(bars, p_values):
    ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.05, 
             f'{p:.2f}', ha='center', va='bottom', fontsize=9)

# 3. GCI值
ax3 = axes[1, 0]
gci_values = [gci_results[case]['GCI32'] for case in case_names]
gci_percent = [gci / gci_results[case]['phi_fine'] * 100 for case, gci in zip(case_names, gci_values)]

bars = ax3.bar(case_names, gci_percent, color='steelblue', edgecolor='black', alpha=0.7)
ax3.axhline(y=5.0, color='r', linestyle='--', linewidth=1.5, label='5% Threshold')
ax3.set_xlabel('Test Case', fontsize=11)
ax3.set_ylabel('GCI (%)', fontsize=11)
ax3.set_title('Grid Convergence Index (Fine Grid)', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, axis='y', linestyle='--', alpha=0.7)

# 添加数值标签
for bar, gci in zip(bars, gci_percent):
    ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.1, 
             f'{gci:.2f}%', ha='center', va='bottom', fontsize=9)

# 4. Richardson外推与数值解对比
ax4 = axes[1, 1]
x_pos = np.arange(len(case_names))
width = 0.35

phi_fine = [gci_results[case]['phi_fine'] for case in case_names]
phi_extrap = [gci_results[case]['phi_extrap'] for case in case_names]

ax4.bar(x_pos - width/2, phi_fine, width, label='Fine Grid Solution', color='skyblue', edgecolor='black')
ax4.bar(x_pos + width/2, phi_extrap, width, label='Richardson Extrapolation', color='lightcoral', edgecolor='black')

ax4.set_xlabel('Test Case', fontsize=11)
ax4.set_ylabel('Temperature (K)', fontsize=11)
ax4.set_title('Fine Grid vs Extrapolated Solution', fontsize=12, fontweight='bold')
ax4.set_xticks(x_pos)
ax4.set_xticklabels(case_names)
ax4.legend()
ax4.grid(True, axis='y', linestyle='--', alpha=0.7)

plt.tight_layout()
plt.savefig(f'{output_dir}/exp5_gci_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n  图5已保存: exp5_gci_analysis.png")

print("\n" + "=" * 60)
print("实验5完成")
print("=" * 60)

6.6 综合V&V流程演示

以下代码演示完整的燃烧模型验证与确认流程。

# -*- coding: utf-8 -*-
"""
实验6: 综合V&V流程演示
完整的燃烧模型验证与确认工作流
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import FancyBboxPatch, FancyArrowPatch
import warnings
warnings.filterwarnings('ignore')
import os

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

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

print("\n" + "=" * 60)
print("实验6: 综合V&V流程演示")
print("=" * 60)

print("\n[6] 完整V&V流程")
print("-" * 50)

# 模拟完整的V&V流程数据
np.random.seed(789)

# 阶段1: 代码验证 - 收敛阶数测试
print("\n  阶段1: 代码验证 (Code Verification)")
print("  " + "-" * 45)

# 模拟不同网格上的误差
N_grids = np.array([20, 40, 80, 160, 320])
h_grids = 1.0 / N_grids
# 二阶格式的理论误差
errors_theory = 0.5 * h_grids**2
# 实际数值误差（带随机扰动）
errors_actual = errors_theory * (1 + np.random.normal(0, 0.1, len(N_grids)))

# 计算观察到的收敛阶数
p_observed = np.polyfit(np.log(h_grids), np.log(errors_actual), 1)[0]
print(f"    理论收敛阶数: p = 2.0")
print(f"    观察收敛阶数: p = {p_observed:.3f}")
print(f"    代码验证状态: {'通过' if 1.8 < p_observed < 2.2 else '需检查'}")

# 阶段2: 计算验证 - 网格收敛性
print("\n  阶段2: 计算验证 (Calculation Verification)")
print("  " + "-" * 45)

# 三套网格的模拟结果
phi_coarse = 1750.5
phi_medium = 1785.3
phi_fine = 1795.8

# Richardson外推
phi_extrap = (4 * phi_fine - phi_medium) / 3
# GCI计算
GCI_fine = 1.25 * abs(phi_fine - phi_medium) / (2**2 - 1)

print(f"    粗网格解: {phi_coarse:.2f} K")
print(f"    中网格解: {phi_medium:.2f} K")
print(f"    细网格解: {phi_fine:.2f} K")
print(f"    Richardson外推: {phi_extrap:.2f} K")
print(f"    GCI (细网格): ±{GCI_fine:.2f} K ({GCI_fine/phi_fine*100:.2f}%)")

# 阶段3: 模型确认 - 实验对比
print("\n  阶段3: 模型确认 (Validation)")
print("  " + "-" * 45)

# 实验数据
T_exp_mean = 1800.0
T_exp_std = 15.0
# 数值结果（考虑数值误差）
T_sim = phi_fine
T_sim_uncertainty = GCI_fine

# 确认度量
d_val = abs(T_sim - T_exp_mean) / np.sqrt(T_sim_uncertainty**2 + T_exp_std**2)
print(f"    实验测量值: {T_exp_mean:.2f} ± {T_exp_std:.2f} K")
print(f"    数值预测值: {T_sim:.2f} ± {T_sim_uncertainty:.2f} K")
print(f"    确认度量 d: {d_val:.3f}")
print(f"    确认状态: {'通过' if d_val < 1.0 else '未通过'} (d < 1.0)")

# 阶段4: 不确定性量化
print("\n  阶段4: 不确定性量化 (Uncertainty Quantification)")
print("  " + "-" * 45)

# 输入参数不确定性
param_uncertainties = {
    'Reaction Rate': 0.08,
    'Heat Transfer': 0.05,
    'Turbulence Model': 0.12,
    'Boundary Conditions': 0.06,
    'Numerical Error': GCI_fine/phi_fine
}

# 总不确定性（平方和开根）
total_uncertainty = np.sqrt(sum([u**2 for u in param_uncertainties.values()]))
print(f"    参数不确定性贡献:")
for param, unc in param_uncertainties.items():
    print(f"      {param:<20s}: {unc*100:.2f}%")
print(f"    总预测不确定度: ±{total_uncertainty*100:.2f}%")

# 阶段5: 敏感性分析
print("\n  阶段5: 敏感性分析 (Sensitivity Analysis)")
print("  " + "-" * 45)

sensitivities = {
    'Activation Energy': 0.45,
    'Pre-exponential Factor': 0.32,
    'Heat Loss Coefficient': 0.28,
    'Inlet Temperature': 0.15,
    'Pressure': 0.08
}

sorted_sens = sorted(sensitivities.items(), key=lambda x: x[1], reverse=True)
print(f"    参数敏感性排序:")
for i, (param, sens) in enumerate(sorted_sens, 1):
    bar = '█' * int(sens * 20)
    print(f"      {i}. {param:<25s}: {sens:.3f} {bar}")

# 综合评估
print("\n  综合V&V评估")
print("  " + "=" * 45)

scores = {
    'Code Verification': 0.95 if 1.8 < p_observed < 2.2 else 0.6,
    'Calculation Verification': 0.90 if GCI_fine/phi_fine < 0.05 else 0.7,
    'Validation': 0.85 if d_val < 1.0 else 0.5,
    'Uncertainty Quantification': 0.88,
    'Sensitivity Analysis': 0.92
}

overall_score = np.mean(list(scores.values()))
print(f"    各阶段评分:")
for stage, score in scores.items():
    status = 'Excellent' if score > 0.9 else 'Good' if score > 0.75 else 'Fair'
    print(f"      {stage:<30s}: {score:.2f} ({status})")

print(f"\n    综合评分: {overall_score:.2f}")
print(f"    模型可信度等级: {'A级(高可信度)' if overall_score > 0.9 else 'B级(中等可信度)' if overall_score > 0.75 else 'C级(需改进)'}")

# 绘制V&V流程图
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. V&V流程各阶段评分
ax1 = axes[0, 0]
stages = list(scores.keys())
score_values = list(scores.values())
colors = ['#2ecc71' if s > 0.9 else '#f39c12' if s > 0.75 else '#e74c3c' for s in score_values]

bars = ax1.barh(stages, score_values, color=colors, edgecolor='black', alpha=0.8)
ax1.axvline(x=0.9, color='green', linestyle='--', linewidth=2, label='Excellent (0.9)')
ax1.axvline(x=0.75, color='orange', linestyle='--', linewidth=2, label='Good (0.75)')
ax1.set_xlabel('Score', fontsize=11)
ax1.set_xlim(0, 1)
ax1.set_title('V&V Process Stage Scores', fontsize=12, fontweight='bold')
ax1.legend(loc='lower right')
ax1.grid(True, axis='x', linestyle='--', alpha=0.7)

# 添加数值标签
for bar, score in zip(bars, score_values):
    ax1.text(score + 0.02, bar.get_y() + bar.get_height()/2, f'{score:.2f}', 
             va='center', fontsize=10, fontweight='bold')

# 2. 收敛阶数验证
ax2 = axes[0, 1]
ax2.loglog(h_grids, errors_actual, 'bo-', linewidth=2, markersize=8, label='Observed Error')
ax2.loglog(h_grids, errors_theory, 'r--', linewidth=2, label='Theoretical (p=2.0)')

# 拟合线
p_fit = np.polyfit(np.log(h_grids), np.log(errors_actual), 1)
h_fit = np.linspace(h_grids.min(), h_grids.max(), 100)
error_fit = np.exp(np.polyval(p_fit, np.log(h_fit)))
ax2.loglog(h_fit, error_fit, 'g:', linewidth=2, label=f'Fit (p={p_observed:.2f})')

ax2.set_xlabel('Grid Spacing h', fontsize=11)
ax2.set_ylabel('Error', fontsize=11)
ax2.set_title('Code Verification: Order of Accuracy', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, which='both', linestyle='--', alpha=0.7)

# 3. 实验对比与确认度量
ax3 = axes[1, 0]
x_pos = np.arange(1)
width = 0.25

ax3.bar(x_pos - width, [T_exp_mean], width, label='Experiment', color='lightblue', 
        edgecolor='black', yerr=[T_exp_std], capsize=5)
ax3.bar(x_pos, [T_sim], width, label='Simulation', color='lightcoral', 
        edgecolor='black', yerr=[T_sim_uncertainty], capsize=5)
ax3.bar(x_pos + width, [phi_extrap], width, label='Extrapolated', color='lightgreen', 
        edgecolor='black', alpha=0.7)

ax3.set_ylabel('Temperature (K)', fontsize=11)
ax3.set_title(f'Validation: d={d_val:.3f} ({"PASS" if d_val < 1.0 else "FAIL"})', 
              fontsize=12, fontweight='bold')
ax3.set_xticks(x_pos)
ax3.set_xticklabels(['Flame Temperature'])
ax3.legend()
ax3.grid(True, axis='y', linestyle='--', alpha=0.7)

# 添加确认度量文本
ax3.text(0.5, max(T_exp_mean, T_sim) + 30, f'Validation Metric d={d_val:.3f}\n({"PASS" if d_val < 1.0 else "FAIL"})', 
         ha='center', fontsize=11, bbox=dict(boxstyle='round', facecolor='lightgreen' if d_val < 1.0 else 'lightyellow', alpha=0.8))

# 4. 不确定性分解
ax4 = axes[1, 1]
param_names = list(param_uncertainties.keys())
unc_values = [v * 100 for v in param_uncertainties.values()]  # 转换为百分比
colors_pie = plt.cm.Set3(np.linspace(0, 1, len(param_names)))

wedges, texts, autotexts = ax4.pie(unc_values, labels=param_names, autopct='%1.1f%%', 
                                    colors=colors_pie, startangle=90)
for autotext in autotexts:
    autotext.set_color('white')
    autotext.set_fontweight('bold')
    autotext.set_fontsize(9)

ax4.set_title(f'Uncertainty Contribution\nTotal: ±{total_uncertainty*100:.1f}%', 
              fontsize=12, fontweight='bold')

plt.tight_layout()
plt.savefig(f'{output_dir}/exp6_vv_process.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n  图6已保存: exp6_vv_process.png")

print("\n" + "=" * 60)
print("实验6完成")
print("=" * 60)

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