何飞龙 许 谦 王 娜 朱春花

(1 新疆大学物理科学与技术学院乌鲁木齐830046)(2 中国科学院新疆天文台乌鲁木齐830011)

(3 中国科学院射电天文重点实验室乌鲁木齐830011)

1 引言

随着望远镜朝着大口径、高工作频段方向发展, 风载荷对望远镜指向精度的影响越来越严重, 天线在进行抗风设计以及保精度观测时更加需要关注台址风场特性. 望远镜台址, 特别是位于山区的台址, 风流动随地形变化产生相对较大差异. 一般为了得到准确的风场信息需要将测风塔布置在离目标区域比较近的位置. 但在望远镜视野范围内不建议有高大物体, 因此, 测风塔位置要与天线区域保持一定距离, 测风塔数据表征天线区域的风场特性就会存在一定误差. 对高性能望远镜指向精度来说, 这种误差不可忽视.因此优选满足台址实际工况条件的测风塔布设位置, 能够将误差尽可能降到最低.

目前, 测风塔的布置大多根据规范进行设计, 需要符合一定要求[1–2], 比如: (1)视野开阔、地势平坦、附近没有障碍物; (2)观测场地附近若有障碍物, 则离障碍物的距离应大于障碍物高度的3倍; (3)测风塔应布置在来风最多的上风口等. 基于规范设置测风塔的目的是为了保障测量数据具有区域代表性. 但当前的规范指标比较宽泛, 存在以下不足: (1)当前的规范主要依据气象观测工程经验总结而来, 气象测量可以在当地选择合适的位置, 但风工程测量目标区域小, 偏重测量小区域风场, 基于规范可选择位置较窄;(2)当测风塔受到障碍物的影响时会产生伪数据, 规范对于远离障碍物的技术指标不够详细和明确, 有的测风塔安装手册指出离障碍物的距离应大于障碍物高度的8–10倍, 若不能满足上述要求, 规范没有给出解决方法, 需依赖设计人员的经验; (3)对于测风塔的布设方案效果也没有给出可行的评估方案, 为了实测可靠性, 技术人员通常会设置多点测量.

随着计算机技术的发展, 国内外不少学者利用数值模拟开展山地风场研究, 并在计算域、边界条件设置、网格结构等方面取得了一些进展[3–5]. 还有学者利用数值模拟分析风速传感器布置的最优位置[6]. 数值模拟的核心是计算流体动力学(computational fluid dynamics, CFD), 即把连续的物理量场离散成一系列有限的变量点, 基于控制方程建立这些离散点的代数方程组, 求解代数方程组获得流场变量的近似解. 它可以仿真接近真实地形上的流场结构, 并借助计算机图形学技术展示流体运动过程[7]. 利用数值模拟优选台址测风塔位置, 能够弥补传统方法难于定量评估的不足, 对提升大口径、高精度望远镜的观测性能起到重要支撑作用. 本研究以新疆奇台110 m口径射电望远镜(QTT)台址为目标, 通过建立台址地理模型开展风场仿真, 借助台址前期布置的测风传感器和气象站所采集的数据, 验证仿真后风场结果的准确性; 依据测风塔布置规范要求布设初选点, 基于风场仿真结果, 提取拟设点和天线位置的风场数据, 通过分析获取测风塔优选位置. 此法还能够为已有望远镜台站精确获取风场特性, 提供补充测风塔或测风传感器布设位置的参考.

2 台址风场数值模拟

2.1 QTT台址介绍

QTT台址位于新疆奇台县半截沟镇石河子村, 坐落在天山东段北麓一个四面环山的盆地内[8], 如图1. 盆地东西长约1.5 km, 南北长约2 km, 海拔约1730–1830 m, 地势南高北低. 盆地四周山体海拔高度1860–2250 m. 该地区常年的主导风向为南风. 在园区内设有一个测风塔, 高60 m. 从地面5 m开始, 每5 m装有一个2D超声波风速风向传感器.园区内还有一个15 m高的气象站. 图1中, F为测风塔、T为天线、Q为气象站.

2.2 地理模型和网格划分

以台址为中心, 划出一个1.4×1.7 km2的矩形区域. 计算域以该矩形为核心区域向四周延伸, 如图2所示. 其中上游边界延伸到盆地外, 左、右和下游边界按上游距离的比例延伸. 计算域顶高设置为5倍地面障碍物的高度[9]. 以正北为0◦方向, 按顺时针方向旋转, 分别模拟0◦、90◦、180◦、270◦方向的风场. 地理数据来自地理云空间(http://www.gscloud.cn)下载的ASTER GDEM (advanced spaceborne thermal emission and reflection radiometer global digital elevation model)数据, 地图的空间分辨率为30 m, 台址区域山体坡度大部分小于45◦, 地图垂直分辨率优于15 m, 可满足本文工况的风场数值模拟需求. 计算域内有许多不规则的山体, 而且核心区域的西侧有一条河谷, 使用结构化网格划分实现比较困难, 因此本文在仿真中将采用自适应比较良好的非结构化网格进行划分, 网格距离最大尺寸设置为32 m. 由于地面对湍流有显著的影响, 地面边界层的处理质量会直接影响近地面风场的精度, 因此采用幂律方式沿地面向上设置边界层网格.第1层高度设置为5 m, 设置8层, 每层增长率指数1.192, 总高为80 m.

图1 QTT台址. F为测风塔、T为天线、Q为气象站.Fig.1 QTT site. F is wind tower, T is antenna, and Q is meteorological station.

图2 风场计算域示意图Fig.2 The diagram of wind field calculation region

2.3 理论模型

本文采用比较成熟的雷诺平均方法进行数值模拟. 该方法对湍流特征量做时均化处理, 将瞬时变量划分为平均变量和脉动变量, 并用雷诺应力项来表示脉动量, 引入湍流模型进行求解, 计算量小, 比较适合工程计算, 个人4核小型工作站基本能够满足计算需求[10]. 在风场模拟中通常把空气假设为低速、不可压的粘性牛顿流体, 本文湍流模型选择适合上述工况且应用比较广泛的可实现k-ε模型. 它是经典的标准k-ε模型的改进模型,推出较新, 与实际流动的物理情况更加符合.

基本控制方程如下:

式中,t是时间;ρ是空气密度;v是速度;Sm是质量源项;p是静压;是应力张量;ρg是重力;F是其他力源项.

2.4 边界条件设置

风场计算域边界条件的设置见表1, 风速剖面一般用对数律或指数律表达式来描述,指数律和对数律两者计算差别很小, 但指数律比对数律计算更加方便, 在工程计算中常采用指数律. 边界入口风速剖面公式如下(风速值为标量):

根据台址地貌, 属于B类, 地面粗糙度α指数取0.15;z0取10 m高度; 设参考速度v0为5 m·s−1.

表1 风场计算域边界条件Table 1 The boundary condition of computing domain

来流湍流强度经验公式采用与本文工况比较接近的日本规范, 取相应的II类地貌参数:

I为湍流强度;zg为梯度风高度, 取350 m.

本文使用Fluent软件, 以上边界条件均可以通过用户自定义函数(User-defined function, UDF)实现.

3 仿真与实测数据比较

实测数据为2017年4—9月的测风塔以及气象站15 m高度的风观测数据. 测风塔数据采集设备是2D超声波风速风向传感器和CR3000微型采集器, 采集频率为每分钟定采一次, 采集时长为1 s; 气象站是ZQZ-A型自动气象站, 其中的风速风向传感器是三杯式风速传感器和单翼式风向传感器, 采集频率为每分钟定采一次, 采集时长为3 s. 将风观测数据处理成风向玫瑰图, 如图3所示, 可知主流风向的来风趋势大致相同, 由于所处位置地形的影响来风频率和大小会有所差异.

图3 风向玫瑰图Fig.3 The wind rose diagram

3.1 实测数据处理

为了使测风塔和气象站的数据能够与数值模拟结果进行比较, 需要将两点的风速归一化处理[11–12]: (1)分别计算测风塔和气象站全部来风的速度平均值i. (2)对测风塔和气象站0◦、90◦、180◦、270◦方向的来风求速度平均值然后用比得到各风向上的风比系数如(5)式. (3)在各风向上用测风塔的风比系数比气象站的风比系数得到两点的风比系数如(6)式. 计算数值模拟中测风塔与气象站的风比系数如(7)式, 与实测风比系数进行比较, 检验数值模拟的可靠性.

式中i= 1、2, 分别对应测风塔、气象站;j= 1、2、3、4, 分别对应0◦、90◦、180◦、270◦; 上标m代表实测、s代表仿真.

3.2 仿真与实测比较

图4为仿真与实测比较结果, 可以看出: 在0◦、90◦和180◦方向上, 数值模拟的风比系数与实测的风比系数符合良好, 变化趋势相同; 但在270◦方向上, 两者的风比系数差值略大, 变化趋势不同. 分析原因, 可能是由于该风向上来风频率极少, 如图3, 因此在数据统计上存在较大的随机性. 由图4中的误差棒可以看到, 误差与来风频率相关, 来风频率越少, 误差越大. 根据理论分析该方向上5 m·s−1来风, 由于测风塔被障碍物遮挡更严重, 实测中测风塔测得的风速要远小于气象站, 而仿真结果符合这一趋势. 综上分析, 本文对QTT台址风场环境的仿真基本符合实际情况, 能够用于测风塔的位置优选.

图4 仿真与实测的风比系数Fig.4 The wind ratio of simulation and measurement

4 天线外围测风塔的布置

假设在距离天线150 m左右的范围设置测风塔不会对天线的视野产生干扰. 拟在天线的东南西北4个方向, 各挑选出1个位置相对较优的设置点. 在图5中可以看到天线外围有许多凹凸不平的小沟壑, 需要将测风塔布置在相对比较平坦且不受周围障碍物遮挡的地方, 初选了拟设点P1、P2、P3、P4.

图5 测风塔位置初选Fig.5 Primary wind tower location

地面气象观测规范指出风传感器布置离地高度一般不小于10 m, 许多台址园区气象站风传感器离地15 m高, 风塔高度也在50 m左右, 所以本文主要分析离地15 m和50 m的风况. 提取仿真结果的数据, 绘制相对地面高度15 m和50 m处的风速云图, 如图6. 由图可知气流受局部地形影响, 相邻区域风速等值线变化比较大. 随着高度上升, 受局部地形影响减弱, 风速等值线逐渐平滑, 气流相对稳定. 因此设置较高的测风塔, 可以减小天线和测风塔之间由于距离所造成的误差.

图6 风速云图Fig.6 The nephogram of wind speed

提取4个拟设点和天线位置的风速值做相对地面15 m和50 m的折线统计图, 由图7中可以看到P3点的风速随风向变化的趋势与天线位置风速变化的趋势最为接近, P4点次之; P2点风速变化趋势与天线位置不同, 但两点的速度值最接近. 随着高度增加, 各点的风速变化趋势大致相同, 但离散程度在减小. 表2为计算拟设点与天线位置风速的均方根误差值(root-mean-square error, RMSE). 随着高度增加, P2、P3和P4点的风速RMSE都在降低, P1点没有变化. 在拟设点中, P3点在15 m高度误差值最大, 随着高度增加, 误差值有所降低, 但也比较大; P2点一直是误差最小的位置. 如果P3点减去固定差值, 则误差最小. 分析P3点误差产生的原因: 从图7中可以看到, 除P3点外, P4点与T点的风速变化趋势也比较接近, 而P4和P3点两者距离比较近, 地形相似, 一定程度上具有地形相关性, 在周围地势的共同作用下导致了P3点产生固定偏差的可能性比较大. 地图精度导致的可能性比较小, 因为所有点的数值模拟结果都在同一精度的地图上产生, 导致单一点偏差的概率很低. 以上分析P2和P3点设置测风塔都比较理想. 测风塔布置中, 如果设置在P3点还需要实际测量P3点和天线位置的风速来消除固定误差. 另外南方来风频率比较高, 测风塔设置在P2点更符合布置在上风口的规范. 综上分析P2点设置测风塔最优.

图7 拟设点与天线位置的风速误差. (a)相对地面15 m; (b)相对地面50 m.Fig.7 The wind speed error between the proposed point and the antenna position. (a) 15 m above ground; (b) 50 m above ground.

表2 拟设点与天线位置的风速RMSETable 2 The wind speed RMSE between the proposed point and the antenna position

5 结论

本文提出了利用数值模拟定量分析射电望远镜台址测风塔最优位置的方法. 基于该方法以QTT台址为研究对象, 开展了台址测风塔布设位置的优选, 初设了4个设置点, 通过仿真与分析, 相对于其他拟设点P2与天线位置的风速误差最小, 最为理想. 本文的数值模拟依据台址地形相应的规范参数设置边界条件, 得到的仿真结果精度满足测风塔布置的需求. 基于规范参数设置, 一定程度上降低研究测风塔布置者的工作量, 方便相关技术人员快速使用.