上海理工大学学报  2017, Vol. 39 Issue (4): 313-319   PDF    
基于OpenFOAM的正弦形沙丘表面流场特性研究
王硕1,2, 凡凤仙1,2, 张盟盟1,2, 胡晓红1,2     
1. 上海理工大学 能源与动力工程学院, 上海 200093;
2. 上海理工大学 上海市动力工程多相流动与传热重点实验室, 上海 200093
摘要: 基于开源计算流体力学软件OpenFOAM,对沙丘表面流场进行数值模拟,将所得到的沙丘表面风速分布与野外观测结果进行对比,验证数值模拟的准确性.在此基础上,获得了距沙丘表面不同高度的速度分布曲线和沙丘表面上的速度矢量场,探讨了数值模拟结果对计算网格的敏感性,考察了沙丘宽高比对流场特性的影响.结果表明:OpenFOAM的标准k-ε模型能够预测低矮正弦形沙丘表面流场的速度转变特性;在贴近沙丘表面的内层,沙丘表面特定高度处的速度分布曲线相对于沙丘轮廓存在相位提前,随着高度的增加,相位提前量减少,速度变化趋于平缓;网格数目对外层速度分布影响较小,但对内层的速度分布则有重要影响;足够多的网格数目是获得精确的沙丘表面流场信息的必要条件,在沙丘流场数值模拟中应当引起重视;保持沙丘宽度不变,随着沙丘宽高比的减小,沙丘表面速度分布曲线的相位提前量变大,背风坡出现回流区,宽高比越小,回流区的长度与强度越大.
关键词: 计算流体力学     正弦形沙丘     湍流     速度转变     网格敏感性    
Flow Behavior over a Sinusoidal Dune Based on OpenFOAM
WANG Shuo1,2, FAN Fengxian1,2, ZHANG Mengmeng1,2, HU Xiaohong1,2     
1. School of Energy and Power Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China;
2. Shanghai Key Laboratory of Multiphase Flow and Heat Transfer in Power Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China
Abstract: Numerical simulations were performed based on the open source computational fluid dynamics (CFD) software OpenFOAM.The numerically obtained velocity profiles over the dune were compared with the field measured results to validate the accuracy of the numerical simulation.On this basis, the velocity profiles at different heights over the dune as well as the velocity vector field were obtained.The sensitivity of the numerical simulation result to the computational mesh was discussed.Furthermore, the effect of the aspect ratio of the dune on the flow behavior was examined.The results show that the standard k-ε turbulence model in OpenFOAM can be used to predict the velocity shift behavior of the flow field over the low sinusoidal dune.In the inner layer close to the dune surface, a phase advance of velocity curve relative to the dune profile is observable within a certain distance over the dune.As the height increases, the phase advance decreases and the velocity begins to flatten.The mesh number has little effect on the velocity profile in the outer layer, while it affects greatly on the velocity profile in the inner layer.It is worth noting that sufficient mesh number is required to attain accuracy information of the flow field over the dune. When the dune width keeps constant, with the decrease in the width-to-height ratio, the phase advance of the velocity profile over the dune increases and the recirculation zone could appear at the leeward.It is found that when the width-to-height ratio becomes smaller, the length and the intensity of the recirculation zone are greater.
Key words: computational fluid dynamics (CFD)     sinusoidal dune     turbulent flow     velocity shift     mesh sensitivity    

沙丘是指由于沙粒在流场中相互作用所形成的堆积体, 常见于海岸、河谷和干旱沙漠地带.沙丘表面流场决定着局地沙粒输运的强度, 从而影响沙丘的形态, 在沙丘演变动力学研究中占据重要地位, 因而备受研究者关注[1-3].目前对沙丘表面流场进行研究的方法可概括为3大类:野外观测[4-6]、风洞实验[7-8]和数值模拟[9-12].野外观测是最为直接的方法, 但是在进行大空间尺度和时间尺度的沙丘研究时, 实验测试和数据处理工作量将会非常复杂庞大, 且难以给出流场结构的详细信息.风洞实验可以较好地反映沙丘环境, 但是在满足相似性条件时, 难以同时保证几何相似和动力学相似.随着计算流体力学的发展和日趋成熟, 数值模拟展示出其方便、灵活的优势, 不仅能够实现长时间、大范围的沙丘流场研究, 还可以灵活调整计算参数, 降低研究的时间代价和成本, 尤其是在沙丘流场细节信息的获取上, 数值模拟成为一种有效途径.

1 国内外研究现状

近年来, 一些学者针对沙丘表面流场特性开展了一系列研究.2004年, Parsons等[13]借助商用软件PHOENICS, 利用RNG k-ε湍流模型对剖面为三角形的风成沙丘表面流场特性进行了二维数值模拟, 准确地预测到沙丘表面流型, 获得了气流在坡脚的停滞区、迎风坡的加速区和背风坡的流动分离和回流区.2005年, Herrmann等[14]采用商用软件FLUENT的k-ε模型对三维新月形沙丘和剖面近似为三角形的二维横向沙丘表面流场开展数值模拟, 展示了沙丘背风坡流动分离的现象.2010年, 江丽娟等[9]基于FLUENT的k-ε模型对新月形沙丘表面流场进行数值模拟, 所得结果能够反映剪应力的定性特征, 但在定量上与实验有一定差异.2012年, Balogh等[15]对湍流模型和壁面函数进行修正, 分别基于FLUENT和开源软件OpenFOAM研究了复杂地形大气边界层流场特性, 结果显示OpenFOAM具有更高的准确性和灵活性.2013年, 周晓斯等[10]基于大涡模拟(LES)方法, 自行编制计算程序, 对小尺度新月形沙丘背风坡流场特性进行二维数值模拟, 所得风速与风洞实验结果吻合良好.同年, Araújo等[11]采用FLUENT的k-ε模型, 对中剖面近似为三角形的横向沙丘表面流场开展数值模拟, 发现了回流区长度随剪切速度的变化关系.2014年, Peralta等[16]基于OpenFOAM的标准k-ε、RNG k-εk-ω SST模型, 采用结构化和非结构化两种网格类型, 模拟了两种复杂地形上大气边界层风场, 并将时均风速和湍动能的模拟结果与实验结果进行对比, 结果显示:第一种地形模拟结果和实验结果符合较好, 在采用RNG k-ε湍流模型和结构化网格时吻合最好; 而对于第二种地形, 模拟结果与实验结果偏差较大.2015年, Michelsen等[12]基于OpenFOAM的RNG k-ε湍流模型, 对新月形沙丘表面流场开展二维和三维模拟, 结果显示, 二维模拟低估了沙丘表面风速.

上述工作主要通过与野外观测和风洞实验数据的对比来验证模型, 进而对沙丘流场动力学开展模拟.然而, 不同的研究者所针对的地形和空间尺度差别很大, 采用的CFD软件及湍流模型也有所不同, 所得的结论缺乏普适性.因此, 有必要结合实验结果, 对沙丘表面流场特性开展更细致的研究, 以寻求适应性强、快速高效的数值模拟方法, 进而获得准确的沙丘表面流场信息.OpenFOAM因其免费、开源、稳定强大的底层类库, 以及能够并行计算的优点受到CFD用户的关注, 并在大气边界层模拟中发挥着日益重要的作用[12, 15-16].最近, Claudin等[4]首次报道了低矮正弦形沙丘表面流场的野外观测结果, 获得了贴近沙丘表面的风速分布, 为沙丘表面流场的精细数值模拟提供了新的实验依据.本文借助OpenFOAM的标准k-ε模型, 采用与文献[4]相同的沙丘轮廓, 对沙丘表面流场开展数值模拟, 将数值模拟结果与文献[4]数据相结合, 分析流场特性.在此基础上, 探讨了网格数目对计算结果的影响, 并对不同宽高比沙丘表面流场进行预测.本文研究可为进一步利用数值模拟方法研究沙丘表面沙粒输运提供基础和参考.

2 计算模型

沙丘表面流场为处于大气边界层内的近地面流场, 其特点是小流速、大雷诺数的湍流.这种湍流速度场往往具有扰动尺度小且变化频率高的特点, 在建立描述流场的纳维-斯托克斯方程时通常采用将速度分解为时均速度和脉动速度, 以降低求解的时间成本.将风速方向设定为x向, 竖直向上方向设定为z向, 与xz垂直的方向设定为y向, 则在笛卡尔坐标系下雷诺平均纳维-斯托克斯(RANS)方程可写为

(1)
(2)

式中:t为时间; u为流速; p为压力; ρ为空气密度; μ为空气动力粘度; δij为克罗内克函数, 当i=j时, δij=1, 当ij时, δij=0.

为使方程封闭, 引入标准k-ε湍流模型计算湍动能k和耗散率ε, 即

(3)
(4)

其中,

(5)

式中:σk, σε, C1ε, C2ε, Cμ为模型常数, 在大气边界层模拟中通常采用σk=1.0, σε=1.30, C1ε=1.21, C2ε=1.92, Cμ=0.03 [17-19]; Eij为平均应变率, 其表达式为

(6)
3 计算条件及方法 3.1 计算区域

为了验证计算模型和计算方法的准确性, 对文献[4]中沙丘表面流场进行数值模拟, 并与其野外观测数据进行对比.文献[4]所研究的沙丘位于平坦的地面上, 沙丘轮廓可以近似为正弦曲线

(7)

式中:zd为沙丘表面的z坐标; zref为参考高度; A为幅值; K为波数, K=2π/λ, λ为波长; x=0对应于沙丘坡顶位置, zd=0对应于沙丘坡脚位置.文献[4]给出, zref=A≈1.1 m, K≈0.14 m-1, λ≈46 m.此外, 为通过数值模拟考察沙丘宽高比对流场特性的影响, 计算时保持沙丘宽度(即波长λ)不变, 通过调整沙丘高度(2A)以改变沙丘的宽高比.

开展数值模拟首先应当确定计算区域, 并划分计算网格.由于沙丘表面速度变化集中在来流方向(x向)和沙丘高度方向(z向), 沙丘长度方向(y向)速度变化很小; 为减少网格数目, 缩短计算时间, 开展二维模拟.为避免边界效应, 进口、出口、顶面与沙丘的距离应足够大.为此, 选择计算区域宽度Δx=200 m, 计算区域高度根据沙丘坡顶高度进行调整, 保持两高度之比不发生变化, 其中沙丘幅值A=1.1 m时, Δz=20 m, 沙丘位于计算区域中心, 如图 1所示.


图 1 计算区域示意图 Fig. 1 Schematic of the simulation domain
3.2 网格划分

计算网格由OpenFOAM中的blockMesh和snappyHexMesh网格生成程序产生.在数值模拟中, 网格越小, 计算精度越高, 但计算耗时也越长.为解决计算精度和计算速度之间的矛盾, 本文分区域设置网格大小等级, 对于速度变化较大的沙丘附近区域采用较密的网格.以幅值A=1.1 m的沙丘表面流场计算区域为例, 图 2给出了不同等级网格区域分布, 在研究不同高度沙丘时, 各区域高度按比例调整.总的来说, 区域Ⅰ网格最大, 区域Ⅱ网格尺寸(长和高)为区域Ⅰ的1/2, 区域Ⅲ网格尺寸为区域Ⅱ的1/2.在计算区域底面(即包含沙丘的地面)添加5层边界层, 在远离下表面的方向上, 边界层网格厚度呈比例增加, 比例系数为1.1.


图 2 不同等级网格区域分布 Fig. 2 Distribution of regions with different mesh levels
3.3 边界条件

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

(8)

式中:κ为冯·卡门常数, κ=0.41;u*为剪切速度, u*=0.5 m/s; z0为空气动力学粗糙长度, z0=0.83 mm.

湍动能和耗散率的表达式为[20-21]

(9)

底面采用壁面边界条件; 顶面设置为移动壁面, 顶面风速由式(8) 确定; 出口设置为压力出口.此外, 对于本文的二维数值模拟, 需要在y向采用单层网格, 并将相应的边界设置为“empty”.

3.4 求解控制

基于SIMPLE算法[22-24]对压力和动量方程进行耦合求解, 该算法利用有限体积方法求解RANS方程, 由于在计算过程中标量储存在网格中心, 矢量储存在网格面中心, 因而需要采用插值方法.本文的数值计算中, 梯度项、拉普拉斯项采用高斯线性插值格式; 速度项采用高斯QUICKV插值格式; 湍动能项以及湍动能扩散项则采用高斯迎风格式.压力、速度、湍动能和耗散率的松弛因子分别为0.3, 0.7, 0.7和0.7.收敛条件为质量和动量通量的残差达到10-8, 压力的残差达到10-7.

4 结果与讨论 4.1 数值模拟结果的验证

为验证数学模型和数值计算方法的合理性, 将数值模拟得到的沙丘表面速度分布与文献[4]的野外观测结果作对比, 如图 3所示.其中, 计算区域网格数目为1.34×105个, 图 3(a)给出了距离沙丘表面0.11 m高度处(内层)的速度分布, 图 3(b)给出了距离沙丘表面0.50 m高度处(外层)的速度分布.图中, u0为最大速度和最小速度的平均值, 连续的曲线为数值模拟结果, 带误差条的离散点为来自文献[4]的野外观测结果.可见, 数值模拟和文献[4]结果吻合较好, 表明二维模拟能够预测正弦形沙丘表面流场特性.但也需要注意到, 在一些位置模拟值比文献[4]略低, 这可能是由于模拟时将三维沙丘视为二维所致, 已有对新月形沙丘表面流场的数值模拟结果显示, 二维模拟低估了沙丘表面风速[12], 后续研究可采用三维模拟进一步优化计算结果.图 3的数值模拟和文献[4]结果均反映出:内层气流在迎风坡加速, 在坡顶上游速度达到最大值, 速度分布相对沙丘轮廓具有相位提前量φ/K≈3.2 m; 外层气流速度相位提前量明显减小.内层和外层速度分布差异是由于流体力学机理不同所致[25-27].在内层, 流体粘性起重要作用, 由于粘性阻力的影响, 风速峰值出现在沙丘坡顶前; 随着高度的增加, 惯性作用增强, 粘性的影响减弱, 风速曲线的相位提前量减小.沙丘表面流场内层速度具有相位提前的现象, 这对沙丘的形成、移动以及演变具有重要影响.气流在迎风侧加速, 背风侧减速, 造成表面剪切应力在迎风侧增加, 背风侧减小, 从而沙丘上游产生吹蚀区, 下游产生沉积区, 进而引发沙丘的演变.


图 3 数值模拟得到的速度分布及与文献[4]的对比 Fig. 3 Velocity profile obtained from numerical simulation and the comparison with document four
4.2 沙丘表面流场速度分布

图 4给出了沙丘表面0.11, 0.50, 1.50, 2.00 m高度处的速度分布, 数值模拟采用的网格数目与图 3相同.可以看出, 速度分布随高度发生显著变化, 随着高度的增加, 最大速度出现的位置向坡顶迁移, 沙丘迎风坡和背风坡速度趋于对称分布, 速度变化变得平缓.


图 4 沙丘表面不同高度处的速度分布 Fig. 4 Velocity profiles at different heights over the dune

图 5给出了沙丘表面流场速度矢量图.由图可见, 受沙丘表面摩擦的影响, 贴近沙丘表面的区域速度较低, 且在迎风坡具有向上的分量, 在背风坡具有向下的分量; 在离沙丘表面较远的区域, 沙丘表面的影响减弱, 流场速度的竖直分量基本消失, 速度方向与来流方向平行.已有研究曾对沙丘背风坡流动分离的现象作过报道[11, 13-14, 28-29], 流动分离的形成与沙丘大小和外形、来流速度均有关系.在本节研究所针对的低矮沙丘及中等来流速度条件下, 无论野外观测还是数值模拟均未发现流动分离现象.由图 5还可以看出, 数值模拟可以全面展现流场丰富的细节信息, 是对野外观测的有效补充.


图 5 沙丘表面流场速度矢量图 Fig. 5 Velocity vectors of the flow field over the dune
4.3 网格数目对流场的影响

为考察数值模拟结果对计算网格的敏感性, 图 6给出了计算区域总网格数目不同时得到的沙丘表面流场速度分布.其中, 图 6(a), (b)分别为沙丘表面0.11, 0.50 m处的结果.对比可知, 网格数目对内层流场速度有重要影响, 而对外层流场速度的影响较小.内层速度分布曲线相对沙丘轮廓曲线具有相位提前量的特征已被野外观测所证实[4], 并能通过理论分析得到合理解释[30].然而, 在网格数目较少, 例如为4.10×104, 7.30×104个时, 数值模拟结果无法反映出速度相位提前这一重要特征.外层速度分布曲线虽然对网格数目敏感性较弱, 但沙丘坡顶附近以及背风坡流速随着网格数目的增加略有增大.同样反映出网格数目较少时, 无法体现速度分布的不对称性.图 6中4种网格数目情况下, 网格数目为1.34×105个时的数值模拟结果与野外观测结果吻合良好, 且继续增大网格数目对计算结果的影响很小, 是数值模拟中适宜的网格选择.该网格数目条件下, 区域Ⅰ, Ⅱ, Ⅲ (见图 2)中网格尺寸分别为0.2 m×0.2 m, 0.1 m×0.1 m, 0.05 m×0.05 m.本文采用的分区划分网格方法, 以及网格相对沙丘宽度和高度的尺寸, 可为正弦形表面流场特性数值模拟研究中的网格划分提供参考.


图 6 不同网格数目时的速度分布 Fig. 6 Velocity profiles for different mesh numbers
4.4 沙丘宽高比对流场的影响

为研究正弦形沙丘几何参数对流场特性的影响, 对不同宽高比的沙丘表面流场进行数值预测.在A=1.1, 1.65, 2.2, 2.75, 3.3 m, 即宽高比λ/(2A)=20.9, 13.9, 10.5, 8.36, 6.97情况下获得了沙丘表面流场信息.图 7给出了距沙丘表面0.11 m与0.50 m处的流场速度分布情况.由图可知, 在距离沙丘表面0.11 m处, 随着沙丘宽高比的减小, 一方面, 速度分布曲线相对于沙丘轮廓线的相位提前量增大, 迎风坡气流速度变化加剧; 另一方面, 当沙丘宽高比λ/(2A)≤10.5时, 背风坡气流速度出现负值, 表明气流流过坡顶后下沉形成回流区.沙丘迎风坡的气流加速和背风坡的回流区与沙丘侵蚀率有重要关联, 沙丘宽高比越小, 侵蚀越显著, 沙丘移动和演变越迅速.在距离沙丘表面0.5 m处, 速度分布曲线与沙丘轮廓曲线的相位差均较小, 尤其是在沙丘宽高比λ/(2A)≥10.5时, 相位差几乎为0;同时, 沙丘宽高比越小, 迎风坡气流加速越快.沙丘宽高比λ/(2A)≥10.5时, 在背风坡气流速度全为正; 而λ/(2A)≤8.36时, 可观察到速度负值, 且高宽比越小, 速度大小越偏离0.这些结果表明, 随着沙丘宽高比的减小, 回流区长度和强度均增大.不同尺度沙丘表面流场特性研究是后续利用气固两相流数值模拟方法考察沙粒输运, 进而认识和掌握沙丘的形成机制、发展规律及变化趋势的基础.


图 7 不同宽高比沙丘表面的速度分布 Fig. 7 Velocity profiles over dunes with different aspect ratios
5 结论

本文基于开源CFD软件OpenFOAM对正弦形沙丘表面流场进行了数值模拟, 并将数值模拟结果与文献[4]的野外观测结果相对比, 验证数学模型和计算方法的正确性.在此基础上, 对沙丘表面流场特性进行分析, 探讨了数值模拟结果的网格敏感性, 考察了沙丘宽高比对流场特性的影响, 得到以下结论:

a.基于OpenFOAM的标准k-ε模型能够得到低矮正弦形沙丘表面流场特性, 由其预测出的内层和外层速度分布与文献[4]中的野外观测结果吻合良好;

b.在贴近沙丘表面的内层, 速度分布曲线相对于沙丘轮廓线存在相位提前, 随着高度的增加, 相位提前量减少, 速度变化趋于平缓, 直至风速的竖直分量不再存在, 沙丘形状对流场的影响消失;

c.计算网格数目对外层速度分布影响较小, 但对内层的速度分布则有重要影响, 足够多的网格数目是获得精确的沙丘表面流场信息的必要条件, 在沙丘流场数值模拟研究中应当引起重视;

d.沙丘宽度不变时, 沙丘宽高比减小, 沙丘表面速度分布曲线的相位提前量增加, 背风坡出现回流区, 并且宽高比越小, 回流区的长度与强度越大.

参考文献
[1] PARTELI E J R, KROY K, TSOAR H, et al. Morphodynamic modeling of aeolian dunes:review and future plans[J]. The European Physical Journal Special Topics, 2014, 223(11): 2269–2283. DOI:10.1140/epjst/e2014-02263-2
[2] CHARRU F, ANDREOTTI B, CLAUDIN P. Sand ripples and dunes[J]. Annual Review of Fluid Mechanics, 2013, 45(1): 469–493. DOI:10.1146/annurev-fluid-011212-140806
[3] 吕萍, 董治宝. 沙丘动力学数值模拟研究方法综述(Ⅰ):沙丘周围流场模拟[J]. 中国沙漠, 2015, 35(3): 526–533. DOI:10.7522/j.issn.1000-694X.2014.00072
[4] CLAUDIN P, WIGGS G F S, ANDREOTTI B. Field evidence for the upwind velocity shift at the crest of low dunes[J]. Boundary-Layer Meteorology, 2013, 148(1): 195–206. DOI:10.1007/s10546-013-9804-3
[5] PELLETIER J D. Deviations from self-similarity in Barchan form and flux:the case of the Salton sea dunes, California[J]. Journal of Geophysical Research:Earth Surface, 2013, 118(4): 2406–2420. DOI:10.1002/2013JF002867
[6] 张正偲, 董治宝. 横向沙丘表面气流湍流特征野外观测[J]. 地理科学, 2015, 35(5): 652–657.
[7] CREYSSELS M, DUPONT P, MOCTAR A, et al. Saltating particles in a turbulent boundary layer:experiment and theory[J]. Journal of Fluid Mechanics, 2009, 625(12): 47–74.
[8] 贾攀, 王元, 张洋. 缩尺沙丘演化特性的风洞实验研究[J]. 西安交通大学学报, 2013, 47(5): 67–71. DOI:10.7652/xjtuxb201305012
[9] 江丽娟, 马高生. 新月形沙丘表面流场的数值模拟[J]. 兰州大学学报(自然科学版), 2010, 46(1): 48–53.
[10] 周晓斯, 王元, 李志强. 小尺度新月形沙丘背风侧流场特性的大涡模拟分析[J]. 西安交通大学学报, 2013, 47(7): 114–119. DOI:10.7652/xjtuxb201307021
[11] ARAÚJO A D, PARTELI E J R, PÖSCHEL T, et al. Numerical modeling of the wind flow over a transversedune[J]. Scientific Reports, 2013, 3: 2858. DOI:10.1038/srep02858
[12] MICHELSEN B, STROBL S, PARTELI E R J, et al. Two-dimensional airflow modeling under predicts the wind velocity over dunes[J]. Scientific Reports, 2015, 5: 16572. DOI:10.1038/srep16572
[13] PARSONS D R, WIGGS G F S, WALKER I J, et al. Numerical modelling of airflow over an idealised transverse dune[J]. Environmental Modelling & Software, 2004, 19(2): 153–162.
[14] HERRMANN H J, ANDRADE JR J S, SCHATZ V, et al. Calculation of the separation streamlines of barchans and transversedunes[J]. Physica A:Statistical Mechanics and Its Applications, 2005, 357(1): 44–49. DOI:10.1016/j.physa.2005.05.057
[15] BALOGH M, PARENTE A, BENOCCI C. RANS simulation of ABL flow over complex terrains applying an Enhanced k-ε model and wall function formulation:implementation and comparison for fluent and OpenFOAM[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2012(104/105/106): 360–368.
[16] PERALTA C, NUGUSSE H, KOKILAVANI S P, et al.Validation of thesimpleFoam (RANS) solver for the atmospheric boundary layer in complex terrain[C]//ITM Web of Conferences.London:EDP Sciences, 2014, 2:01002.
[17] HARGREAVES D M, WRIGHT N G. On the use of the k-ε model in commercial CFD software to model the neutral atmospheric boundary layer[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2007, 95(5): 355–369. DOI:10.1016/j.jweia.2006.08.002
[18] ZHANG X D.CFD simulation of neutral ABL flows[R].Roskilde:Technical University of Denmark, 2009. https://www.mendeley.com/research-papers/cfd-simulation-neutral-abl-flows/
[19] TAPIA X P.Modelling of wind flow over complex terrain using OpenFoam[D].Gävle:University of Gävle, 2009:35-42. https://www.mendeley.com/research-papers/modelling-wind-flow-complex-terrain-using-openfoam/
[20] BLOCKEN B, STATHOPOULOS T, CARMELIET J. CFD simulation of the atmospheric boundary layer:wall function problems[J]. Atmospheric Environment, 2007, 41(2): 238–252. DOI:10.1016/j.atmosenv.2006.08.019
[21] YAN B W, LI Q S, HE Y C, et al. RANS simulation of neutral atmospheric boundary layer flows over complex terrain by proper imposition of boundary conditions and modification on the k-ε model[J]. Environmental Fluid Mechanics, 2016, 16(1): 1–23. DOI:10.1007/s10652-015-9408-1
[22] VERSTEEG H K, MALALASEKERA W. An introduction to computational fluid dynamics:the finite volumemethod[M]. 2nd ed. Harlow: Pearson Education Limited, 2007.
[23] JASAK H.Error analysis and estimation for the finite volume method with applications to fluid flows[D].London:Imperial College, 1996. http://citeseerx.ist.psu.edu/showciting?cid=2508913
[24] 陶文铨. 数值传热学[M]. 2版. 西安: 西安交通大学出版社, 2001.
[25] WEAVER C M, WIGGS G F S. Field measurements of mean and turbulent airflow over abarchan sand dune[J]. Geomorphology, 2011, 128(1/2): 32–41.
[26] BRITTER R E, HUNT J C R, RICHARDS K J. Air flow over a two-dimensional hill:studies of velocity speed-up, roughness effects andturbulence[J]. Quarterly Journal of the Royal Meteorological Society, 1981, 107(451): 91–110. DOI:10.1002/qj.49710745106
[27] BADDOCK M C, WIGGS G F S, LIVINGSTONE I. A field study of mean and turbulent flow characteristics upwind, over and downwind of barchan dunes[J]. Earth Surface Processes and Landforms, 2011, 36(11): 1435–1448. DOI:10.1002/esp.v36.11
[28] WALKER I J, NICKLING W G. Dynamics of secondary airflow and sediment transport over and in the lee of transversedunes[J]. Progress in Physical Geography, 2002, 26(1): 47–75. DOI:10.1191/0309133302pp325ra
[29] GROH C, WIERSCHEM A, AKSEL N, et al. Barchan dunes in two dimensions:experimental tests for minimal models[J]. Physical Review E, 2008, 78(2): 021304. DOI:10.1103/PhysRevE.78.021304
[30] COLOMBINI M. Revisiting the linear theory of sand dune formation[J]. Journal of Fluid Mechanics, 2004, 502: 1–16. DOI:10.1017/S0022112003007201