考虑magnus效应的风沙跃移云模型

考虑magnus效应的风沙跃移云模型

在沙流研究中，床面粗糙度作为表面对空气的影响的基本参数参数被高度重视。这种动态参数与边界层外参考风速或摩阻风速之间的关系长期以来受到了许多研究人员的高度重视。Bagnold通过风洞实验研究了气流与沙粒的相互作用,发现所测得的风速廓线在无沙粒运动时遵循壁面定律,即在半对数纸上风速沿高度变化曲线为一直线,若此直线与y轴的交点为(0,y0),则高度y0即为定床面的粗糙度(图1)。依据其所测得的实验结果,Bagnold推断固定沙面的粗糙度为yo=Dp/30,其中Dp为沙粒的粒径;董治宝等人则通过分析实验结果得出了粗糙度不仅是一个几何参数,而且是一个与流场有关的动力学参数的结论。在起沙后,Bagnold发现实测粗糙度总是大于定床粗糙度;Anderson和Owen则认为,跃移沙粒对风沙跃移层外风速分布的作用类似于定床表面粗糙度的影响。在跃移层内风速受沙粒的阻滞作用而远远偏离了对数分布,但在风沙跃移层以上的高度,风速仍服从对数分布,即在半对数纸上风速沿高度变化曲线仍为一直线,将此直线延长后与y轴相交,此交点处的高度yst的值就是起沙后地表的有效粗糙度。由于在进行实验时,实验条件如沙样的粒径分布、风洞中湍流的发育程度等很难保持一致;同时由于跃移层高度区域不大(一般在离地几厘米到几十厘米量级),且在此区域内沙粒的浓度较高,使得测量风的速度较为困难。因而有关起沙后床面粗糙度的规律仍没有在认识上达成一致。本文以Ungar和Haff提出的将风沙流自平衡机制与沙粒-床面碰撞过程相结合的描述稳态风-沙粒相互耦合的风沙跃移云模型为基础,计算了各种摩阻速度下的风速廓线,并采用Anderson等的方法给出了起沙后各种摩阻速度下的有效粗糙度以及Magnus力对有效粗糙度的影响。1模型的构建在分析风沙跃移运动时,应该考虑:①风对沙粒的加速作用;②沙粒对风的阻滞作用;③下降沙粒与地表的碰撞。1.1沙粒maagus动力的动态特性跃移沙粒在运动过程中受到的作用力主要有:风的拖曳力、Magnus升力及重力。风的拖曳力为:fDx=-ρπD2p8(24vDp√(x-u)2-y2+61+√Dp√(x-u)2-y2/v+0.4)(1a)fDy=-ρπD2p8(24vDp√(x-u)2-y2+61+√Dp√(x-u)2-y2/v+0.4)×y√(x-u)2-y2(1b)其中:fDx、fDy分别表示作用在沙粒上的空气的拖曳力在x、y方向上的分量,x、y为沙粒运动速度在x、y方向上的分量,u为气流速度在y方向上的分量,v为空气的运动粘性系数。许多学者在采用频闪光摄影对沙粒运动的动态过程进行研究时发现,沙粒在运动过程中始终在旋转,因此沙粒还受到Magnus升力的作用,此力可表示为:fLx=18πρD3py(θ-12dudy)(2a)fLy=18πρD3p(x-u)(θ-12dudy)(2b)其中:fLx、fLy分别表示作用在沙粒上的Magnus力在x、y方向上的分量;θ为沙粒旋转的角速度。因此的运动轨迹应满足:mpd2xdt2=fx=fDx+fLx=-ρπD2p8(24vDp√(x-u)2-y2+61+√Dp√(x-u)2-y2/v+0.4)×(x-u)√(x-u)2-y2+18πρD3py(θ-12dudy)(3a)mpd2ydt2=fy=fDy-mg+fLy=-ρπD2p8(24vDp√(x-u)2-y2+61+√Dp√(x-u)2-y2/v+0.4)×y√(x-u)2-y2-mpg+18πρD3p(x-u)(θ-12dudy)(3b)Ιp¨θ=Μ(3c)其中:mp为沙粒质量;fx、fy分别表示作用在沙粒上的合力在x、y方向上的分量;Ip为沙粒的转动惯量;¨θ为沙粒的角加速度。当沙粒被视为球体时,其转动惯量为Ιp=110D2pmp(4)而力矩可表示为:Μ=πμD3p(θ-12dudy)(5)将(4)式和(5)式代入(3c)式得mpd2θd2t=10πμDp(θ-12dudy)(6)1.2混合长度理论受沙粒阻滞力作用的空气运动的Navier-Stokes方程为:其中:u为流体速度;ρ为空气密度;t为时间;p为压力;F为每单位体积的气流所受到的体积力;τ为剪应力张量。记μ为空气的粘性系数,则剪应力张量的分量满足如下形式的牛顿流体的本构关系:τij=μ(∂ui∂xi+∂uj∂xi)(8)对于一无限大平坦沙面上流过的二维稳定气流,Ungar和Haff、Werner等人运用混合长度理论将Navier-Stokes方程简化为封闭形式,其x方向的方程为Fx+ddy[ρk2y2(dudy)2]=0(9)其中;k为vonKarman常数;Fx为跃移沙粒对单位体积风的阻力,其表达式:Fx=-s[fDx↑y↑(y)-fDx↓y↓(y)](10)其中:s为单位时间内从单位面积床面上起跳的沙粒数(假定它们以同一速度起跳);y↑及y↓分别为上升沙粒与下降沙粒在高度y处的铅直速度分量。显然,在沙粒跃移层以上,有Fx=0。将fDx代入(10)式,然后再将(10)式代入(9)式,就得出风沙流中气流的运动方程如下:ddy[ρk2y2(dudy)2=-ρπsD2p8[(24vDp√(x↑-u)2-y↑2+61+√Dp√(x↑-u)2-y↑2/v+0.4)×(x↑-u)y↑√(x↑-u)2+y↑-(24vDp√(x↓-u)2-y↓2+61+√Dp√(x↓-u)2-y↓2/v+0.4)×(x↓-u)y↓√(x↓-u)2+y↓](11)1.3沙粒初始速度的计算在近地层中,地面的摩擦阻力使风降速,且摩擦阻力的影响随高度的增加而减小,因此风速随高度而增大。一般来说,在高度不超过10m的大气边界层内,地面的剪切力(拖曳力)或阻力τb可以近似地看作常数。即在10m以内风速与高度的对数值成正比,亦即风速廓线是随高度呈对数分布的。记u*为摩阻速度,定义为(τ/ρ)1/2,摩阻速度只与来流速度有关;记y0为床面粗糙度。当床面为沙子时,y0正比于沙粒的直径Dp,且在一般情形下,y0可取为Dp/30。于是对应于方程(11)的气流运动方程边界条件可按如下方式给出:y=Dp30时u=0(12a)y\?∞时kydudy=u*(12b)考虑Magnus效应时,对应沙粒运动方程的初始条件为:t=0时x=0,y=0,x=0,y=Voθ=0,θ=θ0(13)其中,Vo是沙粒垂直起跳时的初速度。于是方程(3)、方程(6)与方程(11)在边界条件(12)与初值条件(13)下构成了稳态情形风沙跃移运动的定解问题,显然这是一组风与沙粒运动相互耦合的非线性微分方程构成的初边值问题。2沙粒全沙后冲击韧性的求解由上节给出的风沙跃移运动的基本方程(3)、(6)与(11)可知,当给定单位时间内从单位面积的沙床上起跳的沙粒数即床面起沙率s、摩阻速度u*和沙粒铅直起跳的速度Vo后,通过求解上述方程,一般可以获得沙粒跃移轨迹以及来流风场在跃移层内沿高度的分布情况。然而,在实际计算中,由于在稳定状态下床面起沙率s、摩阻速度u*和沙粒铅直起跳的速度Vo等参数相互之间有一定的内在联系,例如:床面起沙率s与摩阻速度u*是直接相关的,且沙粒铅直起跳的速度Vo与沙床表面性质和沙粒反弹或溅起等有关,因而这些参数并不能预先任意选定。当风沙跃移运动达到稳定状态时,一般可以认为从沙床进入跃移层的沙粒数将与从跃移层落到沙床的沙粒数相同,对此可用击溅函数来描述。最简单的击溅函数是由Ungar和Haff(1987)提出的,即δ-函数形式。为了简化计算,这里我们通过采用这一最简单的击溅函数来给出s、Vo和u*等之间的关系。其δ-击溅函数为:S(V0,Vi)=AV2iDpgδ3(V0-BV3iDpgey)(14)这里S表示被以入射速度Vi下落的一颗沙粒溅起后,以速度V0起跳的沙粒数;ey为y方向的单位矢量;A、B为无量纲常数,它们一般由床面的粗糙度、紧实程度等床面特性决定。显然,(14)式的物理意义为如果一颗沙粒以速度Vi撞击床面,将会有AVi2/(Dpg)颗沙粒以同一起跳速度Vo=BV3i/(Dpg)ey溅起。在跃移运动达到稳定状态时,应有S(Vo,Vi)=1。于是Vim=DpgA(15)其次,为了给出方程(11)中的床面起沙率参数s,我们将边界条件(12)用下式替代y=Dp30时u=0,dudy=C(16)这里c为一选定且可调整的常数。通过对方程(3)、(6)在初始条件(13)下和方程(11)在边值条件(12)下的求解,反算出与给定沙粒起跳速度Vo和床面起沙率s相对应的满足稳定条件(15)的摩阻速度u*。具体计算步骤如下:(1)给出输入参数Vo,s,ν,κ,Dp,mp,ρa的值,计算初始风速分布u0(y)、A以及沙粒的入射速度Vim。(2)选取常数c>0的值,将u=u0(y)代入方程(3)和方程(6)后在初始条件(13)下对其求解以得到沙粒跃移轨迹,并将结果记为x1(t)和y1(t)。(3)将x1(t)(或xi(t))和y1(t)(或yi(t))代入方程(11)后,在边值条件(12)下对其求解以得到新的风速廓线,并将结果记为ui(y)。(4)将ui(y)代入方程(3)和方程(6)后,在初始条件(13)下对其求解以得到沙粒跃移轨迹新的迭代值xi+1(t)和yi+1(t)。(5)用xi+1(t)和yi+1(t)分别代替xi(t)和yi(t),然后进行步骤3的计算过程以得到风速廓线新的迭代值ui+1(y)。重复步骤(3)～(5)的计算过程直到第i+1次迭代算出的ui+1(y)与上一次迭代算出的ui(y)对任意高度y满足:|ui+1(y)-ui(y)|<ε1(17)为止。(6)由xi+1(t)和yi+1(t),可计算出沙粒入射速度Vim的新的值V^im。如果V^im≠Vim,则要通过调整常数Vim并重复步骤(3)～(6)的计算过程直到:|V^im-DpgA|<ε2(18)为止。一般地,V^im大于(小于)由方程(18)得到的Vim,则选取比上一次小(大)的c值。最后,在得到了风速廓线后摩阻速度的值可通过求解方程(12b)得到。同时,人们在研究土壤风蚀及风沙运动时所关注的其他物理量如:单宽输沙率、有效粗糙度以及风速廓线等可以很容易地得到。显然,所有计算结果主要依赖于参数s的值,而s又与摩阻速度密切相关,通过改变s的值我们可以得到包括摩阻速度在内的人们所关注的物理量。因此,通过数值模拟,我们可以得到与摩阻速度相对应的一些宏观物理量。3参数对有效粗糙度的影响图2给出了有效粗糙度随摩阻风速变化曲线。图中各数据点表示用本文模型计算得到的结果,实线为用最小二乘法对实验结果进行拟合后的曲线,从图中可以看出拟合曲线与数值模拟结果吻合得很好。其拟合曲线可用如下方程表示:y0st=2.01622×10-4-0.00207u*+0.00401u*2(19)即有效粗糙度为摩阻速度的二次函数,这与Owen的推测结果定性相似。Owen首先假定运动沙粒对风沙跃移层以上风速的影响类似于床面粗糙度的作用,并且他还进一步假定运动沙粒对风速的影响与风沙跃移层的厚度成正比。在上述假定的基础上,他通过一系列简化和理论推导后指出:起沙后地表的有效粗糙度与摩阻流速u*密切相关,并且它们之间存在着如下关系y0st=cu*2/g(20)比较(19)式与(20)式可以看出,Owen的结论与本文计算结果存在着一定的差异,这说明Owen的简化假设使得其提出的有效粗糙度与摩阻风速之间的关系只是一种简单的定性描述性的关系。yast=-6.2