连续粉体进料与混合建模教程6失重式进料器建模① feed-factor 模型——有效堆积密度与螺杆输送版本声明块工具/软件PharmaPyPurdue CryPTSys2.0 官方已新增 feeder 单元模型为半经验 feed factor 模型本篇的feed_factor函数与本系列PowderFeeder类均为本系列自研示例非官方 API。语言/环境Python 3.9numpy2、matplotlib本篇只依赖 numpy 与 matplotlib。事实边界本篇只复现官方进料器方程的数学结构σ(t)、β、ρ_sat、ρ_in、V_sp 的精确定义与取值以 PharmaPy 官方代码原文为准本系列不编造其物理含义与数值。一句话结论本篇把第 5 篇骨架里的理想跟踪替换成 PharmaPy 2.0 官方进料器方程的结构——F_out(t) ff(t)·ω_screw(t)、ff(t) ρ_eff(t)·V_sp、ρ_eff(t) ρ_sat (ρ_in − ρ_sat)·exp(−σ(t)/β)、ω_imp(t) α·ω_screw(t)——并给出支持标量与数组的feed_factor(t, rho_in, rho_sat, sigma, beta, V_sp)实现当 σ/β 1 时有效堆积密度已完成约 63.2% 的趋饱和行程因 1 − e⁻¹ ≈ 0.632。〇、本篇要解决的认知问题feed factor进料因子这个词在厂商资料里被反复使用却很少有人把它的数学结构摆开看。这一节先钉五个问题PharmaPy 2.0 的进料器方程一共四式每一式在工程上分别回答哪个量由什么决定为什么有效堆积密度不是常数而要写成随某个状态量指数趋近饱和密度的形式这个函数形式在物理上意味着什么V_sp的量纲是什么它与螺距、螺纹深这些螺杆几何是什么关系为什么只能给量纲关系而不能给数值叶轮转速与螺杆转速为什么要写成比例关系助流这件事对流量到底起了什么作用单螺杆、双螺杆、无轴螺旋、振动式、皮带式该怎么选厂商标的 ±0.25% 与 ±3% 到底在比什么为什么容易差一个数量级一、机制解析1.1 官方方程的结构四式各管一件事PharmaPy 2.0 论文给出的进料器方程是半经验模型结构如下此处只复现结构参数语义与取值以官方代码原文为准F_out(t) ff(t) · ω_screw(t) (1) ff(t) ρ_eff(t) · V_sp (2) ρ_eff(t) ρ_sat (ρ_in − ρ_sat) · exp(−σ(t)/β) (3) ω_imp(t) α · ω_screw(t) (4)逐式的工程直觉一式一句话式数学形式工程直觉(1)流量 参数 × 转速质量流量正比于螺杆转速螺杆每转一圈送出一口物料转得越快送得越快(2)每转量 密度 × 体积每转输送的质量 物料的有效堆积密度 × 单转输送体积一口的大小由装得多实和勺子多大共同决定(3)密度随状态量指数趋饱和有效堆积密度随某种状态量 σ 指数趋近饱和密度状态量越大越接近装实的上限越小时越松(4)叶轮转速成比例叶轮impeller转速与螺杆转速成比例助流强度跟着螺杆走保证供料跟得上式 (3) 是全篇最值得琢磨的一式。它是一个一阶松弛形式把当前密度与饱和密度之差记作 Δ ρ_sat − ρ_eff则式 (3) 等价于 Δ(t) (ρ_sat − ρ_in)·exp(−σ(t)/β)——差值按指数衰减衰减速度由 σ(t)/β 决定。这正是一切趋向平衡过程的标准形状离目标越远靠拢的动力越大参数 β 是尺度σ(t) 是进度。必须钉死的边界本系列登记的不可断言项σ(t) 到底代表什么物理量、β 的单位与量级、ρ_sat 与 ρ_in 各自的定义、V_sp 是几何体积还是有效充填体积——这些在可获取的公开摘要中均未给出本系列一律不编造。本篇的教学实现只复现这个函数形状用来跑通结构与练手感把它对上真实物料必须以 PharmaPy 官方代码原文为准。1.2 有效堆积密度为什么不是常数粉体的堆积密度是个被严重低估的变量。同一批粉松装aerated与振实tapped状态下的堆积密度可以差出相当可观的比例料斗料位下降、充填方式改变、振动、含水率变化都会让此刻这堆粉装得实不实发生变化。式 (3) 把这个直觉写成了动力学有效堆积密度 ρ_eff 不是一个常数而是随某种状态量演化的变量。工程上有三个直接后果同一转速不等于同一流量。若 ρ_eff 会变式 (1) 里的 ff 就在变固定转速下流量随之漂移。这正是容积式螺杆给料精度只有 ±1%~±3%而失重式能到 ±0.25%~±1.0%的结构性原因。低速段的流量分辨率更差。转速降下来时每转那一口本身在变小若 ρ_eff 也在变控制灵敏度的衰减比线性预期更快——这为第 8 篇讲调节比 1:10 对控制权限的约束埋下伏笔。补料会把 ρ_eff 扰动一遍。补料瞬间料斗被快速填充堆积状态剧变ρ_eff 随之跳变这是第 7 篇补料期间流量尖峰如何在模型层面产生的机制来源。1.3 V_sp 的量纲关系与螺杆几何式 (2) 决定了V_sp的量纲必须是体积转m³/rev否则式 (1) 的右边量纲对不上升质量流量。把这个量纲与螺杆几何挂钩可以得到一个量纲层面的估算注意这只是几何直觉不是官方V_sp的定义单转输送体积几何估算 螺距 × 螺杆截面积 × 螺杆数 p × (π/4 · D²) × n_screws 单位 m³/rev忽略充填系数几何量进入哪个环节备注螺距 p单转输送体积每转前进的距离螺距越大每转送出的长度越长螺纹深单转输送体积截面积的有效部分与充填系数共同决定实际输送量螺杆直径 D单转输送体积截面积直径加倍截面积四倍螺杆数 n_screws单转输送体积的倍数双螺杆不止两倍还有啮合效应这份估算是上界它假设螺纹内被物料填满而真实螺杆内部有充填系数fill factor——第 5 篇的v_sp_teach_m3_per_rev就明确写了忽略充填系数。把它直接当流量用会系统性高估。正确做法是几何估算法只用来判断量级真正的V_sp要么由标定得到要么以官方定义为准。1.4 叶轮与螺杆为什么成比例式 (4) 说叶轮转速与螺杆转速成比例。这背后是一条很朴素的机械设计逻辑叶轮是助流的不是计量的。失重式进料器loss-in-weight feederLIW的结构里料斗底部有叶轮impeller负责打散、助流把粉体源源不断地喂给旋转螺杆螺杆再把它输出去。叶轮如果比螺杆慢会在螺杆入口形成供料不足螺杆就吃不饱式 (2) 的V_sp名存实亡叶轮如果比螺杆快很多则可能过度搅动、加剧扬尘与偏析。把它设计成与螺杆成比例比例系数 α是用一个操纵量螺杆转速同时控制两件事——计量与供料——的工程妥协。这也是选型时的隐藏判据如果你选的进料器在低转速段必须大幅提升叶轮转速才能保证供料它就不是为宽调节比设计的。1.5 失重式进料器的工作原理把结构放回硬件LIW 的闭环逻辑是四步多来源一致可交叉验证① 整台送料总成料斗 螺杆/振动送料器 驱动安装在精密称重传感器上 │ ② 物料卸出 → 系统连续监测重量下降速率 → 换算为瞬时质量流量 │ ③ 控制器把瞬时流量与目标质量流量比较 → 实时调整电机转速 │ ④ 于是自动补偿堆积密度、含水率、粒度分布的变化关键在②③④ 三步构成一个称重闭环它测的不是转速而是重量在掉多快所以物料性质变了它能自己找回来。这也是为什么 LIW 的精度±0.25%±1.0%能比纯容积式螺杆±1%±3%高一个量级——一个在测结果一个在猜结果。1.6 选型决策表进料器不是一种选错螺杆类型后面的 PID 调得再好也白搭。下表把五类 LIW 的适用材料与关键优势并列分类与厂商一手资料一致类型适用材料关键优势备注单螺杆 LIW自由流动粉体、颗粒最常用结构简单双螺杆 LIW黏聚、易架桥、易压实的粉体需要正位移输送时首选无轴螺旋 LIW黏、纤维状、黏聚材料轴不阻挡强制输送振动式 LIW易碎、易脆、磨蚀性材料零机械剪切适合高价值 API 微量给料皮带式需连续称重测流量的场合称重与输送解耦精度与规格对照全部为厂商转引口径引用时必须标等级指标数值等级LIW 连续给料精度±0.25% ~ ±1.0%相对设定值★ 厂商容积式螺杆给料精度±1% ~ ±3%相对设定值★ 厂商DSL G 连续进料误差典型 ±0.5%对设定值★ 厂商DSL G 批量精度±0.1–0.5%1 kg 批次★ 厂商调节比turndown1:10★ 厂商吞吐范围0.5–25,000 kg/h★ 厂商螺杆螺旋直径18–150 mm★ 厂商称重传感器class C3、不锈钢、IP68/IP69K★ 厂商载荷传感器采样率有的系统达16 次/秒★ 转引某研究数据采集1 Hz并滤波★ 转引铁律级提醒第 4 条报告进料精度时必须写明是对设定值还是对满量程——两者差一个数量级行业宣传口径常混用。本系列统一为相对设定值的连续进料误差。表里所有百分数都已按此口径标注看到别人的 ±0.5% 时先问一句对什么。二、完整代码与逐行剖析第一段代码给出feed_factor的通用实现支持标量与数组并跑一个数值演示。# -*- coding: utf-8 -*- 第 6 篇配套脚本一feed_factor 通用实现支持标量与数组。 本系列自研示例非 PharmaPy 官方 API。 结构ff(t) rho_eff(t) * V_sprho_eff(t) rho_sat (rho_in - rho_sat)*exp(-sigma/beta) 注意sigma、beta、rho_sat、rho_in、V_sp 的精确定义与取值以官方代码/原文为准。 量纲rho_eff [kg/m^3]V_sp [m^3/rev] - ff [kg/rev] ff [kg/rev] * omega_screw [rev/s] 质量流量 [kg/s]。 importnumpyasnpdeffeed_factor(t,rho_in,rho_sat,sigma,beta,V_sp):返回 feed factor ff(t)单位 kg/rev每转输送质量。 参数 ---- t : 时间秒或任意自变量标量或数组 rho_in : 式(3) 的起始密度量kg/m^3语义以官方原文为准 rho_sat : 式(3) 的饱和密度量kg/m^3语义以官方原文为准 sigma : 状态量 sigma(t)标量或数组语义与单位以官方原文为准 beta : 式(3) 的尺度参数语义与单位以官方原文为准 V_sp : 单转输送体积m^3/rev语义以官方原文为准 返回 ---- ff : 与 t / sigma 广播后的同形数组若为标量输入则返回 float t_arrnp.asarray(t,dtypefloat)sigma_arrnp.asarray(sigma,dtypefloat)ifbeta0:raiseValueError(beta 不能为 0式(3) 中 sigma/beta 会溢出)# 广播到同一形状使 t 与 sigma 都可为标量或数组sigma_b,_t_bnp.broadcast_arrays(sigma_arr,t_arr)# 式(3)有效堆积密度 饱和密度 (起始密度 - 饱和密度) * exp(-sigma/beta)rho_effrho_sat(rho_in-rho_sat)*np.exp(-sigma_b/beta)# 式(2)feed factor 有效堆积密度 * 单转输送体积ffrho_eff*V_spreturnfloat(ff)ifff.ndim0elseffdefmain():# 教学示例参数非文献值仅为跑通结构rho_in,rho_sat,beta350.0,600.0,20.0V_sp6.283e-6# m^3/rev约等于 D20 mm、螺距20 mm 的几何估算omega_screw0.5# rev/s即 30 rpm# --- 标量调用sigma 0、beta、2*beta 三个点 ---forsigmain(0.0,beta,2.0*beta):fffeed_factor(0.0,rho_in,rho_sat,sigma,beta,V_sp)print(fsigma{sigma:6.1f}ff{ff:.6e}kg/rev 流量{ff*omega_screw*3600:.3f}kg/h)# --- 数组调用整条 sigma 网格 ---sigma_gridnp.linspace(0.0,120.0,7)ff_arrfeed_factor(np.zeros_like(sigma_grid),rho_in,rho_sat,sigma_grid,beta,V_sp)frac(ff_arr-ff_arr[0])/(ff_arr[-1]-ff_arr[0])# 趋饱和行程占比fors,f,xinzip(sigma_grid,ff_arr,frac):print(fsigma{s:6.1f}ff{f:.6e}kg/rev 行程占比{x*100:5.1f}%)# 理论核对sigma/beta 1 时行程占比应为 1 - exp(-1) 0.6321print(f\nsigmabeta 时理论行程占比 {1-np.exp(-1):.4f})if__name____main__:main()第二段代码画两张图有效堆积密度随 σ 的演化曲线以及不同 β 下的敏感性对比。为便于独立运行此处重复给出feed_factor定义。# -*- coding: utf-8 -*- 第 6 篇配套脚本二有效堆积密度演化曲线与 beta 敏感性对比保存为 PNG。 自研教学脚本非官方 API。图中坐标轴用英文以免缺字注释用中文。 importnumpyasnpimportmatplotlib matplotlib.use(Agg)# 无显示环境下也能保存图片importmatplotlib.pyplotaspltdeffeed_factor(t,rho_in,rho_sat,sigma,beta,V_sp):与配套脚本一相同ff rho_eff * V_sprho_eff 由式(3)给出。t_arrnp.asarray(t,dtypefloat)sigma_arrnp.asarray(sigma,dtypefloat)ifbeta0:raiseValueError(beta 不能为 0)sigma_b,_np.broadcast_arrays(sigma_arr,t_arr)rho_effrho_sat(rho_in-rho_sat)*np.exp(-sigma_b/beta)ffrho_eff*V_spreturnfloat(ff)ifff.ndim0elseffdefrho_eff_curve(sigma,rho_in,rho_sat,beta):直接返回式(3)的 rho_eff便于画图。sigmanp.asarray(sigma,dtypefloat)returnrho_sat(rho_in-rho_sat)*np.exp(-sigma/beta)defmain():rho_in,rho_sat,V_sp350.0,600.0,6.283e-6sigmanp.linspace(0.0,150.0,400)betas[10.0,20.0,40.0]colors[tab:blue,tab:orange,tab:green]fig,axesplt.subplots(1,2,figsize(11,4.2))# --- 左图不同 beta 下 rho_eff 随 sigma 的演化 ---forb,cinzip(betas,colors):axes[0].plot(sigma,rho_eff_curve(sigma,rho_in,rho_sat,b),colorc,labelfbeta {b:g})axes[0].axhline(rho_sat,ls--,colorgray,lw1,labelrho_sat (saturated))axes[0].axhline(rho_in,ls:,colorgray,lw1,labelrho_in (initial))axes[0].set_xlabel(sigma (state variable))axes[0].set_ylabel(rho_eff [kg/m$^3$])axes[0].set_title(Effective bulk density vs sigma)axes[0].legend(fontsize8)axes[0].grid(alpha0.3)# --- 右图b20 时对应的 feed factor 与瞬时流量 ---b0,omega_screw20.0,0.5# rev/sfffeed_factor(np.zeros_like(sigma),rho_in,rho_sat,sigma,b0,V_sp)flow_kg_hff*omega_screw*3600.0axes[1].plot(sigma,flow_kg_h,colortab:red)axes[1].set_xlabel(sigma (state variable))axes[1].set_ylabel(mass flow [kg/h] (at 0.5 rev/s))axes[1].set_title(Resulting mass flow (beta 20))axes[1].grid(alpha0.3)fig.tight_layout()outch06_feed_factor.pngfig.savefig(out,dpi120)print(fsaved:{out})# 数值检查sigma 从 0 增到 5*beta 时还剩多少未走完的行程remainnp.exp(-5.0)print(fsigma 5*beta 时残余行程占比 {remain:.2%})if__name____main__:main()逐段剖析为什么feed_factor要支持数组因为第 7 篇的补料仿真里 σ 是一整条时间序列第 19 篇的敏感性扫描里 β 又是一整条网格。若函数只收 float调用方就得写外层循环既慢又容易把广播写错。让 numpy 广播在函数内部完成是这类纯函数最省事的设计。为什么先校验beta 0因为-sigma/beta会直接除零。更隐蔽的是beta 为负表达式不会报错但会把趋近饱和变成指数发散画出反物理曲线。建议在带标定的工程代码里加beta 0断言本篇为保持短小只拦零。为什么返回float而不是 0 维数组标量输入返回 numpy 0 维数组会让print结果多出array(...)外壳也会让json.dump直接失败。在边界处把 0 维数组转成 Python 标量能省掉下游一大堆类型处理。左图的rho_in与rho_sat两条参考线为什么必要因为式 (3) 的本质就是从 ρ_in 出发、向 ρ_sat 靠拢。把两根渐近线画出来才能一眼看出 β 的作用β 越小曲线拐得越急越早贴到饱和线β 越大越平缓。右图为什么不画无量纲量因为工程上要的是这个 β 对应多少 kg/h 的漂移。把最终流量画出来才能和精度指标相对设定值 ±0.25%~±1.0%直接对照——能对上指标的图才有决策价值。反直觉提醒式 (3) 里的 σ(t) 在所有公开摘要中都没有明确物理定义而它恰恰是有效密度会不会动、动多快的唯一驱动量。很多二次开发项目死在这一点上开发者默认 σ 是常数于是 ρ_eff 是常数、ff 是常数模型在稳态附近看着还行一进补料与低料位区就全线偏离。正确姿势是把 σ 当状态量让它在第 7 篇的料斗模型里被料位驱动。三、常见报错与排查1现象RuntimeWarning: overflow encountered in exp随后ff变成inf或nan。根因sigma/beta过大如 β 取得极小或 σ 量纲与 β 不匹配。exp(-大数)本身不溢出但若写成exp(sigma/beta)或 β 为负指数会正向发散。解法检查 β 的符号与量纲对sigma/beta做区间断言例如-700 sigma/beta 700确认 σ 与 β 同单位——单位不匹配是这里最常见的真凶。2现象ValueError: beta 不能为 0。根因参数未初始化或标定结果为 0。解法β 是式 (3) 的尺度参数为 0 时函数退化为阶跃。要么给出非零初值要么在拟合时加下界约束。3现象数组调用时报ValueError: operands could not be broadcast together。根因t与sigma形状不兼容如长度 400 与长度 7 且都不为 1。解法确保两者形状一致或把其中一个写成标量/长度 1 数组。第 7 篇里t与sigma都是同长度时序注意别把t_grid与sigma_at_refill短数组混传。4现象算出来的ff数量级明显不对差 10⁶。根因V_sp的单位错了。若用 mm³/rev 代入结果会差 10⁻⁹ 倍若用 cm³/rev差 10⁻⁶ 倍。解法全系列统一 SI 派生单位密度 kg/m³、体积 m³/rev、转速 rev/s、流量 kg/s对外展示再换算 kg/h。在注释里写清量纲链比事后 debug 便宜得多。5现象别人质疑你这个 σ/β 是不是编的。根因把教学实现的参数说成了官方定义。σ(t)、β、ρ_sat、ρ_in、V_sp 的精确语义与取值未在公开摘要中给出。解法在所有出现处标注参数语义以 PharmaPy 官方代码原文为准此处为教学示意引用官方方程时只引结构不替官方定义参数。四、动手练习练习 1基础跑通函数并核对三个数值。运行配套脚本一。成功标准量化sigmabeta那行的行程占比 63.2%与理论值 1 − e⁻¹ 相对偏差 1e-6sigma2*beta时行程占比 86.5%1 − e⁻²sigma5*beta时残余行程占比 0.67%e⁻⁵脚本二会打印。练习 2进阶验证 β 的尺度含义。固定rho_in350、rho_sat600、V_sp6.283e-6对beta ∈ {10, 20, 40}分别求rho_eff达到500 kg/m³所需的 σ。成功标准量化三个 σ 值之比接近1:2:4相对偏差 2%并写一句物理解释σ/β 是无量纲进度β 越大需要越大的 σ 才能走到同一进度。练习 3挑战把 feed factor 接回第 5 篇的进料器并做量纲链自检。把feed_factor接进PowderFeeder第 5 篇定义的类的 volumetric 分支用 σ 从 0 线性增到 120、β20、V_sp取几何估算值screw_rpm固定 30 rpm模拟 60 s。成功标准量化输出 60 s 内mass_flow_kg_h的起点与终点同时打印量纲链rho_eff × V_sp → kg/rev、× rev/s → kg/s、× 3600 → kg/h并断言终点流量 / 起点流量落在1.5–2.0之间对应 ρ_eff 从 350 趋近 600 但未走完。五、小结与下一篇预告本篇把 PharmaPy 2.0 官方进料器方程的结构拆开看了式 (1) 说流量正比于转速式 (2) 说每转量等于有效密度乘单转体积式 (3) 说有效堆积密度随状态量指数趋近饱和式 (4) 说叶轮转速跟着螺杆走。feed_factor的实现只有六行但量纲链与广播处理决定了它能不能被后续篇目大量复用。选型部分给出了五类 LIW 的适用材料与厂商口径的精度对照表——请记住那条铁律报精度必须写清对设定值还是对满量程。事实边界再强调一次σ(t)、β、ρ_sat、ρ_in、V_sp 的精确定义与取值本篇没有给出也不允许编造一律以 PharmaPy 官方代码原文为准feed_factor与本系列PowderFeeder均为自研示例。下一篇第 7 篇把 σ 的来历补上建料斗质量平衡 ODEdM/dt −F_out补料时叠加F_refill把补料refill事件接进来写清补料期间临时切到 volumetric 模式、维持最后一次已知转速这一行业实践以及料位下降→堆积密度变化→feed factor 漂移如何与式 (3) 呼应代码给hopper_mass_balance(M0, t_grid, F_out, refill_events)并用scipy.integrate.solve_ivp与手写显式欧拉两种积分方式做对比。第 8 篇再给螺杆转速装上 PID 闭环与抗积分饱和。本篇认知问题回显FAQQ1PharmaPy 2.0 的进料器方程有哪几式每式在工程上管什么A共四式——流量 feed factor × 螺杆转速feed factor 有效堆积密度 × 单转输送体积有效堆积密度随状态量 σ 指数趋近饱和密度叶轮转速与螺杆转速成比例分别对应计量、每转口大小、密度演化与助流供料。Q2为什么有效堆积密度要写成随状态量指数趋近饱和的形式A因为当前密度与饱和密度之差按指数衰减是一切趋平衡过程的标准形状差值越大靠拢越快σ/β 是无量纲进度当 σ/β 1 时有效密度已完成 1 − e⁻¹ ≈ 63.2% 的趋饱和行程。Q3V_sp 是什么量纲与螺距和螺纹深是什么关系AV_sp必须是体积/转m³/rev否则流量式量纲对不上它可用螺距 × 螺杆截面积 × 螺杆数作几何上界估算实际还受充填系数影响其精确定义以 PharmaPy 官方代码或原文为准。Q4叶轮转速为什么写成与螺杆转速成比例A叶轮负责打散助流、把粉体喂给螺杆属供料而非计量环节让它跟随螺杆转速成固定倍数就能用螺杆转速这一个操纵量同时管住计量与供料避免螺杆入口供料不足。Q5单螺杆、双螺杆、无轴螺旋与振动式进料器该怎么选精度怎样比较A单螺杆适合自由流动粉体与颗粒双螺杆适合黏聚、易架桥、需正位移的粉体无轴螺旋适合黏或纤维状物料振动式零机械剪切、适合易碎料与高价值 API 微量给料失重式精度约 ±0.25%~±1.0%容积式螺杆约 ±1%~±3%。