企业官网建设流程全解析

PharmaPy 干法制粒二次开发教程（10）：粉碎模块②——离散破碎 PBM 的数值实现：Hounslow 几何网格与 Kumar–Ramkrishna 固定枢轴法

版本声明块

• 工具/软件：本项目二次开发的离散破碎 PBM 求解函数（本系列二次开发示例，非官方 API），基于 Hounslow 几何网格 + Kumar–Ramkrishna（1996）固定枢轴法
• 语言/环境：Python 3.9 + numpy/scipy（scipy.integrate.solve_ivp）
• 本文目标：把"一粒粒看概率"升格为"整批 PSD 随停留时间演化"，最关键的是每一格算完都验算质量守恒。

一句话结论：在 Hounslow 对数几何网格（固定粒径比）上，破碎 PBE 离散为dN_i/dt = Σ_{j≤i} S_j·B_ij·N_j − S_i·N_i；用 Kumar–Ramkrishna（1996）固定枢轴法构造守恒的B_ij，solve_ivp(method='BDF')积分，每步用np.trapezoid(distrib, x)（或质量份额求和）复验总质量漂移< 相对容差——这是干法制粒破碎仿真可信度的底线。

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

1. 破碎 PBM 的"出生、死亡"两项，各自在离散网格上长什么样？
2. 一个母粒碎成多个子粒，子粒落到不恰好在网格节点上的粒径怎么办？——固定枢轴法解决了什么？
3. 为什么网格要"对数几何"？粒径比恒定（如2^(1/3)）对守恒与B_ij构造为什么重要？
4. solve_ivp积分很顺，但我怎么知道它没"偷偷丢掉质量"？
5. 篇 9 的S(x)概率核，在这里怎么变成一个逐档的速率系数S_i？

一、机制解析

破碎 PBE 的出生-死亡结构（Bilgili 2024，底账 C2）

单个颗粒的"碎与不碎"（篇 9）升级为整批 PSD，核心是一个连续的群体平衡方程：对粒径i，其数量密度随时间的变化 = 出生 − 死亡。

• 死亡项−S_i·N_i：本档颗粒以速率S_i破碎，出去了。
• 出生项Σ_{j≤i} S_j·B_ij·N_j：所有比i大的档j破碎时，有一部分碎片落进i档。B_ij是"j碎的碎片中落进i的比例（破碎函数）"，S_j是j档的选择函数（比破碎速率）。

写成标准形式（Austin/Bilgili 批式，底账 C2）：

dN_i/dt = Σ_{j≤i} S_j·B_ij·N_j − S_i·N_i

这个方程正确的前提是几个假设：磨内近似全混、无团聚/生长、破碎一阶线性（即S_i不随批次变化）。其中S_i受操作条件（转速）强烈影响，而B_ij相对不敏感——这决定了篇 11 参数化的分工：S_i管"多快碎"，B_ij管"碎成什么样"。

连续全混磨（含停留时间与内分级，底账 C2，篇 11/12 铺垫）

真实工业磨是连续流，Bilgili（2024）给连续全混式：d(H·M_i)/dt = f_i − p_i − (S_i/d_i)·p_i + Σ_{j<i} b_ij·(S_j/d_j)·p_j，p_i = d_i·H·M_i。其中d_i是排料速率，若有粒度分级效应，可用归一化排料速率d_i*=τ·d_i充当平均分级函数——这就是"分级进入粉碎模型"的接口，为篇 11 的筛网和篇 12 的 Tromp 埋线。本篇先从更干净、更好验算的批式落地，连续全混式留到篇 20 完整项目。

为什么必须"固定枢轴"：碎片不落在节点上

数值上我们不能连续存粒径，只能用一组离散网格。难点在出生项：母粒j碎出的子粒体积是连续分布的，不可能恰好落在某个网格节点上。若强行把落在两节点间的子粒分给两端，需要既不偷质量、也不伪造数量密度——Kumar–Ramkrishna（1996）Fixed Pivot正是干这个：对每个离散网格（固定枢轴），选择一组约束（通常保留两个矩，如总质量/总体积与总数），把任一落点分配给它两个相邻枢轴，使这两个矩在离散化前后保持。对破碎这种只关心"子粒体积和 = 母粒体积"（零阶+一阶矩守恒即可）的场合，固定枢轴法能构造出精确的质量守恒的B_ij。

序列级术语提示：Hounslow et al.（1988/1990）提出"四倍"粒径几何网格与差分格式；Kumar–Ramkrishna（1996）的 Fixed Pivot 是更一般地把连续 PBE 换成离散 ODE 并保矩的方案——它就是本系列落地破碎 PBM 的首选（底账 G，可靠性优于直接四倍网格的启发式差分）。

网格共轭：对数几何 + bin→守恒的关键

铁律 5 规定"对数等距（粒径比恒定）"。原因有三：①制药 PSD 常跨 1.5 个数量级，对数网格稀疏处不浪费 bin；②固定粒径比让 Bini的平均体积与i+1满足常数比值，出生项的分配向量可预计算；③质量守恒核验用Σ distrib·Δx时，几何网格的Δx随档位变化，正好把"每档占多少体积"体现进来。本系列取x_gr=np.geomspace(10,2000,num=40)（工程选型 30~50 bin，底账 G）。

质量守恒就是验收试金石

数值格式是不是守恒、积分有没有发散，唯一硬判据是总质量不回飘：批式破碎下"母粒子破碎不产生也不消灭质量"，所以积分结束后Σ M_i（质量份额）应与初值一致。这也是"固定枢轴法 + 网格共轭 + 每步验算"三者缺一不可的原因——任一处出错，质量守恒断言立刻暴露。这也是本系列相对"跑出来好看"的普通教程最硬的纪律：一切数值结论都先过守恒这一关。

二、完整代码与逐行剖析

下面用真实可运行的示例实现：Hounslow 几何网格 → 一套S_i（承接篇 9 的 Vogel–Peukert）→ 固定枢轴法构造B_ij→solve_ivp(method='BDF')积分 → 质量守恒核验。S_i/B_ij的具体参数化是篇 11 的正文，这里先用一个幂律S_i与"等质量双子碎"B_ij占位，重点演示数值框架与守恒纪律。

importnumpyasnpfromscipy.integrateimportsolve_ivp# ---- 1) Hounslow 对数几何网格（固定粒径比，铁律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:])# bin 中心（几何中点）dx=x_gr[1:]-x_gr[:-1]# bin 宽nb=numbins-1# bin 数# ---- 2) 选择函数 S_i：承接篇9 概率核，化作逐档速率 (示意参数) ----# 简单幂律：S_i ∝ (x_c / x_ref)^qdefselection(x_c,k0=1.0,q=0.8):returnk0*(x_c/500.0)**q S=selection(x_c)# ---- 3) 破碎函数 B_ij：固定枢轴法构造"双子均分"，保体积/质量守恒 ----# B_ij = b_{i-1,j} 形式：母粒 j 恰好碎成两个体积均分的子粒，# 落在两相邻枢轴 (i, i+1) 间按线性分配，使其体积和=母粒。# 这里做示意实现：子粒等分到最接近的两个 bin，保质量份额。defbreakage_fixed_pivot(x_c):n=x_c.size B=np.zeros((n,n))vol=x_c**3.0# 以体积代表质量份额(密度均匀)forjinrange(n):v_child=0.5*vol[j]# 子粒体积(等分2粒)# 找到 v_child 落点邻近最大枢轴i=int(np.searchsorted(vol,v_child))# 第一个 vol>=v_childi=max(1,min(n-1,i))# 分配比例：使两枢轴体积加权和保持——w1=(vol[i]-v_child)/(vol[i]-vol[i-1])# 分配给 i-1 的份额B[i-1,j]+=w1 B[i,j]+=(1.0-w1)returnB B=breakage_fixed_pivot(x_c)# ---- 4) PBM ODE：dN_i/dt = Σ_j S_j B_ij N_j - S_i N_i ----defpbm_rhs(t,N):birth=np.dot(B,S*N)# Σ_j B_ij * (S_j N_j)death=S*Nreturnbirth-death N0=np.zeros(nb);N0[8]=1.0# 初始：单档(50µm附近)一批母粒M0=np.sum(N0*vol)# 初始质量份额# ---- 5) 积分 + 质量守恒核验 ----t_end=50.0sol=solve_ivp(pbm_rhs,(0.0,t_end),N0,method='BDF',dense_output=True,rtol=1e-8,atol=1e-10)Nf=sol.y[:,-1]Mf=np.sum(Nf*vol)drift=abs(Mf-M0)/M0# 相对质量漂移print(f"初带档: bin8(x_c={x_c[8]:.1f}µm) 末带质量份额相对漂移 ={drift:.2e}")assertdrift<1e-6,"质量守恒失败，请检查 B_ij 构造"# ---- 行为核验：粗粒应当变少、细粒变多（纯破碎）----coarse0=np.sum(N0[x_c>800.0]*vol[x_c>800.0])coarse1=np.sum(Nf[x_c>800.0]*vol[x_c>800.0])print(f">800µm 质量份额: 初{coarse0:.4f}-> 末{coarse1:.4f}")assertcoarse1<coarse0,"纯破碎不应出现粗粒反而增多"print("OK：质量守恒 + 纯破碎右移行为均通过")

逐行要点：

• ratio=(x_hi/x_lo)**(1/(n-1))保证固定粒径比（几何网格），满足铁律 5；x_c用几何中点承载 bin 体积。
• selection用幂律S_i ∝ (x_c/500)^q占位（示意参数，篇 11 才做参数化）；它承接篇 9 的vogel_peukert_breakage_probability——本篇用幂律是"先搭守恒框架"，篇 11 将用物理核替换。
• breakage_fixed_pivot：把v_child=vol[j]/2（双子均分）落在两枢轴间，按体积差w1=(vol[i]−v_child)/(vol[i]−vol[i−1])分配，恰好保体积守恒——B满足Σ_i B_ij = 1（每列和为 1），这是守恒的代数保证。
• pbm_rhs里birth=np.dot(B, S*N)用矩阵乘法一次性算完全部出生项，避免逐 bin 双重循环——性能上这是工程关键（篇 19 会谈向量化批量）。
• solve_ivp(..., method='BDF')：=刚性系统（破碎速率差别可达几个数量级）用隐式 BDF 更稳，rtol/atol收紧保证守恒精度。
• 双断言：质量漂移<1e-6（守恒底线）与粗粒质量减少（破碎右移行为）——把"对错"变成可自动判定的标准。

三、常见报错与排查

质量守恒断言失败，漂移很大。
现象：drift达到1e-2及以上。根因绝大多数在B_ij：breakage_fixed_pivot的分配索引i越界/落点分配不合理，导致Σ_i B_ij ≠ 1。解法：先单独打印B.sum(axis=0)（每列和），必须全为 1；再检查i的 clamp 逻辑（i=max(1,min(n-1,i))）是否把边界档处理对。

用method='RK45'积分发散/NaN。
现象：solve_ivp报步长过小或出现 NaN。根因：破碎系统是刚性，显式 RK45 处理不了跨数量级的速率差，或无足够精度。解法：改用method='BDF'（或LSODA），并把rtol/atol收紧到1e-8/1e-10。

初始条件只在"节点"而不在"边界"，B维度错。
现象：B赋值时报索引错误或N0长度与B不匹配。根因：nb = numbins-1是 bin 数，x_c/dx/N/S/B都应是nb长；若把x_gr（节点长numbins）混入N，维度就炸了。解法：统一所有状态数组长度为nb，把x_gr只作网格坐标，别直接当状态长度用。

固定枢轴法给双子均分太重，子粒晚于 1 步又碎。
现象：B把双子只落在两相邻档，若这两档本身又很大，会反复破碎。根因：这是"双子均分"占位模型的局限（不是守恒问题），真实离散应让子粒分布更广。解法：明确本篇是数值框架演示，真正子粒分布b_ij的物理参数化在篇 11——本篇只保证"结构正确 + 守恒"，不追求子粒形态逼真。

初值只放 bin8，行为断言还通过了吗？
现象：coarse1 < coarse0可能不触发但意义不大。根因：x_c>800的初值只有 bin8 一档，初就极少。解法：练习里建议用一个"宽的平整初值"（如N0=np.ones(nb).copy(); N0[x_c<100]=0之类）再验右移，行为断言才有统计意义。

四、动手练习

1. 守恒自验（必做）：对多个不同N0（单档、平带、偏细、偏粗），各跑一遍求drift。判定：所有情况下drift<1e-6，证明你的B_ij与solve_ivp配置在不同初值下都守恒。
2. 网格收敛：把numbins从40试到80/160，观察同一t_end下 PSD 的d50（用篇 14 的算法）是否收敛、drift是否随网格加密而改善。判定：d50变化<5%且drift持续<1e-6，证明网格密度已够。
3. 把篇 9 概率核接进selection：把幂律S_i替换为vogel_peukert_breakage_probability(x_c*1e-6, W_m_kin)的映射，重跑守恒与行为断言。判定：守恒依旧<1e-6，且随W_m_kin增大 PSD 右移更快——打通"篇 9 概率核 → 篇 10 群体 PBM"。

五、小结与下一篇预告

本篇完成了破碎 PBM 的数值落地：Hounslow 几何网格 + Kumar–Ramkrishna 固定枢轴法构造守恒B_ij，solve_ivp(method='BDF')积分批式dN_i/dt=Σ S_j B_ij N_j − S_i N_i，并用质量漂移<1e-6与粗粒减少两个断言把可信度钉死。你已掌握"离散化守恒格式 + 刚性积分 + 守恒验收"这个能被 PharmaPy 与任何项目复用的内核。

下一篇：篇 11《粉碎③：选择函数 S_i 与破碎函数 B_ij 参数化》——把占位的幂律S_i和双子B_ij换成有物理的模型：S_i的大小依赖、B_ij的碎裂模式与 GGS 型子粒分布、以及筛网内的上限粒径分类d_i。

本篇认知问题回显（FAQ）

Q1：破碎 PBM 的出生、死亡项离散后长什么样？
A1：对档位i，dN_i/dt = Σ_{j≤i} S_j·B_ij·N_j − S_i·N_i——死亡是−S_i·N_i（本档以S_i破碎离开），出生是Σ_{j≤i} S_j·B_ij·N_j（所有更大档j破碎时按B_ij落入i）。B_ij满足Σ_i B_ij=1（守恒代数保证）。

Q2：母粒碎片不在网格节点上怎么处理？
A2：用 Kumar–Ramkrishna（1996）固定枢轴法：把碎片子粒按其体积/质量落点分配给两个相邻枢轴，按体积差线性加权w1=(vol[i]−v_child)/(vol[i]−vol[i−1])，保证离散前后保留两个矩（破碎场景即总质量/总数）。这是"不偷质量、不伪造数密度"的标准保矩离散方案。

Q3：为什么要对数几何网格？
A3：铁律 5。①PSD 跨 ~1.5 个数量级，对数网格在粗端不浪费 bin；②固定粒径比（如2^(1/3)）使相邻 bin 平均体积成常数比，出生分配可预计算；③守恒核验用Σ distrib·Δx时，几何网格的Δx随档变化恰好体现每档体积权重。本系列取geomspace(10,2000,40)（工程选型 30~50）。

Q4：怎么知道solve_ivp没偷偷丢质量？
A4：只信一道质量守恒断言：批式破碎不产生、不消灭质量，故积分结束后Σ N_f·vol / Σ N_0·vol应等于 1 到很高精度（本文<1e-6）。配合B.sum(axis=0)==1的代数检查，以及method='BDF'（刚性）+ 收紧rtol/atol，就能把积分数值误差压到守恒线以下。

Q5：篇 9 的S(x)概率核怎么变成速率系数S_i？
分结束后Σ N_f·vol / Σ N_0·vol应等于 1 到很高精度（本文<1e-6）。配合B.sum(axis=0)==1的代数检查，以及method='BDF'（刚性）+ 收紧rtol/atol，就能把积分数值误差压到守恒线以下。

Q5：篇 9 的S(x)概率核怎么变成速率系数S_i？
A5：本篇第一步用幂律S_i=k0·(x_c/500)^q（示意参数）搭守恒框架；真正物理的做法是把vogel_peukert_breakage_probability(x_c*1e-6, W_m_kin)映射成逐档S_i（练习 3 已打通）。篇 9 给"碎与不碎的概率"，篇 10 把它乘以单位时间碰撞强度变成"单位时间破碎速率"。