具有两个对称平面的刚性颗粒在沉积过程中趋向于稳定、静止方向的实验与数值研究

《Journal of Fluid Mechanics》：Experimental and numerical study of rigid particles with two planes of symmetry approaching a stable, stationary orientation while sedimenting

　　摘要：本研究通过实验和数值方法探讨了在低雷诺数（远小于1）条件下，具有两个正交对称平面的刚性颗粒在高度粘性流体中受重力作用时的运动规律。Joshi和Govindarajan（2025年《物理评论快报》134(1), 014002）从理论上指出，对于这类颗粒，其运动特性会随着两个转动-平移移动性系数的乘积符号的不同而发生变化，这些系数是相对于颗粒质心在一个对称参考系中评估的。然而，仅通过观察颗粒的形状，尚不能立即判断出该乘积是负值、正值还是零。在本文中，我们展示了如何通过实验（利用特定的初始方向）以及基于斯托克斯方程的数值方法来估计这些系数及其乘积的符号。特别值得关注的是“沉降颗粒”——这类颗粒会重新定位并趋向于一个稳定的静止状态，我们的研究重点也放在这一类形状上。实验结果表明，圆锥体、新月形、箭头形状和开口平坦的环状颗粒都属于沉降颗粒，我们通过实验确定了它们的转动-平移移动性系数。随后，我们将每个实验形状重构为许多相互接触的珠子的集合体，并使用精确的Hydromultipole代码来计算该集合体的移动性系数。数值结果与实验结果十分接近，从而确认这些颗粒确实属于沉降颗粒，并能够估算出它们的重新定位时间尺度。我们的发现不仅适用于水基溶液中的非布朗微观物体（基于实验结果和斯托克斯方程的理论推导），也具有普遍适用性。沉降颗粒在相对较短的时间内重新定位到稳定状态的现象，可能在环境科学、生物学、医学或工业领域有应用价值。

1. 引言：在多种系统中，可以观察到复杂形状的小颗粒在粘性流体中的重力沉降现象，这一现象已在文献中得到了广泛研究。例如空气中的雪花（Jayaweera和Mason的研究）、水生环境中的自然絮状沉积物或泥浆（Strom和Keyvani的研究；Spencer等人以及Gu等人）、海洋中的海洋雪（Alldredge和Gotschalk的研究；Diercks和Asper的研究）、密度大于水的微生物的游动（Rizvi等人的研究），以及水和空气中的微塑料纤维的传输（Tatsii等人的研究）。Tatsii等人指出，微塑料颗粒的沉降速度越小，它们被风传输的距离就越远，有时可达数千公里。另一方面，颗粒的沉降速度显著受到其形状的影响。各向异性颗粒在重力作用下会经历复杂的流体力学相互作用，这些相互作用已在许多简化模型系统中得到分析。对于较大的颗粒，惯性效应显得尤为重要（Jayaweera和Mason的研究；Candelier和Mehlig的研究；Ravichandran和Wettlaufer的研究；Maches等人的研究；Tatsii等人的研究）。对于微米级颗粒，采用无惯性描述即可再现其运动规律。实际上，根据具体情境，雷诺数可能小于一（Ravichandran和Wettlaufer的研究），甚至可能小于0.01（Jayaweera和Mason的研究）。对于雷诺数远小于1的粘性流体中不同复杂形状的单个刚性颗粒的沉降现象，已经进行了长期研究（参见Brenner的早期论文）。实验和数值方法（例如Lasso和Weidman、Cichocki和Hinsen、Weber等人的研究）研究了空心圆柱体、球体集合体及刚性结等的沉降行为，目的是确定颗粒的沉降速度与相同质量和体积的球体的斯托克斯速度相比有何差异。一般来说，细长颗粒或盘状颗粒的沉降速度取决于其姿态，并伴随着侧向漂移（Lamb的研究；Taylor的研究）。此外，根据形状的不同，颗粒在沉降过程中可能会发生旋转。对于在不可压缩粘性流体中受重力作用并产生满足静止斯托克斯方程的流体流动的颗粒，其角速度和质心平移速度分别由相对于质心评估的$3\times 3$转动-平移和平移-平移移动性矩阵决定（Happel和Brenner、Durlofsky和Bossis、Felderhof、Kim和Karrila、Graham的研究）。已知对于手性颗粒，存在转动-平移耦合现象，这意味着相对于质心评估的转动-平移移动性矩阵不为零。手性物体的耦合转动和平移动力学已得到广泛研究，包括刚性结（Gonzalez和Maddocks的研究）、类似螺旋桨的手性形状（Doi和Makino的研究）、不对称的刚性纤维（Tozzi等人的研究）、刚性螺旋结构（Palusa等人的研究）以及螺旋形刚性带状物（Huseby等人的研究）。然而，即使对于非手性颗粒，转动-平移耦合也可能是非零的（参见Gonzalez等人的研究；Doi和Makino的研究；Krapf、Witten和Keim的研究）。例如，Kim和Karrila研究了非手性螺旋桨（Kim和Karrila的研究；Makino和Doi的研究；Doi和Makino的研究；Krapf、Witten和Keim的研究）；Ekiel-Je?ewska和Wajnryb通过数值和理论方法研究了弯曲的刚性纤维（Candelier等人的研究）；Miara等人通过实验、数值和理论方法研究了弯曲的刚性圆盘（Miara等人的研究；Vaquero-Stainer等人的研究）。Joshi和Govindarajan最近研究了具有两个正交对称平面的物体的沉降动力学（Joshi和Govindarajan的研究），指出每种形状属于三类不同的运动类型：沉降颗粒（趋向于稳定的静止状态）；漂移颗粒（趋向于倾斜姿态并伴有侧向漂移）；以及振动颗粒（在沉降过程中进行周期性振荡）。本文重点讨论沉降颗粒的动力学特性，并提供了沉降颗粒运动的实验和数值实例，表明它们达到稳定状态的典型时间足够短，具有实际应用潜力（Spencer等人的研究；Tatsii等人的研究）。选择这些实验用的刚性形状是基于Ekiel-Je?ewska和Wajnryb之前关于由三个相同球体组成的等腰三角形（但非等边三角形）的顶部颗粒沉降的研究，其中中间球体与其他球体接触。研究发现，这种颗粒会旋转至中间球体位于底部的稳定姿态。我们选择了几种容易获得的类似铃铛的形状进行实验：圆锥体、箭头形状、新月形和开口不同的开放环形颗粒。实际上，所有这些形状都属于沉降颗粒。然而，最初不确定这两种稳定姿态中哪种是稳定的：箭头形状、新月形和开放环形颗粒会“宽边朝上”重新定位，而圆锥体则“宽边朝下”是稳定姿态。人们可能会直观认为圆锥体的“宽边朝下”是稳定姿态，因为这一侧更重（如Candelier和Mehlig、Nissanka和Burton的研究所示），但这一结论不适用于Ekiel-Je?ewska和Wajnryb研究中的铃铛形状以及我们实验中使用的其他形状。Jayaweera和Mason还指出，在较大的雷诺数（Re ≥ 0.5）下，锥体的尖角小于π/4时，锥体会将尖端向上调整；而在尖角大于π/4时，尖端会向下调整；对于雷诺数远小于1的情况，目前尚无相应的研究结果。Joshi和Govindarajan从理论上指出，对于雷诺数远小于1的颗粒，当两个正交对称平面的转动-平移移动性系数的乘积符号不同时，其运动特性会有所不同（见图1(c)）。本文展示了如何通过实验选择两种不同的初始方向来确定这些系数，并基于斯托克斯方程进行数值计算。我们提出了一种使用珠子模型来完成这一任务的方法，并将其应用于本研究中实验观测的颗粒（沉降颗粒，其转动-平移移动性系数满足μ42 ≤ μ51 ≤ 0）。数值计算得到的转动-平移移动性系数值与实验数据相当吻合，证实这些颗粒属于沉降颗粒，并能预测它们的重新定位时间尺度。本文提出的实验和理论方法也可用于研究振动颗粒（μ42 ≥ μ51）和漂移颗粒（μ42 = μ51 = 0）。本文的结构如下：第2节提供了理论背景并解释了沉降颗粒的转动-平移移动性矩阵的特征。第3节描述了实验装置、所用颗粒及方法。第4节展示了实验结果：初始姿态和实验数据统计；颗粒的平移和旋转典型示例（通过时间序列快照显示）；随时间变化的垂直速度；随时间变化的倾斜角度以及移动性系数的平均值。第5节将实验形状近似为相同球体的集合体，并通过数值方法确定它们的移动性系数。第6节对比了理论和实验结果，第7节对开口可控的颗粒进行了参数化数值研究。最后，第8节总结了研究结果。

2. 理论背景：我们考虑一个在粘性流体中移动的刚性颗粒，其所受的减重力$\boldsymbol{F}$等于重力减去浮力。雷诺数Re远小于1，流体速度$\boldsymbol{v}(\boldsymbol{r})$和压力$p(\boldsymbol{r})$由静止斯托克斯方程描述（2.1）。颗粒的运动是基于（2.1）推导出的，并假设颗粒表面的流体满足无滑移边界条件（Kim和Karrila的研究；Graham的研究）。颗粒的平移速度$\boldsymbol{U}$和角速度$\boldsymbol{\varOmega}$是作用在颗粒上的减重力$\boldsymbol{F}$的笛卡尔分量的线性函数。在实验室参考系$XYZ$中，减重力$\boldsymbol{F}$与$Z$轴方向相反，即$\boldsymbol{F}=-F\boldsymbol{\hat {Z}}$。我们根据质心来计算$\boldsymbol{U}$。这样做的好处是，相对于质心来说，力矩为零。因此，我们有：
$$
(\begin{array}{ll}\boldsymbol{\mu }^{\textit{tt}} \,\boldsymbol{\mu }^{tr} & \boldsymbol{\mu }^{\textit{rt}} \, \boldsymbol{\mu }^{\textit{rr}}\end{array}) \boldsymbol{\cdot } \boldsymbol{F} = 0,
$$
或者，也可以表示为：
$$
\boldsymbol{U} = \boldsymbol{\mu }^{\textit{tt}} \boldsymbol{\cdot }\boldsymbol{F},
$$
以及
$$
\boldsymbol{\varOmega } = \boldsymbol{\mu }^{\textit{rt}} \boldsymbol{\cdot }\boldsymbol{F}.
$$
平移-平移迁移率$\boldsymbol{\mu}^{tt}$和旋转-平移迁移率$\boldsymbol{\mu}^{rt}$是$3 \times 3$的笛卡尔矩阵，它们取决于粒子的朝向。在本文中，如同Joshi与Govindarajan（参考文献Joshi and Govindarajan2025）以及Ekiel-Je?ewska与Wajnryb（参考文献Ekiel-Je?ewska and Wajnryb2009a）在他们的理论研究中所做的那样，我们关注的是那些具有两个垂直对称平面的粒子（但相对于第三个垂直平面的反射并不对称）。我们选择的参考系中，$x$和$y$轴与对称平面垂直，而$x$轴沿着粒子的最大尺寸$L$方向。对于我们大多数粒子来说，它们在$y$方向上的尺寸小于$x$方向；对于锥形粒子来说，$x$和$y$方向上的尺寸是相同的。在这个特殊的参考系中，无量纲迁移率矩阵的形式为：
$$
\begin{array}{ll}
\pi \eta L \boldsymbol{\mu }^{\textit{tt}} &= \left (\begin{matrix}\mu _{11} & 0 & 0\\ 0 & \mu _{22} & 0\\ 0 & 0 & \mu _{33} \end{matrix} \right ),\quad \pi \eta L^2 \boldsymbol{\mu }^{\textit{rt}} = \left (\begin{matrix} 0 & \mu _{42} & 0\\ \mu _{51} & 0 & 0\\ 0 & 0 & 0\end{matrix} \right ).
$$
在下文中，我们将使用$L$作为长度单位，无量纲时间$t$定义为：
$$
t=\frac {F}{\pi \eta L^2}.
$$
其中$F=|\boldsymbol{F}|$，$\tilde {t}$表示有量纲时间。图1展示了两种具有$x\rightarrow -x$和$y\rightarrow -y$对称性的粒子的特殊运动：(a)情况A，粒子绕$y=Y$轴旋转；(b)情况B，粒子绕$x=X$轴旋转。这里的$xyz$是随粒子移动的坐标系，如图(c)所示，$XYZ$是实验室参考系，$\theta$是$z$轴和$Z$轴之间的角度。我们考虑的是那些会旋转并趋向于静止姿态的粒子，不论初始条件如何。这类形状在Joshi与Govindarajan（参考文献Joshi and Govindarajan2025）的理论研究中被称为“settlers”。对于这些settlers，旋转-平移系数$\mu _{42}$和$\mu _{51}$不等于零，并且符号相反，即：
$$
\mu _{42} \,\mu _{51} < 0.
$$
问题是如何通过实验来确定$\mu _{42}$和$\mu _{51}$。Ekiel-Je?ewska与Wajnryb（参考文献Ekiel-Je?ewska and Wajnryb2009a）指出，有两种特殊情况下的动力学非常简单。这些特殊情况A和B在图1中用实验室参考系$XYZ$和随粒子移动的参考系$xyz$进行了说明。对于情况A（如图1a所示），如果我们选择初始条件使得$y$轴和$Y$轴平行，那么由于对称性，$y$轴和$Y$轴将始终保持平行，质心的运动将发生在$XZ$平面上。粒子的姿态由$Z$轴和$z$轴之间的角度$\theta (t)$确定，该角度随时间$t$的变化遵循以下方程：
$$
\frac {\mathrm{d}\theta }{\mathrm{d}t}= \mu _{51} \sin \theta.
$$
对于情况B（如图1b所示），如果我们选择初始条件使得$x$轴和$X$轴平行，那么由于对称性，$x$轴和$X$轴将始终保持平行，质心的运动将发生在$YZ$平面上。粒子的姿态同样由$Z$轴和$z$轴之间的角度$\theta (t)$确定，该角度随时间$t$的变化遵循以下方程：
$$
\frac {\mathrm{d}\theta }{\mathrm{d}t}= -\mu _{42} \sin \theta.
$$
这里，方程(2.7)和(2.8)与Ekiel-Je?ewska与Wajnryb（参考文献Ekiel-Je?ewska and Wajnryb2009a）引入的方程(14)相同，并且在那里作为方程(15)进行了求解。因此，情况A和B的解分别为：
$$
\tan (\theta /2) = A \exp (\mu _{51} t),
$$
和
$$
\tan (\theta /2) = A \exp (-\mu _{42} t).
$$
Ekiel-Je?ewska与Wajnryb（参考文献Ekiel-Je?ewska and Wajnryb2009a）在他们的方程(10)–(12)中给出了settler非平面运动的一般表达式，Joshi与Govindarajan（参考文献Joshi and Govindarajan2025）也有类似的表达式。为了简化分析，在实验中，我们旨在提供粒子如情况A或情况B那样的初始姿态，测量角度$\theta$随时间的变化，然后与方程(2.9)或(2.10)中的相应解进行比较，具体方法将在第4节中描述。在实验中，我们将直接分析有量纲量，特别是定义为以下形式的 time系数：
$$
\tilde {\tau }_{51} \equiv \frac {\pi \eta L^2}{F} \frac {1}{\mu _{51}},\,\, \mbox{ 和 }\,\,\,\, \tilde {\tau }_{42} \equiv \frac {\pi \eta L^2}{F} \frac {1}{\mu _{42}},
$$
基于以下条件：$\mu _{51}t=\tilde {t}/\tilde {\tau }_{51}$和$\mu _{42}t=\tilde {t}/\tilde {\tau }_{42}$。在实验中，我们将展示粒子是如何重新定向并达到静止状态的（见第4节）。通过测量最终的垂直静止质心速度$U_{\!f}$和减重力$F$，我们将确定粒子在静止状态下平移的相关无量纲迁移率系数：
$$
\mu _{33} = \pi \eta L \frac {U_{\!f}}{F}.
$$
我们实验中使用的每种粒子形状随后都被近似为一系列相互接触的珠子，迁移率系数(2.4)是通过Cichocki、Ekiel-Je?ewska与Wajnryb（参考文献Cichocki, Ekiel-Je?ewska and Wajnryb1999）和Ekiel-Je?ewska与Wajnryb（参考文献Ekiel-Je?ewska and Wajnryb2009b）提出的多极方法进行评估的，同时考虑了润滑效应（结果见第5节）。

3. 实验系统、材料和方法
3.1. 实验装置
实验装置几乎与Shashank、Melikhov和Ekiel-Je?ewska（参考文献Shashank, Melikhov and Ekiel-Je?ewska2023）使用的装置相同，但没有使用镜子。实验是在一个内尺寸为200毫米×200毫米×500毫米（分别代表宽度、深度和高度）的玻璃容器中进行的。我们使用了高粘度硅油（由Silikony Polskie制造），其运动粘度为$\nu = 5 \times 10^{-3} \,\text{m}^2 \;\text{s}^{-1}$，密度为$\rho = 970 \,\text{kg m}^{-3}$，温度为25摄氏度。将具有两个正交对称平面的固体物体（以下简称为粒子）放入这个容器中，研究它们是否会重新定向并趋向于静止状态。粒子的运动通过两个分别垂直于容器前后面的摄像机记录下来，摄像机的视线指向容器中心。两个摄像机都直接对准玻璃容器，且两个光源都保持与容器壁平行，并配有散射器，以确保整个容器区域照明均匀，如图2所示。光源、容器和摄像机的布置确保了不透明的沉积粒子在明亮背景下清晰可见。我们使用了两台相同的全帧数码单反相机（Canon 5D Mark IV），分辨率为3000万像素，配备100毫米的定焦镜头。这些摄像机连接到由Esper Ltd.（英国制造）制造的多相机操作触发箱。摄像机通过Esper触发箱的操作系统同时触发，该系统由计算机控制。触发箱同时触发两个摄像机的时间延迟为20毫秒。图像的捕捉速率为一帧每秒。两个摄像机都设置了最高的$f$数$f=32$（对应于最小光圈），曝光时间为1/25秒。物体在给定曝光时间内移动的距离范围为1-7像素。物体在曝光时间内的运动表现为图像模糊。附录B分析了模糊对角度$\theta$测量值的影响，得出结论：在当前研究中可以忽略运动模糊的影响。为了实现足够的亮度并保持适度的噪声水平，两台相机的ISO都设置为160。

3.2. 刚性粒子
本研究涉及的刚性物体（来自波兰的Koniarscy S. C.）具有两个垂直对称平面，但在第三个垂直平面上的反射并不对称：包括玻璃制成的新月形和圆锥形物体、由赤铁矿制成的箭头形状物体以及钢制开口环形物体。我们手动打开了这些环形物体，并调整了它们两端之间的间隙大小。图3展示了这些物体在随粒子移动的$xyz$参考系中的真实图像，其中$z$轴位于对称平面$xz$和$yz$的交点。图3中显示的所有尺寸都是使用高精度数字游标卡尺测量的，测量精度为每次0.01毫米。每个尺寸至少在五个物体上进行了测量，并重复了两到三次，以考虑单个物体尺寸的差异并减少测量误差。物体的尺寸及其标准偏差如下：
(i) 圆锥形物体：直径为$L= 4.39 \pm 0.02$毫米，高度为$H=10.45 \pm 0.05$毫米。一个直径为$0.96\pm 0.05$毫米的圆柱孔穿过圆锥体，与该物体底面平行，并在距离底面$1.70 \pm 0.01$毫米的高度处与中心轴相交。
(ii) 新月形物体：长度为$L=10.18 \pm 0.02$毫米，总高度为$H=4.46 \pm 0.09$毫米，中间高度为$S=3.91 \pm 0.14$毫米，宽度有所不同，中间宽度为$W_M=2.13 \pm 0.02$毫米。两个直径为$0.85\pm 0.02$毫米的孔穿过新月形物体的厚度。这些孔对称地分布在物体底部上方$1.47 \pm 0.06$毫米的位置，中心之间的距离为$2.99 \pm 0.04$毫米。
(iii) 箭头形物体：长度为$L=6.28 \pm 0.08$毫米，总高度为$H=4.93 \pm 0.08$毫米，中间高度为$S=3.81 \pm 0.02$毫米，宽度为$W=2.33 \pm 0.10$毫米。一个与$z$轴同轴的孔直径为$1.06\pm 0.02$毫米，穿过箭头并从顶点伸出。
(iv) 开口环形物体：我们测量了这些环的长度$L$和高度$H$，分别为$L=5.45 \pm 0.01$毫米、$L=5.92 \pm 0.04$毫米、$L=6.48 \pm 0.08$毫米、$L=6.81 \pm 0.02$毫米和$H=4.95 \pm 0.01$毫米、$H=4.90 \pm 0.02$毫米、$H=4.88 \pm 0.03$毫米，开口直径分别为$1.02 \pm 0.05$毫米、$2.02 \pm 0.06$毫米、$3.02 \pm 0.05$毫米和$4.02 \pm 0.05$毫米。所有环的宽度保持恒定为$W=0.596 \pm 0.005$毫米。实验研究中使用的物体从不同角度展示，其尺寸以毫米为单位标明。刻度尺（白线）代表5毫米。物体的物理性质——质量、体积以及减去浮力后的净重$F$——在表1中列出。物体的质量是通过高性能分析天平测量的，每次测量的精度达到0.1毫克。然而，同一类型的单个物体质量存在差异，因此对至少10个相同类型的物体进行了质量测量，并将平均值及其标准差列在表1中。圆锥体、箭头和新月形的体积是通过排水法实验测量的，这是一种用于测量不规则形状固体体积的成熟技术。对于箭头和新月形，我们使用了50个颗粒，并假设排开的水体积等于所有颗粒的总体积。由于无法获得同一批次的材料，我们只使用了12个圆锥形颗粒进行体积测量，从而确定了单个物体的平均体积。对于更小的环形物体，通过测量多个封闭环的尺寸并使用环形体积的几何公式来更精确地计算体积。然后使用平均质量和体积来计算净重$F$（重力减去浮力）。平均值及其不确定性列在表1中。最后，请注意，我们使用的是同一类型的不同物体进行实验，每个物体只使用了一次，并且我们使用平均值来近似它们的质量、体积和净重。

表1. 实验研究中使用的各种物体的物理性质。

3.3 方法论
实验是通过将物体从水箱顶部手动落下，大致位于液体表面的中心位置来进行的。当物体即将进入相机的可视区域时，相机会被触发；当物体离开相机视野时，相机停止记录。在进行每一组实验之前，我们将尺子放在面对两个相机的容器外侧墙壁的中央垂直线上，并调整相机位置，以确保两个相机拍摄的图像具有相同的高度$h$。物体的运动发生在容器中央垂直线附近。由于相机的视角有限（硅油的折射率为1.4），实际图像的高度$h'$大于$h$。为了确定这个修正因子$h'/h=1.077$，我们将尺子放在水箱的中央线上测量实际图像高度$h'$，并将另一个尺子放在容器外侧前壁的中央垂直线上测量图像高度$h$。然后使用这个修正因子将像素位置转换为毫米。大约一个像素对应0.052毫米，但这个校准因子是针对每一组实验专门确定的。我们使用了预先湿润和未预先湿润的颗粒。在某些未预先湿润的圆锥体、箭头和新月形的实验中，我们观察到一个小气泡从物体上的孔中排出，暂时附着在物体上，随后又分离开，有时会分裂成更小的气泡。然而，我们仔细分析后发现，附着在物体上的小气泡并不会显著改变其运动系数（例如，颗粒的沉降速度$U_{\kern-1pt f}$变化不超过2％）。唯一的显著区别是，对于未预先湿润的新月形颗粒，气泡在进入相机视野之前改变了颗粒的方向。物体的运动和重新定向被记录在相机图像的时间序列中。所有图像都被传输到MATLAB中，并进行图像分析以识别其中的物体。图像处理步骤包括对RGB图像进行二值化、背景减除和阈值处理，以获得更清晰的物体轮廓。此外，还采用了噪声去除技术，以避免在高光照条件下由于阴影和物体反射导致的虚假物体检测。

4. 实验结果
4.1 初始方向和实验次数
在实验中，我们控制了颗粒在液体表面的初始方向。如前所述，相机在颗粒从液体自由表面沉降约10厘米后开始记录，并记录它们从初始方向变化后的状态。在初步试验中观察到，所有形状在不同初始配置下都趋向于一个稳定的方向。然后，我们对液面初始方向有三种不同的情况进行了实验：(i)稳定的；(ii)倒立的；(iii)倾斜的。选择“稳定”的初始方向，即$x=X$、$y=Y$和$z=Z$，以确认这一方向确实保持不变。我们进行了27次这样的实验，结果在所有实验中颗粒的方向都没有随时间改变。“倒立”的初始方向，即$x=-X$、$y=Y$和$z=-Z$，目的是观察颗粒围绕$y=Y$轴的旋转（如图1a中的情况A所示），以及其质量中心在$xz$平面（即$XZ$平面）上的移动。这种运动可以用来根据理论关系式(2.9)测量旋转-平移移动系数$\mu _{51}$。“倾斜”的初始方向，即$x=X$，且$z$轴与$Z$轴之间的角度$\theta$大于$\pi /2$且接近$\pi$，目的是观察颗粒围绕$x=X$轴的旋转（如图1b中的情况B所示），从而根据理论关系式(2.10)测量旋转-平移移动系数$\mu _{42}$。在这种情况下，质量中心在$yz$平面（即$YZ$平面）上移动。表2列出了每种类型颗粒的实验次数，包括稳定的、倒立的和倾斜的初始方向。在155次实验中，我们观察到颗粒重新定向到一个稳定的、静止的方向（倒立初始方向有110次，倾斜初始方向有45次）。此外，还有27次实验中各种类型的颗粒在整个实验过程中保持稳定的方向。因此，我们总共进行了182次实验。

表2. 对于所选初始方向的物体，分别列出了实验次数以及用于重新定向定量分析的试验次数。物体名称：C代表圆锥体；A代表箭头；M代表新月形；R$i$代表开口宽度分别为$i = 1, 2, 3, 4$毫米的 ring。对于新月形颗粒，没有星号的数字代表使用未预先湿润颗粒的试验，带星号的数字代表使用预先湿润颗粒的试验，分别用于计算$\mu _{42}$和$\mu _{51}$。所有报告的实验均在Shekhar等人（参考文献Shekhar, Mirajkar, Zdybel, Melikhov and Ekiel-Je?ewska2025）的存储库中。在其中一些实验中，初始条件引发了复杂的运动。附录A中展示了示例。从这些实验中提取旋转-平移移动系数的数值较为困难。幸运的是，与我们预期的情况相符，所选的初始条件通常会导致两种特殊运动之一，分别标记为情况A和情况B，并在图1中说明。对于这两种特殊情况，动力学的解都是解析的且非常简单（如(2.9)和(2.10)所示）。因此，我们从所有实验中选择了情况A和情况B，并用它们来轻松计算$\mu _{51}$和$\mu _{42}$。

4.2 稳定和倒立初始方向的结果
我们记录了在110次实验中倒立初始方向的沉降圆锥体、新月形颗粒、箭头和开口环的重新定向情况。我们使用了预先湿润和未预先湿润的颗粒。对于同一类型的颗粒进行的所有实验，观察到颗粒总是趋向于相同的稳定方向并保持在该方向很长时间。在某些情况下，颗粒的运动比预测的更复杂，涉及角速度的其他分量，有些实验由于只能部分捕捉到重新定向而无法进行定量分析。为了定量分析圆锥体、箭头和开口环的运动，我们选择了颗粒仅围绕$y=Y$轴旋转的实验（如图1a中的情况A所示）。这种选择允许我们应用简单的演化方程(2.9)来计算时间依赖的倾斜角度$\theta$，并确定无量纲的旋转-平移移动系数$\mu _{51}$。然而，对于未预先湿润的新月形颗粒，尽管进行了27次实验，但没有这样的实验。相反，在7次实验中，新月形颗粒仅围绕$x=X$轴旋转（如图1b中的情况B所示），我们选择了这些实验，并将测量到的时间依赖角度$\theta$与(2.10)进行比较，从而得出另一个无量纲的旋转-平移移动系数$\mu _{42}$。关于预先湿润的新月形颗粒的结果将在后面呈现。为了展示形状随时间重新定向的情况，我们在单个图像中结合了两个相机处理后的图像序列，时间间隔$\Delta \tilde {t}$相同（箭头为1秒，新月形为2秒，圆锥体和所有环为4秒）。每种形状的记录时间为：圆锥体68秒；箭头40秒；新月形105秒；开口为1毫米的环180秒；开口为2毫米的环188秒；开口为3毫米的环190秒；开口为4毫米的环196秒。使用MATLAB库的imfuse函数进行了叠加，保持所有图像的相同比例和大小。所有实验中沉降颗粒的快照时间序列都在存储库中呈现（Shekhar等人，参考文献Shekhar, Mirajkar, Zdybel, Melikhov and Ekiel-Je?ewska2025）。在本文中，我们展示了其中的一些结果，将黑白图像反转并裁剪宽度，只保留每幅图像中央垂直线附近的$\pm 3$厘米范围。图4显示了三次实验中的快照，展示了单个沉降物体的时间依赖位置和方向：(a)圆锥体，(b)箭头，(c)新月形。图像是由两台相机（每个面板的上半部分和下半部分）同时记录的。重力方向朝左，颗粒从右向左移动。图4和图5展示了“选定的倒立”实验中各种形状的重新定向行为的快照示例。附录A中展示了从倒立初始方向开始的三个维度（3-D）动力学的两个例子，这些情况具有非零且随时间变化的三组分角速度。图5展示了四个实验中的快照，展示了开口分别为1毫米、(b) 2毫米、(c) 3毫米、(d) 4毫米的沉降环的时间依赖位置和方向。图像是由相机1（每个面板的上半部分）和相机2（每个面板的下半部分）同时记录的。重力方向朝左，颗粒从右向左移动。对于选定的实验，在准确识别图像序列中的物体形状后，对相机1（$XZ$投影）的图像进行了定量分析，除了未预先湿润的新月形颗粒，其图像是由相机2（$YZ$投影）分析的。首先，使用MATLAB的内置函数regionprops在每个时间点$\tilde {t}$识别出物体的几何特性，包括物体的质心位置（以像素为单位），然后通过校准因子转换为毫米，得到$(X_c, Z_c)$，以及物体的倾角$\theta$，这是根据二值图像的二阶中心矩来确定的，具体方法见附录B.1。通过将图像分析步骤中的阈值变化$\pm 5\,\%$，确定了使用regionprops识别物体中心$Z_c$时的不确定性为$\pm 0.2$像素。此外，$Z_c$的不确定性还来自其他因素；例如，校准因子的微小变化。需要注意的是，$Z_c$并不完全等于物体的质心位置的垂直坐标。倾角$\theta$的不确定性是根据附录B.2和B.3中的方法来确定的。时间依赖的垂直速度分量$U_Z(\tilde {t})$是通过计算连续图像之间的时间差$Z_c(\tilde {t}+\Delta \tilde {t})-Z_c(\tilde {t})$得到的，时间间隔约为$\Delta \tilde {t} \, \approx 1$秒。图6显示了实验中物体中心速度的垂直分量$U_Z$与物体垂直位置$Z_c$的关系，这些实验数据来自图4和图5。负值表示测量方向是从罐子顶部到底部。将图6与图4和图5中的快照序列进行比较，我们观察到对于圆锥和新月形状的物体，当它们水平定向时（即$\theta = \pi /2$），$U_Z$会减小并达到最小值，而当$\theta \rightarrow 0$时，$U_Z$会增加。相反，对于箭头形状的物体，当它们水平定向时（即$\theta = \pi /2$），$U_Z$会增加并达到最大值，然后随着$\theta \rightarrow 0$而减小。这与众所周知的规则一致：垂直定向的细长颗粒比水平定向的颗粒沉降得更快。图6还展示了图4和图5中实验实验中物体中心沉降速度的垂直分量$U_Z$。物体名称：C代表圆锥；A代表箭头；M代表新月；R$i$代表开口宽度分别为$i = 1, 2, 3, 4$毫米的环。

图7展示了图4和图5中实验实验中不同形状颗粒的重定向过程。实验数据（符号）的时间依赖性$\tan (\theta /2)$用指数函数（4.1）来近似（实线），其参数$A$和$\tilde {\tau }$的拟合值列在附录C中。不同物体的图表在时间上进行了调整，以便在相同的时间点$\tan (\theta /2)=1$。图7显示了颗粒随时间重定向的情况。在这里，我们将$\tan (\theta /2)$作为维度时间$\tilde {t}$的函数进行绘制，这些实验数据在图4-6中展示和分析了。函数$\tan (\theta /2)$的不确定性是根据附录B.3中描述的标准误差传播方法计算得出的$\theta$的不确定性得出的。为了更好地比较，我们将不同物体的图表在时间上进行了调整，使它们在相同的时间点达到$\tan (\theta /2)=1$。为了确定旋转-平移迁移系数，我们使用了理论推导的时间依赖性（2.9）–（2.10）对$\tan (\theta /2)$进行拟合。因此，我们使用指数衰减函数（4.1）来拟合实验数据：\[ \tan (\theta /2)= A \exp ( - \tilde {t} / \tilde {\tau } ) \]。其中，振幅$A$和特征时间$\tilde {\tau }$是使用MATLAB的非线性数据拟合函数fittype确定的。根据颗粒的不同，特征时间的倒数$1/\tilde {\tau }$要么是系数$1/|\tilde {\tau }_{51}|$，要么是$1/|\tilde {\tau }_{42}|$（见（2.11）。这些拟合参数对应于图7中的实线（每种物体类型的单个实验示例），可以在附录C中找到。图7表明，特征重定向时间取决于颗粒的形状：箭头的最小，新月形的次之，圆锥形的更大，而开口宽度为1毫米的环形的最大。对于初始方向被倒置的新月形颗粒，我们还观察并分析了其质心在$XZ$平面上的运动（见图1中的情况A）。图8展示了具有初始倒置定向的预湿润新月形颗粒在实验中的时间依赖性位置和定向的快照。这种特殊的动态现象在10次具有初始倒置定向的实验中有5次被观察到。

对于初始方向倾斜的物体，我们进行了45次实验。在附录A中，我们提供了一个从初始倾斜方向开始的3D动力学示例，其中角速度的三个分量不为零且随时间变化。图9展示了两个具有初始倾斜方向的实验的快照，显示了单个沉积物体的时间依赖性位置和定向：(a)箭头；(b)开口宽度为1毫米的环。图像是由两个相机（每个面板的顶部和底部）同时记录的。重力方向为左侧，颗粒从右向左移动。如前所述，在某些非预湿润的圆锥形、箭头形和新月形颗粒的实验中，我们观察到有一个小气泡从物体上的孔中释放出来，暂时附着在物体上，然后分离，有时会分裂成更小的气泡。然而，我们仔细分析后发现，附着在物体上的小气泡并不会显著改变迁移系数（例如，颗粒的沉降速度不会减少超过2%）。唯一的显著区别是，对于非预湿润的新月形颗粒，气泡在进入相机视野之前改变了颗粒的定向。

对于倾斜初始方向的物体，我们进行了45次实验。在附录A中，我们提供了一个从初始倾斜方向开始的3D动力学的示例，其中角速度的三个分量不为零且随时间变化。图9展示了两个具有初始倾斜方向的实验的快照，显示了单个沉积物体的时间依赖性位置和定向：(a)箭头；(b)开口宽度为1毫米的环。图像是由两个相机（每个面板的顶部和底部）同时记录的。重力方向为左侧，颗粒从右向左移动。为了进行定量分析，我们选择了颗粒仅绕$x=Y$轴旋转的实验（见图1中的情况B），并使用了来自相机2的图像。图9展示了两个选定实验的快照。我们使用（2.10）来确定旋转-平移迁移系数$\mu _{42}$，方法和之前描述倒置实验的方法相同。

为了提取每种物体类型的平均值，我们从初始方向的倒置类别中选择了颗粒仅绕$x=Y$轴或仅绕$x=X$轴旋转的实验；从初始方向的倾斜类别中选择了颗粒仅绕$x=X$轴旋转的实验（见图1中的情况A和B）。我们称这些实验为“选定的倒置”和“选定的倾斜”。为了计算平均速度，我们考虑了选定的倒置和静止实验（实验数量见表2）。表3展示了平均参数。对于特定颗粒类型的终端垂直速度$U_{f}$，方法是首先计算每个“静止”和每个“选定的倒置”实验的垂直速度$U_{Z,j}(\tilde {t})$的时间平均值及其标准差$s_{j}$。对于“静止”的实验，$U_{Z,j}(\tilde {t})$是在整个时间范围$\tilde {t}$上平均得出的。对于“选定的倒置”实验，$U_{Z,j}(\tilde {t})$是在物体达到稳定静止状态后的时间$\tilde {t}\geqslant \tilde {t}_o$内平均得出的。在这些情况下，每次实验的$\tilde {t}_o$是手动确定的。然后，使用所有选定实验的时间平均垂直终端速度$U_{Z,j}(\tilde {t})$来计算平均终端速度$U_{\kern-1pt f}$及其标准差$s_{ta}$。表3中报告的$U_{\kern-1pt f}$的不确定性是每个实验的标准差$s_{j}$和$s_{ta}$中的最大值。然后根据（2.12）从（2.12）计算了无量纲值$\mu _{33}$，不确定性是使用误差传播方法计算的。

表3展示了平均参数：在静止配置下的颗粒速度$U_{\!f}$、系数$1/\tilde {\tau }_{51}$或$1/\tilde {\tau }_{42}$，以及以质心为参考点的无量纲迁移系数。颗粒类型：C代表圆锥；A代表箭头；M代表新月；R$i$代表开口宽度分别为$i = 1, 2, 3, 4$毫米的环。表4显示了我们实验中的雷诺数（Re）、斯托克斯数（St）和阿基米德数（Ar）都远小于一。通过使用（4.1）拟合$\tan (\theta /2)$的时间依赖性，为每个选定的倒置或倾斜实验确定的无量纲系数$1/\tilde {\tau }_{51}$和$1/\tilde {\tau }_{42}$，随后对所有选定的实验分别进行了平均，平均值及其标准差见表3。然后根据（2.11）计算了无量纲的旋转-平移迁移系数$\mu _{42}$和$\mu _{51}$，不确定性是使用误差传播方法评估的。

图10展示了实验中使用的圆锥形、箭头形和开口环的形状，这些形状被表示为由$N$个接触的相同球形颗粒组成的集合体。颗粒在$\theta =0$时处于静止稳定配置，$H$方向与重力方向一致。(a)圆锥形，$N=6189$；(b)箭头形，$N=4896$；(c)开口宽度为1毫米的环，$N=5182$；(d)开口宽度为3毫米的环，$N=5110$。图11展示了新月形的形状，也被表示为由$N$个接触的相同球形颗粒组成的集合体。图展示了三个辅助表面网格，用于界定模型的边界（见附录D.4）。新月形在$\theta =0$时处于静止稳定配置，$H$方向与重力方向一致。(a)新月形，$N=5513$；(b)新月形，正面和侧面视图。现在使用表1和表3中的值来计算雷诺数、斯托克斯数和阿基米德数，分别为${Re} ={ U_{\kern-1pt f} L}/{\nu }$、${St} = {Re} M/(V {\rho })$和${Ar} ={F}/{(\rho \nu ^2)}$，并将它们列在表4中。这些值都远小于一。

5.理论建模
5.1. 从珠子模型评估迁移系数
在本节中，我们提出了一种确定任意形状颗粒迁移系数的方法。这是基于珠子模型的。具体来说，颗粒形状被近似为由许多接触的球形颗粒组成的集合体，然后对斯托克斯方程进行多极展开以评估迁移矩阵。许多作者已经成功应用珠子模型来确定复杂颗粒形状（包括蛋白质、DNA或聚合物）的迁移系数，例如Cichocki和Hinsen（参考文献Cichocki and Hinsen1995）、Adamczyk等人（参考文献Adamczyk, Cichocki, Ekiel-Je?ewska, S?owicka, Wajnryb and Wasilewska2012）、Zuk, Cichocki和Szymczak（参考文献Zuk, Cichocki and Szymczak2018）以及其中的其他参考文献。现在，我们将把珠子模型应用于实验中使用的颗粒，并证明该模型的准确性足以估计特征重定向时间并唯一确定$\mu _{51}$和$\mu _{42}$的符号，从而对颗粒进行分类。实验中使用的每个刚性颗粒都被近似为大量直径为$d$的相同球形颗粒的集合体。珠子集合体的形状展示在图10和图11中。构建细节在附录D中有说明。对于开口直径为2毫米和4毫米的环状物体，其珠状集合体的构建方式与图10(c,d)中所示的开口直径为1毫米和3毫米的环状物体的构建方式类似。通常，在这项工作中，一个颗粒是由超过4800个球形珠子组成的。关于如何将实验形状表示为珠状集合体的详细信息，请参见附录D。为了观察随着珠子数量的增加，移动系数饱和度的变化情况，也使用了较少数量的珠子。对于这些珠状集合体，(2.4)中定义的平移-平移移动系数和旋转-平移移动系数是通过斯托克斯方程的多极展开来计算的，并且经过润滑修正以加速收敛过程。这种多极方法的版本是由Cichocki等人（参考文献Cichocki, Felderhof, Hinsen, Wajnryb and Bl/awzdziewicz 1994；参考文献Cichocki, Ekiel-Je?ewska and Wajnryb 1999；Ekiel-Je?ewska & Wajnryb 2009；以及Cichocki, Ekiel-Je?ewska & Wajnryb 2014）开发的，并已在Hydromultipole数值代码中实现。我们在多极展开中截断了$LL=2$阶。我们保持较低的截断阶数，以便能够对更多数量的珠子进行计算，并更准确地再现颗粒的形状。在(2.4)中定义的无量纲移动系数的计算结果列在表5中。数值计算确认，对于所有对象，都有$\mu _{51} \,\mu _{42}\lt 0$。因此，基于Joshi & Govindarajan（参考文献Joshi and Govindarajan 2025）和Ekiel-Je?ewska & Wajnryb（参考文献Ekiel-Je?ewska and Wajnryb 2009a）所做的理论稳定性分析，我们得出结论：在本研究中调查的所有颗粒类型中，最终的静止方向是稳定的。在第5.2节中，我们将讨论珠子模型的准确性以及移动系数及其符号对形状变化的敏感性。在第6节中，我们将比较数值和实验结果，并证明两者之间的差异足够小，可以唯一地分类这些颗粒为沉降颗粒。珠子模型可以应用于不同的颗粒形状，并且无需进行实验即可确定移动系数。表5列出了基于质心作为参考点，在(2.4)中定义的无量纲移动系数，这些系数是根据我们实验中使用的颗粒的珠子模型计算得出的。最后，在第7节中，我们将进行一项专门的参数化数值研究，其中我们改变环状开口的半角$\Delta \vartheta$，并计算相应环状碎片的$\mu _{51}$和$\mu _{42}$。结果表明，在整个参数范围内，旋转-平移移动系数的符号保持不变，环状碎片保持了沉降颗粒的动力学特性。

5.2 珠子模型的准确性
在本节中，我们将讨论移动系数的三个不确定性来源：(i) 给定珠状集合体的多极展开的截断；(ii) 实验形状与相应珠状集合体形状之间的差异；(iii) 容器壁和自由表面的影响。(i) 对于刚性颗粒，经过润滑修正后的多极展开的准确性随着多极阶数的增加而迅速提高（参见例如Cichocki等人参考文献Cichocki, Felderhof, Hinsen, Wajnryb and Bl/awzdziewicz 1994；参考文献Cichocki, Ekiel-Je?ewska and Wajnryb 1999）。我们的收敛性测试与先前研究的估计结果一致。对于由741个珠子组成的、近似锥形颗粒的集合体，使用$LL=2$阶和$LL=5$阶计算得到的移动系数之间的差异，对于$\mu _{33}$小于0.2％，对于$\mu _{42}$和$\mu _{51}$小于0.7％。这一结果证明了使用$LL=2$阶是合理的。(ii) 现在我们将讨论珠状集合体表示实验形状的准确性。首先，我们使用了计算集群内存允许的最大珠子数量$N$。根据形状的不同，$N$的值在4870到6189之间变化。随着$N$的增加，数值不确定性通常会降低。作为这种不确定性的一个粗略估计，我们使用了在最高可用分辨率下计算的移动系数与之前较粗糙模型下的移动系数之间的相对差异（详见附录D.5）。计算出的差异不超过15％。此外，珠子模型预测中还有其他不确定性来源。在我们的计算模型中再现精确的实验形状存在实际挑战。特别是，对于新月形颗粒，无法获得精确的边界规格，因此需要基于近似测量尺寸的模型。同样，对于锥形和箭头形颗粒，尽管整体形状较为简单，但诸如角部圆润等细微特征很难精确再现。此外，珠子放置的标准（例如要求珠子中心或整个珠子体积位于给定边界内）对移动系数的值影响较小。增加珠子的数量可以减小单个珠子的半径，从而更接近实际颗粒表面的描述，但数值计算对最小珠子半径有所限制。然而，在这项工作中并不需要高精度。上述准确性足以估计移动系数的值，而且珠子模型适用于本研究的范围。(iii) 最后，让我们估计壁效应，重点讨论壁对位于容器中心的颗粒的移动系数$\mu _{33}$的影响。与无界流体中的移动系数$\mu _{33}^{(\infty )}$相比，我们定义壁修正为$\delta \mu _{33}= (\mu _{33}^{(\infty )}-\mu _{33})/\mu _{33}^{(\infty )}$。根据Goldman等人（参考文献Goldman, Cox and Brenner 1967）以及Lee, Chadwick & Leal（参考文献Lee, Chadwick and Leal 1979）的研究结果，对于半径为$a$的球体，如果它平行或垂直于距离颗粒中心$D$的平面壁移动，那么$\delta \mu _{33}\approx ka/D$，其中$k=9/16$或$k=9/8$；对于向自由表面移动的球体，$k=3/4$。将我们颗粒的流体动力学半径估计为$a\approx 2.5$毫米，并使用界面叠加方法，我们得到$\delta \mu _{33}\approx 7.5\,\%$。或者，我们的容器可以被近似为一个圆柱体，此时可以应用Sano（参考文献Sano 1987）对于沿圆柱体轴线移动的半径为$a$的球体的分析结果，$\delta \mu _{33}\approx ka/D$，其中$D$是圆柱体的半径，$k=2.1$是圆柱体的长宽比。应用这个公式得到壁效应的估计值约为$\delta \mu _{33}\approx 5\,\%$。因此，壁效应小于珠子模型的估计精度。

6. 移动系数的实验值与数值值的比较
在展示所有颗粒的实验测量（表3）和数值估计（表5）的移动系数之间的比较分析之前，我们需要明确为什么确定旋转-平移移动系数$\mu _{51}$和$\mu _{42}$的值及其符号是重要的。原则上，仅仅观察到实验中的特定颗粒形状在某个方向上保持很长时间并不足以推断出该方向是稳定的。颗粒也可能在一个不稳定的静止方向上保持很长时间，直到受到扰动而重新定向。相反，通过实验或数值方法确定$\mu _{51}$和$\mu _{42}$的值及其符号，我们可以判断该颗粒是否为沉降颗粒，如果是沉降颗粒，那么它的稳定静止方向是什么，以及其特征性的重新定向时间尺度是什么。因此，我们的主要目标是估计旋转-平移移动系数的值，而不仅仅是提供精确的值。这些估计有两个目的：(i) 大致预测重新定向时间；(ii) 判断给定形状是否表现为沉降颗粒。关于后者，对于$\mu _{42}$和$\mu _{51}$，理论和实验值之间的绝对差异小于绝对值本身。这确保了$\mu _{42} \,\mu _{51}$的符号得到保留，这足以正确地将颗粒识别为沉降颗粒。此外，如表6所示，实验和理论旋转-平移移动系数之间的相对差异不超过25％，这对于我们的目标来说是足够小的。

7. 参数化数值研究：开口控制的环状物体和移动趋势
为了测试我们的结果在超出之前研究的特定几何形状时的鲁棒性，我们在$xz$平面上对两组具有可调开口的环状碎片进行了形状控制的扫描（见图12(a)）。这些环状物体可以是圆形的（半径为$a$）或椭圆形的（长宽比为$a/b\approx 1.43$，$a$方向沿$x$轴）。圆形和椭圆形环状物体在其整个弧长上都具有均匀的圆形截面，圆形环的半径为$R = 0.122a$，椭圆形环的半径为$R = 0.092a$。每个环状物体都通过接触直径为$d$的球形珠子来离散化，这些珠子的中心位于一个间距为$d$的立方体网格上；只有当整个珠子的体积完全包含在管状环内时才保留该珠子。为了描述不依赖于开口$\Delta \vartheta$的珠子模型的分辨率，我们报告了横截面占据率$N_\perp$，定义为在垂直于局部中心线的厚度为$d$的切片中包含的珠子中心数量（即在半径为$R$的圆盘内）。对于圆形和椭圆形环，当$R=3.6\,d$时，每个横截面有$N_\perp \approx 42$个珠子，横截面取自$x=0$处的$yz$平面。作为一个参考尺度，$15^\circ$的角度段包含圆形环$N_{15^\circ }=325$个珠子，椭圆形环$N_{15^\circ }=451$个珠子。在参数化研究中，我们保持相同的$a,\,b$和$R$；控制参数是附录D.3中定义的开口半角$\Delta \vartheta$，如图12(a)所示。尺寸$L$定义为沿$x$轴的最外层珠子边缘之间的距离。对于圆形环，$L/a$的值范围从1.09到2.21；对于椭圆形环，范围从1.05到2.18。在每种情况下，最小值对应于$\Delta \vartheta = 150^\circ$，最大值对应于$\Delta \vartheta = 90^\circ$。图12显示了开口控制环状碎片的参数化研究。(a) 几何形状和符号。(b–d) 移动系数作为$\Delta \vartheta$的函数：(b) $\mu _{33} \,({R}/{L})$，(c) $\mu _{42}\,( {R}/{L})^2$和(d) $\mu _{51} \,({R}/{L})^2$。蓝色虚线曲线代表圆形环；橙色虚线曲线代表椭圆形环（$a/b\approx 1.43$）。红色星星（R1–R4）代表表5中给出的开口宽度分别为$1,2,3$和$4$毫米的环状物体的数值。对于每个$\Delta \vartheta$，我们使用相同的归一化方法、相同的多极阶数和数值容忍度来计算移动系数。为了比较不同形状，我们报告了乘以$R/L$的移动系数$\mu _{33}$，而$\mu _{51}$和$\mu _{42}$则乘以$(R/L)^2$。通过这种方式，我们比较了相同$R$的环状碎片的移动系数。图12总结了圆形和椭圆形环的开口扫描结果。图12(b)显示，在$60^{\circ } \leqslant \Delta \vartheta \leqslant 150^{\circ }$的范围内，垂直移动系数$\mu _{33}({R}/{L})$呈单调增加。旋转-平移移动系数（图12 c,d）则不具有单调性：在$\Delta \vartheta \approx 120^\circ \!-\!140^\circ$的范围内，$\mu _{42}({R}/{L})^2$和$|\mu _{51}|({R}/{L})^2$出现明显的峰值。这些极值的存在是可以理解的，因为在$\Delta \vartheta \rightarrow 0$（封闭环）和$\Delta \vartheta \rightarrow 180^{\circ }$（完整圆）的极限情况下，旋转-平移移动系数消失。关键的是，$\mu _{42}$ 和 $\mu _{51}$ 的符号保持不变，这表明这些物体在整个家族中仍然保持其沉降特性——它们的优选方向和旋转感得到了保留，尽管它们的振幅有所不同。圆形环（蓝色）系统性地放大了峰值，并略微改变了它们相对于椭圆形环（橙色）的位置，但定性趋势并未改变。红星（R1–R4）标记了表5中的数值结果，这些结果表明环的开口宽度分别为$=1,2,3$毫米，长宽比为$a/b$ = 1.11, 1.22, 1.35和1.43；在相应的$\Delta \vartheta$下绘制时，它们显示出相似的倾向。另一个有用的观察结果是，开口控制环家族的三维方向动力学至少在定性上与之前发现的新月形物体的动力学特征重合。三维方向动力学的基本特征取决于$\mu _{51}/\mu _{42}$；而对$\mu _{51}$的依赖仅通过时间缩放来体现。特别是对于上面讨论的新月形珠子模型，我们得到$\mu _{51}/\mu _{42}=-0.44$，而对于开口宽度为$\Delta \vartheta =120^\circ$的椭圆形环，相应的值为$\mu _{51}/\mu _{42}=-0.45$。这种近乎巧合的现象与观察到的方向动力学的相似性一致，并支持了这样一个解释：一个足够宽的开口环在动力学上可能接近新月形物体的情况。图13显示了不同物体的稳定静止方向，对应于$\theta =0$。这些形状和它们的方向是从实验图像中获得的。同时，不应该过度解读这个比率作为形状描述符：即使对于几何形状不同的物体，$\mu _{51}/\mu _{42}$的相似值和相同的符号也可能是存在的。例如，具有不同高径比（height-to-base-diameter）的圆锥在家族中都具有$\mu _{51}/\mu _{42}=-1$。此外，一个开口非常小的环的$\mu _{51}/\mu _{42}$也非常接近$-1$，但其形状与圆锥不同。因此，比率$\mu _{51}/\mu _{42}$对于分类动力学是有用的，但它并不能唯一确定粒子的几何形状。

8. 结论
在这项工作中，我们通过实验和数值方法研究了在雷诺数远小于1的情况下，刚性各向异性颗粒在粘性流体中的重力沉降行为。我们调查了具有不同形状的颗粒的动力学行为，这些形状包括圆锥、箭头形、新月形和开口宽度为4种不同尺寸的环形物体。所有这些物体都具有两个垂直的对称平面。总的来说，这类形状有三类，它们具有质量不同的动力学特性，Joshi和Govindarajan（参考文献Joshi and Govindarajan2025）对此进行了分类。这项工作的目标是分析其中一类：即那些趋于稳定静止方向的沉降颗粒。Ekiel-Je?ewska和Wajnryb（参考文献Ekiel-Je?ewska and Wajnryb2009a）研究了钟形的颗粒动力学，而Candelier等人（参考文献Candelier, Gustavsson, Sharma, Sundberg, Pumir, Bagheri and Mehlig2025）则研究了半圆形和四分之一圆形颗粒的动力学。实际上，在本研究中进行的实验中，所选物体重新定向并达到了稳定的静止方向，如图13所示，这些结果是通过实验图像得出的。如第2节所述，沉降颗粒的动力学特性由两个旋转-平移移动系数$\mu _{51}$和$\mu _{42}$决定，它们的符号相反。这些系数的值揭示了颗粒重新定向的速度。因此，在这项工作中，我们的目标是在实验中测量这些系数。我们选择了这样的初始方向：由于对称性，这些方向可能导致颗粒绕$y=Y$轴或$X$轴旋转（对于倒置的初始方向），以及绕$X$轴旋转（对于倾斜的初始方向），如图1(a,b)所示。我们为倒置和倾斜的初始方向选择了这样的对称实验，并分别确定了$\mu _{51}$和$\mu _{42}$的值，见表3。同样的实验方法也可以用于研究颤动型和漂移型颗粒。接下来，我们将实验中使用的颗粒形状近似为接触在一起的相同珠子的集合体。使用多极方法，我们评估了这些集合体的移动系数。实验中测得的移动系数（见表3）与数值计算得出的移动系数（见表5）之间的差异不超过25%（详见第6节）。因此，本工作中提出的理论模型可以作为估计具有所需形状的物体的$\mu _{51}$和$\mu _{42}$值的有用工具，而无需进行实验。这些值可以预测动力学的类型，特别是判断一个物体是否为沉降颗粒以及它将以多快的速度旋转到稳定方向。这些信息对于实际应用非常重要。在某些实验中，我们观察到了三维动力学行为，如图14(b,c)的附录A所示。主要原因是，颗粒进入相机视野时的方向不属于A或B情况。形状的轻微不对称可能是某些实验中观察到三维动力学的另一个原因。有趣的是，对于非手性形状（如圆锥、新月形、箭头形和开口环），会发生手性运动。这种效应是由于相对于质心计算出的旋转-平移移动系数不为零造成的。对于具有两个垂直对称平面的非手性形状（称为颤动型颗粒），观察到了更复杂的手性运动（参见Miara等人（参考文献Miara, Vaquero-Stainer, Pihler-Puzovi?, Heil and Juel2024）、Vaquero-Stainer等人（参考文献Vaquero-Stainer, Miara, Juel, Pihler-Puzovi? and Heil2024）、Joshi和Govindarajan（参考文献Joshi and Govindarajan2025）以及Vaquero-Stainer等人（参考文献Vaquero-Stainer, Miara, Juel, Heil and Pihler-Puzovi?2026）的研究）。我们的结果是针对在雷诺数远小于1的高粘性流体中沉降的毫米级颗粒获得的。根据相似性原理（参见例如Homsy等人（参考文献Homsy2008）），所得到的动力学结果可以很容易地重新应用于在相同雷诺数下的水基溶液中几何相似的微粒的运动。最后，我们的结果表明，经过一段时间后，沉降的微沉降颗粒的稀薄悬浮液将形成一个有序的各向异性介质。这样的系统在理论和实践上都具有兴趣。

资金支持
这项工作部分得到了国家科学中心（National Science Centre）的资助，授予编号为UMO-2021/41/B/ST8/04474的资助。

利益声明
作者们声明没有利益冲突。

作者贡献
C.S.、P.Z.和Y.M.对此工作做出了同等程度的贡献。

数据可用性声明
支持本研究发现的数据可以在RepOD（开放数据存储库）中公开获取，网址为https://doi.org/10.18150/MIL6UG。

附录A. 沉降颗粒的三维重新定向
图14展示了最初处于“倒置”或“倾斜”方向的颗粒的三维重新定向示例。在本工作中，这种复杂动力学还没有被定量分析。显然，对于图14中显示的每个实验，角速度的时变非零分量$X$、$Y$和$Z$在某些时间范围内存在。

附录B. 从二维图像测量方向角$\theta$
B.1 使用二阶中心矩估计方向角
我们采用了一种基于二阶中心矩的统计方法来确定二进制图像中物体的方向角，这种方法在Burger和Burge（参考文献Burger and Burge2009）的书中有所描述。阶数为$(p+q)$的归一化中心矩定义为：
$$
\begin{eqnarray}
\mathcal{M}_{pq} = \frac {1}{N_{\mathcal{B}}^{(p+q+2)/2}} \sum _{i \in \mathcal{B}} (x_i - \bar {x})^p (z_i - \bar {z})^q,
\end{eqnarray}
$$
其中$(x_i, z_i)$表示组成二维物体$\mathcal{B}$的像素坐标，$(\bar {x}, \bar {z})=\frac {1}{N_{\mathcal{B}}}\sum _{i\in \mathcal{B}} (x_i, z_i)$是$\mathcal{B}$的质心，$N_{\mathcal{B}}$是$\mathcal{B}$中的像素总数。协方差矩阵$\boldsymbol{\varSigma }$是使用$\mathcal{B}$的二阶矩构建的：
$$
\begin{eqnarray}
\boldsymbol{\varSigma }=\left (\begin{matrix} \mathcal{M}_{20} & \mathcal{M}_{11} \\
\mathcal{M}_{11} & \mathcal{M}_{02}
\end{matrix} \right ).
\end{eqnarray}
$$
上述矩阵通过相似变换$\tilde {\boldsymbol{\varSigma }}= \boldsymbol{R}^{-1}\boldsymbol{\varSigma } \boldsymbol{R}$进行对角化，其中$\boldsymbol{R}$是一个旋转矩阵，它将坐标轴逆时针旋转一个角度$\alpha$，该角度由下式给出：
$$
\alpha = \frac {1}{2}\arctan \left (\frac {2 \mathcal{M}_{11}}{\mathcal{M}_{20}-\mathcal{M}_{02}}\right ),
$$
我们在测量中使用的物体$\mathcal{B}$的方向角$\theta$由下式给出：
$$
\theta = \frac {\pi }{2} + \alpha,
$$
这两个角度都是以弧度为单位的。图15展示了一个示意图，其中标出了计算出的等效椭圆和指示的倾斜角$\alpha$。椭圆以物体的质心为中心，其主要轴也已标出。

B.2 二维物体的等效椭圆
所研究的物体$\mathcal{B}$可以用一个等效椭圆来表示，该椭圆的长半轴为$a$，短半轴为$b$，它具有与物体$\mathcal{B}$相同的协方差矩阵（见图15）。对协方差矩阵进行对角化后，经过角度$\alpha$的旋转后，椭圆的长轴将变为水平方向，协方差矩阵表示为$\tilde {\boldsymbol{\varSigma }}=\mathrm{diag}(a^2,b^2)$。

B.3 方向角的不确定性
我们使用测量误差的传播来估计物体方向角$\theta =\pi /2 + \alpha$的测量不确定性：
$$
\sigma _{\theta }=\sqrt {\left (\frac {\partial \alpha }{\partial \mathcal{M}_{20}}\sigma _{\mathcal{M}_{20}}\right )^2+\left (\frac {\partial \alpha }{\partial \mathcal{M}_{02}}\sigma _{\mathcal{M}_{02}}\right )^2+\left (\frac {\partial \alpha }{\partial \mathcal{M}_{11}}\sigma _{\mathcal{M}_{11}}\right )^2,
$$
其中${\partial \alpha }/{\partial \mathcal{M}_{02}}=-({\partial \alpha }/{\partial \mathcal{M}_{20}})=2\mathcal{M}_{11} D^{-1}$，${\partial \alpha }/{\partial \mathcal{M}_{02}}=(\mathcal{M}_{20}-\mathcal{M}_{02}) D^{-1}$，$D=(\mathcal{M}_{20}-\mathcal{M}_{02})^2+4\mathcal{M}_{11}^2$。我们假设像素坐标$(X, Z)$遵循双变量正态椭圆分布$\mathcal{N}_2((0,0), \boldsymbol{\varSigma })$，其中协方差矩阵与(B3)中的相同，但我们对参数进行了适当的简化，即：
$$
\boldsymbol{\varSigma }=\boldsymbol{R}\tilde {\boldsymbol{\varSigma }}\boldsymbol{R}^{-1}=\left (\begin{matrix} a^2 \cos ^2\alpha +b^2 \sin ^2\alpha & (a^2-b^2)\sin \alpha \cos \alpha \\
(a^2-b^2)\sin \alpha \cos \alpha & a^2 \sin ^2\alpha +b^2 \cos ^2\alpha \end{matrix} \right ),
$$
其中$\alpha =\theta -\pi /2$与测得的物体方向角$\theta$有关。图15展示了这个等效椭圆和指示的倾斜角$\alpha$。我们可以直接写出 $\sigma _{\tan ({\theta}/{2})$ 的表达式，该表达式由 (B14) 给出：
\[
\sigma _{\tan{\kern-1pt} \theta/2}=\left |\frac {\mathrm{d} \tan \frac {\theta }{2}}{\mathrm{d}\theta }\right |\sigma _{\theta , est}=\frac {1}{2\sqrt {N_{\mathcal{B}}}}\frac {ab}{|a^2-b^2|}\left (1+\tan ^2\frac {\theta }{2}\right ).
\]
上述公式用于计算图 7 中的误差条。图像分析是使用 MATLAB 的内置函数 regionprops 进行的，该函数可以直接访问二值图像的属性。值得注意的是，这个函数允许通过 ‘orientation’ 属性获取角度 $\alpha$，以及通过 ‘MajorAxisLength’ 和 ‘MinorAxisLength’ 属性获取等效椭圆的轴长。

B.4 运动模糊的不相关性
由于实验中使用的相机的曝光时间有限，因此由于物体的垂直沉降而产生了运动模糊。我们假设物体的水平速度与其沉降速度相比可以忽略不计，从而可以将运动模糊效应仅在 $z$ 方向上考虑。与这种效应相关的标准差 $\sigma _{blur}$ 是通过将模糊建模为沿 $z$ 轴的均匀分布来获得的，其宽度为 $\textit{UTk}$，从而得到 (B15)：
\[
\sigma ^2_{blur} = \frac {(\textit{UTk})^2}{12}
\]
其中 $U$ 是沉降速度，$T$ 是曝光时间，$k$ 是一个缩放因子，用于根据相机校准将毫米测量距离转换为像素单位。运动模糊改变了二阶中心矩 $\mathcal{M}_{20}$ 的值，通常会导致其由于物体强度沿模糊方向的拉伸而被高估。存在运动模糊的情况下与不存在运动模糊的情况下的二阶矩之间的关系由 (B16) 给出：
\[
\mathcal{M}^{\textit{blur}}_{20}&=&\mathcal{M}_{20}, \\
\mathcal{M}^{\textit{blur}}_{02}&=&\mathcal{M}_{02}+ \sigma ^2_{blur}, \\
\mathcal{M}^{\textit{blur}}_{11}&=&\mathcal{M}_{11}.
\]
因此，物体的方向角会受到影响。我们可以通过以下方程将真实的方向角 $\alpha$ 与从模糊图像中提取的方向角 $\alpha _{blur}$ 关联起来：
\[
\tan 2 \alpha = \frac {\tan 2 \alpha _{blur}}{1+\frac {\sigma ^2_{blur}}{2 \mathcal{M}_{11}} \tan 2 \alpha _{blur}}\approx \frac {\sin 2 \alpha _{blur}}{\cos 2 \alpha _{blur}+\frac {\sigma ^2_{blur}}{a^2_{blur}-b^2_{blur}},
\]
其中 $a_{blur}$ 和 $b_{blur}$ 是从模糊图像中获得的等效椭圆的半轴长度。在没有运动模糊的情况下，估计的方向角与真实的方向角一致，即 $\alpha =\alpha _{blur}$。对于所研究的物体，$\sigma _{blur}^2/ (a_{blur}^2-b_{blur}^2 )\ll 1$，因为 $\sigma ^2_{blur}$ 的数量级是几个平方像素，而 $a_{blur}^2-b_{blur}^2$ 的数量级是几千平方像素。$\tan 2\alpha$ 的相对差异由 (B18) 给出：
\[
\frac {|\tan 2\alpha _{blur}-\tan 2\alpha |}{|\tan 2\alpha _{blur}|}\approx \frac {\sigma _{blur}^2}{a_{blur}^2-b_{blur}^2}\frac {1}{\cos 2\alpha _{blur}}
\]
并且对于几乎所有角度来说都小于 $1\,\%$，除了等于 $45^\circ \pm 2^\circ$ 的角度。因此，在本研究中可以安全地忽略运动模糊的效应。

附录 C. 适用于图 7 中实验数据的拟合参数 (4.1)
表 7 给出了理论函数 (4.1) 对图 7 中显示的数据的拟合参数，其中包括了七个不同形状粒子的示例实验试验。相应的决定系数 $R^2$ 值表明实验结果与理论预测之间有很好的一致性。在我们的研究中，决定系数 $R^2$ 是按以下公式计算的 (C1)：
\[
R^2=1-\frac {\sum _{i=1}^n w_i \left (y_i - \hat {y}_i\right )^2}{\sum _{i=1}^n w_i \left (y_i - \bar {y}_w\right )^2},
\]
其中 $y_i=\tan (\theta _i/2)$ 是在时间点 $\tilde {t}_i$ 处测量的数据，$\hat {y}_i=A \exp (-\tilde {t}_i/\tilde {\tau })$ 是具有拟合参数 $A$ 和 $\tilde {\tau }$ 的模型预测，$w_i = \sigma ^{-2}_{\tan (\theta _i/2)}$ 是与 (B14) 中的不确定性相关的统计权重。$C2$ 中的 $\bar {y}_w=\frac {\sum _{i=1}^n w_i y_i}{\sum _{i=1}^n w_i}$ 是数据点的加权平均值。这里，$n$ 表示测量总数。参数估计是在 MATLAB 中使用加权非线性最小二乘法进行的。

表 7. 对图 7 中七个示例实验试验数据的指数衰减拟合 (4.1) 的拟合参数 $A$ 和 $1/\tilde {\tau }$，以及显示拟合优度的决定系数 $R^2$。注意表 7 中单个试验的 $1/\tilde {\tau }$ 值与表 3 中的相应平均值之间的差异。

附录 D. 实验中使用的粒子珠子模型
D.1 圆锥体
为了构建代表实验中使用的圆锥形粒子的珠子集合，首先直接从物体的实验图像中识别出底角 $\alpha _{i}$：$\alpha _{i}=81.0 ^\circ \pm 0.2 ^\circ$。然后，使用识别的角度 $\alpha _{i}$ 来定义直圆锥的底直径 $L_{i}=4.80 \pm 0.05$ mm 和高度 $H_{i}=15.15 \pm 0.18$ mm。参数方程 $R(z) = ({(H_{i}-z) L_{i} }/{2 H_{i} }$ 通过关联其径向和方位角坐标 $R$ 和 $z$ 来描述圆锥的边界。直径为 $d$ 的球形珠子被放置在一个规则的三维立方晶格中，晶格常数等于球体直径 $d$，并对于球体中心位置的径向和方位角分量 $r$ 和 $z$ 有以下约束：$r \leqslant R(z)$，对于 $0 \leqslant z \leqslant H$，确保所有球体中心严格位于体积内，包括其表面。注意，理想的圆锥体必须在实验测量的高度 $H=10.45 \pm 0.05$ mm 处被截断，这是由 $z$ 的约束指定的。由于理想圆锥的底直径 $L_{i}$ 大于物体的直径 $L$，因此约束 $\sqrt {x^2+y^2} \leqslant L/2$ 移除了那些位于对称轴 $z=0$ 以外距离超过物体半径 $L/2$ 的珠子。约束 $\sqrt {y^2 + (z-H_{\textit{hole}})^2} \gt R_{\textit{hole}}$ 移除了位于通过圆锥体的圆柱孔中的珠子，该孔的半径为 $R_{\textit{hole}}$，并平行于底面并在距离底面 $H_{\textit{hole}}$ 的高度处与中心轴相交。底面的圆滑处理是使用约束 $\sqrt {x^2+y^2} \leqslant L/2 -d$ 实施的。顶部的圆滑处理是使用约束 $\sqrt {x^2+y^2+ ( z- (H-R_{top}) )^2} \leqslant R_{top}$ 实施的，其中参数 $R_{top}=0.25$ mm 代表用于匹配形状的圆滑半径。球体直径被系统地减小，直到流动系数达到近似饱和。结果得到的圆锥形物体的珠子模型表示由 6189 个球体组成，如图 10(a) 所示。

D.2 箭头形物体
为了构建实验中使用的箭头形物体的珠子模型表示，首先创建了一个长为 $L$、高为 $H$、宽为 $W$ 的长方体（立方体），这些尺寸对应于箭头形物体的长度、总高度和宽度（实验值在 § 3.2 中指定）。然后，该长方体被放置以满足以下约束：$-L/2 \leqslant x \leqslant L/2$、$0 \leqslant y \leqslant W$ 和 $0 \leqslant z \leqslant H$。接着，创建了一个晶格常数等于球体直径 $d$ 的规则三维立方晶格，使得其中一个晶格节点的坐标为 $(0,0,0)$。最后，模型被构建为只包含那些中心位于 $(x,y,z)$ 且在物体体积内的球体，这由以下约束定义：
\[
\begin{align}
0 \leqslant |x| \leqslant L/2, \qquad 0 \leqslant y \leqslant W, \qquad k \boldsymbol{\cdot }|x| \leqslant z \leqslant S + k \boldsymbol{\cdot }|x|,
\end{align}
\]
其中斜率 $k$ 是使用箭头形物体的中间高度 $S$ 定义的，$k=( {H-S})/({L/2})$。约束 $\sqrt {x^2 + y^2} \gt R_{\textit{hole}}$ 移除了位于贯穿物体的孔内的珠子，该孔的半径为 $R_{\textit{hole}}$。球体直径被系统地减小，直到所有流动系数的值达到近似饱和。结果得到的箭头形物体的珠子模型表示由 4896 个球体组成，如图 10(b) 所示。

D.3 开环
我们考虑一个形状为半径为 $R > 0$ 的圆形管道的刚性环形体，该管道沿着一个平面椭圆曲线中心排列。管道的准线是由 (D2) 给出的平面椭圆：
\[
\mathcal{C} = \left \{ \boldsymbol{r}(\vartheta ) = \left .\left ( a\cos \vartheta ,\, b\sin \vartheta ,\, 0 \right ) \,\right |\, \vartheta \in [0, 2\pi ) \right \},
\]
其中半长轴为 $a$，半短轴为 $b$。对于所考虑的环，$a$ 和 $b$ 的值是通过实验测量的。椭圆完全位于平面 $z = 0$ 内，参数 $\vartheta$ 逆时针增加。在椭圆上的每个点，我们定义一个局部正交框架 $\{ \boldsymbol{T}(\vartheta ), \boldsymbol{N}(\vartheta ), \boldsymbol{B}(\vartheta ) \}$，其中 $\boldsymbol{T}(\vartheta )$ 是单位切向量，$\boldsymbol{B}(\vartheta ) = \boldsymbol{e}_z$ 是垂直于椭圆平面的单位向量，$\boldsymbol{N}(\vartheta ) = \boldsymbol{B}(\vartheta ) \times \boldsymbol{T}(\vartheta )$ 是平面内的法向量。然后管道表面被参数化为 (D3)：
\[
\mathcal{S} = \left \{ \boldsymbol{X}(\vartheta , \varphi ) = \boldsymbol{r}(\vartheta ) + R \left .\left [ \cos \varphi \, \boldsymbol{N}(\vartheta ) + \sin \varphi \, \boldsymbol{B}(\vartheta ) \right ] \,\right |\, \vartheta \in [0,2\pi ),\, \varphi \in [0,2\pi ) \right \}.
\]
我们使用珠子模型对管道表面所包围的体积进行离散化。珠子中心的位置（通过珠子直径 $d$ 标准化）被放置在一个单位间距的均匀立方网格上，只有那些中心位于 $(x, y, z) \in \mathbb{R}^3$ 且在管道内的珠子被保留。如果一个珠子中心到椭圆曲线 $\boldsymbol{r}(\vartheta )$ 的距离小于或等于有效管道半径 $R-d/2$（只有整个周长位于管道边界内的珠子被包括在内），则认为该珠子在管道内，即 (D4)：
\[
\min _{\vartheta \in [0, 2\pi )} \left \| (x, y, z) - \boldsymbol{r}(\vartheta ) \right \| \leqslant R-\tfrac {d}{2}.
\]
这个条件是通过离散化参数 $\vartheta$ 并找到每个网格点上椭圆上的最近点来数值评估的。为了构建类似于实验中使用的开环，我们移除了位于椭圆上预先指定角度扇区附近的珠子子集。具体来说，我们定义一个中心角度位置 $\vartheta _0 = 90^\circ$，并移除了所有其中心点在椭圆上对应的角度 $\vartheta ^* \in [\vartheta _0 - \Delta \vartheta ,\, \vartheta _0 + \Delta \vartheta ]$ 之外的珠子。为了确保环中的间隙对应于一个物理上有意义的弧长 $\ell _{\textit{gap}}$，我们确定相应的角度宽度 $\Delta \vartheta$，使得 (D5)：
\[
\int _{\vartheta _0 - \Delta \vartheta }^{\vartheta _0 + \Delta \vartheta } \left \| \frac {\mathrm{d}\boldsymbol{r}}{\mathrm{d}\vartheta }(\vartheta ) \right \| \, \mathrm{d}\vartheta = \ell _{\textit{gap}.
\]
这是一个关于 $\Delta \vartheta$ 的非线性方程，因为弧长密度 (D6)：
\[
\left \| \frac {\mathrm{d}\boldsymbol{r}}{\mathrm{d}\vartheta }(\vartheta ) \right \| = \sqrt {a^2 \sin ^2\vartheta + b^2 \cos ^2\vartheta }
\]
取决于 $\vartheta$。我们使用二分法数值解决这个方程。给定一个目标弧长 $\ell _{\textit{gap}}$，确定相应的角度半宽度 $\Delta \vartheta$，以使从 $\vartheta _0 - \Delta \vartheta$ 到 $\vartheta _0 + \Delta \vartheta$ 的椭圆弧长度满足所需长度。这个程序保证了间隙围绕 $\vartheta _0$ 的旋转对称性，并确保了对开口的精确几何控制。主轴 $a$ 和 $b$ 的比率与实验测量值相关，即 $(a+R)/(b+R) \approx L/H$。同样，主轴半径与管道半径的比值与 $(a+R)/R\approx L/W$ 成正相关。在计算过程中，管道表面被重新缩放，从而可以控制该表面上包含的珠子数量 $N$。增加 $R$ 的值（以球体直径 $d$ 为单位）会重新缩放管道表面，从而增加模型中包含的球体数量，进而提高模型的精度。对于不同间隙宽度的环，采用了以下参数来匹配实验中使用的开口环的形状，并在 § 3.2 中进行了规定：(i) $a = 27.6$，$b = 24.9$，$R = 3.2$，$\Delta \vartheta = 10^\circ 47'$，对于间隙宽度为 1 mm 的环，得到 $N = 5182$；(ii) $a = 30.4$，$b = 24.9$，$R = 3.2$，$\Delta \vartheta = 20^\circ 32'$，对于间隙宽度为 2 mm 的环，得到 $N = 5106$；(iii) $a = 33.4$，$b = 24.7$，$R = 3.2$，$\Delta \vartheta = 28^\circ 30'$，对于间隙宽度为 3 mm 的环，得到 $N = 5110$；(iv) $a = 35.0$，$b = 24.4$，$R = 3.2$，$\Delta \vartheta = 36^\circ 35'$，对于间隙宽度为 4 mm 的环，得到 $N = 4870$。所有参数值都以球体直径 $d = 1$ 为单位表示。

D.4 新月形物体 我们在数值模拟中使用的新月形物体是由直径为 $d=1$ 的球体组成，这些球体排列在一个规则的立方晶格上。晶格点对应于球体的中心，其选择基于 $xz$ 平面上一条平面参考曲线的参数描述，以及一个在垂直的 $yz$ 平面上生成平滑变化的椭圆截面的规则。图 16 显示了用于构造新月形珠子模型的表面网格。构造过程从 $xz$ 平面上 $y = 0$ 处定义的一条参考曲线 $\varGamma$ 开始，该曲线由参数方程 (D7) 给出：
\begin{align}
x(p) = \eta _x \left (5\cos p - \sin 2p\right ), \qquad z(p) = \eta _z \left (3 \sin p + 2 \cos 2p\right ),
\end{align}
对于 $p \in [0, 2\pi )$，其中 $\eta _x = 0.84 \,s$ 和 $\eta _z = 0.51 \,s$，$s \gt 1$ 表示一个无量纲缩放因子。选择 $\eta _x$ 和 $\eta _z$ 的值是为了匹配用作参考的物理物体的形状，即实验测量值和数值计算值之间的比例 $L/H$ 拼配得相当。对于曲线 $\varGamma$ 上 $x$ 轴上每个整数倍的 $d$，都确定了 $z$ 的最小值和最大值，分别表示为 $z_{\textit{min}}(x)$ 和 $z_{\textit{max}}(x)$。在每个这样的位置 $x$，我们在 $yz$ 平面上定义了一个以 (D8) 为中心的椭圆截面：
\begin{align}
z_m(x) = \tfrac {1}{2}\left (z_{\textit{min}}(x) + z_{\textit{max}}(x)\right ).
\end{align}
椭圆的长半轴设置为 (D9)：
\begin{align}
a(x) = \tfrac {1}{2}\left (z_{\textit{max}}(x) - z_{\textit{min}}(x)\right ),
\end{align}
而短半轴（沿 $y$ 方向）则根据实验观察结果定义为 (D10)：
\begin{align}
b(x) = \varepsilon a(x),
\end{align}
其中 $\varepsilon = 0.6$，即 $\varepsilon \approx W/S$。如果一个晶格点 $(x, y, z)$ 表示一个球体的中心，并且满足条件 (D11)：
\begin{align}
\left (\frac {y}{b(x)}\right )^2 + \left (\frac {z - z_m(x)}{a(x)}\right )^2 \leqslant 1, \qquad \text{对于 } x \in \varGamma.
\end{align}
则该点被包含在物体中。这个条件产生了一个由橙色网格表示的形状（见图 16）。为了引入所研究几何形状的不对称性，应用了两个额外的过滤规则。首先，排除了所有位于 $xz$ 平面上两个圆柱形孔内的晶格点，这些孔的半径为 $r = 0.4\,s$，中心位于 $(x_c, z_c) = (\pm 1.823\,s,\ -0.548\,s)$。图 16 中用红色网格表示了这些孔的表面。其次，对于每个被接受的截面，$y$ 的允许范围进一步受到不等式 (D12) 的限制：
\begin{align}
y_{-}(x, z) \leqslant y \leqslant y_{+}(x, z), \qquad \text{其中} \qquad y_{\mp }(x, z) = \pm \tfrac {7}{20} \varepsilon z \mp \tfrac {19}{40} b(x).
\end{align}
这个条件对应于那些中心位于图 16 中显示的两个绿色平面之间的球体。图 11 中所示的新月形珠子模型由 5513 个球体组成，对应的缩放因子为 $s = 6.2$。

D.5 珠子模型精度的细节 随着珠子数量 $N$ 的增加，与珠子模型相关的数值不确定性通常会降低。对于所有物体，长度 $L/d$ 作为一个离散化参数；它的增加会导致 $N$ 的增加，从而得到越来越精细的模型。我们定义了两个连续分辨率步骤之间的相对差异 $\Delta \mu$，如下所示 (D13)：
\begin{align}
\Delta \mu = \frac {\mu ^{{(i)}} - \mu ^{{(i-1)}}}{\mu ^{{(i)}},
\end{align}
其中 $\mu ^{{(i)}}$ 和 $\mu ^{{(i-1)}}$ 分别对应于当前（更精细）和之前（更粗糙）分辨率下计算出的迁移系数。珠子模型的分辨率直接由参数对 $(N, L/d)$ 表征。作为对本工作中模拟不确定性的一种粗略估计，我们使用了最高可用分辨率与之前较粗糙模型之间的相对差异 $\Delta \mu$。表 8 总结了所有物体的这些值，需要注意的是，在其他间隙的环的情况下，这些数值与间隙为 1 mm 的环的数值相似。新月形和箭头形状的物体观察到的较大偏差可能与它们的复杂几何形状有关。对于所有粒子，在构建相应的珠子聚集体时，很难准确定义它们的表面。

表 8. 珠子模型的收敛性分析。在本工作中使用的最高分辨率模型（精细）与之前的离散化步骤（粗糙）之间的系数 $\mu$ 进行了比较。分辨率由珠子数量 $N$ 和离散化参数长度 $L/d$ 表示。相对差异 $\Delta \mu$ 计算为 $(\mu ^{\textit{fine}} \! - \!\mu ^{\textit{coarse}}) / \mu ^{\textit{fine}}$。这种困难在一个针对间隙为 1 mm 的环的珠子模型收敛性的系统研究中得到了体现，研究表明迁移系数在分辨率范围 $(824, 39)$ 到 $(5182, 67)$ 内表现出衰减振荡。具体来说，$\Delta \mu _{33}$ 的变化范围在 $\pm 3\,\%$ 之内，而 $\Delta \mu _{42}$ 和 $\Delta \mu _{51}$ 的变化范围在 $\pm 15\,\%$ 之内。这样的振荡不允许将结果外推到 $N\rightarrow \infty$ 的情况，因为在珠子数量不断增加的情况下，如果聚集体被一个包含所有珠子内部的固定表面所限制，这种情况通常是适用的。虽然增加珠子的数量理论上可以减轻这些离散化伪影，但计算成本会显著增加。例如，对于箭头形状的物体，将分辨率提高到 $L/d=33$（下一个可能的增量）将产生一个包含 7488 个珠子的模型。由此产生的矩阵运算超出了我们当前计算基础设施（252 GB 的 RAM）的资源限制。

脚注：?这些作者对这项工作做出了平等的贡献。?目前就职于印度古吉拉特邦艾哈迈达巴德大学的工程与应用科学学院。