PharmaPy 干法制粒二次开发教程4固相与 PBM 基础设施——SolidPhase、粒径离散网格与矩表示版本声明块工具/软件CryPTSys/PharmaPyPhases.SolidPhase / MixedPhases.Slurry语言/环境Python 3.9numpy/scipy/matplotlib本文目标会用 SolidPhase 构造、读写一套颗粒粒径分布PSD并算出其矩与统计量为后续自定义辊压/破碎单元奠定状态容器基础。一句话结论PharmaPy 没有独立 PSD 类粒径分布状态内嵌于SolidPhase/Slurryx_distrib为粒径网格[µm]、distrib为数目密度[#/m³/µm]、moments为矩(um^n)并同时支持离散 bin1D-FVM与矩moments两种表示官方结晶网格用np.geomspace(1, 1500, num35)对数等距。〇、本篇要解决的认知问题为什么 PharmaPy 不单独提供一个PSD 类而是把粒径状态塞进SolidPhase/Slurry里同一个粒径分布凭什么要同时保留离散 bin和矩两种表示它们各自在什么时候用x_distrib、distrib、moments三个参数到底是什么维度、什么单位读错一个会不会导致天差地别的结果官方示例那句np.geomspace(1, 1500, num35)到底在做什么对干法制粒二次开发我该选什么样的网格怎么从一个SolidPhase里取出 d10/d50/d90 这类统计量矩的第零/第一阶为什么要守恒一、机制解析相本身内嵌 PSD而不是独立类PharmaPy 的对象体系跟许多工艺仿真工具不同没有一个孤立的ParticleSizeDistribution类。相反粒径的完整状态网格、数目密度、矩作为成员变量挂在相对象上——PharmaPy.Phases.SolidPhase纯固相、PharmaPy.MixedPhases.Slurry液-固浆料及其流束SlurryStream。SlurryStream是浆料越过单元边界时携带的移动态SolidPhase则是驻留在单元如结晶器里的存档态。这个设计的理由很朴素粒径分布不是独立存在的量它必须有宿主。一个粉碎机关心的是进来的颗粒长什么样进料/给料相所含 PSD一个辊压机关心的是我吐出的料带长什么样固相/j 单元内 Phytomem。把 PSD 绑在相上跨单元传递物料时就只传相对象SlurryStream/LiquidStream/SolidPhase粒径信息跟着相走flowDiagram 拓扑串起来就自然成立——这正是篇 03 flowDiagram 里每个单元只关心挂进来的 Phases/Inlet的原因。对干法制粒二次开发而言这一抉择是友好的我们要新增的辊压机/粉碎机/分级机只需要进料持有一个 PSD、出料重写这个 PSD用现有SolidPhase作宿主即可无需重写状态容器。“双表示”离散 bin有限体积vs 矩momentsSolidPhase同时携带两套互相关联的粒径描述离散表示x_distrib内部网格坐标单位 µmdistrib数目密度分布#/m³/µm。这是把[xmin⁡,xmax⁡][x_{\min},x_{\max}][xmin​,xmax​]用若干bin切分每个 bin 存单位粒径间隔上的颗粒数密度。对 PBE 离散求解method1D-FVM即一维有限体积法而言离散状态就是积分变量——破碎/生长的出生-死亡项逐 bin 累加。官方结晶器BatchCryst的method取值1D-FVM或moments对应的正是选哪一套作为框内的状态变量。矩表示moments一组维度为 um^n 的矩数组n0..N−1n0..N-1n0..N−1默认num_mom4即 μ0、μ1、μ2、μ3。矩是粒径分布的压缩投影μ0 是总颗粒数、μ1/μ0 是平均粒径、μ2 关联比表面积、μ3 关联总体积/质量。矩方法的优势是方程数少、省算不用在几十个 bin 上逐步积分代价是分布形状信息被压缩重建回完整 PSD 需要假设分布族。MOM 方程通常不封闭μ4\mu_4μ4​被截断QMOM/DQMOM 用求积点近似求封闭——这些是篇 14 与篇 10 会深挖的值本期只需认得矩压缩投影。于是设计上就有了明确分工要逐 bin 演算破碎/辊压的完整 PSD —— 用离散表示 1D-FVM只要算平均粒径、总颗粒数等粗指标 —— 用矩表示两者用getMoments/getDistribution/convert_distribution互转。干法制粒的核心 CQA 是完整 PSDd10/d50/d90、细粉我们的粉碎 PBM 走离散表示而跨单元的粗粒度汇总、以及求解稳定性校验会用矩来做。Value anchor矩单位统一为 um^n源码里偷偷乘1e-6最容易被忽略、也最容易踩坑的是单位约定。底账 A4 在PharmaPy.Phases.getPropsPhaseMix内部逐字写着一行mom_metersmom_solid*(1e-6)**conv_exp它把以 µm 为基准的矩换算成以米为基准的矩去算体积分数。换句话说矩moments的维度约定是 um^n单位微米而非米。这与x_distrib的 [µm] 约定、与前面distrib_type的vol_perc/mass_frac语义是自洽的——整个 PSD 链统一走 µm。这个1e-6的幂次conv_exp在实际仿真里非常容易左右结果如果你把moments按米SI喂进去mom_meters会凭空差出(1e−6)n(1e-6)^n(1e−6)n个数量级体积分数随之荒谬。铁律任何传给SolidPhase/Slurry的粒径相关量先统一到 µm/um^n凡是从结果反推物理量先想清楚这个1e-6还在不在。二、完整代码与逐行剖析下面给出一段真实可运行的脚本先构造一个对数等距的粒径网格再据此构造一个SolidPhase写入/读回distrib调用getMoments算矩并反算统计量最后用 matplotlib 画出图 08 风格的 PSD 示意图。importnumpyasnpimportmatplotlib.pyplotaspltfromscipy.statsimportlognorm# ---------------------------------------------------------------# 0. 构造粒径网格对数等距几何网格官方结晶示例同款# 官方 guide 用 np.geomspace(1, 1500, num35)。# 干法制粒工程选型30~50 bin、恒定粒径比 2^(1/3)见下注释。# ---------------------------------------------------------------x_lo,x_hi,num1.0,1500.0,35# 单位 µmx_grnp.geomspace(x_lo,x_hi,numnum)# 对数等距网格节点# 计算每个 bin 的中心几何中点与宽度x_midnp.sqrt(x_gr[:-1]*x_gr[1:])# 几何中点 [µm]dxx_gr[1:]-x_gr[:-1]# bin 宽度 [µm]nbinsdx.sizeprint(网格节点数,x_gr.size,| bin 数,nbins)print(粒径比,x_gr[1]/x_gr[0],数值约 2^(1/3)1.26 量级)# ---------------------------------------------------------------# 1. 构造一条数值密度分布lognormal便于对账# 数密度分布 n(d)d50120 µm几何标准差 1.6单位 #/µm 归一# ---------------------------------------------------------------mu,sigmanp.log(120.0),np.log(1.6)n_normlognorm(ssigma,scalenp.exp(mu)).pdf(x_mid)# #/µm未缩放n_stdn_norm/(n_norm dx)# 归一使总颗粒数1# distrib 的单位约定是 #/m³/µm这里以 1 m³ 取样直接乘个数浓度 C0C01e12# 单位体积内颗粒数 #/m³示例num_densityC0*n_std# [#/m³/µm]喂给 SolidPhase 的就是它# ---------------------------------------------------------------# 2. 构造 SolidPhase签名底账 A4 逐字见文中构造签名# 注意本脚本演示输入直接用 numpy 数组不必真的 import PharmaPy# 真实项目里 from pharmapy.Phases import SolidPhase 传入即可。# ---------------------------------------------------------------# SolidPhase(path_thermo, temp298.15, temp_ref298.15, pres101325,# mass0, mass_fracNone, momentsNone, num_mom4,# distribNone, x_distribNone, distrib_typevol_perc,# moisture0, porosity0, mole_concNone, kv1)distrib_innum_density# 数目密度 [#/m³/µm]moments_inNone# 让 getMoments 从 distrib 现算dummy{x_distrib:x_gr,distrib:distrib_in,kv:1.0}print(─*60)print(数值密度分布已就绪,distrib_in.shape)print(网格覆盖 [%.0f, %.0f] µm%d 个 bin%(x_lo,x_hi,nbins))# ---------------------------------------------------------------# 3. 从离散 distrib 算矩模拟 SolidPhase.getMoments# getMoments(x_distrib, distrib, mom_num)返回维度 um^n 的前 N 个矩# ---------------------------------------------------------------defgetMoments(xd,d,mom_num):从数目密度离散分布算矩 μ_p ∫ x^p n(x) dx。xd 单位 µm → 矩单位 um^p。xnp.sqrt(xd[:-1]*xd[1:])# bin 中心dxxd[1:]-xd[:-1]momsnp.array([np.sum(x**p*d*dx)forpinrange(mom_num)])returnmoms mom_num4momgetMoments(x_gr,distrib_in,mom_num)print(矩 mom[mum] ,mom)# μ0..μ3um^0..um^3mu0mom[0]# 总颗粒数密度print(μ0(总数密度),mu0)print(数均粒径(um),mom[1]/mu0)# μ1/μ0# 反算体积加权 d43 作为粒径量级参考d43mom[3]/mom[2]ifmom[2]0elsefloat(nan)print(体积矩比 d43(um),d43)# ---------------------------------------------------------------# 4. 由体积加权分布求 d10/d50/d90PSD 统计基础篇14 详述# 体积与数目关系 w(d) ∝ n(d)·d³ 用累积体积分位点定义 d10/d50/d90# ---------------------------------------------------------------# 注意分布已归一C01e12 是常数倍数对分位点无影响直接质量化wnum_density*x_mid**3# ∝ 体积权重质量分布cwnp.cumsum(w)/w.sum()# 归一累积质量q10,q50,q90np.interp([0.10,0.50,0.90],cw,x_mid)print(d10 %.1f µm, d50 %.1f µm, d90 %.1f µm%(q10,q50,q90))print(span (d90-d10)/d50 %.3f%((q90-q10)/q50))# ---------------------------------------------------------------# 5. 可视化图 08 风格 PSD对数粒径横轴 双曲线# ---------------------------------------------------------------fig,axplt.subplots(figsize(8,4.6))ax2ax.twinx()# 右轴质量分布ax.plot(x_mid,num_density,ls--,lw1.4,color#3a7bd5,label数目密度 n(d) [1e12·#/m³/µm])ax2.plot(x_mid,w/w.max()*100,ls-,lw1.8,color#e67e22,label质量分布 w(d) [归一化,%])forq,nmin[(q10,d10),(q50,d50),(q90,d90)]:ax2.axvline(q,color#999,ls:,lw1)ax2.text(q,8,nm,rotation90,fontsize9,color#444,haleft,vabottom)ax2.set_yscale(log);ax.set_yscale(log)ax.set_xscale(log)ax.set_xlabel(粒径 d [µm])ax.set_ylabel(数目密度 n(d) [#/m³/µm])ax2.set_ylabel(质量分布 w(d) [归一化,%])ax.set_title(PSD 双表示示意数目密度(n) 与由 n·d³ 导出的质量分布(w))linesax.get_lines()ax2.get_lines()ax.legend(lines,[l.get_label()forlinlines],locupper right,frameonFalse)plt.tight_layout();plt.savefig(figs/fig08_psd.png,dpi150);plt.show()# 质量守恒自验矩由同一 distrib 积分得到∫n·dxμ0 应恒等assertnp.isclose(np.sum(num_density*dx),mu0),discrete int vs moments 不一致print(OK离散积分 ∫n·dx μ0双表示自洽)逐行要点第 0 段建网格给了两种口径官方结晶示例np.geomspace(1, 1500, num35)对数等距、35 bin干法制粒工程选型建议30–50 bin、恒定粒径比如 2^(1/3)——这个 bin 数与粒径比经验在文献里没有权威出处属工程选择正文标注为自选值靠质量守恒自验兜底。distrib_in的维度约定是#/m³/µm除以 bin 宽度dx得到每米³料里、每微米粒径间隔颗粒个数密度与底账 A4 的说明严格对齐。getMoments是SolidPhase.getMoments(x_distrib, distrib, mom_num)的本土化实现返回 um^n 维度的数组这里用 numpy 直写以脱离框架依赖、便于对账。传入SolidPhase时由框架内同名的 getMoments 完成。d10/d50/d90 用体积/质量权重w∝n(d)d³求积分位点与细粉质量分数、体积分位点的 PSD 契约一致span(d90−d10)/d50。末行断言∫n·dx μ0是离散积分 vs 矩双表示的守恒自检呼应系列铁律 4。三、常见报错与排查distrib长度与x_distrib不一致。现象构造SolidPhase时抛维度不匹配/索引错误或 getMoments 算出的 μ0 小得离谱。根因x_distrib是网格节点np.geomspace默认返回节点数组长度为 numbin 数 num−1distrib是每个 bin 一个的值长度应为 num−1。解法要么用np.sqrt(x[:-1]*x[1:])取 bin 中心对齐 distrib要么让两者长度严格匹配写前先print(len(distrib), len(x_gr))。下标越界/IndexError。现象对distrib[i]、x_mid[i]加下标越界。根因把节点数组当 bin 数组用或循环边界取了n。解法统一用bin 中心 节点几何中点这一定义;访问边界一律[:-1]、[1:]。distrib_type混用把质量/体积分数当数目密度喂。现象distrib_type是mass_frac/vol_perc而非数目密度语义若你手里是质量分数却塞进distrib矩系统会整体错比例。根因底账 A4 明确distrib是数目密度#/m³/µm而distrib_type只是声明 mass/vol 口径的属性。解法需要质量/体积分布时用convert_distribution(...)先换算别在distrib里塞质量分数。矩次数num_mom不足。现象想算 d43 用的 μ3 算不出来或重建分布族时形状参数缺失。根因SolidPhase默认num_mom4μ0…μ3某些统计如更高阶、或 needs μ4 的半经验式超出 N−1。解法构造时显式传num_mom或用getMoments(x_distrib, distrib, mom_num)指定需要的阶数。单位 µm·1e-6 换算错体积分数荒谬。现象输出 Cfitness/体积分数大出或小出几个数量级。根因矩的维度约定是 um^n而getPropsPhaseMix内部mom_meters mom_solid * (1e-6)**conv_exp会按 µm 基准换算若你喂进米单位幂次累计失真。解法凡进SolidPhase/Slurry的粒径量先统一到 µm/um^n对结果反推物理量前先确认该1e-6幂次。四、动手练习双表示守恒构造任意 lognormal 数目密度分别用离散 bin 积分和显式矩公式对数正态 3 阶矩的解析式算 μ0、μ1、μ2、μ3比较两者相对误差。判定相对误差 1e-6即通过。PSD 面积归一对给定的distrib/dx验证Σ n_i·dx_i μ0总数目密度守恒并把质量分布w∝n·d³归一化到Σw1。判定两者偏差 1e-9即通过。d 分位数对账把生成的d10/d50/d90用 scipy 逆正态 / 累积积分核对确认与输入 lognormal 参数在 1% 内一致并验证span(d90−d10)/d50随 sigma 增大而单调增大。判定分位数偏差1%、span 单调。五、小结与下一篇预告本篇落地了干法制粒二次开发的状态容器层敲开了SolidPhase/Slurry内嵌 PSD 的机制分清离散 bin1D-FVM与矩moments两种表示及互转方法锁死了 µm/um^n 单位约定与1e-6幂次这个最容易翻车的地方并用np.geomspace给出网格选型与一段真实可运行的对账脚本。从现在起颗粒长什么样这个问题的载体有了——下一站就是把这些容器放进一个真正的、我们亲手写的单元里。下一篇篇 05《如何扩一个自定义单元操作基座》——模仿Crystallizers.py的结构写一个最小的自定义单元类不炸框架让SolidPhase状态的初始化、演变与读取在SimulationExec的 flowsheet 里流转求解。到那时前面的 PSD 容器才真正接入机器。本篇认知问题回显FAQQ1为什么 PharmaPy 不单独做一个 PSD 类而把粒径状态塞进SolidPhase/SlurryA1因为粒径分布必须有宿主——它不属于任何单个公式而属于某段物料的状态。把x_distrib/distrib/moments绑在相对象上跨单元传料时相对象一传PSD 跟着走flowDiagram 拓扑自然成立。纯固相用SolidPhase浆料用MixedPhases.Slurry及其流束SlurryStream。对二次开发也友好新增辊压/粉碎/分级单元只需进料给一个带 PSD 的 SolidPhase、出料重写之不必重做状态容器。Q2为什么同一种 PSD 要保留离散 bin 和矩两种表示A2它们各有所长。离散 bin 是逐 bin 的数目密度distrib用一维有限体积法method1D-FVM逐 bin 演算破碎/生长的出生-死亡项能保留完整分布形状干法制粒 CQA 要完整 d10/d50/d90矩方法是分布的压缩投影μ0 总数、μ1/μ0 平均粒径、μ3/μ2·体积矩比 d43 等状态数少、省算但丢失形状细节方程常不封闭。官方BatchCryst的method就在1D-FVM与moments之间选两者用getMoments/getDistribution/convert_distribution互转。Q3x_distrib、distrib、moments到底是什么维度单位A3x_distrib是内部粒径网格坐标单位 µm通常是网格节点数组distrib是数目密度分布单位#/m³/µm每立方米料内、每微米粒径间隔的颗粒个数逐 bin 一个值moments是矩数组维度 um^nn0…N−1默认num_mom4。三者统一走 µm 约定。Q4官方那句np.geomspace(1, 1500, num35)在做什么干法制粒该选什么网格A4它生成一个从 1 µm 到 1500 µm、35 个节点的对数等距几何网格粒子粒径比恒定来自官方结晶示例。干法制粒工程选型建议30–50 个 bin、恒定粒径比如 2^(1/3)并配质量守恒自检这个 bin 数与粒径比经验在文献里没有权威出处属于工程自选值正文按以官方文档/工程经验为准标注。Q5怎么从SolidPhase取出 d10/d50/d90矩的第零/第一阶为什么要守恒A5先由数目密度distrib按w∝n(d)·d³得质量/体积权重分布再对累积质量曲线求 10%/50%/90% 分位点得 d10/d50/d90span(d90−d10)/d50。矩守恒强调两个硬约束μ0总颗粒数密度在多步演化里只随破碎/生长按物理规律变化、不因数值离散漂移;以及离散积分∫n·dxμ0双表示自洽、质量分布归一Σw1w∝n·d³二者必须成立是校验 PBM 实现是否可信的基石。