0 引　言

复杂流动中的颗粒运动在自然环境和工业生产等领域中普遍存在，比如云层冰晶以及雨雪的形成^[1]、浮游生物的运动与聚集现象^[2]、纸浆纤维与造纸过程^[3-4]、化学制药过程^[5]等。同时，颗粒的形状通常是各向异性的。但是过去大部分的研究用等效的球形点颗粒模型对颗粒运动进行描述，开展颗粒的沉积、聚集以及与湍流相互作用等问题的研究^[6–10]。然而，颗粒的非球形的形状特征在一定程度上会对颗粒的输运产生影响。与此同时，在一些特定的自然与工程问题中，颗粒的形状特征不可忽略。因此由颗粒形状各向异性引起的取向行为就具有十分重要的研究意义，比如造纸过程中细长的纸纤维的取向影响纸张的力学性能^[3-4]；纤维增强材料中纤维的排列决定了材料的微结构，进而影响材料的力学性能；在减阻控制中，纤维的取向行为影响纤维与湍流之间的相互作用^[11]；浮游生物的取向行为会影响生物的聚集与迁移过程^[12]等。

Fan 和 Ahmadi^[13-14]、Zhang等^[15]较早采用欧拉-拉格朗日方法对杆状颗粒在壁湍流中的运动问题进行数值模拟。随后，Mortensen 等^[16-17]、Marchioli等^[6,18-19]、Marchioli & Soldati ^[20] 和Challabotla等^[21-23]对非球形颗粒在壁湍流中的输运过程开展了大量相关研究工作。其中，欧拉-拉格朗日方法是背景流场在欧拉的观点下进行求解，而颗粒则通过拉格朗日追踪的方法进行求解。壁湍流是一类在壁面附近存在强剪切与丰富的拟序结构的湍流，比如槽道、圆管以及边界层湍流等^[24]。因此，壁湍流中的非球形颗粒两相流研究主要关注以下几类问题：1）非球形颗粒在壁湍流中的统计行为（平动、转动以及取向）；2）非球形颗粒与湍流中的相干结构之间的相互作用；3）非球形颗粒对湍流的减阻控制的作用。

本文主要综述了近些年关于微小非球形颗粒在壁湍流中取向行为的研究工作进展，具体结构如下：第1节介绍相关工作使用的理论与数值模拟方法，第2节讨论非球形颗粒在槽道湍流中的取向行为，第3节介绍颗粒倾向性取向行为的机理，第4节讨论颗粒的取向行为会在近壁区呈现“形状敏感性”的现象，最后对本文的内容进行总结和展望。

1 数值方法

前期的壁湍流颗粒数值研究主要方法为欧拉-拉格朗日耦合的直接数值模拟。这里对流场及颗粒求解方法做简单的介绍。

1.1 流体相

直接数值模拟对不可压缩的纳维-斯托克斯方程进行求解，求解的连续方程与动量方程分别为：

$ \frac{\partial {u}_{i}}{\partial {x}_{i}}=0 $

(1)

$ \frac{\partial {u}_{i}}{\partial t}+{u}_{j}\frac{\partial {u}_{i}}{\partial {x}_{j}}=-\frac{1}{{\rho }_{f}}\frac{\partial p}{\partial {x}_{i}}+\nu \frac{{\partial }^{2}{u}_{i}}{\partial {x}_{j}\partial {x}_{j}} $

(2)

其中 $ {u}_{i} $ 为流体速度，p为压力， $\nu $ 是流体的运动学黏度， $ {\rho }_{f} $ 是流体的密度。本文主要介绍的物理模型为两平板间的槽道湍流，其中流向（x）和展向（y）为均匀的周期方向，且在壁面处满足速度不可滑移条件。在流向与展向均使用伪谱方法求解，而在壁面垂直方向（z）采用二阶有限差分格式。在时间推进上采用显示二阶Adams-Bashforth格式。

主要算例的计算域大多为 $ 12h\times 6h\times 2h $ ，其中h为槽道高度的一半。算例的计算网格数为 $ {192}^{3} $ ，在流向与展向为均匀网格，而壁面法向采用双曲正切函数形式的非均匀网格。算例中基于壁面摩擦速度的摩擦雷诺数 ${Re}_{\tau }=180$ 。在本文中上标“+”代表着变量已被黏性尺度进行无量纲化处理。

1.2 颗粒相

本文主要介绍的颗粒尺寸远小于流体的Kolmogorov尺度，且颗粒相为稀疏悬浮，因此不考虑颗粒对流体的反作用以及颗粒与颗粒之间的碰撞作用。同时，颗粒的尺寸非常小（考虑的微粒弛豫时间与流体弛豫时间小于0.1，但颗粒与流体的密度比可以达到1 000左右），颗粒惯性与流体惯性对颗粒运动学的影响可以被忽略^[25]。此时，颗粒可近似认为跟随流体的轨迹线，即

$ \frac{{\rm{d}}{x}_{p,i}}{{\rm{d}}t}={u}_{p,i} $

(3)

其中 $ {u}_{p,i} $ 为颗粒速度。对于跟随流体的颗粒，颗粒速度就近似等于当地的流场速度。本文主要考虑轴对称椭球颗粒，因此颗粒的取向（轴对称椭球颗粒的回转轴方向）可以通过Jeffery方程^[26]得到，即

$ \frac{{\rm{d}}{p}_{i}}{{\rm{d}}t}={O}_{ij}{p}_{j}+\frac{{{\lambda }}^{2}-1}{{\lambda }^{2}+1}\left({S}_{ij}{p}_{j}-{p}_{j}{S}_{jk}{p}_{k}{p}_{i}\right) $

(4)

其中 $ {p}_{i} $ 为颗粒回转轴的方向（单位向量）； $ {O}_{ij} $ 为流体的旋转张量，即 $\dfrac{1}{2}\Bigg(\dfrac{\partial {u}_{i}}{\partial {x}_{j}}-\dfrac{\partial {u}_{j}}{\partial {x}_{i}}\Bigg)$ ； ${S}_{ij}$ 为流体的变形率张量，即 $\dfrac{1}{2}\Bigg(\dfrac{\partial {u}_{i}}{\partial {x}_{j}}+\dfrac{\partial {u}_{j}}{\partial {x}_{i}}\Bigg)$ 。而 $ \mathrm{\lambda } $ 为颗粒的形状参数，即为颗粒的回转轴与赤道轴的长度之比。有时也使用形状因子 $\varLambda =\dfrac{{\mathrm{\lambda }}^{2}-1}{{\lambda }^{2}+1}$ 对颗粒形状进行描述。不同于颗粒参数 $\lambda$ ，形状因子 $\varLambda$ 的取值范围为 $ [-1,1] $ 。当 ${\lambda } > 1$ 时， $0 < \varLambda <1$ ，此时颗粒为杆状颗粒；当 $0 < {\lambda }<1$ 时， $-1 < \varLambda <0$ ，那么颗粒则为碟状颗粒；当 ${\lambda }=1$ 时，颗粒呈球形，此时 $\varLambda =0$ 。通过颗粒的取向方程可知，颗粒的形状参数 ${\lambda }$ （或形状因子 $\varLambda$ ）主要是调节流体变形率在颗粒取向方程中的占比，进而影响颗粒的取向行为。当形状参数 ${\lambda }\to \mathrm{\infty }$ 或0时（形状因子 $\left|\varLambda \right|=1$ ），此时流体变形率的作用最强。

本文中颗粒相的时间推进格式与流体相的时间推进格式保持一致，而颗粒处的流体信息通过二阶三维拉格朗日插值格式从颗粒附近的流体网格中插值得到。同时，颗粒与壁面之间采用完全弹性碰撞模型。取向随机的颗粒在初始时刻被均匀地以当地流场速度放入充分发展的槽道湍流中，待颗粒在流场中充分发展后进行统计与分析。

2 颗粒在壁面附近的取向行为

在槽道湍流中，因为槽道内部湍流非均匀且各向异性，非球形颗粒的取向与转动行为受颗粒形状的影响较大，尤其对于极细长的杆状颗粒和极扁平的碟状颗粒。而且，颗粒的行为在槽道中的不同位置处会呈现出不同的行为特征。

Challabotla等^[22]针对非球形颗粒在槽道中的取向行为展开了研究。图1展示了瞬时细长颗粒与扁平颗粒在近壁面处的分布，可以明显看到颗粒具有倾向性的取向，同时这种取向分布与颗粒形状有关。图2展示了非球形颗粒回转轴与流向、展向与法向方向的夹角余弦绝对值的统计平均值。在槽道中部，因为中部的流动状态趋近于均匀各向同性湍流^[27]，所以非球形颗粒的运动行为与其在均匀各向同性湍流中的结果基本一致，即杆状颗粒倾向性地朝着涡量方向而碟状颗粒垂直于涡量方向。在惯性参考系中看，颗粒的回转轴方向与惯性参考系各个轴的夹角余弦的绝对值的统计值均趋于0.5，这意味着非球形颗粒在槽道中部的取向在空间上是趋于一种随机分布的状态。然而，从槽道的黏性底层区（ ${z}^{+} < 5$ ）延伸到 ${z}^{+}\approx 30$ 的缓冲区，杆状颗粒与碟状颗粒的取向行为呈现出明显的差异，即细长的杆状颗粒在壁面附近倾向性地朝着流向方向，而碟状颗粒的回转轴方向倾向性地朝着壁面法向。颗粒回转轴的方向与流向或法向方向的夹角余弦统计值随颗粒形状参数从0.01至50的变化而单调变化。该变化趋势与非球形颗粒在各向同性湍流中颗粒与涡量夹角随形状参数 ${\lambda }$ 变化趋势类似。

通过Jeffery方程^[26]可知，非球形颗粒的转动行为与取向行为是一种相互影响的关系，即颗粒取向的时间导数 $\dfrac{\mathrm{d}{p}}{\mathrm{d}{t}}={\omega }\times {p}$ ，其中 $ {\omega } $ 为颗粒的转动角速度。在槽道湍流中，壁面附近的涡量场具有较强的各向异性，并且在展向存在非常强的涡量。图3(a)展示了颗粒与流体微团的展向平均旋转角速度。其结果表明在近壁区非球形颗粒的展向转动受到了抑制，其抑制效果随形状偏离球形的程度的增大而增强。但是随着统计位置逐渐远离壁面 （远离缓冲区），非球形颗粒的转动行为受颗粒形状的影响程度会逐渐减弱。非球形颗粒在壁面的转动受抑制也进一步解释了为什么颗粒在壁面存在较强的倾向性的取向行为。前期工作认为非球形颗粒在近壁区的行为（尤其在黏性底层）与颗粒在线性剪切流中的现象类似，并认为颗粒的转动行为类似Jeffery轨迹^[26]，即杆状颗粒回转轴在流向附近（或碟状颗粒回转轴在法向附近）时需要停留较长的时间才会突然进行翻转。因此，在统计的基础上可以得到非球形颗粒在壁面附近具有倾向性的统计学取向行为。图3(b-d)则进一步展示了颗粒在惯性参考系中各个方向上转动角速度的脉动值随颗粒形状参数的复杂变化规律。与此同时，揭育澄等人^[28-29]也研究了中等雷诺数 ${Re}_{\tau }=1\;000$ 中非球形颗粒的取向行为，并发现颗粒的取向行为的统计结果基本与低雷诺数 ${Re}_{\tau}=180$ 接近^[28-29]。该结果也说明了非球形颗粒在近壁处的取向行为具有一定雷诺数无关的特征^[28-29]。

3 颗粒倾向性取向行为的机理

颗粒取向行为的机理一直以来都是关注的热点问题。在均匀各向同性湍流的研究中，目前普遍接受的观点主要有两种：1）非球形颗粒的取向与流体的涡量以及变形率张量的第二特征值方向相关^[30-32]；2）非球形颗粒的取向与流体的拉格朗日拉伸与压缩方向相关^[33-34]。前者主要是基于欧拉的观点去讨论颗粒与流场物理量之间的关系，其中涡量随时间演化的拉格朗日控制方程与杆状颗粒退化的演化方程相似，唯一的不同点在于涡量方程多出了黏性扩散项^[30,32]。同时，由于涡量与变形率张量之间的关联，颗粒的取向同样与变形率张量的第二特征值方向具有较强的相关性。然而，颗粒与涡量相关的解释仅在各向同性湍流或远离壁面区域的湍流有效，但是在具有较强剪切的壁面附近并不成立。因为壁面的存在，流场在壁面会存在较强的展向涡量，但是杆状颗粒与碟状颗粒均不会倾向性地朝着展向^[22]，所以杆状颗粒与碟状颗粒被观察到垂直于涡量方向。同时，在二维流场中，涡量方向永远垂直于流场的平面。因此，第一种解释具有一定的局限性。对于第二种观点，极细长的杆状颗粒可以看成无限小的物质线段，而该物质线段与流体的拉格朗日拉伸方向渐近一致^[30,35-36]，当时间足够长时，杆状颗粒与拉格朗日拉伸方向一致。Ni 等^[34]发现在均匀各向同性湍流中形状参数 $ \mathrm{\lambda }>10 $ 的杆状颗粒与拉格朗日拉伸方向就具有较好的一致性，而且随形状参数 $ \mathrm{\lambda } $ 的变化并不明显。为了进一步验证其在具有各向异性的湍流中是否也同样适用，赵立豪等^[37]在槽道湍流中分析了非球形颗粒与流体拉格朗日拉伸与压缩方向之间的相关性。

3.1 拉格朗日拉伸与压缩方向

流体的拉格朗日拉伸与压缩方向是指初始时刻为球形的流体微团沿着拉格朗日轨迹线发生拉伸与压缩变形的方向^[34]。因此，为了得到流体的拉格朗日拉伸与压缩方向，首先需要沿着流体迹线积分得到对应的变形梯度张量 ${F}_{ij}$ ，即^[33-34]

$ \frac{{\rm{d}}{F}_{ij}}{{\rm{d}}t}={A}_{jk}{F}_{kj} $

(5)

其中 ${A}_{ij}$ 是流体的速度梯度张量，而 $ {F}_{ij} $ 是流体的变形梯度张量。变形梯度张量 $ {F}_{ij} $ 表征了流体微团沿轨迹线的变形情况。 $ {F}_{ij} $ 不是实对称张量，通常采用柯西格林张量来表征流体微团的变形。本文采用左柯西格林张量进行表征，即^[33-34]

$ {M}_{ij}={F}_{ik}{F}_{jk} $

(6)

通常设定初始时刻的 $ {F}_{ij}={\delta }_{ij} $ （ ${\delta }_{ij}$ 为Kronecker符号），即初始时刻为单位球形张量，随时间积分一段时间后，对 $ {M} $ 进行特征分解。其中M的最大特征值对应的特征方向作为拉格朗日拉伸方向，记为 $ {{e}}_{L1} $ 。同理，最小特征值对应的特征方向为拉格朗日压缩方向，记为 $ {{e}}_{L3} $ 。第二特征值对应的特征方向记为 $ {{e}}_{L2} $ 。而在壁面附近，由于流动近似为线性剪切流，所以在统计结果中拉格朗日的拉伸方向指向流向，而压缩方向指向壁面法向^[29]。

3.2 非球形颗粒与拉格朗日拉伸与压缩方向的相关性

图4展示了颗粒回转轴方向 $ {p} $ 与左柯西格林张量的主轴方向 $ {{e}}_{Li} $ 的相关随积分时间的变化，即 $\langle { {({{e}}_{Li}\cdot {p})}^{2} } \rangle$ 。若相关值为1/3，那么颗粒相对于柯西格林张量的主轴方向是随机分布的；而相关值为1意味着颗粒与对应的主轴方向完全一致。图4横轴为相对初始时刻的积分时长^[37]。图4的结果表明了杆状颗粒逐渐趋近于拉格朗日拉伸方向，而碟状颗粒逐渐趋近于拉格朗日的压缩方向。图5为颗粒回转轴方向 $ {p} $ 与左柯西格林张量的主轴方向 $ {{e}}_{Li} $ 的相关沿空间分布，即 $ \langle { {({{e}}_{Li}\cdot {p})}^{2} } \rangle $ 。图5选取积分时间 ${t}^{+}={t}_{0}^{+}$ +136至 ${t}^{+}={t}_{0}^{+}$ +144进行统计，图片来源于文献[37]。通过颗粒回转轴的方向与柯西格林张量的三个主轴方向的夹角的空间分布图，发现杆状颗粒与拉格朗日拉伸方向以及碟状颗粒与拉格朗日压缩方向在近壁区域存在非常强的相关性。但是，细长的杆状颗粒和扁平的碟状颗粒的取向与拉格朗日的拉伸和压缩方向还是存在比较明显的差异，至少与Ni等^[34]在均匀各向同性湍流中提出来 ${\lambda } > 10$ 的形状条件并不相符。而这种差异产生的原因有待进一步研究。

4 颗粒取向行为的形状敏感性问题

在第3节中，本文介绍到杆状颗粒与拉格朗日拉伸方向以及碟状颗粒与拉格朗日的压缩方向具有较强的相关性。但是，颗粒与拉格朗日拉伸与压缩方向的差异随颗粒形状参数 $ {\lambda } $ 变化依然明显，这与Ni 等^[34]观察到 $ {\lambda }>10 $ 的杆状颗粒在各向同性湍流中的取向就基本与拉格朗日拉伸方向一致的结论存在差别。本章主要讨论在近壁面附近这种差异是如何体现，及其原因。

4.1 同一轨迹线上颗粒与拉格朗日拉伸方向的差异

首先，在崔智文等^[38]的工作中，他们选取一条轨迹线，并截取颗粒已经充分发展且在壁面附近有较长时间停留但与壁面无碰撞过程的轨迹段，如图6所示。图6展示了同一条轨迹线上形状参数 ${\lambda }=23.3$ 的杆状颗粒的取向的时间演化图与拉格朗日拉伸方向的差异。图中 $ \mathrm{\phi} $ 是相对于惯性参考系的方位角， $\mathrm{\theta }$ 是俯仰角。槽道的上壁面在 ${z}^{+}=360$ ，当 ${2\;700 < {t}}^{+}<3\;200$ 时，颗粒处在上壁面的粘性底层区 ${z}^{+} > 355$ ，图片来自文献[38]。颗粒的形状参数 ${\lambda }=23.3$ ，对应于 $\varLambda =0.996\;3$ （图中蓝色虚线），此时杆状颗粒的 $ \varLambda $ 参数已经非常接近于1，其与1的差别仅为 $ \mathrm{\delta }\varLambda =1-\varLambda =0.003\;7 $ 。而颗粒的形状因子 $ \varLambda =1 $ （图中红色实线）意味着颗粒的长细比无限大，并且在充分发展的颗粒场中无限长细比的杆状颗粒的取向与流体的拉格朗日拉伸方向一致。而在差别仅为 $ \mathrm{\delta }\varLambda =0.003\;7 $ 的情况下，人们从直观上认为杆状颗粒的取向应该与拉格朗日的拉伸方向基本一致。然而，图6的结果表明，除了远离壁面的区域，杆状颗粒的取向与拉格朗日拉伸方向在近壁区却呈现出迥然不同的运动行为。在非常靠近壁面的位置，拉格朗日的拉伸方向基本朝着流向方向，但杆状颗粒却持续地翻转。图6中的方位角 $ \mathrm{\phi } $ 是杆状颗粒回转轴在 $ x-z $ 平面的投影与 $ x $ 轴的夹角，而俯仰角 $ {\theta} $ 是杆状颗粒回转轴与 $ x-z $ 平面的夹角。

4.2 杆状颗粒在拉格朗日坐标系中的相对分布

为进一步分析，崔智文等^[38]对细长杆状颗粒在拉格朗日坐标系（由柯西格林张量三个主轴方向组成，即 ${{{e}}}_{L1}$ 、 ${{{e}}}_{L2} $ 和 ${{{e}}}_{L3}$ )的相对取向的概率密度分布进行了研究。图7与图8分别展示了该分布在槽道中部（ ${z}^{+}=180$ ）和槽道黏性底层（ ${z}^{+}=4$ ）中的情况。其中方位角 $ \mathrm{\alpha } $ 是杆状颗粒回转轴在 ${{{e}}}_{L1}-{{{e}}}_{L3}$ 平面的投影与 ${{{e}}}_{L1}$ 的夹角，而俯仰角 $\;{\beta } $ 是杆状颗粒回转轴与 ${{{e}}}_{L1}-{{{e}}}_{L3}$ 平面的夹角。结果表明不管是在槽道中部还是壁面的黏性底层，其欧拉角 $ \mathrm{\alpha } $ 与 $\; {\beta } $ 的概率密度分布函数都存在平台区与幂律区这两个明显的区域。其中平台区的宽度随形状差别 $\mathrm{\delta }\varLambda =1-\varLambda$ （即无限长细比与有限长细比颗粒之间的形状因子的差别）的增加而变宽，而幂律区中的分布函数接近 ${P}\left(\mathrm{\alpha }\right)\sim {\mathrm{\alpha }}^{-2}$ 。崔智文等^[38]认为平台区的宽度反映了杆状颗粒取向与拉格朗日拉伸方向的差异。平台越窄意味着杆状颗粒在绝大部分时候基本与拉格朗日拉伸方向一致。但随着平台的变宽，杆状颗粒取向与拉格朗日拉伸方向的差异逐渐增强。通过图7与图8的对比，不难发现，不管在槽道中部还是在黏性底层，欧拉角 $ \mathrm{\alpha }{\text{与}}\mathrm{\beta } $ 的分布是类似的，但是平台的宽度在壁面附近要明显大于槽道中部的情况。如果将槽道中各层中不同形状颗粒的 $ \mathrm{\alpha } $ 分布的平台的宽度 ${\mathrm{\alpha }}_{c}$ 进行提取（图7与图8中的虚线），并将其与颗粒形状 $ \mathrm{\delta }\varLambda =1-\varLambda $ 进行比较（如图9所示），其结果表明在槽道的各个层的平台宽度 ${\mathrm{\alpha }}_{c}$ 与 $ \mathrm{\delta }\varLambda $ 呈较强的线性相关。图9的结果说明了图7与图8中发现的规律不仅存在于壁面附近存在强剪切的区域，而且也适用于槽道的各个区域。

4.3 机理分析与模型验证

4.1节与4.2节分别讨论了 $ \varLambda $ 接近于1的杆状颗粒在近壁面的区域依然会与拉格朗日的拉伸方向存在明显差别。那么产生这种差别的原因是什么？首先，剪切在当前问题中扮演着重要的角色，但仅只有剪切这一种因素并不能解释为什么 ${\text{当}}\varLambda \to 1{\text{但}\varLambda}\ne 1$ 的杆状颗粒在大部分时候会与拉格朗日拉伸方向还具有较好的一致性。因为在只考虑纯剪切的情况，根据Jeffery 方程^[26]，不难发现 $ \varLambda =1 $ 与 $ \varLambda \ne 1 $ 两者颗粒在行为上的迥异差别。对于前者，剪切流中拉格朗日的拉伸方向（ $\varLambda =1$ ）指向流向方向，而压缩方向会朝着法向。但对于后者，杆状颗粒与碟状颗粒均服从Jeffery轨迹在空间内周期转动^[26]。实际上，对于壁湍流，在壁面附近除了强剪切也存在比较明显的湍流脉动，湍流脉动引入的噪声是否会使得 $\varLambda \to 1{\text{但}\varLambda}\ne 1$ 的颗粒在大部分时候还是会与拉格朗日拉伸方向具有较好的一致性？因此，在崔智文等^[38]的研究中，速度梯度的脉动被引入。他们认为槽道不同区域的平均剪切与速度梯度的湍流脉动之比( $s/D$ )是影响颗粒取向行为的主要因素，其中s是平均剪切，D表征速度梯度的脉动影响强度。当 $s/D\gg 1$ 时，剪切的作用强于脉动，此时非常小的形状差异 $ \mathrm{\delta }\varLambda $ 也会使得杆状颗粒与拉格朗日拉伸方向存在明显的差异；当 $s/D\ll 1$ 时，剪切的作用远小于脉动，此时比较大的形状差异 $ \mathrm{\delta }\varLambda $ 也会使得杆状颗粒与拉格朗日拉伸方向基本上表现出较好的一致性。

为了进一步验证该想法，崔智文等^[38]基于Jeffery方程建立了二维简化的颗粒取向模型，并从方程角度出发分析在存在强剪切与速度梯度脉动下颗粒转动周期随形状的变化，即：

$ \frac{\mathrm{d}\mathrm{\phi }}{\mathrm{d}{t}}=-\frac{s}{2}\left(1-\mathrm{cos}2\phi \right)+{\eta }_{\phi } $

(7)

其中， $ \mathrm{\phi } $ 是颗粒回转轴在 $x-z$ 剪切平面内的投影与x轴方向的夹角，s是平均剪切， $ {\mathrm{\eta }}_{\mathrm{\phi }} $ 是随机脉动。利用文献[39]中的方法，结合Fokker-Planck 随机过程模型，最终可以得到颗粒转动周期统计值的表达式为^[38]：

$\begin{split}& \left\langle \tau \right\rangle s = \\&\qquad \sqrt {\text{π}} {12^{\frac{1}{6}}}{\left( {\frac{s}{D}} \right)^{\frac{1}{3}}}{{\rm{\varLambda }}^{ - \frac{2}{3}}}\mathop \int \nolimits_0^\infty \frac{1}{{\sqrt y }}{{{\rm{e}}}^{ - {{\left( {\frac{3}{2}} \right)}^{\frac{1}{3}}}{{\left( {\frac{s}{D}} \right)}^{\frac{2}{3}}}\delta {\rm{\varLambda }}{{\rm{\varLambda }}^{ - \frac{1}{3}}}y - {y^3}}}{\rm{d}}y\end{split}$

(8)

图10为颗粒平均翻转时间 $ \left\langle { \tau } \right\rangle s $ 与颗粒形状 $\mathrm{\delta }\varLambda =1-\varLambda$ 的关系（ $0 < {{z}}^{+}<10$ ），横轴从左到右代表颗粒形状从杆状逐渐变为球形，图中的符号为直接数值模拟结果，虚线为简化模型的结果，图片来自文献[38]。图10展示的是理论模型与直接数值模拟结果的对比，其中直接数值模拟结果通过统计颗粒在所有时刻均停留在范围 $0 < {{z}}^{+}<10$ 的轨迹线片段上的翻滚周期。理论模型中的s由当前层的平均剪切得到，而D通过计算拉格朗日轨迹线上的速度梯度自相关得到。通过图10可知，理论与数值实验结果的符合度高，进而验证了他们的想法，即平均剪切与速度梯度的脉动的比值影响颗粒与拉格朗日拉伸方向的相关性。同时，通过理论模型可知，当 $ \mathrm{\delta }\varLambda \to 0 $ 时，

$ \left\langle \tau \right\rangle s=\sqrt{{\text{π}} }{2}^{\frac{4}{3}}{3}^{\frac{1}{6}}\mathrm{\varGamma }\left(\frac{7}{6}\right){\left(\frac{s}{D}\right)}^{\frac{1}{3}} $

(9)

当 $\mathrm{\delta }\varLambda \gg {\left(\dfrac{D}{s}\right)}^{\frac{2}{3}}$ 时，颗粒的转动周期趋近于Jeffery理论^[26]，

$ \left\langle \tau \right\rangle s=\sqrt{2}{\text{π}} \delta {\varLambda }^{-\frac{1}{2}} $

(10)

图10展示了颗粒的转动周期随形状变化 $ \mathrm{\delta }\varLambda $ 也存在明显的平台区与幂律区，而平台区的含义代表着颗粒的转动行为基本不随形状 $ \mathrm{\delta }\varLambda $ 改变，幂律区则是颗粒周期符合Jeffery理论。这说明当颗粒形状变化进入到幂律区时，颗粒的行为会产生明显的变化。同时，区分平台区与幂律区的临界值满足 $\mathrm{\delta }{\varLambda }_{\mathrm{c}}\sim {\left(\dfrac{s}{D}\right)}^{-\frac{2}{3}}$ 。当 ${\left(\dfrac{s}{D}\right)} \gg 1$ 时， $ \mathrm{\delta }{\varLambda }_{\mathrm{c}}\to 0 $ ，这意味着非常小的 $ \mathrm{\delta }\varLambda $ 也会使得颗粒与 $ \varLambda =1 $ 的情况产生明显的变化，这与之前分析的结果一致。

5 结　论

本文回顾了非球形颗粒在壁湍流中的取向行为的直接数值模拟研究。在忽略颗粒惯性的作用下，近壁区的杆状颗粒倾向性地朝着流向方向而碟状颗粒则倾向性地朝着壁面法向方向，且颗粒的倾向性取向行为的程度随颗粒偏离球形的程度的增加而增强。而非球形颗粒之所以会在壁面具有倾向性取向，是因为偏离球形程度越大的颗粒与流体的拉格朗日拉伸与压缩方向具有非常强的相关性。同时，研究工作也发现即使非球形颗粒的形状因子 $ \left|\varLambda \right|\to 1 $ ，但是该类颗粒的取向与拉格朗日拉伸或压缩方向( $ \left|\varLambda \right|=1 $ )的在壁面附近还是会存在较大的行为差异。在以杆状颗粒为例的研究中发现，杆状颗粒取向与流体拉格朗日拉伸方向的相关性随颗粒形状的影响规律在壁面和槽道中心位置是类似的，但不同点在于壁面附近较大的平均剪切与速度梯度脉动之比s/D会使得颗粒取向与流体拉格朗日拉伸方向产生明显的差异。基于壁面湍流的特点，建立了以平均剪切与速度梯度脉动之比的颗粒转动周期预测的二维简化模型。模型的建立与验证进一步揭示了平均剪切与速度梯度脉动在颗粒取向趋近于拉格朗日拉伸与压缩方向过程中的重要作用。

在未来的研究工作中，可以考虑非球形颗粒取向行为与湍流中的相干结构之间的相互作用关系，并将其应用在纤维减阻控制的研究之中，进而加深对湍流减阻控制作用的理解。此外，真实流动中非球形颗粒随着形状的不规则性增强，易变形特性会逐渐体现出来，而此时模拟计算中颗粒的柔性需要考虑，所以关于柔性颗粒在湍流中的行为研究也是该领域可以进一步探索的问题之一。

致谢：本文涉及到工作是作者赵立豪在挪威科技大学以及回国后在清华大学的研究内容，全部已公开发表。本工作得到了挪威研究理事会、国家自然科学基金以及清华大学国强研究院项目的支持。