CN105608326B - 一种山区复杂地形风场大涡模拟入口边界条件输入方法 - Google Patents

Abstract

本发明公开了一种山区复杂地形风场大涡模拟入口边界条件输入方法，确定目标区域的大涡模拟计算域；生成目标区域的入口速度信息；将计算域入口边界进行分块处理，然后把入口速度信息按多项式插值的方法分别赋给每个分块；用风速监测仪对所述目标区域入口实际位置的脉动风场特性进行监测，并用最小二乘法拟合成风谱；用谐波合成法对拟合的风谱进行脉动风速合成，生成满足实际入口风场特性的脉动风速；将考虑平均和脉动风速后的风速时程数据赋给大涡模拟的入口边界。本发明可以实现平均风和脉动风的输入，较好的解决了山区风场入口输入时的“人为峭壁”问题，给出了与实际情况更为吻合的入口边界。

Description

一种山区复杂地形风场大涡模拟入口边界条件输入方法

技术领域

本发明涉及一种山区复杂地形风场大涡模拟入口边界条件输入方法。

背景技术

山区复杂地形风场的数值研究对风能、大跨度桥梁的选址研究有着广泛的理论与现实意义，而目前对山区复杂地形的风环境数值研究，合理的入口边界条件给定是其首要难题，同时入口边界的合理性是其计算结果正确的重要影响因素。一般情况下入口边界由平均风和脉动风两部分组成，对于入口平均风速的给定，以往研究都对其做了很大程度的近似和简化，Maurizi等采用1/10的斜坡作为气流过渡段，对14km×15km区域范围的山区地形进行了模拟，但其模拟过程中斜坡形式较为简单，其适用性有待验证；胡朋等用曲线过渡段对山区峡谷边界进行处理，但这种处理会引起人为的来流风攻角；李永乐等采用指数率风剖面近似等效入口风速剖面，对深切峡谷桥址区风场空间分布特性进行了数值模拟研究，考虑不同来流风向角和风偏角后得到了峡谷桥址处的风场分布特性，但是其入口边界形式的给定对经验依赖较强，通过对已有文献的调研与分析，发现目前研究对山区风场入口的给定均基于主观判断，其给定结果与实际风场相比偏差较大，同时，由于山区地形起伏多样，模型在截断过程中会出现“人为峭壁”现象(如图2所示)，导致数值模拟入口边界速度场的给定出现地面速度不为零现象，因此，一系列问题的出现使得山区峡谷风场的研究进展缓慢。近年来，随着气象预报模式的快速发展，中尺度气象预报模式与CFD软件的耦合方法近年来得到了广泛应用，这也给实现山区复杂地形入口边界的合理给定带来了曙光，中尺度模式的分辨率可以达到百米量级，但是中尺度模式对于山区陡峭峡谷地区会失效，有可能出现积分溢出的情况，单纯的中尺度模式不能对山区复杂地形处的风场进行详细分析，只有将中尺度模式与CFD软件结合起来才能得到山区复杂地区的详细流场信息。

采用多项式插值的方法可以获得从大尺度到小尺度上的平均风速，虽有误差，但基本能控制在允许范围内。但对于脉动风速的给定，用普通的插值方法误差很大，目前还不能直接通过简单的降尺方法直接获取小尺度的湍流信息，而脉动特性对山区峡谷桥梁的风致振动研究非常重要。Uchida和Ohya采用大涡模拟的方法对9.5km×5km区域范围内空气流动进行数值模拟，使用粗糙方块制造脉动风，获得计算区域内平均风速与脉动风速情况，但其粗糙元的摆放产生的脉动风场可调性差，只能模拟一种特定的脉动特性，同时通过这种方法模拟的脉动风场与实际的边界条件仍存在一些出入。

发明内容

本发明所要解决的技术问题是，针对现有技术不足，提供一种山区复杂地形风场大涡模拟入口边界条件输入方法。

为解决上述技术问题，本发明所采用的技术方案是：一种山区复杂地形风场大涡模拟入口边界条件确定方法，包括以下步骤：

1)确定目标区域的大涡模拟计算域；

2)生成目标区域的入口速度信息；

3)将步骤1)建立的计算域入口边界进行分块处理，然后把步骤2)得到的入口速度信息按多项式插值的方法分别赋给每个分块；

4)用风速监测仪对所述目标区域入口实际位置的脉动风场特性进行监测，并用最小二乘法拟合成风谱；

5)在步骤3)分块后的基础上，用谐波合成法对拟合的风谱进行脉动风速合成，生成满足实际入口风场特性的脉动风速；

6)将考虑脉动风速后的风速时程数据赋给大涡模拟的入口边界。

所述步骤1)中，计算域的确定方法为：对于山区峡谷桥址风场，计算域选取原则以峡谷桥址所在位置为中心，向外延伸N公里，其中N不超过10km。

所述步骤2)中，入口速度信息生成过程包括：

1)输入需要模拟的地形分析资料；

2)以所述计算域为中心，利用多重嵌套手法对WRF模拟区域的计算域进行划分；

3)输入需要模拟时间段的气象分析资料；

4)运行WRF，通过多尺度耦合分析，将计算域分为若干个区域，从最外层区域到最里层区域依次进行计算，将当前层的计算结果作为该当前层的里面一层的计算边界，最终得到最里层计算域的入口速度信息。

所述步骤6)的具体实现过程包括以下步骤：

1)在ICEM中根据网格尺寸和延伸率换算大涡模拟入口边界网格坐标的理论值；

2)将步骤5)生成的脉动风速赋值给步骤1)中入口边界网格坐标的理论值；

3)将附上理论值的的风速时程数据与入口边界实际坐标值进行检验对接；

4)当理论值与实际坐标值之间的误差大于0.1时，返回步骤1)重新换算理论值，当误差小于0.1时，风速时程赋值成功。

与现有技术相比，本发明所具有的有益效果为：本发明可以实现平均风和脉动风的输入，较好的解决了山区风场入口输入时的“人为峭壁”问题，给出了与实际情况更为吻合的入口边界。

附图说明

图1为WRF计算域；

图2为人为峭壁示意图；

图3为入口边界分块示意图；

图4为UDF程序流程图

图5和图6为峡谷地区8月份风玫瑰图；

图7为风剖面实测α值次数分布图；

图8(a)和图8(b)分别为风速与风向角时程；

图9为现场实测拟合功率谱；

图10为计算域示意图；

图11为计算网格图；

图12为数值模拟过程中桥位处风速监测点布置图。

具体实施方式

本发明中，利用WRF模式结合现场监测数据通过自编程序，实现了大涡模拟入口边界条件的输入，具体实施步骤如下：

1)确定大涡模拟计算域，对于山区峡谷桥址风场，计算域选取原则以峡谷桥址所在位置为中心，向周围数公里进行延伸，对于风能研究，计算域选取以风速较大区域为原则，计算域大小一般在10公里左右。

2)通过WRF模式利用多尺度耦合技术，生成目标区域的入口速度信息，具体步骤为：

1)输入需要模拟的地形分析资料(NCEP提供的全球30s地形数据及MODIS下垫面分类资料)。

2)以步骤1的计算域为中心，利用多重嵌套手法对WRF模拟区域的计算域进行划分，如图1所示。

3)输入需要模拟时间段的气象分析资料(美国国家环境预报中心(NCEP1°×1°)再分析资料)。

4)运行WRF，利用多尺度耦合技术，将计算域分为若干个区域(划分原则：需保证目标区域在WRF模拟区域的正中心，划分过程中网格尺寸从最外层依次向最内层减小)，如图1左图所示，从最外层到最里层依次进行计算，将外层的计算结果依次作为下一层的计算边界，最终得到最里层计算域的入口风速。

5)将步骤1)建立的CFD计算域入口边界进行分块处理(如图3)，然后把上述入口风速信息按多项式插值的方法分别赋给大涡模拟入口的每块区域。

6)用风速监测仪对目标区域入口实际位置的脉动风场特性进行监测，并用最小二乘法拟合成风谱。

7)用谐波合成法对拟合的风谱进行脉动风速合成，生成满足实际入口风场特性的脉动风速。

8)自编UDF程序，将考虑平均和脉动风速后的风速时程数据赋给大涡模拟的入口边界，其程序流程图见附图4。

本发明以张花高速澧水大桥所在峡谷为研究背景。桥梁主跨为856m，属典型的山区峡谷大跨径桥梁，桥位所处峡谷谷顶宽420m，谷顶与谷底高差280m。

澧水大桥风速监测系统由观测站、桥塔站和桥跨站组成。观测站根据盛行风方向设立在桥址西南侧，位于本发明数值模拟入口附近，监测站风速仪布置在平坦地区，可以认为其风场特性与数值模拟入口接近，其特性可为数值模拟的入口边界条件取值提供参考，风速监测仪布置高度离地面10m，如图5所示。

桥跨站由三个Young 81000三维超声风速仪组成，布置位置为盛行风方向同侧。桥塔站也沿高度布置了三个Young 81000三维超声风速仪，具体位置如图1所示(其中黄色星标识为风速仪安装位置)，风速仪采样频率为4Hz。为了实现风速时程的实时观测，课题组利用GPRS无线传输系统，将现场风速实时数据远程传输至学校的风速采集中心。

为得到峡谷风场特性，通过对2014年8月的风速时程进行分析，得到了观测站和桥跨站的风玫瑰图，如图5和图6所示。

从图中可以发现8月主导风向为西南风，风向与观测站风速仪布置方位一致。同时，对桥塔不同高度风速时程进行监测，得到了强风作用下风剖面α值出现次数分布情况，如图7所示，从图7中可以发现α值呈正态分布，均值为0.3011，综合考虑澧水大桥桥位地形条件，本发明α值取0.3。

为使数值模拟来流风向角与现场实测保持一致，本发明风速时程取自2014年8月29日观测站实测结果。通过对时程数据的分析可得到监测站平均风速为4.43m/s，平均风向角约为180°，风速和风向角时程如图8(a)和图8(b)所示，通过对风速时程分析可得湍流强度为20.6％。

对监测站风速进行风谱分析，用Kaimal谱形式进行拟合，结果如图9所示，其功率谱拟合公式可表示为：

其中，S_u(n)为顺风向功率谱密度函数(即风谱)，n为风的脉动频率，f＝nZ/U(Z),u^*为气流摩阻系数。

几何模型建立过程中，首先通过空间地理数据云获得大范围的地形高程数据，然后用Global Mapper做进一步处理，从而得到目标区域的地形模型。山体模型采用实际尺寸，计算区域大小取10km×9km×4km，如图10所示。

为保证计算精度，数值模型采用全六面体网格，网格在近地面进行加密，最底层网格高度为1m，高度方向在近地面处网格延伸率为1.05，远离地面网格延伸率为1.15，总网格数为6752495，计算网格通过了无关性测试，如图11所示。

本发明数值模型入口处最低高程为164米，桥跨站风速仪安装高度为505m，桥塔站风速仪安装高度为620m，为使峡谷风场入口速度尽量接近现场实测值，本发明入口边界条件包含了平均风和脉动风两部分，平均风速采用分段函数形式给定，其中，近地面处采用指数率形式，为使风速在峡谷内的变化趋势接近实际情况，本发明指数率变化段高度取836m。平均风剖面的具体表达形式为：当Z≦164m时，速度为0m/s，当164m<Z<1000m时，速度遵循指数率规律，其中α取0.3，当Z≧1000m时，速度为7.5m/s，具体可表示为：

V＝0m/s，Z≦164m

V＝7.5m/s，Z≧1000m

本发明脉动风速是基于谐波合成法进行合成的，合成过程中功率谱根据观测站风速时程等效而来，因此，本发明所用入口边界能够满足观测站的风场特性。通过编制UDF程序对商业软件Fluent进行二次开发，将合成的随机风速时程赋给数值模型入口所对应的网格坐标，其中顺风向功率谱采用公式(1)所述形式，竖向功率谱采用Lumley and Panofsky谱，可表示为：

其中，S_w(n)为脉动风竖向功率谱密度函数，Z为入口实际位置距离地面的高度，K为无量纲常数，本发明取K＝0.4，z₀为地表粗糙高度，本发明取z₀＝1m；z_d为零平面位移(zero plane displacement)。

本发明在数值模拟过程中，除入口边界条件采用用户自定义外，地表采用无滑移边界条件，顶面采用自由滑移边界条件，侧面采用对称边界，出口采用压力出口边界。求解方面，本发明的N-S方程采用PISO方法进行求解，对流项和扩散项均采用二阶中心差分格式，用超松弛方法(SOR)求解压力Poisson方程，压力和动量松弛因子分别取0.3和0.7，在满足柯朗数(CFL)的前提下，时间步长取0.1秒。

数值模拟过程中，在主梁水平方向布置了9个风速监测点，竖向方向在1/2跨，3/5跨，7/10跨处分别布置了20个监测点，监测点具体位置如图12所示。

Claims (2)

1.一种山区复杂地形风场大涡模拟入口边界条件确定方法，其特征在于，包括以下步骤：

1)确定目标区域的大涡模拟计算域；

2)生成目标区域的入口速度信息；入口速度信息生成过程包括：

a)输入需要模拟的地形分析资料；

b)以所述计算域为中心，利用多重嵌套手法对WRF模拟区域的计算域进行划分；

c)输入需要模拟时间段的气象分析资料；

d)运行WRF，通过多尺度耦合分析，将计算域分为若干个区域，从最外层区域到最里层区域依次进行计算，将当前层的计算结果作为该当前层的里面一层的计算边界，最终得到最里层计算域的入口速度信息；

3)将步骤1)确定的计算域的入口边界进行分块处理，然后把步骤2)得到的入口速度信息按多项式插值的方法分别赋给每个分块；

4)用风速监测仪对所述目标区域入口实际位置的脉动风场特性进行监测，并用最小二乘法拟合成风谱：

$\frac{{\mathrm{nS}}_u\left(n\right)}{{\left(u^{*}\right)}^2}=\frac{124f}{{(1+28.5f)}^{5/3}};$

其中，S_u(n)为顺风向功率谱密度函数，即风谱；n为风的脉动频率，f＝nZ/U(Z),u^*为气流摩阻系数；Z和U(Z)分别表示入口实际位置到地面的高度和当地的平均风速；

5)在步骤3)分块后的基础上，用谐波合成法对拟合的风谱进行脉动风速合成，生成满足实际入口风场特性的脉动风速；

6)将考虑脉动风速后的风速时程数据赋给大涡模拟的入口边界。

2.根据权利要求1所述的山区复杂地形风场大涡模拟入口边界条件确定方法，其特征在于，所述步骤6)的具体实现过程包括以下步骤：

ⅰ)在ICEM中根据网格尺寸和延伸率换算大涡模拟入口边界网格坐标的理论值；

ⅱ)将步骤5)生成的脉动风速赋值给步骤ⅰ)中入口边界网格坐标的理论值；

ⅲ)将附上理论值的的风速时程数据与入口边界实际坐标值进行检验对接；

ⅳ)当理论值与实际坐标值之间的误差大于0.1时，返回步骤ⅰ)重新换算理论值，当误差小于0.1时，风速时程赋值成功。