郑德乾， 包延芳， 李亮， 方平治， 马文勇， 汤胜茗

(1.河南工业大学土木工程学院， 郑州 450001；2.中再巨灾风险管理股份有限公司, 重庆 400023; 3.中国气象局上海台风研究所， 上海 200030； 4.石家庄铁道大学土木工程学院， 石家庄 050043)

自然风流经山地时会形成复杂的风场，研究山地地形的风场特性对于山区建筑和桥梁结构抗风设计、风能利用、太阳能支架抗风设计等具有重要的指导意义。

Momomura等[1]对某山区地形进行了为时两年的地形风场实测，发现平均风速和湍流强度受地形影响很大。马如进等[2]通过现场实测对某U形峡谷桥的选址进行了研究，给出了平均风剖面、湍流度和脉动风功率谱的变化范围。DeBray[3]进行了最大坡度为15°山坡的风洞试验，给出了斜坡上下游不同高度处的风速模型，但其研究只针对二维模型；Bowen等[4]进行了更为细致的不同坡度山丘的风洞试验研究，结果发现坡度越陡，山顶风加速效应越大；段静等[5]通过风洞试验研究山体背风面的风场特性时发现背风面的坡度会影响山顶处的地形加速效应，坡度越小，山顶加速效应越明显；牛华伟等[6]采用风洞试验对不同坡度的山丘地形进行了风速和脉动风的研究，发现坡度越大，竖向风速分量极大值反而减小，山体表面脉动风速越大。

与现场实测和风洞试验方法相比，CFD(computational fluid dynamics)数值模拟方法便于参数分析。根据对N-S(Naiver-Stokes)方程处理方式的不同，CFD数值模拟方法主要有[7-13]：基于时间平均的RANS(Reynolds averaged Naiver-Stokes)方法和基于空间平均的大涡模拟(Large eddy simulation，LES)方法。王凯等[7]在研究近地表风场受坡度影响时，坡度比在0.3°以下对地形风场影响较小，坡度比在0.7°以上，再附涡更加明显。闫姝等[8]研究了山丘的绕流场分布，发现当山丘高宽比较大时，背风面的绕流稳定性较差，易发生流动分离；但由于以上研究所采用的RANS方法仅能得到时间平均值，所以仅研究了平均风，而未涉及脉动风特性。与RANS方法相比，LES大涡模拟方法基于空间平均，对大尺度涡直接求解，对小尺度涡则采用亚格子模型模拟，因而能够模拟瞬态流场，能够同时获取平均和脉动风速场信息[9]。通过山丘地形风场特性的实测结果与数值模拟结果的对比，发现与RANS方法相比，LES方法能够更加准确地模拟流场[10-11]。杨祥生等[12]在坡度对不同地形风场的数值模拟研究中发现，背风面坡度越大湍流强度越大，但该研究中未描述分析山体周围的具体绕流情况。许昌等[13]通过LBM (lattice Boltzmann method)和LES结合的方式计算分析了不同坡度山丘地形风场，发现坡度越大，速度衰减越大，背风面的流动分离越剧烈，尾迹区越长，能量耗散越慢。目前，在坡度对山丘地形风场特性、地形加速效应的影响，以及影响机理还有待进一步研究。

为此，以理想三维山丘为研究对象，采用基于空间平均的大涡模拟方法，入流脉动的合成采用基于自保持边界条件[14]的改进涡方法[15]。进行了15°、21.8°、30°和45°4种坡度三维山丘非定常绕流计算，在将坡度为21.8°的三维山丘LES计算结果与相应风洞试验结果的对比，验证了数值模拟方法及参数的有效性基础上，分析了坡度对山丘不同位置处的平均风剖面和脉动风剖面的影响，从时均和瞬态流场角度对不同坡度山丘尤其是背风面进行了流场机理分析；最后，结合各国规范[16]，研究分析了坡度对三维山丘地形风场加速效应的影响范围。

1 计算模型及参数设置

研究对象为余弦型三维对称山丘，模型示意图如图1所示，相应的数学表达式为[17]

(1)

式(1)中：z为离地高度；H为不同坡度模型的高度；L为迎风面或背风面水平长度，L=1 000 mm；山体坡度定义为tanα=H/L，考虑的4个坡度α=15°、21.8°、30°和45°，通过改变模型高度H来实现。

为捕捉山体周围的瞬态流场信息，CFD数值模拟计算采用基于空间平均的LES大涡模拟方法，模型几何缩尺比为1∶1 000。

数值模拟中，计算域大小取为23L(流向)×6L(展向)×6L(竖向)，阻塞比不超过3%。采用非均匀结构化网格进行离散，在壁面关心区域进行网格加密，如图2所示，考虑两种网格布置的影响，网格伸展率均为1.15，其中Mesh-1网格总数为976 500，最小网格尺寸为0.000 3L；Mesh-2网格总数为1 904 000，最小网格尺寸为0.000 1L，两套网格对应壁面y+≤5.5。

图1 山体形状示意图Fig.1 Schematic diagram of mountain shape

图2 网格划分示意图Fig.2 Sketches of the mesh discretization

入流面采用速度入口，根据三维山丘所处的B类地貌，采用基于自保持边界条件的涡方法[15]合成大涡模拟入流脉动，涡量值取95。出流面为压力出口，计算域两侧及顶部均设置为对称边界条件，山体表面及地面均设置为无滑移壁面。压力速度耦合采用SIMPLEC算法，空间离散格式采用具有二阶精度的Bounded central differencing格式；时间离散格式为二阶隐式，时间步长为0.004 s；亚格子模型选用Dynamic Smagorinsky-Lilly动态亚格子模型。

2 大涡模拟方法及参数验证

为说明所提数值模拟方法及参数的有效性，首先以坡度21.8°三维山丘为对象，分别将大涡模拟所得B类湍流边界层风场，以及三维山丘绕流风场与试验结果进行了对比分析。

2.1 湍流边界层风场大涡模拟结果验证

图3为不同网格分辨率情况下，大涡模拟所得三维山丘迎风面山脚前方水平距离2.5L位置处的平均风剖面、湍流度剖面，以及该位置距地面H高度处的脉动风速谱与试验风场[18]结果(图3中“EXP”)的对比。

U为流向平均风速；Uref为距离地面4H高度处来流流向平均风速； Iu为流向湍流度，Iu=σu/U，其中σu为流向脉动风速根方差， 表示其脉动风速；Lu为流向湍流积分尺度；f为频率；S(f)为脉动风速 自功率谱密度函数；σ2为脉动风速的方差；z/H为无量纲高度； fS(f)/σ2为脉动风速功率谱；U/Uref为流向平均风速剖面； fLu/Uref为湍流积分尺度图3 湍流边界层风场大涡模拟与试验结果比较Fig.3 Comparisons of the turbulent boundary wind field between LES and wind tunnel test

由图3可知：两种网格工况的平均风剖面和湍流度剖面，大涡模拟结果与风洞试验结果基本一致，进一步提高网格分辨率，对风剖面模拟结果精度的提高不明显。不同网格模型的脉动风速谱在低频处均能与Von Karman谱较好地吻合，而在高频处由于网格过滤效应出现了一定程度的衰减，但基本能够反映惯性子区的能量分布。可见，所提的大涡模拟方法及参数能较好地重现B类湍流边界层风场。

2.2 三维山丘绕流风场大涡模拟结果验证

图4为不同网格分辨率情况下，三维山丘y=0纵剖面不同位置处的平均风速剖面和脉动风速剖面与风洞试验结果[18]的对比。

可以看出大涡模拟所得不同位置处的流向和竖向平均风速剖面均与风洞试验有较好的一致性，两种分辨率网格结果无明显差异，说明在Mesh-1网格布置情况下，进一步增加网格分辨率对平均风剖面模拟精度的提高不明显。对于不同位置处的流向和竖向脉动风速剖面来说，与风洞试验结果的差异主要体现在x/H=2.5和3.75位置近壁区(其中x为水平坐标，H为山体高度)，其中Mesh-2网格的流向脉动风速值[图4(c)]略高于试验值，而Mesh-1网格的竖向脉动风速值[图4(d)]则略低于试验结果；x/H=-2.5位置远离地面处的竖向脉动风速值[图4(d)]与试验值相比，Mesh-1网格结果偏低而Mesh-2网格结果偏高。总体上，两套网格的大涡模拟结果均能够反映山丘周围及其尾流场的脉动风速变化趋势。

综上，所提大涡模拟方法及参数能够较好地重现来流紊流边界层风场，同时能够较有效地模拟三维山丘地形非定常绕流风场特性。综合大涡模拟精度和计算效率，坡度对三维山丘风场的影响分析均基于Mesh-1网格。

x/H为无量纲水平长度；W/Uref为竖向平均风速剖面； σu/Uref为流向脉动风速剖面；σw/Uref为竖向脉动风速剖面； W为竖向平均风速；σw为竖向脉动风速根方差图4 三维山丘模型风场大涡模拟与试验结果对比 (坡度21.8°)Fig.4 Comparisons of wind field over 3D hilly model between LES and experimental data(slope 21.8°)

3 坡度对地形风场的影响及机理分析

基于Mesh-1网格，采用所提的大涡模拟方法及参数，进行了4个坡度三维山丘的非定常绕流数值模拟，分析坡度对三维山丘风场的影响，并结合流场分析结果讨论其影响机理。

3.1 平均风速和脉动风速剖面比较

图5为不同坡度三维对称山丘的流向和竖向平均风速剖面和脉动风速剖面的对比。

对于流向平均风剖面[图5(a)]，总体上，不同坡度的流向平均风速剖面在迎风面区域的变化趋势基本一致；流向平均风速均随着离地高度的增加而增大，受山体的阻挡，迎风面山脚处产生了回流现象而形成驻涡，使得山脚处的流向平均风发生了突变，随着坡度的增大，驻涡尺寸逐渐增大，对流向平均风起到了一定的抑制作用，使得流向平均风逐渐减小。具体分析如下。

(1)沿流向的风加速效应在山顶近壁区位置达到最大，数值比入流面近壁区相应位置增大约1.12倍，随着坡度的增大，山顶近壁区的流向平均风呈现先减小后增大的趋势，45°流向平均风最大，较15°增大了约14%。

(2)在背风面处，以坡度30°山丘模拟结果为分界线，当坡度小于30°时，随着坡度的增大，旋涡尺寸逐渐增大，且逐渐远离壁面，流向平均风也随之逐渐减小；当坡度大于30°时，山脚处的回流现象使得旋涡靠近壁面且远离山顶，流向平均风呈现增大趋势，这也是45°山顶处风速增大的原因。在尾流区，随着流向距离的增加，不同坡度的变化趋势逐渐接近并趋于一致。

对于竖向平均风速剖面[图5(b)]，整体上不同坡度的结果无明显差异，只是迎风面山腰处和山顶近壁区域有较小差别，以及背风面及尾流区有一定突变。在迎风面处，不同坡度的竖向平均风速剖面均随高度呈现先增大后减小的趋势，在山腰位置近壁区达到最大，且随坡度的增大逐渐增大；受地形影响，背风面的流向平均风剖面表现出的流动分离现象更为剧烈，大致在x/L=0.44～0.50位置[图5(a)]发生流动分离并在近壁区出现明显的回流现象，在下游x/L=2.0位置处逐渐恢复正常。

从图5(c)、图5(d)所示脉动风速剖面的比较结果来看，总体上，脉动风速值均随着测点高度的增加有所减小，在山脚处，随着坡度的增大，驻涡尺寸逐渐增大，能量逐渐集中，对应的脉动风速逐渐增大；在背风面和尾流区，随着坡度的增大，湍流流动更为复杂，产生不同尺度的涡，可以明显观察到流向和竖向脉动风速剖面形状在较大高度范围内均存在较明显的变形，且脉动风速值也明显增加，其中在x/L=0.5位置处增幅最大。局部分析如下。

(1)对于流向脉动风速剖面[图5(c)]，在山顶近壁区，随着坡度的增大，山顶处的涡量逐渐丰富且靠近壁面，涡结构更加集中，流向脉动风速值逐渐增大且影响高度逐渐靠近壁面。在背风面山腰位置处，随着坡度的增加，45°背风面的尾流较宽，能量更为集中，从而使得45°山丘的变化趋势明显大于其他坡度结果，数值约为15°模型的19%。

(2)对于竖向脉动风速剖面[图5(d)]，在山顶近壁区，坡度45°时脉动风速值最大，数值约为其他坡度的1.05倍。在背风面及尾流区域0.6H高度范围内，随着坡度的增大，竖向脉动风速值均呈现增大的趋势，这是因为45°背风面涡道较宽，能量更为集中。

图5 不同坡度三维山丘风剖面对比Fig.5 Comparisons of wind profiles over 3D hills with different slopes

3.2 瞬态流场分析

图6和图7分别为不同坡度的三维对称山丘在4个相同时刻，在y=0纵剖面和z/H=0.5水平截面的瞬态涡量图。可以看出：整体上，不同坡度下三维山丘周围涡量分布较为丰富，正涡与负涡同时存在。受山丘的阻挡，气流在迎风面前方形成了不同尺寸的驻涡(图6)，随着涡旋沿山丘的逐渐演变，在山顶处(图6)和山丘两侧(图7)出现漩涡脱落现象。在背风面，受山丘三维效应的影响，随着坡度的增加，背风面的涡脱更加复杂，能量更为分散，近壁区的涡尺度也相对增大，一定数量的小尺度涡在尾流区混合后演变成为大尺度涡(图6、图7)。

从不同时刻y=0纵剖面的瞬态涡量分布情况来看(图6)，在山顶处，随着坡度的增加，旋涡尺度逐渐增大，脉动风速逐渐增大(图5)。在背风面山腰处，不同坡度的三维山丘背风面均形成了较明显的旋涡尾迹，大尺度涡也相对增加，随着坡度的增大，涡结构逐渐变宽，脉动风逐渐增大(图5)。

从不同时刻的z/H=0.5水平截面的瞬态涡量分布情况来看(图7)，随着坡度的增加，背风面的大尺度涡更加集中，脉动风速值逐渐增大(图5)，旋涡尾迹逐渐变宽，漩涡发展距离变长，能量更强。在背风面山脚位置及尾流区域，随着坡度的增加，漩涡脱落强度逐渐增强，旋涡尾迹逐渐变宽，尾迹区更长。大坡度山丘背风面的能量变化更为复杂。

图6 不同坡度的y=0纵剖面瞬态涡演化规律比较Fig.6 Comparisons of the transient vortex evolution law at y=0 longitudinal section of the hill with different slopes

图7 不同坡度山丘z/H=0.5水平截面瞬态涡量演变比较Fig.7 Comparisons of the transient vortex evolution at z/H=0.5 horizontal section of the hill with different slopes

图8 不同坡度山丘y=0纵剖面时均流线图Fig.8 Time averaged streamline of y=0 longitudinal section of the hill with different slopes

3.3 时均流场分析

图8为不同坡度的三维对称山丘流向y=0纵剖面时均流线图，可以看出：对于不同坡度的三维对称山丘，在迎风面均产生了驻涡，随着坡度的增大，驻涡尺寸逐渐变大，这是由于受三维山丘的影响产生了回流现象且坡度越高阻挡越显著所致。在山丘顶部，随着坡度的增大，风速逐渐增大，导致流动分离点逐渐下移，在山丘顶部发生流动分离后，在背风面均产生了较大尺度的尾涡并附有漩涡脱落的现象产生。

随着坡度的增大，漩涡尺寸逐渐增大，附着涡范围也逐渐增大，4个坡度山丘的涡核心坐标分别为：15°山丘(x/L=0.40，z/L=0.80)、21.8°山丘(x/L=0.50，z/H=0.30)、30°山丘(x/L=0.62，z/L=0.48)、45°山丘(x/L=0.60，z/L=0.70)。可见，随着坡度的增大，涡核心逐渐下移，其中45°的涡核心更靠近下游，这是山顶处近壁区风速较大(图5)所致。

4 不同坡度的地形加速效应

地形加速效应是山地地形风场的重要指标，对比分析不同坡度三维对称山丘的迎风面山腰(x/L=-0.5)，山顶(x/L=0)，以及背风面山腰(x/L=0.5)3个位置处的地形加速效应对比。地形加速效应可表示为

(2)

式(2)中：Ui(x,z′)为数值模拟所得的平均风速；U0(z′)为对应高度处来流平均风速。

(3)

式(3)中：σui(x,z′)为大涡模拟所得脉动风速根方差。

图9为4种坡度三维山丘的地形加速效应随距离地面高度的变化曲线，可以看出：在迎风面山腰位置处[图9(a)]，当z/H<0.5，地形加速效应整体上随着坡度的增加呈逐渐减小趋势，在z/H=0.7位置处，15°模型的地形加速效应最大(值为0.98)，30°模型的地形加速效应最小(值为0.96)；随着测点高度的增加，不同坡度的地形加速效应逐渐趋于一致，数值趋近于1.0；不同坡度的脉动风速加速效应在迎风面的变化趋势基本一致，数值在0.01左右，随着坡度的增大，呈现减小的趋势，不同坡度之间差值最大在0.004左右。

在山顶位置处[图9(b)]，不同坡度的地形加速效应均随测点高度的增加呈现先增大后减小的趋势，且存在明显的拐点，均在z/H=0.05左右时地形加速效应达到最大值，其中45°数值最大(值为1.40)。不同坡度的脉动风速加速效应数值均接近0.1，以坡度30°为界，当坡度小于30°时，随着坡度的增大，不同坡度脉动风的地形加速效应变化逐渐减小；坡度大于30°时趋势则相反，当z/H>1.0时，坡度45°的脉动风速加速效应最大，但不同坡度的最大差值在0.04左右。

在背风面[图9(c)]，整体上，不同坡度在背风面山腰处平均风的地形加速效应变化趋势一致，随着测点高度的增加逐渐趋向于1.0。在背风面近壁区山腰位置处，当z/H<0.7时，以坡度30°为界，当坡度小于30°时，地形加速效应值随坡度的增大逐渐减小；当坡度大于30°时则相反。背风面的脉动风速加速效应随坡度的变化较为复杂，随着测点距离地面高度的增加呈现先减小后增大的趋势，其中45°的数值最大，较其他坡度增大约7倍。

由以上分析可以看出，在迎风面、山顶和背风面处，坡度对地形加速效应的最大影响高度分别在z/H=0.8、0.1、0.5拐点高度处，为此将该位置的迎风面、山顶和背风面地形加速效应数值模拟结果与不同国家的规范值[16]进行对比分析，如图10所示，限于篇幅，仅给出了《中国建筑结构荷载规范》(GB 50009—2012)(简称中国规范)、StructuralDesignActionsPart2：WindActions(AS/NZS 1170.2—2011)(简称澳大利亚规范)、AIJRecommendationsforLoadsonBuildings(AIJ—2004)(简称日本规范)，山地修正系数值，详细的计算方法见文献[16]。由图10可知：在迎风面[图10(a)]，地形加速效应大涡模拟结果随坡度的增加逐渐增大，而各国规范值随着坡度的增加逐渐趋于定值，各国规范的取值范围在1.0～3.0，《中国建筑结构荷载规范》 (GB 50009—2012)取值最大。与规范相比，不同坡度的大涡模拟结果均小于规范取值，当坡度小于30°时，大涡模拟结果与澳大利亚规范值更为接近；当坡度大于30°时，大涡模拟结果更接近日本规范值。

在山顶位置[图10(b)]处，随着坡度的增加，大涡模拟结果呈现先减小后增大的趋势，总体上与日本规范更为接近。当坡度α=15°时，大涡模拟结果为1.26，稍高于日本规范值1.25；当坡度大于15°时，其他规范值均较大，其中中国规范值最大(值为3.62)。

在背风面位置[图10(c)]处，随着坡度的增大，大涡模拟结果呈现先减小后增大的趋势，整体均较小于规范取值，其中坡度15°和45°大涡模拟结果与澳大利亚规范取值较为接近，坡度21.8°和30°结果均小于各国规范取值。

由以上分析可见，在迎风面、山顶和背风面处，中国、澳大利亚、日本规范取值离散性较大，随坡度变化规律不尽一致，规范值均较为保守，本文结果可为后续地形加速效应值的修正提供参考。

图10 大涡模拟所得地形加速效应最大影响高度处地形 加速效应随坡度变化曲线与规范结果比较Fig.10 Comparison of LES results for terrain acceleration effects at different slopes with different national standard values at maximum impact altitude

5 结论

(1)坡度对三维山丘绕流风场的影响主要集中在山顶近壁区处以及背风面和尾流区一定高度范围内。当坡度小于30°时，山顶位置处的流向平均风和脉动风速随坡度的增加呈减小趋势，当坡度大于30°时则呈增大趋势；背风面近壁区的流向平均风速有类似现象，但脉动风速在近壁区一定高度范围内随坡度增大呈现增大趋势。

(2)从流场角度来看，随着坡度的增加，山顶处旋涡尺度逐渐减小，背风面流动分离点逐渐向下游移动，引起涡旋中心位置逐渐向下游移动且远离山丘壁面，其脉动风逐渐增大，这是地面摩擦效应所致。在背风面近壁区处，当坡度小于30°时，随着坡度的增大其旋涡尺度不断增大，且旋涡尾迹逐渐变宽，旋涡脱落强度逐渐增强，脉动风速逐渐减小；坡度大于30°时则相反。

(3)不同坡度的地形加速效应趋势基本类似，均随着高度的增加逐渐趋近于1.0，主要差异体现在近壁区，以坡度30°为界，当坡度小于30°时，迎风面、山顶和背风面处的地形加速效应均随坡度的增加逐渐减小，而坡度大于30°时则呈相反趋势；在地面附近地形加速效应值较大的位置，各国规范取值离散性较大且取值偏保守。