考虑粗糙地表的近地层风场仿真模拟研究及误差分析

2020-11-04刘震卿吴晓波李秋明

空气动力学学报 2020年4期

刘震卿,张冲,吴晓波,李秋明

(1.华中科技大学土木工程与力学学院,武汉 430074;2.中建国际投资(河南)有限公司,郑州 450000)

0 引言

微地形上空紊流流场预测应用于许多重要工程,例如结构安全、农业风害、航空安全等,在绘制复杂地形风能分布中也发挥重要作用,其决定风力机排布位置。风资源在一些多山地区丰富,但风速分布复杂,因此精确预测复杂地形的风速和紊流强度非常重要。地表粗糙度是反映近地层动力特征的一个重要指标,是研究风场竖向剖面分布的重要参数,因而地表粗糙度对复杂地形风场研究具有重要意义。实际地表大多覆盖植被或者村落城镇,若需准确预测这些区域的风能分布情况,首先需要模拟地表粗糙对流场的影响[1-4]。

地表粗糙度长度是边界层气象学中常用的一个重要参数,是指在边界层大气中,近地层风速向下递减到零时的高度。根据Pattanapol等[5]对复杂地表粗糙度模拟方法的总结:对于复杂地表粗糙度的模拟方法包括两种,一种方法是基于湍流模型中的壁面函数法;另一种方法是在湍流模型输送方程中添加源项加以考虑。

对于壁面函数法,严格说来当且仅当计算区域入口处的来流边界条件与湍流模型及地表粗糙特性相协调时,才能得到水平均匀性的平衡大气边界层,但由于湍流的复杂性,在数值模拟中精确实现平衡大气边界层常存在困难。Richards和Hoxey[6]基于k-ε湍流模型提出了满足平衡大气边界层的几个原则,给出了一组来流边界条件,并对湍流模型参数的取值进行了探讨,但在其研究中所给出的湍动能边界条件为一常数,这与实测以及试验结果不符。Yang等[7]从标准k-ε 湍流模型方程出发,假设湍流在整个边界层内处于平衡态,推导出一组随高度变化的来流边界条件,并考虑壁面粗糙度,实现了边界层流场的平衡,但仅满足湍流能量在边界层内层范围才处于平衡。此外使用壁面函数法模拟地表粗糙度的还有Daniel[8]等、唐煜等[9]、邓院昌等[10]。对于在湍流模型输送方程中添加源项的粗糙度模拟方法,胡朋等[11]提出通过在湍动能输送方程与耗散率输送方程中添加源项的方法,使来流边界条件与湍流模型相协调,但SST k-ω 湍流模型中各系数很难确定。

另外Fluent提供了多种湍流模型,包括k-ε模型、k-ω 模型和LES等等。但是不同的湍流模型对于求解的适用性也不同,湍流模型的选择同样对风场数值模拟结果的精度和资源耗费有着十分重要的影响。

本文基于Fluent流体计算平台,使用自编UDF程序模拟精细化地表粗糙度,该程序实现了根据由Wind Pro中下载得到的目标区域地表粗糙度长度文件自动添加覆盖植被阻力源项来模拟精细化地表粗糙度,考虑了计算域内湍流结构信息以及覆盖植被类型、叶面积密度等相关参数的影响。以Askervein山为例建立数值模型选择了不同的湍流模型进行风场模拟,与实测数据对比,研究了湍流模型对复杂地形风场CFD 模拟的影响,分析了湍流模型的适用性,并进一步探讨了单一地表粗糙度长度的影响。

1 Askervein山实测试验

Askervein项目[12]是在国际能源署主持下开展的一项关于低山丘陵大气边界层的合作研究。Askervein位于苏格兰外赫布里底群岛南端的南尤伊斯特的西海岸,该山经纬度坐标为(57°11′N,7°22′W),相对高度为116 m,海拔高度126 m。Askervein山体示意图如图1 所示,A 线与B 线相交于Askervein山的最高点(HT-Hilltop),AA 线与B 线相交于中心点(CP-Centre Point)。此外还在距离此处3 km 的达利堡附近设置了一个参考点(RSReference Site),用于对与Askervein相遇之前未受干扰的流体运动进行详细测量。

图1 Askervein山体示意图Fig.1 Contour map of Askervein Mountain

在试验过程中,沿着A、B 和AA 三条线布置了50座以上的测风塔,大多数测风塔高10 m,还有多座更高的测风塔。其中记录代号为TU03-B 的测风数据发生在1983年10月3日的14:00～17:00风流稳定,风向约210°,由于该山实测数据丰富、便于将数值模拟结果与实测值对比,故本文以Askervein山为研究对象建立地形模型。

PA Taylor等[12]表示,除了附近的一些建筑物以及位于海岸线附近的沙丘之外,Askervein山和周围地区的地表粗糙度长度基本一致。地表粗糙度长度的初始估测值在0.01～0.05 m,根据参考点RS处的风剖面测量数据表明,地表粗糙度长度的范围大多在0.01～0.03m 的范围内。若作为单一的代表值,建议z0取值0.03 m。

2 数值模型与方法

2.1 湍流模型

数值模拟方法主要分为三种[13],包括直接数值模拟DNS、雷诺平均法RANS和大涡模拟LES。作为最准确的湍流模拟方式,DNS 可以直接求解三维N-S方程组,而不需要求平均或近似值来估计和控制误差,但其对网格精度要求很高,计算量很大,比较耗时;雷诺平均法不需要计算各尺度的湍流脉动,降低了计算量,但无法解析具有较强脉动特征的小尺度涡;大涡模拟直接模拟大尺度紊流运动,利用次网格尺度模型模拟小尺度紊流运动对大尺度紊流运动的影响,但较雷诺平均法更耗时。所以本文选取了后两种数值模拟方法——雷诺平均和大涡模拟,因而选取了基于雷诺平均法而建立的k-ε 模型、RNG k-ε 模型以及基于大涡模拟的LES模型。

(1)标准k-ε 模型

二维不可压缩流场质量方程和动量方程为:

式中,ui为平均速度,为雷诺应力,μ 为黏性系数,ρ 为空气密度,p 为压强,代表粗糙度阻力源项。湍动能k 和湍动耗散率ε 的运输方程分别如下式,其中湍流脉动能k 方程为(3)式,湍流耗散率ε 方程为(4)式:

平均速度梯度引起湍动能k 产生了Gk项;标准k-ε湍流模型中的经验常数C1ε、C2ε、σk、σε一般取值如下:C1ε＝1.44、C2ε＝1.92、σk＝1.0,σε＝1.3,μt为湍流黏度,表示为μt ＝ρCμk2/ε。Cμ为经验常数Cμ＝0.09。式中fk和fε为阻力源项,以模拟地表粗糙度对风场的影响。

(2)RNG k-ε 模型

RNG 模型考虑了湍流漩涡的影响,提供了一个考虑低雷诺数流动黏性的解析公式,这些均使得RNG k-ε 模型有更高的计算精度和稳定性,其质量方程和动量方程如(1)、(2)式,湍动能方程和耗散率方程如下所示:

式中,C1ε、C2ε、C3ε为经验常数,一般取值C1ε＝1.42、C2ε＝1.68、C3ε＝1.72,αk＝αε＝1.393,μeff＝μ+μt,Gb为由浮动引起的湍动能。Rε代表平均应变率对ε影响的附加项。式中fk和fε为阻力源项,以模拟地表粗糙度对风场的影响。

(3)LES模型

湍流运动可以看做是由很多尺度大小不同的涡组成的,大尺度涡比小尺度涡运动更剧烈,具有各向异性的特点,对于湍流的平均运动影响更明显,可以实现动量、质量和能量等的交换;小尺度涡相对通常较弱,主要通过非线性的作用对大尺度涡的运动产生影响,作用主要表现为耗散[13],因为共性更多,近似于各向同性。不可压常黏性系数的连续性方程和N-S方程为:

式中,Sij为拉伸率张量,Sij＝(∂ui/∂xj+∂ui/∂xj)/2;ρ为为流体密度;γ 为分子黏性系数。滤波后的动量方程如(10式),其中f~u,i为考虑地表粗糙度阻力源项。

亚格子模型采用Smagorinsky-Lilly,亚格子应力(SGS)的表达式如下:

式中,Ls—亚格子的混合长度;κ—卡门常数,κ＝0.42;d—壁面网格到壁面的距离;V—控制体体积。

三种湍流模型使用Fluent求解三维非稳态Navier-Stokes方程,空间离散方式采用有限体积法,二阶中心差分格式用于对流项与黏性项,二阶隐式格式用于非稳态项的时间推进;压力的插值格式将选取标准格式;SIMPLE算法用于压强速度解耦。

2.2 粗糙度模拟

本文采用在计算域内添加阻力源项的方法来模拟地表粗糙度。为了考虑地表粗糙区域对流场的阻碍作用,通过在Navier-Stokes动量方程中加入阻力源项加以考虑,在阻力源项的计算中,需要首先确定粗糙元的阻力系数,针对此参数已开展过相关试验研究,并给出了经验取值。由于Askervein山地表粗糙元为覆盖植被,根据Enoki和Ishihara[10]的研究,当粗糙区域粗糙元为覆盖植被时,采用表达式:

上式中,z 为距离地表高度,h 为植被高度。对于不同地表覆盖情况,粗糙度长度z0与覆盖植被高度h满足h＝a×z0,a 为两种粗糙度表征参数的转换系数。本文根据既有研究,选取转换系数a＝20.0[12],抗力系数CD,t＝0.4[10]。

对于k-ε 模型湍动能和湍动耗散率的运输方程中的阻力源项fk和fε由Enoki与Ishihara[10]提出的阻力源项计算方法求得:

式中,u 为速度大小,βp＝1.0,βd＝4.0,Cpε1＝1.5,Cpε2＝0.6。

2.3 网格模型

考虑到后续边界条件的设置以及最高点HT 的高度,计算域长宽高设置为6 km×6 km×1.5 km,底面中心点为Askervein山的最高点HT 如图2(b)所示,为使地形平稳过渡到光滑地区,设置宽0.4 km 的过渡区圆环如图2(a)。建立网格模型时,参考Askervein山的尺寸,以最高点HT 所在垂线为中心,半径1 km 的圆柱体区域内进行网格细化加密,网格分辨率为15 m 如图2(c),以精确显示此地的地形地势。山体周围逐渐过渡到平坦地形,设置网格尺寸平稳的从加密区过渡到边界处,最大网格分辨率为80 m。垂直方向上在近地面进行网格加密最小网格尺寸0.05 m,以准确捕捉近地面流场,采用相邻网格尺寸比值为一定值的σ 网格;为减少由于网格尺寸急剧变化而带来的数值误差,计算域内网格最大生长率为1.2;网格总量约为200万。采用非结构化三棱柱网格,以充分拟合几何边界,同时也可以保证每层单元拓扑结构的一致性。图2为CFD 三维实际地形模型示意图。

图2 CFD三维实际地形模型及网格分布示意图Fig.2 CFD three-dimensional actual terrain model and grid distribution diagram

2.4 边界条件

计算域左侧采用速度入口,右侧采用压力出口,四周壁面采用对称边界条件(∂u/∂n＝0,∂w/∂n＝0,v＝0),底面采用无滑移边界条件(u＝0,v＝0,w＝0,∂p/∂n＝0)。根据Richards等[6]的大气边界层理论,入口处速度选择参考点RS处210°风向角时测风塔的风速,采用以下模型:

式中,u*为摩擦速度,κ＝0.4[6]为卡门常数,z0为地表粗糙度长度。根据参考点RS 测风数据得到:u*＝0.618 m/s,z0＝0.03[16]。使用标准k-ε 模型和RNG k-ε 模型时,入口边界的湍动能、耗散率依次设置为:

式中,hg＝891.82 m,当高度超过界限值后,湍动能趋于零,取值0.0001。

董喜阳，1986年生于吉林九台。文学硕士。诗人、作家，兼事文学、美术评论。鲁迅文学院第三十四届中青年作家高级研讨班（青年作家班）学员。中国作协会员。作品散见于《诗刊》《人民日报》《大家》《诗选刊》等刊物。现居长春。

3 湍流模型分析

3.1 风速比

Askervein试验包括不同的来流风速,为了便于分析数值模拟的计算结果,在此便不能使用实际风速进行比较,通常使用一个无量纲参数风速比S 来描述风加速效应。

上式中,U(z)表示距离地面高度z 处的风速,U0(z)为入口处距离地面z 高度处的风速,取z＝10 m。定义风速比标准偏差衡量计算结果与实测结果的偏差,风速比标准偏差由下式计算:

式中,n 为测风点个数,Sei表示第i 个测风点距离地面10 m 高度的风速实测值与入口基准点RS相同高度的风速实测值的风速比,Sai表示第i 个测风点离地面10 m 高度的CFD 数值模拟风速值与入口基准点RS相同高度的风速实测值的风速比。

3.2 精细化粗糙地表模拟

从WindPro中下载包含Askervein山的计算域范围内的地表粗糙度长度文件,通过自编程序插值得到计算域内网格节点处粗糙度长度。图3显示了计算域内地表粗糙度长度分布情况,可以看出地表粗糙度长度分布在0.01～0.05 m 范围内。

图3 计算域内地表粗糙度长度分布图(单位:m)Fig.3 Roughness length distribution in computational domain(unit:m)

3.3 10 m 高度处风速比

基于精细化地表粗糙度模拟,选用三种不同湍流模型进行风场模拟。沿AA 线和B 线布置的10 m高测风塔,均可以获得比较全面的数据,由于A 线和AA 线之间具有几何地形的相似性,因此仅对贯穿山体中心的AA 线和B 线进行10 m 高度处风速比分析。图4、5分别显示了沿AA 线和B 线10 m 高度处的风速比与实测数据[17]的对比,图中柱状图分别代表各湍流模型相比实测值的风速比标准偏差。

图4 沿A线10 m 高度的风速比Fig.4 Wind speed ratio at 10 m along line A

从图4中可以看出,沿AA 线10 m 高度处风速在遇到山体阻碍后风速降低,到达山顶产生较大的风加速,随后在背风面处风速迅速降低。三种湍流模型风速比计算结果在迎风面趋势一致,均与实测值吻合较好,其中LES模型吻合更好;在背风面则显示出一定误差,主要由于在背风面发生流动分离,导致误差产生,但LES模型可以捕捉背风面风速迅速下降的现象。整体上LES计算结果标准偏差为0.0564,相比于标准k-ε 模型结果0.0733 和RNG k-ε 模型0.0681误差均要小,精度更高。

图5对比了三种湍流模型计算结果沿B线10 m高度处风速比和实测数据。B线为山脊线,可以发现随着山脊升高,风速比也随之增大,这是由来流方向垂直于B线,山脊越高风加速效应越明显,导致风速比越大。整体上LES计算结果相比于实测值偏差最小为0.0758。而标准k-ε 模型和RNG k-ε 模型结果偏差分别为0.1412、0.1217,两者计算结果差别不大,但误差均大于LES模型。

图5 沿B线10 m 高度的风速比Fig.5 Wind speed ratio at 10 m along line B

表1 三种湍流模型风速比标准偏差分析Table 1 Analysis of wind speed ratio standard deviation under

3.4 风剖面

三种湍流模型在点HT 处风剖面风速比计算结果对比如图6所示。从图中可以看出三种湍流模型在HT 处的风剖面风速比曲线与实测值基本一致,但是在近地面,LES的峰值大于其余两种湍流模型,更接近实测值,误差更小;标准k-ε 模型和RNG k-ε模型在5～20 m 高度时接近风速比实测值,但在高处风速计算结果明显偏大。LES模型计算结果与实测值偏差11.7%,相比于标准k-ε 模型和RNG k-ε 模型偏差分别下降了7.7%和2.7%,模拟误差更小。

图6 HT处的垂直风剖面风速比Fig.6 Vertical wind profile wind speed ratio at HT

3.5 湍动能

Askervein山布置的测风塔中沿A 线分布的测风塔主要用于采集湍流数据,故将三种湍流模型模拟得到沿A 线10 m 高度处的湍动能分布于实测结果对比,如图7 所示,其中湍动能标准偏差定义同式(17)。由图可以看出湍动能在山顶减小,随后在距离山顶大约200 m 的背风面附近达到极大值。三种湍流模型计算结果中LES与实测结果吻合最好,标准偏差为0.4191,相比于标准k-ε 模型和RNG k-ε 模型0.9376、0.8860均要小。

图7 A线10 m 高度处湍动能分布Fig.7 Turbulent energy of line A at 10 m hight

图8 HT处垂直方向的湍动能Fig.8 Vertical turbulent energy at HT

三种湍流模型在点HT 处垂直方向上湍动能分布如图8所示。图8显示,虽然三者的趋势基本一致,但是标准k-ε 和RNG k-ε 模型计算结果明显比实测值大,由于标准k-ε 和RNG k-ε 无法解析具有较强脉动特征的小尺度涡,尤其在近地面处,表明在此处并不适用,而LES模型能较好捕捉近地面的小尺度涡进而模拟值与实测值基本吻合,误差更小。LES、RNG k-ε 和标准k-ε 对点HT 垂直方向的湍动能模拟结果的标准差分别为0.3951、1.3467和0.8702,相比之下LES更为精确。

4 单一地表粗糙度

上一小节讨论了基于精细化地表粗糙度模拟下不同湍流模型的性能,结果表明,相比于其他几种湍流模型,LES模型计算结果精度更高。但通常难以获得较精确的精细化地表粗糙度分布,本节为了探明模拟复杂地形风场时,单一地表粗糙度能否获得理想结果,选择LES湍流模型,设置地表粗糙度长度z0分别为0.00,0.01,0.02,0.03,0.04,0.05,0.06,工况设置如表2所示。

模拟不同工况下的风场,将距离地面10 m 高度的风速转换为风速比,与实测值[17]进行比较。图9、图10和图11分别显示了沿A 线、AA 线和B线距地10 m 高度处模拟风速比与实测数据的对比,图中柱状图分别代表各湍流模型相比于实测值的风速比标准偏差。

从图9～图11可以看出,地表粗糙度长度对实际地形的CFD 数值模拟结果影响比较大。当地表粗糙度长度设置为0.00时,即采用光滑地表不考虑粗糙度长度对风场的阻碍作用,在A 线、AA 线的所有测风点以及B线大部分测风点处的数值模拟结果明显远大于实测风速值,同时也大于各个粗糙度工况的模拟结果。对于粗糙度长度的设置,从图9～图11均可看出随着地表粗糙度长度的增大,各测点的数值模拟风速值逐渐降低。当地表粗糙度长度设置为0.06 m 和0.05 m 时,整体的风速模拟值远小于实测值,标准偏差较大。而当地表粗糙度长度设置为0.02～0.04 m 时,整体的风速比曲线与实测值的吻合度较高,表明该地形的地表粗糙度长度的平均值接近0.03 m,与实际地表植被平均分布情况一致[17]。

表2 工况汇总表Table 2 Summary of cases

图9 沿A线10 m 高度的风速比Fig.9 Wind speed ratio at 10 m along line A

图10 沿AA线10 m 高度的风速比Fig.10 Wind speed ratio at 10 m along line AA

图11 沿B线10 m 高度的风速比Fig.11 Wind speed ratio at 10 m along line B

另外从图9和图10中,可以看到在迎风面山坡,地表粗糙度长度为0.02 m 时风速比与实测值吻合较好;在背风面山坡,地表粗糙度长度为0.04 m 时吻合较好,表明该地区地表粗糙度长度分布的差异性。

同时与采用单一地表粗糙度长度工况结果对比,统计各工况所有测点风速比相比于实测值的标准偏差如表3所示。由表3可以看出,当地表粗糙度长度设置为0.02～0.03时,模拟风速比误差为0.0663～0.0784相比于精细化地表粗糙下模拟误差0.0620十分接近,而当单一地表粗糙度长度小于0.02或大于0.03时,模拟风速比误差则过大,证明了模拟复杂地形风场时合理设置单一地表粗糙度能有效减小模拟误差,满足精度要求。