2. 一维对流扩散方程迎风格式与中心格式对比

问题设定与离散格式

考虑一维稳态对流扩散方程：

u·dφ/dx = Γ·d²φ/dx²,  0 < x < L

取 L=1，边界条件 φ(0)=0、φ(1)=1，速度 u 与扩散系数 Γ 为常数。定义网格 Péclet 数 Pe_h = u·Δx/Γ，它衡量单个网格间距内对流与扩散的相对强度。均匀网格 N 个单元，Δx = L/N，节点 x_i = i·Δx（i=0..N）。

两种格式对一阶导数项的处理不同：

• 中心差分（CDS）：dφ/dx ≈ (φ_{i+1} − φ_{i−1})/(2Δx)，二阶精度，但系数矩阵在 Pe_h > 2 时失去对角占优，产生非物理振荡。
• 迎风差分（UDS）：u>0 时 dφ/dx ≈ (φ_i − φ_{i−1})/Δx，一阶精度，无条件满足对角占优，但引入数值耗散，边界层被抹平。

解析解为 φ(x) = (exp(Pe·x/L) − 1)/(exp(Pe) − 1)，其中 Pe = u·L/Γ 为全局 Péclet 数。

离散方程与实现

对内部节点 i=1..N−1 写出离散方程。记 a = Γ/Δx²，b = u/(2Δx)（CDS）或 b = u/Δx（UDS，u>0）。

中心格式：(a + b)·φ_{i−1} − 2a·φ_i + (a − b)·φ_{i+1} = 0

迎风格式（u>0）：(a + u/Δx)·φ_{i−1} − (2a + u/Δx)·φ_i + a·φ_{i+1} = 0

组装为三对角系统 A·φ = d，用 Thomas 算法求解。以下 Python 代码完整实现两种格式并计算与解析解的误差：

import numpy as np

def solve_cds_uds(N, Pe_global, L=1.0):
    dx = L / N
    Gamma = 1.0
    u = Pe_global * Gamma / L
    x = np.linspace(0, L, N + 1)
    a = Gamma / dx**2

    def thomas(lower, diag, upper, rhs):
        n = len(diag)
        c = np.zeros(n); d = np.zeros(n)
        c[0] = upper[0] / diag[0]; d[0] = rhs[0] / diag[0]
        for i in range(1, n):
            m = diag[i] - lower[i] * c[i-1]
            c[i] = upper[i] / m if i < n-1 else 0.0
            d[i] = (rhs[i] - lower[i] * d[i-1]) / m
        phi = np.zeros(n)
        phi[-1] = d[-1]
        for i in range(n-2, -1, -1):
            phi[i] = d[i] - c[i] * phi[i+1]
        return phi

    # 中心差分
    b = u / (2 * dx)
    n = N - 1
    lo = np.full(n, a + b); di = np.full(n, -2*a); up = np.full(n, a - b)
    rhs = np.zeros(n); rhs[0] = -(a + b) * 0.0; rhs[-1] = -(a - b) * 1.0
    phi_cds = np.concatenate(([0.0], thomas(lo, di, up, rhs), [1.0]))

    # 迎风格式 (u > 0)
    c = u / dx
    lo = np.full(n, a + c); di = np.full(n, -(2*a + c)); up = np.full(n, a)
    rhs = np.zeros(n); rhs[-1] = -a * 1.0
    phi_uds = np.concatenate(([0.0], thomas(lo, di, up, rhs), [1.0]))

    phi_ex = (np.exp(Pe_global * x / L) - 1) / (np.exp(Pe_global) - 1)
    return x, phi_cds, phi_uds, phi_ex

for Pe in [1, 5, 10, 20]:
    x, pc, pu, pe = solve_cds_uds(N=20, Pe_global=Pe)
    print(f"Pe={Pe:2d}  CDS最大误差={np.max(np.abs(pc-pe)):.4f}  "
          f"UDS最大误差={np.max(np.abs(pu-pe)):.4f}  "
          f"CDS最小值={pc.min():.4f}")

运行后可观察到：Pe=1 时两种格式都接近解析解；Pe=5 时 CDS 开始出现负值（非物理振荡）；Pe=20 时 CDS 振荡剧烈，而 UDS 虽无振荡但边界层被显著抹平。

网格 Péclet 数与格式行为

关键判据是 Pe_h = u·Δx/Γ = Pe/N。CDS 的系数 (a − b) 在 Pe_h > 2 时变号，导致矩阵非对角占优，解出现振荡。UDS 的系数始终同号，保证单调性，但一阶截断误差相当于附加了数值扩散 Γ_num = u·Δx/2，使有效 Péclet 数降低。

结果判读与格式选择

对比两种格式的误差与形态，可归纳出以下判读逻辑：

若必须使用 CDS 且 Pe_h > 2，可加密网格使 Pe_h ≤ 2，或采用混合格式（Pe_h > 2 时切换为迎风）。验证时检查解是否出现负值或超出边界值范围——这是非物理振荡的直接信号；同时比较不同 N 下的最大误差，CDS 在 Pe_h < 2 时应呈二阶收敛，UDS 呈一阶收敛。