正弦形沙丘背风坡回流区特性研究

祁紫薇，王 硕，凡凤仙，胡晓红

（1.上海理工大学 能源与动力工程学院，上海 200093；2.上海理工大学 上海市动力工程多相流动与传热重点实验室，上海 200093）

沙丘是不同粒径的颗粒在风力的长期作用下堆积而形成的自然地形地貌的一种，其形态、大小各不相同。沙丘的形成和演变过程强烈受制于风沙两相流运动规律，与气流特征、颗粒运动紧密联系，是一种极其复杂的物理过程。在研究沙丘形成和演变动力学时应充分考虑沙丘外形、沙丘表面流场、沙粒粒径等多方面因素。其中，沙丘表面流场决定着局地沙粒输运的强度，在沙丘演变动力学研究中占据重要地位。在沙丘表面流场的3 类研究方法中，野外观测研究对象单一，而且受到恶劣的观测环境以及仪器精度的限制，加之背风坡的湍流特性十分复杂，依靠野外观测难以对背风坡流场结构进行精确描述；风洞实验存在同时保证几何相似和动力学相似的困难；数值模拟能够方便、灵活地调整计算参数，便于获得流场结构的细节信息，是研究沙丘表面流场特性的一种有效路径[1-2]。

在沙丘表面流场数值模拟研究中，计算流体力学软件PHOENICS[3]、FLUENT[4-6]、OpenFOAM[2,5,7-8]等以及自编程序[9]先后得到应用，所采用的描述流场的数学模型包括标准k-ε、RNGk-ε、k-ωSST 湍流模型、大涡模拟等，文献[2]对此作了较为详尽的评述。由于不同的研究者所针对的沙丘形貌和尺度差别很大，采用的计算软件及湍流模型也有所不同，所得的结论缺乏普适性。此外，已有数值模拟研究往往基于野外观测或风洞实验数据考量模型的准确性，对沙丘表面流场特性，特别是正弦形沙丘背风坡回流区特性的研究仍不够。笔者所在课题组曾借助OpenFOAM 的标准k-ε模型对正弦形沙丘表面流场开展数值模拟，将模拟结果与实验结果相结合验证了数值模拟的正确性和网格敏感性，给出了沙丘表面内层和外层的速度分布曲线。在此基础上，本文采用与文献[2]一致的k-ε模型，通过改变沙丘的高度和宽度，开展数值模拟，考察沙丘表面流场的结构与背风坡回流区特性，以确定正弦形沙丘几何参数对其表面背风坡回流区的定量影响，为进一步研究沙丘演变动力学提供基础和参考。

1 计算模型

描述湍流的模型，如标准k-ε、RNGk-ε、k-ωSST 湍流模型、大涡模拟等各有特点。标准k-ε模型是工程流场计算中主要的工具，其适用范围广，具有合理的计算成本和较高的计算精度，在流场模拟中具有广泛的应用。RNGk-ε模型相比标准k-ε模型修正了耗散率方程，对于复杂的剪切流，具有大应变率、旋涡、分离等的流动问题，具有比标准k-ε模型更高的精度。k-ωSST 模型包含了来自ω方程的交叉扩散，在湍流粘度中考虑了湍流剪应力的传播，使得其用于流场模拟具有更高的精度。大涡模拟介于直接数值模拟（DNS）与雷诺平均纳维-斯托克斯（RANS）方法之间，大尺度漩涡直接模拟，小尺度漩涡用RANS 方法求解，在实际工程应用中需要很大的计算代价。由于在前期研究中，利用标准k-ε模型取得了与实验吻合良好的结果[2]，而标准k-ε模型计算成本相对较低，因此本文模拟中采用标准k-ε模型。

沙丘表面流场可用笛卡尔坐标系下的雷诺平均纳维-斯托克斯（RANS）方程进行描述，可写为

式中：t为时间；u为流速；p为压力；ρ为空气密度；μ为空气动力粘度；δij为克罗内克函数，当i=j时，δij=1，当i≠j时，δij=0。

为使方程封闭，标准k-ε湍流模型给出了湍动能k和耗散率ε的方程如下：

式中：σk，σε，C1ε，C2ε为模型常数，在大气边界层模拟中通常采用σk=1.0，σε=1.30，C1ε=1.21，C2ε=1.92[10-12]；μt为湍流粘度；Eij为平均应变率张量。

式中，模型常数Cμ=0.03[10-12]。

2 计算条件及方法

2.1 计算区域

本文的研究是针对文献[13]的沙丘形貌进行的，即沙丘为狭长型沙丘（长30 m），如果按实际尺寸进行模拟，计算网格数目庞大，导致计算成本很高。由于沙丘表面速度变化主要集中在来流方向（x向）和竖直方向（z向），因而可忽略y向的速度分量，开展二维模拟。事实上，已有研究也表明二维模拟能够得到和实验吻合良好的结果[2]。图1 给出了二维计算区域示意图，沙丘位于计算区域底面中心，其轮廓曲线为

式中：zd为沙丘表面的z坐标；A为正弦曲线的振幅；k为波数，k=2π/λ，λ为波长。

图1 计算区域示意图Fig.1 Schematic diagram of the simulation domain

数值模拟时，为避免边界效应，进口、出口、顶面与沙丘的距离应足够大。为此，计算区域设置时取Δx=4.35λ，Δz=18.20A，研究中发现，如果增大Δx和Δz，计算结果几乎不发生变化。

2.2 网格划分

为减少计算量，在划分网格时，可以对于速度变化较大的区域采用较密的网格，对于速度变化较小的区域采用较疏的网格。基于此，将由Δx和Δz确定的沙丘表面流场计算区域划分为3 个子区域，如图2 所示。其中，区域III 的大小由=λ×（4.55A）确定，区域II 的大小由（1.09λ）×（9.10A）确定。确定网格大小的原则为：区域I 的网格最大，为0.4 m×0.4 m；区域Ⅱ的网格边长为区域I 的1/2，区域Ⅲ的网格边长为区域Ⅱ的1/2。此外，在计算区域底面（即包含沙丘的地面）添加5 层边界层，在远离下表面的方向上，边界层网格厚度成比例增加，比例系数为1.1。

图2 网格划分所采用的3 个区域Fig.2 Three regions used for mesh generation

2.3 边界条件

进口采用速度边界条件，速度值可由大气边界层风速的对数分布率给出[13]

式中：κ为冯·卡门常数，κ=0.41；u*为剪切速度；z0为空气动力学粗糙长度。数值模拟中，选择u*=0.5 m/s，z0=0.83 mm。

进口湍动能k和耗散率ε的表达式为[14-15]

底面采用壁面边界条件，顶面设置为移动壁面，顶面风速由式（8）确定，出口设置为压力出口。此外，基于OpenFOAM 进行二维数值模拟时，需要在y向采用单层网格，并将相应的边界设置为“empty”。

2.4 求解控制

数值模拟中，采用SIMPLE 算法[16-18]对压力和动量方程进行耦合求解，收敛条件为质量和动量通量的残差达到10-8，压力的残差达到10-7；同时，压力的松弛因子设置为0.3，速度、湍动能和耗散率的松弛因子均设置为0.7；所采用的插值格式见表1。

表1 数值模拟采用的插值格式Tab.1 Interpolation schemes used in numerical simulations

3 结果与讨论

3.1 沙丘高度对流场的影响

正弦形沙丘的振幅反映了沙丘的高度，为探讨沙丘高度对沙丘表面流场的影响，保持波长λ=46 m 不变，在振幅A=1.65，2.20，2.75，3.30 m情况下开展数值模拟。模拟得到的沙丘表面不同高度处的速度曲线如图3 所示，其中图（a）～（d）分别对应于上述4 种振幅情况，u0为对应曲线上最大速度和最小速度的平均值。由图3（a）可以看出：在沙丘高度较小时，如A=1.65 时，在距沙丘表面较高的高度范围内（1.5～5.0 m），速度的分布与沙丘外形类似，即呈对称形，速度最大值出现在沙丘坡顶位置；在贴近沙丘表面的区域（沙丘表面0.5 m 高度以内区域），速度最大值出现在沙丘迎风坡，即坡顶上游，并且受沙丘轮廓影响，背风坡流速先降低后增加，但速度值始终为正值，这表明在沙丘的背风坡无明显回流区。对比图3（a）～（d）可知，随着沙丘高度的增加，在A=2.20 m 时，贴近沙丘表面的区域中背风坡流速出现负值，表明背风坡流场出现了回流区。图中，流速为负值的沙丘宽度（x）和高度区间反映了回流区的范围，负的流速在数值上反映了回流区的强度。可见，随着沙丘高度的增加，背风坡回流区范围变大、强度增强。图3 也反映出与沙丘迎风坡相比，背风坡流速的变化更为剧烈，因此背风坡沙粒受到的剪切力作用也更大；在沙丘背风坡出现回流区时，沙丘表面沙粒还将受到涡流影响，因此背风坡沙粒较迎风坡沙粒更易于发生迁移。与此同时，贴近沙丘表面的气流经过沙丘坡顶后，在背风坡转而向下，同时水平方向速度迅速降低并出现回流，可以推测迎风坡加速气流挟带的沙粒在输运到背风坡时，倾向于降落在回流区的起始位置。

图3 沙丘高度对沙丘表面流场的影响Fig.3 Effect of dune height on the flow field over the dune surface

3.2 沙丘表面回流区尺度的确定

沙丘背风坡回流区的存在对沙丘的移动以及演变有着重要影响，在对沙丘表面流场进行定量研究时，回流区尺度（x向尺度用回流区长度表示，z向尺度用回流区高度表示）的确定是一个关键的问题。为确定回流区尺度，对数值模拟结果进行处理时，首先绘制沙丘表面流场的流线图，根据图中气流的回流与分离来确定回流区尺度。图4 给出了不同高度沙丘的表面流场回流区。从图4 可以看出，随着沙丘高度的增加，回流区的长度和高度以及位置变化明显，图中示出了回流区的长度。结果表明：在振幅增加过程中，回流区长度增大；沙丘高度较小时，在沙丘背风坡的底部贴近沙丘表面处形成扁平的连续涡流；随着沙丘高度的增加，回流区中心从贴近沙丘表面的位置逐渐上移，回流区范围逐渐增大。图中所示3 种情况下，回流区长度与沙丘高度（2A）之比依次为1.02，3.33，4.23，图4（b）和（c）中回流区长度范围与Wiggs[19]的野外观测数据相符，即回流区长度通常为沙丘高度的1.6～5.4 倍。但是，需要说明的是，野外观测是对高度为9.6 m、坡顶与坡脚距离为86 m、宽度为130 m 的新月形沙丘进行的。此外，江丽娟[6]对高度3～15 m，背风坡坡度为32°的新月形沙丘表面流场进行数值模拟研究，也得到回流区长度为沙丘高度的2～8 倍。本文针对的沙丘形状与这些研究不同，但得出的回流区长度却与之相符。

图4 不同高度沙丘的回流区Fig.4 Recirculation zones of dunes with different heights

3.3 沙丘宽高比对回流区尺度的影响

图5 给出了不同计算条件下回流区长度随宽高比的变化关系。数值模拟时，采用了不同的计算区域设置，其中“原计算区域”表示按Δx=4.35λ、Δz=18.20A（见图1）的原则确定的计算区域，“区域放大后”表示对计算区域长度和高度分别放大1.5 倍后的计算区域。对比相同沙丘几何参数下“原计算区域”和“区域放大后”对应的沙丘背风坡回流区长度的计算结果，可知两种计算区域下的计算结果吻合良好，显示出本文数值模拟所选择的计算区域对于沙丘表面流场回流区特性研究是合理的。沙丘轮廓线振幅A恒定时的曲线代表了回流区长度随沙丘宽度的变化情况，可见：随着沙丘宽度的增加，回流区长度减小；并且沙丘高度越大，回流区长度的减小越迅速。沙丘轮廓线波长λ恒定时的曲线代表了回流区长度随沙丘高度的变化情况，可见，随着沙丘高度的减小，回流区长度迅速减小。总地说来，随着宽高比的增加，回流区长度减小。然而，回流区长度不仅取决于宽高比，还取决于沙丘宽度与高度的数值。图6 给出了回流区高度随沙丘宽高比的变化情况。可见，沙丘宽高比越大，回流区高度越小。这表明在沙丘宽高比发生变化时，回流区在横向和纵向的尺度同时发生变化，沙丘背风坡回流区对沙丘几何参数很敏感。

图5 回流区长度随沙丘宽高比的变化关系Fig.5 Relationship between the length of recirculation zone and the aspect ratio of the dune

图6 回流区高度随沙丘宽高比的变化关系Fig.6 Relationship between the height of recirculation zone and the aspect ratio of the dune

3.4 回流区成因分析

图7 给出了λ=23 m、A=1.1 m 情况下，沙丘背风坡的涡流分布情况。由于在沙丘背风坡气流受到的剪切应力不均匀，造成气流流动不稳定，容易形成涡流。沙丘背风坡的回流区刚出现时极不稳定，由一系列附着于沙丘表面的小涡组成，如图7 所示。对于理想的平坦表面，定常流动条件下，气流为边界层流动，不发生流动分离和回流。沙丘的存在，改变了气流的流动特性。对于正弦形沙丘，沙丘高度增加或宽度减小造成沙丘坡度增大，对风的阻碍作用增强，使气流通过沙丘坡顶和背风坡坡脚时，产生更大的压力差，从而在背风坡更易产生漩涡，形成回流区，进而影响沙粒的迁移。需要说明的是，本文模拟着重从流场的角度探讨了沙丘表面回流区特性，同时考虑沙粒的运动，利用气固两相流理论研究沙丘的演变过程将是下一步需要重点开展的工作。

图7 沙丘（λ=23 m，A=1.1 m）背风坡的涡流分布情况Fig.7 Eddy distribution at leeward side of the dune（λ=23 m,A=1.1 m）

4 结论

基于OpenFOAM 的标准k-ε模型，利用数值模拟方法研究了不同几何参数条件下正弦形沙丘表面流场特性，特别是背风坡回流区特性，给出了不同条件下回流区的特点，以及回流区长度和高度随沙丘宽高比的变化关系，分析了回流区的成因，得到以下结论：

a.沙丘高度较低时，在背风坡贴近沙丘表面的区域气流速度先降低后增加；随着沙丘高度的增加，背风坡气流速度降低更为迅速，气流不稳定，出现回流区。

b.在沙丘宽度或高度一定的情况下，随着宽高比的增加，回流区尺度减小；但是回流区尺度不仅受沙丘宽高比的影响，还取决于沙丘宽度与高度的数值。

c.沙丘背风坡回流区萌生时贴近沙丘表面出现一系列小涡，且分布不规律，随着沙丘高度的增加或宽度的减小，背风坡气流不稳定性加剧，涡的尺度增大，强度增加，演变为回流区。