典型株型沙生灌丛对风沙流场影响的数值模拟

闫 晴，李菊艳，尹忠东，刘金苗，柳宏才

（1.北京林业大学水土保持学院，北京 100083；2.新疆维吾尔自治区水土保持生态环境监测总站，新疆 乌鲁木齐 830002）

植被对地表起到保护作用，增加植被覆盖率是防风固沙的主要手段［1］。然而，在干旱、半干旱地区，由于水分条件的限制，植被主要以低盖度沙旱生灌木为主，许多灌木以单株或单丛的形式存在，其与风蚀关系的研究可简化为单植株的防风固沙研究［2］。但不同类型植株对周围气流场的影响存在差异［3］。因此，研究不同株型单株植被周围的风速与输沙特征，对干旱、半干旱区合理选择风蚀防治中的植被类型具有重要意义，同时也对建设不同形态的稀疏灌丛防护林工程具有一定的参考意义。

株型会显著影响植株的防风固沙效益，目前有关植株形态与风蚀关系的研究主要存在两种方法，一类是采用传统的风洞试验或野外试验探究植株形态对风蚀影响的差异，如亢力强等［4］通过风洞试验测量了细长状植株和上大下小形状植株对空气动力学粗糙度的影响；Abbas等［5］对两种株型植株进行研究，指出在不同风速条件下，女贞（Ligustrum lu⁃cidum）对强风的抵抗力大于秋英（Cosmos bipinna⁃tus）。此类研究主要探究了单株植物对周围风场特征的影响，对于沙粒的运动及植株周围输沙特征关注较少。另一类是利用数值模拟探究单一植被要素周围的流场变化，相比野外试验或风洞试验，数值模拟在数据获取及工作量等方面具有明显优势。李正农等［6］将豆瓣黄杨（Buxus sinica）简化为3D对称模型，有效模拟了树木对于风场的影响；Liu等［7］通过CFD软件研究了不同形状单株植物周围的流场，指出底部大而顶部小的植物具有更好的防风效率。以上研究增强了人们对植株形态影响单株植被防风效果的认识，但在现有单一植被要素数值模拟中，更侧重探究乔木等高大植株周围的风场变化，对于灌丛及不同株型灌丛周围风速与输沙特征的研究较少。

沙生灌丛在沙质荒漠生态系统中占有重要地位，维持了干旱、半干旱地区生态环境的稳定［8］。沙生灌丛根系发达，生命力顽强，在生长季枝叶密度较大，区域植被盖度相应增加，可有效降低风动量，从而起到防风阻沙作用。植被形态和结构是控制风蚀的优先参数［9-11］，不同株型灌丛的防风阻沙效益存在显著差异。本文基于数值模拟，以3 种典型株型的荒漠灌丛为研究对象，通过分析不同类型灌丛前后风速、空气动力学粗糙度、沙粒体积分数等变化规律，探讨不同株型植被的防风阻沙效益，对比得出优势灌丛株型，并利用已有风洞试验进行验证，旨在揭示不同植株形态对风沙流活动的影响存在差异，促进干旱、半干旱地区防风固沙工程设计中对灌丛株型的重视，为防风固沙工程及生态恢复工程中合理选择不同株型的防风沙植被提供参考。

1 数据与方法

1.1 研究对象

本文研究区位于乌兰布和沙漠东北部，区域内生长着不同株型的荒漠沙生灌丛，以坛形、梭形、帚形居多，如白刺（Nitraria sphaerocarpa）、柠条（Cara⁃gana Korshinskii）、霸王（Sarcozygium Bunge）、杨柴（Hedysarum mongolicum）、沙拐枣（Calligonum mon⁃golicum）、梭梭（Haloxylon ammodendron）、花棒（He⁃dysarum scoparium）等［12］，对比选取白刺、梭梭以及沙拐枣作为3种典型株型参照。模拟所需植被特征参数源于已有野外调查试验［13］，调查于灌丛生长季进行，在林场灌丛均匀分布地带设置3 个100 m×100 m 的样方，抽取20 株高度在1 m 左右的同龄灌木，记录高度、冠幅等植被特征参数，计算平均值作为确定3 类株型模型的参考依据，结果如表1 所示。本文主要探究二维尺度下植株周围的流场特征，为便于植株建模，对植株高度、冠幅进行均值处理及整数化处理。经过试算与模拟，确定3 类株型模型植被参数如表2所示。

表1 已有野外调查结果Tab.1 Field survey results available

表2 植株模型参数Tab.2 Model plant parameters

由表1可知，白刺的最大侧影面积出现在0.10～0.30 m 处，呈现底部大、顶部较窄的坛形，枝下高度较低，整体贴近地面，重心偏低；梭梭整体呈现中间宽两头窄的梭形，形态较为高大，重心在植株中间；沙拐枣的枝下高度最大，最大侧影面积出现在植株上部，呈现帚形，重心偏高。

1.2 数值模拟方法

1.2.1 植株建模 由于沙粒在风沙流运动中所受的拽力、重力基本在同一平面内，故本文采用二维简化模型。根据表2 及前人建模经验［14］，采用GAM⁃BIT 软件，充分利用灌丛形态的对称性，建立了白刺、梭梭、沙拐枣的孔隙几何模型，依次为坛形、梭形、帚形。简化植株模型是植株形态的近似和概化，通过较小的植株孔隙度以及各类实测植被特征参数来表现生长季灌丛的植株形态，因此，模拟结果具有一定的可靠性。通过多次试算和模拟，流域长度设置为30 m，计算高度为10 m，沙床厚度为0.05 m；植株简化模型距入口5 m，计算域简化图如图1所示。

图1 流场简化图Fig.1 Model diagram

模型形状较为规则，故网格划分类型采用Quad形式(四边形网格)，网格划分方法采用结构化网格划分。本文主要讨论沙粒的流场特征以及植株前后沙粒和沙床表面的相互作用，因此，对地面及植株附近的网格进行加密，经过多次模拟和调整，确定网格尺寸增长率为1.04，计算域网格总数为44647个，最小正交质量远远大于0.99，最大正交歪斜率远远小于0.001，网格质量良好，满足计算要求。

1.2.2 参数设置 结合前人［14］模拟参数设置，本文将风沙流视为不可压缩流体。模型左侧入口处采用速度边界条件，右侧采用压力出口边界条件，上壁面采用对称边界，下壁面及植株模型采用固体无滑移边界。风沙流粒径一般为0.075～0.25 mm，沙粒粒径设置为ds=0.1 mm，沙粒密度ρs=2650 kg·m-3，初始沙粒体积分数为0.02%，沙床初始沙床堆积率α=0.625，黏度μs=0.047 Pas，空气密度ρ=1.225 kg·m-3，空气运动黏度μ=1.7894×10-5Pas。入口边界速度为典型风速廓线流如下式所示：

式中：v0为摩阻速度（m·s-1）；y0为粗糙长度（m）；k为冯卡门系数，取值为0.4；y为高度（m）；v(y)为y处的风速值（m·s-1）。

求解模型采用欧拉双流体模型，附加湍流模型。由于需要考虑流场中速度变化、求解精度问题，故采用非定常瞬态求解方法来模拟，流场求解采用Phase Coupled SIMPLE 算法，空间离散格式采用二阶迎风格式，时间步长取0.001 s。

1.2.3 基本方程 本文模拟气流为不可压缩流，主要包含的控制方程有连续方程、动量方程和湍流模型方程，根据模拟需要，湍流模型选择适合高雷诺数的标准k-ε模型。

连续性方程为：

式中：ux、uy分别为两个方向的速度矢量；ρ为密度（kg·m-3）；t为时间（s）。

动量方程为：

式中：p为流体微元体上的压强（Pa）；ū为速度矢量；τxx、τxy、τyx、τyy都是因分子黏性作用而产生的作用在微元体表面上的黏性应力τ的分量（Pa）；g为重力加速度（m·s2）。

标准k-ε模型的湍动能k以及耗散率方程为：

1.3 分析方法

空气动力学粗糙度是近地表风速为零的高度，是反映地表对风阻抗的重要参数［15］，粗糙度越大，表明地面对风的削弱能力越强。其计算公式如下：

式中：Z0为空气动力学粗糙度（cm）；V1、V2分别是Z1、Z2两个高度的风速（m·s-1）。

植株防风效能计算公式如下：

式中：ηxz为(x，z)处的防风效率；x为距植株的水平距离（m）；z为距地面高度（m）；vxz为植株前后(x，z)处的风速值（m·s-1）；vz为裸沙地在z高度处的初始风速值（m·s-1）。

1.4 风洞试验概述

为验证模拟结果的可靠性，本文利用已有风洞试验［13］验证相关结论。该试验在中国林业科学院沙漠林业实验中心进行，采用直流开口吹气式风洞，宽、高均为2 m，试验段为30 m。风洞通过变频器调节风速，本文设计风速为6 m·s-1。由于本文主要利用风洞试验测定的各株型灌丛阻沙效益来验证相关模拟结果，因此，以下内容主要详述了有关风洞试验中测定植株阻沙效益的试验布置。采用阶梯式集沙仪测定灌丛周围不同高度层的输沙量。选择具有代表性的白刺、梭梭、沙拐枣灌丛，在样方内筛选高1 m 左右的同龄单株灌木，并将具有代表性的植株取回试验地，精确测定其植被特征参数。将采集到的3类生长季灌丛单株固定在风洞试验区的中轴线处，在株后1 H处布设一个1 m高、50层的阶梯式积沙仪，每层入口截面为2 cm×2 cm。为确保沙源充足，在试验段布设5 cm 厚的沙床，沙样经70 目筛网过滤。每种灌丛单株在6 m·s-1风速下进行吹蚀试验，吹蚀时间为2 min。吹蚀结束后取出各高度层的集沙槽，用精度为千分之一的天平测定各个高度层的积沙量，并计算输沙率。每类灌丛试验重复3 次，并在每次吹蚀试验结束后补充沙量并平整沙面。同时设置空白对照组，即在试验段不放任何植株的条件下，测得输沙率作为对照，对比分析各灌丛的阻沙效应。

2 结果与分析

2.1 风沙两相流场合理性验证

为验证本文风沙流场设置的合理性，将模拟结果与已有风洞试验结果［13］进行对比。研究表明风沙运动中沙粒体积分数与输沙率的分布状况一致［16］。因此，本文将风洞试验中6 m·s-1（即摩阻风速为0.3752 m·s-1）风速下，株后1 H处白刺和裸沙地的输沙率数据作为对照，分析同条件下数值模拟所得白刺沙粒体积分数随高度的分布值，结果如图2 所示。株后1 H处输沙率与沙粒体积分数随高度变化规律大体一致，均表现为随高度层的增加逐渐减小，总体呈现对数减小趋势，此规律符合沙粒运移特征和风沙流结构特点，说明模拟结果具有较高的可靠性，可以较为准确地反映植株周围的流场特征。

图2 输沙率和沙粒体积分数随高度分布Fig.2 Distribution of sand transport rate and sand volume fraction with height

2.2 单株灌木的防风效果

2.2.1 植株周围风速流场特征比较 为了对比3种株型植被对流场的影响，选择2 m 高度处风速为6 m·s-1的风速廓线进行模拟（下文速度均与此相同，不再赘述），图3展示了3种株型灌丛的风速分布特征。由图3 可知，在相同来流条件下，3 种植株上方均出现了加速区，迎风侧出现遇阻减速区，同时背风侧出现了紊流减速区、涡旋加速区以及恢复区。气流在株前-2 H 处出现扰动，形成遇阻减速区，其原因为灌丛枝型较为密集，枝条丛生增加受风面积，加速了风能的损耗。气流到达植株前方时，一部分气流汇聚抬升并在植株上方形成加速区；另一部分气流穿过植株，经过涡旋作用，在植株背风侧风速急剧下降，从而形成了紊流减速区，该区是植株的主要防护区，其大小直接影响植株的防护效果。株后，随着株距的增加，风速逐渐恢复至正常水平，从而植株失去防护能力。值得注意的是，植株周围出现了明显的涡旋加速区，3 类植株涡旋加速区的范围大小表现为沙拐枣＞梭梭＞白刺，而3者对应枝下高度依次为0.12 m、0.08 m、0.03 m，表明涡旋加速区呈现随枝下高度增加而明显增大的趋势，与已有研究结论相似［13］。

图3 灌木周围风速云图Fig.3 Cloud chart of wind speed around shrubs

2.2.2 株后涡流分析 已有研究表明，植株防沙的气动性原理主要是局部风速的降低和气流的循环泡（即涡旋流）［17］。图4展示了不同株型灌丛周围的风速矢量。从图中可以看出背风侧主要存在3个涡流，其中涡流I、II是主要防护区域，且受灌丛株型影响，两涡流的位置不同。白刺、梭梭及沙拐枣的涡流Ⅰ分别位于株后0.5～3.5 H、0.5～5 H、0.5～5.5 H处，涡流II 的位置依次为4～5 H、6～7 H、6.5～7.5 H。分析3 个涡流的形成原因如下：由于气流穿过植株时受到枝叶挤压，越过植株后过流断面突然增大，从而形成向上、向下运动的分流，向下运动的部分气流在距离植株不远处抬升，又返回植株方向，从而形成涡流I，并成为植被积沙的主要区域；涡流II形态虽较小，但强度较弱，也会导致大量沙粒在此沉积，其形成与涡流I 的下沉气流与涡旋加速区的气流相碰撞有关，气流挤压后形成涡流II；受涡流I、II 末端影响，气流再次被挤压形成涡流Ⅲ。植株阻沙能力与涡流的强度及尺度有直接关系，涡流范围越大、强度越弱，沙粒越易沉积，从而植株的阻沙效益越好。3 类植株涡流I 范围大小依次为梭梭＞沙拐枣＞白刺，强度依次为沙拐枣>梭梭>白刺。梭梭与沙拐枣受涡旋加速区的影响，涡流I虽范围大，但强度也较大。因此，在积沙初始阶段，流经梭梭、沙拐枣的挟沙气流主要在株后6～7 H 处积沙，而白刺主要沉积在株后3 H附近。

图4 灌木周围风速矢量图Fig.4 Vector diagram of wind speed around shrubs

2.2.3 植株后方风速廓线对比分析 为进一步探究3 类株型灌丛在垂直方向上的防护效果，在6 m·s-1风速条件下，提取3类株型株后1 H处垂直截面，绘制风速廓线图，并以裸沙地为对照。根据图5可知，3类株型的风速轮廓线存在明显差异。首先表现在近地表起始风速不同，白刺的起始风速最小，为1.57 m·s-1，梭梭与沙拐枣的起始风速均大于裸沙地，表明在株后1 H处的近地表区域，白刺具有抗风蚀作用，而梭梭、沙拐枣反而会促进植株周围发生风蚀，此现象值得在生态恢复工程中进一步研究。此外，3类株型的防护高度也存在明显差异，白刺、梭梭、沙拐枣的防护高度依次为1.1 m、1.5 m、1.65 m，沙拐枣的防护高度最高，与其植株高度相呼应。

图5 单株灌木后1 H处风速廓线图Fig.5 Wind speed profile of a single shrub 1 H later

进一步分析3类株型灌丛与裸沙地的风速廓线还可以发现，裸沙地风速轮廓线基本呈现对数形式增长，而灌丛株后风速轮廓线并没有形成规律的对数函数曲线，而是随高度增加呈现竖向“W”变化的趋势，气流垂向分布可依据其变化趋势划为低空减速区、涡旋区、高空加速区。受植株的阻挡作用，白刺风速廓线在0.05～0.3 m 高度范围内逐渐减小，并在0.3 m 高度处达到风速极小值，风速由1.57 m·s-1降低至0.18 m·s-1，风速降幅为88.53%，该区域即为低空减速区；随后气流在0.3～0.8 m高度层内进入涡旋区，风速随高度增加呈现先小幅度增大后减小的趋势；在0.8 m高度处以后，风速以较大的加速度呈现快速增长的趋势，植株也逐渐失去防护能力，此区即为高空加速区。与白刺不同，梭梭的低空减速区出现在0.1～0.4 m 高度层内，风速从3.81 m·s-1降低至0.18 m·s-1，风速降幅为95.27%，并在0.4 m高度出现风速极小值；随后出现涡旋区，风速在0.4～0.8 m 出现小范围波动，并在0.8 m 处风速达到次小值；气流在0.8 m 以上呈现随高度增加而增大，即高空加速区。沙拐枣具有相似的变化规律，其低空减速区为0.05～0.4 m，风速降幅为64%，0.4～0.8 m为涡旋区，风速出现先增长后降低的特征，并在0.8 m高度出现风速极小值，随后进入高空加速区。总体而言，不同株型的风速轮廓线存在差异，风速极小值出现高度不同，各植株最优防护高度也存在明显区别，白刺、梭梭、沙拐枣的最优防护高度依次为0.2～0.4 m、0.3～0.6 m、0.8～1 m。值得注意的是，0.8 m 以上高度，各植株均呈现风速随高度增加而增大的趋势，表明0.8 m处是风速发展的重要分界点，此现象值得进一步探究。

2.2.4 株后空气动力学粗糙度分析 由公式（7）可知，观测高度和风速的选择是计算空气动力学粗糙度的关键。根据前文分析可知，0.5～2 m间风速轮廓线的变化可以代表整个灌丛层的风速特征，故选择3类株型植被背风侧和裸沙地0.5 m和2 m两个高度的风速值，并计算得出3 类株型株后不同株距处的空气动力学粗糙度，计算结果如图6所示。根据图6可知，随株距增大，3 类株型株后的粗糙度逐渐减小，并在15 H 处与裸沙地粗糙度基本持平。3 类株型在株后1 H 处的粗糙度最大，白刺、梭梭、沙拐枣的粗糙度依次是裸沙地的26倍、20倍、13倍。白刺的空气动力学粗糙度明显高于其他两种株型，沙拐枣的粗糙度最小，表明白刺削弱近地表风力的能力更强，具有更好的防风作用。

图6 灌丛株后空气动力学粗糙度Fig.6 Aerodynamic roughness of shrub plants

2.2.5 植株防风效能分析 分别提取6 m·s-1风速条件下各高度层不同株距下的风速值，通过公式（8）计算单株灌木的防风效能。分析图7，株前水平防风效益随株距减小而增大，尤其在-2～0 H 范围内，防风效益剧增；垂向防风效益随高度增加而减小，但不存在负值，表明各灌丛在其高度范围内对风速均具有明显的削弱作用。株后根据各高度层防风效益的变化趋势，可将其分为近地表区（0.1 m、0.2 m、0.3 m）和中高空区（0.5 m、0.7 m、0.9 m）。在近地表区，株后水平防风效率随株距增加呈现波浪式变化；垂向防风效益在株后5 H范围内，呈现随高度增加而增大的趋势。在中高空区的垂向防风效益上，白刺、梭梭防风效益随高度增加呈现波动式降低，而沙拐枣呈现波动式增大；在水平方向上整体呈现随株距增加而降低的趋势，白刺、梭梭、沙拐枣的防风效益分别在12 H、11 H、13 H处基本下降至50%，且在该区内，株后1 H 处防风效果最好。以沙拐枣为例，3 个高度层最大防风效益依次为89.8%、98.1%、94.4%，均出现在株后1 H处。

图7 灌木植株防风效益Fig.7 Protection benefit of shrub plants

进一步分析图7 可知，除梭梭、沙拐枣在0.1 m高度层防风效益出现负值外，各株型在低于其高度内均产生了有效的防护，但各灌丛的最优防护区域存在明显差异。白刺最优防护区域出现在0.3 m高度层、株后5 H范围内，此区域内植株防风效益均大于90%；梭梭最优防护区域出现在0.5 m 高度层、株后4 H 范围内，而沙拐枣则出现在0.9 m 高度层、株后2～5 H 范围内。总体而言，各株型在-2～10 H 范围内均可有效降低风速，但3 类株型灌丛的区域防风效益显著不同，白刺在近地表区具有较好的防风效果，而梭梭、沙拐枣在植株中高区域具有较好的防风效益。

2.3 单株灌木周围沙粒运动特征分析

2.3.1 植株前后输沙特征分析 当风速达到起沙风速时，地表沙粒就会在风力作用下开始以蠕移、跃移、悬移的形式运动，进而形成风沙流。灌丛能够拦蓄风沙流中的沙粒，使其沉积在植株周围，从而形成相应的灌丛沙丘，但不同株型灌丛的阻沙能力不同。为了探究3 种典型株型荒漠灌丛的阻沙作用，图8给出了3类株型植株在6 m·s-1风速条件下，T=4 s、T=10 s的积沙云图。白刺前期迎风坡积沙量较少，沙粒主要沉积在背风坡，但与植株保持一定距离；随时间增加，背风坡积沙点逐渐向植株方向靠近，T=10 s 时，株后积沙点已从0.5 H 处蔓延至植株底部，且背风坡积沙长度由3 H 变化至5 H，株前积沙长度由0.5 H增加至3.5 H，此时株前、株后总积沙长度为8.5 H。梭梭迎风坡整体积沙量较少且主要集中在植株附近，背风坡积沙量较多，积沙距离同样表现出随时间增加逐渐靠近植株的趋势，T=10 s 时，株后积沙点已从2 H 处到达1 H 处，株后积沙长度大致从4 H 变化至5 H，最大积沙高度从3 cm增长至6 cm，10 s时总积沙长度为6 H。沙拐枣迎风坡积沙最少，背风坡在4 s 时存在较大的积沙距离，表现为在距植株2.5 H 处积沙，随时间增加至10 s，株后积沙点到达1 H 处，总积沙长度增加量为1.5 H，积沙高度增加量为5 cm，10 s 时总积沙长度为4.5 H，相比其他株型，沙拐枣的阻沙效益相对较差。

图8 灌丛周围积沙云图Fig.8 Sand accumulation cloud atlas around shrubs

2.3.2 植株前后地表蚀积特征分析 在6 m·s-1的风速条件下，气流与沙床相互作用，撞击和挤压使得沙床上的沙粒以多种形式运动，从而导致沙床上的沙粒发生重排列，进而使沙床堆积率发生变化。在实际工作中，沙床堆积率不易精准检测，但数值模拟可以解决这一难点。为了探究不同形态灌丛植株前后的地表蚀积状况，本文初始条件设置了厚度0.05 m 且均匀的沙床，单一粒径初始堆积率为0.625，通过提取T=10 s 时刻的沿程沙床堆积率，反映不同株型对区域蚀积状况的影响。根据植株周围风沙运移规律以及前人相关研究［18］，文中提出假设：若0.05 m 以下高度层沙粒体积分数低于0.625，且同时满足0.05 m以上高度层沙粒体积分数在0.02附近，则表明该区域发生风蚀；若0.05 m 以上高度层体积分数均大于0.02，且此刻0.05 cm以下高度层的沙粒体积分数在0.625 附近，则说明该区域出现风积现象。

从图9a 可以明显看出，风沙流在距入口3 m 左右处受到白刺影响，0.05 m 及以下高度层沙粒体积分数始终保持在0.625 左右，0.05 cm 以上高度层沙粒堆积率随株距缩短而增加，表明株前发生了风积作用。在白刺周围，0.05 cm以上高度层沙床堆积率发生轻微扰动，但其堆积率始终高于0.5，表明植株附近发生了风积作用。随后气流逐渐恢复，但受涡流区II、Ⅲ的影响形成弱风蚀区、风积区。梭梭和沙拐枣在株前一定范围内也存在不同程度的风积作用，但与白刺不同的是，两者在植株下方或株后一定范围内，高空沙粒堆积率虽增加，但梭梭在株下0.04 m 高度层、沙拐枣在株下0.04 m、0.05 m 高度层，沙粒体积分数均出现低于0.02 的状况，表明该区域发生风蚀作用，风蚀区域范围分别距入口5～5.5 m、4.5～6 m；气流穿越植株后，两者分别在株后距入口约6～10 m、6～9 m发生风积作用。总体而言，3类株型灌丛对地表蚀积微地貌的影响存在显著差异，株前、株后方均存在一定程度的风积状况，但除白刺外，其他两种植株周围出现不同程度的风蚀现象。

图9 灌木周围蚀积状况Fig.9 Erosion around shrubs

3 讨论

3.1 不同株型灌丛防风能力差异的原因

沙生灌丛的形态差异是不同灌丛对空间、光、水等资源的适应策略，灌丛在适应不同生存环境时的反应特征主要体现在其外部形态上［19］。生存环境影响植株形态的发育进程，而植株形态也会影响其自身的防风阻沙效益［20］。本文通过模拟发现，单株灌丛周围的流场可以划分为5 个区域：遇阻减速区、上方加速区、紊流减速区、涡旋加速区、恢复区。其中，涡旋加速区的形成与植株枝下高度相关。较大的枝下高度会使风沙流穿过时，附近流场产生“狭管效应”［21］。各株型灌丛的枝下高度不同，导致涡旋加速区范围存在明显差异，并呈现随枝下高度增加而增大的趋势。各株型背风侧均存在3个涡流，涡流强度影响沙粒的沉积位置。由于坛形灌丛的枝叶密集，迎风面透风系数较小，导致涡流Ⅰ强度较小，即在风影区风速降幅较大，并在风影区形成风影堆积。因此，在积沙初始阶段，白刺主要在株后3 H 处积沙。梭形、帚形植株底部枝叶分布稀疏，疏透度较大，气流经过植株孔隙时受到挤压而加速，形成强度较大的涡流Ⅰ。相比涡流Ⅰ，涡流Ⅱ的强度较小，从而导致梭形、帚形植株最初在株后6～7 H处积沙。

3 类灌丛株后1 H 处的风速廓线图呈现竖向“W”变化趋势，气流垂向分布可划分为低空减速区、涡旋区、高空加速区，其分布规律与植株形态相关。受枝下高度的影响，起始风速较大，后受植株、沙粒蠕移、跃移等因素影响，风速随高度增加而减小，从而出现低空减速区；涡旋区的出现与植株本身的复杂枝系结构相关，使气流穿过复杂孔隙时出现小尺度湍流；高空加速区则是受植株阻挡影响，风速随高度增加而增大。3 类株型灌丛株后1 H 处的风速极小值依次出现在0.3 m、0.4 m、0.8 m 高度处，最优防护高度依次为0.2～0.4 m、0.3～0.6 m、0.8～1 m，不同株型的防护范围及最优防护高度均存在显著差异。研究表明，防风效应取决于迎风面侧影面积的大小，植株形态决定其侧影面积大小及分布，从而影响其防护范围及最优防护高度［22-23］。白刺形态低矮，呈现上窄下宽的坛形，近地表区域枝系繁茂，最大侧影面积出现在0.1～0.3 m 高度层，因此白刺可以有效降低近地表区风速，其最优防护高度为0.2～0.4 m。梭梭整体形态呈现中间宽两头窄的梭形，具有一定的枝下高度，植株底部形成涡旋加速区，从而容易导致局部风速变大，但由于其迎风面侧影面积较大，且最大侧影面积出现在0.3～0.6 m 高度层，因此，其具有较好的防风效果，最优防风高度出现在植株中间高度处。沙拐枣呈现上宽下窄的帚形，具有较大的枝下高度，导致其在近地表的防风能力较弱，且容易引起植株根部周围的侵蚀，但由于它形态高大，最大侧影面积出现在0.6～1 m高度，导致其在0.6 m以上高度层具有较好的防风效果，并在0.8～1 m表现出最优防护高度。

总体而言，各株型灌丛的防风能力存在显著差异。坛形植株具有较为紧密的枝叶结构，且其底部侧影面积较大，可以有效降低近地表风速；梭形、帚形植株虽然底部枝叶相对稀疏，容易形成“狭管效应”，但两者形态高大，冠层中高部分透风系数小，因此，其在中高区域具有较好的防风效益。

3.2 不同株型灌丛阻沙能力差异的原因

3类株型灌丛的阻沙能力呈现白刺＞梭梭＞沙拐枣，总积沙长度依次为8.5 H、6 H、4.5 H。与防风效益相似，株型也会影响植株的阻沙能力。已有研究表明，半球形冠形植株比锥体冠形及梭形冠形植株的截沙能力更强［24］。本文通过研究发现白刺的积沙范围最大，梭梭次之，沙拐枣最小，且梭梭、沙拐枣在植株周围存在明显风蚀区，因此，相较于梭梭、沙拐枣，白刺具有更好的防风阻沙能力，此结论与已有试验结果相似［25-26］。其原因为风沙活动主要发生在近地表0～0.3 m 处，而白刺呈现下宽上窄的坛形，贴近地面生长，当风沙流经白刺时，底部枝系结构紧密，对沙粒有明显的阻挡作用，导致越过和透过植株的沙粒减少，起到良好的阻沙作用，与其容易形成沙堆的特点一致［27］。沉积在3类灌丛周围的积沙量始终表现出背风侧高于迎风侧，产生上述积沙迁移规律的主要原因为：植株前后对风沙流的削减能力存在差异，风沙流流至植株时，迎风侧风速虽有所降低，但大部分气流受挤压抬升，此时部分动能转化为重力势能，加之风沙流流经植株内部孔隙时与植株相撞，损失的动能较大，气流挟沙能力进一步降低，因此，更多沙粒堆积在背风坡，甚至形成小沙堆。

就灌丛对地表蚀积微地貌影响的研究中，本研究发现3 类株型灌丛株前、株后方均存在一定程度的风积状况，但梭梭、沙拐枣分别在距入口5～5.5 m、4.5～6 m处存在不同程度的风蚀现象，这与两者存在较大的枝下高度有关。白刺株前部分风沙流从植株周围加速通过，带走部分沙粒，但由于白刺近地表区冠层侧影面积大，导致植株下方沙粒沉积量远大于吹蚀量，因此，其株下并未发生风蚀而是出现沙粒沉积。而流经梭梭、沙拐枣的风沙流受“狭管效应”的影响，在植株下方加速通过并带走大部分沙粒，说明两类灌丛根部及周围存在一定的风蚀现象。

在土壤风蚀防治研究中，研究单一灌木或灌丛防护效益的同时，还要考虑灌木林带的群体防护效应和机制［28-29］。合理配置灌丛可以使水分条件限制下的有限植被达到良好的防治风蚀效果。利用数值模拟研究单一植株周围流场的目的之一是设计更有效的稀疏灌丛防护林，如Guo 等［30］通过数值模拟发现行距在2.5～3 H 的沙蒿防护林防护效益最好，而柳树防风林的最优间距为1～1.5 H。本研究发现白刺、梭梭、沙拐枣的最优防护高度依次为0.2～0.4 m、0.3～0.6 m、0.8～1 m，而3者的阻沙能力依次为白刺＞梭梭＞沙拐枣。综上所述，在实际风蚀防治工程中，可采用白刺与梭梭、沙拐枣相结合的配置方式，既能发挥白刺的阻沙作用，又能利用梭梭、沙拐枣在中高空区较大的防风效应。

4 结论

本文利用欧拉双流体非定常模型，对3 种典型株型荒漠灌丛周围的风沙流场进行了研究，并通过已有风洞试验验证了流场合理性，得出以下结论：

（1）单株灌丛周围的流场可分为5 个区域，其中，涡旋加速区呈现随枝下高度增加而增大的趋势。株后均出现3 个涡流，但同一涡流的强度存在差异，导致在积沙初始阶段，梭形、帚形植株主要在株后6～7 H处积沙，而坛形主要在株后3 H处积沙。

（2）3 类株型株后1 H 处的风速轮廓线呈现竖向“W”变化趋势。受植株最大侧影面积高度层的影响，坛形、梭形、帚形的风速极小值依次出现在0.3 m、0.4 m、0.8 m 高度处，最优防护高度依次为0.2～0.4 m、0.3～0.6 m、0.8～1 m。3 类株型株后的空气动力学粗糙度逐渐减小，且坛形植株的粗糙度明显高于其他两种株型。

（3）3 类株型在-2～10 H 范围内均可有效降低风速，株后近地表区防风效益表现为坛形＞梭形＞帚形；中高空区防风效益均随株距增加而减小。3类株型最优防护范围存在明显差异，坛形植株最优防护区域出现在0.1～0.3 m 高度、株后5 H 范围内；梭形最优防护区域在0.3～0.6 m高度、株后4 H范围内；而在0.6～1 m、株后2～5 H 范围内帚形防风效益均高于90%。

（4）在T=10 s时，灌丛株前、株后方均存在一定程度的风积状况，白刺、梭梭、沙拐枣周围总积沙长度分别为8.5 H、6 H、4.5 H。受“狭管效应”的影响，梭梭、沙拐枣分别在距入口5～5.5 m、4.5～6 m处存在不同程度的风蚀现象。相比而言，白刺的阻沙效果较好，但在实际防风固沙工程建设中应使其与梭梭、沙拐枣结合，既能发挥白刺的阻沙作用，又能利用梭梭、沙拐枣在中高区较大的防风效果。