CFD-DEM耦合方法已越来越广泛地应用于多相流研究，如推移质输移、流化床反应等。在CFD-DEM耦合方法中，流体运动由Euler框架下的连续方法描述，固体颗粒在考虑粒间碰撞及颗粒-流体作用的基础上由Lagrange方法追踪。由于Lagrange框架下的离散方法仅用来描述系统的微观特性，在空间和时间尺度上远小于连续性方法所描述的宏观特性[1]，因此需要准确、高效地建立宏观尺度(连续介质)与微观尺度(离散介质)间的桥梁，即流体体积分数场。

目前，按照颗粒粒径D与CFD网格尺寸L的相对关系，常用颗粒解析耦合方法(D≫L，如浸没边界法等)和颗粒不解析耦合方法(L≫D)。前者解析颗粒周围的流动细节，但计算代价高昂，仅用于极少数颗粒的计算；后者只考虑质心受力，在大规模颗粒计算中应用广泛。依据颗粒位置、尺寸与所属CFD(computational fluid dynamics)网格间的相对关系，颗粒不解析耦合方法已有PCM(particle centroid method)、DPVM(divided particle volume method)和SKM(statistical kernel method)等不同复杂程度的算法[2]。

PCM是一种简单的体积分数场算法。它将颗粒的全部体积归属于颗粒中心所在的CFD网格，当颗粒中心临近CFD网格边缘时，引起的误差可能达到50%[3]。同时，当一个CFD网格中包含多个颗粒时，全部颗粒体积累加可能大于CFD网格体积；当一个CFD网格中仅包含一个颗粒时，该颗粒体积也可能大于CFD网格体积。这种非物理的不光滑结果将导致整个CFD-DEM模拟计算的数值不稳定。因此，PCM常用于CFD网格尺寸显著大于颗粒尺寸(如L约为3~10倍D)的情形[4]。

DPVM最早由Wu等[5]提出，用于解决当颗粒覆盖多个CFD网格时PCM所引起的较大误差。DPVM将颗粒体积按其实际位置在所覆盖的CFD网格中分配。相对于PCM，DPVM能得到相对光滑的体积分数场，模拟结果的准确性有显著提升[3]。但是，由于颗粒体积在不同网格中的分配与网格形状、拓扑关系等有关，DPVM在非均匀网格中较难实现，也无法解决D＞L时非物理计算结果所引起的数值不稳定问题。

SKM将颗粒体积依据内核函数在整个计算域内进行加权分配，内核函数包含盒式分布函数(此时，SKM退化为PCM)、Gauss分布函数[6]和Johnson's SB分布函数[7]等。计算域内某一点的体积分数由不同颗粒在该点处的体积权重叠加得到[8]。在D＞L时，仍能够得到光滑的体积分数场，保证耦合计算的数值稳定性。但由于颗粒体积在整个计算域内分配，需要对每个颗粒所影响的CFD网格进行循环搜索，计算代价高昂。当实施并行计算时，因独立计算的子区域之间要通过子区域交界面的信息交互实现多个处理器的协同并发计算，颗粒影响范围往往穿越不同处理器负责区域之间的内部边界，在内部边界处截断颗粒影响范围将导致颗粒体积不满足守恒条件、产生误差可达50%。这为SKM的广泛应用，特别是在大规模并行计算中的应用带来挑战[9-10]。

综上所述，体积分数场的现有算法还存在不同类型的局限，在粒径D与网格L同量级的流动情形下表现突出[2]，为诸如湍流与泥沙颗粒相互作用研究带来明显限制。特别是为充分解析湍流结构而难以改变CFD网格尺寸L时，现行算法将对颗粒粒径范围提出严格限制。为此，本文以CFD-DEM耦合方法在湍流-泥沙问题中的应用为导向，以大规模并行计算的高效实现为研究目标，提出一种改进型SKM算法，以期适用于D~L的情形，解决简单算法(如PCM、DPVM)在D＞L时所引起的较大误差和数值不稳定性、克服传统SKM计算耗时长和并行算法较难实现等问题。

本文所采用的SKM算法的内核函数为三维Gaussian分布函数[6]：

其中: T为向量转置符号；b为Gaussian分布函数参数，表征颗粒所影响的空间范围；xi为计算域内第i个颗粒的中心位置坐标，x为空间坐标。在三维空间上，此表达式满足归一化条件：

在计算域内，空间某一点的体积分数值由不同颗粒在该点处的体积权重叠加得到：

其中: α(x)为空间位置x处的颗粒体积分数值，Np为所有影响范围包含该空间位置的颗粒数量，Vp, i为计算域内第i个颗粒的体积。

为降低SKM循环搜索每个颗粒所覆盖CFD网格导致的高昂计算代价、解决并行计算中颗粒影响范围穿越计算核心内部边界等问题，本文提出采用特征点剖分方法代替传统的网格搜索方法。

以图 1所示的均匀网格为例，特征点剖分方法的实现过程如下。

1) 确定颗粒影响范围。以颗粒i的球心(坐标xi)为中心、以搜索半径R为参数，确立颗粒空间影响范围(即搜索半径R覆盖区域的外接立方体)。

2) 布置特征点。在立方体内三维方向均匀布置特征点，每个方向上的特征点间距为该方向CFD网格的最小间距(以位于x、y平面上的特征点为例，如图 1b中的Lx、Ly)，若计算域在x、y方向的长度分别为LX、LY，则特征点在此平面上的坐标为

其中，int[·]为取整运算。

3) 确定有效特征点。特征点中心到颗粒中心距离小于搜索半径R时为有效特征点，反之为无效特征点，即有效特征点满足：

4) 标记CFD网格。通过搜索有效特征点所属的CFD网格，标记该颗粒所影响的所有CFD网格。因特征点间距为CFD网格间距最小值，因而所有被影响的CFD网格都会被标记。当不同特征点属于同一网格时，只标记一次，如图 1b所示，从而使颗粒体积守恒条件满足。以位于x、y平面上的特征点为例，若当前网格中心坐标为xmi=(xmi, ymi)，且网格在x、y方向上长度分别为Lmx、Lmy，则被标记网格满足：

其中，j，k如式(4) 所示。

5) 计算CFD网格体积权重。完成网格标记后，以颗粒中心坐标为xi，以CFD网格中心坐标为xmi，由式(1) 计算该颗粒在空间位置xmi处产生的体积权重值：

6) 对所有颗粒进行循环。按照式(3) 将所有对空间位置x处有影响的颗粒体积权重值叠加，得到该时刻空间位置x处的体积分数值，实现从Lagrange场到Euler场的信息转化。

与传统SKM方法相比，特征点剖分法具有2个优势：1) 并行计算中，特征点剖分法将网格的迭代搜索转化为空间一系列位置点的搜索，因空间位置信息不受计算区域限制，故可以穿越计算核心间的内部边界，实现跨区域并行计算，解决内部边界区域因误差而产生的数值不稳定问题；2) 显著提升计算效率，特征点剖分法在物理存储量相同的条件下可以大幅减少计算耗时，特别是在颗粒粒径大于CFD网格尺寸且计算域内颗粒数量较多时更为显著，下面从算法上对计算效率的提升进行分析。

以均匀网格为例，当一个颗粒的影响半径刚好覆盖N个CFD网格时，传统SKM算法的网格搜索过程为[6]：首先搜索颗粒中心所在CFD网格(1次循环)，然后搜索该网格的6个相邻网格(61次循环)，进而再搜索6个相邻网格各自的相邻网格(62次循环)，如此循环直至超出颗粒影响半径，则循环总次数T1为

在相同条件下，特征点剖分方法则只需要在边长为2N-1个网格的立方区域内进行搜索，则循环总次数为

由此可见，在传统SKM算法中，计算耗时随颗粒覆盖的CFD网格数成指数变化，当D＞L时，计算耗时将显著大于特征点剖分方法的三次方关系，且两者之间的差距将随颗粒尺寸增大而进一步扩大(如图 1b所示)。特别是当计算域内包含大量颗粒时(假设包含M个颗粒)，两者之间的计算耗时差距将进一步扩大为单颗粒计算耗时差距的M倍：

当M较大时，特征点剖分方法节省的计算耗时将是十分显著的。该方法提高了SKM的计算效率，也可以根据计算精度要求调整特征点分布数量。当特征点间距小于CFD网格最小尺寸时，可以保证识别全部CFD网格；当精度要求不高时，也可以减少特征点数量，以进一步提高计算效率。

在特征点剖分方法的基础上，为进一步优化算法，对颗粒进行预分类，将计算域中的颗粒分为3种类型：内部颗粒、并行颗粒和边界颗粒，如图 2所示。内部颗粒是指影响范围不与任何计算域边界或计算核心之间内部边界相交的颗粒(到任何边界距离大于搜索半径R)；并行颗粒是指影响范围与计算核心之间内部边界相交的颗粒；边界颗粒是指影响范围与计算域边界存在相交的颗粒。对不同类型的颗粒在算法上进行单独处理，将避免大量的重复运算，减少计算耗时，如边界颗粒可由Zhu和Yu[7]提出的镜像颗粒法或其他方法[11]单独处理。

在上述特征点剖分法中，R和b是关键参数，其中参数b决定算法在空间上的理论搜索范围，直接影响算法的计算效率，而参数R决定空间搜索的实际截断范围，直接影响颗粒体积的守恒性和计算误差。因此，为保证算法的效率与精度，下文结合典型案例，给出R和b的选取原则并进行检验。

选择静态颗粒排列实验进行参数R和b的率定和验证。之所以选择该实验，是因为颗粒排列问题作为数学上的经典问题存在理论上的解析解[12-13]，便于参数的率定与验证；另外算例实现相对简单，配合CFD网格尺寸及形状变化，满足对改进型SKM算法在不同场景下的应用检验要求。通过不同的算例设置，主要考察颗粒尺寸与CFD网格尺寸发生相对变化及应用不同CFD网格形态时不同算法的适用性。

静态颗粒排列数值实验包含颗粒立方排列与颗粒紧密立方排列2个算例(排列形式如图 3a和3b所示)，计算整个域内的空间体积分数场，取计算域中线的数值计算结果与理论解析解[12-13]进行对比，前者用于率定模型参数，后者用于验证。其中：1) CFD网格尺寸L与颗粒粒径D的特征比值L/D在0.27~5.00之间变化(具体为0.27、0.33、0.67、1.00、1.60、3.08、5.00)，用于分析不同算法对颗粒粒径范围的适用性；2) CFD网格的分布形式分为均匀网格与非均匀网格(如图 3c和3d所示)，用于分析不同算法对网格形态的适用性；3) 在计算过程中，记录部分算例的计算耗时，用于检验不同算法的效率。

关于SKM的2个模型参数：决定颗粒在空间影响范围的搜索半径R和Gaussian分布函数的参数b，研究成果并不多。虽然已知搜索半径R为参数b的函数[14]，但是b的确定目前尚未有明确结论，如Sun和Xiao[14]提出R=3b和b=6D，Sun等[15]提出R=b和b=0.45等。本文通过颗粒立方排布算例，确定搜索半径R与参数b的取值。

数值实验发现，计算结果与解析解之间的相对误差随R/b值呈单调递减趋势(如图 4a所示)，即R/b值越大，计算结果越准确。当R/b＞2.5时，与解析解之间的相对误差小于1%。其次，当R/b取值一定时，相对误差随b/D值的增大呈收敛趋势。若能取得相对误差收敛时的最小b/D值，则能够在保证满足精度要求的情况下，实现最高的计算效率。如图 4b所示，最小b/D值随R/b值呈线性关系，相关系数0.9884，通过拟合得到：

由此可以给出改进型SKM算法参数的确定方法：先根据精度要求确定R/b值(如要求相对误差小于1%时，则可取满足R/b > 2.5的任意值)，再根据R/b值通过式(11) 得到b/D值。根据颗粒粒径D值，即可确定模型参数。

如图 5所示，对于2种不同的颗粒排布形式(颗粒立方排列算例结果如图 5a—5c所示，颗粒紧密立方排列算例结果如图 5d—5f所示)，PCM与DPVM均呈现较差的适用性，当且仅当D≪L时，计算结果接近理论解；当D＞L时，PCM与DPVM均无法计算出准确结果。相较于PCM，DPVM计算结果在整体趋势上更接近理论解，说明体积分割方法可以一定程度上减小计算误差，避免颗粒体积全部归属某一个CFD网格所造成该处体积分数出现非物理性负值的情况，但当D＞L时，DPVM的误差仍然达到50%~100%，说明其在此种情况下仍然不适用。由此可见，如果将存在如此大误差的计算结果代入CFD-DEM耦合计算中，必然会产生较大误差并引起严重的数值不稳定性。

与PCM、DPVM相比，改进型SKM(使用节2.2中方法确定的模型参数)给出了准确的计算结果，如图 6a所示，当CFD网格为非均匀分布时，改进型SKM也同样给出了准确结果，证明其在不同应用场景下的适用性和准确性。当CFD网格尺寸与颗粒粒径的特征比值L/D在0.27~5.00之间变化时，不会对改进型SKM的计算结果产生明显影响，说明该算法可用于颗粒粒径与CFD网格尺寸相对变化范围较广的不同场景。L/D在0.27~5.00之间变化时计算结果的稳定性说明，该算法可以作为当D≪L时使用PCM或DPVM、与当D≫L时使用颗粒直接解析方法之间的过渡性算法，从而解决当颗粒粒径与CFD网格尺寸量级相当时传统算法不适用的问题。

同时，对于颗粒立方排列实验中颗粒粒径D与CFD网格尺寸L存在相对变化的4个算例，在模拟过程中同步记录计算耗时，如图 6b所示，PCM与DPVM计算效率基本相同(两者曲线重叠)，改进型SKM计算耗时虽然略高于PCM与DPVM，但与传统SKM(参数选用Sun和Xiao[14]所用参数)相比，特别是在颗粒粒径相对较大时(如颗粒粒径与CFD网格尺寸相同时)，改进型SKM大大降低了计算耗时，能够显著提高计算效率。

为进一步验证算法的可靠性和参数的适用性，采用静水颗粒沉速实验，将数值计算结果与沉速经验公式进行对比分析。

Concha[16]提出采用如下公式计算静水中的颗粒沉速：

其中：u*为无量纲颗粒沉速，d*为无量纲颗粒粒径。由如下表达式给出：

其中，P和Q分别为

其中：μf为流体黏度，ρf为流体密度，Δρ为颗粒与流体密度差，g为重力加速度值。

如图 7所示，在长、宽、高分别为100 mm×100 mm×200 mm的计算域内进行单颗粒沉降数值实验。初始状态下，颗粒被放置于距离底面190 mm处开始自由沉降，计算过程中仅考虑拖曳力、浮力和重力作用。颗粒粒径分别选用1、2、5、10和15 mm 5种工况，CFD网格尺寸与颗粒粒径比值在0.3~4.5之间变化，颗粒密度ρs为1 120 kg/m3。流体环境选用Cate等[17]使用的4组流体工况，分别使用PCM、DPVM和改进型SKM进行CFD-DEM的耦合模拟计算。算例参数见表 1。

如图 8a所示，相比于Concha[16]的经验公式，PCM和DPVM均无法给出准确的计算结果，沉速计算值系统偏大，这主要是由于在PCM和DPVM计算中，颗粒中心处的体积分数值往往产生非物理性负值或极小值，导致在计算过程中该处的流体速度为0或极小，由于颗粒的拖曳力与颗粒和流体间速度差成正比，因此会导致计算得到的拖曳力系统偏大，从而得到系统性偏大的沉速计算结果。在颗粒粒径较小时，即L≫D时，2种算法的计算结果接近准确值；但当颗粒粒径较大时(如d*＞100时)，2种算法所产生的误差最大可能接近一个数量级。就2种算法本身而言，DPVM的计算结果相对优于PCM，但颗粒粒径相对于CFD网格尺寸同量级或更大时，其计算结果仍不可用。

相比之下，使用节2.2中方法所确定的参数模型，传统型和改进型SKM在所有算例中均能得到准确的计算结果(两者计算精度相同，仅列出改进型SKM结果)，特别是在颗粒粒径相对于CFD网格尺寸同量级或更大时，也能得到准确沉速值，证明了该算法及其参数的适用性。同时，如图 8b所示，改进型SKM在保证相同计算精度的情况下，大大提升了模型的计算效率，以算例1中L/D=4.5时为例，其计算耗时仅为传统型SKM计算耗时的1/5；当计算域内颗粒数量进一步增加时，改进型SKM的计算效率优势将更加明显。

综上所述，本文所提出的改进型SKM算法在充分保障计算效率的情况下，能够适用于CFD网格尺寸与颗粒粒径同量级时的场景，能够覆盖颗粒粒径较小时使用传统算法(如PCM、DPVM)，颗粒粒径较大时使用颗粒直接解析模拟方法[18]间的过渡区域，有助于CFD-DEM耦合方法在多相流数值模拟计算，特别是在大规模颗粒固液耦合模拟计算中的应用。

本文提出了CFD-DEM耦合计算中的一种改进型SKM算法。通过特征点剖分、颗粒预分类及关键参数确定等技术，解决了传统SKM等算法在实际应用中遇到的非物理性失真、计算效率低等问题，提升了流体体积分数场的计算效率。通过静态颗粒排列实验与静水沉速实验对所提出的算法及其参数进行检验，结果表明：改进型SKM算法使计算颗粒粒径不再受到CFD网格尺寸的限制，覆盖了颗粒粒径较小时使用传统算法(如PCM、DPVM)，颗粒粒径较大时使用颗粒直接解析模拟方法间的过渡区域，有助于CFD-DEM耦合方法在多相流数值模拟计算，特别是在大规模颗粒固液耦合模拟计算中的应用。