CN108534982A - 一种用于山区风场模拟的过渡曲线 - Google Patents

Abstract

本发明公开一种用于山区风场模拟的过渡曲线，其推导过程如下：首先根据风洞收缩段曲线：双三次曲线、双三次曲线和维氏曲线推导出地形过渡段曲线，再通过二维CFD模拟对得到的过渡段曲线进行评估，对比不同过渡曲线的平均加速比、沿不同过渡曲线的平均风攻角和沿不同过渡曲线的平均额外湍流强度，并引入评估综合指标，提出新的过渡曲线。本发明能够很好的模拟山区深切峡谷风场，具有很好的过渡性。

Description

一种用于山区风场模拟的过渡曲线

技术领域

本发明涉及山区风场模拟技术领域，具体为一种用于山区风场模拟的过渡曲线。

背景技术

在利用CFD(Computational Fluid Dynamics计算流体动力学)或者风洞实验进行山区风场模拟时，地形模型的有限范围会导致山区地形在合适的位置需要被截断，这就会造成“人造悬崖”，如图1所示。特别是深切峡谷这类复杂的地形，落差高达几百上千米，因而这个在风洞地板(或者数值计算域底部)和地形边缘之间的“人造悬崖”也会形成巨大的高差，并且引起气流在模型边缘分离，从而可能会造成气流加速，改变风攻角甚至产生不真实的流型。为了减小“人造悬崖”对流场造成的变异，一个连接风洞地板(或者数值计算域底部)和地形边缘的过渡段是必需的。Hu et al.(2015)基于无粘流体圆柱绕流理论提出了一种理想曲线，并且借助风洞实验比较了二维理想曲线和二维传统斜坡分别对气流的过渡效果。结果显示理想曲线有更好的过渡表现。并且，Li et al.(2017)将这类二维理想曲线应用到一个真实的深切峡谷风洞实验当中。尽管这个理想曲线比斜坡在过渡气流方面有更好的表现，但是它的推导过程并未考虑空气粘度。因此将它应用于模拟真实地形的三维风场还需要深入研究。

发明内容

针对上述问题，本发明的目的在于提供一种过渡效果更好，可以用于模拟山区深切峡谷风场的过渡曲线的推导方法。技术方案如下：

一种用于山区风场模拟的过渡曲线的推导方法，包括以下步骤：

步骤1：根据风洞收缩段曲线推导出地形过渡段曲线；

a)引入风洞收缩段曲线：双三次曲线、双三次曲线和维氏曲线：

双三次曲线是由两条三次曲线由拐点连接起来的，其表达式如下：

其中，h是在轴向距离x处的横截面的半高，H_i是入口截面的半高，H_o是出口截面的半高，L是收缩段的长度，x_f是拐点的轴向距离；

五次曲线的表达式如下：

h＝[-10(x/L)³+15(x/L)⁴-6(x/L)⁵][H_i-H_o]+H_i

维氏曲线的表达式如下：

其中，a为调整参数；

b)定义过渡段曲线起点和终端分别为S(0,0)和E(L,H)；其中H是过渡曲线的高度，相当于收缩段中的H_i-H_o，L是过渡段的长度，相当于收缩段的长度L；

相应地，定义等效斜率为K＝H/L，等效倾角为β＝tan^-1(K)；横坐标x对应的纵坐标值由 H_i-h决定；

则将风洞收缩段曲线表达方式用于地形过渡段，双三次过渡曲线由下式确定：

五次过渡曲线表示为：

y＝H[10(x/L)³-15(x/L)⁴+6(x/L)⁵]

对维氏曲线收缩段，假设H_i→∞，得到适用于地形模拟的维氏过渡曲线：

其中，α为评估变量，上式取α＝1/3；

步骤2：通过二维CFD模拟对得到的过渡段曲线进行评估：

a)定义通过过渡曲线后的风速在x处的平均加速比S_x为：

其中，N是地面至临界高度H_w的网格数量，临界高度H_w是风速变化小于5％的高度；v_i是x位置距地面第i个网格的水平风速；v_0i是入口处距地面第i个网格的水平风速；Δh_i是第i个网格和第i-1个网格的高差；

对比x_f取不同值的双三次过渡曲线、五次过渡曲线及α取不同值的维氏过渡曲线中各不同过渡曲线的平均加速比，得出：且当α值较大时，维氏过渡曲线的平均加速比上升最快且也最早开始下降至更低值；

c)定义通过过渡曲线后在x处的平均风攻角A_x为：

其中，a_i是x位置处距地面第i个网格风攻角值；

对比沿所述各不同过渡曲线的平均风攻角，得出：当α值较大时，维氏过渡曲线在风攻角方面表现更好；

d)定义通过过渡曲线后在x位置处超过入流处的平均额外湍流强度T_x为：

其中，T_i为x位置处距离地面第i个网格的湍流强度；T_0i为距离地面第i个网格的入流湍流强度；

对比沿所述各不同过渡曲线的平均额外湍流强度，得出：当α达到一个较大值时，维氏过渡曲线的湍流强度突然增大；

步骤3：通过对过渡曲线的评估提出新的过渡曲线：

引入评估综合指标CI：

其中，S_j,A_j和T_j分别是所述各不同过渡曲线的平均加速比、平均风攻角和平均额外湍流强度值；S_r,A_r和T_r分别是参考过渡曲线的平均加速比，平均风攻角和平均额外湍流强度值；

选取CI值最小的曲线为新的过渡曲线。

进一步的，所述CI值最小的曲线为α取50时的维氏过渡曲线，其平均加速比、平均风攻角和平均额外湍流强度值在x＝0.6位置处均小于5％，则过渡曲线在x＝0.6位置处被截断，在L取单位长度为1的情况下，为了保证截断后曲线的形状，取新的等效斜率K_o＝K/0.6，原等效斜率K＝0.5774，则K_o＝0.9623，从而新的过渡曲线表达式如下：

本发明的有益效果是：本发明考虑了空气粘度，能够用于模拟山区深切峡谷风场，具有更好的过渡性。

附图说明

图1为过渡曲线用于连接截断后的山区地形。

图2为风洞收缩段曲线示意图。

图3为过渡曲线。

图4为几种过渡曲线的对比。

图5为BCTC曲线(图5a)和WTC曲线(图5b)的变形。

图6为计算域和边界条件。

图7a、7b为G1网格的示意图。

图8为沿W50过渡曲线的沿程风剖面。

图9为不同过渡曲线的平均加速比。

图10为沿不同过渡曲线的平均风攻角分布。

图11为沿不同过渡曲线的平均额外湍流强度分布。

图12为桥梁、山谷和测风塔在山区的具体位置。

图13为测风塔示意图。

图14为风向风速玫瑰图。

图15为截断后的实际数值模型。

图16为实际地形过渡段曲面。

图17为计算域网格。

图18为地形模型网格示意图；18a地形表面网格，18b面网格细化示意图。

图19为CFD入流方向。

图20为50m高度处实测和CFD模拟得到的50m处平均竖向风攻角值。

图21为风剖面指数的概率。

具体实施方式

下面结合附图和具体实施例对本发明做进一步详细说明。本实施例中用于山区风场模拟的过渡曲线的推导方法，具体步骤如下：

步骤1：基于风洞收缩段的地形过渡曲线。

风洞收缩段是风洞实验室最重要的组成部分之一，如图2所示。它的主要作用是均匀加速空气至试验段需要的风速。此外，需要保证在出口位置由于逆压产生的湍流尽量小。目前广泛使用的风洞收缩段包括双三次曲线、五次曲线和维氏曲线。

双三次曲线是由两条三次曲线由拐点连接起来的，其表达式如下，

其中h是在轴向距离x处的横截面的半高(或半径)，H_i是入口截面的半高(或半径)，H_o是出口截面的半高(或半径)，L是收缩段的长度，x_f是拐点的轴向距离，详见图3所示。五次曲线表达式如下：

h＝[-10(x/L)³+15(x/L)⁴-6(x/L)⁵][H_i-H_o]+H_i (2)

维氏曲线表达式：

其中a为调整参数。

一般来说，模拟山区深切峡谷地形时，合适的过渡段需要尽量减小由“人造悬崖”引起的不利影响并且使得在地形边缘处的来流风跟目标风场尽可能一致。根据Hu et al.(2015)，设计地形过渡曲线时应该满足两个准则：第一，主要的风场特性(例如平均风速和风攻角)应该与“未受打扰”的来流风一致；第二，过渡段的水平长度尽可能短。明显过渡曲线的设计需要在这两个准则之间寻求一个合理的平衡。

如前所述，收缩段能够在风洞实验段中成功地过渡气流。受此启发，地形过渡段可以由风洞收缩段推导而来。为了更好地阐述推演过程，典型的过渡段曲线被给出如图3，其起点和终点分别是S(0,0)和E(L,H)，其中H是过渡曲线的高度等于图2中收缩段的(H_i-H_o)，L是过渡段的长度。相应地，定义等效斜率和等效倾角分别为K＝H/L andα＝tan^-1(K)。根据图2，图3中横坐标x对应的纵坐标值由H_i-h决定。

由图2和图3的比较，风洞收缩段是可以被用于地形过渡段的。由方程(1)可知，双三次过渡曲线(BCTC)可由下式确定：

其中曲线拐点横坐标值x_f为0.5。类似地，五次过渡曲线被给出：

y＝H[10(x/L)³-15(x/L)⁴+6(x/L)⁵] (5)

因为风场模拟时，天空高度是足够大的，故对维氏曲线收缩段，假设H_i→∞，这样可以得到适用于地形模拟的维氏过渡曲线(WTC)：

其中α为评估变量，上式取α＝1/3。注意到BCTC和WTC曲线的形状分别是可以随x_f和α的值而改变。

为了与由风洞收缩段推导得到的过渡曲线比较，由Hu et al.(2015)提出的理想过渡曲线 (STC)的表达式也给出如下：

其中m和r由以下两个式子获取：

注意到求解过程需要进行迭代。然后这四条曲线形状被显示在图4。

根据上述的讨论，将改变x_f和α的值后，BCTC和WTC的形状被展示于图5。其中Wα表示拥有α值的WTC曲线，Bx代表拥有拐点横坐标值为x_f＝0.1x的BCTC曲线。

步骤2：二维过渡曲线的过渡效果对比。

通过二维CFD模拟对上步骤中得到的过渡曲线进行评估，其计算域和边界条件显示于图 6。对比时，所有的过渡曲线均采用等效斜率K＝0.5774。

CFD数值模拟使用的卡迪尔坐标下的RANS方程如下：

其中这一项需要被数值模拟来封闭上述方程。数值模拟采用的是由Menter(1993)提出的 SST k-ω湍流模型。而k和ω由下式(Fluent 2016)确定：

其中G_k表示湍动能k的产生；G_ω表示ω的生成；Y_k和Y_ω分别代表由于湍流导致的k和ω的耗散；S_k和S_ω是用户自定义源项；D_ω代表扩散项。

数值模拟采用二阶压力插值，k,ω和动量项均采用二阶迎风格式，Chorin(1968)提出的 SIMPLEC(Semi Implicit Method for Pressure Linked Equation Consistent)算法被使用于压力- 速度耦合中

图6所示的计算域中，过渡曲线的起点被定为坐标原点。计算域的高度和长度分别为50H 和71L。其他具体的尺寸可详见图6。相应地，可以得到其阻塞比为2％，这小于Bekele(2002) 建议的3％。出口是压力为0的压力出口；计算域顶面是对称边界；入口的平均风剖面由下式给出：

其中U_z是z高度处的纵向风速，U_b＝10m/s是参考高度Z_b＝L/100所对应的风速。此外，地面边界条件设置为无滑移壁面。

近壁面的首层网格高度设为L/n(n＝1000)，壁面垂直方向的网格高度增长率为1.05。主计算域的总网格数量约为2.5×10⁶，这种网格设置定义为G1。针对网格独立性的检验，一种细化的网格方案定义为G2。而G2和G1唯一的区别在于n＝1500，其总网格数量为3.8×10⁶。而计算结果显示G1和G2的区别十分小，故后面的计算均采用G1网格方案。图7所示为是网格示意图。

后面的讨论可得到W50拥有更为良好的过渡效果，故图8显示了W50曲线在不同位置处的水平风速剖面与入口处的对比。其中，实线为不同位置处的水平风速剖面，虚线为入口处的风速剖面。要注意的是L在下面的讨论中设为单位长度。在x＝-2处，近地面的水平风速开始减小，并且在x＝0处这种趋势更为明显。然而，在x＝1处水平风速开始有明显的增大，且随着x的增大，这种增大的幅值在高度方向上趋于均匀。相应地，水平风速剖面的形状逐渐变得跟入口风速剖面一致，除了数值上稍微有点增大。

通过上述讨论可知，通过过渡曲线后的风速在近壁面有所增长。为了评估这种加速，在 x处的平均加速比被定义为：

其中N是地面至临界高度H_w的网格数量，临界高度H_w是风速变化小于5％的高度；v_i是x位置距地面第i个网格的水平风速；v_0i是入口处距地面第i个网格的水平风速；Δh_i是第i个网格和第i-1个网格的高差。图9是不同过渡曲线的平均加速比对比。可以看出，平均加速比的变化趋势均为先增长到最大值再逐渐减小。具体来看，当α值较大时，Wα过渡曲线的平均加速比上升最快且也最早开始下降至更低值。

与水平风速类似，沿程风攻角也有所变化。在x处的平均风攻角定义如下：

其中a_i是x位置处距地面第i个网格风攻角值。图10是沿着过渡曲线的平均风攻角分布。同平均加速比类似，当α值较大时，Wα过渡曲线在风攻角方面也表现更好。

总的来说，就稳定加速和减小风攻角这两方面，当α值较大时，Wα过渡曲线有更好过渡效果。然而，如果α值过大，气流可能在曲率过大的区域发生分离，影响流场稳定性。为了将这一因素考虑进去，在x位置处超过入流处的平均额外湍流强度定义为：

其中T_i为x位置处距离地面第i个网格的湍流强度；T_0i为距离地面第i个网格的入流湍流强度。平均额外湍流强度沿不同过渡曲线的分布展示在图11中。可以看出，当α达到一个较大值(例如80)时，湍流强度突然增大。其原因可能是在过渡曲线某个位置发生了分离，导致了湍流强度的增大。

步骤3：通过对过渡曲线的评估提出新的过渡曲线。

通过步骤2利用二维的CFD模拟对于一系列过渡曲线进行比较，拥有较大α值的过渡曲线表现较好，然而α不能过大。为了综合评估这些过渡曲线的过渡效果。一个综合指标(CI) 被提出，如下所示：

其中S_j,A_j和T_j分别是j过渡曲线的平均加速比，平均风攻角和平均额外湍流强度值；S_r,A_r和 T_r分别是参考过渡曲线的对应的值。这里将STC曲线作为参考过渡曲线。表1总结了在x＝1 处每条过渡曲线的CI值。能够看出W50曲线的CI是最小的，这就说明其在地形过渡方面表现相对更好。

表1不同过渡曲线在x＝1位置处的CI值

通过综合分析，W50过渡曲线有着更好过渡效果。此外，W50的三项指标的变化率在x＝0.6 位置处均小于5％，详见图9～11。而根据过渡段设计的第二条准则，过渡段的水平长度应该尽可能短。综上，W50过渡曲线在x＝0.6位置处被截断。在L取单位长度为1的情况下，为了保证截断后曲线的形状，取新的等效斜率K_o＝K/0.6＝0.9623，故更短的过渡曲线表达式如下：

注意的是，因为等效斜率被维持为常数，故过渡段的水平长度将随着地形边缘的高度而变化。本是实施例的建模方法即采用这种等效斜率不变的形式。

下面结合实例，对本发明的实施进行详细的说明。

实施实例：深切峡谷的三维风场CFD模拟

1)现场实测概况

普立大桥是一座主跨628m的悬索桥，坐落于云南省山区，如图12所示。三个高程大致相等的山峰分别位于起东北方，西南方和东南方，并且形成一个“Y”字形的山谷。大桥位置处深切峡谷的高差大约200m左右。显然，桥梁设计使用的风参数是不能直接由设计规范得到的。因此进行了桥位处的风场实测(从2013年2月到2015年10月)

一个50m高的测风塔被直立在大桥附近一个没有障碍的小山上，如图12所示。在测风塔上，以下这些测风计被用于测量风速和方向：5个NRG测风计从10m高的位置每隔10m安装1 个；3个NRG风标从10m高的位置每隔20m安装1个。测风计的测量范围为0～96m/s，采样频率为1Hz。风速和风向每隔10min传输一次。同时，2个超声测风计(Young 81000)安装在测风塔30m和50m高的位置，被用于测量脉动风速。其测量范围是0～40m/s，测量频率为4Hz。图13为测风塔的示意图。图14为50m高度处的风向玫瑰图。

2)建模和求解参数设置

首先，从CGIAR-CSI(Consultative Group on International AgriculturalResearch Consortium for Spatial Information国际农业空间信息研究联合会协商小组)下载到山区的STRM (Shuttle Radar Topography Mission航天飞机雷达地形测量任务)高程数据，利用global mapper 软件将其转化成三维坐标值。接着，利用建模软件将点云建模成地形曲线，数值地形见图15。如图16，将地形模型的边缘分成许多块，a1是地形边界线，a2和a4是所提出的过渡曲线， a3是拟合曲线。因为所提的过渡曲线等效斜率K₀是个常数，等于0.9623，故过渡段的水平长度是变化的。最后，使用a1，a2，a3和a4将三维过渡曲面拟合出来，详见图16所示。

本模拟的计算域高度为10H_d，来流畅度L₁＝5D，尾流长度10D，其中H_d是计算域底部到地形最高点的竖向距离，D为地形模型的最大尺寸。故计算域的尺寸为11D×16D×10H_d，如图17所示。阻塞率近似为3％。

考虑到山区地形的复杂性，在地形附近采用非结构化网格，而外部区域采用结构化网格。为了获得更精确的结果，在地形壁面采用了6层棱柱体网格。并且为了确保数值稳定，在结构化网格和非结构化网格之间也采用了棱柱体网格。这里仍然采用了两套网格方案M1和M2。其中在M1中，地形壁面的首层网格高度取为D/6000，而平行于壁面的首层网格尺寸设为 D/600。棱柱体网格竖向增长率为1.05，而非结构化网格各方向的增长率均为1.2。在M2中，地形壁面的首层网格高度取为D/10000，而平行于壁面的首层网格尺寸设为D/1000，而其他设置同M1。另外，两套网格方案结构化网格的竖向增长率均为1.05，横向增长率为1.2。由于两套网格得到结果相差很小，故本模拟网格方案采用M1，网格总数为7.9×10⁷，如图18 所示。图18(b)为测风塔附近的表面网格，其网格被细化且最小网格为1m。

本模拟的计算域高度为10H_d，来流畅度L₁＝5D，尾流长度10D，其中H_d是计算域底部到地形最高点的竖向距离，D为地形模型的最大尺寸。故计算域的尺寸为11D×16D×10H_d，如图17所示。阻塞率近似为3％。

考虑到山区地形的复杂性，在地形附近采用非结构化网格，而外部区域采用结构化网格。为了获得更精确的结果，在地形壁面采用了6层棱柱体网格。并且为了确保数值稳定，在结构化网格和非结构化网格之间也采用了棱柱体网格。这里仍然采用了两套网格方案M1和M2。其中在M1中，地形壁面的首层网格高度取为D/6000，而平行于壁面的首层网格尺寸设为 D/600。棱柱体网格竖向增长率为1.05，而非结构化网格各方向的增长率均为1.2。在M2中，地形壁面的首层网格高度取为D/10000，而平行于壁面的首层网格尺寸设为D/1000，而其他设置同M1。另外，两套网格方案结构化网格的竖向增长率均为1.05，横向增长率为1.2。由于两套网格得到结果相差很小，故本模拟网格方案采用M1，网格总数为7.9×10⁷，如图18 所示。图18(b)为测风塔附近的表面网格，其网格被细化且最小网格为1m。

3)结果与分析

通过分析超声波风速计测到的数据可得到在30m和50m高度处的竖向风攻角值。同样地，这些数据也可以由CFD模拟获取。图20给出了实测和CFD模拟得到的50m处平均竖向风攻角值(对于实测数据是取10分钟求一次平均值)的概率。可以看出CFD模拟得到的结果跟实测数据是比较符合的。总的来说，平均竖向风攻角主要集中在-10到20度之间，而分布在-5到0 度之间的数据概率是最大的，可以占到将近50％，并且其他竖向风攻角基本上呈均匀分布，分别都在10％左右。

基于5个不同高度处NRG风速计测得的平均风速，平均风速剖面能够被得到。这里采用公式(18)对实测数据进行拟合：

其中U_z是在z高度处的平均风速，是U₁₀在10m高度处的平均风速，α是需要拟合的风速剖面的指数。当然，CFD模拟得到的结果也可以用同样的办法得到风剖面指数值。图21 给出了CFD模拟和实测数据得到的平均风速剖面指数值的概率。对比可知两者之间是非常接近的。

Claims (2)

1.一种用于山区风场模拟的过渡曲线的推导方法，其特征在于，包括以下步骤：

步骤1：根据风洞收缩段曲线推导出地形过渡段曲线；

a)引入风洞收缩段曲线：双三次曲线、双三次曲线和维氏曲线：

双三次曲线是由两条三次曲线由拐点连接起来的，其表达式如下：

其中，h是在轴向距离x处的横截面的半高，H_i是入口截面的半高，H_o是出口截面的半高，L是收缩段的长度，x_f是拐点的轴向距离；

五次曲线的表达式如下：

h＝[-10(x/L)³+15(x/L)⁴-6(x/L)⁵][H_i-H_o]+H_i

维氏曲线的表达式如下：

其中，a为调整参数；

b)定义过渡段曲线起点和终端分别为S(0,0)和E(L,H)；其中H是过渡曲线的高度，相当于收缩段中的H_i-H_o，L是过渡段的长度，相当于收缩段的长度L；

相应地，定义等效斜率为K＝H/L，等效倾角为β＝tan^-1(K)；横坐标x对应的纵坐标值由H_i-h决定；

则将风洞收缩段曲线表达方式用于地形过渡段，双三次过渡曲线由下式确定：

五次过渡曲线表示为：

y＝H[10(x/L)³-15(x/L)⁴+6(x/L)⁵]

对维氏曲线收缩段，假设H_i→∞，得到适用于地形模拟的维氏过渡曲线：

其中，α为评估变量，上式取α＝1/3；

步骤2：通过二维CFD模拟对得到的过渡段曲线进行评估：

a)定义通过过渡曲线后的风速在x处的平均加速比S_x为：

其中，N是地面至临界高度H_w的网格数量，临界高度H_w是风速变化小于5％的高度；v_i是x位置距地面第i个网格的水平风速；v_0i是入口处距地面第i个网格的水平风速；Δh_i是第i个网格和第i-1个网格的高差；

对比x_f取不同值的双三次过渡曲线、五次过渡曲线及α取不同值的维氏过渡曲线中各不同过渡曲线的平均加速比，得出：且当α值较大时，维氏过渡曲线的平均加速比上升最快且也最早开始下降至更低值；

c)定义通过过渡曲线后在x处的平均风攻角A_x为：

其中，a_i是x位置处距地面第i个网格风攻角值；

对比沿所述各不同过渡曲线的平均风攻角，得出：当α值较大时，维氏过渡曲线在风攻角方面表现更好；

d)定义通过过渡曲线后在x位置处超过入流处的平均额外湍流强度T_x为：

其中，T_i为x位置处距离地面第i个网格的湍流强度；T_0i为距离地面第i个网格的入流湍流强度；

对比沿所述各不同过渡曲线的平均额外湍流强度，得出：当α达到一个较大值时，维氏过渡曲线的湍流强度突然增大；

步骤3：通过对过渡曲线的评估提出新的过渡曲线：

引入评估综合指标CI：

其中，S_j,A_j和T_j分别是所述各不同过渡曲线的平均加速比、平均风攻角和平均额外湍流强度值；S_r,A_r和T_r分别是参考过渡曲线的平均加速比，平均风攻角和平均额外湍流强度值；

选取CI值最小的曲线为新的过渡曲线。

2.根据权利要求1所述的用于山区风场模拟的过渡曲线的推导方法，其特征在于，所述CI值最小的曲线为α取50时的维氏过渡曲线，其平均加速比、平均风攻角和平均额外湍流强度值在x＝0.6位置处均小于5％，则过渡曲线在x＝0.6位置处被截断，在L取单位长度为1的情况下，为了保证截断后曲线的形状，取新的等效斜率K_o＝K/0.6，原等效斜率K＝0.5774，则K_o＝0.9623，从而新的过渡曲线表达式如下：