CN106875488A - 一种反射面天线面板风压系数数值模拟方法 - Google Patents

Abstract

本发明公开了一种反射面天线面板风压系数数值模拟方法，包括：首先采用ProE软件建立天线面板几何模型；然后将几何模型导入前置划分网格软件ICEM‑CFD，建立计算流场区域范围，进行流场区域网格划分；其次将网格文件导入流场求解软件FLUENT，通过一系列设置，计算求解，得到天线面板风压分布；最后将FLUENT的计算求解结果导出，编制MATLAB程序求出分块平均风压系数值。本发明解决了目前反射面天线缺少相关风洞试验数据的问题，通过软件仿真方法，确定多个工况下天线面板的分块风压系数值。方法可提供可视化的流场信息，具有成本低、可重复操作、便于控制的改变模拟环境条件、结果采集全面的特点。

Description

一种反射面天线面板风压系数数值模拟方法

技术领域

本发明涉及CFD数值仿真领域，特别是一种反射面天线面板风压系数数值模拟方法。

背景技术

在大型高精度天线系统设计中，天线反射面的表面精度是主要技术性能指标之一，它直接影响着天线的电性能。随着天线口径的增大，工作频率的提高，对天线反射面精度要求也越来越高。对于露天工作的天线，风荷载对天线工作有着较大影响。在风荷载作用下，天线反射面会偏离的设计曲面，从而使天线的系统性能，如天线增益、天线副瓣和天线交叉极化等主要指标下降，严重失去数据真实性。对于巨型射电望远镜结构，由于其反射面面积巨大、实际工作工况多样(俯仰角、风向角的变化等等)，风荷载的作用更是不可忽略。故对巨型射电望远镜结构在风场中风压分布规律进行了解十分重要。对于大口径射电望远镜，在进行抗风设计时，缺少相关的风洞试验数据，且在现行规范中找不到合适的风载体型系数。因此，研究大口径射电望远镜结构的风场绕流特性和风压分布规律非常重要。

研究射电望远镜天线结构风载荷特性主要有以下四种方法：

1.理论分析

对于大口径巨型射电望远镜天线结构，不仅受风环境的来流湍流影响，由于体型巨大，巨型反射面结构的钝体绕流现象尤为显著，撞击、分离、再附、尾流作用的特征湍流影响不可忽视，无法将如此复杂多样的空气流动完全用解析方法来确定。

2.现场实测

a.实测周期长、受仪器观测精度的影响；

b.由于不能对气候、地形等环境条件加以改变和控制，故有很大的局限性；

c.安装数量庞大的传感器，对射电望远镜天线结构本身和射电望远镜的观测精度的影响已不容忽视。

3.风洞试验

由于望远镜结构的反射面相对较薄，布置测压管很不方面，导致增大模型厚度造成尺寸效应；此外风洞试验技术中的动力相似条件很难满足，也必将带来试验误差。

4.数值模拟

以CFD为基础的数值模拟方法具有成本低、可重复操作、便于控制的改变模拟环境条件、结果采集全面、可提供可视化的流场信息；与风洞试验方法相比，CFD数值模拟方法可不受模型尺度影响，从而克服了风洞试验中雷诺数相似难满足的问题。

发明内容

基于上述问题，本发明的目的在于解决目前反射面天线缺少相关的风洞试验数据的问题，通过数值模拟方法，确定多个工况下天线面板的分块风压系数值。

实现本发明目的的技术解决方案是:

一种反射面天线面板风压系数数值模拟方法，包括以下步骤：

(1)按照反射面天线面板原型尺寸建立几何模型，采用Proe三维绘图软件建立指定反射面天线面板口径和焦径比的旋转抛物面；

(2)将Proe三维绘图软件建立的几何模型导入ICEM-CFD软件中，建立流场区域；采用非结构化网格划分流场区域，导出为网格文件；

(3)检查非结构化网格质量，若非结构化网格质量大于0.2，最大偏斜率低于0.9，则进入下一步；否则跳转至步骤(2)重新进行网格划分；

(4)把步骤(2)生成的网格文件导入到FLUENT软件中，通过选择湍流模型、定义边界条件，进行计算求解得到反射面天线面板风压节点的风压系数；

(5)将反射面面板节点风压系数数据以ASCII形式导出，再导入到MATLAB数据处理软件中；根据已有风洞风压系数的面板分块样式，通过编制MATLAB程序，求出面板分块区域的平均风压系数；

(6)重复步骤(1)～(5)，得出俯仰角E＝0°、方位角A＝30°，俯仰角E＝0°、方位角A＝60°这两个工况下的软件仿真风压系数值；采用软件仿真得出的风压系数值与已有风洞风压系数值进行对比，验证软件数值仿真天线面板风压系数的可行性。

进一步，所述步骤(1)中，采用Proe三维绘图软件建立指定反射面天线面板口径和焦径比的旋转抛物面的俯仰角E＝0°，方位角A＝0°。

进一步，所述步骤(2)中，建立流场区域通过下述方法实现：

2a)计算流场区域的设置包括上游尺寸L₁，下游尺寸L₂，计算域迎风面宽度B及高度H，计算模型的尺寸为l×b×h，确定计算域设置参数；

2b)设定流场区域阻塞率ε不大于3％：

式中，D是天线口径，B是计算域迎风面宽度，H是计算域迎风面高度。

进一步，所述步骤(2)中，采用非结构化网格划分流场区域，在近壁处生成边界层加密区域，模型表面近壁处至少生成10层边界层网格，网格由密到疏的过渡比率为1.05-1.2。

进一步，所述步骤(4)中，天线面板风压节点分布，通过下述方法实现：

4a)选择湍流模型为雷诺应力模型：

其中，ρ为密度，u'_i和u'_j为速度脉动分量，t为时间，C_ij为对流项，D_T,ij为湍动扩散项，D_L,ij为分子粘性扩散项，P_ij为剪应力产生项，G_ij为浮力产生项，φ_ij为压力应变项，ε_ij为粘性耗散项，F_ij为系统旋转产生项；

4b)定义流场区域边界条件，包括流域入流处采用的FLUENT中的速度进口条件、流域出口边界采用的完全发展出流边界条件、区域顶部和两侧采用的对称边界条件和区域底部和面板模型表面采用的无滑移的壁面条件；

4c)选用SIMPLEC算法，风压系数C_p的计算公式如下：

其中p为在FLUENT中成为表压，p_ref为在软件中设为0，v为参考位置处的风速。

进一步，所述步骤4b)中，计算流域入流处采用FLUENT中的速度进口条件(velocity-inlet)，包括计算域入口处的湍流强度I、湍流动能k和湍流耗散率ε，具体表达式如下：

l＝100(z/30)^0.5

其中z、分别是流域中任意高度和对应的平均风速，l为湍流积分尺度，α为地面粗糙度。

进一步，所述步骤4b)中，出口边界条件采用流场任意物理量ψ沿出口法向梯度为零，的完全发展出流边界条件。

进一步，风压系数C_p为望远镜结构上下表面的平均压力系数的差值：

C_p＝C_{p_up}-C_{p_down}

式中C_{p_up}为迎风面风压系数，C_{p_down}为背风面风压系数，C_p为风压系数。

进一步，所述步骤(6)中，风压系数C_p值与已有风洞风压系数值进行对比仿真误差通过下式得到：

式中C_Pi为第i块软件仿真的风压系数值，c_pi为第i块风洞实验得出的风压系数值。

本发明与现有技术相比，具有以下特点：

本发明采用数值模拟方法，得出天线面板分块风压系数值，以CFD为基础的数值模拟，可提供可视化的流场信息，具有成本低、可重复操作、便于控制的改变模拟环境条件、结果采集全面的特点。与风洞试验方法相比，CFD数值模拟方法可不受模型尺度影响，从而克服了风洞试验中雷诺数相似难满足的问题。

通过对比三个工况下数值模拟方法和风洞试验得出的分块风压系数值，得出软件仿真得到的结果基本与试验数据吻合的结论，进而可以丰富多个工况下的风压系数值，解决了除风洞试验数据试验的工况外缺少风压系数的问题，为天线抗风设计提供数据支持。

附图说明

图1是本发明方法的流程图；

图2是本发明使用的软件之间的关系流程框图；

图3是计算域尺寸设置项目；

图4是ICEM-CFD建立计算域示意图；

图5是面板附近网格示意图；

图6是面板分块样式；

图7是风洞试验(E＝0°，A＝0°)得出的风压系数；

图8是软件仿真与风洞试验(E＝0°，A＝0°)风压系数对比；

图9是软件仿真与风洞试验(E＝30°，A＝0°)风压系数对比；

图10是软件仿真与风洞试验(E＝60°，A＝0°)风压系数对比。

具体实施方式

下面结合附图及具体实施案例对本发明做进一步说明。

参照图1，一种反射面天线面板风压系数数值模拟方法，其包括以下步骤：

(1)按照反射面天线面板原型尺寸建立几何模型。

采用原型尺寸能避免尺寸效应可能对结果所产生的影响。由于目前的网格生成技术，对于被复杂的空间桁架结构分割的具有尖角位置的复杂空间流场区域，无法进行有效的网格划分进行数值求解，而且背架结构并非对风场造成影响的主要因素，因此不对其进行考虑。

采用Proe三维绘图软件建立指定反射面天线面板口径、焦径比的旋转抛物面(俯仰角E＝0°，方位角A＝0°)，因为已有风洞数据中天线的焦径比是0.5，为了将软件仿真的结果和已有风洞数据进行对比，因此采用Proe三维绘图软件建立口径26m，焦径比f/D＝0.5的天线旋转抛物面模型，俯仰角E和方位角A均为0°，即风正吹天线。

(2)计算流域建立和网格划分。

1)计算流域的确定

位于大气边界层中的天线面板结构风对其绕流时，处于一个完全开口的流动风场中，但风对天线面板作用的影响具有一定的范围，在数值模拟时可给定有限的三维计算区域，并确定计算区域各边界的边界条件以模拟实际的流动风场。计算流域的确定需要考虑，随着计算区域的增大，网格数量增加，计算量加大，计算时间增长；然而计算区域过小，则可能导致模拟失真。合理的选择计算区域有助于模拟结果准确性的提高和计算量的减少。

将Proe三维绘图软件建立的几何模型导入ICEM-CFD软件中，建立流场区域，使用的软件之间的关系流程框图见图2所示。

建立流场区域具体步骤包括：

①计算流场区域的设置包括上游尺寸L₁，下游尺寸L₂，计算域迎风面宽度B及高度H，计算模型的尺寸为l×b×h，计算流场区域尺寸如图3所示，确定计算域设置参数；确定计算域设置参数如下：

L₁＝5l,L₂＝15l,B＝12b,H＝8h

考虑到巨型射电望远镜的多种工况情况(不同俯仰角、方位角的变化),为实现对计算域的统一设置，选取巨型射电望远镜旋转抛物面反射面的口径D作为特征尺寸，并令l＝b＝h＝D＝26m，则：

L₁＝130m,L₂＝390m,B＝312m,H＝208m

②由于阻塞效应应对结构周围流场造成干扰，得不到正确的模拟结果，因此通常要求阻塞率不大于3％:

其中D为天线口径，B为迎风面宽度，H为迎风面高度。

经计算阻塞率ε＝0.818％≤3％满足要求；天线反射面结构底部最靠近地面的高度在3m-8m，本例取4m，就是说天线反射面中心高度为17m，天线反射面最高处为30m。建立的计算域如图4所示。

2)网格单元的选择

网格的类型和网格的数量对计算结果的准确性具有很大的影响。非结构化网格因其对模型的自适应性好，可解决较为复杂流场的网格划分，得到了广泛的应用。本例选用非平衡壁面函数方法来处理近壁区域的流体运动。在近壁处生成边界层加密区域。模型表面近壁处至少生成10层边界层网格。网格由密到疏的过渡比率以1.05-1.2之间为宜。面板模型附近生成网格的样式如图5所示。

(3)检查非结构化网格。

若网格质量大于0.2，网格最大偏斜率低于0.9，则进入下一步；否则跳转至步骤(2)重新进行网格划分，最终网格质量大于0.3，网格最大偏斜率为0.73。

(4)若步骤(3)满足要求，则导出ICEM-CFD生成的网格文件到FLUENT软件中，进行FLUENT软件设置。通过选择湍流模型、定义边界条件、设置求解控制参数和初始化流场进行计算求解，得到反射面天线面板风压节点的风压系数。具体步骤如下：

4a)选择湍流模型：湍流模型选为求解精度最高的雷诺应力(RSM)模型，选用非平衡壁面函数以模拟壁面附近的复杂流动现象。

雷诺应力输运方程：

其中，ρ为密度，u'_i和u'_j为速度脉动分量，t为时间，C_ij为对流项，D_T,ij为湍动扩散项，D_L,ij为分子粘性扩散项，P_ij为剪应力产生项，G_ij为浮力产生项，φ_ij为压力应变项，ε_ij为粘性耗散项，F_ij为系统旋转产生项；

4b)定义流场区域边界条件：边界条件用于定义在流动进口处的流动速度及相关其它标量型流动变量。

①入口边界条件

计算流域入流处采用FLUENT中的速度进口边界条件(velocity-inlet)。该边界条件适用于不可压缩流动，对于可压缩流动问题时会使得入口处的总温度或总压有一定的波动，导致非物理结果，所以可压缩流问题不适合采用速度进口边界条件。此例为不可压缩流，可采用以velocity-inlet边界。需对流动速度v、湍动能k和耗散率ε定义。计算流域入流处采用FLUENT中的速度进口边界条件(velocity-inlet)。入口处的湍流强度I、湍流动能k和湍流耗散率ε的具体表达式如下:

l＝100(z/30)^0.5

其中z，分别是流域中任意高度和对应的平均风速，z由模型底部开始算起，α＝0.20；l为湍流积分尺度；

本例为B类地貌，平均风速剖面、湍流动能和湍流耗散率值采用FLUENT提供UDF(User-Defined Function，用户自定义函数)编制C语言程序与FLUENT作接口实现。

②出口边界条件

流域出口边界条件采用的完全发展出流边界条件(outflow)：采用流场任意物理量ψ沿出口法向梯度为零，的完全发展出流边界条件。

③计算区域顶部和两侧为对称边界条件(symmetry)。

④计算区域底部和面板模型表面采用无滑移的壁面条件(wall)。

不考虑热交换，屏蔽能量方程，流体介质是空气，具有不可压缩性，密度为常数，材料参数使用缺省值；本例研究天线面板周围的风场，关心的是风对面板表面的压力作用，以天线面板迎风面表面的平均压力值为监测对象，在收敛标准判断上，除使控制方程的迭代残余量达到设定值要求外，还要注意对研究对象表面的压力的监测，当其表面压力基本不发生变化可认为计算已经进入稳定状态。

本例需要将风对天线面板的作用进行定量的分析，所以采用高精度二阶迎风格式配合较小的欠松弛系数，选用SIMPLEC算法。

4c)参考压力位置的确定

对于不可压缩流动，由于边界条件中不包括压力边界条件，必须设置一个参考压力位置。在计算流域中选一个不受模型影响，并且远离流域边界的位置作为参考压力位置，FLUENT在每次迭代后，都将整个压力场减去这个参考压力位置的压力值，从而使得所有的压力均按照参考压力位置的值来度量，压力绝对值常比相对值高几个数量级，若采用压力绝对值进行计算，则会导致计算存在较大的误差。确定参考压力位置根据所建模型和坐标系，选择点(100m，160m，60m)为参考压力点，将此点的参考压力设为0，使表压和绝对压力相等；FLUENT中风压系数C_p的计算公式如下：

其中p在FLUENT中成为表压，因为将参考压力设为0，所以此处也为绝对压力(absolute pressure)，p_ref在软件中设为0，v为参考位置处的风速。因此C_p的取值与参考位置的选取有关，在工程应用中如参考位置选取不同，需要进行换算，换算公式如下：

其中C_Ph为参考高度为h时的风压系数，C_PH为参考高度为H时的风压系数，α为地面粗糙度。

本例在reference values中将速度v设为23.7m/s；本例计算的风压系数是以参考高度为10m处的风速计算的，实际工程上应该以天线反射面最高处为参考高度，本例为30m，采用换算公式将参考高度为10m的风压系数转换成参考高度为30m的风压系数：

C_P30为参考高度为30m时的风压系数，C_P10为参考高度为10m时的风压系数；

从计算域入口设定的初始值进行全局初始化；设定残差值为10^-4，迭代步数为3000步。

(5)将面板节点风压系数数据以ASCII形式导出，转换成excel文件，再导入到MATLAB数据处理软件中。

由于望远镜天线结构反射面各点位置的风荷载作用实际上是前后表面的压力差，故风压系数C_p应为望远镜结构上下表面的平均压力系数的差值：

C_p＝C_{p_up}-C_{p_down}

式中C_{p_up}为迎风面风压系数，C_{p_down}为背风面风压系数，C_p为风压系数。

图7是已有风洞试验得出的天线面板风压系数分布图(正吹)，将天线面板分成5圈×24块＝120块，分块样式如图6所示，编制matlab程序，对每个分块内的所有节点的风压系数值求平均值，以这个值来表示该分块区域的风压系数值，图8是软件仿真和风洞数据对比情况。

(6)已经得出了俯仰角E＝0°，方位角A＝0°，这个工况下的软件仿真风压系数值。重复步骤(1)～(6)得出俯仰角E＝0°，方位角A＝30°、俯仰角E＝0°，方位角A＝60°这两个工况下的软件仿真风压系数值，如图9和图10所示。

表1软件仿真误差

工况	RMSE	MRE
			E＝0°，A＝0°	0.0727	3.99％
E＝0°，A＝30°	0.1230	6.71％
			E＝0°，A＝60°	0.1899	13.34％

风压系数C_p值与已有风洞风压系数值进行对比仿真误差通过下式得到：

式中C_Pi为第i块软件仿真的风压系数值，c_pi为第i块风洞实验得出的风压系数值。

由于已有风洞实验并没有给出地貌类型和计算平均风压系数的参考高度，而且本例简化了方位架、馈源等结构，故数值模拟结果与风洞试验数据之间必然存在一定差异。通过计算这两者之间的差异，得出软件仿真可行的结论。采用数值模拟方法，可以得到较好的望远镜反射面板风压分布规律。可以进一步丰富天线各个工况下的面板风压分布情况。

以上所述，仅是本发明的实施例，并非对本发明做任何限制，凡是根据本发明技术对以上实施例所做的任何简单修改，变更以及等效结构变化，均仍属于本发明技术方案的保护范围内。

Claims (9)

1.一种反射面天线面板风压系数数值模拟方法，其特征在于，包括以下步骤：

(1)按照反射面天线面板原型尺寸建立几何模型，采用Proe三维绘图软件建立指定反射面天线面板口径和焦径比的旋转抛物面；

(2)将Proe三维绘图软件建立的几何模型导入ICEM-CFD软件中，建立流场区域；采用非结构化网格划分流场区域，导出为网格文件；

(3)检查非结构化网格质量，若非结构化网格质量大于0.2，最大偏斜率低于0.9，则进入下一步；否则跳转至步骤(2)重新进行网格划分；

(4)将步骤(2)生成的网格文件导入到FLUENT软件中，通过选择湍流模型、定义边界条件，计算求解得到反射面天线面板风压节点的风压系数；

(5)将反射面面板节点风压系数数据以ASCII形式导出，再导入到MATLAB数据处理软件中；根据已有风洞风压系数的面板分块样式，通过编制MATLAB程序，求出面板分块区域的平均风压系数；

(6)重复步骤(1)～(5)，得出俯仰角E＝0°、方位角A＝30°，俯仰角E＝0°、方位角A＝60°这两个工况下的软件仿真风压系数值；采用软件仿真得出的风压系数值与已有风洞风压系数值进行对比，验证软件数值仿真天线面板风压系数的可行性。

2.根据权利要求1所述的一种反射面天线面板风压系数数值模拟方法，其特征在于，所述步骤(1)中，采用Proe三维绘图软件建立指定反射面天线面板口径和焦径比的旋转抛物面的俯仰角E＝0°，方位角A＝0°。

3.根据权利要求1所述的一种反射面天线面板风压系数数值模拟方法，其特征在于，所述步骤(2)中，建立流场区域通过下述方法实现：

2a)计算流场区域的设置包括上游尺寸L₁，下游尺寸L₂，计算域迎风面宽度B及高度H，计算模型的尺寸为l×b×h，确定计算域设置参数；

2b)设定流场区域阻塞率ε不大于3％：

$\mathrm{\ϵ}=\frac{\mathrm{\π}{\left(\frac{D}{2}\right)}^2}{B\·H}=\frac{{\mathrm{\πD}}^2}{4B\·H}\≤3\%$

式中，D是天线口径，B是计算域迎风面宽度，H是计算域迎风面高度。

4.根据权利要求1所述的一种反射面天线面板风压系数数值模拟方法，其特征在于，所述步骤(2)中，采用非结构化网格划分流场区域，在近壁处生成边界层加密区域，模型表面近壁处至少生成10层边界层网格，网格由密到疏的过渡比率为1.05-1.2。

5.根据权利要求1所述的一种反射面天线面板风压系数数值模拟方法，其特征在于，所述步骤(4)中，天线面板风压节点分布，通过下述方法实现：

4a)选择湍流模型为雷诺应力模型：

$\frac{\∂\left(\mathrm{\ρ}\stackrel{\‾}{u_i^{\′}{u}_j^{\′}}\right)}{\∂t}+C_{ij}=D_{T,ij}+D_{L,ij}+P_{ij}+G_{ij}+{\mathrm{\φ}}_{ij}+{\mathrm{\ϵ}}_{ij}+F_{ij}$

其中，ρ为密度，u′_i和u′_j为速度脉动分量，t为时间，C_ij为对流项，D_T,ij为湍动扩散项，D_L,ij为分子粘性扩散项，P_ij为剪应力产生项，G_ij为浮力产生项，φ_ij为压力应变项，ε_ij为粘性耗散项，F_ij为系统旋转产生项；

4b)定义流场区域边界条件，包括流域入流处采用的FLUENT中的速度进口条件、流域出口边界采用的完全发展出流边界条件、区域顶部和两侧采用的对称边界条件和区域底部和面板模型表面采用的无滑移的壁面条件；

4c)风压系数C_p的计算公式如下：

$C_p=\frac{p-p_{ref}}{0.5{\mathrm{\ρv}}^2}$

其中p为在FLUENT中成为表压，p_ref为在软件中设为0，v为参考位置处的风速。

6.根据权利要求5所述的一种反射面天线面板风压系数数值模拟方法，其特征在于，所述步骤4b)中，计算流域入流处采用FLUENT中的速度进口条件(velocity-inlet)，包括计算域入口处的湍流强度I、湍流动能k和湍流耗散率ε，具体表达式如下：

$k=1.5{(\stackrel{\‾}{u}\·I)}^2$

$\mathrm{\ϵ}=\frac{1}{l}\·{0.09}^{\frac{3}{4}}{k}^{\frac{3}{2}}$

$I=\left\{\begin{array}{cc}0.23& z\≤5\\ 0.1{\left(\frac{z}{350}\right)}^{-\mathrm{\α}-0.05}& 5\≤z\≤350\end{array}\right.$

l＝100(z/30)^0.5

其中z、分别是流域中任意高度和对应的平均风速，l为湍流积分尺度，α为地面粗糙度。

7.根据权利要求5所述的一种反射面天线面板风压系数数值模拟方法，其特征在于，所述步骤4b)中，出口边界条件采用流场任意物理量ψ沿出口法向梯度为零，的完全发展出流边界条件。

8.根据权利要求5所述的一种反射面天线面板风压系数数值模拟方法，其特征在于，风压系数C_p为望远镜结构上下表面的平均压力系数的差值：

C_p＝C_{p_up}-C_{p_down}

式中C_{p_up}为迎风面风压系数，C_{p_down}为背风面风压系数，C_p为风压系数。

9.根据权利要求1所述的一种反射面天线面板风压系数数值模拟方法，其特征在于，所述步骤(6)中，风压系数C_p值与已有风洞风压系数值进行对比仿真误差通过下式得到：

$RMSE=\sqrt{\frac{{\mathrm{\Σ}}_{i=1}^{96}{(C_{Pi}-c_{pi})}^2}{96}}$

$MRE=\frac{1}{96}\underset{i=1}{\overset{96}{\Σ}}\left|\frac{C_{Pi}-c_{pi}}{c_{pi}}\right|$

式中C_Pi为第i块软件仿真的风压系数值，c_pi为第i块风洞实验得出的风压系数值。