用于揭示均匀冷却状态下粗糙颗粒流控制方程的非线性动力系统规范化稀疏识别方法

《Powder Technology》：Normalized sparse identification of nonlinear dynamics system for discovering the governing equation of rough granular flow in the homogeneous cooling state

　　•提出了一种可解释的机器学习方法——归一化SINDy。•该方法旨在揭示颗粒流的控制方程。•该途径有效解决了字典矩阵的病态问题。•通过数据得出了HCS中粗糙颗粒流的控制方程，这些数据是通过离散元法生成的。•所发现的控制方程在多种工况下与颗粒流动力学理论相符。•这些控制方程能够有效描述T(t)和R(t)的变化规律。

引言
自然界中的泥石流、制药生产以及食品的生产与运输都与颗粒流系统密切相关[1]、[2]。随着计算资源的日益丰富以及先进算法的发展，计算模拟已成为研究颗粒流流体动力学的主流方法，其中主要包括基于局部连续体假设的欧拉方法（如连续体模型）[3]以及拉格朗日方法（如离散元法，DEM）[4]、[5]。拉格朗日方法以牛顿运动定律为基础，通过分析作用在每个粒子上的力来研究其运动情况。在工业或自然界的颗粒流系统中，由于粒子数量庞大，计算需求极高，往往会导致成本过高。为在计算效率与精确度之间取得平衡，人们通常采用基于局部连续体假设的N-S方程的欧拉方法[6]来模拟大规模颗粒流系统，因此这种方法得到了广泛应用。然而，这些方法的精确度完全依赖于合理的粒子-粒子相互作用模型，比如固体-固体应力、固体-固体阻力以及其他出现在控制方程中的本构关系。用于确定粒子间相互作用的传统方法通常包括基于高精度拉格朗日模拟数据的直接统计平均方法[7]、[8]，以及基于统计力学原理的颗粒流动力学理论（KTGF）[9]、[10]。KTGF的理论基础是硬球模型假设，但该假设无法解释粒子间的长程相互作用[11]。此外，KTGF还假定所有粒子间的相互作用都仅通过二元碰撞发生[12]，因此传统的KTGF[9]无法体现多体碰撞效应对固体本构关系的影响[13]。统计平均方法[12]则侧重于动力学部分和碰撞部分这两种力学机制，与动力学理论类似，进而推导出一些极限情况。随着机器学习和人工智能技术的快速发展，这些方法被广泛用作发现颗粒流中各种本构关系的工具[14]、[15]、[16]。学习过程的可解释性是对这些方法进行分类时的重要因素。最普遍的分类方式是将它们分为黑箱模型和白箱模型。黑箱模型的主要特征是运用各种神经网络[17]、[18]、[19]，这类模型依赖大量数据集，并通过多种神经网络架构来学习一系列本构关系，比如固体应力张量与应变之间的关系。这类神经网络方法的学习和训练目标仅仅基于已知数据中的数值或概率相关性，而通过自动微分在神经网络中获取参数的过程仅仅是一种数学运算，无法与物理定律相联系，因此神经网络的训练过程仍然不具备可解释性。为提高学习和训练过程的可解释性，人们提出了物理信息神经网络（PINN）方法[20]、[21]，该方法将物理定律作为损失函数整合到最终优化的参数目标中，从而确保通过自动微分在神经网络中确定的参数优化过程必然符合这些物理定律。例如，在探索颗粒剪切流中固体应力张量与应变之间关系的最终目标时[22]，损失函数就会包含这种关系，神经网络中的参数优化过程必须满足这一关系，从而显著提升学习结果的精确度和泛化性能。除了黑箱模型，白箱模型也被用来揭示颗粒流系统中的本构关系，这类方法可以直接以特定形式输出本构关系，同时学习过程中考虑的其他替代方案也是已知的，这类方法主要属于回归技术范畴。通过回归方法获取本构关系的常见方法包括符号回归和稀疏回归。符号回归方法[23]、[24]、[25]能够从数据中独立提取物理定律，并能使用常见的数学符号自动生成替代项，这种方法的优点是替代项较为完整，但也存在三个缺点：首先，学习过程的计算需求很高；其次，候选模型的复杂性往往会导致过拟合现象；第三，学习得到的结果可能包含一些难以用理论解释的项。为在计算复杂度与方程发现的有效性之间取得平衡，非线性动力学的稀疏识别方法（SINDy）[26]、[27]、[28]、[29]能够结合基于物理定律的专家经验，将控制方程或本构关系中常见的项作为字典矩阵中的候选项。该方法有STRidage和TrainSTRidage算法作为支撑，从而确保学习结果的准确性和稀疏性。稀疏采样技术已被证明能够大幅减少计算工作量并提升学习速度，不过这类方法目前主要应用于气体或液体流动系统[27]、[29]、[30]，在颗粒流系统领域的应用研究相对较少。尽管如此，这些方法的基本原理仍有望被应用于颗粒流系统。早期的研究尝试将SINDy方法应用于具有光滑颗粒的均匀冷却状态下的最简单颗粒流案例[28]，该系统中的控制方程仅为描述颗粒温度随时间变化的常微分方程，主要的挑战在于确定控制方程中能量耗散源项的形式，而这项工作已经证明了能够准确得到不同情况下能量耗散源项的完整形式，同时还能推导出颗粒温度的时间变化关系。因此，本文的目标是将该方法扩展到具有粗糙颗粒的均匀冷却状态下的颗粒流系统[31]、[32]，该系统的控制方程被扩展为一组常微分方程，这些方程分别描述了颗粒平动温度/转动温度随时间的变化情况，其中平动颗粒温度与转动颗粒温度之间存在巨大差异，这增加了学习的复杂性，而且在实际颗粒流系统中，不同的基本流体动力学变量之间往往存在数量级的差异，因此这项工作有望探讨新的SINDy方法是否能够发现更复杂的颗粒流系统，从而为建立颗粒流或多相流系统的本构关系提供新的思路。

第2节的目的是全面概述用于校准具有粗糙颗粒的均匀冷却状态下颗粒流中单个粒子属性的精确方法，随后会通过高精度模拟方法结合合理的颗粒材料来生成数据集。第3节介绍了为HCS中粗糙颗粒流的控件方程而发现的归一化SINDy方法，重点关注涉及颗粒平动温度/转动温度的场景。第4节详细介绍了该新型方法在该系统不同工况下的学习结果。第5节则总结了研究成果，并阐述了未来的研究方向。

章节片段
HCS中具有粗糙颗粒的颗粒流
HCS中具有粗糙颗粒的颗粒流是指在三周期边界条件下，单个颗粒的速度最初遵循麦克斯韦分布的一种场景。由于颗粒之间的非弹性碰撞和接触摩擦，颗粒的平动颗粒温度（T）和转动颗粒温度（R）会随时间逐渐发生变化。据推测，在特定的时间内，系统内的颗粒会保持某种状态。

方法：SINDy
根据HCS中具有粗糙颗粒的颗粒流系统的特点，只有颗粒的平动温度/转动温度会随时间变化，这种变化可以用两个常微分方程来表示，其中的非稳态项为(dT/dt, dR/dt)，而唯一的不确定源项则是颗粒的平动能量耗散和转动能量耗散。参考颗粒流动力学理论[31]、[32]的知识，这些源项可能与能量耗散项存在关联。

结果与讨论
从数学角度来看，归一化SINDy方法能够完全规避字典矩阵的病态问题，通过利用伪稀疏向量??′=(ξ′,ζ′)与真正稀疏向量??=(ξ,ζ)之间的线性变换关系，可以精确地找到HCS中具有粗糙颗粒且平动温度与转动温度存在数量级差异的颗粒流的控制方程，因此本研究能够实现这一目标。

结论
本研究成功提出了一种可解释的机器学习方法——归一化SINDy，可用于精确发现具有显著不同基本流体动力学特性的流动系统的控制方程，本研究主要聚焦于具有粗糙颗粒的均匀冷却状态下的颗粒流，为此采用了颗粒流动力学理论来校准数据样本，从而在离散模拟中得到与实际情况相符的颗粒属性。

CRediT作者贡献声明
赵碧丹：写作——审稿与编辑，写作——初稿撰写，可视化，验证，监督，软件，资源，项目管理，方法论，研究，资金获取，正式分析，数据整理，概念构建。宋佳超：写作——初稿撰写，可视化，软件，方法论，正式分析，数据整理。荆文正：写作——审稿与编辑。王俊武：写作——审稿与编辑，项目管理，方法论，资金获取。

利益冲突声明
作者声明不存在任何可能影响本文研究的已知财务利益或个人关系。

致谢
本研究得到了中国国家自然科学基金（22578485、22378399、22478421）、中国石油大学（北京）科学基金（2462024YJRC012、2462024YJRC008）以及中国国家重点研发计划（煤炭重大专项2025ZD1701403）的资助。

赵碧丹|宋佳超|荆文正|王俊武