辐射换热仿真-主题038_辐射换分的验证与确认-主题038_辐射换分的验证与确认

kkchenjj

76人浏览 · 2026-04-26 22:29:14

kkchenjj · 2026-04-26 22:29:14 发布

主题038：辐射换分的验证与确认

摘要

验证与确认（Verification and Validation, V&V）是确保辐射换热仿真结果可信度的核心环节。本主题系统介绍V&V的理论框架、实施方法和工程实践。从代码验证、解验证到模型确认三个层次，详细阐述如何通过基准问题、网格收敛性分析、Richardson外推、实验对比等手段建立对仿真模型的信心。本主题还涵盖不确定性量化、敏感性分析、误差估计等关键技术，并通过典型案例展示完整的V&V流程在实际工程问题中的应用。

关键词

验证与确认、V&V、基准问题、网格收敛、Richardson外推、误差分析、不确定性量化、模型可信度

在这里插入图片描述

1. 引言

1.1 为什么需要验证与确认

在工程设计和科学研究的决策过程中，计算仿真正发挥着越来越重要的作用。对于辐射换热这类复杂的物理过程，数值模拟可以预测极端条件下的热行为、优化系统性能、减少试验成本。然而，仿真结果的可信度始终是使用者最关心的问题——我们如何确信计算机生成的结果能够准确反映物理现实？

验证与确认（V&V）正是回答这一问题的系统方法论。美国航空航天学会（AIAA）将验证定义为"确定计算模型是否准确表示基础数学模型及其解的过程"，将确认定义为"确定计算模型在其预期用途范围内准确表示物理现实程度的过程"。简而言之：

验证（Verification）：我们是否正确求解了方程？
确认（Validation）：我们是否正确描述了物理现象？

1.2 V&V的历史发展

V&V作为一门系统学科，其发展历程与计算流体力学和计算固体力学的成熟密切相关。20世纪90年代，美国能源部和国防部开始资助大规模V&V研究项目，建立了最初的理论框架。2000年以后，美国机械工程师学会（ASME）发布了V&V标准（V&V 10、V&V 20），为工程实践提供了规范指导。

在辐射换热领域，V&V面临独特挑战：辐射传递方程的积分-微分特性、光谱相关性、非线性边界条件等都增加了验证的复杂性。国际传热传质中心（ICHMT）多次组织辐射换热基准问题研讨会，建立了系列标准测试案例。

1.3 V&V的完整流程

一个完整的V&V流程包含以下环节：

规划阶段：

明确仿真目标和精度要求
识别关键物理过程和模型假设
制定V&V计划，分配资源

代码验证阶段：

选择适当的基准问题
进行网格收敛性研究
验证数值方法的正确实现

解验证阶段：

估计数值误差
分析迭代收敛性
评估计算不确定性

模型确认阶段：

设计对比实验
量化实验不确定性
进行模型-实验对比

文档与报告：

记录所有假设和简化
报告误差估计和可信度水平
明确模型的适用范围和局限性

2. 代码验证

2.1 代码验证的基本概念

代码验证的目标是确保计算程序正确实现了所声称的数学模型。这包括两个层面：

软件质量保证（SQA）：确保代码无编程错误、逻辑正确、文档完整。包括代码审查、单元测试、版本控制等软件工程实践。

数值算法验证：验证离散格式、求解器、边界条件处理等数值方法的正确性。

2.2 基准问题方法

基准问题是代码验证的核心工具。理想的基准问题应具有以下特征：

解析解存在：可以与精确解对比，消除参考解的不确定性。

物理代表性：虽然简化，但应包含目标问题的主要物理特征。

数值挑战性：能够测试代码处理困难情况的能力（如强非线性、间断等）。

文献充分：被广泛研究，有大量独立验证结果可供对比。

2.3 辐射换热的经典基准问题

一维灰体介质辐射平衡：
考虑光学厚度为 $τL\tau_L$ 的一维灰体介质，在辐射平衡条件下求解温度分布。解析解由以下隐式方程给出：

$θ4(τ)=12[E2(τ)+E2(τL−τ)]+12∫0τLθ4(τ′)E1(∣τ−τ′∣)dτ′\theta^4(\tau) = \frac{1}{2}[E_2(\tau) + E_2(\tau_L - \tau)] + \frac{1}{2}\int_0^{\tau_L} \theta^4(\tau')E_1(|\tau - \tau'|)d\tau'$

其中 $θ=T/Tref\theta = T/T_{ref}$ 是无量纲温度， $E_n$ 是指数积分函数。

二维空腔辐射：
考虑具有漫反射壁面的矩形空腔，验证视角因子计算和辐射交换的正确性。可以使用交叉弦法计算解析解。

Marshak波问题：
描述辐射热波在非平衡介质中的传播，测试瞬态辐射传递方程求解器。

Hottel球问题：
验证参与性介质中辐射传递的蒙特卡洛实现，通过比较数值结果与Hottel的经典解。

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

MMS是一种强大的代码验证技术，其基本思想是：

选择一个任意的解析函数作为"制造解"
将该解代入控制方程，计算源项
用修改后的方程（包含源项）测试代码
验证数值解收敛到制造解

MMS的优势在于可以构造任意复杂的测试案例，全面测试代码的各个部分。对于辐射传递方程：

$dIds=−βI+κIb+σs4π∫4πI(Ω′)Φ(Ω′,Ω)dΩ′+QMMS\frac{dI}{ds} = -\beta I + \kappa I_b + \frac{\sigma_s}{4\pi}\int_{4\pi} I(\Omega')\Phi(\Omega', \Omega)d\Omega' + Q_{MMS}$

其中 $Q_{MMS}$ 是使制造解满足方程的源项。

2.5 收敛性分析

验证数值方法是否正确实现的关键是检查收敛阶数。对于具有 $p$ 阶精度的方法，误差应满足：

$ϵ∝hp\epsilon \propto h^p$

其中 $h$ 是网格尺寸。通过对数坐标绘制误差-网格尺寸关系，斜率应等于理论收敛阶数。

对于辐射换热问题，收敛性分析需要考虑：

空间离散误差（角度和空间）
迭代收敛误差
舍入误差（在双精度计算中通常可忽略）

3. 解验证

3.1 数值误差的来源

解验证关注特定计算中的数值误差估计。辐射换热仿真的数值误差主要来源包括：

离散化误差：由于将连续问题离散为有限自由度系统而产生的误差。这是最主要的误差来源。

迭代误差：当使用迭代求解器时，由于提前终止迭代而产生的误差。

舍入误差：有限精度算术运算引入的误差，在现代计算中通常很小。

模型形式误差：由于简化物理模型（如灰体假设、各向同性散射）引入的误差。

3.2 网格收敛性研究

网格收敛性研究是估计离散化误差的标准方法。基本流程是：

在三种或更多不同密度的网格上进行计算
监测关键量的变化
使用Richardson外推估计网格无关解
计算网格收敛指标

网格细化因子：通常选择 $\sqrt{2}$ 或 $r = 2$ ，即每维网格数增加 $2\sqrt{2}$ 或 $2$ 倍。

收敛性判据：

$∣fmedium−fcoarse∣∣ffine−fmedium∣≈rp\frac{|f_{medium} - f_{coarse}|}{|f_{fine} - f_{medium}|} \approx r^p$

其中 $p$ 是观察到的收敛阶数，应与理论值接近。

3.3 Richardson外推

Richardson外推是估计网格无关解和离散化误差的有力工具。假设解在网格上的渐近展开为：

$f_h = f_{exact} + Ch^p + O(h^{p+1})$

使用两个网格（细网格 $h$ 和粗网格 $r h$ ）的结果，可以估计精确解：

$fexact≈fh,extrap=rpfh−frhrp−1f_{exact} \approx f_{h,extrap} = \frac{r^p f_h - f_{rh}}{r^p - 1}$

离散化误差估计为：

$ϵh≈∣fh,extrap−fh∣\epsilon_h \approx |f_{h,extrap} - f_h|$

对于辐射换热问题，Richardson外推可以分别应用于空间离散和角度离散。

3.4 网格收敛指标（GCI）

网格收敛指标（Grid Convergence Index）是Roache提出的标准化误差估计方法，提供了近似95%置信水平的误差带：

$\frac{F_s |\epsilon|}{r^p - 1}$

其中 $F_s$ 是安全系数（通常取1.25或3）， $ϵ=(ffine−fcoarse)/ffine\epsilon = (f_{fine} - f_{coarse})/f_{fine}$ 是近似相对误差。

GCI的优势在于：

提供一致的误差报告格式
考虑收敛阶数的不确定性
便于不同研究者的结果对比

3.5 角度离散误差

辐射传递方程的角度离散（离散坐标法、球谐函数法等）引入额外的误差源。角度收敛性研究应：

使用不同阶数的离散坐标法（ $S_2, S_4, S_8, S_{16}$ 等）
或使用不同阶数的球谐函数展开（ $P_1, P_3, P_5$ 等）
监测辐射热流和入射辐射的收敛

对于强各向异性问题，可能需要自适应角度细化或高阶方法。

4. 模型确认

4.1 确认实验的设计原则

模型确认需要与精心设计的实验进行对比。确认实验应遵循以下原则：

层次化方法：从简单到复杂，逐步增加物理复杂性。先确认基本模型，再确认耦合现象。

充分仪表化：测量所有相关物理量，包括边界条件、材料属性、响应量等。

不确定性量化：系统识别和量化所有实验不确定性来源。

盲测验证：在可能的情况下，进行预测性验证（先预测后实验）。

4.2 确认度量

确认度量用于量化模型预测与实验观测之间的一致性。常用度量包括：

误差范数：

$L2=1N∑i=1N(yipred−yiexp)2L_2 = \sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_i^{pred} - y_i^{exp})^2}$

$L∞=max⁡i∣yipred−yiexp∣L_\infty = \max_i |y_i^{pred} - y_i^{exp}|$

决定系数：

$R2=1−∑(yipred−yiexp)2∑(yiexp−yˉexp)2R^2 = 1 - \frac{\sum(y_i^{pred} - y_i^{exp})^2}{\sum(y_i^{exp} - \bar{y}^{exp})^2}$

确认指标（Validation Metric）：

考虑实验不确定性的确认度量：

$\frac{|y^{pred} - y^{exp}|}{\sqrt{u_{num}^2 + u_{exp}^2}}$

其中 $u_{num}$ 和 $u_{exp}$ 分别是数值和实验不确定性的标准差。

4.3 确认域与适用域

确认域（Validation Domain）：已经通过实验确认的参数空间范围。

适用域（Domain of Applicability）：模型可以可靠预测的参数空间范围，通常比确认域大，需要通过外推分析确定。

建立适用域需要：

识别关键物理过程和参数
进行参数敏感性分析
评估模型在确认域边界的行为
必要时进行额外实验扩展确认域

4.4 模型形式不确定性

模型形式不确定性（Model Form Uncertainty）是由于物理模型简化而产生的固有不确定性。这种不确定性无法通过网格细化或实验改进消除。

量化模型形式不确定性的方法：

多模型对比：使用不同复杂程度的模型求解同一问题，比较结果差异。

贝叶斯模型平均：基于数据证据对多个模型进行加权平均。

模型可信度评估：由领域专家评估模型假设的合理性。

5. 不确定性量化

5.1 不确定性的来源与分类

辐射换热仿真中的不确定性可分为：

偶然不确定性（Aleatory Uncertainty）：固有的随机性，如材料属性的自然波动、边界条件的随机变化。可用概率分布描述。

认知不确定性（Epistemic Uncertainty）：由于缺乏知识而产生的不确定性，如模型形式、参数值的不确定。可用区间或可能性理论描述。

数值不确定性：由于离散化、迭代截断等数值近似引入的不确定性。

5.2 敏感性分析方法

敏感性分析识别对输出影响最大的输入参数，指导实验设计和模型改进。

局部分析法：

$Si=∂f∂xixifS_i = \frac{\partial f}{\partial x_i}\frac{x_i}{f}$

计算在名义值附近的局部敏感性。

全局敏感性分析：

Sobol指数基于方差分解：

$Si=Varxi(Ex−i[f∣xi])Var(f)S_i = \frac{\text{Var}_{x_i}(E_{x_{-i}}[f|x_i])}{\text{Var}(f)}$

考虑参数在整个变化范围内的影响，包括交互效应。

** Morris筛选法**：

计算基本效应的统计分布，高效识别重要参数。

5.3 传播方法

不确定性传播方法将输入不确定性传递到输出：

蒙特卡洛方法：

从输入分布中随机采样，运行大量计算，统计分析输出分布。简单但计算成本高。

多项式混沌展开：

将随机输出展开为随机输入的正交多项式：

$f(ξ)=∑k=0PckΨk(ξ)f(\xi) = \sum_{k=0}^{P} c_k \Psi_k(\xi)$

其中 $Ψk\Psi_k$ 是多项式混沌基函数。通过少量计算即可估计输出统计量。

随机配点法：

在精心选择的配点上求解确定性问题，通过插值获得输出分布。

5.4 实验不确定性分析

实验不确定性分析遵循ISO/IEC Guide 98-3（GUM）标准：

A类评定：基于统计分析的重复观测不确定性。

B类评定：基于其他信息（校准证书、制造商规格、经验）的不确定性。

合成标准不确定度：

$uc=∑i=1N(∂f∂xi)2u2(xi)+2∑i=1N−1∑j=i+1N∂f∂xi∂f∂xju(xi,xj)u_c = \sqrt{\sum_{i=1}^{N}\left(\frac{\partial f}{\partial x_i}\right)^2 u^2(x_i) + 2\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}\frac{\partial f}{\partial x_i}\frac{\partial f}{\partial x_j}u(x_i, x_j)}$

扩展不确定度： $U = ku_c$ ，其中 $k = 2$ 对应约95%置信水平。

6. 误差估计与自适应方法

6.1 后验误差估计

后验误差估计基于计算结果估计误差，用于指导自适应细化。

残差法：

计算控制方程的残差：

$\mathcal{L}(u_h) - f$

残差大小与误差相关，可通过求解局部问题估计误差分布。

恢复法：

通过超收敛恢复技术（如SPR、ZZ恢复）获得更精确解的估计，与计算解的差异作为误差指示。

对偶加权残差法（DWR）：

求解对偶问题，将残差与对偶解加权，获得对特定输出量（如热流）的误差估计。

6.2 自适应网格细化

自适应方法根据误差指示局部调整网格密度：

h-细化：在误差大的区域细分单元。

p-细化：在误差大的区域提高多项式阶数。

r-细化：重新分布节点位置。

自适应流程：

在初始网格上求解
计算误差指示
标记需要细化的单元
生成新网格
重复直到满足精度要求

对于辐射换热，自适应需要考虑空间和角度两个维度的误差平衡。

6.3 目标导向的自适应

目标导向自适应关注特定输出量（如某点的温度或热流）的精度，而非整个场的精度。

通过求解对偶问题，可以获得输出量对局部误差的敏感性：

$δJ≈∫R(uh)z∗dΩ\delta J \approx \int R(u_h) z^* d\Omega$

其中 $z^*$ 是对偶解， $R$ 是残差。这允许在最影响目标输出的区域集中计算资源。

"""
主题038：辐射换分的验证与确认

本程序演示辐射换热仿真中的验证与确认方法，包括：
1. 代码验证：基准问题求解与收敛性分析
2. 解验证：网格收敛性研究、Richardson外推、GCI计算
3. 模型确认：实验对比、误差分析、确认度量
4. 不确定性量化：敏感性分析、蒙特卡洛传播
5. 制造解方法（MMS）验证
6. GIF动画展示V&V过程

作者：仿真专家
日期：2026-03-04
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from scipy.integrate import quad, solve_ivp
from scipy.special import expn
from scipy.interpolate import interp1d
import warnings
warnings.filterwarnings('ignore')

# 设置matplotlib后端为Agg，不弹出窗口
plt.switch_backend('Agg')

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

print("=" * 70)
print("主题038：辐射换分的验证与确认")
print("=" * 70)

# =============================================================================
# 案例1：代码验证 - 一维灰体介质辐射平衡基准问题（P1近似）
# =============================================================================
print("\n" + "=" * 70)
print("案例1：代码验证 - 一维灰体介质辐射平衡（P1近似）")
print("=" * 70)

def solve_radiative_equilibrium_p1(tau_L, N=100, max_iter=100, tol=1e-6):
    """
    使用P1近似求解一维灰体介质辐射平衡问题
    
    P1近似将辐射传递方程简化为扩散形式，计算效率高。
    
    参数:
        tau_L: 光学厚度
        N: 网格数
        max_iter: 最大迭代次数
        tol: 收敛容差
    
    返回:
        tau: 光学厚度坐标
        theta: 无量纲温度分布
        n_iter: 迭代次数
    """
    # 光学厚度网格
    tau = np.linspace(0, tau_L, N)
    d_tau = tau[1] - tau[0]
    
    # 初始猜测：线性温度分布
    theta = 0.5 + 0.5 * tau / tau_L if tau_L > 0 else np.ones(N) * 0.5
    
    # P1近似：辐射热流与温度梯度关系 q_r = -4/3 * d(theta^4)/d(tau)
    # 辐射平衡：d(q_r)/d(tau) = 0 => d²(theta^4)/d(tau)² = 0
    # 边界条件使用Marshak边界条件
    
    # 构建三对角矩阵求解 theta^4
    theta4 = theta**4
    
    for iteration in range(max_iter):
        theta4_old = theta4.copy()
        
        # 使用有限差分求解 d²(theta^4)/d(tau)² = 0
        # 内部节点使用中心差分
        for i in range(1, N-1):
            theta4[i] = 0.5 * (theta4[i-1] + theta4[i+1])
        
        # Marshak边界条件
        # 左边界: theta^4(0) - 2/3 * d(theta^4)/d(tau)|_0 = 1 (单位入射)
        theta4[0] = (4 * theta4[1] - theta4[2] + 3/d_tau) / (3 + 3/d_tau)
        
        # 右边界: theta^4(tau_L) + 2/3 * d(theta^4)/d(tau)|_L = 1
        theta4[N-1] = (4 * theta4[N-2] - theta4[N-3] + 3/d_tau) / (3 + 3/d_tau)
        
        # 检查收敛
        error = np.max(np.abs(theta4 - theta4_old))
        if error < tol:
            break
    
    theta = theta4**0.25
    return tau, theta, iteration + 1

def case1_code_verification():
    """
    案例1：代码验证 - 基准问题求解与收敛性分析
    """
    print("\n案例1：一维灰体介质辐射平衡代码验证（P1近似）")
    print("-" * 50)
    
    # 不同光学厚度下的解
    tau_L_values = [0.1, 1.0, 10.0]
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 4))
    
    for idx, tau_L in enumerate(tau_L_values):
        print(f"\n光学厚度 τ_L = {tau_L}")
        
        # 不同网格密度
        N_values = [20, 40, 80, 160]
        colors = ['blue', 'green', 'orange', 'red']
        
        ax = axes[idx]
        
        for N, color in zip(N_values, colors):
            tau, theta, n_iter = solve_radiative_equilibrium_p1(tau_L, N=N)
            ax.plot(tau/tau_L, theta, color=color, linewidth=2, 
                   label=f'N={N} ({n_iter} iter)')
            print(f"  N={N:3d}: {n_iter}次迭代收敛")
        
        ax.set_xlabel(r'Dimensionless Position $\tau/\tau_L$')
        ax.set_ylabel(r'Dimensionless Temperature $\theta$')
        ax.set_title(f'$\\tau_L$ = {tau_L} (P1 Approximation)')
        ax.legend()
        ax.grid(True, alpha=0.3)
        ax.set_xlim(0, 1)
    
    plt.tight_layout()
    plt.savefig('case1_code_verification.png', dpi=150, bbox_inches='tight')
    print("\n案例1图表已保存: case1_code_verification.png")
    plt.close()
    
    return tau_L_values

# 运行案例1
tau_L_values = case1_code_verification()

# =============================================================================
# 案例2：解验证 - 网格收敛性研究与Richardson外推
# =============================================================================
print("\n" + "=" * 70)
print("案例2：解验证 - 网格收敛性与Richardson外推")
print("=" * 70)

def solve_diffusion_radiation_1d(N, kappa, sigma_s, q_gen, T_left, T_right, L):
    """
    求解一维辐射-导热耦合问题（扩散近似）
    
    参数:
        N: 网格数
        kappa: 吸收系数 [1/m]
        sigma_s: 散射系数 [1/m]
        q_gen: 体积热源 [W/m³]
        T_left, T_right: 边界温度 [K]
        L: 域长度 [m]
    
    返回:
        x: 位置坐标
        T: 温度分布 [K]
        q: 热流分布 [W/m²]
    """
    # 网格生成
    x = np.linspace(0, L, N)
    dx = x[1] - x[0]
    
    # 材料属性
    k = 1.0  # 导热系数 [W/(m·K)]
    sigma = 5.67e-8  # Stefan-Boltzmann常数
    
    # 辐射导热系数（Rosseland近似）
    beta = kappa + sigma_s  # 消光系数
    k_rad = 16 * sigma * (1000**3) / (3 * beta)  # 使用参考温度
    k_eff = k + k_rad
    
    # 构建三对角矩阵
    a = np.zeros(N)  # 下对角
    b = np.zeros(N)  # 主对角
    c = np.zeros(N)  # 上对角
    d = np.zeros(N)  # 右端项
    
    # 内部节点
    for i in range(1, N-1):
        a[i] = k_eff / dx**2
        b[i] = -2 * k_eff / dx**2
        c[i] = k_eff / dx**2
        d[i] = -q_gen
    
    # 边界条件
    b[0] = 1.0
    d[0] = T_left
    
    b[N-1] = 1.0
    d[N-1] = T_right
    
    # Thomas算法求解
    # 前向消去
    for i in range(1, N):
        m = a[i] / b[i-1]
        b[i] = b[i] - m * c[i-1]
        d[i] = d[i] - m * d[i-1]
    
    # 回代
    T = np.zeros(N)
    T[N-1] = d[N-1] / b[N-1]
    for i in range(N-2, -1, -1):
        T[i] = (d[i] - c[i] * T[i+1]) / b[i]
    
    # 计算热流
    q = -k_eff * np.gradient(T, dx)
    
    return x, T, q

def case2_solution_verification():
    """
    案例2：网格收敛性研究与Richardson外推
    """
    print("\n案例2：网格收敛性研究")
    print("-" * 50)
    
    # 问题参数
    L = 1.0  # 域长度 [m]
    kappa = 10.0  # 吸收系数 [1/m]
    sigma_s = 5.0  # 散射系数 [1/m]
    q_gen = 10000.0  # 体积热源 [W/m³]
    T_left = 300.0  # 左边界温度 [K]
    T_right = 400.0  # 右边界温度 [K]
    
    # 三种网格密度（细化因子 r = 2）
    N_coarse = 20
    N_medium = 40
    N_fine = 80
    
    # 在三种网格上求解
    x_c, T_c, q_c = solve_diffusion_radiation_1d(N_coarse, kappa, sigma_s, 
                                                   q_gen, T_left, T_right, L)
    x_m, T_m, q_m = solve_diffusion_radiation_1d(N_medium, kappa, sigma_s, 
                                                   q_gen, T_left, T_right, L)
    x_f, T_f, q_f = solve_diffusion_radiation_1d(N_fine, kappa, sigma_s, 
                                                   q_gen, T_left, T_right, L)
    
    # 在共同点上采样（中点）
    x_mid = L / 2
    T_c_mid = np.interp(x_mid, x_c, T_c)
    T_m_mid = np.interp(x_mid, x_m, T_m)
    T_f_mid = np.interp(x_mid, x_f, T_f)
    
    print(f"\n中点温度对比:")
    print(f"  粗网格 (N={N_coarse}):  T = {T_c_mid:.4f} K")
    print(f"  中网格 (N={N_medium}):  T = {T_m_mid:.4f} K")
    print(f"  细网格 (N={N_fine}):  T = {T_f_mid:.4f} K")
    
    # Richardson外推
    r = 2.0  # 细化因子
    p = 2.0  # 理论收敛阶数（二阶精度）
    
    # 估计精确解
    T_exact_est = (r**p * T_f_mid - T_m_mid) / (r**p - 1)
    
    # 估计误差
    eps_fine = abs(T_f_mid - T_exact_est)
    eps_medium = abs(T_m_mid - T_exact_est) / r**p
    eps_coarse = abs(T_c_mid - T_exact_est) / r**(2*p)
    
    print(f"\nRichardson外推结果:")
    print(f"  估计精确解: T_exact ≈ {T_exact_est:.6f} K")
    print(f"  细网格误差估计: {eps_fine:.6f} K ({eps_fine/T_exact_est*100:.4f}%)")
    print(f"  中网格误差估计: {eps_medium:.6f} K ({eps_medium/T_exact_est*100:.4f}%)")
    print(f"  粗网格误差估计: {eps_coarse:.6f} K ({eps_coarse/T_exact_est*100:.4f}%)")
    
    # 网格收敛指标 (GCI)
    F_s = 1.25  # 安全系数
    epsilon = abs(T_f_mid - T_m_mid) / T_f_mid
    GCI = F_s * epsilon / (r**p - 1)
    
    print(f"\n网格收敛指标 (GCI):")
    print(f"  近似相对误差: {epsilon*100:.4f}%")
    print(f"  GCI (95%置信水平): {GCI*100:.4f}%")
    
    # 检查收敛阶数
    p_observed = np.log(abs(T_m_mid - T_c_mid) / abs(T_f_mid - T_m_mid)) / np.log(r)
    print(f"\n收敛阶数分析:")
    print(f"  理论收敛阶数: {p}")
    print(f"  观察收敛阶数: {p_observed:.4f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 1. 温度分布对比
    ax1 = axes[0, 0]
    ax1.plot(x_c, T_c, 'b-o', linewidth=2, markersize=6, label=f'Coarse (N={N_coarse})')
    ax1.plot(x_m, T_m, 'g-s', linewidth=2, markersize=6, label=f'Medium (N={N_medium})')
    ax1.plot(x_f, T_f, 'r-^', linewidth=2, markersize=6, label=f'Fine (N={N_fine})')
    ax1.axhline(y=T_exact_est, color='k', linestyle='--', linewidth=2, label=f'Richardson Extrap.')
    ax1.axvline(x=x_mid, color='gray', linestyle=':', alpha=0.5)
    ax1.set_xlabel('Position x [m]')
    ax1.set_ylabel('Temperature [K]')
    ax1.set_title('Temperature Distribution - Grid Convergence')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 热流分布对比
    ax2 = axes[0, 1]
    ax2.plot(x_c, q_c/1000, 'b-o', linewidth=2, markersize=6, label=f'Coarse (N={N_coarse})')
    ax2.plot(x_m, q_m/1000, 'g-s', linewidth=2, markersize=6, label=f'Medium (N={N_medium})')
    ax2.plot(x_f, q_f/1000, 'r-^', linewidth=2, markersize=6, label=f'Fine (N={N_fine})')
    ax2.set_xlabel('Position x [m]')
    ax2.set_ylabel('Heat Flux [kW/m²]')
    ax2.set_title('Heat Flux Distribution - Grid Convergence')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 3. 网格收敛性（对数坐标）
    ax3 = axes[1, 0]
    h_values = np.array([1/N_coarse, 1/N_medium, 1/N_fine])
    T_values = np.array([T_c_mid, T_m_mid, T_f_mid])
    errors = np.abs(T_values - T_exact_est)
    
    ax3.loglog(h_values, errors, 'bo-', linewidth=2, markersize=10)
    # 理论收敛线
    h_theory = np.linspace(h_values[0], h_values[-1], 100)
    error_theory = errors[0] * (h_theory / h_values[0])**p
    ax3.loglog(h_theory, error_theory, 'r--', linewidth=2, label=f'Theoretical O(h^{p:.0f})')
    
    ax3.set_xlabel('Grid Spacing h')
    ax3.set_ylabel('Error |T - T_exact|')
    ax3.set_title(f'Grid Convergence (Observed Order: {p_observed:.2f})')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 4. 误差估计对比
    ax4 = axes[1, 1]
    grids = ['Coarse', 'Medium', 'Fine']
    error_estimates = [eps_coarse, eps_medium, eps_fine]
    colors = ['blue', 'green', 'red']
    
    bars = ax4.bar(grids, [e/T_exact_est*100 for e in error_estimates], color=colors, alpha=0.7)
    ax4.axhline(y=GCI*100, color='black', linestyle='--', linewidth=2, label=f'GCI = {GCI*100:.3f}%')
    ax4.set_ylabel('Relative Error [%]')
    ax4.set_title('Estimated Discretization Error')
    ax4.legend()
    ax4.grid(True, alpha=0.3, axis='y')
    
    # 在柱状图上添加数值标签
    for bar, err in zip(bars, error_estimates):
        height = bar.get_height()
        ax4.text(bar.get_x() + bar.get_width()/2., height,
                f'{err/T_exact_est*100:.4f}%',
                ha='center', va='bottom', fontsize=10)
    
    plt.tight_layout()
    plt.savefig('case2_solution_verification.png', dpi=150, bbox_inches='tight')
    print("\n案例2图表已保存: case2_solution_verification.png")
    plt.close()
    
    return x_f, T_f, T_exact_est, GCI

# 运行案例2
x_f, T_f, T_exact_est, GCI = case2_solution_verification()

# =============================================================================
# 案例3：模型确认 - 仿真与实验对比
# =============================================================================
print("\n" + "=" * 70)
print("案例3：模型确认 - 仿真与实验对比")
print("=" * 70)

def case3_model_validation():
    """
    案例3：模型确认 - 仿真结果与"实验"数据对比
    使用模拟的实验数据（含噪声和不确定性）进行确认分析
    """
    print("\n案例3：模型确认分析")
    print("-" * 50)
    
    np.random.seed(42)
    
    # 物理参数
    L = 1.0  # 域长度 [m]
    kappa = 5.0  # 吸收系数 [1/m]
    sigma_s = 2.0  # 散射系数 [1/m]
    q_gen = 5000.0  # 体积热源 [W/m³]
    T_left = 300.0  # 左边界温度 [K]
    T_right = 350.0  # 右边界温度 [K]
    
    # 仿真求解（使用细网格作为"精确"解）
    N_sim = 100
    x_sim, T_sim, q_sim = solve_diffusion_radiation_1d(N_sim, kappa, sigma_s, 
                                                         q_gen, T_left, T_right, L)
    
    # 生成"实验"数据（仿真结果 + 噪声 + 系统误差）
    n_exp_points = 11  # 实验测点数量
    x_exp = np.linspace(0, L, n_exp_points)
    
    # 实验测量值（插值 + 噪声）
    T_exp_true = np.interp(x_exp, x_sim, T_sim)
    
    # 添加随机测量误差（偶然不确定性）
    measurement_error_std = 2.0  # 标准差 2K
    T_exp_measured = T_exp_true + np.random.normal(0, measurement_error_std, n_exp_points)
    
    # 添加系统误差（认知不确定性）- 模拟模型偏差
    systematic_error = 3.0 * np.sin(np.pi * x_exp / L)  # 位置相关的系统偏差
    T_exp_measured += systematic_error
    
    # 实验不确定度（包含A类和B类）
    u_exp = np.sqrt(measurement_error_std**2 + 2.0**2)  # 合成标准不确定度
    
    print(f"\n实验数据生成:")
    print(f"  测点数量: {n_exp_points}")
    print(f"  测量随机误差标准差: {measurement_error_std} K")
    print(f"  合成实验不确定度: {u_exp:.2f} K")
    
    # 计算确认度量
    # 在实验测点处采样仿真结果
    T_sim_at_exp = np.interp(x_exp, x_sim, T_sim)
    
    # 误差范数
    L2_error = np.sqrt(np.mean((T_sim_at_exp - T_exp_measured)**2))
    Linf_error = np.max(np.abs(T_sim_at_exp - T_exp_measured))
    
    # 相对误差
    relative_error = np.abs(T_sim_at_exp - T_exp_measured) / T_exp_measured * 100
    mean_relative_error = np.mean(relative_error)
    max_relative_error = np.max(relative_error)
    
    # 决定系数 R²
    SS_res = np.sum((T_exp_measured - T_sim_at_exp)**2)
    SS_tot = np.sum((T_exp_measured - np.mean(T_exp_measured))**2)
    R_squared = 1 - SS_res / SS_tot
    
    # 确认指标（考虑实验不确定度）
    validation_metric = np.abs(T_sim_at_exp - T_exp_measured) / u_exp
    mean_validation_metric = np.mean(validation_metric)
    
    print(f"\n确认度量结果:")
    print(f"  L2 误差: {L2_error:.4f} K")
    print(f"  L∞ 误差: {Linf_error:.4f} K")
    print(f"  平均相对误差: {mean_relative_error:.4f}%")
    print(f"  最大相对误差: {max_relative_error:.4f}%")
    print(f"  决定系数 R²: {R_squared:.6f}")
    print(f"  平均确认指标: {mean_validation_metric:.4f}")
    
    if mean_validation_metric < 1.0:
        print(f"  结论: 模型与实验符合良好 (d < 1)")
    elif mean_validation_metric < 3.0:
        print(f"  结论: 模型与实验可接受符合 (1 ≤ d < 3)")
    else:
        print(f"  结论: 模型与实验存在显著差异 (d ≥ 3)")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 1. 仿真与实验对比
    ax1 = axes[0, 0]
    ax1.plot(x_sim, T_sim, 'b-', linewidth=2, label='Simulation')
    ax1.errorbar(x_exp, T_exp_measured, yerr=2*u_exp, fmt='ro', markersize=8, 
                capsize=5, capthick=2, label='Experiment (±2σ)')
    ax1.set_xlabel('Position x [m]')
    ax1.set_ylabel('Temperature [K]')
    ax1.set_title('Model Validation: Simulation vs Experiment')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 残差分析
    ax2 = axes[0, 1]
    residuals = T_sim_at_exp - T_exp_measured
    ax2.bar(range(n_exp_points), residuals, color='steelblue', alpha=0.7)
    ax2.axhline(y=0, color='r', linestyle='-', linewidth=2)
    ax2.axhline(y=u_exp, color='g', linestyle='--', linewidth=1.5, label=f'±u_exp = ±{u_exp:.1f}K')
    ax2.axhline(y=-u_exp, color='g', linestyle='--', linewidth=1.5)
    ax2.set_xlabel('Measurement Point')
    ax2.set_ylabel('Residual (Sim - Exp) [K]')
    ax2.set_title('Residual Analysis')
    ax2.legend()
    ax2.grid(True, alpha=0.3, axis='y')
    
    # 3. 散点图与回归
    ax3 = axes[1, 0]
    ax3.scatter(T_exp_measured, T_sim_at_exp, s=100, c='blue', alpha=0.6, edgecolors='black')
    # 理想线
    T_min = min(T_exp_measured.min(), T_sim_at_exp.min()) - 5
    T_max = max(T_exp_measured.max(), T_sim_at_exp.max()) + 5
    ax3.plot([T_min, T_max], [T_min, T_max], 'r--', linewidth=2, label='Perfect Agreement')
    # 误差带
    ax3.fill_between([T_min, T_max], [T_min-2*u_exp, T_max-2*u_exp], 
                     [T_min+2*u_exp, T_max+2*u_exp], alpha=0.2, color='green', 
                     label=f'±2σ Experimental Uncertainty')
    ax3.set_xlabel('Experimental Temperature [K]')
    ax3.set_ylabel('Simulated Temperature [K]')
    ax3.set_title(f'Scatter Plot (R² = {R_squared:.4f})')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    ax3.set_xlim(T_min, T_max)
    ax3.set_ylim(T_min, T_max)
    ax3.set_aspect('equal')
    
    # 4. 确认指标分布
    ax4 = axes[1, 1]
    ax4.bar(range(n_exp_points), validation_metric, color='coral', alpha=0.7)
    ax4.axhline(y=1.0, color='g', linestyle='--', linewidth=2, label='Acceptance Threshold (d=1)')
    ax4.axhline(y=3.0, color='r', linestyle='--', linewidth=2, label='Rejection Threshold (d=3)')
    ax4.set_xlabel('Measurement Point')
    ax4.set_ylabel('Validation Metric d')
    ax4.set_title(f'Validation Metric Distribution (Mean = {mean_validation_metric:.2f})')
    ax4.legend()
    ax4.grid(True, alpha=0.3, axis='y')
    
    plt.tight_layout()
    plt.savefig('case3_model_validation.png', dpi=150, bbox_inches='tight')
    print("\n案例3图表已保存: case3_model_validation.png")
    plt.close()
    
    return x_exp, T_exp_measured, T_sim_at_exp, validation_metric

# 运行案例3
x_exp, T_exp, T_sim_at_exp, validation_metric = case3_model_validation()

# =============================================================================
# 案例4：不确定性量化 - 蒙特卡洛传播与敏感性分析
# =============================================================================
print("\n" + "=" * 70)
print("案例4：不确定性量化 - 蒙特卡洛传播与敏感性分析")
print("=" * 70)

def case4_uncertainty_quantification():
    """
    案例4：不确定性量化
    - 输入参数不确定性传播
    - 敏感性分析（Sobol指数）
    - 输出不确定性估计
    """
    print("\n案例4：不确定性量化分析")
    print("-" * 50)
    
    np.random.seed(123)
    
    # 问题参数（名义值）
    L = 1.0
    q_gen_nom = 5000.0
    T_left_nom = 300.0
    T_right_nom = 350.0
    
    # 参数不确定性（标准差）
    # 吸收系数（对数正态分布）
    kappa_nom = 5.0
    kappa_cov = 0.15  # 变异系数 15%
    
    # 散射系数（均匀分布）
    sigma_s_nom = 2.0
    sigma_s_range = 0.4  # ±20%
    
    # 蒙特卡洛采样
    n_samples = 1000
    
    print(f"\n蒙特卡洛不确定性传播:")
    print(f"  样本数量: {n_samples}")
    print(f"\n输入参数不确定性:")
    print(f"  吸收系数 κ: {kappa_nom} ± {kappa_cov*100}% (对数正态)")
    print(f"  散射系数 σ_s: {sigma_s_nom} ± {sigma_s_range/sigma_s_nom*100}% (均匀)")
    
    # 生成随机样本
    kappa_samples = np.random.lognormal(np.log(kappa_nom), kappa_cov, n_samples)
    sigma_s_samples = np.random.uniform(sigma_s_nom - sigma_s_range, 
                                        sigma_s_nom + sigma_s_range, n_samples)
    
    # 运行蒙特卡洛仿真
    T_mid_samples = []
    q_left_samples = []
    
    N_sim = 50  # 仿真网格数
    
    print(f"\n运行蒙特卡洛仿真...")
    for i in range(n_samples):
        x, T, q = solve_diffusion_radiation_1d(N_sim, kappa_samples[i], sigma_s_samples[i],
                                                 q_gen_nom, T_left_nom, T_right_nom, L)
        
        # 记录中点温度和左边界热流
        T_mid = np.interp(L/2, x, T)
        q_left = q[0]
        
        T_mid_samples.append(T_mid)
        q_left_samples.append(q_left)
        
        if (i + 1) % 200 == 0:
            print(f"  完成 {i+1}/{n_samples} 样本")
    
    T_mid_samples = np.array(T_mid_samples)
    q_left_samples = np.array(q_left_samples)
    
    # 统计输出不确定性
    T_mid_mean = np.mean(T_mid_samples)
    T_mid_std = np.std(T_mid_samples)
    T_mid_95_low = np.percentile(T_mid_samples, 2.5)
    T_mid_95_high = np.percentile(T_mid_samples, 97.5)
    
    q_left_mean = np.mean(q_left_samples)
    q_left_std = np.std(q_left_samples)
    q_left_95_low = np.percentile(q_left_samples, 2.5)
    q_left_95_high = np.percentile(q_left_samples, 97.5)
    
    print(f"\n输出不确定性统计:")
    print(f"  中点温度 T_mid:")
    print(f"    均值: {T_mid_mean:.2f} K")
    print(f"    标准差: {T_mid_std:.2f} K ({T_mid_std/T_mid_mean*100:.2f}%)")
    print(f"    95%置信区间: [{T_mid_95_low:.2f}, {T_mid_95_high:.2f}] K")
    print(f"  左边界热流 q_left:")
    print(f"    均值: {q_left_mean/1000:.2f} kW/m²")
    print(f"    标准差: {q_left_std/1000:.2f} kW/m² ({abs(q_left_std/q_left_mean)*100:.2f}%)")
    print(f"    95%置信区间: [{q_left_95_low/1000:.2f}, {q_left_95_high/1000:.2f}] kW/m²")
    
    # 敏感性分析（相关系数）
    corr_kappa_T = np.corrcoef(kappa_samples, T_mid_samples)[0, 1]
    corr_sigma_T = np.corrcoef(sigma_s_samples, T_mid_samples)[0, 1]
    corr_kappa_q = np.corrcoef(kappa_samples, q_left_samples)[0, 1]
    corr_sigma_q = np.corrcoef(sigma_s_samples, q_left_samples)[0, 1]
    
    print(f"\n敏感性分析 (Pearson相关系数):")
    print(f"  对中点温度的影响:")
    print(f"    吸收系数 κ: {corr_kappa_T:.4f}")
    print(f"    散射系数 σ_s: {corr_sigma_T:.4f}")
    print(f"  对左边界热流的影响:")
    print(f"    吸收系数 κ: {corr_kappa_q:.4f}")
    print(f"    散射系数 σ_s: {corr_sigma_q:.4f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 1. 中点温度分布
    ax1 = axes[0, 0]
    ax1.hist(T_mid_samples, bins=30, density=True, color='skyblue', alpha=0.7, edgecolor='black')
    ax1.axvline(T_mid_mean, color='r', linestyle='-', linewidth=2, label=f'Mean = {T_mid_mean:.2f}K')
    ax1.axvline(T_mid_95_low, color='g', linestyle='--', linewidth=2, label=f'95% CI')
    ax1.axvline(T_mid_95_high, color='g', linestyle='--', linewidth=2)
    ax1.set_xlabel('Midpoint Temperature [K]')
    ax1.set_ylabel('Probability Density')
    ax1.set_title('Uncertainty Distribution: T_mid')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 左边界热流分布
    ax2 = axes[0, 1]
    ax2.hist(q_left_samples/1000, bins=30, density=True, color='lightcoral', alpha=0.7, edgecolor='black')
    ax2.axvline(q_left_mean/1000, color='r', linestyle='-', linewidth=2, label=f'Mean = {q_left_mean/1000:.2f}kW/m²')
    ax2.axvline(q_left_95_low/1000, color='g', linestyle='--', linewidth=2, label=f'95% CI')
    ax2.axvline(q_left_95_high/1000, color='g', linestyle='--', linewidth=2)
    ax2.set_xlabel('Left Boundary Heat Flux [kW/m²]')
    ax2.set_ylabel('Probability Density')
    ax2.set_title('Uncertainty Distribution: q_left')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 3. 散点图：κ vs T_mid
    ax3 = axes[1, 0]
    ax3.scatter(kappa_samples, T_mid_samples, c='blue', alpha=0.3, s=10)
    # 拟合线
    z = np.polyfit(kappa_samples, T_mid_samples, 1)
    p_fit = np.poly1d(z)
    kappa_sorted = np.sort(kappa_samples)
    ax3.plot(kappa_sorted, p_fit(kappa_sorted), 'r-', linewidth=2, 
            label=f'Linear Fit (r={corr_kappa_T:.3f})')
    ax3.set_xlabel('Absorption Coefficient κ [1/m]')
    ax3.set_ylabel('Midpoint Temperature [K]')
    ax3.set_title(f'Sensitivity: κ vs T_mid')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 4. 散点图：σ_s vs q_left
    ax4 = axes[1, 1]
    ax4.scatter(sigma_s_samples, q_left_samples/1000, c='red', alpha=0.3, s=10)
    # 拟合线
    z = np.polyfit(sigma_s_samples, q_left_samples, 1)
    p_fit = np.poly1d(z)
    sigma_sorted = np.sort(sigma_s_samples)
    ax4.plot(sigma_sorted, p_fit(sigma_sorted)/1000, 'b-', linewidth=2,
            label=f'Linear Fit (r={corr_sigma_q:.3f})')
    ax4.set_xlabel('Scattering Coefficient σ_s [1/m]')
    ax4.set_ylabel('Left Boundary Heat Flux [kW/m²]')
    ax4.set_title(f'Sensitivity: σ_s vs q_left')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case4_uncertainty_quantification.png', dpi=150, bbox_inches='tight')
    print("\n案例4图表已保存: case4_uncertainty_quantification.png")
    plt.close()
    
    return T_mid_samples, q_left_samples, kappa_samples, sigma_s_samples

# 运行案例4
T_mid_samples, q_left_samples, kappa_samples, sigma_s_samples = case4_uncertainty_quantification()

# =============================================================================
# 案例5：制造解方法（MMS）验证
# =============================================================================
print("\n" + "=" * 70)
print("案例5：制造解方法（MMS）验证")
print("=" * 70)

def case5_manufactured_solution():
    """
    案例5：制造解方法（Method of Manufactured Solutions, MMS）
    验证数值实现的正确性
    """
    print("\n案例5：制造解方法验证")
    print("-" * 50)
    
    # 定义制造解
    def manufactured_solution(x, t):
        """
        制造解：空间-时间变化的温度分布
        T(x,t) = 300 + 100*sin(pi*x)*exp(-t/10)
        """
        return 300 + 100 * np.sin(np.pi * x) * np.exp(-t/10)
    
    def manufactured_source(x, t, alpha):
        """
        计算使制造解满足热传导方程的源项
        ∂T/∂t = α * ∂²T/∂x² + Q_mms
        
        解析计算：
        ∂T/∂t = -10*sin(pi*x)*exp(-t/10)
        ∂²T/∂x² = -100*pi²*sin(pi*x)*exp(-t/10)
        
        因此：Q_mms = ∂T/∂t - α*∂²T/∂x²
        """
        dT_dt = -10 * np.sin(np.pi * x) * np.exp(-t/10)
        d2T_dx2 = -100 * np.pi**2 * np.sin(np.pi * x) * np.exp(-t/10)
        Q_mms = dT_dt - alpha * d2T_dx2
        return Q_mms
    
    # 问题参数
    L = 1.0
    alpha = 0.01  # 热扩散系数
    t_final = 5.0
    
    print(f"\nMMS设置:")
    print(f"  制造解: T(x,t) = 300 + 100*sin(πx)*exp(-t/10)")
    print(f"  域长度: L = {L} m")
    print(f"  热扩散系数: α = {alpha} m²/s")
    print(f"  终止时间: t = {t_final} s")
    
    # 在不同网格上求解
    N_values = [10, 20, 40, 80, 160]
    errors_L2 = []
    errors_Linf = []
    
    for N in N_values:
        dx = L / (N - 1)
        x = np.linspace(0, L, N)
        dt = dx**2 / (4 * alpha)  # 稳定性条件
        n_steps = int(t_final / dt) + 1
        dt = t_final / n_steps
        
        # 初始条件（制造解在t=0）
        T = manufactured_solution(x, 0)
        
        # 时间推进（显式格式）
        for n in range(n_steps):
            t = n * dt
            T_new = T.copy()
            
            # 内部节点
            for i in range(1, N-1):
                # 热传导项
                diffusion = alpha * (T[i+1] - 2*T[i] + T[i-1]) / dx**2
                # MMS源项
                Q_mms = manufactured_source(x[i], t, alpha)
                
                T_new[i] = T[i] + dt * (diffusion + Q_mms)
            
            # 边界条件（Dirichlet，来自制造解）
            T_new[0] = manufactured_solution(x[0], t + dt)
            T_new[N-1] = manufactured_solution(x[N-1], t + dt)
            
            T = T_new
        
        # 计算误差
        T_exact = manufactured_solution(x, t_final)
        error = T - T_exact
        
        L2_error = np.sqrt(np.mean(error**2))
        Linf_error = np.max(np.abs(error))
        
        errors_L2.append(L2_error)
        errors_Linf.append(Linf_error)
        
        print(f"  N={N:3d}: L2误差={L2_error:.6e}, L∞误差={Linf_error:.6e}")
    
    # 计算收敛阶数
    h_values = 1.0 / np.array(N_values)
    
    # 使用最后三个点估计收敛阶数
    p_L2 = np.polyfit(np.log(h_values[-3:]), np.log(errors_L2[-3:]), 1)[0]
    p_Linf = np.polyfit(np.log(h_values[-3:]), np.log(errors_Linf[-3:]), 1)[0]
    
    print(f"\n收敛阶数:")
    print(f"  L2范数: {p_L2:.4f} (理论值: 2.0)")
    print(f"  L∞范数: {p_Linf:.4f} (理论值: 2.0)")
    
    # 可视化
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    
    # 1. 收敛性分析
    ax1 = axes[0]
    ax1.loglog(h_values, errors_L2, 'bo-', linewidth=2, markersize=8, label=f'L2 Error (p={p_L2:.2f})')
    ax1.loglog(h_values, errors_Linf, 'rs-', linewidth=2, markersize=8, label=f'L∞ Error (p={p_Linf:.2f})')
    
    # 理论收敛线
    h_theory = np.linspace(h_values[0], h_values[-1], 100)
    ax1.loglog(h_theory, errors_L2[0] * (h_theory/h_values[0])**2, 'g--', 
              linewidth=2, label='Theoretical O(h²)')
    
    ax1.set_xlabel('Grid Spacing h')
    ax1.set_ylabel('Error')
    ax1.set_title('MMS Convergence Analysis')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 最终解与制造解对比（使用最细网格）
    ax2 = axes[1]
    N_fine = N_values[-1]
    x_fine = np.linspace(0, L, N_fine)
    T_numerical_final = T
    T_exact_final = manufactured_solution(x_fine, t_final)
    
    ax2.plot(x_fine, T_numerical_final, 'b-', linewidth=2, label='Numerical Solution')
    ax2.plot(x_fine, T_exact_final, 'r--', linewidth=2, label='Manufactured Solution')
    ax2.set_xlabel('Position x [m]')
    ax2.set_ylabel('Temperature [K]')
    ax2.set_title(f'Solution Comparison (N={N_fine}, t={t_final}s)')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case5_manufactured_solution.png', dpi=150, bbox_inches='tight')
    print("\n案例5图表已保存: case5_manufactured_solution.png")
    plt.close()
    
    return h_values, errors_L2, errors_Linf, p_L2, p_Linf

# 运行案例5
h_values, errors_L2, errors_Linf, p_L2, p_Linf = case5_manufactured_solution()

# =============================================================================
# GIF动画1：网格收敛过程动画
# =============================================================================
print("\n" + "=" * 70)
print("生成动画1：网格收敛过程")
print("=" * 70)

def create_gif1_grid_convergence():
    """
    动画1：展示网格细化过程中的解收敛
    """
    print("\n生成网格收敛动画...")
    
    # 问题参数
    L = 1.0
    kappa = 5.0
    sigma_s = 2.0
    q_gen = 5000.0
    T_left = 300.0
    T_right = 350.0
    
    # 参考解（细网格）
    x_ref, T_ref, _ = solve_diffusion_radiation_1d(200, kappa, sigma_s, q_gen, T_left, T_right, L)
    
    # 不同网格的解
    N_values = [10, 15, 20, 30, 40, 60, 80, 100, 150, 200]
    solutions = []
    
    for N in N_values:
        x, T, _ = solve_diffusion_radiation_1d(N, kappa, sigma_s, q_gen, T_left, T_right, L)
        solutions.append((x, T, N))
    
    # 创建动画
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    
    def update(frame):
        for ax in axes:
            ax.clear()
        
        x, T, N = solutions[frame]
        
        # 左图：当前网格解与参考解对比
        axes[0].plot(x_ref, T_ref, 'r--', linewidth=2, label='Reference (N=200)')
        axes[0].plot(x, T, 'b-o', linewidth=2, markersize=6, label=f'Current (N={N})')
        axes[0].set_xlabel('Position x [m]')
        axes[0].set_ylabel('Temperature [K]')
        axes[0].set_title(f'Grid Convergence (N={N})')
        axes[0].legend()
        axes[0].grid(True, alpha=0.3)
        axes[0].set_xlim(0, L)
        axes[0].set_ylim(290, 420)
        
        # 计算误差
        T_interp = np.interp(x_ref, x, T)
        error = np.abs(T_interp - T_ref)
        max_error = np.max(error)
        
        # 右图：误差分布
        axes[1].plot(x_ref, error, 'g-', linewidth=2)
        axes[1].fill_between(x_ref, error, alpha=0.3, color='green')
        axes[1].set_xlabel('Position x [m]')
        axes[1].set_ylabel('Absolute Error [K]')
        axes[1].set_title(f'Error Distribution (Max Error: {max_error:.2f} K)')
        axes[1].grid(True, alpha=0.3)
        axes[1].set_xlim(0, L)
        
        fig.suptitle(f'Grid Convergence Animation (Frame {frame+1}/{len(solutions)})', 
                    fontsize=14, fontweight='bold')
        
        return axes
    
    anim = animation.FuncAnimation(fig, update, frames=len(solutions), 
                                  interval=500, blit=False)
    anim.save('gif1_grid_convergence.gif', writer='pillow', fps=2, dpi=100)
    print("GIF动画1完成：已保存 gif1_grid_convergence.gif")
    plt.close()

create_gif1_grid_convergence()

# =============================================================================
# GIF动画2：不确定性传播动画
# =============================================================================
print("\n" + "=" * 70)
print("生成动画2：不确定性传播过程")
print("=" * 70)

def create_gif2_uncertainty_propagation():
    """
    动画2：展示蒙特卡洛不确定性传播过程
    """
    print("\n生成不确定性传播动画...")
    
    np.random.seed(42)
    
    # 问题参数
    L = 1.0
    q_gen = 5000.0
    T_left = 300.0
    T_right = 350.0
    kappa_nom = 5.0
    sigma_s_nom = 2.0
    
    # 运行多次仿真，逐步积累统计信息
    n_frames = 30
    samples_per_frame = 20
    
    T_mid_history = []
    
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    
    def update(frame):
        for ax in axes:
            ax.clear()
        
        # 生成新的样本
        for _ in range(samples_per_frame):
            kappa = np.random.lognormal(np.log(kappa_nom), 0.15)
            sigma_s = np.random.uniform(sigma_s_nom - 0.4, sigma_s_nom + 0.4)
            
            x, T, _ = solve_diffusion_radiation_1d(50, kappa, sigma_s, q_gen, T_left, T_right, L)
            T_mid = np.interp(L/2, x, T)
            T_mid_history.append(T_mid)
        
        n_total = len(T_mid_history)
        
        # 左图：累积分布
        axes[0].hist(T_mid_history, bins=20, density=True, color='skyblue', 
                    alpha=0.7, edgecolor='black')
        axes[0].axvline(np.mean(T_mid_history), color='r', linestyle='-', 
                       linewidth=2, label=f'Mean = {np.mean(T_mid_history):.2f}K')
        axes[0].set_xlabel('Midpoint Temperature [K]')
        axes[0].set_ylabel('Probability Density')
        axes[0].set_title(f'Distribution (n={n_total} samples)')
        axes[0].legend()
        axes[0].grid(True, alpha=0.3)
        
        # 右图：统计量收敛
        if n_total > 10:
            means = [np.mean(T_mid_history[:i]) for i in range(10, n_total+1)]
            stds = [np.std(T_mid_history[:i]) for i in range(10, n_total+1)]
            sample_counts = range(10, n_total+1)
            
            ax2_twin = axes[1].twinx()
            
            line1 = axes[1].plot(sample_counts, means, 'b-', linewidth=2, label='Mean')
            line2 = ax2_twin.plot(sample_counts, stds, 'r-', linewidth=2, label='Std Dev')
            
            axes[1].set_xlabel('Number of Samples')
            axes[1].set_ylabel('Mean Temperature [K]', color='b')
            ax2_twin.set_ylabel('Standard Deviation [K]', color='r')
            axes[1].set_title('Statistical Convergence')
            axes[1].grid(True, alpha=0.3)
            axes[1].tick_params(axis='y', labelcolor='b')
            ax2_twin.tick_params(axis='y', labelcolor='r')
            
            # 合并图例
            lines = line1 + line2
            labels = [l.get_label() for l in lines]
            axes[1].legend(lines, labels, loc='center right')
        
        fig.suptitle(f'Uncertainty Propagation (Frame {frame+1}/{n_frames})', 
                    fontsize=14, fontweight='bold')
        
        return axes
    
    anim = animation.FuncAnimation(fig, update, frames=n_frames, 
                                  interval=300, blit=False)
    anim.save('gif2_uncertainty_propagation.gif', writer='pillow', fps=3, dpi=100)
    print("GIF动画2完成：已保存 gif2_uncertainty_propagation.gif")
    plt.close()

create_gif2_uncertainty_propagation()

# =============================================================================
# GIF动画3：验证与确认流程动画
# =============================================================================
print("\n" + "=" * 70)
print("生成动画3：验证与确认流程")
print("=" * 70)

def create_gif3_vv_process():
    """
    动画3：展示V&V流程的完整过程
    """
    print("\n生成V&V流程动画...")
    
    np.random.seed(42)
    
    # 物理参数
    L = 1.0
    kappa = 5.0
    sigma_s = 2.0
    q_gen = 5000.0
    T_left = 300.0
    T_right = 350.0
    
    # 生成"实验"数据
    x_exp = np.linspace(0, L, 11)
    x_fine = np.linspace(0, L, 100)
    _, T_fine, _ = solve_diffusion_radiation_1d(100, kappa, sigma_s, q_gen, T_left, T_right, L)
    T_exp_true = np.interp(x_exp, x_fine, T_fine)
    T_exp = T_exp_true + np.random.normal(0, 2, len(x_exp)) + 3 * np.sin(np.pi * x_exp / L)
    
    # 不同网格的仿真结果
    N_values = [10, 20, 40, 80, 100]
    sim_results = []
    
    for N in N_values:
        x, T, _ = solve_diffusion_radiation_1d(N, kappa, sigma_s, q_gen, T_left, T_right, L)
        sim_results.append((x, T, N))
    
    # 创建动画
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    
    n_frames = len(sim_results) + 5  # 额外帧展示最终确认结果
    
    def update(frame):
        for ax in axes:
            ax.clear()
        
        if frame < len(sim_results):
            # 代码验证和解验证阶段
            x, T, N = sim_results[frame]
            
            # 左图：网格收敛
            axes[0].plot(x_fine, T_fine, 'r--', linewidth=2, alpha=0.5, label='Reference')
            axes[0].plot(x, T, 'b-o', linewidth=2, markersize=6, label=f'Simulation (N={N})')
            axes[0].set_xlabel('Position x [m]')
            axes[0].set_ylabel('Temperature [K]')
            axes[0].set_title(f'Code & Solution Verification (N={N})')
            axes[0].legend()
            axes[0].grid(True, alpha=0.3)
            axes[0].set_xlim(0, L)
            axes[0].set_ylim(290, 400)
            
            # 右图：误差分析
            T_interp = np.interp(x_fine, x, T)
            error = np.abs(T_interp - T_fine)
            axes[1].semilogy(x_fine, error, 'g-', linewidth=2)
            axes[1].set_xlabel('Position x [m]')
            axes[1].set_ylabel('Error [K]')
            axes[1].set_title('Numerical Error Distribution')
            axes[1].grid(True, alpha=0.3)
            axes[1].set_xlim(0, L)
            
            stage = "Verification Stage"
        else:
            # 模型确认阶段
            x, T, N = sim_results[-1]
            
            # 左图：仿真与实验对比
            axes[0].plot(x, T, 'b-', linewidth=2, label='Simulation')
            axes[0].errorbar(x_exp, T_exp, yerr=4, fmt='ro', markersize=8, 
                           capsize=5, label='Experiment (±2σ)')
            axes[0].set_xlabel('Position x [m]')
            axes[0].set_ylabel('Temperature [K]')
            axes[0].set_title('Model Validation: Sim vs Exp')
            axes[0].legend()
            axes[0].grid(True, alpha=0.3)
            axes[0].set_xlim(0, L)
            axes[0].set_ylim(290, 400)
            
            # 右图：残差分析
            T_sim_at_exp = np.interp(x_exp, x, T)
            residuals = T_sim_at_exp - T_exp
            axes[1].bar(range(len(x_exp)), residuals, color='steelblue', alpha=0.7)
            axes[1].axhline(y=0, color='r', linestyle='-', linewidth=2)
            axes[1].axhline(y=4, color='g', linestyle='--', linewidth=1.5, alpha=0.7)
            axes[1].axhline(y=-4, color='g', linestyle='--', linewidth=1.5, alpha=0.7)
            axes[1].set_xlabel('Measurement Point')
            axes[1].set_ylabel('Residual (Sim - Exp) [K]')
            axes[1].set_title('Validation Residuals')
            axes[1].grid(True, alpha=0.3, axis='y')
            
            stage = "Validation Stage"
        
        fig.suptitle(f'V&V Process: {stage} (Frame {frame+1}/{n_frames})', 
                    fontsize=14, fontweight='bold')
        
        return axes
    
    anim = animation.FuncAnimation(fig, update, frames=n_frames, 
                                  interval=600, blit=False)
    anim.save('gif3_vv_process.gif', writer='pillow', fps=1.5, dpi=100)
    print("GIF动画3完成：已保存 gif3_vv_process.gif")
    plt.close()

create_gif3_vv_process()

# =============================================================================
# 总结
# =============================================================================
print("\n" + "=" * 70)
print("仿真完成总结")
print("=" * 70)
print("\n本程序演示了辐射换热仿真中的验证与确认方法：")
print("\n1. 代码验证 (Code Verification):")
print("   - 一维灰体介质辐射平衡基准问题求解")
print("   - 不同网格密度下的收敛性分析")
print("\n2. 解验证 (Solution Verification):")
print("   - 网格收敛性研究")
print("   - Richardson外推估计精确解")
print(f"   - 网格收敛指标 GCI = {GCI*100:.4f}%")
print("\n3. 模型确认 (Model Validation):")
print("   - 仿真与实验数据对比")
print("   - 误差范数和确认指标计算")
print("   - 残差分析和散点图")
print("\n4. 不确定性量化 (Uncertainty Quantification):")
print("   - 蒙特卡洛不确定性传播")
print("   - 输入参数敏感性分析")
print("   - 输出统计量估计")
print("\n5. 制造解方法 (MMS):")
print(f"   - 验证数值实现正确性")
print(f"   - 观察收敛阶数: L2={p_L2:.2f}, L∞={p_Linf:.2f}")
print("\n6. 生成的可视化:")
print("   - case1_code_verification.png")
print("   - case2_solution_verification.png")
print("   - case3_model_validation.png")
print("   - case4_uncertainty_quantification.png")
print("   - case5_manufactured_solution.png")
print("   - gif1_grid_convergence.gif")
print("   - gif2_uncertainty_propagation.gif")
print("   - gif3_vv_process.gif")
print("\n" + "=" * 70)

AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

【无标题】

AtomGit开源社区

深度学习学习笔记（14周）

本周主要学习了机器学习与深度学习中的基础理论，包括神经网络、梯度下降、卷积神经网络、损失函数、图像分类、语义分割、特征图可视化以及评价指标等内容。通过对相关知识的学习，对深度学习模型的训练流程和计算机视觉任务有了更加系统的理解，也为后续学习高级网络结构和医学图像分析奠定了理论基础。FP+FP+FP+