主题072：仿真结果验证与确认（V&V）

一、引言

1.1 V&V的重要性

仿真结果验证与确认（Verification and Validation, V&V）是确保仿真模型可靠性和可信度的核心过程。随着仿真技术在工程决策中扮演越来越重要的角色，V&V已成为仿真工作流程中不可或缺的环节。

验证（Verification）：

• 回答"我们是否正确地求解了方程？"
• 关注数值解与数学模型之间的一致性
• 检查代码实现、算法正确性、数值精度

确认（Validation）：

• 回答"我们是否求解了正确的方程？"
• 关注数学模型与真实物理现象之间的一致性
• 评估模型的适用范围和预测能力

不确定性量化（Uncertainty Quantification, UQ）：

• 识别和量化仿真中的各种不确定性来源
• 评估不确定性对结果的影响
• 为决策提供置信区间
  

1.2 V&V的标准与指南

国际标准：

• ASME V&V 10-2006：计算固体力学V&V指南
• ASME V&V 20-2009：计算流体力学和传热学V&V标准
• AIAA G-077-1998：计算流体力学V&V指南
• OECD/NEA：核工程仿真V&V最佳实践

关键概念：

• 确认度量（Validation Metrics）：定量评估模型与实验数据的一致性
• 确认层次（Validation Hierarchy）：从单元问题到全系统确认的递进过程
• 预测能力（Predictive Capability）：模型在未测试条件下的预测可信度

二、代码验证（Code Verification）

2.1 软件质量保证

代码审查：

• 静态代码分析：检查语法错误、潜在bug、代码规范
• 同行评审：由其他开发者审查代码逻辑和实现
• 单元测试：对每个函数/模块进行独立测试

版本控制：

• 使用Git等版本控制系统
• 记录代码变更历史
• 支持代码回溯和分支管理

文档化：

• 用户手册：软件功能和使用方法
• 理论手册：数学模型和数值方法
• 验证报告：验证案例和结果

2.2 数值算法验证

精确解验证：
对于存在解析解的简单问题，直接比较数值解与精确解：
$$\text{误差}=|u_{\text{num}}-u_{\text{exact}}|$$

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

1. 选择一个光滑的解析函数作为"制造解"
2. 将制造解代入控制方程，得到源项
3. 修改代码以包含源项
4. 比较数值解与制造解
5. 验证收敛阶数是否符合理论预期

示例：验证一维热传导方程求解器

# 制造解：u(x,t) = sin(πx) * exp(-π²t)
# 代入方程 ∂u/∂t = ∂²u/∂x² 得到源项 f = 0
# 验证数值解的收敛阶数

2.3 网格收敛性分析

Richardson外推法：
利用不同网格上的计算结果估计精确解：
$$f_{\text{exact}}\approx f_h+\frac{f_h-f_{rh}}{r^p-1}$$

其中：

• $f_h$：细网格解
• $f_{rh}$：粗网格解（网格尺寸为 $rh$）
• $r$：网格细化比例
• $p$：收敛阶数

网格收敛指标（Grid Convergence Index, GCI）：
$$GCI=\frac{F_s|ϵ|}{r^p-1}$$

其中：

• $F_s$：安全系数（通常取1.25）
• $ϵ=\frac{f_{rh}-f_h}{f_h}$：相对变化

收敛阶数估计：
$$p=\frac{\ln \left(\frac{f_{2h}-f_{4h}}{f_h-f_{2h}}\right)}{\ln (r)}$$

2.4 时间步长收敛性

对于瞬态问题，需要验证时间步长的影响：

时间收敛分析：

1. 固定空间网格，改变时间步长
2. 计算不同时间步长下的解
3. 验证时间收敛阶数
4. 确定合适的时间步长

稳定性限制：

• CFL条件：$\Delta t\le C\frac{\Delta x}{|u|}$
• 扩散限制：$\Delta t\le \frac{\Delta x^2}{2D}$

三、计算验证（Calculation Verification）

3.1 迭代收敛性

残差监控：

• 记录每次迭代的残差范数
• 验证残差是否单调下降
• 确定收敛标准

收敛判据：

• 绝对残差：$\Vert R\Vert <{ϵ}_{\text{abs}}$
• 相对残差：$\frac{\Vert R^{(k)}\Vert }{\Vert R^{(0)}\Vert }<{ϵ}_{\text{rel}}$
• 解的变化：$\Vert u^{(k+1)}-u^{(k)}\Vert <{ϵ}_{\text{sol}}$

迭代次数分析：

• 记录达到收敛所需的迭代次数
• 分析迭代次数与问题规模的关系
• 评估求解器效率

3.2 数值误差估计

截断误差：
$$\tau =L_h(u)-L(u)=O(h^p)$$

舍入误差：

• 浮点数精度限制
• 误差累积和传播
• 条件数分析

数值扩散与耗散：

• 人工粘性影响
• 色散误差分析
• 相位误差

3.3 敏感性分析

输入参数敏感性：

• 改变模型参数，观察结果变化
• 计算敏感性系数：$S_i=\frac{\partial y}{\partial x_i}$
• 识别关键参数

边界条件敏感性：

• 测试不同边界条件的影响
• 评估边界条件的不确定性
• 验证边界处理的正确性

初始条件敏感性：

• 对于瞬态问题，测试初始条件的影响
• 评估系统的稳定性
• 识别混沌行为

四、模型确认（Model Validation）

4.1 实验设计

确认实验层次：

1. 单元问题（Unit Problems）：隔离单个物理现象
2. 基准问题（Benchmark Problems）：综合多个物理过程
3. 子系统（Subsystems）：部分真实系统
4. 全系统（Complete System）：完整真实系统

实验设计原则：

• 覆盖模型的适用范围
• 包含边界条件变化
• 考虑参数敏感性
• 重复性验证

测量不确定性：

• 校准测量设备
• 估计测量误差
• 记录环境条件

4.2 确认度量

定性比较：

• 可视化对比（云图、曲线）
• 流动特征识别
• 趋势一致性

定量度量：

1. 误差范数：
$$L_2\text{误差}=\sqrt{\frac{{\sum }_{i=1}^N(y_{\text{sim},i}-y_{\text{exp},i}{)}^2}{N}}$$

2. 相对误差：
$${ϵ}_{\text{rel}}=\frac{|y_{\text{sim}}-y_{\text{exp}}|}{|y_{\text{exp}}|}$$

3. 决定系数（R²）：
$$R^2=1-\frac{\sum (y_{\text{exp}}-y_{\text{sim}}{)}^2}{\sum (y_{\text{exp}}-{\bar{y}}_{\text{exp}}{)}^2}$$

4. 确认度量（Validation Metric）：
$$d=\frac{|y_{\text{sim}}-y_{\text{exp}}|}{\sqrt{U_{\text{sim}}^2+U_{\text{exp}}^2}}$$

其中 $U$ 表示不确定性。

5. 特征比较：

• 峰值误差
• 相位误差
• 频率响应

4.3 确认矩阵

建立仿真结果与实验数据的系统对比：

物理量	实验值	仿真值	误差	接受标准	状态
最大应力	250 MPa	245 MPa	2%	<5%	✓
固有频率	100 Hz	98 Hz	2%	❤️%	✓
温度分布	-	-	-	-	-

4.4 模型适用性评估

适用范围确定：

• 参数空间边界
• 几何限制
• 物理条件限制

外推风险评估：

• 模型在验证范围外的可信度
• 外推的不确定性增长
• 保守性设计建议

模型改进方向：

• 识别模型缺陷
• 提出改进建议
• 优先级排序

五、不确定性量化

5.1 不确定性来源

偶然不确定性（Aleatory Uncertainty）：

• 固有的随机性
• 不可减少
• 用概率分布描述

认知不确定性（Epistemic Uncertainty）：

• 知识不足导致
• 可通过更多数据减少
• 用区间或概率分布描述

数值不确定性：

• 离散化误差
• 迭代误差
• 舍入误差

5.2 不确定性传播

蒙特卡洛方法：

1. 从输入分布中随机采样
2. 运行仿真
3. 重复多次
4. 统计分析输出

泰勒级数方法：
$$y\approx y_0+\sum _{i=1}^n\frac{\partial y}{\partial x_i}\Delta x_i$$

多项式混沌展开：
$$y(\xi )=\sum _{k=0}^P{y}_k{\Psi }_k(\xi )$$

5.3 验证与确认的置信度

置信度评估：

• 验证置信度：基于代码验证的完整性
• 确认置信度：基于实验对比的充分性
• 综合置信度：考虑所有不确定性来源

置信度等级：

• Level 1：初步验证，有限确认
• Level 2：部分验证，部分确认
• Level 3：全面验证，充分确认
• Level 4：严格验证，广泛确认

六、V&V最佳实践

6.1 V&V计划

计划内容：

• V&V目标和范围
• 验证策略和方法
• 确认实验设计
• 资源和时间安排
• 接受标准

文档要求：

• 问题描述
• 数学模型
• 数值方法
• 验证案例
• 确认实验
• 结果分析
• 结论和建议

6.2 自动化V&V工具

代码验证工具：

• 代码覆盖率分析
• 静态分析工具
• 回归测试框架
• MMS自动生成

网格收敛工具：

• 自动网格细化
• Richardson外推
• GCI计算
• 收敛阶数验证

确认分析工具：

• 实验数据管理
• 自动对比分析
• 确认度量计算
• 报告生成

6.3 持续V&V

版本控制集成：

• 每次代码变更自动运行测试
• 回归测试防止引入错误
• 性能基准监控

回归测试：

• 建立验证案例库
• 定期运行全部案例
• 对比历史结果

知识管理：

• 积累V&V经验
• 建立最佳实践库
• 培训新成员

"""
主题072：区域法辐射换热
一维区域法求解参与性介质辐射
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve
import os

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

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

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

sigma = 5.670374419e-8

def zone_method_1d(L, N, kappa, eps_wall, T_left, T_right, q_vol=0):
    """
    一维区域法求解参与性介质辐射换热
    L: 介质厚度
    N: 区域数
    kappa: 吸收系数
    eps_wall: 壁面发射率
    T_left, T_right: 左右壁面温度
    q_vol: 体积热源
    """
    dx = L / N
    x = np.linspace(dx/2, L - dx/2, N)  # 区域中心位置
    
    # 直接交换面积计算（简化模型）
    # 气体-气体交换面积
    gg = np.zeros((N, N))
    for i in range(N):
        for j in range(N):
            if i == j:
                # 自交换面积
                gg[i, j] = 4 * kappa * dx * (1 - 0.5 * kappa * dx)
            else:
                # 区域间交换面积
                r = abs(i - j) * dx
                gg[i, j] = kappa**2 * dx**2 * np.exp(-kappa * r) / (4 * np.pi * r**2) if r > 0 else 0
    
    # 气体-壁面交换面积
    gs_left = np.zeros(N)
    gs_right = np.zeros(N)
    for i in range(N):
        r_left = x[i]
        r_right = L - x[i]
        gs_left[i] = kappa * dx * np.exp(-kappa * r_left) / (4 * np.pi * r_left**2)
        gs_right[i] = kappa * dx * np.exp(-kappa * r_right) / (4 * np.pi * r_right**2)
    
    # 构建能量平衡方程
    # A * T^4 = B
    A = np.zeros((N, N))
    B = np.zeros(N)
    
    for i in range(N):
        # 自身发射项
        A[i, i] = 4 * kappa * dx
        
        # 与其他气体区域的交换
        for j in range(N):
            if i != j:
                A[i, i] += gg[i, j]
                A[i, j] = -gg[i, j]
        
        # 与壁面的交换
        A[i, i] += gs_left[i] + gs_right[i]
        
        # 右端项
        B[i] = q_vol * dx + gs_left[i] * sigma * T_left**4 + gs_right[i] * sigma * T_right**4
    
    # 求解
    E_b = solve(A, B)
    T_gas = (E_b / sigma)**0.25
    
    return x, T_gas

print("="*60)
print("仿真1：一维区域法")
print("="*60)

# 参数
L = 1.0  # m
N = 20
kappa = 1.0  # 1/m
eps_wall = 0.8
T_left = 1000  # K
T_right = 300  # K

# 求解
x, T_gas = zone_method_1d(L, N, kappa, eps_wall, T_left, T_right)

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

# 1. 温度分布
axes[0].plot(x*100, T_gas, 'b-o', linewidth=2, markersize=8)
axes[0].axhline(y=T_left, color='r', linestyle='--', linewidth=2, label=f'T_left = {T_left}K')
axes[0].axhline(y=T_right, color='g', linestyle='--', linewidth=2, label=f'T_right = {T_right}K')
axes[0].set_xlabel('Position (cm)', fontsize=11)
axes[0].set_ylabel('Gas Temperature (K)', fontsize=11)
axes[0].set_title(f'1D Zone Method (κ={kappa} m⁻¹)', fontsize=12)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 2. 不同吸收系数的影响
kappa_values = [0.1, 0.5, 1.0, 2.0, 5.0]
colors = ['blue', 'green', 'orange', 'red', 'purple']

for kappa_val, color in zip(kappa_values, colors):
    _, T_kappa = zone_method_1d(L, N, kappa_val, eps_wall, T_left, T_right)
    axes[1].plot(x*100, T_kappa, color=color, linewidth=2, label=f'κ = {kappa_val} m⁻¹')

axes[1].set_xlabel('Position (cm)', fontsize=11)
axes[1].set_ylabel('Gas Temperature (K)', fontsize=11)
axes[1].set_title('Effect of Absorption Coefficient', fontsize=12)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/zone_method_1d.png', dpi=150, bbox_inches='tight')
plt.close()
print("图1：一维区域法图已保存")

print("\n" + "="*60)
print("仿真2：二维区域法")
print("="*60)

def zone_method_2d(Lx, Ly, Nx, Ny, kappa, T_boundary):
    """
    二维区域法简化实现
    """
    dx = Lx / Nx
    dy = Ly / Ny
    
    x = np.linspace(dx/2, Lx - dx/2, Nx)
    y = np.linspace(dy/2, Ly - dy/2, Ny)
    X, Y = np.meshgrid(x, y)
    
    # 简化：使用迭代法求解
    T = np.ones((Ny, Nx)) * T_boundary
    
    for iteration in range(100):
        T_old = T.copy()
        
        for i in range(Ny):
            for j in range(Nx):
                # 简化能量平衡
                # 考虑与相邻区域的交换和辐射损失
                neighbors = []
                weights = []
                
                # 四个方向的邻居
                if i > 0:
                    neighbors.append(T_old[i-1, j])
                    weights.append(1.0)
                if i < Ny-1:
                    neighbors.append(T_old[i+1, j])
                    weights.append(1.0)
                if j > 0:
                    neighbors.append(T_old[i, j-1])
                    weights.append(1.0)
                if j < Nx-1:
                    neighbors.append(T_old[i, j+1])
                    weights.append(1.0)
                
                # 加权平均
                if neighbors:
                    T[i, j] = np.average(neighbors, weights=weights)
                
                # 辐射损失（简化）
                T[i, j] -= 0.01 * kappa * dx * (T[i, j] - T_boundary)
        
        # 边界条件
        T[0, :] = T_boundary
        T[-1, :] = T_boundary
        T[:, 0] = T_boundary
        T[:, -1] = T_boundary
        
        # 检查收敛
        if np.max(np.abs(T - T_old)) < 0.1:
            break
    
    return X, Y, T

# 参数
Lx, Ly = 1.0, 1.0
Nx, Ny = 20, 20
kappa = 1.0
T_boundary = 1000  # K

# 求解
X, Y, T_2d = zone_method_2d(Lx, Ly, Nx, Ny, kappa, T_boundary)

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

# 1. 温度分布云图
contour = axes[0].contourf(X*100, Y*100, T_2d, levels=20, cmap='jet')
axes[0].set_xlabel('x (cm)', fontsize=11)
axes[0].set_ylabel('y (cm)', fontsize=11)
axes[0].set_title('2D Zone Method Temperature Distribution', fontsize=12)
plt.colorbar(contour, ax=axes[0], label='Temperature (K)')

# 2. 中心线温度分布
center_idx = Ny // 2
axes[1].plot(X[center_idx, :]*100, T_2d[center_idx, :], 'b-', linewidth=2, label='Horizontal Center')
axes[1].plot(Y[:, Nx//2]*100, T_2d[:, Nx//2], 'r--', linewidth=2, label='Vertical Center')
axes[1].set_xlabel('Position (cm)', fontsize=11)
axes[1].set_ylabel('Temperature (K)', fontsize=11)
axes[1].set_title('Centerline Temperature', fontsize=12)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/zone_method_2d.png', dpi=150, bbox_inches='tight')
plt.close()
print("图2：二维区域法图已保存")

print("\n所有仿真完成！")