3.1 计算机的难处

第 2 讲得到的方程里全是 ∂T∂t∂t∂T​、∂T∂r∂r∂T​ 这种“连续变化率”。计算机不会处理连续的东西——它只认识一个个孤立的数字。

所以要做一个“降维”操作：

• 空间上切：把半径 0~2 cm 切成薄片。题目要求“每隔 0.1 cm”输出，那干脆就按 0.1 cm 切，得到 21 个位置点：0, 0.1, 0.2, …, 2.0 cm；

• 时间上切：题目要求“每隔 1 s”输出，那就按 1 s 一步往前推，1800 s 就是 1800 步。

这样，T(r,t)T(r,t) 就变成了一张大表格：

初始条件（28℃）填好了第一列。我们的任务：用方程推出"从第 t 列填出第 t+1 列"的规则，然后让计算机重复 1800 次。

先约定记号，避免混乱：

Tin≡位置 i 在时刻 n 的温度Tin​≡位置 i 在时刻 n 的温度

• ii：空间格子编号，i=0i=0 是中心，i=N=20i=N=20 是表面，Δr=0.001Δr=0.001 m（0.1 cm）；

• nn：时间步编号，n=0n=0 是开始，Δt=1Δt=1 s。

水分浓度 CinCin​ 完全同理。

3.2 导数变成“差”

核心思想：导数的定义本来就是"差 ÷ 距离"的极限，那我们干脆不取极限，直接用差：

∂T∂t∣i≈Tin+1−TinΔt（用下一时刻和当前时刻的差）∂t∂T​​i​≈ΔtTin+1​−Tin​​（用下一时刻和当前时刻的差）

空间导数同理，用相邻格子的差。但方程右边那个 1r∂∂r(kr∂T∂r)r1​∂r∂​(kr∂r∂T​) 有两层导数嵌套，直接展开容易错。我们换一个更聪明、更不容易错的办法——有限体积法（就是第 2 讲“洋葱圈记账”的直接翻版）。

3.3 有限体积法：把记账直接写下来

回忆第 2 讲的记账规则：每个格子存热的变化 = 流进 − 流出。

内部格子 ii（1 ≤ i ≤ N−1）

格子 ii 是介于 ri−1/2ri−1/2​ 和 ri+1/2ri+1/2​ 之间的一个薄环（ri+1/2=(i+12)Δrri+1/2​=(i+21​)Δr 是格子两端的界面）。

界面 i+12i+21​ 处的热流量（傅里叶定律，导数用温差近似）：

qi+1/2=−k⋅2πri+1/2L⏟圆筒面积⋅Ti+1−TiΔrqi+1/2​=−k⋅圆筒面积2πri+1/2​L​​⋅ΔrTi+1​−Ti​​

界面 i−12i−21​ 处同理：qi−1/2=−k⋅2πri−1/2L⋅Ti−Ti−1Δrqi−1/2​=−k⋅2πri−1/2​L⋅ΔrTi​−Ti−1​​

记账（格子体积 Vi=2πriΔr LVi​=2πri​ΔrL）：

ρcpViTin+1−TinΔt=qi−1/2−qi+1/2ρcp​Vi​ΔtTin+1​−Tin​​=qi−1/2​−qi+1/2​

两边把公因子 2πL2πL 约掉，整理得（这一步是纯代数，你自己拿草稿纸验证一遍会非常有收获）：

  Tin+1−TinΔt=αriΔr2[ ri+1/2(Ti+1n+1−Tin+1)−ri−1/2(Tin+1−Ti−1n+1)]  ΔtTin+1​−Tin​​=ri​Δr2α​[ri+1/2​(Ti+1n+1​−Tin+1​)−ri−1/2​(Tin+1​−Ti−1n+1​)]​

其中 α=k/(ρcp)α=k/(ρcp​)。你看到了什么？方程右边用的全是 n+1n+1 时刻（下一时刻）的温度——这叫隐式格式，马上解释为什么用它。

水分方程：一模一样

把 TT 换成 CC、αα 换成 DD（注意 DD 依赖 CC，每个格子取当前值，见 3.6 节）：

Cin+1−CinΔt=DiriΔr2[ ri+1/2(Ci+1n+1−Cin+1)−ri−1/2(Cin+1−Ci−1n+1)]ΔtCin+1​−Cin​​=ri​Δr2Di​​[ri+1/2​(Ci+1n+1​−Cin+1​)−ri−1/2​(Cin+1​−Ci−1n+1​)]

结构完全相同——再次印证热和水是“同一种数学”。

3.4 两个边界格子怎么办？

① 中心格子 i=0i=0：处理“除以零”

方程里有 1rr1​，而 r0=0r0​=0，直接算会爆炸。用对称性解决：

中心格子是一个半径 Δr22Δr​ 的小圆柱。它只有外表面（r=Δr/2r=Δr/2）跟外界交换热量，中心点没有热流出（对称）。所以记账变成：

ρcpπ(Δr2)2L⏟小圆柱体积T0n+1−T0nΔt=−k⋅2πΔr2L⋅T1n+1−T0n+1Δr⏟从格子1流进来的热ρcp​小圆柱体积π(2Δr​)2L​​ΔtT0n+1​−T0n​​=从格子1流进来的热−k⋅2π2Δr​L⋅ΔrT1n+1​−T0n+1​​​​

整理后：

T0n+1−T0nΔt=4αΔr2(T1n+1−T0n+1)ΔtT0n+1​−T0n​​=Δr24α​(T1n+1​−T0n+1​)

水分版把 αα 换成 D0D0​ 即可。记住这个结果：中心格子的公式比内部格子简单，只是系数换成了 4α/Δr24α/Δr2 乘以 (T1−T0)(T1​−T0​)。

② 表面格子 i=Ni=N：把“对流”接进来

表面格子向外的一面不是药材，而是热空气。所以它的记账是：

存热变化=qN−1/2⏟从里面格子流进的+q对流⏟空气交给表面的存热变化=从里面格子流进的qN−1/2​​​+空气交给表面的q对流​​​

对流热（牛顿冷却定律）：q对流=h⋅2πrNL⋅(Tairn−TN)q对流​=h⋅2πrN​L⋅(Tairn​−TN​)

其中 TairnTairn​ 就是附件1在 nΔtnΔt 秒时的烘房温度（用线性插值取值，附件1只有每60s一个点）。

同样体积用薄环体积 VN=π(rN2−rN−1/22)LVN​=π(rN2​−rN−1/22​)L，整理即可。水分版：

存水变化=(从格子N−1扩散进来的)−hm⋅2πrNL⋅(CN−Cairn)存水变化=(从格子N−1扩散进来的)−hm​⋅2πrN​L⋅(CN​−Cairn​)

注意：CairnCairn​ 也是附件1查表插值。

3.5 为什么用“隐式”（用下一时刻的值）？

对比一下两种写法：

• 显式格式：右边用 nn 时刻（已知）的值 → 直接算出 Tin+1Tin+1​，一步到位，很简单；

• 隐式格式：右边用 n+1n+1 时刻（未知）的值 → 21 个未知数联立，要解方程组。

显式这么简单，为什么不用？因为它有个致命限制——稳定性条件：

Δt≤Δr22αΔt≤2αΔr2​

代入数字：Δr=0.001Δr=0.001 m（0.1 cm），α≈1.7×10−7α≈1.7×10−7，得

Δt≤10−63.4×10−7≈3 sΔt≤3.4×10−710−6​≈3 s

问题1要求输出间隔 1 s，勉强能过；但问题2的水分扩散用附录3公式，D=2.4×10−3e−0.45Ce−3850/TD=2.4×10−3e−0.45Ce−3850/T，在高温干燥时 DD 可以达到 2.4×10−3×e−0.45×0.15×e−3850/323≈2×10−92.4×10−3×e−0.45×0.15×e−3850/323≈2×10−9 m²/s 量级——注意热扩散率 αα 那边问题不大，但水分的稳定性要求 Δt≤Δr2/(2D)Δt≤Δr2/(2D)，只要 DD 不超过 5×10−75×10−7 都没问题。真正的问题在于：显式格式在边界系数大时（本题 hh 较大、格子薄）容易出现数值振荡。

结论：用隐式格式最稳妥——它无条件稳定，Δt 取多大都不会发散，多花的那点计算量（每步解一个 21 元三对角方程组）对计算机来说九牛一毛。写论文时这也是个可以论述的选型理由。

3.6 变系数怎么办？（为问题2做铺垫）

问题1里 αα 是常数，一切好办。问题2里 cp,k,Dcp​,k,D 每个格子、每个时刻都不同。处理办法简单粗暴但有效——“冻结”系数：

在第 nn 步计算第 n+1n+1 步时，所有系数（ρ,cp,k,Dρ,cp​,k,D）都用 nn 时刻的已知值算出来，当作常数代入方程组。

因为 Δt=1Δt=1 s 很短，物性在一秒内变化极小，这样引入的误差完全可以忽略。如果论文想更严谨，可以在每步内做 1~2 次迭代：先解一遍得到 n+1n+1 的值 → 用新值重算系数 → 再解一遍 → 直到变化可忽略。数值实验中通常一次就够。

3.7 完整算法流程（伪代码）

到这里，问题1的算法已经完全成型。把它写成流程：

输入：Δr=0.001 m，N=20（21个格子），Δt=1 s，共1800步
初始化：T[i] = 28.0，C[i] = 2.55（所有 i）

for n = 0 to 1799:            # 每步推进1秒
    T_air = 插值(附件1温度表, (n+1)*Δt)
    C_air = 插值(附件1浓度表, (n+1)*Δt)
    D[i]  = 7e-9 * exp(-0.89*C[i])   # 每个格子的扩散系数

    # 组装 21×21 三对角方程组：
    #   格子0：中心公式（系数 4D/Δr²）
    #   格子1~19：内部公式（含 r_{i±1/2}）
    #   格子20：表面公式（含 h 和 h_m 对流项）
    # 解出 T_new[i]，解出 C_new[i]

    T = T_new
    C = C_new

    # 顺手把 (n+1, 各位置) 的结果写进 result1.xlsx

每步的计算量：解一个 21×21 三对角方程组，用Thomas 算法（追赶法）只需约 60 次乘加，1800 步总共不到 0.1 秒。问题2要算 259200 步，也就几秒钟。

3.8 本讲小结

1. 离散化 = 把连续的 T(r,t)T(r,t) 变成表格 TinTin​；

2. 导数 = 相邻格子的差；圆柱的 1rr1​ 几何因子通过有限体积法的 ri±1/2ri±1/2​ 自然体现，不用死记；

3. 三个特殊位置各有各的公式：中心（对称、系数×4）、内部（标准式）、表面（把对流项 h(Tair−TN)h(Tair​−TN​)、hm(CN−Cair)hm​(CN​−Cair​) 接进记账）；

4. 用隐式格式保证无条件稳定；变系数用“冻结”处理；

5. 整个求解 = 一个三对角方程组 × 每秒一次 × 几十万次。

自测两题

1. 如果把 Δr 减半（0.05 cm），稳定性允许的显式格式最大 Δt 变为多少？（提示：Δr2Δr2 变成 1/4 → 允许 Δt 变为 1/4，约 0.75 s——这就不够 1 s 输出间隔了，更说明隐式格式的必要）

2. 表面格子 i=Ni=N 的记账里，为什么没有“流出到 rN+1/2rN+1/2​”这一项？（提示：那外面已经是空气了，热量交接改由对流项负责，边界即“墙”）

下一讲预告

理论全通了！第 4 讲我给出问题1的完整可运行 Python 代码（逐行注释，你不需要会编程也能看懂结构），运行后自动生成 result1.xlsx 和论文中表1、表2 的数值，并画出温度/水分的时空分布图。之后第 5 讲在它基础上扩展问题2、3、4（变物性、烘干终点判断、移动边界），你会发现 90% 的代码是复用的。