近床面区挟沙气流流速垂线分布特征（倪晋仁，李振山）
倪晋仁1,2，李振山2,3
(1.北京大学 环境科学中心, 2. 水沙科学教育部重点实验室; 3. 清华大学 水利工程系)
摘要：以测试参数较为齐全的挟沙气流风洞实验资料为基础，结合已有各家资料重点对近床区挟沙气流流速垂线分布进行了系统的分析。在对近床面区沙粒跃移层内气流剪切应力分布的理论分析基础上，指出了该区的剪切应力由沙粒临界起动风速、沙粒垂直速度的垂线变化和输沙量垂线分布三者决定，并进一步导出了近床面区的挟沙气流流速垂线分布廓线方程，由该方程计算得到的结果与实测资料吻合很好。此外，研究还表明，Bagnold 提出的不同风力条件下挟沙气流流速廓线聚焦于一点的观点基本正确，但这一点不宜作为跃移层内外风速廓线的分界点，实际的分界点高度是沙粒跃移层外气流摩阻速度平方的0.015倍；Owen 提出的跃移层外风速廓线方程基本合理，但其中有效粗糙度公式中的系数C取值为 0.04 时更符合实测资料；来流速度相同的情况下，挟沙气流的摩阻速度大于清洁气流的摩阻速度。
关键词：风速廓线; 近床面区; 能量损耗; 摩阻速度; 风洞实验
作者简介：倪晋仁（1962-），男，山西山阴人，北京大学教授。
1 引言
挟沙气流流速垂线分布因受运动沙粒的影响而不同于遵循对数规律的清洁气流的流速垂线分布。早在1941年，Bagnold［１］就通过风洞实验发现挟沙气流的速度廓线在半对数坐标中不再完全是一条直线，且不同风力下的气流速度廓线在离床面 0.3cm处重新聚焦于一点，并指出在来流流速一致的情况下，挟沙气流的摩阻速度与清洁气流的摩阻速度相同。然而，后来不少学者的实验结果与此有所不同［２～6］。Owen［７］从理论上对挟沙气流流速分布进行分析后认为，沙粒跃移层以上气流速度垂线分布符合半对数规律，但是他以沙粒平均跃移层的高度作为动力粗糙度，这也与后来其他学者关于粗糙度取值变化很大的研究结论不符［8，9］。关于沙粒跃移层内的风速分布，Owen 认为可用“颗粒剪切应力”的概念进行分析，但至今这方面的研究较少［10，11］。
通过分析各家实验条件及测试数据后发现，挟沙气流实测结果的分歧和矛盾主要源于各家实验数据往往只反映了挟沙气流流速分布的某一侧面，具体表现为实验的风速相对较小、近床面测点较少、实验时铺沙床面长度太短、气流中沙粒浓度相对较低，所获资料受测量精度所限未能充分反映出沙粒运动对气流速度分布的影响。此外，只测量挟沙气流风速分布，而无对应的清洁气流风速分布资料，就无法有效显示二者详细区别之点。实测资料的代表性不强导致了理论分析中的参数确定较为粗糙，如上下层风速廓线的分界高度不明确、上层风速廓线中粗糙度取值变动大和沙粒跃移层内气流剪切应力分布难以确定等。本文针对上述问题，从实验和理论两方面重点对近床面区挟沙气流流速分布进行较为深入的探讨。
2 风洞实验方法
为了充分揭示挟沙气流的速度分布特征，本实验在详细比较了以往多家风洞实验的优缺点后，针对已有风洞实验的不足，在实验中增大了风速范围，增加了铺沙床面长度，采用了两种不同级配的沙粒，尤其是增加了近床面区的风速测点，并同时测量了挟沙气流流速分布和相应的清洁气流流速分布。
实验是在中国科学院兰州沙漠研究所沙坡头沙漠科学实验研究站室内、外两用风洞中进行的。该风洞为直流闭口吹气式活动风洞，实验段宽1.2m、高1.2m、长21m、扩散段长3m。实验时距实验段前缘1m处开始铺沙，直至风洞出口处，沙床宽1.2m，厚6cm，长达23m，这样就确保了风沙流的充分发展。实验时床面刮平后都由水平尺进行了校准。实验分A组和B组，A组实验沙粒较细，沙源为腾格里沙漠东南缘沙坡头沙丘沙，平均粒径为0.17mm；B组实验沙粒较粗，沙源为同一地区孟家湾附近沙丘沙，平均粒径为0.35mm。
风速测量采用防沙毕托管连接数字压力仪采集，并通过微机自动记录。防沙毕托管由合金管制成，并经过了专门的率定，它的内径为2.4mm，外径为3mm，其动压孔迎向气流，静压孔背向气流。数字压力仪为国产BPY-I型数字压力仪，测量误差小于0.15%。通过预备实验中观测床面蚀积情况发现，风沙流在距实验段前缘10m之内就已经得到了充分发展。因而实验确定的风速测量位置应位于10m之后，具体设置在距实验段前缘13m的断面中央。风速测点布置时特意在近床面安排了多个测点。沿垂线共测量11个点，距床面高度分别为0.25、0.5、0.75、1、2、3.5、5、10、20、35和60cm。其中前六个测点通过防沙毕托管提升实现，后五个测点各置一个毕托管，共设6个防沙毕托管。采样时间间隔为1秒钟，每一测点采集约15个数值，然后取平均值，11个测点测风时间共需5分钟左右。以往各家风洞实验的风洞轴线风速多低于15m.s-1，不能有效揭示大风速范围内沙粒输移规律的状况。本实验在确定沙粒临界起沙风后，又设定了五档风速，这五档风速以风洞实验段前缘风洞轴线风速（简称进口风速）为标准，8.5、11.5、13.5、16.5和22.5m.s-1。将床面喷湿，然后分别以这五个风速吹风，便得这五档风速对应的清洁气流。清洁气流流速的测量位置及测量方法与挟沙气流的情形相同。
实验步骤为：首先测量沙粒临界起沙风速，接着测量五档风速下的挟沙气流流速分布，然后喷湿床面，重复这五个风速，测量清洁气流流速。每次实验后都重新布置床面，以确保各实验中床面的一致性。
3 风洞实验结果
3.1 风速分布特点
两组实验的气流流速测量结果见图1。气流的摩阻速度是通过对风速测点值进行对数曲线拟合求得。对于清洁气流，所有测点都用于摩阻速度的计算；对于挟沙气流，只计算上段廓线的摩阻速度，用于计算的测点限制在3cm 以上。
| 由图1可见：（1）在各种风速条件下清洁气流的速度分布廓线在半对数图上为一条直线, 说明流速分布都较好地遵循半对数规律；（2）挟沙气流中，风速廓线上段在半对数图上仍为一条直线，说明继续遵循半对数规律，廓线下段呈逐渐下凹趋势，偏离半对数规律，并且随着风速的增大，偏离点有上升趋势，证实了Bagnold［１］的实验结果；（3）不同风力下的挟沙气流速度廓线确实存在焦点，其高度随沙粒粒径的增大而升高（A 组实验中，风速廓线的焦点约在1cm 左右；B组实验中，廓线焦点在 2cm 左右），焦点的高度远大于Bag［１］提出的 0.3cm，也不符合 Zingg［１３］提出的焦点高度是沙粒粒径10倍的函数关系；（4）焦点以下风速值变化情况与沙粒粒径大小密切相关，沙粒较细时（A组实验），风速基本不随主流速度改变而变化，这 |
Fig.1 Wind velocity profiles over beds with and without sand movement fortwo test sands |
| 与Gerety 和 Slingerland［５］所得的结果相同；沙粒较粗时（B组实验），风速随主流速度的变化有减小趋势，这又与 Bagnold［１］的结果相同，说明以往实验结果各反映了挟沙气流流速分布的一种情况；（5）由于挟沙气流在近床面区的风速分布开始偏离半对数规律的范围随主流速度的增大而增大，因此在风速较小时该偏离点高度小于焦点的高度；在风速较大时，偏离点高度往往又大于焦点的高度；就多数情况而言，将焦点作为风速分布开始偏离半对数规律分界点会产生很大的误差，也带来许多不便，因而它不宜作为挟沙气流流速廓线上下段之间的分界点；（6）来流速度相同的情况下，挟沙气流流速分布廓线上段与所对应的清洁气流流速分布廓线并不象 Bagnold［１］实验所指出的那样相互平行，而是前者的流速梯度大于后者的流速梯度，这表明挟沙气流上段的摩阻速度要大于清洁气流的摩阻速度。这一结果与 Neuman 和 Nickling［６］的实验结果相似。实际上，挟沙气流的流速随着高度的增加会逐渐与所对应的清洁气流流速趋向一致，表现在风速分布廓线上就是挟沙气流的流速廓线与清洁气流的流速廓线随着高度的增加逐渐靠近直至最后相交；（7）沙粒较细时，挟沙气流的流速都大于同高度的临界起沙风速；沙粒较粗时，挟沙气流的下部流速却都小于同高度的临界起沙风速；说明沙粒大小对近床面区风速大小有强烈影响。 |
3.2 气流能量损耗及损耗率特征
与清洁气流相比，挟沙气流由于运动沙粒的影响会损耗一部分能量，这部分能量就是气流用于输移沙粒的能量，也可近似地认为是沙粒从气流中所获得的能量。气流在单位体积内损耗能量（ΔE）可由同一高度清洁气流流速（U′）和挟沙气流流速（U）的平方差与空气密度的乘积表示，即
ΔE＝ρ(U′2－U2) | (1) |
η＝ρ(U′2-U2)/ρU′2 | (2) |
气流能量的损耗率分布见图3。可以看出：（1）气流能量的损耗率随高度的增加呈单调减小趋势，并在近床面区取得极大值，这一点与气流的能量损耗分布不尽相同；（2）来流速度不大时，能量损耗率分布也近似遵循半对数规律，分布曲线在半对数坐标图中呈直线；（3）随着来流速度的增大，上层气流能量损耗率分布仍然近似遵循半对数规律，下层气流的能量损耗率分布廓线呈下凹状，已偏离半对数规律，且偏离点随来流速度增大逐渐升高；（4）随着来流速度的加大，气流能量的损耗率也增大。
图2 气流能量损耗垂线分布 |
图3 气流能量损耗率垂线分布 |
4 跃移层外气流流速分布方程
Owen 提出了跃移层外气流流速分布方程，即
U/U＝1/κln(z/z′0) | (3) |
式中z′0(＝CU2*2g)为有效粗糙度，C为系数，取0.02。但是，后来不同研究者［８，９］提出的系数C值的变化范围较大，约在0.019～0.16之间。考虑到Owen得出C=0.02时所依赖的是Bagnold［１］、Chepil［１４］和Zingg［１３］三家的实验数据，其中Chepil 实验时的风洞长度仅为3.6m，Zingg 实验时只在实验段中央 46cm 宽的底板上铺上了沙粒（风洞实际宽 90cm ），挟沙气流或者处于不饱和状态，或者受侧向二次流影响严重［１５］，因而由这两组数据来反映挟沙气流流速分布的特点可能不具有太强的代表性。由于C值的大小与挟沙气流的饱和状态关系密切，所以在确定系数C的值时就要充分考虑风沙流的饱和状态。通常情况下，当风沙流流经的距离大于10倍的风洞高度时才能得到充分发展，达到饱和状态。我们从以往风洞实验数据中，选择了符合这一条件的实验数据，并示于图4。这些数据来源于Bagnold［１］、Horikawa 和 Shen［１６］及 Neuman 和 Aljaars［１７］。可以看出大部分数据点落在了（3）式所表示的风速廓线左上侧，表明系数C应大于0.02。
|
图5 本文实验资料与（3）式的比较 |
τ＝τa＋τp | (4) |
大量的研究［１０，１１，１８］表明，沙粒跃移层内气流的速度分布可表示为
| (5) |
上式表明挟沙气流流速U的变化规律取决于τa 的变化规律，由（4）式又知，在给定τ的情况下，τa 的大小完全依赖于τp 的大小。这样，研究U的变化规律问题最终归结为寻求τp 的变化规律的问题。刘大有等［１９］从剪切应力的定义出发，探讨了τp 的表达形式，指出沙粒的剪切应力正比于向上颗粒的垂直质量通量以及向下颗粒与向上颗粒的水平平均速度之差，即
| (6) |
式中 m、n、w、u分别表示单一沙粒质量、向上或向下运动的沙粒数、沙粒的垂直速度和沙粒的水平速度，下标↑和↓分别代表向上运动和向下运动，字母上的“－”表示平均。将（6）式转换为相对形式，即
| (7) |
下标0表示床面处各参数。
以（7）式为基础，可以进一步探讨沙粒剪切应力的具体表达式，（7）式右边第一项可表示为
| (8) |
式中 和 分别为床面和高度 z 处的沙粒水平平均速度。鉴于平衡输沙状态下向上运动和向下运动的沙粒数相等，所以2mn↑＝q，2mn↑0＝q0,q 为某一高度的输沙量，q0为床面处输沙量。因此 |
| (9) |
研究表明［２０，２１］，输沙量的垂线分布按指数规律递减，根据Nalpanis等［２２］的研究，输沙量垂线分布函数可用下式表达
q/q0＝e-zH | (10) |
式中H为近床面区沙粒跃层有效厚度, 即跃移层内外风速廓线的分界高度。将（10）式代入（9）式得
| (11) |
如果忽略沙粒垂直运动方向上的空气阻力，根据能量守恒定律有
| (12) |
由（12）式得
| (13) |
将（11）式和（13）式代入（7）式并整理得
| (14) |
近似地认为＝0.5(
↓)，＝0.5(↑0+↓0)，所以（14）式右边第三项可变换为
| (15) |
沙粒起跃进入气流中，由于气流速度大于沙粒运动水平速度，因而沙粒在气流中是逐渐加速的。尽管当下降沙粒接近床面时，沙粒的速度可能会超过同高度的气流速度，进而此时沙粒运动速度会有所减小，但从沙粒整个运动过程看，下降沙粒在接近床面时速度减小的过程显然处于次要地位，与沙粒加速过程的重要性不属于同一数量级，可以忽略。假设沙粒在气流中的水平平均加速度为a，即可得到
| (16) |
| (17) |
| (18) |
t↑+ t↑＝t | (19) |
式中t↑和C分别表示上升颗粒和下降颗粒在气流中的运动时间，t为沙粒在气流中运动的总时间。将（16）至（19）式代入（15）得
| (20) |
根据刘贤万［２３］及 McEwan 和 Willetts［２４］的研究，t可表示为
| (21) |
β为系数。同理，上升沙粒和下降沙粒在气流中运动的时间分别为
| (22) |
| (23) |
将（21）式至（23）式代入（20）式，并结合（13）式，不难得出
| (24) |
由于床面处的气流剪切应力等于沙粒的临界剪切应力，所以
τp0＝τ-τa0＝ρU2*-ρU2*t | (25) |
将 (24) 式和（25）式代入 (14)式，并整理得
τp=ρU2*(1－r2)（1－z/H）e-z/H | (26) |
式中r＝U*t/U*o由（4）式得
τa＝ρU2*[1-(1－r2)(1－z/H)e-z/H] | (27) |
上式即为近床面区跃移层内气流剪切应力方程。可见该剪切应力由沙粒临界起动风速、沙粒垂直速度的垂线变化和输沙量垂线分布三者所决定。将（27）式代入（5）式便得
| (28) |
上式即为近床面区跃移层内气流速度分布的微分形式。上式中根式内的各项参数形式简单、意义明确。等号右边根式内第一项（1－r2）表现了沙粒临界起动风速对气流速度分布的影响；第二项(1-z/H)反映了沙粒跃移过程中沙粒垂直速度变化对气流速度分布的影响，第三项e-z/H*表达了输沙量垂直分布对气流速度的影响。对（28）式进行数值计算，并将计算结果与本文实测资料相比较发现H应满足下式
2gH/U2*＝0.29 | (29) |
此时U＝5U*（见图6）。上式可以变换为
H＝U2/6.7g≈0.015U2* | (30) |
将(30) 式代入 (28) 式，且将（28）式中的r还原为 U*t与U*之比，并将系数C＝0.04代入（3）就得到了有关挟沙气流流速分布的方程，即
| (31) |
上式与实测数据吻合很好（图6）。
| 6 结论 (1)近床面区沙粒跃移层内的气流剪切应力由沙粒临界起动风速、沙粒垂直速度的垂线变化和输沙量垂线分布三者共同决定，由此导出了沙粒跃移层内气流速度廓线方程，该方程与实测资料吻合很好。 (2)通过对精测的风洞实验资料分析，发现 Bagnold 提出的不同风力条件下挟沙气流流速廓线交于一点的观点基本正确，但这一点不宜作为跃移层内外风速廓线的分界点，实际的分界点高度是沙粒跃移层外气流摩阻速度平方的0.015倍；来流速度相同的情况下，挟沙气流的摩阻速度大于清洁气流的摩阻速度。 (3)由沙粒运动引起的气流能量损耗随距床面的高度之增加开始时逐渐增大，到达某一高度达到极大值后又逐渐减小。气流能量损耗及损耗的极大值高度都随来流速度的变大而增加，基本上也随沙粒粒径的增大而增加。气流能量的损耗率随高度增加呈单调减小趋势，并随来流速度的加大而增大。 |
图6 挟沙气流流速廓线方程（(31)式与多家实验资料的比较 |
(4)Owen 提出的跃移层外风速廓线方程中的系数C=0.02 所依赖的多数实测挟沙气流流速分布数据不具代表性，本文实验数据及多家典型实测数共同表明C的取值为 0.04。
致谢：本项研究中的风洞实验得到了中国科学院兰州沙漠研究所刘贤万研究员、王国昌工程师和沙坡头沙漠研究站及王涛研究员、肖洪浪研究员、高鹏恂同志的大力支持。此外，本文还得到了北京大学吴世亮、曲轶众和钱征寒同志的帮助。作者在此一并致谢。
参考文献［1］Bagnold R A. The physics of blown sand and desert dunes. London: Methuen, 1941. 57～65.［2］河村龙马。飞砂の的研究。东京大学理工学研究所报告。1951, 5 (3/4):95-112.
［3］Kadib A A. A function of sand movement by wind. University of California, Berkeley, Tech. Rept. HEL-2-12. 1965. 68 - 69.
［4］White B R, Mounla H. An experimental study of Froude number effect on wind～tunnel saltation. Acta Mechanica, 1991， Suppl.1: 145 -157.
［5］Gerety K M, Slingerland R. Nature of the saltating population inwind tunnel experiments with heterogeneous size-density sands. In: Brookfield M E, Ahlbrandt T S, (eds). Eolian Sediments and Processes. Developments in Sedimentology, 38. Amsterdam: Elsevier, 1983. 115-131.
［6］Neuman C M, Nickling W G. Momentum extraction with saltation: implications for experimental evaluation of wind profile parameters. Boundary-Layer Meteorology, 1994, 68: 35-50.
［7］Owen, P R. Saltation of uniform grains in air. Journal of Fluid Mechanics, 1964, 20: 225-242.
［8］Sherman D J. An equilibrium relationship for shear velocity and apparent roughness length in aeolian saltation. Geomorphology, 1992, 5: 419-431.
［9］McEwan I K, Willetts B B. Adaptation of the near-surface wind to the development of sand transport. Journal of Fluid Mechanics, 1993, 252: 99-115.
［10］Raupach M R. Saltation layers, vegetation canopies and roughness lengths. Acta Mechanica 1991, Suppl. 1: 83-96.
［11］McEwan I K. Bagnold's kink: a physical feature of a wind velocity profile modified by blown sand? Earth Surface Processes and Landforms, 1993, 18: 145-156.
［12］McEwan I K, Willetts B B. Numerical model of the saltation cloud. Acta Mechanica, 1991 Suppl. 1: 53-66.
［13］Zingg A W. Wind tunnel studies of the movement of sedimentary material.Proc. 5th Hydraulics Conf., Iowa University Studies in Engineering Bulletin , 1953, 34: 111-135.
［14］Chepil W S. Dynamics of wind erosion. Soil Science 1945, 60: 397-411.
［15］Gerety K M. Problems with determination of U* from wind-velocity profiles measured in experiments with saltation. In: Barndorff-Nielsen O E, Moller JT, Willetts B B (eds.) Proceedings of International Workshop on the Physics of Blown Sand. Department of Theoretical Statistics University of Aarhus, 1985. 271-300.
［16］Horikawa K, Shen H W. Sand movement by wind action (on the characteristics of sand traps). U S Beach Erosion Board, Tech. Mem, No. 119, 1960. 51.
［17］Neuman M C, Aljaars M. Wind tunnel measurement of boundary-layer response to sediment transports. Boundary-Layer Meteorology, 1997, 84: 67-83.
［18］Anderson R S, Haff P K. Wind modification and bed response during saltation of sand in air. Acta Mechanic, 1991, Suppl. 1: 21-52.
［19］刘大有, 董飞, 贺大良。风沙二相流运动特点的分析。地理学报, 1996, 51(5): 434-444.
［20］吴正。风沙地貌学。北京: 科学出版社, 1987. 61.
［21］Livingstone I, Warren A. Aeolian Geomorphology: An Introduction. Harlow: Addison Wesley Longman, 1996. 22.
［22］Nalpanis P, Hunt J C R, Barrett C F. Saltating particles over flat beds. Journal of Fluid Mechanics, 1993, 251: 661-685.
［23］刘贤万。实验风沙物理与风沙工程学。北京: 科学出版社, 1995. 31.
［24］McEwan I K, Willetts B B. On the prediction of bed-load sand transportrate in air. Sedimentology, 1994, 41: 1241-1251.
