稳定床面风沙起动的碰撞耗散机制：CFD-DEM耦合模拟与非均匀粒径响应

Collisional Dissipation Mechanism of Aeolian Sand Entrainment over Stable Beds: CFD-DEM Coupled Simulation and Non-Uniform Particle Size Response

摘要: 风沙起动是风蚀与沙尘输送的核心环节，其细观机制尚不明确。文章基于CFD-DEM双向耦合框架，建立可解析颗粒离散运动、相间动量交换及湍流调制作用的数值模型，系统研究稳定床面单颗粒起动与非均匀床面的粒径选择性起动行为。研究发现，稳定床面单颗粒起动并非传统认知的碰撞溅射主导，而是流体动量累积与碰撞耗散相互竞争的动力学过程；碰撞在此主要扮演耗散角色，颗粒成功起动需跨越由碰撞累积决定的临界能量门槛。非均匀床面模拟进一步揭示粒径对起动概率的非单调调控：细颗粒因深陷孔隙、高频碰撞耗散，起动概率反低于中等颗粒；粗颗粒呈现为碰撞驱动、流体辅助的模式；中等颗粒体现由接触主导向流体主导的过渡。

Abstract: Wind-blown sand initiation is the core link of wind erosion and dust transport, yet its mesoscopic mechanism remains unclear. Based on the CFD-DEM two-way coupling framework, this study establishes a numerical model capable of resolving discrete particle motion, interphase momentum exchange, and turbulence modulation, systematically investigating single-particle initiation on stable bed surfaces and size-selective initiation behaviors on non-uniform beds. The study finds that single-particle initiation on stable beds is not dominated by collisional splashing as traditionally believed, but rather a dynamic process involving competition between fluid momentum accumulation and collisional dissipation; collisions here primarily play a dissipative role, and successful particle initiation requires crossing a critical energy threshold determined by cumulative collisions. Non-uniform bed simulations further reveal the non-monotonic regulation of particle size on initiation probability: fine particles exhibit lower initiation probability than medium-sized particles due to deep pore embedding and high-frequency collisional dissipation; coarse particles demonstrate a collision-driven, fluid-assisted mode; while medium-sized particles represent a transition from contact-dominated to fluid-dominated regimes.

1. 引言

风沙运动作为典型的气固两相流，广泛存在于自然环境(如沙漠化、沙尘暴)与工业过程(如气力输送、流化床)中[1]-[4]。因此，对这些输运机制的基本了解对物理学家和工程师来说是重要的，其中一个关键作用往往是确定泥沙运动开始的门槛条件，也称为初始运动。准确理解颗粒在风场作用下的起动机理与运动规律，是预测风蚀、评估输沙通量及优化相关工业设备的基础。由于风沙流动的时空多尺度性、相间强耦合性及颗粒运动的离散性，单纯依靠理论分析或物理实验难以全面、精细地捕捉其内在动力学机制，为应对这一挑战，数值模拟已成为研究气固两相流动不可或缺的重要手段[4]-[6]。

本研究旨在建立一种改进的CFD-DEM双向耦合数值模型，该模型基于欧拉–拉格朗日框架，通过计算流体力学(CFD)求解连续流体相，并利用离散单元法(DEM)显式追踪每个固体颗粒的运动。模型充分考虑了颗粒间的碰撞接触力，以及颗粒体积分数对局部流场的影响，能够更真实地模拟从稳定床面条件下颗粒的蠕移、碰撞到最终跃起的完整动力学过程。本文旨在通过该耦合模型，深入研究均匀与非均匀沙床中颗粒的起动机理，分析粒径、风速及床面结构对起动行为的影响，以揭示风沙起动的细观物理机制。

2. 数值方法

对于CFD-DEM耦合算法，本研究采用了“未解析法”。该方法将流体网格单元尺寸设定为大于颗粒直径，故固体颗粒轨迹无法在单元内完全解析。具体而言，流体流动通过连续介质的局部平均连续性和动量方程来求解，并且使用离散元法直接跟踪单个颗粒的运动。颗粒–流体相互作用力主要包括阻力、压力梯度、流体剪应力张量引起的粘性力。

考虑到DEM的大量计算需求，本文的流场采用了考虑流固耦合的区域平均风速剖面模型。在模拟开始时，在流向方向(Y方向)上采用了对数风速分布，其表示见式(1)：

$u\left(z\right)=\frac{u_{*}}{\kappa }\ln \frac{z-z_0}{z_0}$ (1)

其中u_*表示摩阻风速，k = 0.4是冯卡门常数，z和z₀ = d_p/30分别是风速为零点以上的高度(在计算设置中定义)和沙床粗糙度[2] [7]。

2.1. 流体相控制方程

任何流动问题都满足质量守恒方程，即单位时间内流体微团质量增加的部分，等于流入该微团的质量流量，该定律的数学表达为连续性方程，其表达式见式(2)。类似地，动量守恒定律在流体力学中的数学表述为Navier-Stokes方程，它描述了流体微团的动量变化率等于其所受外力之和，具体表达式见式(3)：

$\frac{\partial \left({\alpha }_f{\rho }_f\right)}{\partial t}+\nabla \cdot \left({\alpha }_f{\rho }_f{u}_f\right)=0$ (2)

$\frac{\partial \left({\alpha }_f{\rho }_f{u}_f\right)}{\partial t}+\nabla \cdot \left({\alpha }_f{\rho }_f{u}_f{u}_f\right)=-{\alpha }_f\nabla P+\nabla \cdot \left({\alpha }_f{\tau }_f\right)+{\alpha }_f{\rho }_{f}g-F_{pf}$ (3)

2.2. 固相控制方程

固相颗粒的运动主要由牛顿第二运动定律求得，颗粒在流场中主要受到重力、流体对颗粒的作用力以及接触力三类作用[8]。颗粒的平动与转动方程如下：

$m_p\frac{d{u}_p}{dt}=m_{p}g+F_p+F_c$ (4)

$I_p\frac{d{\omega }_p}{dt}={\displaystyle \sum T_c}+T_f$ (5)

式中m_p和𝑢_𝑝分别是颗粒的质量(kg)和速度(m∙s⁻¹)，ｇ为重力加速(m∙s⁻¹)，F_c为颗粒之间或者颗粒–壁面之间的接触力(N)；ω_P是颗粒的角速度(r∙ads⁻¹)，Tc以及T_f是接触扭矩和流体扭矩(Nm)，F_p为流体对颗粒的作用为(N)。

颗粒在流场中运动时，会受到流体的阻碍作用，它是影响颗粒在流场中运动主要的力，阻力与许多因素有关，包括颗粒表面的粗糙程度，以及颗粒的旋转速度，颗粒形状，流体是否可压缩以及温度的变化等有关。阻力的计算公式如下[9]：

$F_d=\frac{1}{2}{C}_{d}A{\rho }_f\left|u_f-u_p\right|\left(u_f-u_p\right)$ (6)

$C_d={\left[0.6+\frac{4.8}{\sqrt{{\mathrm{Re}}_p}}\right]}^2{\epsilon }_i{}^{-\xi }$ (7)

$\xi =3.7-0.65{\text{e}}^{\left[-\frac{1}{2}{\left(1.5-\log 10{\mathrm{Re}}_p\right)}^2\right]}$ (8)

式中，A是颗粒的投影面积，C_d是阻力系数，Re是颗粒雷诺数，$\epsilon$ 是孔隙率，$\xi$ 是考虑颗粒周围颗粒的阻碍作用。连续相与离散相间的耦合主要是通过交换动量、质量等完成的，属于双向耦合，示意图如图1所示。

本研究聚焦于高颗粒密度床面条件下的风沙起动机理，此时床面颗粒紧密堆积，相互作用显著，如图2所示。准确解析近壁区流动结构以获取可靠的壁面剪切应力，仍是量化起动阈值的关键。SST k-ω湍流模型继承了标准k-ω模型优异的近壁解析能力，可在y⁺ ≈ 1的条件下直接积分至壁面，无需求助壁面函数便能精确捕捉床面附近的剪切应力分布。此外，SST模型通过在边界层内部切换到k-ω格式、在外区采用k-ε格式，在保持近壁分辨率的同时，增强了对复杂流动分离与混合的预测鲁棒性，并维持了适中的计算成本[10] [11]。综合考虑本研究高颗粒浓度床面所涉及的复杂颗粒–流体–湍流相互作用，以及精确捕捉起动阈值的核心需求，采用SST k-ω湍流模型是兼顾精度、稳定性与计算效率的合适选择。在50 d_p × 10 d_p × 1000 d_p区域内随机自然落下颗粒将此区域填满形成床面，稳定床层的判据是最大颗粒速度V_max必须小1.0e⁻⁴ m∙s⁻¹，床层顶部定义为风速归零位置以下0.46 d_p处[6]。床面前方预留20 d_p的通道使湍流能够充分发展。此外，为避免初始流场不稳定，模拟初始阶段将颗粒固定0.05 s以稳定流场。本文所用的物理参数在表1中呈现。

Figure 1. Schematic diagram of gas-solid two-phase coupling

图1. 气固两相间的耦合示意图

Figure 2. Schematic diagram of the model

图2. 模型示意图

Table 1. Model parameters and constants

表1. 模型中运用的参数和常数

	g	E	${\nu }_2$	e	${\mu }_s$	$\nu$	$\mu$	${\rho }_f$
颗粒密度	重力加速度	弹性模量	泊松比	恢复系数	滑动摩擦系数	空气运动黏度	静摩擦系数	空气密度
2650 (kg∙m^−^3)	9.8 (m∙s^−^2)	1.0 × 10^6 (Pa)	0.3	0.2	0.05	1.5 × 10^−^5 (m^2∙s^−^1)	0.3	1.225 (kg∙m^−^3)

3. 结果讨论

本研究通过CFD-DEM双向耦合数值模拟，成功再现了风沙输运过程中的关键物理现象，包括颗粒的蠕移、跃移以及复杂的颗粒–床面碰撞动力学。模型通过引入耦合对数风廓线并合理设定零风速点高度，构建了贴近真实湍流边界层条件的、具有活跃床面蠕移特征的流动环境。模型验证结果表明，模拟获得的流体起动摩阻风速u_*t在0.19~0.22 m∙s⁻¹之间，平均值为0.20 m∙s⁻¹，与Bagnold、Chepil等人的实验结果高度一致[1] [12]-[15]；同时也在Iversen以及Shao理论预测的200 um颗粒起动阈值范围之内。此外，Fu在模型中引入的湍流脉动效应亦使预测结果与上述区间吻合(如图3所示) [16]。上述一致性充分证实了本模型在湍流条件下颗粒起动模拟方面的准确性与可靠性。

Figure 3. Threshold diagram for 200 um particle entrainment

图3. 200 um颗粒起动阈值图

对颗粒运动下随机选取的颗粒轨迹分析(图4)揭示了与传统认知不同的颗粒起跳模式。模拟结果显示，床面颗粒并非由单一、剧烈的越移颗粒冲击瞬时“溅射”至空中，而是经历了一个复杂的动量逐步累积过程。颗粒在起跳前，其运动轨迹表现为在床面附近长时间的滚动、蠕移并伴随着一系列频繁的颗粒间碰撞；有的颗粒可能会在经历初始的滚动后因能量被床面吸收而停止运动。在颗粒速度演化曲线上，平缓的上升段对应于流体拖曳力与摩擦阻力的持续作用，而尖锐的突变点则标志着有效的碰撞事件。许多颗粒在经历初始滚动后，其能量被床面吸收而停止运动；而另一些颗粒则能在经历数次至数十次有效碰撞后，移动超过自身直径10倍或更长的距离，最终才获得足够动量脱离蠕移层，转变为跃移颗粒。因此，床面颗粒在稳态输送下的起动应该是一个动量积累的过程，而不是由单一跳跃冲击过程以及流体引起的。

图4展示了目标颗粒在流体起动过程中的运动轨迹与其速度演变。从轨迹图可以清晰地看出，颗粒的运动并非一蹴而就，而是经历了明显的阶段性：初始的轻微晃动、沿床面的持续滚动，并最终在轨迹末端脱离床面实现起跳。对应的速度曲线则更为细致地揭示了这一过程的动力学特征：速度并非单调增长，而是呈现出阶梯式上升与瞬时骤降交替的复杂模式。

Figure 4. The relationship between the starting trajectory of a 200 μm particle and its velocity as the particle rolls.

图4. 200 um颗粒起动轨迹以及其颗粒速度随着颗粒滚动变化关系

速度曲线中的瞬时骤降点(如图4所示)与颗粒在轨迹中发生方向突变或速度减缓的位置精确对应，这些被明确标记为碰撞事件。每一次碰撞都导致了颗粒动能的显著损失。这与Wang等人在活跃蠕移层上发生的碰撞是促进颗粒起跳能量的主要来源有本质的不同[6]。在本研究的单颗粒起动场景中，颗粒冲击的是一个静态或准静态的稳定床面，碰撞更像是一个撞击–反弹/停止的过程，大部分动能在此过程中通过摩擦、阻尼和非弹性变形转化为内能而耗散掉。

在图5中对颗粒动量的精细分析揭示出单颗粒流体起动与群体输运状态在碰撞动力学上的本质差异。在单颗粒起动场景中，颗粒与床面固定颗粒之间的碰撞主要起能量耗散作用，而非动量传递与起跳激励。每次碰撞均导致颗粒速度骤降，动量损失显著。这一机制差异清晰界定了稀疏颗粒流与饱和输沙状态下碰撞行为的物理本质。

动量定量分析表明(图5)，在单颗粒起动过程中，流体拖曳累积动量I_flu是颗粒能量增加的正贡献源。而碰撞累积动量的总和为负，其绝对值随碰撞次数增加而增大。这意味着，颗粒必须依靠流体持续做功，克服一系列碰撞带来的能量损失，才能最终积累起足够的动能以脱离床面。因此，在此场景下，碰撞是起动过程的主要阻力机制。

颗粒能否成功起动，取决于流体拖曳力所做的功(I_flu)在扣除了所有碰撞损失(I_dmp)后，是否仍能使颗粒获得大于起跳阈值的净起动。我们的统计发现，成功起动的颗粒其I_flu通常远大于I_imp。这表明，存在一个流体动力必须超越的临界能量耗散门槛。

Figure 5. Momentum distribution of incipiently moving particles

图5. 起动颗粒动量分布图

本章构建了对数正态分布床面以研究天然沙地的起动选择性。首先在计算域中创建了10,000个位置和直径随机的粒子，颗粒直径服从对数正态分布，平均值d_p = 200 μm，标准差σp = 40 μm，颗粒在重力作用下经历自由落体运动，直至形成平静的沉积层。图6显示了模拟生成的床面不同直径颗粒的体积分数分布，细颗粒(135~180 μm)主要填充在粗颗粒间隙中，形成致密的屏蔽层。

Figure 6. Volume ratio of particles of each size fraction to bed surface

图6. 各粒径颗粒体积与床面体积比

如图2所示，在床面上方施加对数风廓线以启动颗粒运动。根据暴露在床面颗粒施加流体阻力。我们将从零速点到一个平均直径d_p的高度的层视为蠕动层，因为这一层内的颗粒主要进行滚动或爬行运动模式，携带或抛射粒子的标准是粒子的底部离开爬行层，本文将蠕动层视为颗粒可能起动的床层，分别统计8种直径颗粒在该起动层的颗粒数。

Figure 7. Number of entrained particles versus friction velocity for different particle diameters

图7. 不同直径颗粒随着摩阻风速变化时颗粒起动数

本文将摩阻风速u_*从0.16 m∙s⁻¹起，以每0.2 s增加0.01 m∙s⁻¹的速率递增。由图7可知，新形成的床面由于并没有经历风选过程，床面比较不稳定：在u_* = 0.16 m∙s⁻¹时的总体颗粒起动数反而比u_* = 0.17 m∙s⁻¹时更高；经过摩阻风速从0.16 m∙s⁻¹到0.18 m∙s⁻¹的风选作用过后，床面趋于稳定，此时方呈现出随着摩阻风速增高而起动颗粒数递增的现象。为精准量化粒径非均匀性对颗粒起动行为的调控机制，本研究选取具有代表性的典型颗粒，追踪其从静止、蠕移至最终离床的完整力学演化过程。图8~10分别呈现了粒径D =135 μm、200 μm及300 μm三类特征颗粒在恒定摩阻风速u_∗ = 0.4条件下的受力–速度–轨迹耦合响应时程，系统揭示了不同尺寸颗粒在相同流场中的动力学异质性。

图8揭示了D = 135 μm细颗粒的起动特征。其轨迹长度约12 d_p，远超越200 um的5 d_p (图8)，且短暂起动过程中经历了高达12次高频碰撞。受力图表明，流体作用力始终占据主导，但由于细颗粒深陷于粗颗粒构成的孔隙网络中，碰撞耗散导致其净动量积累效率显著降低。速度曲线显示，速度在t = 0.01 s后呈现平稳增长，此时才离开蠕移层起动完成。该现象证实，细颗粒虽与流体耦合更紧密，但孔隙约束与高频耗散的双重抑制作用使其起动概率并非随粒径减小而单调递增，挑战了传统风沙理论中粒径越小越易蚀的经典假设。

图9所示的200 μm颗粒因床层不稳定且暴露度更大而率先起动，初始速度也快于更小颗粒；然而，正因其前置位置，前方颗粒的阻碍作用使其初始所受总力反而低于流体力，导致需要约0.01 s的延迟后，流体力才逐渐接近总力并满足起动条件。由图9可见，颗粒经历持续蠕移后最终脱离床面。而受力曲线显示，起动初期流体力显著高于总力，二者在t = 0.055~0.06 s区间内发生交叉后趋于稳定，这一转变反映了颗粒–床面相互作用由接触主导向流体主导过渡；速度曲线则表明颗粒速度在波动中总体上升，碰撞导致明显的瞬时减速，但随后在流体作用下恢复加速，最终经过约5 d_p后完成起跳。

Figure 8. Trajectory, acting forces, fluid forces, and velocity variation of 135 μm particles during incipient rolling motion

图8. 135 um颗粒起动轨迹以及其受到的力、流体力、颗粒速度随着颗粒滚动变化关系

Figure 9. Trajectory, acting forces, fluid forces, and velocity variation of 200 μm particles during incipient rolling motion

图9. 200 um颗粒起动轨迹以及其受到的力、流体力、颗粒速度随着颗粒滚动变化关系

图10中d_p = 300 μm粗颗粒在相同流场条件下，滚动距离超过23 d_p后最终才离地。受力曲线显示，碰撞事件后总力出现骤增，而流体力因颗粒瞬时暴露于气流中也上升。尽管流体力有所增强，但其绝对值仍显著低于碰撞转移的冲击力。颗粒获得的瞬时加速动能中，碰撞提供的能量占比远超流体直接做功的贡献，形成碰撞驱动、流体辅助的模式。这揭示了粗颗粒起动能量主要来源于碰撞解锁机制，流体仅起辅助作用，其起动阈值对床面结构敏感性更强。

Figure 10. Trajectory, acting forces, fluid forces, and velocity variation of 300 μm particles during incipient rolling motion

图10. 300 um颗粒起动轨迹以及其受到的力、流体力、颗粒速度随着颗粒滚动变化关系

4. 结论

本研究基于CFD-DEM双向耦合数值模拟方法，深入探讨了风沙颗粒的起动机理。主要结论如下：

首先，所建立的数值模型具有良好的可靠性。通过引入对数风廓线及湍流–颗粒相互作用模型，模拟成功再现了贴近真实条件的风沙边界层。获得的摩阻风速起动阈值与经典实验及理论预测高度一致，验证了模型在捕捉颗粒起动关键物理过程方面的准确性。

其次，研究揭示了稳定床面上单颗粒起动的“累积–耗散”动力学本质。颗粒起动并非瞬时完成的单一事件，而是一个流体拖曳力持续做功(动量累积)与颗粒–床面碰撞不断耗能相互竞争的复杂过程。在此场景下，碰撞主要扮演能量耗散者的角色，颗粒成功起动的关键在于流体输入的能量能否超越由碰撞决定的临界耗散阈值。

再者，在非均匀床面中，粒径显著影响颗粒的起动机理，挑战了“粒径越小越易起动”的简单认知。细颗粒深陷孔隙网络，受高频碰撞耗散抑制，起动概率并未随粒径减小而单调增加；中等颗粒的起动体现了流体与接触作用的动态过渡；而粗颗粒的起动则更多地依赖于碰撞传递的冲击动量，呈现“碰撞驱动、流体辅助”的模式。

最后，本研究为理解风沙起动的细观机理提供了新视角，但仍有扩展空间。未来工作可聚焦于多因素耦合影响、高浓度稠密相流动的模拟、复杂地形应用以及通过更高精度的实验数据对模型进行细化和验证。