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

    典型株型沙生灌丛对风沙流场影响的数值模拟    闫 晴，李菊艳，尹忠东，刘金苗，柳宏才（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 s1.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 height2.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 shrubs2.2.2 株后涡流分析 已有研究表明，植株防沙的气动性原理主要是局部风速的降低和气流的循环泡（即涡旋流）［17］图4展示了不同株型灌丛周围的风速矢量从图中可以看出背风侧主要存在3个涡流，其中涡流I、II是主要防护区域，且受灌。