4. 含内热源的一维热传导建模与源项处理

问题设定与解析解基准

本节处理一维稳态热传导方程在含空间分布内热源时的离散与求解。控制方程为：

d/dx( k dT/dx ) + S(x) = 0,  0 < x < L
T(0) = T_L,  T(L) = T_R

取导热系数 k = 1 W/(m·K)，长度 L = 1 m，两端温度 T_L = T_R = 0 °C，内热源采用抛物线分布 S(x) = S0 · 4x(L−x)/L²，峰值 S0 = 100 W/m³ 位于 x = 0.5 m。该源项在两端为零、中间最大，便于观察源项对温度剖面的“鼓包”效应。

对上述 S(x) 积分两次并代入边界条件，可得解析解：

T(x) = (S0/(3kL²)) · (L³x − 2Lx³ + x⁴) / ... 简化后：
T(x) = (S0/(3k)) · (x·L − 2x³/L + x⁴/L²) / L ... 直接给出：
T(x) = (S0/(3k)) · ( L·x − 2x³/L + x⁴/L² ) / L
     = (S0/(3kL)) · ( L²x − 2x³ + x⁴/L )

为避免公式歧义，直接给出可验证形式：T(x) = (S0/(3k)) · (x·L − 2x³/L + x⁴/L²) / L，在 x = 0.5 处峰值 T_max = (S0·L²)/(3k) · (1/2 − 1/4 + 1/16) = (100/3)·(5/16) ≈ 10.4167 °C。该值将作为数值解的校验基准。

有限体积离散与源项处理

采用均匀网格，节点 i = 0..N，间距 Δx = L/N，控制体积宽度 ΔV = Δx（一维）。对控制体积 i 积分控制方程：

[k dT/dx]_{i+1/2} − [k dT/dx]_{i−1/2} + ∫_{i−1/2}^{i+1/2} S(x) dx = 0

扩散项用中心差分，源项用控制体积平均。关键在源项积分 ∫ S dx 的近似方式，这直接决定离散一致性：

• 节点采样（中点法）：∫ S dx ≈ S(x_i)·Δx，二阶精度，实现最简单。
• 控制体积精确积分：对抛物线 S(x) 在 [x_{i−1/2}, x_{i+1/2}] 上解析积分，得到与节点采样不同的等效源强，可检验源项离散误差。

离散方程（内部节点）：

(k/Δx)(T_{i−1} − 2T_i + T_{i+1}) + S_i·Δx = 0

整理为三对角形式 a_i T_{i−1} + b_i T_i + c_i T_{i+1} = d_i，其中 a_i = c_i = k/Δx，b_i = −2k/Δx，d_i = −S_i·Δx。边界节点直接置 T_0 = T_L、T_N = T_R。

Python 实现：源项离散与三对角求解

下面给出完整可运行脚本，包含两种源项离散方式、三对角矩阵求解（Thomas 算法）以及与解析解的误差对比。

import numpy as np

def solve_heat_source(N, L=1.0, k=1.0, S0=100.0, TL=0.0, TR=0.0, source_mode="midpoint"):
    dx = L / N
    x = np.linspace(0.0, L, N + 1)
    # 节点源强
    S_node = S0 * 4.0 * x * (L - x) / L**2

    if source_mode == "midpoint":
        # 中点法：∫S dx ≈ S(x_i)·dx
        S_int = S_node * dx
    elif source_mode == "exact":
        # 抛物线在 [x-dx/2, x+dx/2] 上解析积分
        # S(x) = 4S0/L² (Lx - x²)
        # ∫ = 4S0/L² [L(x²/2) - x³/3] 在区间端点求差
        def F(xx):
            return 4.0 * S0 / L**2 * (L * xx**2 / 2.0 - xx**3 / 3.0)
        S_int = np.array([F(xi + dx/2) - F(xi - dx/2) for xi in x])
    else:
        raise ValueError("source_mode 必须为 'midpoint' 或 'exact'")

    # 组装三对角系统（内部节点 1..N-1）
    n = N - 1
    a = np.full(n, k / dx)      # 下对角
    b = np.full(n, -2.0 * k / dx)  # 主对角
    c = np.full(n, k / dx)      # 上对角
    d = -S_int[1:N]             # 右端项

    # 边界贡献
    d[0]  -= a[0] * TL
    d[-1] -= c[-1] * TR

    # Thomas 算法
    cp = np.zeros(n)
    dp = np.zeros(n)
    cp[0] = c[0] / b[0]
    dp[0] = d[0] / b[0]
    for i in range(1, n):
        denom = b[i] - a[i] * cp[i-1]
        cp[i] = c[i] / denom
        dp[i] = (d[i] - a[i] * dp[i-1]) / denom
    T = np.zeros(N + 1)
    T[0], T[N] = TL, TR
    T[n] = dp[-1]
    for i in range(n - 2, -1, -1):
        T[i+1] = dp[i] - cp[i] * T[i+2]
    return x, T, S_node

def exact_T(x, L=1.0, k=1.0, S0=100.0):
    return (S0 / (3.0 * k)) * (x * L - 2.0 * x**3 / L + x**4 / L**2) / L

# 网格收敛性对比
for N in [10, 20, 40, 80]:
    x, T, _ = solve_heat_source(N, source_mode="midpoint")
    Te = exact_T(x)
    err = np.max(np.abs(T - Te))
    print(f"N={N:3d}  dx={1.0/N:.4f}  max|T-Texact|={err:.3e}  T_mid={T[N//2]:.6f}")

# 两种源项离散方式对比（N=40）
x, Tm, _ = solve_heat_source(40, source_mode="midpoint")
_, Te_, _ = solve_heat_source(40, source_mode="exact")
print("midpoint vs exact 源项离散最大差:", np.max(np.abs(Tm - Te_)))

运行后预期观察：中点法误差随 N 增大按 O(Δx²) 下降；T_mid 收敛到约 10.4167 °C。两种源项离散方式在细网格下差异趋于零，说明源项积分近似与扩散项离散精度匹配时整体保持二阶一致性。

结果判读与一致性检查

数值解与解析解对比时，重点检查三件事：

1. 峰值位置与幅值：源项峰值在 x = 0.5，温度峰值也应位于中点，幅值收敛到 10.4167 °C。若峰值偏移，说明源项采样位置或控制体积划分有误。
2. 边界条件：两端温度必须严格为 0。若边界节点被错误地纳入内部方程求解，会出现非零边界值。
3. 收敛阶：中点法源项离散下，最大误差应随 N 加倍约降为 1/4。若只降为 1/2，说明源项积分精度低于扩散项，整体降为一阶。

下表给出示例网格下的预期误差量级（由上述脚本可直接复现，数值为教学示例的典型量级，非实测承诺）：

源项离散一致性要点

源项处理的本质是保证 ∫S dx 的近似精度不低于扩散项。中点法对光滑源项是二阶精度，与中心差分扩散项匹配；若源项在控制体积内剧烈变化（如集中热源、界面阶跃），中点法会退化，此时应改用控制体积精确积分或对源项做分段处理。本例的抛物线源项是检验这一点的良好基准：两种离散方式在细网格下结果一致，说明实现正确。