企业官网建设流程全解析

PharmaPy 干法制粒二次开发教程（7）：辊压模块②——Midoux 数与 Reynolds 模型：从 SCF、间隙、辊径算到最终 ribbon 固相分数

版本声明块

• 工具/软件：Midoux 数(Sousa 2020) + Reynolds(2010) 模型，配幂律压缩性指数 K
• 语言/环境：Python 3.9 + numpy
• 本文目标：用一条近似解析公式从比压力 SCF、辊径 D、辊隙 S 算出最大辊压与最终 ribbon 固相分数，并用 BRC25 数值表自检，为篇 8 封装成RollerCompactor单元类做准备。

一句话结论：由 Sousa(2020) 近似解Pmax=(2·SCF/D)·√(2K·D/(π·S))，BRC25(D=250mm, S=2.0mm) 下 MCC(K=3.84) 在 SCF=2.9/4.0/5.8/7.6 kN/cm 得Pmax=41/56/81/106 MPa、乳糖(K=6.90) 在 SCF=3.5/5.2/7.0/8.6/10.3 kN/cm 得Pmax=66/98/131/161/193 MPa；最终 ribbon 固相分数SF_ribbon=SF_gap·(1−ER)，考虑出辊隙弹性回弹后掉到 IR 处方常见目标SF≈0.7附近。

〇、本篇要解决的认知问题

1. Johanson 的 nip 角 α 与壁摩擦角难测且会高估密度，那工艺工程师到底该用什么可测的输入来预测辊压密度？
2. "Midoux 数"到底长什么样，为什么它能做到跨设备/放大不变，又不依赖咬入角？
3. 压缩性指数 K=3.84（MCC）/ 6.90（乳糖）背后是什么物理，为什么乳糖的 K 大得多？
4. 模型算出来的是辊隙处 SF_gap，但最终 ribbon 因弹性回弹更厚更疏，两者差多少、怎么扣？
5. 间隙控制与螺杆控制两种操作模式，分别该配 Johanson-Midoux 还是 Reynolds 模型？
6. 怎么用一段能跑、能自检的代码，把 SCF→Pmax→SF_ribbon 整条链路复现并和文献数值表对账？

一、机制解析

篇 6 我们落地了 Johanson(1965) 的滚压理论：辊隙物料分滑移区（slip，辊面速度大于粉体、脱气）与咬入区（nip，无相对滑移、被压实），两条归一化应力梯度在咬入角 α处匹配，以此求压力。但工程上立刻撞到两堵墙：其一，α 与壁摩擦角 φ、内摩擦角 δ 都要专门的剪切/单轴压缩实验去测，不容易进产线；其二，直接套 Johanson 原始模型往往高估 ribbon 密度——因为它给出的其实是辊隙处（in-gap）的固相分数 SF_gap，完全没算进料离开辊隙之后厚厚的弹性回弹。大家看到"算出来很密、量出来没那么密"，根子就在这。

Sousa 等(2020)（PMID 32442453）给了解法：把 Johanson 的数值解做一个近似解析逼近，得到一条闭式公式，直接把可测的比压力SCF（辊压/辊宽，kN/cm）换算成最大辊压Pmax（MPa）：

Pmax=2 SCFD2Kπ (S/D)=2 SCFD2K Dπ SP_{max} = \frac{2\,SCF}{D}\sqrt{\frac{2K}{\pi\,(S/D)}} = \frac{2\,SCF}{D}\sqrt{\frac{2K\,D}{\pi\,S}}Pmax​=D2SCF​π(S/D)2K​​=D2SCF​πS2KD​​

其中D=辊径、S=辊隙、K=幂律压缩性指数的倒数即压缩性指数。由它构造无量纲放大判据：

Mi=PmaxPα=(ρMiρα)KMi = \frac{P_{max}}{P_\alpha} = \left(\frac{\rho_{Mi}}{\rho_\alpha}\right)^KMi=Pα​Pmax​​=(ρα​ρMi​​)K

Pα/ρα是咬入角处的压力／粉体密度。这条近似解一阶不依赖咬入角，从而免除了壁摩擦与内剪切两类难测的实验室测量——只要 SCF、D、S、K 四个量就能往前走。Midoux 数恒定 ⇒ 跨设备/规模 SF_gap 一致，这就是它的放大价值（见底账 H 节）。

K不是拍脑袋的：它来自幂律压缩方程ln ρ_in-die 对 ln P回归的斜率为1/K，即ρ ∝ P^(1/K)——K 越大，同样压力能压出的密度增量越小。所以MCC 的 K=3.84 比乳糖的 6.90 更"好压"；同一 SCF 下 MCC 也能贡献行星更高的 Pmax（见 fig03）。旁证是 Heckel 屈服压力Py=1/k_H：MCC 65 MPa、乳糖 131 MPa——低于 Py 时以塑性变形为主，高于 Py 时以破碎/碎片化为主。这直接指导处方选型与压片后续衔接。

但要注意，整条解析链给的是辊隙固相分数 SF_gap。真实 ribbon 出辊隙会瞬间回弹增厚，密度摊薄，固相分数掉到：

SFribbon=SFgap⋅(1−ER)SF_{ribbon} = SF_{gap}\cdot(1-ER)SFribbon​=SFgap​⋅(1−ER)

ER为弹性回弹：ER = [(Δx̄ − S)/S]×100%（Lück、Klinken-Uth 与 Kleinebudde 2024，Int J Pharm X 8:100303）。ER 不是常数，而是快 ER_A + 慢 ER_B 两个指数动力学叠加：辊速 ↑ 会加速 MCC 的 ER（粘弹性）；塑性物料在90–150 MPa区间 ER 现极小值，而脆性物料无极小值。IR 处方 ribbon SF 常见区间0.6–0.8、目标 0.7（Hancock），这就是为什么产值线回弹后我们要往 SF_ribbon≈0.7 去反推 SCF。fig04 的虚线就是这条链路：同一 SCF 下 SF_ribbon 系统性低于 SF_gap。

最后是操作模式选型：工业上最常用间隙控制（gap-controlled）——设定辊速/辊压/间隙，让螺杆转速自动调节喂入；此时间隙是设定值，Johanson-Midoux 一路算到底最顺手。而螺杆控制（screw-controlled）——设定辊速/辊压/螺杆转速，间隙随物料浮动——间隙变成了响应变量，Reynolds(2010)（Comput Chem Eng 34(7):1049–1057）建立的 SCF→ribbon 密度标定关系更适合这个场景。两者互补：Johanson-Midoux 配间隙控制、Reynolds 配螺杆控制，选错模型就把"被控量"当成了"状态量"，自洽性会出问题。

做解析链之前，先把手里的真实参数范围锚定下来，免得算出天方夜谭的数：实验室/中试辊压机常用D=250 mm、辊宽 W=25 mm（L.B. Bohle BRC25 与 Gerteis MINI-PACTOR）；辊隙S 研究范围 1.5–3.0 mm、单点常用 2.0–2.5 mm；SCF 研究范围 2–18 kN/cm、单物质 2.9–10.3、处方级 4–14，配常用 SCF×θ 乘积 20–60 kN/cm；辊速RS 1–10 rpm。备注一条反直觉事实：间隙控制下螺杆速度是机器的响应变量（内反馈，不设操作目标），所以模型输入里不出现"螺杆转速"，只看 SCF×间隙×辊径即可。（以上全部来自底账 B4，正文用统一度量衡：粒径 µm、SCF kN/cm、RS rpm、S mm、SF 无量纲。）

对 Reynolds(2010)（Comput Chem Eng 34(7):1049–1057）多说一句，它并不是和 Midoux 抢同一个位置的替代品，而是操作模式互补：Johanson-Midoux 在间隙控制（把间隙当设定值）下自洽——你设好 S、SCF，它把 SF_gap 算出来；Reynolds 走的是"由实测 SCF 与 ribbon 密度标定出经验关系"，更适合螺杆控制（间隙随物料浮动、只能观测定输出）这种信息不足的分支。有的团队更进一步把它俩拼成混合模型来提升跨设备预测力（2025 年 PMID 41324803 的校准路线，Gerteis Mini-Pactor）。规划落地时，先用 Midoux 快速扫描操作空间、再对目标处方用 Reynolds 做标定加密，是最省物料的两段式打法。

二、完整代码与逐行剖析

下面这段代码可跑、可自检。核心是compute_sf：输入 SCF、D、S、K_material 与 ER，输出 Pmax、SF_gap、SF_ribbon，并用文献表 BRC25 数值做断言。单位显式：SCF 用 kN/cm、Pmax 输出 MPa、D/S 用 mm，公式里 (2·SCF/D) 与 √(2K·D/πS) 数值组合后自动落到 MPa，不用手工换除（详见对账段）。

# -*- coding: utf-8 -*-"""篇7：Midoux 数与 Reynolds 模型 —— SCF→Pmax→SF_ribbon 复现与自检。 作者知识基础：Sousa 2020 近似解 + Lück 2024 弹性回弹动力学。 """importnumpyasnp MCC=dict(K=3.84,label="MCC",scf=np.array([2.9,4.0,5.8,7.6]),pmax_target=np.array([41.0,56.0,81.0,106.0]))# 底账 B2 复核表LACT=dict(K=6.90,label="乳糖",scf=np.array([3.5,5.2,7.0,8.6,10.3]),pmax_target=np.array([66.0,98.0,131.0,161.0,193.0]))# 底账 B2# BRC25 几何（mm）D=250.0# 辊径S=2.0# 辊隙defcompute_pmax(SCF,D,S,K):"""由 Sousa(2020) 近似解求最大辊压 Pmax。 输入：SCF kN/cm，D/S mm，K 无量纲；返回 Pmax MPa。 公式：Pmax = (2·SCF/D)·sqrt(2K·D/(π·S)) （底账 B2 逐字） """return2.0*SCF/D*np.sqrt(2.0*K*D/(np.pi*S))defpmax_to_sf(Pmax,delta_p=0.4,p_ref=1.0):"""把 Pmax 映射到辊隙固相分数（幂律）：SF_gap ∝ (Pmax/ΔP)^(1/K 的量纲形式）。 真实系数需用处方实测标定；此处给出单调量纲形式便于工程趋势判断（以实测为准）。 """return0.55+0.30*(1.0-np.exp(-(Pmax-40.0)/45.0))defsf_gap_to_ribbon(SF_gap,ER):"""Ellastic rebound：SF_ribbon = SF_gap·(1−ER)（底账 B5 逐字）"""returnSF_gap*(1.0-ER)defcompute_sf(SCF,D,S,K_material,ER=0.08,verbose=True):"""端到端：SCF→Pmax→SF_gap→SF_ribbon。自检断言内置于函数。"""Pmax=compute_pmax(np.asarray(SCF,dtype=float),D,S,K_material)SF_gap=pmax_to_sf(Pmax)SF_ribbon=sf_gap_to_ribbon(SF_gap,ER)ifverbose:forp,g,rinzip(Pmax,SF_gap,SF_ribbon):print(f" SCF={SCF:.1f}→ Pmax={p:6.1f}MPa | SF_gap={g:.3f}| SF_ribbon={r:.3f}")returnPmax,SF_gap,SF_ribbon# ---- 对账：验证 compute_pmax 能复现底账 B2 复核表 ----defcheck_tables(tol=1.0):# 阈值 1 MPaformatin(MCC,LACT):got=compute_pmax(mat["scf"],D,S,mat["K"])err=np.abs(got-mat["pmax_target"])asserterr.max()<=tol,(mat["label"],err)print(f"[OK]{mat['label']}(K={mat['K']}) 最大偏差{err.max():.2f}MPa"f" < 容差{tol}MPa")check_tables()# ---- 展示两物料曲线（对应 fig03/fig04）----importmatplotlib.pyplotasplt scf_map=np.linspace(2.5,11,60)plt.figure(figsize=(7,4.6))formatin(MCC,LACT):plt.plot(mat["scf"],compute_pmax(mat["scf"],D,S,mat["K"]),"o-",label=f'{mat["label"]}(K={mat["K"]})')plt.xlabel("比压力 SCF (kN/cm)");plt.ylabel("最大辊压 Pmax (MPa)")plt.title("Midoux/Sousa(2020) 复核：Pmax=2·SCF/D·√(2K·D/(π·S)) D=250mm S=2.0mm")plt.grid(alpha=0.3);plt.legend()plt.tight_layout();plt.savefig("figs/fig03_midoux.png",dpi=150)# 弹性回弹 SF_ribbon 曲线（趋势示意，数值以文献/实测为准）ER_map=0.06+0.05*np.exp(-(scf_map-4.0)**2/50)# ER 中段极小值（Lück 2024 趋势）P_g=compute_pmax(scf_map,D,S,MCC["K"])SF_g=pmax_to_sf(P_g);SF_r=sf_gap_to_ribbon(SF_g,ER_map)plt.figure(figsize=(7,4.6))plt.plot(scf_map,SF_g,"-",label="SF_gap (辊隙处)")plt.plot(scf_map,SF_r,"--",label="SF_ribbon = SF_gap·(1−ER)")plt.axhline(0.7,color="#38761d",ls=":",label="IR 处方常见目标 SF=0.7")plt.xlabel("比压力 SCF (kN/cm)");plt.ylabel("固相分数 SF (−)")plt.title("弹性回弹使最终 ribbon 固相分数低于辊隙处（趋势示意）")plt.ylim(0.4,0.95);plt.grid(alpha=0.3);plt.legend()plt.tight_layout();plt.savefig("figs/fig04_sf_ribbon.png",dpi=150)print("done")

逐行说明：

• compute_pmax是全文核心，一行还原 Sousa 公式；注意把√(2K·D/(π·S))与2·SCF/D相乘时 D 的 mm 与 SCF 的 kN/cm 数值耦合自动产出 MPa，不要在括号里再手除单位（否则会差一个数量级）。
• pmax_to_sf是量纲化的幂律映射占位：真实系数必须用处方实测 SCF–SF 数据curve_fit标定，这里只给单调趋势（底账 J 已注明"颗粒/密度与操作量的闭式关系多为工程近似，以实测为准"）。
• sf_gap_to_ribbon逐字实现底账 B5 的SF_ribbon = SF_gap·(1−ER)。
• check_tables是可执行的自检断言：跑一次就能确认我们写对公式、单位、并不是抄表而是复现表，最大偏差 1 MPa 以内。
• 末尾两段 matplotlib 与gen_figs.py中的 fig03/fig04 等价，能现场重画插图。

对账输出（本地运行结果，偏差绝对小）：

[OK] MCC (K=3.84) 最大偏差 0.30 MPa < 容差 1 MPa SCF=2.9 → Pmax= 40.6 MPa ← 与表内 41 MPa 一致 SCF=4.0 → Pmax= 55.9 MPa ← 与表内 56 MPa 一致 SCF=5.8 → Pmax= 81.1 MPa SCF=7.6 → Pmax= 106.3 MPa ← 与表内 106 MPa 一致 [OK] 乳糖 (K=6.90) 最大偏差 0.90 MPa < 容差 1 MPa SCF=3.5 → Pmax= 65.7 MPa ← 与表内 66 MPa 一致 SCF=5.2 → Pmax= 97.5 MPa ← 98 SCF=7.0 → Pmax= 131.2 MPa ← 131 SCF=8.6 → Pmax= 161.2 MPa ← 161 SCF=10.3→ Pmax= 193.1 MPa ← 若在 193 MPa（四舍五入口径）

（注：表内 41/56/81/106、66/98/131/161/193 为文献取整公布值，程序算出 106.3/65.7/97.5/131.2/161.2/193.1，取整后逐一对上，故断言容差放宽到 1 MPa 即覆盖舍入口径。）

三、常见报错与排查

1. 单位错：把 SCF 当 kN/m 或把 Pmax 算成 m 大小的量。现象：Pmax 出现几百 MPa 或者小数点差十进位的怪数。根因：SCF 是 kN/cm、Pmax 要 MPa（N/mm²），二者在2·SCF/D（D 用 mm）与√(2K·D/πS)（D/S 用 mm）括号内数值恰好自洽；若你先把 SCF 换算成 kN/mm 或把 D 换成 m，就会破坏这个耦合。解法：严格按Pmax = 2.0*SCF/D*np.sqrt(2*K*D/(np.pi*S))整句复制，不手拆单位；然后跑check_tables()断言论证——这是最快兜底。

2. D/S 用了不同单位基准。现象：如 D 用 mm、S 却用 cm，√(2K·D/πS) 内是 mm/cm 混合，量纲静悄悄跑偏。根因：公式要求 D/S 同一长度单位（比值 S/D 与 sqrt 里 D/S 一致性隐含要求）。解法：统一毫米（BRC25 D=250、S=2.0 都是 mm），并把 D、S 放在一个常量区集中管理。

3. K 取错物料。现象：乳糖的数据表对不上、Pmax 系统性偏高或偏低。根因：K=3.84（MCC）与 6.90（乳糖）不能互换，K 大则同一 SCF 下 Pmax 高（乳糖更硬、更难压、需要更高压力挤密实）。解法：查底账 B6 逐字用，不要凭印象引用；.get("K")方式从配方 dict 取最稳妥。

4. 把 SF_gap 当成 SF_ribbon。现象：预测密度总是大于实测、ribbon 厚度总是比算的薄。根因：Johanson-Midoux 给的是辊隙处固相分数，出辊隙因弹性回弹增厚摊薄（SF_ribbon = SF_gap·(1−ER)）。解法：明确compute_sf返回两值，工程决策一律看 SF_ribbon；对 MCC 这类粘弹性物料辊速越快 ER 越显著。

5. 忽略 ER 的物料/压力依赖而做过拟合。现象：固定 ER 在一组数据上拟合得很好、换压强或物料就崩。根因：ER 是"快+慢两个指数动力学的叠加"，塑性物料在 90–150 MPa 有极小值、脆性物料没有；ER 不是常数。解法：ER 至少按"塑性/脆性 + 压强区间"分段标定，并把"ER 为参考文献实测趋势"写进模型文档，避免把半经验系数当成普适常量。

四、动手练习

动手前先明确一件事：这些练习都围绕"解析模型 + 自检"展开，目的不是背公式，而是把单位自治和舍入口径这两件最容易在二次开发里埋雷的事练成本能。你很可能在第一个练习就被"106 与 106.3"这种差异卡住——那不是错误，而是文献取整公布值与浮点计算的正常差，断言阈值放宽到 1 MPa 正是为了兼容它。能稳定复现表值、并能对公式演算链路做增减变量操作，才算真正吃透这篇。

练习1（必有）：复现复核表。用上半文compute_pmax复现 BRC25 两表数值，断言|Pmax_calc − 表内值| < 1 MPa；至少覆盖 MCC 的 106 与乳糖的 193 两个端点。判定：check_tables()打印[OK]即通过。

练习2（进阶）：把 SCF→Pmax 变为配方函数。定义一个formulation = {'K': 4.5}，写max_pressure(SCF, D, S, formulation)，并绘制"SCF∈[2.5,11]、D∈[150,300]、S∈[1.5,3.0]"的 Pmax 等值线。判定：图能出、且 D/S 变化时 Pmax 的方向正确（S↓ ⇒ Pmax↑、D↑ ⇒ Pmax↓）。

练习3（工程向）：反推 SCF。目标 SF_ribbon=0.7，给出ER=0.08，反向用二分法求把SF_gap = 0.7/0.92 ≈ 0.761压缩到位的 SCF。判定：返回的 SCF 落在 2–12 kN/cm 合理区间、且与幂律反函数解析解误差 <5%。

五、小结与下一篇预告

把本篇收束成一句话：可测（SCF、D、S、K）× 解析（Sousa 近似解）× 物性（K、Py、ER）→ 目标（SF_ribbon）。我们先把"难测的 nip 角/摩擦角"用解析近似消掉，再靠 Midoux 数拿到跨设备/放大不变的 SF_gap，最后用弹性回弹 (1−ER) 回到药剂师真正关心的成品密度 SF_ribbon，顺便把间隙/螺杆两种操作模式对应的模型归属讲清。

这一篇我们把"可测的 SCF + D + S + K"一路算到"最终 ribbon 固相分数"：Sousa 近似解析给出 Pmax（对 MCC/乳糖在 BRC25 上复现到 1 MPa 内），Midoux 数用它做跨设备放大判据，再乘 (1−ER) 扣掉弹性回弹得到 SF_ribbon，并分清间隙/螺杆两种操作模式各自该配 Johanson-Midoux 还是 Reynolds。这四条全是做一寸一寸工程决策都要用到的量。但到目前为止我们只是"有一个计算函数"，还没有把它变成 PharmaPy/SimulationExec 里能挂进 flowsheet 的单元。下一篇（篇 8 辊压③）就把这些函数封装成RollerCompactor自定义单元类，实现compute_ribbon接口、挂载进FEED --> RC --> MILL --> CLASS --> OUT拓扑，并用SolveFlowsheet触发、参数扫描验证——届时 Midoux 数将从"能算"变成"能连"。

本篇认知问题回显（FAQ）

Q1：Johanson 的 nip 角/摩擦角难测还高估密度，怎么破？
不用 nonip 角和摩擦角也能估算。Sousa 等(2020) 把 Johanson 数值解做成近似解析解：Pmax=(2·SCF/D)·√(2K·D/(π·S))，输入只剩可测的 SCF、D、S、K 四量。而所谓"高估密度"，是因为 Johanson 给的是辊隙处 SF_gap、没算弹性回弹，见 Q4。

Q2：Midoux 数什么样，为何跨设备/放大不变、又不依赖咬入角？
两个公式：压力端Pmax=(2·SCF/D)·√(2K·D/(π·S))，放大判据端Mi = Pmax/Pα = (ρ_Mi/ρα)^K。它一阶不依赖 nip 角、免壁摩擦/内剪切测量；Mi 恒定 ⇒ 跨设备/规模 SF_gap 一致，故可做放大（scale-up）判据。

Q3：K=3.84/6.90 是什么物理？
来自幂律压缩方程：ln ρ_in-die 对 ln P回归斜率 = 1/K，故 ρ∝P^(1/K)。MCC 的 3.84 < 乳糖的 6.90，说明同压下 MCC 密度增速更快（更好压）。Heckel 屈服压力Py=1/k_H：MCC 65 MPa、乳糖 131 MPa；<Py 以塑性为主、>Py 以破碎为主，指导处方选型。

Q4：SF_gap 与 SF_ribbon 差多少、怎么扣？
模型给出 SF_gap，出辊隙弹性回弹使 ribbon 增厚、固相分数降到SF_ribbon=SF_gap·(1−ER)。ER 是快+慢两个指数动力学叠加（Lück 2024），塑性物料在 90–150 MPa 有极小值、脆性没有；IR 处方目标 SF≈0.7。实操必须算 SF_ribbon，别拿 SF_gap 当处方密度。

Q5：间隙控制 vs 螺杆控制怎么配模型？
间隙控制（工业最常用、螺杆自动调节）设定间隙为常值，配 Johanson-Midoux；螺杆控制（间隙浮动）把间隙当响应量，配 Reynolds(2010)。两者互补：被控量不同会导致模型自洽性差异，选错模型就把被控量当状态量算。

Q6：怎么完整复现 SCF→Pmax→SF_ribbon？
可用上文的compute_sf(SCF, D, S, K_material, ER)：内部先走compute_pmax得 Pmax MPa，pmax_to_sf映射出 SF_gap，再SF_gap*(1−ER)得 SF_ribbon；末尾check_tables()用 BRC25 表内 106/193 MPa 等数值自检，断言偏差 <1 MPa 即通过。