📢 แจ้งเตือน: qingfengxd/article/details/

A 题药材的烘干问题 — 题目分析报告

目标竞赛：2026 高教社杯全国大学生数学建模竞赛（CUMCM 2026）
题目：A 题药材的烘干问题

一、问题背景与目标

热风烘干是中药材加工的关键工序，包含预热平衡与恒温干燥两个阶段。题目要求建立药材内部温度场 $T (r, t)$ 与水分浓度场 $C (r, t)$ （干基含水率，kg/kg）的数学模型，并完成四个子问题：

子问题	任务	关键输出
问题 1	建立预热平衡阶段温湿度变化模型（常物性）	表 1、表 2；result1.xlsx
问题 2	建立整个烘干过程（两阶段）模型（变物性，附录 3 公式）	表 3、表 4；result2.xlsx
问题 3	由问题 2 模型确定烘干所需时间（各处 $C < 0.15$ ）	表 5；result3.xlsx
问题 4	考虑尺寸收缩（附件 2 半径数据，附录 4 公式）确定烘干时长	表 6；result4.xlsx

二、几何与坐标

药材近似为圆柱体，长 $L = 25$ cm，半径 $R = 2$ cm。由于 $\gg R$ （长度是直径的 6 倍以上），轴向与径向的耦合较弱，且题目要求的结果均以到药材中心（轴线）的距离 $r∈[0,R]r\in[0,R]$ 给出，故采用无限长圆柱的轴对称一维径向模型。

空间坐标： $r∈[0,R]r\in[0,R]$ ， $r = 0$ 为药材中心（轴线）， $r = R$ 为药材表面。
问题 4 中 $R$ 随含水率流失而收缩， $R = R (t)$ （附件 2 给出），构成移动边界问题。

三、控制方程

设 $T (r, t)$ （°C 或 K）为温度， $C (r, t)$ （kg/kg 干基）为水分浓度。

能量守恒（径向导热）

$r∂T∂r)\rho c_p \frac{\partial T}{\partial t}=\frac{1}{r}\frac{\partial}{\partial r}\!\left(k\,r\frac{\partial T}{\partial r}\right)$

质量守恒（径向水分扩散）

$r∂C∂r)\frac{\partial C}{\partial t}=\frac{1}{r}\frac{\partial}{\partial r}\!\left(D\,r\frac{\partial C}{\partial r}\right)$

其中 $ρ\rho$ 密度、 $c_p$ 比热容、 $k$ 导热系数、 $D$ 水分扩散系数；问题 1 取常数，问题 2–4 取附录给出的经验公式（随 $C, T$ 变化，构成非线性抛物型方程组）。

四、边界条件与初始条件

对称条件（ $r = 0$ ）

$∂T∂r∣r=0=0,∂C∂r∣r=0=0\left.\frac{\partial T}{\partial r}\right|_{r=0}=0,\qquad \left.\frac{\partial C}{\partial r}\right|_{r=0}=0$

第三类边界（ $r = R$ ，对流换热与对流传质）

$−k∂T∂r∣r=R=h(T∣R−T∞(t)),−D∂C∂r∣r=R=km(C∣R−C∞(t))-k\left.\frac{\partial T}{\partial r}\right|_{r=R}=h\left(T\big|_{R}-T_\infty(t)\right),\qquad -D\left.\frac{\partial C}{\partial r}\right|_{r=R}=k_m\left(C\big|_R-C_\infty(t)\right)$

其中 $W/(m2⋅K)h=25\ \mathrm{W/(m^2\cdot K)}$ 对流换热系数， $m/sk_m=8\times10^{-7}\ \mathrm{m/s}$ 对流传质系数， $C∞(t)T_\infty(t),\,C_\infty(t)$ 为烘房（环境）温度与水分浓度，由附件 1 给定。

初始条件

$kg/kgT(r,0)=28\ ^\circ\mathrm{C},\qquad C(r,0)=2.55\ \mathrm{kg/kg}$

五、物性参数

问题 1（常物性，附录 2）

$W/(m⋅K)\rho=820\ \mathrm{kg/m^3},\quad c_p=2600\ \mathrm{J/(kg\cdot K)},\quad k=0.36\ \mathrm{W/(m\cdot K)}$
$m/sh=25\ \mathrm{W/(m^2\cdot K)},\quad k_m=8\times10^{-7}\ \mathrm{m/s}$
$为干基含水率)D=7\times10^{-9}\,e^{-0.89/C}\quad(\mathrm{m^2/s},\ C\ \text{为干基含水率})$

注：经 PDF 版面几何核验（分式横线上下方字符定位），三处扩散系数公式的水分指数均为分式 $a / C$ ，而非乘积 $a C$ 。

问题 2、3（变物性，附录 3）

$ρ=650+128C,cp=1450+2736CC+1,k=0.21+0.38CC+1\rho=650+128C,\qquad c_p=1450+\frac{2736C}{C+1},\qquad k=0.21+\frac{0.38C}{C+1}$
$K)D=2.4\times10^{-3}\,e^{-0.45/C}\,e^{-3850/T}\quad(\mathrm{m^2/s},\ T\ \text{取 K})$

问题 4（变物性 + 收缩，附录 4）

$ρ=760+90C,cp=1850+2150CC+1,k=0.12+0.20CC+1\rho=760+90C,\qquad c_p=1850+\frac{2150C}{C+1},\qquad k=0.12+\frac{0.20C}{C+1}$
$e−3850/T(m2/s)D=4.2\times10^{-4}\,e^{-0.30/C}\,e^{-3850/T}\quad(\mathrm{m^2/s})$

六、烘房环境条件（附件 1）与外推假设

附件 1 给出 $0≤t≤144000\le t\le 14400$ s（4 h）内烘房温度与水分浓度，采样间隔 60 s：

时刻	环境温度 $T∞T_\infty$	环境水分浓度 $C∞C_\infty$
0 s	28.0 °C	0.01963 kg/kg
1200 s	37.9 °C	0.02891 kg/kg
3600 s	47.5 °C	0.04272 kg/kg
6000 s	49.7 °C	0.04841 kg/kg
14400 s	50.2 °C	0.04986 kg/kg

可见预热平衡阶段为 $0∼0\sim$ 约 6000 s：环境温度由 28 °C 快速升至约 50 °C 并趋于平台；环境湿度因药材放湿同步升高至约 0.05 kg/kg。

外推假设（已记录依据）：附件 1 仅覆盖前 4 h，而问题 2–4 需覆盖 2–3 天的全过程。假设：

$0≤t≤144000\le t\le 14400$ s：直接采用附件 1 的 $T∞(t),C∞(t)T_\infty(t),C_\infty(t)$ （线性插值）；
$t > 14400$ s（恒温干燥阶段）：环境温度稳定在平台值 $∘CT_\infty=50\ ^\circ\mathrm{C}$ ；环境水分浓度维持末值 $C∞=0.0499C_\infty=0.0499$ kg/kg。

依据：附件 1 显示两级参数在 $t≈6000t\approx 6000$ s 后已进入平台，恒温阶段的"参数不同"体现为环境条件由上升转为恒定，而物性公式按题目要求统一采用附录 3 / 附录 4。

七、移动边界处理（问题 4）

问题 4 的求解域 $[0, R (t)]$ 随时间收缩。令 $ξ=r/R(t)∈[0,1]\xi=r/R(t)\in[0,1]$ 为归一化径向坐标，由均匀收缩假设，材料点在运动中保持 $ξ\xi$ 不变，即 $ξ\xi$ 就是随材料运动的参考坐标。

在该参考坐标下，控制方程保持为纯扩散形式，收缩的影响完全体现在度量因子 $1/R^2(t)$ 上：

$⁣(ξ∂T∂ξ)\frac{\partial T}{\partial t}=\frac{k}{\rho c_p R^2(t)}\frac{1}{\xi}\frac{\partial}{\partial\xi}\!\left(\xi\frac{\partial T}{\partial\xi}\right)$
$⁣(ξ∂C∂ξ)\frac{\partial C}{\partial t}=\frac{D}{R^2(t)}\frac{1}{\xi}\frac{\partial}{\partial\xi}\!\left(\xi\frac{\partial C}{\partial\xi}\right)$

表面边界（ $ξ=1\xi=1$ ）相应变为 $∂ξT=h(T−T∞)-k/R\,\partial_\xi T=h(T-T_\infty)$ ， $∂ξC=km(C−C∞)-D/R\,\partial_\xi C=k_m(C-C_\infty)$ 。 $R (t)$ 由附件 2 的半径数据三次样条插值得到。

关于"平流项"的说明（口径统一）：若把收缩问题直接以空间坐标 $r$ 表述并作坐标变换，形式上会出现附加的平流项 $ξR˙/R⋅∂ξ(⋅)\xi\dot R/R\cdot\partial_\xi(\cdot)$ 。但如前所述，均匀收缩下 $ξ\xi$ 恰为材料参考坐标，在该参考系中物理上不存在相对干物质的净对流，故本文不引入附加平流项，直接求解上述纯扩散形式。程序实现中通过令 $R˙\dot R$ 项不生效（dR_func=None）来落实这一口径，与论文正文表述一致。

数值对照（用于说明该处理的影响幅度）：若改为在空间坐标下保留附加平流项求解，问题 4 的烘干时间为 187167 s（约 51.99 h），与本文采用的参考坐标纯扩散结果 181944 s（约 50.54 h）相差约 2.9%。该差异已在灵敏度讨论的误差范围内，不改变烘干时长的量级结论。

八、数值方法

采用有限体积法在 $N$ 个控制体上离散（径向网格均匀， $r$ 网格间距 0.1 cm，即 21 个节点覆盖 $0∼20\sim2$ cm），时间方向采用隐式（后退 Euler）— 非线性系数取上一层迭代值的格式，空间离散得到三对角线性方程组，用 Thomas 算法求解。

通量形式保证柱坐标下 $r = 0$ 处对称条件的自然处理（ $r$ 加权的面通量在轴心趋于零）；
采用隐式格式使 $D$ 变化范围大（约 $m2/s10^{-10}\sim10^{-8}\ \mathrm{m^2/s}$ ）时仍可稳定取 $Δt=1\Delta t=1$ s；
网格与时间步长通过时间步/空间步减半的收敛性检验验证（见编程手报告）。

输出网格：论文表格按题目指定距离（0, 0.5, 1, 1.5, 2 cm）与时间点取值；result 文件按每秒、每 0.1 cm 输出完整结果。

九、四个子问题的建模要点

问题 1：常物性预热模型， $t∈[0,1800]t\in[0,1800]$ s，输出表 1/表 2 及 result1.xlsx。
问题 2：变物性两阶段全过程模型， $t∈[0,3]t\in[0,3]$ h 输出表 3/表 4，全过程 $1$ s 分辨率写入 result2.xlsx。
问题 3：以问题 2 模型推进至各点 $C < 0.15$ 全部满足，取该时刻为烘干时长；输出表 5 与 result3.xlsx。
问题 4：移动边界 + 附录 4 物性，以附件 2 半径轨迹驱动，求烘干时长；输出表 6 与 result4.xlsx。

口径与边界说明（防歧义）：

独立模型体系计数：问题 1 = 1 套（常物性预热模型）；问题 2、3 共用 1 套（变物性两阶段全过程模型）；问题 4 = 1 套（移动边界变物性模型）。均满足"每子问题最多两个独立模型体系"约束，灵敏度/误差检验/可视化不另计。
result2.xlsx 输出时段：问题 2 的模型覆盖"整个烘干过程"（预热 + 恒温两阶段），故 result2.xlsx 保存全过程每秒、每 0.1 cm 的完整结果，论文表 3/表 4 仅摘录前 3 h；全过程时长与问题 3 判定的烘干时长一致。
$R (t)$ 超出附件 2 覆盖范围的处理：附件 2 覆盖 $0≤t≤2592000\le t\le 259200$ s（72 h），此后半径已平台于 1.198 cm。若模拟时长超过 72 h，取 $R (t) = R (259200) = 1.198$ cm、 $R˙=0\dot R=0$ 恒定外推（对结果影响可忽略）。
时间网格：取 $Δt=1\Delta t=1$ s 与题目要求的输出分辨率一致，避免输出插值引入额外误差。

十、模型评价与验证方案

验证方案：

网格无关性： $Δr\Delta r$ 减半、 $Δt\Delta t$ 减半，检查关键点温度/水分浓度变化 $10^{-3}$ ；
守恒性检验：全柱内水分总量减少量应等于表面累计传质量（数值积分闭合误差检验）；
极限检验：令 $D→∞D\to\infty$ 时浓度场应趋于均匀并与集中参数模型一致；令 $h,km→0h,k_m\to0$ 时应退化为绝热/绝质（初值保持）情形；
对比文献/解析：常物性绝热圆柱的解析级数解可用于问题 1 边界简化情形的量级校验。