欧拉颗粒流模型下的近床面风沙跃移速度分布特征

刘 芳, 金阿芳, 热依汗古丽·木沙

(新疆大学机械工程学院, 乌鲁木齐 830047)

风沙颗粒在沙床表面的运动是一类复杂的多相流过程,随着风速的增加,粒径小于500 μm的沙粒最先被风吹起。当气流对沙层表面施加的剪应力超过阈值0.05 N/m2时,沙粒就会发生跃移运动,相较于沙粒的蠕移和悬移运动,跃移运动所需风速最小,且是风沙输移中沙量的主要来源,因此沙粒的跃移运动是风沙流研究中最重要的物理运动状态。

从风沙边界层动力学的研究角度来看,在大气边界层中,风沙流是固体颗粒与气流的相互运动,能在大气边界层中形成明显受风沙运动干扰的反馈两相耦合作用的贴地气流层。在风沙流中,由于受到气相裹挟和湍动、沙粒运动和碰撞等因素的影响,风沙流动行为具有高度的随机性、无规律性和结构不稳定性,即风沙流在跃移阶段的运动具备稠密颗粒气固两相耦合的所有非线性系统特征。邹学勇等[1]利用随机概率方法推导出起跃沙粒垂直初速度的分布函数,进一步解释了风沙流结构呈对数分布。亢力强等[2]采用相位多普勒粒子分析仪(phase Doppler particle analyzer，PDPA) 测量了风沙两相流动中沙床面上沙粒起跳和碰撞速度的概率分布以及不同高度处沙粒速度的概率分布,实验结果表明沙床面上沙粒起跳和碰撞速度概率分布均呈现对数正态分布规律, 沙粒垂直速度概率分布在不同高度处都可用正态分布函数来描述,碰撞和起跳角度均可表示为指数分布函数。肖锋军等[3]同样使用PDPA测量了近床面沙粒的跃移速度,并给出了沙粒冲击速度和角度的联合概率分布函数。程旭[4]在对沙粒的法向速度、水平速度和起动角度的分布研究发现,各曲线均表现出正态分布的特征,且沙粒起动速度的频率分布与风速无关。何丽红等[5]对已有的沙粒起跳垂直速度分布函数进行了非参数检验,最后提出了跃移沙粒起跳垂直速度分布函数服从Weibull分布的结论。罗生虎等[6]对风沙流中沙粒的起跳初速度进行分析发现该速度分布函数服从瑞利分布。金阿芳等[7]运用光滑粒子流体动力学法(smooth particle hydrodynamics，SPH)方法模拟了平坦沙床上沙粒的起跳至回落过程,并对随机抽取的1 000个沙粒进行了统计分析后发现,跃移沙粒起跳速度和起跳角的概率密度函数分别服从对数正态分布和指数分布。

为了完整考虑颗粒相的湍流输运过程,并考虑现实环境中计算机硬件及CPU速度的限制,运用欧拉-欧拉两相颗粒流模型对风沙流运动进行模拟。由于欧拉模型涉及的参数众多,眼下有关欧拉-欧拉模型对风沙流的模拟研究还不充分,张伟等[8]利用欧拉双流体模型研究了该模型中风速和碰撞恢复系数对风沙运动的影响。在此基础上,运用Fluent软件建立欧拉两相流中的颗粒流模型,对风沙流的起动、发展及衰退过程进行模拟。沙粒的速度分布作为研究输沙量分布的重要参数,主要包含水平和铅垂两个方向的运动信息。对风沙流整个运动流程中的物理数据进行统计分析,从沙相速度分布特征中研究沙粒跃移运动中的普遍规律。

1 风沙流数学模型

1.1 近地表层风沙运动的控制方程

在大气边界层里,风沙运动最活跃的区域是下部低于2 m的近地表层,受地表条件的影响,贴地层的速度梯度大并伴有明显的湍流特性,流场受下垫面影响大。风沙活动取决于地表沙粒和大气之间的动力和热力作用,在非极端情况的天气条件下,热量因素与动力因素相比影响较小,可忽略不计。用欧拉-欧拉模型模拟风沙流运动时,气体和沙粒相具有各自的质量和动量平衡输运方程。各相质量平衡方程为

(1)

式(1)中:αk为体积分数；ρk为密度；Uk为k相的平均速度，下标k为g表示气体,k为p表示颗粒相。

各相的动量平衡方程为

(2)

式(2)中:Pg为平均气体压力;gi为重力加速度;∑k,ij为有效应力张量;Ik,i为不考虑平均气体压力影响的相间平均动量传递速率。当只考虑相之间的拖曳力和平均粒子弛豫时间标度τp时,Ik,i为

(3)

式(3)中:

(4)

(5)

(6)

式中：CD为阻力系数,由于沙粒起动不仅受气流平均速度的影响,湍流效应也会影响沙粒的扩散,引入颗粒雷诺数Rep来说明沙粒对气流的影响；dp为颗粒直径；vg为气体分子运动黏度；局部瞬时相对速度vr,i等于局部粒子速度up,i与瞬时气体速度ug,i之差；Vr,i为平均颗粒流相对速度,Vr,i=Up,i-Ug,i。颗粒应力张量∑p,ij定义为

(7)

式(7)中:碰撞压力Pp和体积黏度λp决定了粒子间碰撞时的能量损失。

1.2 两相的体积分数方程

多相流各相间的相互贯穿和融合用相体积分数的概念来描述,这里沙相和气相体积分数分别表示为αp和αg,体积分数能表示各相所占的空间大小,其中每一相都分别满足质量守恒定律和动量守恒定律。多相流中k相的体积Vk定义为

(8)

1.3 湍流模型

由于气流表观雷诺数的范围为15 000～50 000,且底部沙粒形成的粗糙表面会导致气流速度呈明显的湍流分布。这进一步说明了在近地面高度的沙面流场属于湍流边界层。现运用标准κ-ε湍流模型对方程进行封闭。

2 数值计算

计算区域的设置如图1所示。根据风沙流跃移运动中的参数变化特点,将计算区域设置为1 m×1 m的二维矩形区域,网格为均匀四边形网格,数量为111 556个。沙粒密度ρp=2 650 kg/m3,空气密度ρg=1.255 kg/m3,沙床厚度为0.05 m。

图1 数值计算平坦沙床示意图Fig.1 Schematic diagram of numerical calculation of flat sand bed

2.1 欧拉颗粒流模型设置

在欧拉模型中,将沙相设置为颗粒相,颗粒黏度选用syamlal-obrien模型,颗粒体积黏度选用lun-et-al模型,最大堆积率为0.63,曳力模型设置为gidaspow。

2.2 沙相设置

在矩形计算域中标记出0.05 m×1 m的沙相区域,沙相区域网格数为5 678个,流场初始化后利用Patch功能将沙相区域的体积分数设为1。

2.3 边界条件

算例中边界条件的设置为左边界为空气相净风速度入口,下边界无滑移壁面边界,为使风场发展充分,且符合自然风条件,模型上边界设置为对称边界,右边界设置为压力出口边界。

2.4 速度入口边界条件

3 计算结果

在粒径为0.000 1 m球形沙粒组成的平坦静止沙床上进行模拟,摩阻风速为0.265 m/s,该风速下,不同时刻在沙床中心沿高度y分布的沙相体积分数αp，如图2所示。由图2可见,除沙床及其表面附近区域外,从沙粒起动到稳定的不同时刻αp在lnαp-h的半对数坐标系中均呈直线分布,且αp在10-3～10-7范围内的数据与实验结果[9-10]吻合良好,这说明该模型能有效模拟风沙流的跃移运动,且相比于实验结果其数据范围更为广泛,即该模型有一定的预测能力。在0.3 m高度处(距初始床面0.25 m),沙相浓度出现饱和,其中饱和高度的αp=10-7,随着时间的推移,沙相在气流中的浓度逐渐保持稳定,而在饱和层上方,沙相体积分数衰减加快。不同时刻沙相体积分数云图如图3所示。

图2 不同时刻沙相体积分数随高度变化Fig.2 The volume fraction of sand phase changes with height at different times

图3 不同时刻沙相体积分数云图Fig.3 Cloud maps of sand phase volume fraction at different times

沙粒的速度分布是研究输沙量分布的重要参数,与拉格朗日法不同,欧拉方法不能得到单颗沙粒的运动速度,但通过对网格内沙相速度的平均就可以得到不同时刻沙相的速度信息,在此涉及的沙粒速度均为网格节点处沙相的平均速度。在对全流程不同高度处沙粒垂向速度进行统计发现,沙粒垂向速度在不同高度处均呈正态分布,该结论与文献[2]结果相一致。由图4可以看出, 沙粒垂向速度分布在不同高度均存在高度集中区域,且中心随高度向右偏移,随着高度的增加,沙粒垂向速度分布越发松散,这是由于沙粒浓度及沙粒与气流速度两相发生耦合引起的反馈。

图4 全流程不同高度处沙粒垂向速度分布Fig.4 Vertical velocity distribution of sand grains at different heights in the whole process

现有研究的沙粒速度分布是通过统计若干颗沙粒的速度信息后得到的,统计结果受沙粒数目影响大,且常用的实验方法在统计沙粒垂向速度时,往往会遗漏掉速度为0的沙粒,这样得到的沙粒速度分布是不准确的。根据不同时刻沙粒体积分数αp的变化,对αp≥10-7的沙粒速度信息进行统计,得到不同时刻沙床附近所有沙粒在跃移运动中的速度分布,如图5所示。

由图5可知，不同时刻αp≥10-7区域的沙相水平速度呈对数正态分布。进一步分析得到,沙相水平速度在沿流场方向均是正值,且速度分布除0时刻外整体向左偏移,分布并不对称,因此沙相水平速度并不符合普通正态分布,而用对数正态分布拟合较好。即lnvx～N(μ,σ2),其中沙相水平速度vx为取值为正的连续随机变量,vx的概率密度为

f(vx,μ,σ)=

(9)

式(9)中:μ、σ2分别为原始数据的均值和方差,它们的值与风沙颗粒的位置、形状及尺度参数相关。由图5可知,在不同时刻沙粒的水平速度分布并不相同,随时间推移,vx的均值变小,且整体分布趋于平缓。

图5 不同时刻沙相水平速度分布Fig.5 Horizontal velocity distribution of sand phase at different times

图6所示为在不同时刻,αp≥10-7区域的沙相垂向速度vy分布。图6中的分布曲线为拉普拉斯分布(双指数分布)。可以看出,沙粒的上升段(vy>0)和下降段(vy<0)速度都服从指数分布,在沙粒初始运动阶段,起跳沙粒数量逐渐增多,随后沙粒开始降落,参与运动的沙粒数目逐渐增大,最终保持稳定的跳跃状态,但沙粒的上升和下降速度分布并不对称(分布向下降段偏移),尽管沙粒对床面的冲击和碰撞会激发更多沙粒参与运动,但沙粒对风速的阻碍会导致风速减小,若没有风力和沙源的补充,可以预测沙面最终能趋于平静,即再无沙粒运动。沙相垂向速度分布可以表示为vy～La(μ,b),vy的概率密度为

图6 不同时刻沙相垂向速度分布Fig.6 Vertical velocity distribution of sand phase at different times

(10)

式(10)中：μ为位置参数；b>0为尺度参数。

4 结论

采用二维欧拉颗粒两相流模型对风沙的流动过程进行模拟,实现了风沙流沙相体积分数及跃移层沙粒速度分布的研究,得到以下结论。

(1)风沙流从起动到稳定的过程中沙相体积分数在lnαp-h的半对数坐标系中均呈直线分布,该模型不仅能有效模拟风沙流的跃移运动,还具有一定的预测能力。沙相浓度在0.3 m高度处出现饱和,其中饱和高度的沙相体积分数为10-7,在饱和层上方,沙相体积分数衰减加快。

(2)通过对体积分数大于等于10-7的沙粒速度信息进行统计,得到不同时刻沙相水平速度呈对数正态分布、沙相垂向速度分布为拉普拉斯分布(双指数分布),并给出了相应的速度分布概率密度函数。

(3)尽管风沙流是一种多耦合、非线性的随机运动,但利用欧拉颗粒流模型及数理统计的方法能够有效地揭示风沙颗粒流体的动力学运动特征。由于没有考虑混合粒径及沙粒形状对速度分布的影响,这些因素会影响概率密度函数中参数的确定,这也将是后期研究的一个方向。