6.1 问题4在问什么？

前面三个问题都把药材当成半径固定 2 cm 的圆柱。但附件2给了实测数据：药材在烘干过程中会缩水——半径从 2 cm 缩到约 1.2 cm，缩掉 40%！

问题4要求：把这个收缩考虑进模型，重新确定烘干时长。

用附录4的经验公式（物性公式换了一组系数），加上附件2的半径数据 R(t)R(t)，重新解整个烘干问题。

表面上看“只是换个参数”，实际上有一个根本性的困难：

计算区域本身在随时间变化：一开始方程定义在 0≤r≤20≤r≤2 cm，后来定义在 0≤r≤1.20≤r≤1.2 cm。你的网格点（0.1 cm 间隔）会不断有“点跑到边界外面去”，差分格式直接崩掉。

这就是移动边界问题（Stefan 类问题的简化版）。解决办法是数学里最常用的套路之一：换个坐标，让移动的边界“看起来不动”。

6.2 核心技巧：相对坐标变换

① 定义新坐标

x=rR(t)x=R(t)r​

xx 是“相对位置”：x=0x=0 永远是中心，x=1x=1 永远是当前时刻的表面——不管表面实际在哪。

这样定义的函数：

u(x,t)≡C(x R(t),  t)u(x,t)≡C(xR(t),t)

读作：“相对坐标 x 处的水分浓度”。它和原来的 C(r,t)C(r,t) 描述的是同一件事，只是换了个“视角”。温度同理定义 θ(x,t)=T(xR(t),t)θ(x,t)=T(xR(t),t)。

好处立竿见影：新方程永远定义在固定的 0≤x≤10≤x≤1 上，网格切 21 个点不动，边界永远在 x=1x=1，再也不会有网格点“跑出去”。

② 代价：方程会多出一项

天下没有免费的午餐。坐标换着容易，但“网格在动”这个物理事实不会消失，它会以附加项的形式回到方程里。下面把推导一步步做出来——这是本讲唯一的硬数学，每一步都给你展开。

6.3 完整推导（跟着链式法则走）

我们要把原方程

∂C∂t=1r∂∂r(D r ∂C∂r)∂t∂C​=r1​∂r∂​(Dr∂r∂C​)

翻译成 u(x,t)u(x,t) 的方程。翻译规则就是链式法则。

第一步：空间导数（简单）

r=xR(t)r=xR(t)，固定 tt 时 RR 是常数，所以：

∂C∂r=∂u∂x⋅∂x∂r=1R∂u∂x∂r∂C​=∂x∂u​⋅∂r∂x​=R1​∂x∂u​

导数被“放大”了 1/R1/R 倍——几何上 obvious：同样的浓度变化发生在更小的半径里，梯度当然更陡。

第二步：时间导数（关键一步，注意陷阱）

求 ∂C/∂t∂C/∂t 时，rr 是固定的。但 x=r/R(t)x=r/R(t) 里 RR 在变——固定 rr 时，xx 自己也在随时间变！所以链式法则有两项：

∂C∂t∣r=∂u∂x⋅∂x∂t∣r⏟因为 x 在动+∂u∂t∣x⏟真正的局部变化∂t∂C​​r​=因为 x 在动∂x∂u​⋅∂t∂x​​r​​​+真正的局部变化∂t∂u​​x​​​

先算 ∂x∂t∣r∂t∂x​​r​：

x=rR(t)  ⇒  ∂x∂t∣r=−r R˙R2=−x R˙Rx=R(t)r​⇒∂t∂x​​r​=−R2rR˙​=−RxR˙​

其中 R˙=dR/dtR˙=dR/dt 是收缩速率（负值，药材在变小）。代回去：

∂C∂t∣r=∂u∂t∣x−R˙ xR ∂u∂x∂t∂C​​r​=∂t∂u​​x​−RR˙x​∂x∂u​

这一项就是“新数学”的全部：多出来的 −R˙xR∂u∂x−RR˙x​∂x∂u​ 叫对流附加项。

它的物理意义非常直观，值得单独说：

药材收缩时，每个物质点被“拖”着往中心移动。你站在固定相对位置 xx 上观察，会发现浓度在变化——不是因为水扩散了，而是因为你脚下的地面在移动（边界向你压过来）。附加项描述的就是这个“地面移动”效应。

类比：你在匀速行驶的火车上往前扔球，地面上的人看到球速 = 车速 + 你扔的速度。R˙R˙ 就是“车速”。

第三步：拼装完整方程

把第一、二步的结果代回原方程。右边的扩散算子：∂∂r=1R∂∂x∂r∂​=R1​∂x∂​，1r=1xRr1​=xR1​，两层嵌套后（每层出一个 1/R1/R，外面再除 xRxR，共 1/R21/R2）：

1r∂∂r(Dr∂C∂r)=1xR2∂∂x(D x ∂u∂x)r1​∂r∂​(Dr∂r∂C​)=xR21​∂x∂​(Dx∂x∂u​)

于是得到移动坐标系下的水分方程：

  ∂u∂t=1xR2∂∂x(D x ∂u∂x)+R˙ xR ∂u∂x  ∂t∂u​=xR21​∂x∂​(Dx∂x∂u​)+RR˙x​∂x∂u​​

（我调整了符号写成 +R˙xR+RR˙x​：R˙<0R˙<0，所以这一项实际效果是把高浓度的中心值往里“拉”。）

温度方程完全同构——做一模一样的变换：

∂θ∂t=1ρcp⋅1xR2∂∂x(k x ∂θ∂x)+R˙ xR ∂θ∂x∂t∂θ​=ρcp​1​⋅xR21​∂x∂​(kx∂x∂θ​)+RR˙x​∂x∂θ​

结构解读（帮你记忆）

∂u∂t⏟局部变化=1xR2∂∂x(Dx∂u∂x)⏟扩散（同老方程，多了1/R2）+R˙xR∂u∂x⏟收缩对流局部变化∂t∂u​​​=扩散（同老方程，多了1/R2）xR21​∂x∂​(Dx∂x∂u​)​​+收缩对流RR˙x​∂x∂u​​​

• 扩散项：和原方程长得一样，只是 rr 换成 xx，且系数多了 1/R21/R2（区域越小，扩散“跨过全场”越快，RR 在分母合理）；

• 对流项：唯一的 newcomer，系数 R˙/RR˙/R 量级多大？看附件2：前期 1800 s 内 RR 从 2.0 → 1.873，R˙≈−7×10−5R˙≈−7×10−5 cm/s，R˙/R≈−3.5×10−4 s−1R˙/R≈−3.5×10−4 s−1；后期几乎为零。而扩散项量级 D/R2∼10−9/(0.02)2≈2.5×10−6D/R2∼10−9/(0.02)2≈2.5×10−6……等等，看似对流项 10−410−4 反而更大？不对——对流项是 R˙xR∂u∂xRR˙x​∂x∂u​，还要乘梯度 ∂u/∂x∂u/∂x。做个正经的比较：对流项与扩散项之比 ∼∣R˙∣RD≈7×10−7×0.025×10−9≈3∼D∣R˙∣R​≈5×10−97×10−7×0.02​≈3——收缩项和扩散项同量级，绝不能扔掉！（这正是出题人专门设置问题4的原因：这道题里收缩不是小修正，是主效应之一。）注意这个估算用了 R˙≈7×10−7R˙≈7×10−7 m/s（前期收缩最快的阶段），可以说“前期收缩效应显著，后期 RR 趋稳后衰减”。把这个量级比较写进论文，直接说明问题4的必要性——非常好的论证素材。

6.4 定解条件也要搬家

初始条件：u(x,0)=2.55u(x,0)=2.55，θ(x,0)=28θ(x,0)=28——不变。

中心：∂u∂x∣x=0=0∂x∂u​​x=0​=0（对称性在新坐标下保持）——不变。

表面 x=1x=1：原条件是

−D∂C∂r∣r=R=hm[CR−Cair]−D∂r∂C​​r=R​=hm​[CR​−Cair​]

翻译过来：∂C∂r=1R∂u∂x∂r∂C​=R1​∂x∂u​，所以

  −DR∂u∂x∣x=1=hm[u(1,t)−Cair(t)]  −RD​∂x∂u​​x=1​=hm​[u(1,t)−Cair​(t)]​

比原来的条件多除了一个 RR。这个细节最容易漏：收缩后表面面积变小、梯度变陡，综合效果通过这个 1/R1/R 体现。温度边界同理：−kR∂θ∂x∣x=1=h[θ(1,t)−Tair]−Rk​∂x∂θ​​x=1​=h[θ(1,t)−Tair​]。

对流传质总量的变化：跑掉的水分 = hm(CR−Cair)×hm​(CR​−Cair​)× 当前表面积 2πRL2πRL。注意两点：

1. 表面积在变小（RR 缩小）→ 后期蒸发总面积下降，这是收缩“拖慢干燥”的一面；

2. 严格写能量/质量守恒时，被排出的水分体积 ≈ 收缩的体积（药材失水所以变细）——附件2的 R(t)R(t) 数据本质上就是累计失水量的记录。可以用它做一个漂亮的质量守恒校验：

药材当前含水量初始含水量=∫0R(t)C⋅2πrL dr+骨架质量∫0R0C0⋅2πrL dr+骨架质量初始含水量药材当前含水量​=∫0R0​​C0​⋅2πrLdr+骨架质量∫0R(t)​C⋅2πrLdr+骨架质量​

用附录4的 ρ=760+90Cρ=760+90C 可以从半径反推平均含水率，与数值解的体积平均含水率对比——两条独立途径算同一个量，对上了就是模型正确性的强证据。这是问题4最有含金量的验证，强烈建议做。（这里有一个建模自由度：R(t)R(t) 直接取附件2插值，还是从失水量守恒反推？建议主线用附件2插值 + 守恒校验，改进讨论里可以提“由失水动力学反推 R(t)R(t)”的封闭模型。）

6.5 R˙R˙ 怎么取？——数值细节

方程里需要 R˙(t)=dR/dtR˙(t)=dR/dt，但附件2只给离散点（1800 s 间隔）。三个方案：

重要观察：附件2数据到 259200 s（72 h）时 R≈1.198R≈1.198 cm 且几乎不再变化。这说明：

• 收缩主要发生在前 10 小时（2.0 → 1.3），之后迅速趋缓；

• 收缩的“终点”约 1.2 cm 而不是 0——因为干药材骨架还在（干基含水率 0.15 ≠ 0，仍有少量水）。

如果按问题3算出的烘干时间超过 72 h，附件2数据就不够用了——需要用拟合方案 B 外推（R→R∞≈1.2R→R∞​≈1.2），或者用守恒关系延拓。这是问题4特有的数据边界问题，务必检查你的烘干时长和附件2的时间覆盖范围是否匹配，并在论文中说明处理方法。

6.6 数值求解：改动比想象的小

好消息：求解框架完全复用。改动清单：

1. 网格：固定 x∈[0,1]x∈[0,1]，21 个点（对应输出“0.1 cm”的要求——注意！见 6.7 的讨论）；

2. 每步组装方程组时：扩散项的系数从 DiriΔr2 ri±1/2ri​Δr2Di​​ri±1/2​ 换成 DixiR2Δx2 xi±1/2xi​R2Δx2Di​​xi±1/2​，RR 取当前值；

3. 新增对流项：R˙ xiRui+1−ui−12ΔxRR˙xi​​2Δxui+1​−ui−1​​（中心差分）——把它作为已知系数项加入三对角方程组的右端或系数矩阵。隐式处理扩散项 + 显式处理对流项（对流系数用当前时刻值）在这类问题上足够稳定，因为 R˙R˙ 变化平缓；

4. 边界条件：表面式里多除一个 R(t)R(t)；

5. 物性：用附录4的公式（ρ=760+90Cρ=760+90C 等，D=4.2×10−4e−0.30Ce−3850/TD=4.2×10−4e−0.30Ce−3850/T——注意 DD 的系数比附录3小很多：4.2×10−44.2×10−4 vs 2.4×10−32.4×10−3，但 e−0.30Ce−0.30C 衰减更慢。量级自检：C=2、T=323K 时 D≈4.2×10−4×0.55×6.6×10−6≈1.5×10−9D≈4.2×10−4×0.55×6.6×10−6≈1.5×10−9 m²/s，比问题2/3 的 DD 小约 3 倍——提示问题4的烘干时间会比问题3长，这是新经验公式的效果，不是模型错误）；

6. 停止条件：与问题3相同，max⁡xu(x,t)<0.15maxx​u(x,t)<0.15。

稳定性再检查：对流项显式处理的稳定性要求 Courant 数 ∣R˙∣RRΔtΔxR∣R˙∣R​ΔxΔt​……准确说以 xx 坐标计的对流速度是 ∣R˙∣x/R≤∣R˙∣/R∣R˙∣x/R≤∣R˙∣/R，Courant 数 =∣R˙∣RΔtΔx≈3.5×10−4×11/20=7×10−3⋘1=R∣R˙∣​ΔxΔt​≈3.5×10−4×1/201​=7×10−3⋘1，安全得不能再安全。扩散项仍是隐式，无条件稳定。论文里可以加一句这个检验。

6.7 输出格式的坑：表6的列头变了！

仔细读题：

“在论文中按表6的格式给出每隔 6 h、到药材中心距离每隔 0.5 cm 的水分浓度”，表6列头是“0, 0.5, „ , 药材表面”。

和表5对比，表6有两个变化：

1. 最后一列叫“药材表面”而不是 2 cm——因为表面在动！最早表面在 2 cm，后来缩到 1.2 cm，“表面”是一个随时间的量。所以最后一列填的是 u(1,t)u(1,t)，即每个时刻当前表面的浓度；

2. 中间列的物理意义也要想清楚：0.5 cm、1 cm 这些是绝对距离还是相对位置？题目写“到药材中心距离每隔 0.5 cm”——是绝对距离。那么这些固定绝对距离的点，在新坐标里的位置是 x=0.5/R(t)x=0.5/R(t)，随时间变！这就麻烦了：我们的解存在固定 xx 网格上，要取固定 rr 的值就得插值。

   两种处理：

   • 方案甲（推荐）：解完后，对每个输出时刻 tt，把固定 rr = 0, 0.5, 1.0, 1.5 cm 换算成 x=r/R(t)x=r/R(t)，在 u(x,t)u(x,t) 上插值取值；若 r>R(t)r>R(t)（该绝对距离已在表面之外，比如后期 R=1.2R=1.2 cm 时 1.5 cm 处已不存在），该格填“—”或取表面值并加注。这是更忠实题意的做法；

   • 方案乙：直接输出固定 xx = 0, 0.25, 0.5, 0.75, 1 的值，并在论文中说明“距离按相对半径给出”。省事但可能被扣题意分。

   建议：论文表格用方案甲，并在表注里写明“绝对距离超过当前半径的位置以表面值计”。

3. result4.xlsx：“每隔 60 s 时间、0.1 cm 距离”——同样面临绝对/相对的解释问题。同样用方案甲：固定绝对距离 0, 0.1, …，超出当前半径取表面值。或者更稳妥地，在说明文档里写清楚你的取法。这类格式歧义，主动声明你的处理方式就不丢分，闷头不说不声明才丢分。

烘干结束时间的判定同问题3：max⁡xu<0.15maxx​u<0.15，由中心点（x=0x=0）决定。

6.8 问题4预期结果的定性画像 + 与问题3的对比

跑完后论文里最值钱的一张表就是问题3 vs 问题4 对比：

注意附录4的 DD 系数明显小于附录3，所以 t4>t3t4​>t3​ 的主要贡献很可能来自经验公式切换，而收缩本身的净效应（加速 vs 拖慢的竞争）才是“问题4模型”想要回答的科学问题。论文里要把这两个效应拆开讨论（可以先算“附录4公式 + 固定半径”的中间版本做归因），这一段归因分析是问题4的精华，能写出来说明你真正理解了模型，而不是跑了代码交差。

6.9 问题4检查清单

四个问题总览（全系列收官）

一条暗线贯穿始终：每一讲我们都先做量级估算、再跑数值、再拿结果对照预估。这套“估算 → 计算 → 验证”的循环，比任何单点技巧都重要——它让你在 300 h 和 3 h 的荒谬结果面前不掉坑里，也让论文的每一张表都自带可信度。

自测两题

1. 如果药材是吸湿膨胀（R˙>0R˙>0）而不是收缩，附加对流项的方向如何改变？它会把中心的水“推”出去还是“拉”住？（提示：R˙>0R˙>0 时边界往外跑，物质点被拖向外——看对流项符号）

2. 附件2末尾 RR 稳定在 1.2 cm 而不是 0，如何用“干基含水率 0.15 + 骨架密度”粗算验证这个收缩终点是否合理？（提示：湿药材总质量 ∝ ρ(C0)R02ρ(C0​)R02​，干骨架质量 ∝ ρ(0.15)R∞2ρ(0.15)R∞2​，解出 R∞/R0R∞​/R0​ 与 1.2/2.0 = 0.6 对比）

收官。