企业官网建设流程全解析

PharmaPy 干法制粒二次开发教程（11）：粉碎模块③——选择函数 S_i 与破碎函数 B_ij 参数化与筛网内分级

版本声明块

• 工具/软件：本项目二次开发的破碎参数化（本系列二次开发示例，非官方 API）
• 语言/环境：Python 3.9 + numpy/scipy（scipy.special与np.interp筛网截断）
• 本文目标：把篇 10 的骨架换成有物理血肉的S_i/B_ij，并引入"筛网孔径截断可破碎档"。

一句话结论：选择函数S_i = S_ref·(x_i/x_ref)^q·ω_scale(转子速)随粒径幂律上升、随叶轮转速放大（操作依赖）；破碎函数B_ij用"局部崩解 + 多重碎裂"的分支子粒分布，且其累积服从 GGS 型Y=100(x/k)^α；并令筛网孔径把x_i>筛孔的档直接送破碎、≤筛孔档保留——三者合成后在篇 10 的solve_ivp里跑出"粗粒右移、被筛上限封顶"的合理 PSD。

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

1. S_i（多快碎）与B_ij（碎成什么样）为什么该分开参数化？各自的操作依赖是什么？
2. "选择函数随粒径幂律上升"直觉上对不对？叶轮转速是怎么进来的？
3. 破碎函数的分支子粒分布：局部崩解与多重碎裂（泛碎裂）区别在哪？GGS 子粒分布怎么写？
4. 筛网在破碎 PBM 里到底怎么"截断"？是挡住"大档"还是把大档强制变碎？
5. 参数化完之后，我怎么判断 PSD 输出是"合理的右移 + 上限封顶"而不是乱七八糟？

一、机制解析

S_i与B_ij的分工：速率 vs 形态

前面两条铁律性事实（底账 C2）：S_i受操作条件（转速）强烈影响，B_ij相对不敏感。这是把两者分开建模的根本动机——也是篇 17 参数标定时"把S_i各转速分别标、B_ij全转速共用一个"的依据。分工一句话：S_i告诉你"哪些档、多快被碎"，B_ij告诉你"碎了之后碎成什么尺度分布"。

选择函数S_i：粒径幂律 + 操作放大

物理直觉：大的颗粒更脆、更容易被撞碎→S_i随x_i增大而增大，普遍用幂律S_i = S_ref·(x_i/x_ref)^q（q>0，典型 0.6~1.2）。这里的x_ref/S_ref是"参考粒径处的比破碎速率"，是标定量。而操作依赖（尤其叶轮尖端速度，前述锥磨决定细粉的最大因素，底账 C1）用一个ω_scale乘子进来：ω_scale = (u_tip/u_ref)^p，叶轮转越快，单位时间穿透、撞击次数越多，破碎速率同比放大。Dickson 把S_ref标成转速的函数即可。这样S_i就同时含"粒径依赖"（幂律）与"操作依赖"（转速乘子）两条轴，与真实粉碎一致——同一配方，转速越高碎得越多、细粉越多。

破碎函数B_ij：分支子粒分布（局部崩解 vs 多重碎裂）

真实破碎不是"只变成两颗"，而是多种碎裂模式并存。Reynolds（2010，底账 C5）对制药锥磨提出通用碎片分布函数，可看成两种极端模式的叠加：

• 局部崩解（surface/edge chipping）：母粒只崩掉几个小碎片，剩余"核心"仍接近母粒大小 → 产生细粉 + 保留粗核。
• 多重碎裂 / 泛碎裂（massive breakage）：母粒整体碎成较均匀的若干子粒 → 产生拉宽的中间峰。

工程上常用一个分支系数把两者按比例混合：母粒以β份额整体碎裂成人n个体积均分子粒（多重碎裂），以(1−β)份额只崩出少量小碎片（局部崩解）。子粒的累积分布尤其常用GGS（Gaudin-Schuhmann）表述：Y=100(x/k)^α（底账 C4），k是粒度模数（≈最大子粒尺寸/母粒尺寸的函数），α是分布陡度。GGS 的好处是一条对数轴的直线、斜率即α，标定非常直观：测出子粒 PSD，取对数线性回归就得k、α。

序列提示：GGS 的姊妹 RRB 形F(d)=1−exp(−(d/λ)^k)在篇 14（PSD 报告层）与 lognormal 一起作为拟合家族的候选；本篇先用 GGS 描述子粒分布，因为它直接定义了"母粒碎成多大的子粒"这个破碎函数本体。

筛网：截断可破碎档，而非"滤掉大粒"

关键澄清：在粉碎机本体里，筛网是破碎几何约束——x_i大于筛孔的那几档，不属于"能直接通过的稳定粒"，它们会在转子撞击 + 过筛挤压下继续破碎（或至少被判定为"必碎"）；而x_i不超过筛孔的档是"可稳定存活"的。所以模型里我们用筛网孔径mesh把S_i对x_i>mesh的档截断放大（工程近似的做法是把这些档设为"必碎"、即S_i置大），并把B_ij的子粒上限封在mesh附近。这样最终 PSD 的右上端被筛孔压住——这正是"锥磨给出一致上限粒径、收窄 PSD"（底账 C1）的数值体现。

边界（真实性铁律）：S_i幂律指数q、B_ij分支比β、GGS 的k/α都由物料与磨机标定（篇 17 反演），本篇用示意参数演示结构和趋势，不宣称通用数值。

二、完整代码与逐行剖析

下面是真实可运行的参数化示例：定义选择函数（粒径幂律 + 转速放大 + 筛网截断）与破碎函数（分支 + GGS 子粒分布），再灌进篇 10 的solve_ivp求 PSD，最后做"合理右移 + 上限封顶"行为断言。

importnumpyasnpfromscipy.integrateimportsolve_ivp# ---- 网格（同篇10，铁律5）----x_lo,x_hi,numbins=10.0,2000.0,40ratio=(x_hi/x_lo)**(1.0/(numbins-1))x_gr=x_lo*ratio**np.arange(numbins)x_c=np.sqrt(x_gr[:-1]*x_gr[1:])dx=x_gr[1:]-x_gr[:-1]nb=numbins-1vol=x_c**3.0# ---- 1) 选择函数：粒径幂律 + 叶轮转速放大 + 筛网截断 ----defselection(x_c,x_ref=500.0,S_ref=0.5,q=1.0,rpm_tip=2000.0,rpm_ref=2000.0,p=1.2,mesh=1000.0):S=S_ref*(x_c/x_ref)**q# 粒径幂律S=S*(rpm_tip/rpm_ref)**p# 操作依赖(rpm)S=np.where(x_c>mesh,10.0,S)# 筛网截断：>mesh 必碎(示意)returnS# ---- 2) 破碎函数：分支(多重碎裂+局部崩解) + GGS 子粒分布 ----defbreakage_ggs(x_c,beta=0.6,alpha=0.7,mesh=1000.0):n=x_c.size B=np.zeros((n,n))forjinrange(n):xm=x_c[j]# GGS 累计子粒分布：Y(<x)=100(x/k)^α；k 与母粒尺寸成正比k=0.5*xm# 多重碎裂：子粒均匀撒在 [小端:k]，用 GGS 累积求每档份额frac=np.zeros(n)foriinrange(j+1):# 只到母粒档(i<=j)Y_i=(x_c[i]/k)**alphaifk>0else0.0upper=(x_gr[i+1]/k)**alpha frac[i]=min(1.0,upper)-min(1.0,Y_i)# 档内 GGS 累计增量frac=frac/frac.sum()# 归一保证 ΣB_ij=1# 有份额的档被"筛网截断"压到 <= meshallowed=x_c<=mesh frac=np.where(allowed,frac,0.0)frac=frac/frac.sum()B[:j+1,j]=frac[:j+1]returnB S=selection(x_c)B=breakage_ggs(x_c)# ---- 3) PBM 积分（复用篇10）----defpbm_rhs(t,N):returnnp.dot(B,S*N)-S*N N0=np.zeros(nb)N0[x_c<=600.0]=1.0# 较宽的进料带：≤600µmN0[x_c>600.0]*=0.0# (示意) 初值以中细档为主M0=np.sum(N0*vol)sol=solve_ivp(pbm_rhs,(0.0,20.0),N0,method='BDF',rtol=1e-8,atol=1e-10)Nf=sol.y[:,-1]drift=abs(np.sum(Nf*vol)-M0)/M0print(f"守恒漂移 ={drift:.2e}")assertdrift<1e-6# ---- 4) 行为断言：合理右移 + 上限封顶 ----above=x_c>1000.0below=x_c<=100.0print(f">1000µm 质量份额初/末:{np.sum(N0[above]*vol[above])/M0:.3f}-> "f"{np.sum(Nf[above]*vol[above])/np.sum(Nf*vol):.3f}")print(f"<100µm 质量份额末:{np.sum(Nf[below]*vol[below])/np.sum(Nf*vol):.3f}")assertnp.sum(Nf[above]*vol[above])/np.sum(Nf*vol)<1e-9,"筛网应把>1000µm压没"print("OK：守恒 + 筛上限封顶 + 右移，破碎参数化合理")

逐行要点：

• selection三段的顺序即两条轴：先粒径幂律，再转速乘子，最后筛网截断把>mesh档S_i置10.0（"必碎"示意）。mesh就是筛网孔径（如 Comil 的 G043 或 1.0 mm 锉纹筛，底账 C1 给G040 锉纹筛≈1.0mm等），工程上把它当 PSD 右上限。
• breakage_ggs用 GGS 累积分布Y=(x/k)^α求每档份额，k=0.5·x_m让最大子粒约为母粒一半（示意），frac归一保证Σ_i B_ij=1（守恒代数，篇 10 铁律）。
• 筛网在破碎函数里的第二个作用：allowed = x_c<=mesh把B_ij中"子粒仍>筛孔"的档清零再归一，保证碎产物也被封在筛下。这是"筛网同时管上限粒径与破碎程度"的数值落地。
• 行为断言分两半：守恒<1e-6；上限封顶>1000µm份额压到<1e-9。至于"右移"看输出的<100µm/粗档对比（受初值影响，练习 3 给平带的严格右移验收）。
• 这是"参数化"演示，S_ref/q/p/β/α/k全为示意参数，正文与注释反复标注"以标定/文献为准"（真实性铁律）。

三、常见报错与排查

筛网截断后守恒崩了。
现象：加S=np.where(x_c>mesh,10,B)后drift变大。根因：把S_i置10.0是"快速破碎"而非"直接消失"——S_i只影响速率，不消灭质量；但若你的B_ij没同步把子粒压到筛下，母粒可能"碎出的也是>筛孔"，最终仍有粗粒存留导致守恒看着怪。解法：确认breakage_ggs里allowed的frac清零再归一（第 2 步），让破碎产物也被封顶；守恒应该仍<1e-6。

S_i幂律指数q给太大，小档数字退化。
现象：q=2.0时x_c<50的S_i趋近 0，几乎不碎，把 PSD"锁死"在细端。根因：幂律对细粒衰减过快，比破碎速率被低估。解法：把q控制在典型范围（0.6~1.2），或用带下限的S_i=S_ref·max((x/x_ref)^q, S_floor)，避免细档完全"冻住"。

子粒全落在母粒所在档导致"碎了个寂寞"。
现象：B_ij对角近 1，播完Nf≈N0，行为断言不触发。根因：breakage_ggs的k给太大/α给太小，GGS 累计在母粒邻档就饱和，子粒几乎留在母粒档。解法：检查k=0.5·x_m这类"最大子粒≈母粒一半"的设定，并让frac覆盖到更小的档（必要时k更小、α更大拉陡分布）。

mesh和x_c交错，封顶判断失败。
现象：>1000µm断言没压到<1e-9。根因：mesh用了筛网孔径的名义值，但x_c是几何中心——x_c[-1]（~1576µm）与mesh=1000之间还有档，若allowed只处理恰好<=mesh会留尾巴。解法：明确"筛孔对应 bin 上界"，用x_gr节点而非x_c中心判断，或在断言里留1e-3量级的合理残存（工程余量）。

转速p乘子把S_i顶到指数爆炸、ODE 刚性加剧。
现象：rpm_tip很大时S_i有些档 >>1，solve_ivp步长骤小变慢。根因：操作放大无上限，速率不在同一量级。解法：给S_i设上限（如clip到50.0或 mesh 区间的"必碎"值），再收紧atol；物理上无限快的破碎也无意义。

四、动手练习

1. 守恒回归：把mesh换800/1200 µm、β换0.3/0.9各跑一遍，确认drift<1e-6均成立。判定：守恒与参数值无关，只与B_ij归一的代数结构有关——证明你的破碎函数"对任何参数都不偷质量"。
2. 筛网敏感性：固定其余参数，扫描mesh=[800,1000,1200]并打印 PSD 的右端质量份额。判定：mesh越大右端越有余量、越小越压越紧，三者单调——还原"筛网孔径定上限粒径"，可作篇 16 工艺窗口的先验。
3. 转速敏感性（右移验收）：用平带初值（N0[x_c<=1200]=1，>1200为 0），rpm_tip=[1000,2000,4000]各跑一遍，打印<125µm细粉份额。判定：转速越高细粉份额越大（右移/细化增强），且全工况drift<1e-6——这正是"锥磨叶轮尖端速度决定细粉"的模型化体现，也为篇 16 细粉 CQA 埋线。

五、小结与下一篇预告

本篇把破碎骨架填上物理血肉：S_i=粒径幂律×转速放大×筛网截断，B_ij=分支（多重碎裂+局部崩解）×GGS 子粒分布×筛下限制，并在守恒与"上限封顶/右移"断言下验证。你已掌握分离、“速率 vs 形态”、且"筛网同时管上限粒径与破碎程度"这两条主线——这套S_i/B_ij/mesh三件套就是今后参数标定（篇 17）的待估对象。

下一篇：篇 12《分级模块：Tromp 分级效率与切割粒径》——破碎之后进入分级：用 Tromp 曲线量化"某个粒径进粗流/细流的比例"，算 d50 切割粒径，并看清开路/闭路（返料）循环物流在 flowsheet 里怎么搭。

本篇认知问题回显（FAQ）

Q1：为什么 S_i 和 B_ij 要分开参数化？各自依赖什么？
A1：铁律性事实是S_i受操作条件（转速）强烈影响、B_ij相对不敏感（底账 C2）。所以S_i把粒径（幂律q）与操作（叶轮转速乘子）都放进去，B_ij只需描述"碎了之后的子粒尺度分布"，跨转速可共用。这直接决定篇 17 标定时"转速分开标S_i、全转速共用B_ij"。

Q2：选择函数随粒径幂律上升直觉对吗？转速怎么进来？
A2：对的——大颗粒更脆更易被撞碎，S_i = S_ref·(x_i/x_ref)^q（q典型 0.6~1.2）。转速通过乘子(u_tip/u_ref)^p进来：叶轮尖端速度越高，单位时间撞击/过筛挤压越多，破碎速率同比放大。两轴合起来即"同配方转速越高碎得越多、细粉越多"（=锥磨叶轮尖端速度定细粉，底账 C1）。

Q3：局部崩解 vs 多重碎裂区别？GGS 子粒分布怎么写？
A3：局部崩解只崩掉少量小碎片、核仍近母粒（产细粉+保留粗核）；多重碎裂整体碎成多个较均匀子粒（拉宽中间峰）。用分支系数β混合两者。子粒累积分布用 GGS：Y=100(x/k)^α（k粒度模数≈最大子粒、α陡度），取对数成直线、斜率即α，标定直观。

Q4：筛网在破碎模型里怎么"截断"？
A4：不是"滤掉大粒"，而是几何约束——x_i>筛孔的档判定为"必碎"，把其S_i截断置大（本篇置10.0），同时把B_ij中"子粒仍>筛孔"的档清零，使碎产物也被封在筛下。两处共同作用，让最终 PSD 右上端被筛孔压住 = 锥磨"一致上限粒径、收窄 PSD"。

Q5：怎么判断 PSD 输出"合理右移 + 上限封顶"？
A5：两套断言：①守恒drift<1e-6恒成立（只与B_ij归一结构有关，与参数无关）；②上限封顶>mesh质量份额<1e-9（筛网压头）；③右移/细化用平带初值扫rpm_tip，转速越高<125µm细粉份额越大。三者同时满足才算参数化合理。