第 4 讲：问题2——全程烘干模型
4.0 先看题目到底要什么
“请建立整个烘干过程药材温度和水分浓度变化规律的数学模型……在论文中按表3、表4给出 3 h 内每隔 0.5 h 的结果，并将每隔 1 s、距离每隔 0.1 cm 的完整结果存入 result2.xlsx。”
注意一个容易看漏的点：论文表格只要求前 3 小时，但模型必须能算全程（2~3 天）——因为问题3要求烘干时间，那个模型就是问题2的模型多跑几十万步。所以：
问题2的真正交付物 = 一个“全程可跑”的模型。3 小时的表只是展示窗口，result2.xlsx 只存前 3 小时（0~10800 s）。
4.1 升级点一：物性从“常数”变成“状态的函数”
对比一下，问题1（附录2）和问题2（附录3）的参数：
关键变化：现在 ρ,cp,kρ,cp,k 都依赖水分浓度 CC，DD 同时依赖 CC 和 TT。
4.2 升级点二：方程真正“耦合”了
第一题讲过，问题1里两个方程是解耦的。现在不一样了，看依赖链：
T→ e−3850/T D→ C,C→ cp(C), k(C) TT e−3850/T D C,C cp(C),k(C) T
温度高 → 水分扩散快（Arrhenius 型，TT 是开尔文！马上讲这个坑）；
水分多 → 热容大、导热快。
谁也离不开谁，必须每个时间步同时解两个方程。这就是第 3 讲“冻结系数”策略的用武之地：
每一步开始时，用当前已知的 Tin,CinTin,Cin 算出每个格子的 ρi,cp,i,ki,Diρi,cp,i,ki,Di，当作常数；然后温度方程和水分方程各解一个三对角方程组（两个方程组之间这一步内不迭代也行，迭代一次更严谨）；解完更新 T,CT,C，进入下一步。
伪代码只需要在第一题的循环里加三行：
for n in ...:
T_air, C_air = 插值(附件1, (n+1)*dt)
# ---- 新增：更新每个格子的物性 ----
ρ[i], cp[i], k[i] = 附录3公式(C[i]) # 用当前 C
D[i] = 2.4e-3 * exp(-0.45*C[i]) * exp(-3850/(T[i]+273.15))
# ---- 解两个三对角方程组（同第一题）----
T_new = solve_tridiagonal(温度组装(ρ,cp,k,h,T_air))
C_new = solve_tridiagonal(水分组装(D,hm,C_air))
T, C = T_new, C_new
计算量：全程 259200 s = 259200 步 × 2 个 21 元三对角求解，纯 Python 也就十几秒，加一版 numpy 向量化更快。完全不是负担。
4.3 三个必须写进论文的“坑”
坑①：开尔文与摄氏度（全军覆没高发区）
DD 公式里的 e−3850/Te−3850/T，题目明确写了 T 为药材温度（单位：K）。这是 Arrhenius 指数的标准形式。代入时必须：
TK=T℃+273.15TK=T℃+273.15
验证一下量级就知道为什么必须是开尔文：T=323T=323 K（50℃）时 e−3850/323=e−11.92≈6.6×10−6e−3850/323=e−11.92≈6.6×10−6；如果错用 50 代入，e−3850/50=e−77≈3×10−34e−3850/50=e−77≈3×10−34——水分扩散系数直接归零，永远烘不干，程序一跑结果就荒谬。这个错连程序都会“看起来正常”地跑完，结果全错，所以必须靠量级检查抓出来。
顺带验证 DD 的量级（323 K、C=2 时）：
D≈2.4×10−3×e−0.9×e−11.92≈2.4×10−3×0.4×6.6×10−6≈6.4×10−9 m2/sD≈2.4×10−3×e−0.9×e−11.92≈2.4×10−3×0.4×6.6×10−6≈6.4×10−9 m2/s
比问题1的 DD（~7×10⁻¹⁰）大约一个数量级——高温让扩散快了 10 倍，这正是“预热之后才能有效干燥”的物理原因，也是烘房要先预热平衡的工程逻辑。这个对比可以写进论文，把问题1和问题2串成有因果的故事。
坑②：附件1只给到 14400 s，之后烘房条件怎么办？
附件1的烘房温度/浓度数据只到 4 小时，但烘干要 2~3 天。必须外推，这是问题2特有的建模决策。观察附件1末段：温度稳定在 50℃ 附近小幅波动，浓度稳定在 0.05 附近——典型的“恒温干燥阶段工况”。
合理假设（写进假设一节）：
恒温干燥阶段（t > 14400 s）烘房工况恒定：Tair=50Tair=50 ℃，Cair=0.05Cair=0.05 kg/kg（取附件1末段均值）。
可以在论文里画一张烘房工况图：实线（附件1数据）+ 虚线（外推段），一目了然。这一条是问题2相对问题1唯一的模型外假设，务必显式写出。
坑③：物性公式的量级自检
代入典型值检查附录3公式是否合理（防抄错系数）：
C=2.55（湿药材）：ρ = 650+128×2.55 ≈ 976 kg/m³，c_p ≈ 1450+2736×0.718 ≈ 3415 J/(kg·K)，k ≈ 0.21+0.38×0.718 ≈ 0.48 W/(m·K)——比问题1的常数物性偏大，符合“含水多→又重又吸热又导热”的直觉；
C=0.1（干药材）：ρ ≈ 663，c_p ≈ 1698，k ≈ 0.245，D 还受温度加成——药材越烘越“轻、省热、难导热”。
两条物理趋势都对，说明公式抄对了。抄系数是这类题最常见的失分点，自检 30 秒，值得。
4.4 3 天里会发生什么？——先做定性预测，再验证数值结果
跑代码之前，先预测结果长什么样（这也是评审爱看的“解的定性分析”）：
温度：
0~1 h：跟随问题1的预热，整根药材从 28℃ 快速爬升；
1~3 h：中心也逐渐逼近 50℃，径向温差从 ~8℃ 收窄到 1~2℃；
之后：整体贴着 50℃（题给框架无潜热，会略高于真实物理的湿球温度——这是上一讲讨论过的框架局限，此处不重复展开）。
水分：
0~0.5 h：几乎不动（DD 只有 7×10⁻¹⁰，预热期特征时间 150 h 的结论仍然支配前期）；
0.5~3 h：温度上来后 DD 涨了 10 倍，表面水分开始明显下降，径向形成“外干内湿”的梯度——干燥前沿从表面向中心推进；
之后几十小时：中心水分被缓慢“拖”下来。中心是瓶颈，烘干时间由 C(0,t)C(0,t) 穿过 0.15 的时刻决定（这正是问题3）。
热-质耦合的体现：表面先干 → 表面 cp,kcp,k 下降、DD 上升 → 表面升温更快、失水更快，梯度自我强化。你的数值结果应能看到表面与中心的 DD 相差数倍——如果全场的 DD 几乎相同，说明耦合没生效，查代码。
4.5 输出要求核对表
另外，模型验证一节建议放一张全程总览图：横轴 0~72 h（对数或分段坐标），画出中心温度、表面温度、中心水分、表面水分四条曲线。问题3的“烘干结束时间”在这张图上就是一个交点——问题2和问题3用一张图打通，结构非常漂亮。
4.6 问题2小结
相对问题1，只多了三件事：
物性公式换成附录3（每步用冻结系数重算，注意 T 用开尔文）；
两方程联立（每步各解一个三对角方程组，可迭代一次）；
烘房工况外推（t>14400 s 取恒温 50℃ / 0.05，写明假设）。
求解器、网格、隐式格式、插值方法全部复用问题1。所以代码结构上应该把问题1的求解器写成函数 solve(T_air_func, C_air_func, prop_func, t_end, ...)，问题2、3、4 都只是换参数调用——这是工程上最重要的设计决策。
自测两题
若忘记把温度换成开尔文，数值结果会出现什么“荒谬特征”，你能据此反查出 bug 吗？（提示：C 永远 ≈ 2.55，表面也不干——因为 e−3850/50e−3850/50 约 10−3410−34）
为什么预测“3 h 末中心温度接近 50℃、但中心水分几乎没降到 0.15”不矛盾？（提示：热扩散率 α∼10−7α∼10−7 与水分 D∼10−9∼10−8D∼10−9∼10−8 差 1~2 个数量级——热走得快，水走得慢，这正是“预热平衡”存在的意义）