非均匀风作用下的沙粒跃移运动研究

非均匀风作用下的沙粒跃移运动研究

1沙粒结构及流场分析跃移运动是沙害物理的中心主题。本研究包括沙害层中气流速度、运动状态、输沙量等物理因素的空间分布，以及单宽输沙量的测量和计算。目前的研究方法主要有两种———实验观测和数值模拟。许多学者在这一领域,特别是数值模拟上做出了开拓性的工作。Owen［１］提出了沙粒跃移轨迹的一般形式;White等［２］研究了马格努斯力对沙粒跃移轨迹的影响;Ungar等［３］模拟了沙粒与空气之间的耦合作用,得出了跃移达到稳定之后的风速廓线和输沙通量分布;Zheng等［４］研究了电场力对输沙率的影响。现有研究绝大多数是针对平坦床面上的沙粒跃移运动,很少有涉足倾斜床面的。Iverson等［５］研究了倾斜床面对沙粒起跳的影响。White等［６］对沙丘坡面上的跃移进行了初步的探讨,给出了倾斜床面上沙粒运动控制方程,研究了沙丘表面颗粒粒径分布情况。蒋红等［７］研究了沙丘迎风坡面沙粒跃移及跃移对风场的反馈作用。Sauermann等［８］研究了新月型沙丘迎风坡面上单宽输沙率分布情况。李建国等［９］研究了不同角度和不同高度防沙堤周围的流场分布。李得禄等［１０］研究了丘间低地不同位置处的风速廓线和输沙率分布。李双权等［１１］研究了沙质风蚀坑中起伏地形与风场之间的相互作用。在已有的沙丘跃移运动的数值模拟中,初始风场一般根据普朗特边界层理论用对数廓线近似表示,这与实际沙丘周围的风场分布存在较大的差别。本文结合风工程和空气动力学的研究成果［１２－１５］,提出了风沙流层中平均风速分布的计算公式,在此基础上,利用颗粒轨道模型,模拟了倾斜坡面上颗粒跃移运动状态,统计分析了输沙浓度的空间分布以及在单位时间内通过单位宽度面积所输运的沙量,即单宽输沙率。所得计算结果与已有的研究成果吻合良好。其中当地起沙率根据风洞实验［１６］所得的拟合公式估算。在保证精度的前提下,本文的计算简单易行,这为风沙流工程应用提供了一个可行的模型和方法。2理论模型本节分为两个部分,第一部分是沙丘周围的流场分布,第二部分是坡面沙粒的跃移计算。2.1垂直风速的连续方程对于倾斜坡面风场的研究,通常用一无量纲参数加速比ΔS来定量描述:式中:u(z)为倾斜坡面以上z高度处的风速;u０(z)为平坦地面以上z高度处的风速。大量学者针对山地的加速效应进行了风洞试验和数值模拟研究,其中以Lemelin等［１４］提出的算法(简称LSD模型)最为典型,可以用于定量描述三维对称小山、二维沙丘、山坡和路堤周围的风场。其平均风速的加速比ΔS(x,z)计算公式为式中:H为山体最大高度;x为某一点与坡脚点之间的水平距离,在坡脚处x=0;L为特征长度,是坡顶点到H/2高度处的水平距离;z为倾斜坡面上的局地高度。Kaimal等［１３］和Taylor等［１５］研究指出,在二维山体结构中最大加速比ΔSｍａｘ为在本文计算中,利用(4)式进行流场计算,这样可以充分近似模拟沙粒对风场的反馈作用。如图1所示,根据LSD模型,当给定对应平坦地面上风速的水平分量u０(z)和速度加速比ΔS,则倾斜地面上任一点(x,z)处水平风速u(x,z)为式中:z′=z-xtanθ,θ为坡面倾斜角。为了得到垂直风速,必须应用连续性方程:于是,对于倾斜坡面,有由式(5)、式(6)和式(7)可知,在沙丘坡面上风速方向平行于坡面方向。由图1,地面坐标系oxz、局地坐标系o′x′z′、坡面坐标系ox１z１之间的关系为假定参考水平地面上的风速分布由对数廓线给出,即式中:κ=0.41,为VonKarman常数。则空间任一点(x,z)处的风速分布为式中:z′=z-xtanθ在坡面坐标系中,任一点(x１,z１)的风速分布为进一步,可由下式求得坡面上摩阻风速分布:2.2空气动力学参数re仅考虑空气拖曳力和重力,则颗粒运动方程可写成式中:ρ为空气密度;ρｐ为沙粒密度;Dｐ为颗粒直径;为颗粒相对运动速度;为阻力系数;Re=ρVrDp/μ为颗粒雷诺数,μ为空气动粘性系数;g为重力加速度,θ为沙丘坡面倾斜角。假设沙粒垂直于坡面起跳,且起跳速度采用Anderson等［１７］提出的Gamma分布:当坡面上起沙率分布函数N(x１)已知时,各断面处的单宽输沙率可由下式算出:3单宽输沙率与数值结果比较为了探讨本文所提算法的可行性,我们对风洞实验［７］与野外观测［８］的风沙单宽输沙率等进行了数值模拟。算例1:选取迎风坡倾角θ=10°,与风洞实验［７］中所选取的沙丘迎风坡倾角相同,沙丘高度H=4m,则L=11.3m,ΔSｍａｘ=0.705,上游参考摩阻风速u＊=0.5m·s－１,z０=0.5mm。计算时,在迎风坡面上从坡脚到坡顶均匀选取13个断面进行跃移计算。此外,坡脚、坡中和坡顶3处的起沙率［７］,分别为2.152×10６、3.209×10６和6.037×10６,其他各坡面处的起沙率,由线性插值得出。图2给出了倾斜坡面上单宽输沙率的实验测量值与数值模拟结果的比较。从图2可以看出,从坡脚到坡顶,随着坡面位置的升高,各截面处的单宽输沙率逐渐增加,且变化率也逐渐加大。这是由于气流的质量守恒,沿坡面向上,流体不断加速,同时地表摩阻风速也逐渐变大,引起地表起沙率的增加,进而提高了风沙流单宽输沙率。另外,本文的数值模拟结果与蒋红等［７］的数值模拟值吻合良好。但与实验结果比较,在坡顶处相差较大,这可能是在实际流动中,在坡顶处产生分离和涡流,降低了气流的输沙能力所致。图3给出了坡脚、坡中和坡顶3处的输沙浓度沿高度的分布,依次标为A、B、C。从图3发现,数值模拟与实验测量结果吻合良好,本文算法可以定量表征倾斜坡面不同位置处输沙浓度沿高度的变化。在给定截面上,随着高度的增加,单位面积的输沙率逐渐减小,大致呈负指数分布,与文献中的报道相一致［４，７］。算例2:为了与文献比较,取沙丘高度H=34m,长度L=84m,则ΔSｍａｘ=0.81,迎风坡倾角θ=11.44。上游摩阻风速为u＊=0.36m·s－１,z０=0.25mm。这与文献的研究条件相一致。从坡脚到坡顶均匀选取17个截面进行计算。由于文献中未给出起沙率的数据,本文参考风洞实验［１６］,获得如下的线性拟合关系式首先利用式(13)、(14)可得各截面处的摩阻风速,之后根据式(19)可得沙丘迎风坡面各截面处的起沙率,然后再利用起沙率数据进行坡面不同位置处的单宽输沙率的计算。由图4可知,本文的计算结果与Sauermann等［８］的野外观测相一致。在坡脚和坡顶处偏大,在坡中位置偏小。这一方面,如同算例1的情形,实际流动在坡脚和坡顶可能有分离或回流产生,同时自然界沙丘顶部一般是曲面,是三维流动,而文中平面假设只是其近似;另一方面,可能是由(19)式所得的起沙率和实际起沙率存在一定差别,导致偏差较大。4不同风速条件下从动态风速分布特征本文给出了沙丘迎风面颗粒跃移运动的一种较为简捷的算法。模拟结果表明:沙丘迎风坡面上的沙粒跃移情况较平坦床面要复杂许多,起伏地形对风场的影响较大,风速的变化使得不同位置处的沙粒起跳率和初速度都不尽相同;从坡脚到坡顶,随着高度的增大,坡面摩阻风速、起沙率、单宽输沙率等物理量也将随之增加。比较模拟结果与已有文献可得,总体变化趋势一致,但是在某些位置上存在偏差。分析原因可能有:沙丘坡面性质复