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

PharmaPy 干法制粒二次开发教程10粉碎模块②——离散破碎 PBM 的数值实现Hounslow 几何网格与 Kumar–Ramkrishna 固定枢轴法版本声明块工具/软件本项目二次开发的离散破碎 PBM 求解函数本系列二次开发示例非官方 API基于 Hounslow 几何网格 Kumar–Ramkrishna1996固定枢轴法语言/环境Python 3.9 numpy/scipyscipy.integrate.solve_ivp本文目标把一粒粒看概率升格为整批 PSD 随停留时间演化最关键的是每一格算完都验算质量守恒。一句话结论在 Hounslow 对数几何网格固定粒径比上破碎 PBE 离散为dN_i/dt Σ_{j≤i} S_j·B_ij·N_j − S_i·N_i用 Kumar–Ramkrishna1996固定枢轴法构造守恒的B_ijsolve_ivp(methodBDF)积分每步用np.trapezoid(distrib, x)或质量份额求和复验总质量漂移 相对容差——这是干法制粒破碎仿真可信度的底线。〇、本篇要解决的认知问题破碎 PBM 的出生、死亡两项各自在离散网格上长什么样一个母粒碎成多个子粒子粒落到不恰好在网格节点上的粒径怎么办——固定枢轴法解决了什么为什么网格要对数几何粒径比恒定如2^(1/3)对守恒与B_ij构造为什么重要solve_ivp积分很顺但我怎么知道它没偷偷丢掉质量篇 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 批式底账 C2dN_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 铺垫真实工业磨是连续流Bilgili2024给连续全混式d(H·M_i)/dt f_i − p_i − (S_i/d_i)·p_i Σ_{ji} b_ij·(S_j/d_j)·p_jp_i d_i·H·M_i。其中d_i是排料速率若有粒度分级效应可用归一化排料速率d_i*τ·d_i充当平均分级函数——这就是分级进入粉碎模型的接口为篇 11 的筛网和篇 12 的 Tromp 埋线。本篇先从更干净、更好验算的批式落地连续全混式留到篇 20 完整项目。为什么必须固定枢轴碎片不落在节点上数值上我们不能连续存粒径只能用一组离散网格。难点在出生项母粒j碎出的子粒体积是连续分布的不可能恰好落在某个网格节点上。若强行把落在两节点间的子粒分给两端需要既不偷质量、也不伪造数量密度——Kumar–Ramkrishna1996Fixed Pivot正是干这个对每个离散网格固定枢轴选择一组约束通常保留两个矩如总质量/总体积与总数把任一落点分配给它两个相邻枢轴使这两个矩在离散化前后保持。对破碎这种只关心子粒体积和 母粒体积零阶一阶矩守恒即可的场合固定枢轴法能构造出精确的质量守恒的B_ij。序列级术语提示Hounslow et al.1988/1990提出四倍粒径几何网格与差分格式Kumar–Ramkrishna1996的 Fixed Pivot 是更一般地把连续 PBE 换成离散 ODE 并保矩的方案——它就是本系列落地破碎 PBM 的首选底账 G可靠性优于直接四倍网格的启发式差分。网格共轭对数几何 bin→守恒的关键铁律 5 规定对数等距粒径比恒定。原因有三①制药 PSD 常跨 1.5 个数量级对数网格稀疏处不浪费 bin②固定粒径比让 Bini的平均体积与i1满足常数比值出生项的分配向量可预计算③质量守恒核验用Σ distrib·Δx时几何网格的Δx随档位变化正好把每档占多少体积体现进来。本系列取x_grnp.geomspace(10,2000,num40)工程选型 30~50 bin底账 G。质量守恒就是验收试金石数值格式是不是守恒、积分有没有发散唯一硬判据是总质量不回飘批式破碎下母粒子破碎不产生也不消灭质量所以积分结束后Σ M_i质量份额应与初值一致。这也是固定枢轴法 网格共轭 每步验算三者缺一不可的原因——任一处出错质量守恒断言立刻暴露。这也是本系列相对跑出来好看的普通教程最硬的纪律一切数值结论都先过守恒这一关。二、完整代码与逐行剖析下面用真实可运行的示例实现Hounslow 几何网格 → 一套S_i承接篇 9 的 Vogel–Peukert→ 固定枢轴法构造B_ij→solve_ivp(methodBDF)积分 → 质量守恒核验。S_i/B_ij的具体参数化是篇 11 的正文这里先用一个幂律S_i与等质量双子碎B_ij占位重点演示数值框架与守恒纪律。importnumpyasnpfromscipy.integrateimportsolve_ivp# ---- 1) Hounslow 对数几何网格固定粒径比铁律5----x_lo,x_hi,numbins10.0,2000.0,40ratio(x_hi/x_lo)**(1.0/(numbins-1))# 固定粒径比x_grx_lo*ratio**np.arange(numbins)# 节点x_cnp.sqrt(x_gr[:-1]*x_gr[1:])# bin 中心几何中点dxx_gr[1:]-x_gr[:-1]# bin 宽nbnumbins-1# bin 数# ---- 2) 选择函数 S_i承接篇9 概率核化作逐档速率 (示意参数) ----# 简单幂律S_i ∝ (x_c / x_ref)^qdefselection(x_c,k01.0,q0.8):returnk0*(x_c/500.0)**q Sselection(x_c)# ---- 3) 破碎函数 B_ij固定枢轴法构造双子均分保体积/质量守恒 ----# B_ij b_{i-1,j} 形式母粒 j 恰好碎成两个体积均分的子粒# 落在两相邻枢轴 (i, i1) 间按线性分配使其体积和母粒。# 这里做示意实现子粒等分到最接近的两个 bin保质量份额。defbreakage_fixed_pivot(x_c):nx_c.size Bnp.zeros((n,n))volx_c**3.0# 以体积代表质量份额(密度均匀)forjinrange(n):v_child0.5*vol[j]# 子粒体积(等分2粒)# 找到 v_child 落点邻近最大枢轴iint(np.searchsorted(vol,v_child))# 第一个 volv_childimax(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 Bbreakage_fixed_pivot(x_c)# ---- 4) PBM ODEdN_i/dt Σ_j S_j B_ij N_j - S_i N_i ----defpbm_rhs(t,N):birthnp.dot(B,S*N)# Σ_j B_ij * (S_j N_j)deathS*Nreturnbirth-death N0np.zeros(nb);N0[8]1.0# 初始单档(50µm附近)一批母粒M0np.sum(N0*vol)# 初始质量份额# ---- 5) 积分 质量守恒核验 ----t_end50.0solsolve_ivp(pbm_rhs,(0.0,t_end),N0,methodBDF,dense_outputTrue,rtol1e-8,atol1e-10)Nfsol.y[:,-1]Mfnp.sum(Nf*vol)driftabs(Mf-M0)/M0# 相对质量漂移print(f初带档: bin8(x_c{x_c[8]:.1f}µm) 末带质量份额相对漂移 {drift:.2e})assertdrift1e-6,质量守恒失败请检查 B_ij 构造# ---- 行为核验粗粒应当变少、细粒变多纯破碎----coarse0np.sum(N0[x_c800.0]*vol[x_c800.0])coarse1np.sum(Nf[x_c800.0]*vol[x_c800.0])print(f800µm 质量份额: 初{coarse0:.4f}- 末{coarse1:.4f})assertcoarse1coarse0,纯破碎不应出现粗粒反而增多print(OK质量守恒 纯破碎右移行为均通过)逐行要点ratio(x_hi/x_lo)**(1/(n-1))保证固定粒径比几何网格满足铁律 5x_c用几何中点承载 bin 体积。selection用幂律S_i ∝ (x_c/500)^q占位示意参数篇 11 才做参数化它承接篇 9 的vogel_peukert_breakage_probability——本篇用幂律是先搭守恒框架篇 11 将用物理核替换。breakage_fixed_pivot把v_childvol[j]/2双子均分落在两枢轴间按体积差w1(vol[i]−v_child)/(vol[i]−vol[i−1])分配恰好保体积守恒——B满足Σ_i B_ij 1每列和为 1这是守恒的代数保证。pbm_rhs里birthnp.dot(B, S*N)用矩阵乘法一次性算完全部出生项避免逐 bin 双重循环——性能上这是工程关键篇 19 会谈向量化批量。solve_ivp(..., methodBDF)刚性系统破碎速率差别可达几个数量级用隐式 BDF 更稳rtol/atol收紧保证守恒精度。双断言质量漂移1e-6守恒底线与粗粒质量减少破碎右移行为——把对错变成可自动判定的标准。三、常见报错与排查质量守恒断言失败漂移很大。现象drift达到1e-2及以上。根因绝大多数在B_ijbreakage_fixed_pivot的分配索引i越界/落点分配不合理导致Σ_i B_ij ≠ 1。解法先单独打印B.sum(axis0)每列和必须全为 1再检查i的 clamp 逻辑imax(1,min(n-1,i))是否把边界档处理对。用methodRK45积分发散/NaN。现象solve_ivp报步长过小或出现 NaN。根因破碎系统是刚性显式 RK45 处理不了跨数量级的速率差或无足够精度。解法改用methodBDF或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_c800的初值只有 bin8 一档初就极少。解法练习里建议用一个宽的平整初值如N0np.ones(nb).copy(); N0[x_c100]0之类再验右移行为断言才有统计意义。四、动手练习守恒自验必做对多个不同N0单档、平带、偏细、偏粗各跑一遍求drift。判定所有情况下drift1e-6证明你的B_ij与solve_ivp配置在不同初值下都守恒。网格收敛把numbins从40试到80/160观察同一t_end下 PSD 的d50用篇 14 的算法是否收敛、drift是否随网格加密而改善。判定d50变化5%且drift持续1e-6证明网格密度已够。把篇 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_ijsolve_ivp(methodBDF)积分批式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。本篇认知问题回显FAQQ1破碎 PBM 的出生、死亡项离散后长什么样A1对档位idN_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_ij1守恒代数保证。Q2母粒碎片不在网格节点上怎么处理A2用 Kumar–Ramkrishna1996固定枢轴法把碎片子粒按其体积/质量落点分配给两个相邻枢轴按体积差线性加权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(axis0)1的代数检查以及methodBDF刚性 收紧rtol/atol就能把积分数值误差压到守恒线以下。Q5篇 9 的S(x)概率核怎么变成速率系数S_i分结束后Σ N_f·vol / Σ N_0·vol应等于 1 到很高精度本文1e-6。配合B.sum(axis0)1的代数检查以及methodBDF刚性 收紧rtol/atol就能把积分数值误差压到守恒线以下。Q5篇 9 的S(x)概率核怎么变成速率系数S_iA5本篇第一步用幂律S_ik0·(x_c/500)^q示意参数搭守恒框架真正物理的做法是把vogel_peukert_breakage_probability(x_c*1e-6, W_m_kin)映射成逐档S_i练习 3 已打通。篇 9 给碎与不碎的概率篇 10 把它乘以单位时间碰撞强度变成单位时间破碎速率。