0 引言

桥位风特性参数是大跨度桥梁抗风设计的基础^[1]。对于开阔江面或河道，其桥位风特性参数已有大量实测资料和研究成果可供参考，对应桥位风特性参数比较容易确定。对于复杂地形桥位，由于桥位地形差别较大，其风特性研究难度相对较大。目前复杂地形桥位风特性研究主要有现场实测、风洞试验和CFD数值模拟。现场实测是目前桥位风特性参数研究最为可靠的方法^{[2, 3]}，但现场实测需要耗费大量的人力物力和时间，且仅能获得某几个固定位置的风特性参数，而对于跨江河、峡谷大型桥梁则无法直接得到河中央桥面处风速特性。风洞试验通过地形缩尺模型来研究桥位处风特性参数，也是桥位风特性参数研究的重要方法之一，但地形模型风洞试验方法试验周期较长、试验经费较高^{[4, 5]}。近年来复杂地形桥位风特性CFD数值模拟越来越受到关注。Yamaguchi和Ishihara在进行风洞试验同时采用数值模拟方法某区域风场进行了研究，数值模拟结果与试验结果吻合较好^[6]。Uchida和Ohya采用大涡模拟方法对9.5 km×5 km区域范围内空气流动进行了数值模拟，使用粗糙方块自造脉动风，获得计算区域内平均风速与脉动风速情况，并研究了地形因素对风场影响^[7]。Bitsuamlak采用面向对象编程语言C++开发了CFD程序来求解雷诺应力时均N-S方程，采用k－ε模型分别对不同地形风特性进行了数值模拟，数值模拟结果与已有风洞试验结果吻合较好^[8]。肖仪清等人针对我国南部地区某岛屿三维地形模型进行了CFD数值模拟，实际地形范围为10.1 km×6.9 km，竖直方向取2 000 m，分别采用40 m×40 m和20 m×20 m两种大小网格，分别采用RNG k－ε模型和SST k－ω模型进行计算。研究表明：SST k－ω模型计算结果优于RNG k－ε模型计算结果；20 m×20 m网格计算结果优于40 m×40 m网格计算结果^[9]。祝志文、张士宁等人针对乌江大桥桥位地形进行了CFD数值模拟，实际地形范围为10 km×6 km，竖直方向取3 000 m，近地面采用三角形网格自由划分，第一层网格距离地面5 m，网格总数约为110万，采用Standard k－ε模型进行计算，研究给出了桥址附近风场特征以及相关参数^[10]。周志勇等以1∶1 000比例建立了实际大小为23 km×27 km 的某山区三维地形模型，网格划分方式有120 m×120 m，80 m×80 m，60 m×60 m等3种形式，网格数量分别为110万、400万、520万。采用Realizable k－ε湍流模型，分别设置了13个监测点进行数据监测。计算表明：采用多面体网格划分方式在较少网格数量下也可以得到较高计算精度，可较显著地节省计算网格数和内存资源^[11]。李永乐等选用Laminar层流模型简化模拟具有复杂地形地貌区域(8 km×8 km)平均风场，分析了桥址区风速沿高度方向和沿主梁方向变化特点，并讨论了不同攻角情况下桥面风速标准及其与梯度风速的比值关系^[12]。

综上所述，复杂地形桥位风特性CFD数值模拟越来越受到关注。本文在综述以上研究成果的基础上，以拟建的山西临猗黄河大桥为工程背景，根据大桥桥位处风观测塔的实测数据对桥位附近风场空间分布特性进行数值模拟研究，以获得桥轴线不同位置风速沿高度分布规律，为大桥抗风设计提供参考。

1 工程背景及桥位地形数据 1.1 工程概况

拟建的山西临猗黄河大桥及引线工程是山西省高速公路网规划“三纵十一横十一环”第十一横的重要组成部分。临猗黄河大桥是跨越黄河小北干流禹门口至潼关河段的特大型桥梁，推荐桥位起点位于山西省临猗县孙吉镇西里村西侧，终点位于陕西省合阳县百良镇三汲村附近。山西临猗黄河大桥位于运城附近黄河段，河面较宽，大约有6.4 km，桥位轴线大致垂直于水流方向，两岸整体上“西高东低”，是一个巨大的“U”形河谷，两岸地形陡峻，西岸陕西侧基本上呈黄土台塬地貌，黄土高崖耸立，崖顶至河床高差达200 m左右；东岸山西侧地势较为平缓，基本上属典型的黄土高原地貌，沟壑纵横、大型冲沟两侧有大量小冲沟呈典型的“鸡爪”地形，桥位附近地形如图1所示。根据《公路桥梁抗风设计规范》(JTG/T D60—01—2004)可知，山西运城市百年一遇10 m高度10 min平均时距风速为27.6 m/s。大桥目前处于方案设计阶段，本桥具有工程规模大、墩高、桥位风速大等特点，大桥抗风设计应予以重视。

1.2 桥位地形数据

考虑到临猗黄河大桥东西方向横跨黄河，大桥总长约6 km左右，结合计算机内存大小等因素，选定桥位地形计算域为东西向边长10.2 km，南北向边长为9.0 km，高度方向为2.0 km。计算域下边界地(河)面采用航天飞机雷达地形测量系统(Shuttle Radar Topography Mission,SRTM)所提供的90 m分辨率数字地形高程模型^[13-15]。桥址附近风观测塔处经纬度为E: 110°24′34″，N: 35°14′43″。通过位置确定在中科院网上下载对应地形数据，对应压缩包文件名为srtm_59_05。通过Global Mapper导出桥址处附近地形文件，桥位附近地形如图2所示；再由逆向工程软件Imageware将数字高程地形点云图生成自由曲面地形图，如图3所示。

2 数值模拟方法 2.1 网格划分

临猗黄河大桥桥位附近地形风特性CFD数值模拟计算域大小为10.2 km×9 km×2 km。进行网格划分时，近地面网格采用80 m×80 m的结构化网格，即整个计算域在东西向和南北向分别划分125份和112份，近地面第一层网格厚度分别取0.5，1.0，2.0 m和4.0 m，沿高度方向网格厚度渐变，一共划分40层，以进行不同网格划分计算结果的比较。计算域均采用六面体网格进行划分，以确保计算精度和效率，整个计算域网格总数约为56万。计算域近地面网格及入口边界网格划分如图4所示，整体计算域网格划分如图5所示，具体计算网格划分参数如表1所示。

2.2 边界条件

计算域边界条件设置如下：入口(风从正北方吹)边界条件设置为速度入口边界，根据桥位处风实测结果，取桥位处10 m高度处年平均风速为2.7 m/s，桥位处地表为B类地表，对应风剖面指数取为α=0.16，梯度风高度取为350 m。采用UDF方式设置入口风剖面，如图6所示。出口边界条件设置为出流边界条件(Outflow)。粗糙度设置：粗糙度设置与地表情况密切相关，由于缺乏该地区的地表统计资料(或数字地面模型)，粗糙度取为0.03 m。近壁条件的处理：近壁条件采用FLUENT中的标准壁面函数处理；计算域顶面以及两侧面均采用对称边界条件(Symmetry)。计算中分别在桥位附近风观测塔、主桥总长的中央、四分点分别由低至高设置了50个监测点，共设置200个监测点进行风速监测，监测点位置如图7所示。

2.3 湍流模型

目前湍流数值计算主要有3种方法，即直接数值模拟(DNS)方法、大涡模拟方法(LES)和雷诺平均N-S方程(RANS)方法。综合考虑选用RANS进行湍流模拟，RANS方法首先将满足动力学方程的湍流瞬时运动分解为平均运动和脉动运动两部分，然后把脉动运动的贡献通过雷诺应力项来模拟，即通过湍流模型来封闭雷诺平均N-S方程使之可以求解。综合比较确定采用RNG k－ε模型对临猗黄河大桥桥位地形风特性进行CFD计算。

RNG k－ε模型源于严格统计技术，与标准k－ε模型相似但有如下改进：在ε方程中增加了一个条件，可以有效提高计算精度；考虑了湍流旋涡，提高了在这方面的精度；RNG理论为湍流Prandtl数提供了一个解析表达式；在RNG k－ε模型中，通过在大尺度运动和修正以后的黏度项来体现小尺度的影响，从而使小尺度运动系统地从控制方程中去除。RNG k－ε模型能模拟射流撞击、分离流、二次流和旋流等中等复杂流动。

(1)

(2)

式中，G_k为由平均速度梯度所引起的湍动能，其严格定义为${G_k} = - \rho \overline {u{'_i}u{'_j}} \partial {u_j}/\partial {x_i}$，采用Boussinesq 假设，即G_k=μ_tS²，其中S为平均速度应变率张量，即$S \equiv \sqrt {2{S_{ij}}{S_{ij}}} $；G_b为由浮力所产生的湍动能； Y_M为可压缩湍流脉动扩胀对总耗散率的影响。这些参数与标准k－ε模型相同。α_k与α_ε分别为湍动能k和耗散率ε的有效普朗特数的倒数；S_k和S_ε分别为用户定义的源项。

采用FLUENT中的RNG k－ε模型进行桥位地形风特性分布流场计算。具体求解设置为：采用SIMPLEC(Semi implicit method for pressure linked equation consistent)算法求解动量方程(即N-S方程)中速度分量和压力的耦合问题；空间离散采用二阶中心差分，时间离散采用二阶迎风格式差分，计算时间步长为1 s，共计算10 000步。限于篇幅，本文仅给出风从正北方吹时的计算结果。

3 计算结果 3.1 不同网格计算结果

为了比较不同网格计算精度，分别针对表1中给出的4种网格进行了流场计算。表2及图8分别给出了不同网格对应的1^#监测点位置10，30，50 m和80 m高度处风速结果和时程曲线。从图8中可以看出：当计算到约1 500 s时，不同模型对应的各计算点风速大小趋于稳定，由于地形的影响，风速会有一定的波动；测风塔位置处10，30 m及50 m高度处网格4的计算结果与其余3种网格计算结果相比偏小，80 m高度处4种网格的计算结果总体相差不大。由表2可知，4种网格计算结果总体上与实测值较为接近，最大相对误差约为-14.5%。产生这一误差的主要原因是地形的影响，由于计算时采用了桥位风观测塔10 m高度处的实测平均风速作为入口风速边界，实际上入口边界与风观测塔之间存在较大的距离；另一方面入口风速较小，由地表而产生的脉动风对计算结果影响较大。在实际计算中可采用桥位基本风速作为速度入口边界，以减小地形对风速的影响。综合考虑计算量与精度，选择网格2进行实际桥位地形风特性计算。

图9所示为不同网格计算得到的各监测点位置140 m高度处风速时程曲线。从图9中可以看出，1^#、2^#及3^#监测点位置140 m高度处不同网格计算得到的风速时程在计算步达到1 500步后总体较为平稳，且风速大小基本相同；4^#监测点处各网格计算得到的风速时程存在较大的脉动，主要原因是4^#监测点附近的局部地形较为复杂，当风从北侧吹时，4^#监测点位于前方山丘尾流中，即该区域可能存在较大的瞬时风速。这对大桥抗风设计或施工期的抗风安全措施设置具有一定的参考意义。

3.2 桥位风空间分布特性

由临猗黄河大桥的设计方案可知，大桥跨中处桥面距离水面高度约为140 m，故分别给出近地(河)面以及桥中央处距离河面140 m的高度处风速等值线分布图，如图10所示。从图10中可以看出，桥跨中央距离河面140 m高度处(桥面高度附近)河道上方风速分布总体较为均匀，靠近陕西侧风速变化较为剧烈，而靠近山西侧风速变化相对平缓，这主要是由于黄河东岸山西侧地势较为平缓而陕西侧则地势较为起伏引起的。图11分别给出了桥跨中央、四分点以及风观测塔位置处沿高度方向的风速分布图。从图11中可以看出，风从正北方吹时，各监测点处平均风速沿高度增加而增大。

4 结论

以拟建的山西临猗黄河大桥为工程背景，采用RNG k－ε模型对复杂山区地形桥位风特性进行了CFD数值模拟研究，得到如下主要研究结论：

(1) 采用SRTM所提供的90 m分辨率数字地形高程模型，借助Global Mapper、Imageware等软件可实现复杂地形模型CFD建模，可满足复杂桥位地形风特性数值模拟要求。

(2) 测风塔位置10,30 m及50 m高度处网格4计算结果与其余3种网格计算结果相比偏小，80 m高度处4种网格计算结果总体相差不大；4种网格计算结果总体上与实测值较为接近，最大相对误差约为-14.5%。

(3) 风从正北方吹时，各监测点处平均风速沿高度的增加而增大；1^#,2^#及3^#监测点位置140 m高度处且风速大小基本相同；4^#监测点位于前方山丘尾流，风速脉动较大。