沙漠运动

模拟风沙运动的离散颗粒动力学模型

孙其诚，王光谦

（清华大学水利水电工程系)

摘要：根据风沙运动的特点，建立了基于单沙粒动力学的离散颗粒动力学模型。依据该模型，模拟了9 000个沙粒在风力作用下由静止不动到充分发展的稳定态的全过程，细致研究了沙床处于稳定态时床面形态、跃移轨迹以及跃移质/蠕移质间交换规律，所得结果有助于理解风沙运动的内部机理，这是目前常用的模型所做不到的。研究表明：离散颗粒动力学模型能够真实的模拟风沙运动现象，研究范围涉及风沙运动规律和床面形态的发生与发展等关键的问题，因而比其它模型更具优势。

关键词：风沙运动；风沙流结构；模型；计算机模拟

基金项目：国家杰出青年科学基金(59525914)；中国博士后科学基金资助项目。

作者简介：孙其诚（1970-），男，清华大学水利系博士后，主要从事气固两相流的实验研究和数值计算。

收稿日期：2000-09-06

近几年春季，我国华北地区沙尘天气频频发生, 并涉及部分南方省市, 对民航交通及日常生活带来严重影响, 直接经济损失数以亿计, 因此加强风沙运动规律的研究已经迫在眉睫。国内外的众多学者对此进行了大量的室内外实验研究和理论上的探索[1，2]，但是由于风沙运动过程的内在复杂性和影响因素的多样性，至今仍有许多方面没有搞清楚，比如沙粒的起动机制等。我国则更是存在应用技术多、基础研究少、过程研究少的突出问题。

20世纪80年代以来，越来越多的研究者开始对风沙运动过程进行数值模拟，提出的模型形式各异。研究发现，几乎所有的模型首先假设风沙运动处于稳定态，所涉及到的空气含沙量、跃移沙粒的速度分布函数和轨迹、跃移沙粒与

蠕移沙粒的能量、动量交换、蠕移沙粒的溅射系数等反映风沙运动内在机理的及其关键的因素无不是通过随机采样产生[3，4]。这些模型只注重最终的统计结果，忽略了中间过程中每一时刻沙粒的位置、速度、能量耗散等非常重要的信息。事实上，风沙运动过程受外界因素的影响总是处于非稳定态，那么模型中采用的参数就是一种不精确的平均，模拟结果相对粗糙。目前，风沙运动的研究方法逐渐倾向于微观分析，并与宏观综合相结合；研究深度也不仅限于对现象的描述，而是更着重于揭示内在原因及其演化历程，并力图尽可能准确的预测其未来的发展。因而，现有的模型对于全面反映风沙运动的非均匀特征和动态变化规律具有一定的局限性，不适合风沙运动过程基础研究发展的需要。

随着CPU速度的提高和高效率计算方法的开发，从20世纪90年代中期，从单颗粒的受力分析入手，并跟踪颗粒轨迹的离散颗粒模型[5，6](discrete element model; granular dynamics model)首先在气固流态化、颗粒技术等领域发展起来，应用范围越来越广，为此国际著名期刊《Powder Technology》出版了介绍离散颗粒模型最新研究进展的专辑(2000年第4期)。在风沙物理学领域，在对沙粒的基本规律有了深刻认识的基础上，我们就可以采用越来越精细的离散颗粒模型，越来越准确的模拟风沙运动的全过程。本文建立了适用于风沙运动本身特点的离散颗粒动力学模型，模拟了9 000个沙粒在风力作用下从静止到充分发展的全过程，对推移质的含沙量分布和输沙量、跃移质与蠕移质的交换以及床面形态的发生与发展等关键问题做了详尽的研究，所得结果与已知观测结果符合较好。

1 离散颗粒动力学模型

沙粒在运动过程中不断地与其它沙粒发生碰撞，由于碰撞时动量的交换瞬时完成，其冲力远大于重力、气流曳力和沙粒间的摩擦力，因此在这一瞬间仅考虑碰撞引起的速度的变化；在两次碰撞之间，每个沙粒在曳力和重力的作用下运动。根据沙粒在上述不同过程的受力特点，其运动分解为受冲力支配的瞬时碰撞过程和受曳力和重力控制的非瞬时漂移过程，从而建立了对沙粒-沙粒以及沙粒-气流两种作用分别处理的离散颗粒动力学模型（细节见文献7，8）。

1.1 沙粒区域分布与近地层流场

任意堆积的沙粒的空间位置都可以分为图1所示的3个区域：下部的沙粒密实排布，处于两虚线之间的沙粒则排布疏松（厚度为一个沙粒粒径），而上方的沙粒漂浮在空中。空间位置的不同对应着本质不同的受力情况：上部的沙粒主要受到重力和气流曳力的作用，称之为近地气流层；中部的沙粒则涉及到起动问题，需考虑沙粒间的摩擦力、重力和曳力，称之为地表层；下部的沙粒由于地表层的掩盖，不受气流的曳力，称之为屏蔽层。

野外观测和风洞实验结果表明，在发生风沙运动以前，近地层风速遵循对数分布[1，9]

||=u0lny/y0; y0=D/30

(1)

其中y, y

分别为床面高程和粗糙度，D为沙粒粒径，均以mm计。在本模型中，按照式(1)给定近地层流场，并假定流线与床面平行，如图1实线所示。

1.2 沙粒的起动

地表层中沙粒的受力情况及作用力表达式已经清楚，但是不同的研究者采用的起动条件不尽相同[10]。本文采用的起动条件为

F D=μ·mg

(2)

其中F

D 为沙粒受到的曳力：F

D

=1/8πD2C

D

ρ

g

||,μ是沙粒的摩擦系数。

图2 沙粒碰撞示意图

图1 沙粒分层与风速分布示意图

Fig.1 Layer configuration and assumption

of local wind velocity

Fig.2 Collision diagram of sand particles

1.3 沙粒运动分解

1.3.1 碰撞过程

假设沙粒是球形的准刚体，不考虑沙粒的变形；碰撞为两体碰撞，且发生在质心平面上。颗粒的碰撞首先假设为两体碰撞，碰撞点仅为两颗粒的接触点。

如图2所示, 在发生碰撞的两沙粒质心平面建立碰撞坐标系, 原点选在碰撞点

上, y轴和x轴分别沿法向和切向。碰撞过程中产生瞬时法向冲量和切向冲量,

两沙粒动量和角动量守恒。设碰撞结束后平均法向冲量和切向冲量分别为P

y , P

x

,

则速度和角速度改变可表示为下列方程

沙粒(3)

沙粒(4)

其中 m

a ,m

b

是沙粒a,b的质量;(υ

x,a0

,υ

y,a0

,ω

a0

,(υ

x,b0

,υ

y,b0

,ω

b0

)分别为沙

粒a,b碰前的速度和角速度;(υ

x,a ,υ

y,a

,ω

a

),(υ

x,b

,υ

y,b

,ω

b

)分别为沙粒a,b碰

后的速度和角速度;I

a ,I

b

分别为沙粒a,b的转动惯量,I

a

=1/2m

a

R2

a

,I

b

=1/2m

b

R2

b

。

Wang和Mason[11]把颗粒的碰撞划分为5种模式,得到了相应的P

y 和P

x

的

分析解。本文采用的二维圆片形颗粒仅涉及其中的两种: 滑动碰撞(sliding collision,物理意义为接触点的切向相对速度始终不为零), 粘滞碰撞(sticking collision, 物理意义为摩擦力相对很大, 以至于刚刚碰撞, 接触点