关于分析方法及其在求解多变量断裂种群平衡方程中的效率的说明

《Advanced Powder Technology》：A note on the analytical approach and its efficiency for solving multivariate breakage population balance equations

　　普拉克拉蒂·库什瓦|吉特拉吉·萨哈|斯特凡·海因里希印度阿萨姆邦古瓦哈提印度理工学院数学系，邮编781 039摘要在粉末技术中，破碎平衡方程具有极其重要的意义，尤其适用于那些涉及颗粒尺寸减小、研磨、造粒以及材料精细处理的工艺过程。破碎过程的动态特性可以通过多变量公式很好地描述，但由于强烈的非线性以及内部和外部坐标的影响，这类方程很难求解。为了解决这些难题，本文采用基于函数同伦的半解析技术，在多变量框架下求解线性及碰撞型破碎方程。所提出方法的准确性和计算效率已通过与阿多米安分解法和有限体积法的比较得到验证。数值结果表明，相较于现有方法，基于同伦的方法在计算精度方面更具优势。因此，同伦分析方法可被视为解决科学研究和工业应用中复杂多变量颗粒问题的强大工具。引言破碎平衡方程为描述在各种物理相互作用过程中因碎裂而形成的颗粒尺寸分布提供了数学框架[1]、[2]。在破碎过程中，颗粒会不断受到机械应力或颗粒间相互作用的影响而发生断裂或分裂。这一理论在粉末技术中有着广泛的应用，因为控制颗粒的尺寸和分布对于保证产品质量和工艺效率至关重要[3]、[4]、[5]。例如，在化工和食品行业中，产品的反应性和溶解度都取决于颗粒尺寸[6]；而在制药行业，均匀的颗粒尺寸则能确保药物剂量的稳定性[7]、[8]。传统上，破碎平衡方程多为单变量形式，仅追踪单一颗粒属性（如尺寸、体积或质量）。此类模型已在文献中被广泛研究并应用[9]、[10]。然而，颗粒的尺寸、形状、成分、含水量等属性往往同时发生变化[11]，而这些属性又会影响破碎的概率以及颗粒的破碎方式。比如，在结晶过程中，颗粒的形状和纯度都很重要[12]、[13]；在造粒过程中，颗粒的孔隙率和成分也会影响其强度和破碎行为[14]、[15]。多变量破碎平衡方程的应用对工业领域有着深远影响，它使得人们对许多制造系统中的颗粒过程进行精确建模、模拟和控制成为可能。制药、化工、聚合物、食品等行业都能利用这一能力来预测和控制多种颗粒属性，从而改善产品设计、提升产品质量，同时还能显著降低生产成本，减少因昂贵且耗时的实验而产生的产品浪费。因此，开发能够纳入多种颗粒属性的多变量平衡模型就显得极为重要，这样才能更准确地描述相关工艺过程。这类模型有助于研究人员和科学家优化控制工艺，进而提升产品的整体性能。根据其背后的物理机制，破碎过程大致可分为两类：(i)外部诱导破碎，即颗粒因流体动力、热力或机械力等外部应力而发生碎裂[16]；(ii)碰撞型破碎，即颗粒因彼此间的相互作用而发生碎裂[17]、[18]。这两种机制的根本区别在于引发碎裂过程的驱动力的性质。在前者中，颗粒在外部力量的作用下独立破碎；而在后者中，破碎则是由于颗粒之间的碰撞所导致的。因此，外部诱导破碎通常用线性平衡方程来建模，而碰撞型破碎则因涉及颗粒间的相互作用项而需用非线性公式来描述。这类破碎指的是颗粒受到流体剪切力、湍流、温度梯度或机械载荷等外部应力作用而发生的碎裂现象。在这种条件下，每个颗粒的碎裂都是独立于其他颗粒的，其碎裂过程完全由该颗粒的基本属性以及所受的外部场决定。发生碎裂的颗粒被称为母颗粒，而由此产生的碎片则称为子颗粒。由于每次碎裂过程只涉及一个颗粒，因此相应的平衡方程是关于颗粒数密度函数的线性方程[19]。这类破碎过程因此被称为外部诱导型或线性破碎（见图1）。对于任意整数d?1，设?R+d?0,∞d为d变量实空间，其元素为?x→?{x1,x2,…,xd}，且对所有r=1,2,…,d都有xr?0。在多变量框架中，x→的各个分量被称为内部坐标，它们代表颗粒的尺寸、含水量、形状因子、熵值等属性。在此背景下，数密度函数n(x→,t)表示时间t?0时颗粒属性x→的分布情况。多变量线性破碎平衡方程可表示为[19]、[17]（见图2）：?n(x→,t)?t=∫x→∞→b(x→,y→)S(y→)n(y→,t)dy→-S(x→)n(x→,t)，该方程还包含初始条件n(x→,0)=n0(x→)?0，对所有x→∈R+d均成立。此处，我们用∫x→∞→dy→表示∫x1∞∫x2∞…∫xd∞dyd…dy2dy1。在方程(1.1)中，破碎函数b(x→,y→)表示由于母颗粒具有属性y→而产生的、具有属性x→的子颗粒的分布情况。当x→0，而当x→?y→时，b(x→,y→)=0，这意味着只有当母颗粒y→的各属性都大于子颗粒x→的对应属性时，才会发生破碎现象。显然，如果x→的各属性大于y→，则不会发生破碎。函数S(y→)则表示母颗粒y→被选中进行破碎的速率。此处，n、b和S都是关于其自变量的非负函数。一般来说，破碎函数具备以下性质：1.表达式ν(y→)表示在一次破碎过程中产生的子颗粒总数，即ν(y→)=∫0→y→b(x→,y→)dx→。这一条件体现了破碎过程的守恒定律，即粒子数量保持不变。积分域[0,y→]反映了物理上的约束，即每个子颗粒在各个坐标上的数值都必须小于母颗粒。此外，ν?2意味着母颗粒会碎裂成两个或更多颗粒。2.多变量建模中的守恒定律有两种类型[11]，即：a)属性守恒定律：定义f(x→)=∑r=1dxr，那么f(y→)=∫0→∞→f(x→)b(x→,y→)dx→；b)超体积守恒定律：定义f?(x→)=∏r=1dxr，那么f?(y→)=∫0→∞→f?(x→)b(x→,y→)dx→。为更清晰地理解这些守恒定律的物理含义，以干燥造粒过程为例。在该过程中，如果蒸发可以忽略不计，那么母颗粒的总质量和含水量在碎裂后依然保持不变。在平衡方程中，可通过将含水量设定为一个内部坐标，并确保子颗粒的总含水量与母颗粒相等来实现这一效果。需要注意的是，当d=1时，上述方程就简化为标准的体积守恒定律。当两个颗粒发生碰撞时，任一颗粒的破碎概率以及碎片的属性分布都会取决于碰撞颗粒的特定属性，如尺寸、结构、材料特性以及相对速度。这类破碎过程可能导致三种不同的结果：(i)破坏性碎裂，即两个碰撞颗粒都会发生碎裂；(ii)属性传递，即碰撞产生的碎片可能会附着在另一个颗粒上，从而导致该颗粒的部分生长或聚集；(iii)单颗粒碎裂，即仅有一个颗粒因受到另一个颗粒的撞击而发生碎裂。由于颗粒间的相互作用直接决定了碎裂的动态过程，因此这种情况比“线性”情况更为复杂，相关内容可在文献[18]、[20]中找到。在本研究中，我们考虑的是其中一个碰撞颗粒因受到撞击而发生碎裂，而另一个颗粒则保持不变。这样一来，其中一个碰撞颗粒就是产生子碎片的母颗粒，而另一个颗粒则起到催化剂的作用，提供必要的碰撞冲击力，自身却不发生碎裂。碰撞诱导型破碎在数学上可通过以下积分-偏微分方程来表示[17]、[18]：?n(x→,t)?t=∫0→∞→∫x→∞→K(y→,z→)B(x→|y→;z→)n(y→,t)n(z→,t)dy→dz→-n(x→,t)∫0→∞→K(x→,y→)n(y→,t)dy→，该方程还包含初始条件(1.2)。此处，n(x→,t)的含义与线性破碎情况相同。K(x→,y→)表示具有属性x→和y→的颗粒之间的碰撞率。一般来说，碰撞核函数是对其自变量对称且非负的。B(x→|y→;z→)是碰撞破碎分布函数，它表示由于颗粒y→与颗粒z→发生碰撞，从而产生具有属性x→的颗粒的速率。与线性破碎情况类似，碰撞破碎函数B(x→|y→;z→)也具备以下性质：1.如果对于所有i=1,2,…,d都有yi>xi，那么B(x→|y→;z→)>0；2.碎裂后产生的子颗粒的平均数量为ν?(y→,z→)=∫0→y→B(x→|y→;z→)dx→；3.守恒定律与线性碎裂情况的扩展形式类似，分别为：a) f(y→)=∫0→∞→f(x→)B(x→|y→;z→)dx→，b) f?(y→)=∫0→∞→f?(x→)B(x→|y→;z→)dx→。除了颗粒属性分布之外，不同阶的矩也在描述平衡系统的动态特性方面起着重要作用，因为它们代表了不同的物理量。为定义多变量平衡方程的矩，可设m=∑r=1dmr，其中mr为非负整数。n(x→,t)的m阶矩可表示为Mm1m2…md(t)=∫0→∞→∏r=1dxrmrn(x→,t)dx→。与单变量情况类似，零阶矩M00…0表示颗粒的总数。第一阶矩M00…1…0(t)（第r位为1）则表示xr属性的总含量。因此，共有d个一阶矩，为保持颗粒属性总量的守恒，就需要保持所有一阶矩的守恒。此外，M11…1（每个分量均为1）表示系统的总超体积，它用于描述内部属性的累积量。例如，在二维情况下，如果内部坐标为高度和宽度，那么它们的超体积就代表面积；在三维情况下，当高度为第三个内部坐标时，超体积则代表颗粒的总体积。关于超体积的更多细节可参考文献[11]。根据破碎平衡方程的应用场景，不同阶矩的守恒也有不同的要求。在某些情况下，超体积守恒至关重要，而在其他情况下，整体属性的守恒则对成功应用尤为关键。设?φ(x→)?∏r=1dxrmr，那么线性破碎情况下m阶矩随时间变化的一般方程可表示为dMm1m2…md(t)dt=∫0→∞→∫0→y→φ(x→)b(x→,y→)dx→-φ(y→)S(y→)n(y→,t)dy→。线性破碎情况下零阶矩随时间的变化则可表示为dM00…0(t)dt=∫0→∞→S(y→)n(y→,t)[ν(y→)-1]dy→。类似地，对于碰撞型破碎，m阶矩随时间的变化则为dMm1m2…md(t)dt=∫0→∞→∫0→∞→∫0→y→φ(x→)B(x→|y→;z→)dx→-φ(y→)K(y→,z→)n(y→,t)n(z→,t)dy→dz→。对于零阶矩，上述方程可简化为dM00…0(t)dt=∫0→∞→∫0→∞→[ν?(y→,z→)-1]K(y→,z→)n(y→,t)n(z→,t)dy→dz→。章节节选文献综述与研究动机在现有文献中，已有不少文章致力于求解单变量线性及碰撞型破碎方程。由于积分-偏微分方程具有非局部特性，其精确解的求得通常仅限于选择特定的简单核函数和初始条件。然而，实际的自然和工业应用中所面临的模型则往往涉及更为复杂的核函数和初始条件。为更深入地理解和掌握这类系统，人们通常会采用数值方法。单变量方法论通常，我们认为问题的初始解估计值应与初始条件一致，即?0(x→,t)=n(x→,0)。设q∈[0,1]为嵌入参数，借助该参数，我们可以通过同伦函数H构建从n(x→,t)到?(x→,t)的连续映射，该同伦函数的定义为?H[?(x→,t);q,c0,H]?(1-q)L[?(x→,t)-?0(x→,t)]-qc0ξ(x→,t)N[?(x→,t)]=0，如此一来，随着嵌入参数q从0变化到1，解函数?(x→,t)也会从初始估计值?0(x→,t)逐渐变为精确解。数值讨论在本节中，我们根据破碎模型和变量数量对测试案例进行了分类讨论，共分析了六种案例，其中前两种属于线性破碎模型，其余四种则为碰撞型破碎方程模型。为验证所提出的方法，我们将基于HAM的解法与ADM、FVS以及精确解进行了比较。对于ADM和HAM算法，我们在数值验证时最多考虑级数解的前七项。精确解结论我们成功设计并应用了同伦分析方法，用于求解多变量线性和碰撞型破碎平衡方程。我们还推导出了级数解的通用递推公式。为验证HAM方法的准确性，我们将颗粒密度函数及其积分性质与成熟的ADM和FVS方法进行了对比。为了进行全面的数值分析，我们选择了多个变量数量不同的测试案例。本研究为数据可用性问题提供了一种新的解决方案。本研究所得结果所依赖的数据，可在合理请求下向相应作者获取。利益冲突声明所有作者均声明不存在任何利益冲突。CRediT作者贡献说明Stefan Heinrich：写作——审阅与编辑、验证、监督、项目管理。Jitraj Saha：写作——审阅与编辑、可视化、监督、资源协调、项目管理、研究工作、形式分析、概念构建。Prakrati Kushwah：写作——审阅与编辑、初稿撰写、验证、软件应用、方法设计、研究工作、形式分析、数据整理、概念构建。致谢PK感谢印度政府Anusandhan国家研究基金会提供的国家博士后奖学金（PDF/2025/005253）带来的财政支持。JS也感谢印度政府Anusandhan国家研究基金会在本研究期间通过核心研究资助（CRG/2023/001483）以及ARG-MATRICS计划（ANRF/ARGM/2025/000442/MTR）给予的支持。