基于完全解析的浸没边界方法（Immersed Boundary Method），得出了关于球形圆柱体的新的阻力与升力之间的关联关系

《AIChE Journal?AIChE》：A new drag and lift correlation for spherocylinders from fully resolved Immersed Boundary Method

　　摘要

许多工业过程涉及非球形颗粒的处理，例如矿物开采和生物质转化。理解颗粒的流体动力学对于控制和优化这些过程至关重要。为了将现有的研究从球形颗粒扩展到椭球圆柱形颗粒，我们进行了大量的颗粒分辨计算流体动力学模拟，以获得静态组装体上的阻力和升力。在1500次模拟中，改变了长宽比（）、雷诺数（）和固体分数（），以进行研究。模拟数据被用来建立新的阻力相关性（球体的误差为1.3%）和升力相关性（球体的误差为9.7%）。此外，随机排列的椭球圆柱形颗粒在阻力上表现出较大的方差，这与局部流动现象有关，这种方差在新的相关性中得到了量化。

1 引言

固体悬浮液在许多（工业）应用中经常出现，包括血液悬浮液（例如血浆中的红细胞）、矿物开采、生物质燃烧和化学反应器。控制这些过程需要理解这些多相系统的动力学。在许多工业应用中，悬浮液中的颗粒是非球形的。在各种应用中，颗粒形状显著影响了颗粒的动力学行为。将这些颗粒近似为球体可能会导致与观察到的颗粒行为不符。因此，需要准确描述这些颗粒的动力学。这些动力学可以通过最重要的表面力来描述，即阻力和升力，对于非球形颗粒来说，这些力取决于颗粒相对于流动的方向。一种常见的方法是在阻力或升力中加入一个球形系数，而不考虑颗粒的方向。虽然加入球形系数可以近似阻力，但众所周知，非球形颗粒的阻力强烈依赖于颗粒在流动中的方向，因此应该考虑这一点。H?lzer和Sommerfeld开发了一个通用相关性，将阻力与颗粒的纵向和横截面积联系起来。除了这些实验得出的相关性外，还进行了颗粒分辨模拟，以获得特定形状颗粒的相关性，包括椭球体、立方体和椭球圆柱体。在蠕动流动条件下，阻力主要由粘性阻力控制，因为流线通常遵循颗粒的形状。因此，对于与流动方向一致的颗粒和表面光滑的颗粒（如椭球体），粘性阻力会增加，因为在这两种情况下，流动都覆盖了颗粒的更大面积。在更高的雷诺数下，压力分布控制阻力，因为形状阻力变得占主导地位。在这种情况下，与流动方向一致的颗粒的阻力会降低，而与流动方向垂直的、具有较大平坦面积的颗粒的阻力会增加。此外，非球形颗粒的尾流过渡到不稳定流动的起始雷诺数会更低，对于具有尖锐边缘的颗粒（长方体和圆柱体）和倾斜角度较大的颗粒，这种不稳定性会更加明显。这种不稳定的尾流增强了形状阻力，同时也增强了升力和扭矩的振荡。这些现象对颗粒几何形状和雷诺数非常敏感，使得难以用单一的相关性来描述它们。除了这些形状外，椭球圆柱体需要特别关注，因为它们在工业应用中很常见，例如生物质气化和催化过程，其中它们是由圆柱体形成的，其尖锐的顶端会迅速磨损成更圆润的形状。椭球圆柱体在蠕动流动条件下的阻力较大，而在高雷诺数下阻力较小。因此，需要一个单独的相关性。考虑到实际的（工业）应用，颗粒通常以密集排列的方式存在。上述相关性能够预测孤立非球形物体上的力；然而，需要对其进行修改以考虑邻近颗粒的影响。通常，这是通过将之前讨论的非球形颗粒阻力相关性与密集球体系统的相关性结合起来实现的。一些例子包括Beetstra等人、Hill等人、Tang等人、Tenneti等人以及Ergun参考文献中的经验相关性。这些相关性将已知球体密集堆积的阻力从斯托克斯流动区域（Re1）扩展到更高的雷诺数。这些修正也用于非球形颗粒，结合了形状校正的孤立颗粒阻力模型。另一种方法是由Felice和Wen、Yu的经验工作提出的，它校正了单个孤立颗粒的阻力，以考虑其他颗粒的存在。第一种方法在高固体分数下提供了更高的准确性，而后者在稀薄系统中更准确。在密集堆积中，尾流的发展受到邻近颗粒的干扰，改变了局部压力分布，从而增加了形状阻力。此外，下游颗粒由于上游颗粒的尾流而经历较低的阻力，降低了局部速度。这种屏蔽效应在低雷诺数和稀薄系统中最为明显。随着固体分数的增加，狭窄间隙中的润滑应力增加，导致更大的粘性阻力。在稀薄堆积中，由于邻居效应最小，阻力接近单个颗粒的阻力。在更密集的系统中，阻力的时间波动加剧，这是由邻居尾流相互作用引起的。这种从稳定尾流到不稳定尾流的转变发生在显著较低的雷诺数下。

关于密集堆积的非球形颗粒的模拟研究非常有限。对于椭球体，已经推导出几个相关性。在所有这些相关性中，孤立椭球体的阻力都经过了密集系统的调整。这些研究表明，密集非球形颗粒系统的阻力不仅取决于单个颗粒的局部排列，还取决于整个颗粒组装在域中的集体方向。这导致了显著的流动模式（例如，对齐颗粒的通道化或错位床层的更曲折路径），这显著改变了相对于球体的阻力。对于非球形颗粒，观察到的阻力波动被放大，这是由于角度依赖的形状阻力对压力场局部扰动的敏感性。此外，在密集的非球形系统中，升力的影响相对于孤立椭球体变得更加显著。对于椭球圆柱体，Sanjeevi和Padding使用格子-玻尔兹曼方法推导出了一个阻力相关性，适用于长宽比为4的密集堆积椭球圆柱体。这个相关性基于同一作者获得的单个椭球圆柱体的阻力相关性。为了考虑长宽比对密集椭球圆柱体系统阻力的影响，我们进行了浸没边界模拟，以推导出新的阻力 and 升力相关性，这些相关性取决于长宽比、雷诺数、固体体积分数和颗粒方向。与当前的研究水平相比，这些相关性适用于不同的长宽比。此外，与Sanjeevi和Padding的工作相比，这项工作为完全随机配置提供了更具统计相关性的数据。这是在广泛验证和确认模拟方法之后完成的。尽管最初是为固定床配置推导出的，但该相关性不仅限于固定床反应，还用于流化床反应器的CFD–DEM模拟。Mahajan等人、Bo和Jianzhong的研究中使用的流化床雷诺数为300-600，长宽比为4-5，固体分数覆盖了当前研究的雷诺数范围（0.1-500）、长宽比范围（1-6）和固体分数范围（0-0.5）。

2 方法

采用颗粒分辨计算流体动力学模拟来解决这个多相系统。流动被认为是牛顿的、不可压缩的和等温的。该系统可以用方程（1）和（2）中的连续性和纳维-斯托克斯方程来描述。

（1）

（2）

其中 和 是系统中的流体速度和压力， 和 是流体的粘度和密度。使用结构化且交错的笛卡尔网格进行空间离散化。对流项使用MINMOD通量限制器通过延迟校正进行离散化。粘性应力项使用Dijkhuizen等人的公式进行解耦，并使用中心差分方案半隐式求解。这些非线性方程使用Biconjugate梯度稳定方法（Bi-CGStab(2)）求解。Chorin的投影方法用于高效求解纳维-斯托克斯方程和连续性方程。三个解耦的纳维-斯托克斯方程使用Trilinos IfPack2包中的不完全LU预处理器求解，以获得初步的速度场。为了满足连续性方程（即具有无散度场），使用Trilinos MueLu包中的代数MultiGrid预处理器求解压力校正方程。最终更新的速度是通过根据校正方程中的更新压力来校正近似速度得到的。有关更多信息，请参考Dijkhuizen等人的工作。

在固体表面上施加无滑移边界条件，使用Deen等人描述的二阶隐式幽灵单元浸没边界方法（IBM）。这种IBM应用于离散化的动量方程。对于每个与固体单元相邻的流体节点，使用插值幽灵值代替固体内部的邻居值。这个幽灵值基于使用表面值和两个相邻流体单元的二阶一维拟合，或者当第二个节点也在固体中时的固体速度。椭球圆柱体在固体表面的交点被视为半球形顶端和圆柱形中部与流体节点之间的最近交点，从而在两者之间实现逐渐过渡。如果表面与第一个流体节点之间的距离接近零，则在这种二阶拟合中会出现奇异性，在这种情况下，使用一阶拟合来确定幽灵值。颗粒上的总力使用方程（3）计算，作为粘性阻力和形状阻力的总和：

（3）

粘性力（也称为摩擦阻力）可以根据表面上的速度梯度使用二阶拟合来计算。投影面积可以近似为单元面的面积。如果梯度是垂直于力方向的，则使用Baltussen等人描述的方法对投影面积进行校正。这可以通过图1来说明。形状阻力（或压力阻力）可以通过将两个外部流体节点的压力外推到表面来计算，以获得表面的压力梯度。这是表面积分公式，与体积积分公式相比提供了更高的准确性。有兴趣的读者可以参考Deen等人的原始工作。如果第二个流体节点位于物体内部，则表面压力取为第一个流体节点的压力。

为了在阻力公式的表述上保持一致性，下面澄清了本工作中使用的定义。对于通过颗粒组装体的流动，颗粒主要经历两种力，即由于压力梯度产生的浮力（其中 是颗粒的体积）和（摩擦）阻力。这些力的总和是颗粒上的总力，有些出版物中称之为阻力。假设系统中有 个颗粒，颗粒的体积与系统的体积（）通过固体分数（）相关。根据牛顿第三定律，力等于反作用力，并且通过方程（4）与压降相关联。

（4）总力（F）与（摩擦）阻力（D）之间的关系可以推导如下：

（5）在接下来的章节中，将报告阻力（D）的计算方法，该方法与Beetstra等人19年和Tang等人21年的工作类似。阻力将基于斯托克斯流动条件下相同体积的孤立球体的解析阻力进行归一化，如方程（6）所示。

（6）其中Dv是体积等效直径，v是系统的入口速度。由于颗粒是非球形的，颗粒的方向很重要。在图2中显示了颗粒的局部参考框架，这与Sanjeevi和Padding32的工作类似。作用在颗粒上的力可以分解为全局参考框架的三个单位向量U、V和W。阻力（D）的方向与U相同，并且平行于入口速度。升力（L）的方向与W相同，垂直于U，并且与颗粒的主轴在同一平面内。

3 验证与确认
首先对所描述的数值方法进行验证和确认。由于使用了与Huijgen等人41年的工作相同的数值框架，因此建议读者参考该工作，以验证斯托克斯流动条件下阻力和扭矩的实现，并根据Schiller和Naumann的相关性42年验证孤立球体周围的流动。在接下来的章节中，验证将扩展到包括密集的球体系统和孤立的椭球柱体。

3.1 简单立方体排列的球体
为了确保该方法能够准确捕捉密集堆积中的流动，考虑了通过结构化球体堆积的流动。Sangani和Acrivos43年报告了斯托克斯流动条件下简单立方体、体心立方体和面心立方体晶格的阻力解析解。Zick和Homsy44年进行的模拟也作为这些相关性的参考数据。在这项工作中，模拟了一个由N个球体组成的简单立方体晶格，分辨率为Δx和Δy。图3A显示了不同固体分数下中心球体的阻力。随着固体分数的增加，阻力也随之增加。图3B表明，随着固体分数的增加，流体通过颗粒之间的更狭窄空隙以更高的速度流动。颗粒表面的插值压力显示出从前到后的压力梯度，从而产生形状阻力。此外，在颗粒表面附近观察到较大的速度梯度，导致皮肤摩擦增加。随着颗粒分数的增加，这些效应被放大，导致颗粒上的阻力增大。与解析解的小偏差可以由模拟设置的差异解释。当前工作采用了一个产生压力梯度的入口边界条件，代表了最终模拟的设置，与文献中由体力驱动的完全周期性域相反。

3.2 随机排列的球体
验证进一步扩展到随机排列球体的阻力，这与我们最终的模拟设置相当。对于这种配置，文献中有几种相关性，其中球体上的阻力与不同的雷诺数和固体分数有关19-23, 25。表2中列出的模拟设置与用于椭球柱体静态排列的模拟设置相同。对于当前的模拟，使用了0.2的固体分数。颗粒位置是使用Sanjeevi和Padding32描述的暴力方法生成的。图4A比较了与文献中报告的几种相关性19-23, 25计算出的阻力。尽管可用相关性之间存在差异，但计算出的阻力仍在相关结果范围内。Beetstra等人19年的相关性提供了模拟与模型之间最接近的相似性。对于不同的固体分数，差异会增加。此外，这还证明了阻力计算和求解器对于随机排列的球体也是准确的。图4B可视化了颗粒周围的流场，其中流线表示平均流动路径。在某些区域，流动主要是层流的，尽管在局部区域观察到流动分离和小涡流。这些局部湍流效应也会导致阻力计算的时间变化。这是最终分析中需要考虑的一个重要因素。

3.3 孤立椭球柱体上的阻力
浸没边界法可以应用于任意表面；然而，之前的验证仅包括球体。为了将验证扩展到非球形物体的实现，模拟了孤立椭球柱体周围的流动。对于孤立椭球柱体，需要较大的计算域，类似于Huijgen等人41年中对孤立球体阻力的研究。表1列出了模拟设置。使用的域和分辨率取决于模拟设置，即在低雷诺数下厚边界层需要较大的计算域，但允许较低的网格分辨率；而在高雷诺数下边界层变薄，需要更高的分辨率。所选的分辨率与文献中常用的值一致。

表1. 模拟设置：孤立椭球柱体上的流动。参数
Re < 3
计算域（L,H,D）：25a,29a,25a；24a,17a,14a；18a,22a,9a
分辨率：10；20
颗粒位置：0.3L, 0.5H, 0.5D
Courant数：0.1
长宽比（AR）：4
粘度（μ）：1.48
密度（ρ）：1.293

由于相对于流动的方向对于非球形颗粒来说是一个重要因素6, 9，因此通过改变方向角（θ）来进行模拟。如图5A所示，随着角度的增加，颗粒上的阻力也随之增加。模拟结果显示，不同方向角下的阻力与Sanjeevi等人9年和Feng和Michaelides13年的相关性吻合良好。这一结果证明了IBM适合模拟该系统。这种验证是在整个相关雷诺数范围内进行的，如图5B所示。在低雷诺数下，相对阻力趋于趋于1，对应于斯托克斯流动条件。随着雷诺数的增加，颗粒尾流中出现了流动分离，这与球形颗粒类似。与文献的轻微偏差可能是由于数值方法及其校准的差异（Sanjeevi等人9年采用了格子-玻尔兹曼方法，而Feng和Michaelides13年使用了连续驱动的浸没边界方法）。

3.3 孤立椭球柱体上的阻力
浸没边界法可以应用于任意表面；然而，之前的验证仅包括球体。为了将验证扩展到非球形物体的实现，模拟了孤立椭球柱体周围的流动。对于孤立椭球柱体，需要较大的计算域，类似于Huijgen等人41年中对孤立球体阻力的研究。表1列出了模拟设置。使用的域和分辨率取决于模拟设置，即在低雷诺数下厚边界层需要较大的计算域，但允许较低的网格分辨率；而在高雷诺数下边界层变薄，需要更高的分辨率。所选的分辨率与文献中常用的值一致。

由于相对于流动的方向对于非球形颗粒来说是一个重要因素6, 9，因此通过改变方向角（θ）来进行模拟。如图5A所示，随着角度的增加，颗粒上的阻力也随之增加。模拟结果显示，不同方向角下的阻力与Sanjeevi等人9年和Feng和Michaelides13年的相关性吻合良好。这一结果证明了IBM适合模拟该系统。这种验证是在整个相关雷诺数范围内进行的，如图5B所示。在低雷诺数下，相对阻力趋于趋于1，对应于斯托克斯流动条件。随着雷诺数的增加，颗粒尾流中出现了流动分离，这与球形颗粒类似。与文献的轻微偏差可能是由于数值方法及其校准的差异（Sanjeevi等人9年采用了格子-玻尔兹曼方法，而Feng和Michaelides13年使用了连续驱动的浸没边界方法）。

4 模拟设置
为了获得代表所有可能的椭球柱体配置的独立数据点，生成了完全随机的配置。生成过程遵循Sanjeevi和Padding32的方法。对于椭球柱体的模拟，使用了PACKMOL软件包。PACKMOL是为基于粘接球体表示的分子动力学模拟生成初始配置而开发的。每个椭球柱体由N个球体表示。表2总结了模拟设置。

在x方向的下部施加了入口条件，而在x方向的上部施加了1巴的压力。通过这种方式，可以精确地实现所需的雷诺数。在y和z方向施加周期性边界条件。为了防止入口和出口边界的影响，因为颗粒会穿透边界，所以应用了一个长度为L的入口和出口区域（见图2）。图7提供了设置的更详细表示。为了消除边界对计算阻力的影响，只包括位于域中心一半的颗粒，即图7中的紫色区域，这些颗粒被纳入阻力计算中。

图7显示了模拟设置的示意图。红色区域表示颗粒堆积，紫色区域表示在阻力计算中考虑的颗粒。为了确保结果与网格无关，使用50个颗粒进行了网格分辨率研究。颗粒受到雷诺数ρ的影响，这些雷诺数涵盖了本工作中使用的条件。图8A显示了相对于最高分辨率，阻力误差的减少。在ρ=1时，误差已经低于2%，并且在ρ=2时进一步减少到1%。考虑到这一点相对于计算成本的增加，选择了一个网格分辨率用于本研究的剩余部分。阻力被认为是进行分析的最适用指标，而本研究的主要发现是推导出的阻力相关性。

图8使用50个颗粒进行了网格收敛性研究（A），并使用ρ和L表示的域大小研究了粒子数量（B）（A）网格收敛性研究，（B）粒子数量测试。最后，考虑了域中的粒子数量。Sanjeevi和Padding32使用了200个颗粒，而用于球体堆积的数值研究的粒子数量各不相同（16,20；20–160,22；54,19；10821）。基于使用5到20种独立配置的分析，我们对本工作中所有图表中的每个数据点进行了20次独立模拟（因为使用超过15种配置时对平均阻力没有影响）。图8B显示了在不同粒子数量下的结果。结果清楚地表明，当粒子数量超过75个时，不会产生显著变化。因此，在本研究的剩余部分中使用了99个粒子。通过调整域的大小同时保持粒子大小不变来达到所需的固体分数。图9显示了本研究中调查的参数空间。选择参数空间中的间隙是为了确保最佳利用计算资源，同时仍然可以获得构建相关性所需的数据。这种方法确保了所得到的相关性在研究的雷诺数范围（0.1到500）、固体分数（0.1到0.5）和长宽比（1-6）内都是有效的。

5 结果
在本节中，我们将介绍对（密集）球体和椭球圆柱体堆积体进行的模拟。关于椭球圆柱体的结果将分为两部分：一部分是关于阻力的，另一部分是关于升力的。图10显示了在不同长宽比下通过堆积体的流动情况。该图表明，粒子会干扰流动，且每个粒子周围的局部流动不同，导致测量到的阻力和升力存在偏差。一个重要的观察结果是，由于施加的压力梯度，整体流动是从入口边界到出口边界的，尽管可能存在局部偏差。这一点已经通过压力场得到了证实。因此，根据方程（5）给出的浮力修正应该仅应用于x方向上的力，即阻力。

5.1 球体的阻力
在3.2节中观察到，尽管模拟结果落在现有相关性的范围内，但没有一个能够完全再现观察到的结果。因此，基于本工作的模拟结果推导出了一个新的球体阵列相关性。相关性的形状与Tang等人的类似。方程（8）显示了相关性，相关系数列在表3中。

图11A显示了模拟结果、推导出的相关性以及参考文献Feng和Michaelides中给出的长宽比为1的孤立椭球圆柱体的相关性之间的比较，后者应该得到Schiller和Naumann的相关性。模拟结果的再现非常准确，如图11B中的奇偶性图所示。

5.2 椭球圆柱体的阻力
一般来说，非球形粒子的阻力可以用正弦关系来描述。在这种公式中，阻力取决于在和处的阻力，以及角度，如方程（9）所示。
需要注意的是，指数的值存在争议，或者Feng和Michaelides提出了更高阶的多项式拟合。Sanjeevi和Padding使用了对于椭球圆柱体阵列，这也很好地拟合了当前的数据，因此在本研究中也使用了这个值。图12显示了得到的阻力结果，不同的颜色表示不同的粒子配置。数据点的分布证实了独立模拟的数量足以得出总体趋势。这些结果也证实了正弦关系的有效性。在密集区域（），阻力的偏差比稀疏区域（）更大，在稀疏区域，阻力趋近于孤立物体的渐近值。这是由于粒子之间的相互作用和相对位置导致的，从而产生了更明显的速度波动，进而导致阻力波动更大。这些增加的波动也降低了阻力比。Sanjeevi和Padding的工作中也观察到了这些效应。由于局部循环模式和场中的压力分布，个别粒子甚至可能出现负阻力。在图12中已经可以看到稀疏系统的一些局部循环模式。

图13显示了新相关性在描述当前模拟数据集方面相比参考文献Sanjeevi和Padding的相关性有显著改进。

5.3 方差
图12显示了阻力的显著方差，Sanjeevi和Padding、Li等人和Cao等人也观察到了这一点。相邻粒子的存在引入了非均匀的流动模式，即由于不同的局部粒子配置，观察到了不同的局部流动模式，从而导致阻力方差。为了在Euler-Lagrange模型中再现这些效应，应该包括这一点。有多种方法可以包含这些波动。可以从使用已知平均阻力和标准差在高斯分布上生成随机值开始。在Lattanzi等人的工作中，阻力是使用Ornstein–Uhlenbeck过程来建模的。时间波动来自动力学理论和颗粒温度。在自由移动的粒子集合中，这可以进一步研究。由于当前设置中只能研究空间波动，因此标准差将从观察到的阻力波动中得出。图16显示了空间波动的高斯分布。

图14显示了不同雷诺数和固体分数下每个粒子（点）的阻力与粒子角度之间的关系，并按角度分组（带误差棒的框符号）与文献（红线）和推导出的拟合（绿线）进行了比较。此外，图12还显示，对于低雷诺数（），阻力比趋于减小。在斯托克斯流动条件下，形状阻力减小，阻力主要由摩擦力主导，这较少依赖于粒子的形状。这也得到了Sanjeevi和Padding的证实，并在他们的相关性中显示出来。与Sanjeevi和Padding的相关性之间的差异可能源于使用的方法（格子-Boltzmann vs. IBM）和模拟设置的差异，例如边界条件和单向粒子堆积的添加。基于本工作中进行的验证，作者相信当前数据集准确地代表了系统，尽管没有独立的测量来验证绝对准确性。因此，即使对于长宽比为4的情况，也需要推导出新的相关性。方程（9）被拟合到数据中，以确定图9中每次模拟的和。这些值被收集起来，以推导出一个统一模型，其中阻力取决于倾斜角度、雷诺数、固体分数和长宽比。为了推导出这种相关性，需要考虑几种极限情况。对于稀疏系统（），应该使用参考文献Feng和Michaelides的相关性，因为它适用于所有长宽比。在球形系统的极限情况下（），阻力相关性预计将趋近于方程（8）中推导出的球体阵列的相关性。椭球圆柱体的推导阻力相关性显示在方程（10）中。这里，是方程（8），和是来自Feng和Michaelides的工作，A和B是拟合到数据中的修正因子。

图15显示，与参考文献Sanjeevi和Padding的相关性相比，新相关性在描述当前模拟数据集方面有显著改进。

图14显示了不同长宽比对阻力影响的效果。总体而言，随着长宽比的增加，阻力增加，这是由于随着长宽比的增加，流动中的扰动增强所致。对于更细长的粒子，粘性阻力增加，而在较低雷诺数下，流动附着在粒子的更大面积上。对于倾斜角度为的粒子，随着长宽比的增加，阻力减小，因为在以形状阻力为主导的情况下，即高雷诺数时，正面面积减小。对于垂直流动，这种效应是相反的，图14中也说明了这一点。Tang等人和Tenneti等人还观察到了反屏蔽效应，即由于狭窄的粒子间隙中的润滑应力增加，阻力随着固体分数的增加而增加。

图15显示了整个模拟集的奇偶性图，与Sanjeevi和Padding的相关性（SP）以及Feng和Michaelides的孤立椭球圆柱体相关性（不同空隙分数）进行了比较：方程（8）、Beetstra等人（BHK）、Tang等人（TPKKH）、Felice（DF）。在较高的固体分数下，颗粒间的距离减小，导致更多的颗粒聚集在颗粒的附近，从而产生更多的局部排列，因此方差也增大。总之，归一化方差与某些参数有很强的相关性。图17C显示了长宽比与方差之间的弱相关性。与球体的方差相比，球柱体的方差系统性地更大，这一点在He和Tafti15的研究中也观察到了，其中标准差从15%增加到了32%。这是由于颗粒周围流场的对称性丧失，导致尾流不对称以及局部流动现象的波动增大。基于观察到的关系，根据图17中的趋势提出了以下相关性，这也表明数据与相关性之间有很好的一致性：

（此处应插入具体的相关性公式或表达式）

5.4 升力

除了阻力之外，升力在这些非球形系统中也很重要，特别是在高雷诺数下，升力的大小可以与阻力相当。Sanjeevi和Padding的研究中显示了不同长宽比下的升力与倾角之间的关系。图18展示了在Re=10和AR=0.2时不同长宽比的升力情况。与Sanjeevi和Padding的研究类似，升力的分布更广，因为升力对局部流动动力学更为敏感。在45°附近可以看到一个最佳值，这与孤立非球形颗粒的结果相似。此外，升力的大小随着长宽比的增加而增加。在Re趋于无穷大时，升力的大小为零，这对于非旋转球体来说是预期的。这种趋势可以通过非球形程度的增加导致流体垂直于流动方向的偏转增大，从而产生更大的流体-颗粒力来解释。这一点也得到了图10的证实。由于沿流动方向的压强较高，而流动浅区域的压强较低，这种力来源于形状阻力。尽管Sanjeevi等人9和Zastawny等人10观察到升力存在轻微的偏斜，但在图18的数据中无法清晰地观察到这一点，因此根据奥卡姆剃刀原理忽略了这一因素。为了推导出量化这些颗粒升力的相关性，需要考虑一些边界条件。对于孤立颗粒，升力应减少到孤立颗粒的升力，为此使用了Feng和Michaelides13的相关性，并结合了方程（7）的简化模型，因为该相关性适用于整个研究范围。此外，对于长宽比为1的静止颗粒，由于假设颗粒周围的流动是对称的，因此没有升力，这也符合Feng和Michaelides的相关性。基于这些要求，升力的相关性被拟合到方程（12）中。系数C是通过使用PySR软件包来最佳地表示数据的。

5.5 对齐

从结果来看，与Sanjeevi和Padding之前的工作有显著差异。之前工作的主要区别之一是考虑了对齐的阵列。因此，本文探讨了对齐的效果。图19比较了通过随机堆积和通过单向堆积（无论是顺流还是横流对齐）的流动情况。在单向堆积中，流动遵循优先路径，导致堆积中的明显通道效应。相比之下，随机堆积中的流动更加均匀。这种通道效应使得推导统一的阻力相关性变得复杂。此外，工业应用中（如流化床反应器）是否会发生对齐也是值得怀疑的。在Mahajan等人的实验工作中，只观察到小范围的对齐现象。因此，只有随机堆积被用于推导阻力相关性。

6 结论

使用二阶幽灵单元浸没边界法，数值研究了不同长宽比（AR）的球柱体的阻力和升力。该方法已经经过了彻底的验证和确认。为了确保相关性与配置的独立性，每个模拟结果包括99个颗粒和20次独立模拟。模拟显示，随着长宽比、雷诺数、固体分数和倾角的增加，阻力也会增加。为了推导球柱体的相关性，首先推导了球体阵列的相关性，以提供当长宽比接近1时的极限情况（方程8）。这种球体阵列的阻力相关性的平均绝对百分误差为1.5%。球柱体阻力相关的另一个极限是孤立球体的阻力。基于这些极限，推导出了新的球柱体阻力相关性。这种相关性是基于0°和90°倾角下的阻力构建的，取决于Re、AR和AR。结果通过正弦关系与倾角相关联，正如文献中所证明的有效方法。新推导的相关性（方程10）的平均绝对百分误差为9.7%，这比使用孤立球柱体相关性和空隙率校正的组合有了显著改进。此外，还推导了相对方差的相关性，该相关性仅依赖于固体分数和倾角。此外，还研究了升力。对于这些非球形颗粒，升力是显著的。升力在某个特定角度达到最大值。对于孤立球柱体观察到的轻微偏斜，在球柱体阵列中并未出现。随着固体分数、雷诺数和长宽比的增加，升力也随之增加。基于模拟，为球柱体阵列推导出了新的升力相关性。获得的阻力和升力相关性可以用于涉及非球形颗粒的大规模过程的模拟，例如流化床反应器的CFD-DEM模拟。这将比文献中的相关性更好地预测这些非球形颗粒的流体动力学特性。为了进一步提高阻力相关性的准确性，应该使用自由移动的球柱体进行模拟，以包括颗粒移动性对阻力的影响。这些模拟将更真实地反映流化床反应器的一部分。

作者贡献

A. H. Huijgen：概念化、数据整理、形式分析、研究方法、软件、监督、验证、可视化、写作——原始草稿准备。
Y. Sterkenburgh：概念化；研究方法；软件；写作——审阅和编辑。
Y. Tang：概念化、形式分析、资金获取、监督、写作——审阅和编辑。
N. G. Deen：形式分析；监督；项目管理；写作——审阅和编辑；资金获取。
M. W. Baltussen：资金获取；资源；项目管理；写作——审阅和编辑；形式分析；概念化；研究方法；软件；研究；验证；监督。

致谢

本出版物是Open Technology Programme项目“Understanding of wet agglomeration of non-spherical particles in fluidized beds”（项目编号19951）的一部分，该项目由荷兰研究委员会（NWO）资助。我们感谢SURF（www.surf.nl）在使用国家超级计算机Snellius方面的支持。作者感谢prof.dr.ir. Johan Padding就他们过去工作的解释和一般建议所提供的宝贵讨论，以及prof.dr.ir. Hans Kuipers在研究初期阶段提供的宝贵讨论。

资金信息

本工作得到了荷兰科学研究组织（NWO）OPT项目的支持，项目编号为“Understanding of wet agglomeration of non-spherical particles in fluidized beds/19951”。

数据可用性声明

原始数据和处理脚本已放在appendix graphs.zip（一个zip文件）中。获取图3A的脚本位于graphs.zip/dragOnSCLatticeSphere/dragOnSCLatticeSphere.py，而原始数据位于graphs.zip/dragOnSCLatticeSphere/rawdata。图3B的原始数据位于graphs.zip/dragOnSCLatticeSphere/raw_simulation_data_es_028。获取图4A的脚本位于graphs.zip/dragOnArrayOfParticles/plotSpheres.py，而数据位于graphs.zip/dragOnArrayOfParticles/processeddata/drag99_acsi_by_ar_phi_re_iter_res_025_025.json。图4B的原始数据位于graph.zip/dragOnArrayOfParticles/raw_simulation_data_AR_1_es_02_re_100。获取图5A和6A的脚本位于graph.zip/dragOnSingleSpheroCylinder/plotSpheroCylinderRe30.py，而获取图5B和6B的脚本位于graph.zip/dragOnSingleSpheroCylinder/plotSpheroCylinderAllRe.py。两者的原始数据都位于graph.zip/dragOnSingleSpheroCylinder/rawdata。获取图8A的脚本和原始数据位于graph.zip/resolutionStudy/PlotResolution.py和graph.zip/resolutionStudy/rawdata。获取图8B的脚本和原始数据位于graph.zip/plotNumberOfParticles.py和graph.zip/numberOfParticles/rawdata。获取图9的脚本位于graph.zip/parametricspace/plotParamSpace.py。获取图13-20的脚本位于graph.zip/dragOnArrayOfParticles/plotSpheroCylinders.py。数据位于graphs.zip/dragOnArrayOfParticles/processeddata/drag99_acsi_by_ar_phi_re_iter_res_025_025.json和graphs.zip/dragOnArrayOfParticles/processeddata/lift99_acsi_by_ar_phi_re_iter_res_025_025.json。图12的原始数据位于graph.zip/dragOnArrayOfParticles/raw_simulation_data_AR_2_es_01_re_10、graph.zip/dragOnArrayOfParticles/raw_simulation_data_AR_4_es_01_re_10和graph.zip/dragOnArrayOfParticles/raw_simulation_data_AR_6_es_01_re_10。图19的原始数据位于graph.zip/dragOnArrayOfParticles/raw_simulation_data_AR_4_es_02_re_10-aligned_0、graph.zip/dragOnArrayOfParticles/raw_simulation_data_AR_4_es_02_re_10-aligned_90和graph.zip/dragOnArrayOfParticles/raw_simulation_data_AR_4_es_02_re_10。获取图17的脚本位于graph.zip/dragOnArrayOfParticles/plotVariance.py。数据位于graphs.zip/dragOnArrayOfParticles/processeddata/drag99_acsi_by_ar_phi_re_iter_res_025_025.json。获取图16的脚本位于graph.zip/dragOnArrayOfParticles/plotPDF.py。数据位于graphs.zip/dragOnArrayOfParticles/processeddata/drag99_acsi_by_ar_phi_re_iter_res_025_025.json。使用PySR进行升力相关性的拟合在graphs.zip/dragOnArrayOfParticles/findCorrelationLift.py中完成，经过大量迭代后，将适当的相关性打印到graphs.zip/dragOnArrayOfParticles/processeddata/correlationlift.txt中。其中最佳的相关性被包含在本工作中。支持本研究发现的数据可以在4TU研究数据中公开获取，网址为https://doi.org/10.4121/0710ce66-e74a-4938-93d8-7ac129a13768。