PharmaPy 干法制粒二次开发教程6辊压模块①——Johanson 滚压理论与咬入角版本声明块工具/软件Johanson(1965) 滚压理论干法制粒辊压建模基础语言/环境Python 3.9 numpy/scipy本文目标讲透两区滑移/咬入机理与咬入角 α 的物理含义写一个求 α、画归一化应力梯度匹配的数值实现并诚实呈现该模型高估 ribbon 密度、未计弹性回弹的边界。数值来源纪律本篇所有概念、公式、数值一律来自本系列事实底账_调研素材参考.md的 B1/B4 条目凡原文逐字公式未能在底账核实的均在正文/图/码中明确标注以 Johanson(1965) 原文/文献为准。咬入角 26°→9°、临界 17° 等数值均出自底账 B1Al-Asady 转引。一句话结论Johanson(1965) 把对辊间的粉床按辊面线速度与粉体速度的关系分为滑移区slip辊面快、脱气与咬入区nip无相对滑移、被压实咬入角 α 就是两区交界角且大部分压实发生在咬入区α 由滑移区与咬入区两条归一化应力梯度在 α 处匹配唯一确定。模型忽略辊速准静态假设直接给出的是辊隙处固相分数 SF_gap因未计入间隙之后的弹性回弹ER而系统性高估最终 ribbon 密度即有SF_ribbon SF_gap——这一局限正是下一篇引入 Midoux 数与 Reynolds 修正的伏笔。〇、本篇要解决的认知问题为什么对辊里看似一直在压机制上却要拆成滑移区和咬入区两段咬入区凭什么说大部分压实发生在那里咬入角 α 到底是怎么从数学上求出来的为什么是两条归一化应力梯度曲线相交而不是别的判据为什么模型允许忽略辊速和螺旋转速这对二次开发意味着什么约束为什么直接套用 Johanson 会高估ribbon 密度SF_gap和SF_ribbon差在哪一步咬入角 26°→9°、临界 17° 这些数字在什么条件下测/推导的我在自己的设备上能直接复用吗一、机制解析1.1 对辊到底在做什么从挤压到两区辊压机roller compactor的本质是两台相对旋转的压辊把疏松粉体卷进去、压成致密的料带ribbon。想象两个圆柱相向转动粉体在辊隙上方堆积被带进越来越窄的辊隙最终在最小间隙辊隙 S处被压成最密的状态。听起来全流程都在压但对辊粉床的力学结构并不均一——Johanson(1965) 的理论贡献之一就是把这个全在压的朴素画面精确切成了两个力学性质截然不同的区滑移区slip位于最外侧辊面线速度大于粉体的整体运动速度粉体相对辊面打滑。这个区的主要功能是脱气——松散粉体中残存的空气在这里被挤出去同时粉体开始被轻微压缩并积聚。因为相对辊面滑动摩擦力主导应力增长比较平缓。咬入区nip位于辊隙最窄段辊面线速度与粉体速度锁定、无相对滑移粉体被咬住并一同被刚性压实。大部分压实密度提升恰恰发生在这段——只有到了咬入区粉体才在两点之间真正承压被压密压力随角度快速攀升。把压实主要集中在咬入区这句话记牢它是后面模型只描述咬入区压力-密度关系的原因也是为什么咬入角 α 是全局最关键的几何量——α 一旦确定咬入区的起始位置就定了压缩压力/密度的积分区间也就定了。此外压缩应力方向有一个值得注意的约定Johanson 把压实应力看作沿水平方向垂直于垂直进料压力 P₀即水平推力把粉体从辊隙上方往下压入咬入区。图figs/fig02_johanson.png直观画出了这两区以咬入角为界的几何位置读者应把外→内 滑移区→咬入区、界角叫 α这个空间感钉进脑子里。1.2 咬入角两条归一化应力梯度的相交点比起画图Johanson 最漂亮的地方在如何定量求 α。Théo滑移区里粉体相对辊面滑移其压力沿角度 θ 的增长受壁摩擦角 φ支配咬入区里粉体被锁死压实压力增长受粉体本体内摩擦角 δ更准确说受粉体可压缩性与塑性力学支配。于是两类区的压力随角度的相对变化率即归一化应力梯度1σdσdθ\frac{1}{\sigma}\frac{d\sigma}{d\theta}σ1​dθdσ​各自是 θ 的函数一条由壁摩擦 φ 驱动、一条由内摩擦 δ 驱动两条梯度曲线在某一个角度上数值相等——这个交点Johanson 定义为咬入角 α。换句话说α 的定义性条件是一条梯度连续性方程(1σdσdθ)slip∣θα (1σdσdθ)nip∣θα \left.\left(\frac{1}{\sigma}\frac{d\sigma}{d\theta}\right)_{\rm slip}\right|_{\theta\alpha} \;\; \left.\left(\frac{1}{\sigma}\frac{d\sigma}{d\theta}\right)_{\rm nip}\right|_{\theta\alpha}(σ1​dθdσ​)slip​​θα​(σ1​dθdσ​)nip​​θα​它在物理上说的是跨界处压力不能出现斜率突变——如果滑移区的压力增长突然跳变为咬入区的更快增长那分界点处压力/密度计算会撕裂这不符合粉体是连续介质的基本假设。所以两区在数学上必须光滑衔接而 α 就是那个让衔接成立的角度。需要诚实指出底账 B1Johanson 原文求 α 的那两条逐字导数表达式以及滑移区与咬入区的具体归一化应力梯度公式在本系列底账中未能逐字核实见底账 J未尽核实项。因此本篇给出的实现采用结构 数值求根的方式——我们把两条梯度在 α 处匹配这一结构忠实落地用一组物理合理的梯度函数形式滑移区随壁摩擦 φ、咬入区随内摩擦 δ 单调变化但明确声明梯度函数逐字形式以 Johanson(1965) 原文/权威文献为准这里的版本是工程实现用于复现原理与理解行为不冒充原文。这符合本系列真实性铁律 6数值来源标注。1.3 模型输入与准静态假设Johanson 模型的输入底账 B1/B4粉体压缩曲线/可压缩性即压力→密度关系压缩性指数 K下一篇用到内摩擦角 δ、壁摩擦角 φ辊几何辊径 D、辊宽 W、辊隙 S目标辊压比压力 SCF。而更重要的纪律是忽略什么Johanson(1965) 采用准静态假设——忽略辊速与螺旋转速。也就是说在你写代码、跑仿真时Johanson 这一层根本不接收辊速/螺旋转速作为参数压力-密度关系被当成与时间无关的纯几何摩擦平衡问题。这既是优点解析、省算也是边界下一篇 Reynolds 模型正是为螺杆/辊速也成为影响因素而存在两者互补Johanson 适合间隙控制模式、Reynolds 适合螺杆控制模式。对二次开发的心智模型提示Johanson 层是空间/几何层它回答应力沿角度怎么分布、α 多大操作-响应的层留给 Midoux/Reynolds。1.4 局限SF_gap ≠ SF_ribbon底账 B1 用一个小节专门记录局限重要这也是专家诚实的一部分nip 角与摩擦角难测δ、φ 是力学单元固有权参数实验室直接测很麻烦工程上常靠历史数据/压片压缩数据反推误差不可忽视直接套用会高估 ribbon 密度Dec 2003、Bi 2014 转引——理由见下一条模型输出的是在辊隙gap位置计算的固相分数SF_gap它描述辊隙最窄处、仍受约束时粉体可达的密实程度。但 ribbon 一旦流出间隙构件不再受辊面约束弹性应力释放、体积回弹弹性回弹 ERelastic recovery密度下降。所以链条是KaTeX parse error: Expected EOF, got _ at position 86: …n}}_{\text{ SF_̲gap − ER 的净结果}}…即模型给出的SF_gap是上限/乐观值最终成型 ribbon 的固相分数SF_ribbon必然低于它。这正是为什么篇 7 要引入 Midoux 无量纲数做解析放大、并用 Reynolds 与弹性回弹 ER 修正把SF_gap换算到SF_ribbon——Johanson 本身止步于SF_gap缺口留给下一篇补。1.5 咬入角典型数值与可调杠杆底账 B1 给了两组可引用的关键数值制药临界咬入角一般 17°并且提高壁摩擦可增大咬入角壁摩擦 φ 大脱气滑移段更长、咬入更充分这是处方/设备设计时可用的杠杆MCC 恒定压力/间隙下辊速从 3 rpm 提到 7 rpm实测咬入角从 26° 降到 9°Al-Asady 数据转引——注意这条 26°→9° 是实测咬入角随辊速的响应它揭示了一个 Johanson 准静态模型并不内置的现象辊速实际会改变咬入角快辊压实不如慢辊充分。这与 1.3 的准静态忽略辊速是不同层面的事实建模上 Johanson 理论上不显式含辊速但现实测量里辊速是有影响的。二次开发时要有这种模型假设 vs 实测响应的分辨力。需要强调这些数值是特定设备/物料MCC的实测/推导结论不是普适常数。工程上要按底账纪律标出来源均可溯源底账 B1不可当成万能默认。数值复现时把17°、26°→9°当作数量级与行为趋势参考而非直接塞进自己的仿真。二、完整代码与逐行剖析下面实现 Jessica/工程师版的两区应力梯度匹配 数值求根。声明在代码注释中再强调一次逐字梯度公式以 Johanson(1965) 原文/文献为准本脚本是结构正确、数值求根的工程实现用于复现梯度匹配求 α的原理与行为。importnumpyasnpfromscipy.optimizeimportbrentqimportmatplotlib.pyplotasplt# ---------------------------------------------------------------# 0. 模型输入底账 B4/B1单位约定m, rad# D辊径, S辊隙, W辊宽W 本篇只做派生需无量纲时不直接用# δ内摩擦角, φ壁摩擦角单位弧度度→弧度换算见下# ---------------------------------------------------------------D_roll250e-3# 辊径 250 mmB4L.B. Bohle BRC25 / Gerteis MINI-PACTORS_gap2.5e-3# 辊隙 2.5 mmB4常用 2.0–2.5 mmdelta_deg40.0# 内摩擦角度示例phi_deg15.0# 壁摩擦角度示例deltanp.radians(delta_deg)phinp.radians(phi_deg)# ---------------------------------------------------------------# 1. 两区归一化应力梯度的结构实现# 声明真实性纪律Johanson 原文逐字梯度表达式未在本系列事实# 底账核对到见 _调研素材参考.md 条目 J故以下用物理合理的形式# - 滑移区梯度 G_slip(θ)受壁摩擦 φ 驱动随 θ 增大而增加水平线# 纯壁摩擦决定的量对 θ 不敏感物理上即打滑段压力增长由摩擦定# - 咬入区梯度 G_nip(θ)受粉体可压缩/内摩擦驱动随接近辊隙(θ→0)# 而单调增大越窄压力爬坡越陡# 逐字公式请以 Johanson(1965)[ASME J Appl Mech 32(4):842–848] 原文/文献为准。# 本实现用于复现两梯度在 α 处匹配 → 求根这一结构不冒充原文。# ---------------------------------------------------------------defG_slip(theta,phi,D,S):滑移区归一化应力梯度壁摩擦主导取物理合理的量级。# 单位检查theta 为弧度此处记为“结构示意”绝不宣称逐字等于 Johanson 原式return4.0*phi# 量级上与壁摩擦角相关工程近似见声明defG_nip(theta,delta,D,S):咬入区归一化应力梯度越靠近辊隙(θ→0)压力爬坡越陡。# 结构随 D/S 增大辊隙相对更窄整体梯度增大随 θ 减小而增大。# 声明逐字式以 Johanson(1965) 原文为准。K_geoD/S# 几何压缩比无量纲returnnp.log1p(K_geo)*(delta/(deltatheta1e-9))# 2. 匹配方程在奇数角 θ 处使两梯度相等 → 根即 αdefmatch(theta):returnG_slip(theta,phi,D_roll,S_gap)-G_nip(theta,delta,D_roll,S_gap)# 数值求根brentq 需要两端异号否则换区间lo,hi1e-3,np.radians(40.0)try:alphabrentq(match,lo,hi)exceptValueError:alphanp.nanprint(注意当前 δ/φ 下 [lo,hi] 内无根——需检查摩擦角取值或搜索区间)alpha_degnp.degrees(alpha)print(咬入角 α %.2f rad %.2f deg%(alpha,alpha_deg))# 验证α 处两条梯度相等print(|G_slip(α) − G_nip(α)| ,abs(match(alpha)))# ---------------------------------------------------------------# 3. 画出两区归一化应力梯度随角度变化标注交点 α呼应 figs/fig02# ---------------------------------------------------------------thetanp.linspace(1e-3,np.radians(30),400)gsnp.array([G_slip(t,phi,D_roll,S_gap)fortintheta])gnnp.array([G_nip(t,delta,D_roll,S_gap)fortintheta])plt.figure(figsize(8,4.8))plt.plot(np.degrees(theta),gs,lw2,color#e67e22,label滑移区 σ/σ (壁摩擦 φ))plt.plot(np.degrees(theta),gn,lw2,color#3a7bd5,label咬入区 σ/σ (内摩擦 δ))ifnp.isfinite(alpha):plt.axvline(alpha_deg,ls--,color#555)plt.text(alpha_deg,max(gs.max(),gn.max())*0.7,α%.2f°%alpha_deg,rotation90,fontsize10)plt.xlabel(自辊隙中平面的角位置 θ [deg]);plt.ylabel(归一化应力梯度 1/σ dσ/dθ)plt.title(两区归一化应力梯度匹配求解咬入角 α结构示意逐字式以 Johanson 1965 为准)plt.legend(frameonFalse);plt.tight_layout()plt.savefig(figs/fig02_johanson_demo.png,dpi150);plt.show()逐行要点单位纪律铁律 3角度一律np.radians转弧度再进三角函数与求根避免度/弧度混用。φ、δ是弧度制theta也是弧度。这是本篇最容易翻车的点。结构 vs 逐字式G_slip/G_nip用壁摩擦主导内摩擦几何压缩比主导的物理合理结构代码与图片上都标注逐字式以 Johanson(1965) 原文/文献为准。真做工程仿真时把这两个函数体替换成原文推导的精确表达式或从文献照抄调用骨架与求根逻辑不用动——这就是结构 数值求根的可复用价值。**brentq/brentq 双侧异号**求根前确认match(lo)与match(hi)异号否则抛ValueError代码用 try/except 拦截并给出诊断信息。图输出相对路径figs/fig02_johanson_demo.png主文引用的示意图是figs/fig02_johanson.png两区示意本脚本额外画梯度匹配求 α的数学演示图二者互补但一个是概念示意、一个是数值复现。三、常见报错与排查角度单位混用度当弧度 / 弧度当度α 求出来不是负就是 NAN。现象α离奇几百上千度或求根失败梯度曲线形态异常。根因np.radians只转换了你记得的那些而φ/δ/theta之间单位不一致match无法在合理区间异号。解法在脚本开头统一把全部角度转弧度并用一处assert np.degrees(phi)phi_deg自检凡是进入三角函数/求根的量都写弧度注释。梯度匹配无根 / 初始猜测差ValueError或算出毛刺。现象brentq抛f(a) and f(b) must have different signs。根因当前的 δ/φ/D/S 组合下两条梯度在[lo,hi]内不相交例如壁摩擦太小滑移区梯度始终低于咬入区梯度。解法这不是 bug 而是参数空间问题——扩大搜索区间、或调大 φ1.4 节提到提高壁摩擦可增大咬入角正是这种参数敏感性也可以在求根前打印match(lo), match(hi)确认异号。把 SGAP 当 SFR把SF_gap直接当成最终 ribbon 密度。现象仿真结果的 ribbon 密度系统性偏高且与实测对不上。根因Johanson 只给SF_gap辊隙处、受约束未计入出间隙后的弹性回弹ERSF_ribbon SF_gap − ER一定更低底账 B1。解法不要在本篇就把SF_gap断言为 ribbon 实现密度为 ribbon 密度留一个后续 ER 修正的减项篇 7 落地并把Johanson 高估 ribbon 密度Dec 2003、Bi 2014 转引写进模型文档。忽略摩擦角测量的过度拟合。现象调 φ/δ 能把某一组辊压数据拟合得异常漂亮但换一组设备/处方就崩。根因nip 角与摩擦角难实验室测量底账 B1δ/φ 往往靠反求自由度大可拟合但泛化差。解法把 δ/φ 当作需跨工况标定并留置信区间的力学参数不要只针对单点数据过拟合用物理趋势咬入角 17°、随 φ 增而增做合理性约束。辊速/螺旋转速硬塞进 Johanson。现象想研究转速对密度影响于是给药 JT模型加了rpm参数结果与其他文献对不上。根因Johanson 是准静态模型忽略辊速与螺旋转速底账 B1辊速对咬入角的实测影响26°→9°是不在此模型内的真实响应。解法别把 rpm 硬塞进 Johanson要含辊速效应改用下一篇 Reynolds/螺杆控制模型Johanson 管几何应力、Reynolds 管操作响应。四、动手练习梯度匹配求 α 的基线用底账 B4 参数D250 mm、S2.5 mm选一组 δ、φ如上述示例跑脚本求出 α。判定脚本|G_slip(α)−G_nip(α)|1e-6且 α 为正值量级落在数度到十几度把 φ 增大 5°α 单调增大呼应提高壁摩擦可增大咬入角。辊隙敏感性固定 D、δ、φ把 S 从 3.0 mm 扫到 1.5 mm底账 B4 范围观察 α 是否随 S 减小而略增几何压缩比变陡。判定画出 α–S 曲线α 平滑变化、无奇点并能解释间隙越窄、咬入越提前。诚实边界练习写一个断言assert SF_ribbon SF_gap把SF_gap想成一个正标量、SF_ribbon由它减去 ER 比例得到。判定该断言在给定 ER∈(0,1) 下恒成立并能在注释里写出Johanson 高估 ribbon 密度Dec 2003、Bi 2014 转引这一本文核心局限的出处。五、小结与下一篇预告本篇把干法制粒辊压的力学与几何底层钉住了理解了对辊粉床 滑移区脱气、辊面快 咬入区无滑移、被压实大部分压实发生在这里明白了咬入角 α 是两区交界角、由滑移/咬入两条归一化应力梯度在 α 处匹配唯一确定吃透了 Johanson 的输入压缩性、δ、φ、D/W/S、辊压与准静态假设忽略辊速/螺旋转速并诚实记下了它**止步于SF_gap、因未计弹性回弹而高估 ribbon 密度SF_ribbon SF_gap**这一关键边界。这段机理是自定义辊压单元的第一块砖。下一篇篇 07《辊压②Midoux 数与 Reynolds 模型》——用 Midoux 无量纲数做 SCF→Pmax 的解析放大Pmax 2·SCF/D·√(2K·D/(π·S))并用 Reynolds 模型补上 Johanson 欠缺的螺杆/辊速响应最后接上弹性回弹 ER 修正把SF_gap换成最终SF_ribbon真正补齐 Johanson 高估的缺口。本篇认知问题回显FAQQ1为什么对辊里要拆成滑移区和咬入区为什么大部分压实发生在咬入区A1按辊面线速度与粉体速度的关系Johanson(1965) 把粉床分两段滑移区slip最外侧辊面快于粉体、粉体打滑、主要作用是脱气与轻微压缩和咬入区nip辊隙最窄处粉体被咬住无相对滑移、被刚性压实。大部分压实密度提升发生在咬入区——只有那里粉体才真正两点受约束被压密。这是一个用力学真实还原看似全程在压的收敛答案也是咬入区压力-密度关系成为模型核心的原因。Q2咬入角 α 怎么从数学上求为什么是两条梯度相交A2α 定义为滑移区与咬入区两条归一化应力梯度1σdσdθ\frac{1}{\sigma}\frac{d\sigma}{d\theta}σ1​dθdσ​在 α 处数值相等的角度梯度连续性。道理是跨两区分界处压力随角度的斜率不能突变否则粉体介质被撕开滑移支线由壁摩擦 φ 驱动、咬入支线由内摩擦 δ 驱动二者交点的角度就是咬入角 α。求解方式即数值求根求G_slip(θ)−G_nip(θ)0。Q3为什么模型可以忽略辊速和螺旋转速A3因为 Johanson(1965) 采用准静态假设——把压力-密度平衡当作与时间/转速无关的纯几何摩擦问题。所以模型输入里没有辊速、螺旋转速。代价是它管不了辊速对压实的影响这类动力现象下一篇 Reynolds 模型正是为螺杆/辊速响应而存在两者互补Johanson 适合间隙控制、Reynolds 适合螺杆控制。Q4为什么直接套用 Johanson 会高估 ribbon 密度SF_gap 与 SF_ribbon 差在哪A4模型输出的固相分数是辊隙gap处粉体仍受约束计算的SF_gap是一个乐观上限ribbon 一旦流出间隙、辊面约束消失弹性应力释放引起弹性回弹ER密度下降最终SF_ribbon SF_gap。Johanson 未计入这一步所以直接套用高估 ribbon 密度Dec 2003、Bi 2014 转引。SF_gap到SF_ribbon的换算由下一篇 Midoux/Reynolds ER 修正完成。Q5咬入角 26°→9°、临界 17° 在什么条件下得到能直接复用吗A5这些不出自 Johanson 解析结论而是文献/实测与推导值底账 B1制药临界咬入角一般 17°且提高壁摩擦可增大咬入角设计杠杆MCC 恒定压力/间隙下把辊速从 3 提到 7 rpm实测咬入角从 26° 降到 9°Al-Asady 转引。二者都不能直接当普适常数塞进自己的仿真——应作为数量级与行为趋势的参考并结合本设备物料/几何重新求取。