山区桥址处CFD计算域的选取方法

张亮亮，吴波，杨阳，刘连杰

(重庆大学 土木工程学院，重庆 400045)

山区桥址处CFD计算域的选取方法

张亮亮，吴波，杨阳，刘连杰

(重庆大学 土木工程学院，重庆 400045)

山区桥址处风场具有较强的随机性与不确定性，在选取其CFD计算域时，盲目参考已有的工程可能会造成较大的模型误差，或增加大量的计算开销。为解决该问题，提出了山区复杂地形CFD计算域选取的方法与步骤，验证了其准确性。具体方法为：设置一个大范围的基准计算域进行初算，通过后处理软件绘制平均风压系数极差随高度的变化曲线、壁面附近水平面上的静压偏差等值云图以及横风向各平面静压值与对应边界面的均方根差值曲线，分别筛选出基准计算域在高度方向、顺风向及横风向上对整体风场贡献可忽略的区域，余下部分则可用于该区风场的求解。

计算域设置；山区地形；计算流体力学；模型误差；求解效率

在山区建设大跨度桥梁时，风环境参数对桥梁结构的设计与安全评价具有重要的参考价值，是必须首要解决的问题之一[1-2]。现有规范对山区风特性的描述较少，风洞实验难以模拟大范围区域的风场，现场实测又难以捕捉大范围的风特性数据，而计算风工程(CFD)具有周期短、数据全面、费用低等优点[3-4]。随着计算机技术的发展，CFD在实际工程中的应用越来越广泛，其模拟精度也得到了进一步的提高[5]。

利用CFD技术进行风场模拟时，首先要进行计算域的设置。CFD数值模拟的误差主要包括模型误差、离散误差与迭代误差[6]。其中，计算域设置不合理所引起的模型误差是影响计算精度的首要原因，应尽量消除。若计算域设置过小，则模型误差较大，即使将离散误差、迭代误差控制到最低，也无法准确反映实际的风场特性[7]。Fujiwara等[8]发现，不同大小的计算域对同一位置的计算结果差异很大，认为模型边界应远离计算点足够远。但是，若将计算域设置太大，虽可降低模型误差，但却大大增加了计算开销，降低了求解效率。尤其在进行高精度、多工况的实际模拟时，过大的计算域对于计算周期的耗费是难以估量的[8]。如何设置大小合适的计算域，兼顾计算精度与求解效率，是笔者主要解决的问题。目前已有一些关于CFD计算域设置的研究或建议[9-12]。崔利民等[9]针对一个双向对称、孤立的正弦山丘提出了山体地形下低矮房屋数值风洞模拟的计算域设定方法；Franke等[12]给出了群体建筑物CFD风场的竖向、侧向与顺风向长度建议值。然而，这些研究均是针对孤立山丘或规则建筑物的风场，其风场分布对称、较为规则且有迹可循。而山地地形在大范围区域内包含不计其数的山丘、沟壑、河流、房屋等，风场具有较强的随机性与不确定性，不能简单地参照已有工程的经验进行设置，针对规则风场所提出的方法或建议值也不适用。有经验的学者在进行山区CFD模拟时，往往先进行试算[13-14]。笔者以地处长江交汇口山地地貌的大宁河特大桥为依托，提出山区桥址处复杂地形CFD计算域设置的一般方法与步骤。

1 工程概况

重庆大宁河特大桥是G42重庆段二期高速公路的重要工程，桥型为上承式无铰钢桁拱桥，跨度400 m，矢高80 m，桥轴线垂直穿过大宁河和两岸山坡。桥址区位于大宁河与长江交汇口附近，地表为缓、陡相间的折线型斜坡，坡度55°～75°，属于典型的山地地貌。桥位范围内最大地面标高为518 m左右，河底标高仅为90 m，其相对高差达428 m，切割深度大。

现场为期2 a的实测数据表明，桥址处的主导风向为东北风和东南风，100 a重现期的基本风速为26.8 m/s。

2 计算域的设置方法

2.1 基准计算域

2.1.1 计算域大小 定义如下参数：H为壁面最高点离计算域顶面的高度，LE、LW、LS、LN分别为桥址中心距计算域东、西、南、北边界面的距离(图1)。

山区桥址处的CFD计算域一般为长/宽5～15 km、高1～3 km[13-15]的长方体区域。先设置一个较大范围的基准计算域进行初算，将桥位附近有代表性的地貌全部包含在内，其大小为25 km(桥轴方向)×20 km(桥轴法向)×6 km(竖直方向)，即LE=LW=12.5 km，LS=LN=10 km，H=6.0 km。基准计算域的作用为：从其计算结果中提取相关指标，作为设置H、LE、LW、LS、LN的依据；将其解作为标准解，计算各检验计算域的模型误差。

图1 基准计算域及定义的参数Fig.1 Reference computational domain and its parameters

使用Google Earth获取地形底面的高程数据，取样间隔30 m，共计获得584 714个离散高程点。将其导入逆向工程软件Imageware，拟合四阶地形曲面。然后将曲面导入网格划分平台Gambit，形成计算域。设置边界条件：入流面取速度入口Velocity-inlet，出流面为自由流Outflow，地形底面为Wall壁面，其余为对称边界Symmetry[13]。

2.1.2 网格划分与求解 为便于基准计算域与后续工作中检验计算域的对比，所有计算域均分块为内部区域与外部区域。内部区域是风场计算所重点关心的区域，覆盖桥址中心附近5 km×4 km的范围，壁面划分60 m尺度的三角形网格；外部区域的壁面划分90 m尺度的非结构化网格。高度方向上，靠近壁面的第一层网格厚度10 m，逐渐加大，增长因子1.1，最大尺度150 m。基准计算域共划分棱柱体网格单元3 542 798个。

将网格文件导入Ansys Fluent 15.0，选用全隐式分离求解器，时间、空间离散均采用二阶差分格式，压力与速度耦合选用SIMPLE算法。本工作的目的在于选择计算域大小，尚未进入风场特性的高精度模拟的阶段。因此，为提高工作效率，采用层流(laminar)模型，来流设置为均匀流，风速取本区的基本风速26.8 m/s。

2.2 高度H的设置

在模拟孤立建筑物的风场时，一些学者[10-11]将平均风压系数极差ΔCp(z)作为设定计算域高度H的依据，认为ΔCp(z)越大，该高度的水平平面对建筑物周围风场的影响越大；反之，影响越小，即该平面已不在建筑物的风场区。山地地貌的风场虽复杂得多，但在竖直方向上具有和孤立建筑物风场同样的特征[10]。因此，将ΔCp(z)引入到为山地风场，作为计算域高度设置的依据。基准计算域的ΔCp(z)计算式为

(1)

式中：Cp为测点平均风压系数;pi为测点静压；p∞为参考点(计算域顶面)静压；vH为参考点风速；ρ为空气密度；Cp,max(z),Cp,min(z)分别为基准计算域Z高度水平面上平均风压系数的最大值、最小值，二者差值即该平面的平均风压系数极差ΔCp(z)。

图2给出了基准计算域的ΔCp(z)曲线。可见，ΔCp随离地高度Z的增加而逐渐减小，即，离壁面越高的平面受地形的影响越小。对于Z=3.5 km以上的平面，其ΔCp均低于0.01，远小于Z=3.5 km以下的部分，而且其曲线斜率接近于0，几乎不再变化。说明基准计算域Z=3.5 km以上的区域对整体风场的贡献微弱，已不在壁面风场区范围内。

图2 基准计算域的平均风压系数极差ΔCp(z)曲线Fig.2 ΔCp(z)curve of reference computational domain

2.3 顺风向长度LE、LW的设置

崔利民等[9]以Y轴中心线(顺风向)上的相对静压变化曲线作为设置计算域长度的指标。但其计算模型是一个双向对称的、孤立的正弦山丘，其静压分布关于Y轴对称。然而对于山地地貌的风场而言，Y轴中心线附近的风场与边界附近风场相去甚远，崔利民的方法难以全面反应水平面上风场的分布。

CFD模拟风场时，顺风向的两个边界(x方向，本例对应LE、LW)分别为速度入口与压力出流。为评价计算域东、西面部分对整体风场的贡献，经试算与筛选，发现静压偏差δP可综合反映风场在x方向上的分布，其计算式为

(2)

为保持出流面的一致，避免额外误差，在评价东面部分对整体的影响时，以东面边界为入流面；同样，在评价西面部分的影响时，以西面边界为入流面。

由图2可知，离壁面越高的平面对风场的影响越小。偏于保守地考虑，选择壁面最高点附近的水平面(Z=0.5 km)为代表，在后处理软件Tecplot中绘出其δp(z)的等值云图(图3)。

图3 基准计算域Z=0.5 km平面的静压偏差δp等值云图(大于均方根值S(δp)者)Fig.3 δp-contour of reference computational domain at Z=0.5km plane

图3给出的是大于均方根值S(δp)的部分——贡献大的区域。就图(a)中x=8.5～12.5 km、图(b)中x=-8～-12.5 km的区域而言，其δp几乎均低于S(δp)，说明基准计算域以东8.5～12.5 km、以西8～12.5 km的区域对整体的贡献很小。

2.4 横风向宽度LN、LS的设置

CFD模拟风场时，横风向的两个边界面通常设置为对称边界Symmetry(y方向，本例对应LN、LS)。其上所有物理量的梯度(grad(φn))均为0。若要以其他面代替现有边界，则该面上各物理量的解必须与现有边界的解十分接近。

以0.5 km的增量，分别在基准计算域y=2.0～9.5 km的各平面上以200 m的间距设置监控点，每个面共计2 500个；考察这些监控点的静压与北边界面(y=10.0 km)上对应点的差值，差值的均方根计算式为

式中：pi,y为y平面上第i个监控点的静压；pi,0为北边界面上对应点的静压。

同样地，在y=-2.0～-9.5 km的各平面上设置监控点，考察其与南边界面(y=-10.0 km)上对应点的差值，均方根差值的计算方法同式(3)，不过pi,0对应于南边界面。

Sp(y)曲线见图4。由图4(a)可见，y=-8.0～-9.5 km各平面的Sp均低于0.5%，曲线斜率接近于0。说明将南面边界移动至计算域中心以南8.0～9.5 km的区域所造成的误差是可以忽略不计的。同样，由图4(b)，在不显著增加计算误差的情况下，可将北面边界移动至计算域中心以北6.0～9.5 km的区域。

图4 基准计算域各y平面与南/北边界的静压均方根差值Sp(y)Fig.4 Sp(y) curve of reference computational domain

实际上，y=-6.5～-8.0 km的各平面的Sp也低于0.5%，在实际工程中基本上可忽略不计。但偏于保守地认为，则以接近于0的斜率趋近于0的曲线段才是对整体贡献可以忽略的部分。显然，y=-6.5～-8.0 km曲线段的斜率并不接近于0，因此，不认为该部分对整体的贡献可以忽略。

按这3种方法分别进行高度、顺风向长度、横风向宽度的设置，舍弃掉对整体贡献可被忽略的区域，所选定的用于实际求解的计算域为：H=3.5 km、LE=8.5 km、LW=8.0 km、LS=8.0 km、LN=6.0 km，其体积仅为基准计算域的0.27倍。

2.5 计算域设置步骤

山区桥址处CFD计算域的设置步骤为：

1)设置一个大范围的基准计算域，将本区域内所有具代表性的地貌包含在内；获取本区气象数据，以本区的盛行风向为入口，以相对于实际求解较疏的网格、较低阶的计算方法进行快速初算；

2)绘制基准计算域的平均风压系数极差随离地高度Z的变化曲线ΔCp(z)，以曲线上趋近于零值的拐点高度作为实际计算域的高度；

3)分别以顺风向(盛行风)上的两个边界面为速度入口进行求解，绘制其壁面附近水平面的静压偏差云图δp，舍弃掉低于其均方根值S(δp)的部分，余下部分的长度即为实际计算域的顺风向长度；

4)在横风向上的各平面设置监控点，考察各监控点静压与对应两个边界面的差值，分别绘制其均方根误差曲线Sp(y)，以趋近于零值的拐点作为实际计算域的横风向长度。

3 算法验证

计算域高度、顺风向长度、横风向宽度的设置方法是在不断试算、反复筛选的基础上而提出的，必须对其进一步的验证。验证思路为：按上述方法已经评价出对整体风场贡献可忽略不计的区域，分别舍弃掉这些区域，建立若干个检验计算域。对比检验计算域的解与标准解(基准计算域的解)的差异。若差异可忽略不计，说明舍弃掉区域不会增强显著的模型误差，即验证了上述方法的准确性；若差异不可忽略，则说明上述方法并不适用。

3.1 检验计算域

以0.5 km为增量，分别将H、LE、LW、LS、LN作为唯一变量，设置H=1.0～5.5 km的10个、LE=2.5～12.0 km的20个、LW=2.5～12.0 km的20个、LN=2.0～9.5 km的18个、LS=2.0～9.5 km的18个(共计86个)检验计算域。

实测数据表明，该区的主导风向为东北风与东南风。因此，除LE=2.5～12.0 km的20个检验计算域外，其余检验计算域均以西面边界入流，形成由西向东的流向。与之对比的基准计算域也取西面边界为入流面进行求解，其网格划分、求解方法均与各检验计算域完全一致，从而保证了H、LW、LS、LN分别作为变量的唯一性；对于LE=2.5～12.0 km的20个检验计算域，若以西面边界入流，则各计算域与基准计算域的出流面不一致，可能导致额外的误差。因此，这20个检验计算域均以东面边界入流，形成由东向西的流向，与之对比的基准计算域也取东面边界入流。

3.2 五参数对计算误差的影响

在各检验计算域内部区域Z=1/4Hj、2/4Hj、3/4Hj、Hj(Hj为检验计算域的高度)处的水平面上以100 m的间隔设置2 000个监控点，每个检验计算域共计4×2 000=8 000个监控点。对比检验计算域各监控点的顺风向、横风向、竖向风速与基准计算域对应点的差值，其均方根误差的计算式为

(4)

式中：vi,j为检验计算域第i个监控点的风速；vi,0为基准计算域对应点的风速。若某检验计算域的顺风向速度vx、横风向速度vy、竖向速度vz的误差(Svx、Svy、Svz)均趋近于0，且所在曲线段斜率接近于0，则认为其模型误差可忽略不计。

Svx、Svy、Svz随H、LE、LW、LS、LN的变化曲线如图5。由图5可见，模型误差可忽略不计的检验计算域有：H=3.5～5.5 km的5个，LE=8.5～12.0 km的8个，LW=8～12.0 km的9个，LN=6.0～9.5 km的8个以及LS=7.5～9.5 km的5个。说明H=3.5～6.0 km、LE=8.5～12.5 km、LW=8～12.5 km、LN=6.0～10.0 km、LS=7.5～10.0 km这些区域是可以被舍弃的。

而按照该设置方法，选定的计算域为：H=3.5 km，LE=8.5 km，LW=8.0 km，LS=8.0 km，LN=6.0 km。可见，除LW偏于保守外(多估计了0.5 km)，其余参数均与与Svx、Svy、Svz误差曲线的结果准确吻合，从而验证了该方法的准确性。

图5 高度(图(a))、顺风向(图(b)、(c))及横风向长度(图(d)、(e))对计算误差Svx、Svy、Svz的影响曲线Fig.5 Impact of H、LE、LW、LSand LNon calculation error

4 结 论

1)山区风场具有较强的随机性与不确定性，在选取其计算域时，不可盲目参照已有的工程经验。

2)提出了山区桥址处CFD计算域的选取方法，验证了其准确性。

3)进行了一些偏于保守的处理：认为曲线上以近于0斜率趋近于零值者才是对整体贡献可忽略的部分；工程运用时可视精度需要选择贡献大小的分界；以受壁面影响最大的水平面代表计算域在顺风向的分布，工程运用时可视精度需要选择离壁面稍远一些的平面。

4)计算域的选取工作尚未涉及实际求解的高精度需求，因此,可以较疏的网格、较低阶算法进行快速初算；再通过简单处理即可准确、直观地筛选出实际求解的计算域，并不会增加额外的计算负担。

5)本例所选定的实际计算域体积仅为基准计算域的0.27倍，在进行密网格、高阶算法、复杂湍流模型、多工况的实际求解时，可大大减少工作量，缩短计算周期。

[1] 周志勇，肖亮，丁泉顺，等.大范围区域复杂地形风场数值模拟研究[J].力学季刊，2010，31(1)：101-107. Zhou Z Y,Xiao L,Ding Q S,et al. Numerical simulation study of wind environment for the flow around large region with complex terrain [J]. Chinese Quarterly of Mechanics,2010,31(1): 101-107. (in Chinese)

[2] 贺德馨.我国风工程研究现状与展望[J].力学与实践，2002，24(4)：10-19. He D X. The status and vista of wind engineering studies in China [J]. Mechanics in Engineering,2002,24(4): 10-19. (in Chinese)

[3] 项海帆.进入21世纪的桥梁风工程研究[J].同济大学学报:自然科学版,2002,30(5): 529-532. Xiang H F. Study on bridge wind engineering into the 21st century [J]. Journal of Tongji University: Natural Science,2002,30(5): 529-532. (in Chinese)

[4] Ding Y L,Zhou G D,Li A Q,et al. Statistical characteristics of sustained wind environment for a long-span bridge based on long-term field measurement data [J]. Wind and Structures,2013,17(1): 43-68.

[5] Blocken B,Stathopoulos T,Carmeliet J,et al. Application of computational fluid dynamics in building performance simulation for the outdoor environment: An overview[J]. Journal of Building Performance Simulation,2011,4(2): 157-184.

[6] 曾锴，汪丛军，黄本才，等.计算风工程中几个关键影响因素的分析与建议[J].空气动力学学报，2007，25(4)：504-508. Zeng K,Wang C J,Huang B C,et al.Suggestion and analysis of several key factors in computational wind engineering [J]. Acta Aerodynamica Sinica,2007,25(4): 504-508. (in Chinese)

[7] Zhu Z W,Liu Z Q. CFD prediction of local scour hole around bridge piers [J]. Journal of Central South University of Technology: English Edition,2012,19(1): 273-281.

[8] Fujiwara T. Application of statistical system identification method to the analysis of wind-induced current and wind set-up [J]. Coastal Engineering Journal,1985,28: 117-123.

[9] 崔利民，彭兴黔，时凌琳，等.山体地形下低矮房屋数值风洞模拟的计算域设定[J].华侨大学学报:自然科学版，2010，31(4)：463-467. Cui L M,Peng X Q,Shi L L,et al. Computational domain setting about numerical wind tunnel in the simulation of the low-rise housing in the mountain terrain [J]. Journal of Huaqiao University: Natural Science,2010,31(4): 463-467. (in Chinese)

[10] Jakob N,Gurevych I. Extracting opinion targets in asingle-and cross-domain setting with conditional random fields [C]// Conference on Empirical Methods in Natural Language Processing,Cambridge,MA,United States,2010: 1035-1045.

[11] Hu W L,Zhang Y H,Wang L B. Numerical simulation on turbulent fluid flowand heat transfer enhancement of a tube bank fin heat exchanger with mounted vortex generators on the fins [J]. Journal of Enhanced Heat Transfer,2010,18(5): 361-374.

[12] Franke J,Hellsten A,Heinke K,et al. The COST 732 best practice guideline for CFD simulation of flows in the urban environment: A summary [J]. International Journal of Environment and Pollution,2011,44: 419-427.

[13] 祝志文，张士宁，刘震卿，等.桥址峡谷地貌风场特性的CFD模拟[J].湖南大学学报:自然科学版，2011，38(10)：13-17. Zhu Z W,Zhang S N,Liu Z Q,et al. CFD simulation of wind field at bridge site on gorge terrain [J]. Journal of Hunan University: Natural Science,2011,38(10): 13-17. (in Chinese)

[14] 李永乐，蔡宪棠，唐康，等.深切峡谷桥址区风场空间分布特性的数值模拟研究[J].土木工程学报，2011，44(2)：116-122. Li Y L,Cai X T,Tang K,et al. Study of spatial distribution feature of wind fields over bridge site with a deep-cutting gorge using numerical simulation [J]. China Civil Engineering Journal,2011,44(2): 116-122. (in Chinese)

[15] Huang W F,Xu Y L. Prediction of typhoon design wind speed and profile over complex terrain [J]. Structural Engineering and Mechanics,2013,45(1): 1-18.

(编辑 胡英奎)

Method of setting up the wind field of a mountain bridge site

ZhangLiangliang,WuBo,YangYang,LiuLianjie

(School of Civil Engineering,Chongqing University,Chongqing400045,P. R. China)

The wind field of a mountain bridge site always shows a strong randomness and uncertainty. As a result,setting up its computational domain in CFD simulation by a simply reference to the existing experience causes a large amount of model errors and huge additional computational overhead. To solve this problem,a method to the selection of CFD computational domain of complex terrains is explored and verified. To start the procedure,a large reference computational domain should be chosen and preliminary solved. Then,its solution should be post-processed to draw the following curves or contour images: range of mean pressure coefficient align with height,static pressure deviation contour of the horizontal plane near the top of the bottom surface,the root-mean-square error of the static pressure of the crosswind planes compared to the corresponding edge surface. According to these curves and contour images,areas parts which have slight contribution on the overall wind field are figured out and abandoned,and the selected computational domain is constitutive of the rest parts of the reference domain.

computational domain settings; mountain terrain; CFD; model error;computational efficiency

10.11835/j.issn.1674-4764.2015.05.002

2015-07-05 基金项目：国家自然科学基金(51578098)

张亮亮(1956-)，教授，博士生导师，主要从事桥梁力学性能分析、桥梁抗风研究，(E-mail)zll200510@126.com。

Foundation item：National Natural Science Foundation of China(No.51578098)

V211.3;O368

A

1674-4764(2015)05-0011-07

Received:2015-07-05

Author brief：Zhang Liangliang (1956-),professor,doctoral supervisor,main research interests:bridge mechanics performance analysis & wind resistance of bridge,(E-mail) zll200510@126.com.