企业官网建设流程全解析

PharmaPy 干法制粒二次开发教程（4）：固相与 PBM 基础设施——SolidPhase、粒径离散网格与矩表示

版本声明块

• 工具/软件：CryPTSys/PharmaPy（Phases.SolidPhase / MixedPhases.Slurry）
• 语言/环境：Python 3.9，numpy/scipy/matplotlib
• 本文目标：会用 SolidPhase 构造、读写一套颗粒粒径分布（PSD），并算出其矩与统计量，为后续自定义辊压/破碎单元奠定状态容器基础。

一句话结论：PharmaPy 没有独立 PSD 类，粒径分布状态内嵌于SolidPhase/Slurry：x_distrib为粒径网格[µm]、distrib为数目密度[#/m³/µm]、moments为矩(um^n)，并同时支持离散 bin（1D-FVM）与矩（moments）两种表示；官方结晶网格用np.geomspace(1, 1500, num=35)对数等距。

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

• 为什么 PharmaPy 不单独提供一个"PSD 类"，而是把粒径状态塞进SolidPhase/Slurry里？
• 同一个粒径分布，凭什么要同时保留"离散 bin"和"矩"两种表示？它们各自在什么时候用？
• x_distrib、distrib、moments三个参数到底是什么维度、什么单位？读错一个会不会导致天差地别的结果？
• 官方示例那句np.geomspace(1, 1500, num=35)到底在做什么，对干法制粒二次开发我该选什么样的网格？
• 怎么从一个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 矩（moments）

SolidPhase同时携带两套互相关联的粒径描述：

1. 离散表示：x_distrib（内部网格坐标，单位 µm）+distrib（数目密度分布#/m³/µm）。这是把[xmin⁡,xmax⁡][x_{\min},x_{\max}][xmin​,xmax​]用若干bin切分，每个 bin 存单位粒径间隔上的颗粒数密度。对 PBE 离散求解（method='1D-FVM'，即一维有限体积法）而言，离散状态就是积分变量——破碎/生长的出生-死亡项逐 bin 累加。官方结晶器BatchCryst的method取值'1D-FVM'或'moments'，对应的正是选哪一套作为框内的状态变量。

2. 矩表示：moments，一组维度为 um^n 的矩数组，n=0..N−1n=0..N-1n=0..N−1，默认num_mom=4（即 μ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 是完整 PSD（d10/d50/d90、细粉），我们的粉碎 PBM 走离散表示；而跨单元的粗粒度汇总、以及求解稳定性校验，会用矩来做。

Value anchor：矩单位统一为 um^n，源码里偷偷乘1e-6

最容易被忽略、也最容易踩坑的是单位约定。底账 A4 在PharmaPy.Phases.getPropsPhaseMix内部逐字写着一行：

mom_meters=mom_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, num=35)。# 干法制粒工程选型：30~50 bin、恒定粒径比 2^(1/3)，见下注释。# ---------------------------------------------------------------x_lo,x_hi,num=1.0,1500.0,35# 单位 µmx_gr=np.geomspace(x_lo,x_hi,num=num)# 对数等距网格节点# 计算每个 bin 的中心（几何中点）与宽度x_mid=np.sqrt(x_gr[:-1]*x_gr[1:])# 几何中点 [µm]dx=x_gr[1:]-x_gr[:-1]# bin 宽度 [µm]nbins=dx.sizeprint("网格节点数",x_gr.size,"| bin 数",nbins)print("粒径比",x_gr[1]/x_gr[0],"（数值约 2^(1/3)=1.26 量级）")# ---------------------------------------------------------------# 1. 构造一条"数值密度分布"：lognormal，便于对账# 数密度分布 n(d)，d50=120 µm，几何标准差 1.6（单位 #/µm 归一）# ---------------------------------------------------------------mu,sigma=np.log(120.0),np.log(1.6)n_norm=lognorm(s=sigma,scale=np.exp(mu)).pdf(x_mid)# #/µm（未缩放）n_std=n_norm/(n_norm @ dx)# 归一，使总颗粒数=1# distrib 的单位约定是 #/m³/µm：这里以 1 m³ 取样，直接乘个数浓度 C0C0=1e12# 单位体积内颗粒数 #/m³（示例）num_density=C0*n_std# [#/m³/µm]，喂给 SolidPhase 的就是它# ---------------------------------------------------------------# 2. 构造 SolidPhase（签名底账 A4 逐字：见文中"构造签名"）# 注意：本脚本演示输入，直接用 numpy 数组，不必真的 import PharmaPy；# 真实项目里 `from pharmapy.Phases import SolidPhase` 传入即可。# ---------------------------------------------------------------# SolidPhase(path_thermo, temp=298.15, temp_ref=298.15, pres=101325,# mass=0, mass_frac=None, moments=None, num_mom=4,# distrib=None, x_distrib=None, distrib_type='vol_perc',# moisture=0, porosity=0, mole_conc=None, kv=1)distrib_in=num_density# 数目密度 [#/m³/µm]moments_in=None# 让 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。"""x=np.sqrt(xd[:-1]*xd[1:])# bin 中心dx=xd[1:]-xd[:-1]moms=np.array([np.sum(x**p*d*dx)forpinrange(mom_num)])returnmoms mom_num=4mom=getMoments(x_gr,distrib_in,mom_num)print("矩 mom[mum] =",mom)# μ0..μ3，um^0..um^3mu0=mom[0]# 总颗粒数密度print("μ0(总数密度)=",mu0)print("数均粒径(um)=",mom[1]/mu0)# μ1/μ0# 反算体积加权 d43 作为粒径量级参考d43=mom[3]/mom[2]ifmom[2]>0elsefloat("nan")print("体积矩比 d43(um)=",d43)# ---------------------------------------------------------------# 4. 由体积加权分布求 d10/d50/d90（PSD 统计基础，篇14 详述）# 体积与数目关系 w(d) ∝ n(d)·d³ ；用累积体积分位点定义 d10/d50/d90# ---------------------------------------------------------------# 注意：分布已归一（C0=1e12 是常数倍数，对分位点无影响），直接质量化w=num_density*x_mid**3# ∝ 体积权重（质量分布）cw=np.cumsum(w)/w.sum()# 归一累积质量q10,q50,q90=np.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,ax=plt.subplots(figsize=(8,4.6))ax2=ax.twinx()# 右轴：质量分布ax.plot(x_mid,num_density,ls="--",lw=1.4,color="#3a7bd5",label="数目密度 n(d) [1e12·#/m³/µm]")ax2.plot(x_mid,w/w.max()*100,ls="-",lw=1.8,color="#e67e22",label="质量分布 w(d) [归一化,%]")forq,nmin[(q10,"d10"),(q50,"d50"),(q90,"d90")]:ax2.axvline(q,color="#999",ls=":",lw=1)ax2.text(q,8,nm,rotation=90,fontsize=9,color="#444",ha="left",va="bottom")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)")lines=ax.get_lines()+ax2.get_lines()ax.legend(lines,[l.get_label()forlinlines],loc="upper right",frameon=False)plt.tight_layout();plt.savefig("figs/fig08_psd.png",dpi=150);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, num=35)（对数等距、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默认返回节点数组，长度为 num，bin 数 = num−1），distrib是每个 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_mom=4（μ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幂次。

四、动手练习

1. 双表示守恒：构造任意 lognormal 数目密度，分别用"离散 bin 积分"和显式矩公式（对数正态 3 阶矩的解析式）算 μ0、μ1、μ2、μ3，比较两者相对误差。判定：相对误差< 1e-6即通过。
2. PSD 面积归一：对给定的distrib/dx，验证Σ n_i·dx_i = μ0（总数目密度守恒），并把质量分布w∝n·d³归一化到Σw=1。判定：两者偏差< 1e-9即通过。
3. d 分位数对账：把生成的d10/d50/d90用 scipy 逆正态 / 累积积分核对，确认与输入 lognormal 参数在 1% 内一致；并验证span=(d90−d10)/d50随 sigma 增大而单调增大。判定：分位数偏差<1%、span 单调。

五、小结与下一篇预告

本篇落地了干法制粒二次开发的状态容器层：敲开了SolidPhase/Slurry内嵌 PSD 的机制，分清离散 bin（1D-FVM）与矩（moments）两种表示及互转方法，锁死了 µm/um^n 单位约定与1e-6幂次这个最容易翻车的地方，并用np.geomspace给出网格选型与一段真实可运行的对账脚本。从现在起，"颗粒长什么样"这个问题的载体有了——下一站就是把这些容器放进一个真正的、我们亲手写的单元里。

下一篇：篇 05《如何扩一个自定义单元操作（基座）》——模仿Crystallizers.py的结构，写一个最小的自定义单元类（不炸框架），让SolidPhase状态的初始化、演变与读取在SimulationExec的 flowsheet 里流转求解。到那时，前面的 PSD 容器才真正接入机器。

本篇认知问题回显（FAQ）

Q1：为什么 PharmaPy 不单独做一个 PSD 类，而把粒径状态塞进SolidPhase/Slurry？
A1：因为粒径分布必须有宿主——它不属于任何单个公式，而属于"某段物料的状态"。把x_distrib/distrib/moments绑在相对象上，跨单元传料时相对象一传，PSD 跟着走，flowDiagram 拓扑自然成立。纯固相用SolidPhase，浆料用MixedPhases.Slurry及其流束SlurryStream。对二次开发也友好：新增辊压/粉碎/分级单元只需"进料给一个带 PSD 的 SolidPhase、出料重写之"，不必重做状态容器。

Q2：为什么同一种 PSD 要保留离散 bin 和矩两种表示？
A2：它们各有所长。离散 bin 是逐 bin 的数目密度distrib，用一维有限体积法（method='1D-FVM'）逐 bin 演算破碎/生长的出生-死亡项，能保留完整分布形状（干法制粒 CQA 要完整 d10/d50/d90）；矩方法是分布的压缩投影（μ0 总数、μ1/μ0 平均粒径、μ3/μ2·体积矩比 d43 等），状态数少、省算，但丢失形状细节，方程常不封闭。官方BatchCryst的method就在'1D-FVM'与'moments'之间选，两者用getMoments/getDistribution/convert_distribution互转。

Q3：x_distrib、distrib、moments到底是什么维度单位？
A3：x_distrib是内部粒径网格坐标，单位 µm（通常是网格节点数组）；distrib是数目密度分布，单位#/m³/µm（每立方米料内、每微米粒径间隔的颗粒个数，逐 bin 一个值）；moments是矩数组，维度 um^n（n=0…N−1），默认num_mom=4。三者统一走 µm 约定。

Q4：官方那句np.geomspace(1, 1500, num=35)在做什么？干法制粒该选什么网格？
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/d90，span=(d90−d10)/d50。矩守恒强调两个硬约束：μ0（总颗粒数密度）在多步演化里只随破碎/生长按物理规律变化、不因数值离散漂移;以及离散积分∫n·dx=μ0（双表示自洽）、质量分布归一Σw=1（w∝n·d³）二者必须成立，是校验 PBM 实现是否可信的基石。