大型射电望远镜日照热误差及其补偿的仿真研究

2024-02-02陈浩祥赵武林项斌斌李东伟

工程科学与技术 2024年1期

雷震，宁亮，陈浩祥，赵武林，项斌斌，李东伟

(1.长安大学工程机械学院，陕西西安 710006；2.中国电子科技集团公司第三十九研究所，陕西西安 710016；3.中国科学院新疆天文台，新疆乌鲁木齐 830011)

射电望远镜广泛应用于射电天文、测控导航等领域，可通过特定形状的反射面实现电磁波向馈源的汇聚。然而环境载荷会使其结构产生机械形变，引起反射面形状误差进而降低望远镜的性能[1]。重力、风载对天线性能的影响及其控制方法已较为成熟[2–6]，钱宏亮等[3]研究了上海65m天线在重力及风载作用下，主反射面精度的变化情况。Antebi等[4]分析了天线如何在重力情况下使用副面变形调整保证天线的精度。热对低频、中小口径天线性能的影响并不明显，且可以通过采用天线罩、温控系统等进行控制，但目前筹建的新疆110m大口径全可动射电望远镜（QTT）口径大、工作频率高，且新疆地区日照强烈，太阳辐射所引入的结构热误差问题亟须重视[7–8]。

目前，国内外学者对大口径全可动射电望远镜热行为的研究较少，主要对中小型口径的射电望远镜进行了实测研究和数值模拟。Greve等[9]通过大量实验分析，改进了IRAM–30m天线的温控系统，使天线性能得到明显改善。Bremer等[10]进一步构建了IRAM–30m天线热变形的简化数据模型，实现了热误差实时计算。Ambrosini等[11]将实测温度数据加载到VLBI–32 m天线有限元模型上对天线热变形进行了预测。MacDonald[12]和DiCarlo等[13]针对带天线罩的37 m天线(HUSIR)研究了其太阳/风/温度及电磁性能变化情况。Tsela等[14]针对南非的LLR望远镜建立了其数据驱动的热梯度及精度预测模型。连培园等[15]针对天线日照温度场分布特征提出了温度场实时预估方法。常文文[16–17]、易乐天[18]等建立了南山站25m天线反射面的有限元模型，探讨了天线面板、座架等结构由于太阳辐照引起的非均匀温度场及热致变形。刘泽鑫等[19]对武清70m天线座架也进行了类似研究。钱宏亮等[3]通过数值模拟的方法，研究了上海65m天线在日照温度场作用下，主反射面精度的变化趋势。

日照热在天线结构上产生的温度场具有非均匀性、时变性等特点，使天线热变形呈现出复杂的时空分布以及多因素相关性。热环境、天线姿态等不确定性因素，使得依据经验公式评估天线热变形精度的方法难以满足现代天线的设计要求[20–21]。为保证QTT的正常运行，研究其在日照作用下的热变形及其变迁规律一方面可为其热误差补偿提供参考，另一方面可用于预测未来一段时间内天线热变形的趋势，提前规划调整方案。

本文首先对QTT的热问题进行了建模，对天线在典型夏至日的温度场及变形情况进行了数值模拟，分析了不同情况下反射面精度的变化情况及其规律；最后给出了几种工况下副面的位置补偿量，通过分析补偿前后主反射面的光程差，探讨了天线热误差的时空分布规律及其影响机理，有助于预估、理解天线日照热误差问题，并为热补偿方案提供依据。

1 天线分析模型

1.1 热环境建模

天线与外界的热量传递途径如图1所示，主要包括天线接收太阳辐射、天线与空气的对流换热、天线与地面和远空的辐射换热。换热过程与风速、大气清洁度、地面温度、天空温度等因素密切相关。天线结构在上述因素的共同作用下形成瞬态时变的温度场。

图1 射电望远镜换热简图Fig. 1 Heat transfer diagram of a radio telescope

1.1.1太阳辐射热模型

太阳辐射热模型是用于计算地球上的物体接收到的日照热功率的数学模型，常用的有Dilger模型[9]和ASHRAE晴空模型[22–23]，本文采用ASHRAE晴空模型计算天线接收到的太阳热辐射功率。

ASHRAE晴空模型[23]中地面物体接收的太阳辐射由直接辐射、散射辐射和反射辐射3部分构成，其中穿过大气层直射物体的部分为直接辐射，被周围环境反射后间接照向物体的部分称为反射辐射，被大气中的灰尘等分子散射的部分为散射辐射。

直接辐射强度GND计算方法如下：

式中：A为大气质量为零时的太阳辐射强度；B为大气消光系数； βs为太阳高度角；CN为大气清洁度；θ为太阳对构件表面的入射角，若cos θ小于0则表示构件处于阴影中。

散射辐射强度计算方法如下：

式中，C为散射辐射强度， α为物体表面相对于地面的倾角。

反射辐射强度计算方法如下：

式中：pg为地面反射率，一般取0.3[24]；GtH为落在壁面之前的水平面或地面上的总辐射量,GtH=(C+sin βs)GND。

入射到天线表面的太阳辐射热流密度为：

天线实际得到的太阳热流密度为：

式中， γ为天线表面的太阳辐射吸收率，天线表面一般涂有白漆，故 γ取值为0.3[24]。

ASHRAE模型中A、B、C这3个系数可通过中国辐射强度观测数据拟合得到[25]，CN、 βs、 α 、 θ可根据天线地理位置及构件空间方位计算得到[16]。

1.1.2对流换热模型

天线与环境的对流换热可用牛顿换热定律[26]表述如下：

式中：hc为对流换热系数，其取值Saetta[27]建议使用公式hc=4.0V+5.6，其中V为风速；Tx为构件表面温度；Ts为环境温度，可按下式计算[26]：

式中，Tsmax、Tsmin分别为环境最高温度、最低温度，t为时间。

1.1.3净长波辐射模型

净长波辐射指天线与地面、远空的辐射换热。天线表面得到的净长波辐射可表述为[22]：

式中：ε为发射率；σ为黑体辐射系数；Fwg、Fws分别为构件表面对地面、天空的辐射角系数；Tsky为远空温度；Tg为地面温度，可按下式计算：

式中，E、F、M为常量，取值可见文献[28–29]。

1.1.4瞬态温度场模型

天线结构的热分析遵循热力学第一定律，瞬态热分析的能量平衡方程为[30]：

式中，C为比热矩阵，考虑系统内能的增加，T和T˙分别为温度向量以及其对时间的导数，K为传导矩阵,包含导热系数、对流系数及辐射率和形状系数等，Q为节点热流率向量。

1.2 阴影遮挡与辐射传热系数

1.2.1天线主反射面阴影遮挡

QTT主反射面直径110m，面积巨大，在太阳转动过程中，面板会对自身产生遮挡，明暗区及分界点示意如图2所示。根据光线追踪原理，可以找到反射面边缘点在反射面上的投影点，进而得到反射面阴影区与太阳照射区的分界线，计算公式如下：

式中：N(t,AZ,EL)为太阳入射光线方向向量，是关于时间t、天线方位角AZ与俯仰角EL的函数；X′为反射面边缘点；X0为反射面与太阳入射光线的交点。

图2 反射面太阳照射明暗区及分界点示意图Fig. 2 Self shading effect of reflector

1.2.2辐射角系数

物体之间的辐射换热除了与物体绝对温度有关，还与换热表面之间的相对位置有关，辐射角系数[31]描述了这种关系，两微元表面之间的辐射角系数如图3所示，计算公式如下：

式中：Ib1为常数；Eb1为微元表面相对于黑体的辐射力；θ1、θ2分别为dA1、dA2的维度角；dΩ1为微元表面dA1对dA2所张的立体角；r为dA1到dA2的半径。

图3 微元表面辐射角系数Fig. 3 Surface radiation shape factor

1.3 主反射面精度分析及副反射面补偿方法

图4为天线变形分析示意图，变形后其反射面形状不再是理想抛物面，电磁波不再完全汇聚于原设计焦点[32]（图4中F点），为缓解上述问题，大口径射电望远镜基本都采用“保型”设计，即要求各个姿态下重力变形后的反射面形状皆尽可能接近参数不同的抛物面[1]（图4中虚线），这样将馈源移动到新的抛物面的焦点上（图4中F'点）后，能够较好恢复原来的天线性能。图4中uA与wA为顶点的两个方向位移，上述虚拟抛物面被称为最佳吻合抛物面[32]，变形后的反射面相对于最佳吻合抛物面对应点之间半光程差的均方根误差值（RRMS）作为衡量反射面精度的指标。

图4 天线变形分析示意图Fig. 4 Diagram of the antenna deformation analysis

图5为设计抛物面经过位姿刚性变换及焦距变换得到最佳吻合抛物面的过程。OXYZ及O1X1Y1Z1分别为原设计抛物面、最佳吻合抛物面的坐标系，原点O与O1分别为它们的顶点，OZ与O1Z1分别为它们的焦轴，任意点在两个坐标系下坐标分别为（x,y,z）和（x1,y1,z1），设计抛物面方程为：

最佳吻合抛物面方程为：

式（14）～（15）中，f、h分别为设计抛物面的焦距及焦距的变化量。

图5 刚体位姿变换示意Fig. 5 Rigid body pose transformation

由图5可知，在旋转量十分微小情况下，其线性化坐标转换方程为：

式中，uA、vA、wA为最佳吻合抛物面相对于设计面的平移量， ϕx、 ϕy、 ϕz为相对各轴的旋转量，上述6个参数为待定值。

将式（16）带入式（15），略去其2阶微量，可得到最佳吻合抛物面在原设计面坐标系下的方程：

因为抛物面为圆周对称结构，故反射面绕焦轴转动量 ϕz对天线性能无影响，可以省略，最终得到含6个参数的最佳吻合抛物面方程。

变形后的抛物面相对于最佳吻合抛物面的误差如图4所示，略去2阶微量，得到逐点误差关于上述待定参数的线性表达式如下：

式中，x、y、z为反射面节点坐标，u、v、w分别为节点对应的位移，上述6个参数为已知量。

最佳吻合抛物面的参数应使 ∆2达到最小值，上述问题为线性最小二乘问题；对 ∆2求关于6个待定参数的偏导，并令其为零，即得到关于最佳吻合抛物面参数的线性方程组：

式中，A为节点的初始位置信息矩阵，X为最佳吻合抛物面的参数向量，d为变形后反射面的节点位移信息矩阵。

则以均方根值为度量的反射面精度的计算方式如下：

式中，N为背架上弦节点数量[1]，RRMS反映了反射面的形状偏离抛物面的程度。

天线主面热变形的精度评价及副面补偿计算流程如图6所示，即根据反射面的热变形数据，由式（19）（20）确定与之对应的最佳吻合抛物面参数及天线热变形型面精度（RRMS），将副反射面的焦点调整到上述新的主反射面焦点位置（见图4）即可部分实现天线误差补偿，RRMS越小则说明热变形后反射面形状越接近于抛物面，补偿效果越好。

图6 精度评价及副面补偿量计算流程Fig.6 Accuracy evaluation and subreflector repositioning

副反射面的位置调整量与最佳吻合参数的关系如下：

式中，UX、UY、UZ分别为副反射面3个方向的位置调整量。

2 结果与分析

2.1 天线有限元模型及分析流程

QTT反射体的有限元模型如图7所示。天线桁架材料为钢，反射面材料为铝，材料的热特性见表1。天线热力耦合分析流程如图8所示。

图7 天线反射体模型Fig. 7 Antenna reflector model

表1 材料参数Tab.1 Material parameters

图8 天线热力耦合分析流程Fig. 8 Thermal-mechanical coupling analysis process

2.2 热误差分析及补偿

日照温度变化具有一定的时间尺度，本文研究天线在季节性均匀温度变化与逐日温度变化两种情况下的响应情况。

2.2.1季节性热误差分析

天线反射体通过俯仰轴与方位座架连接（如图7所示），根据轴承连接处是否允许反射体沿轴承轴向移动，本文选取两种极端的俯仰轴约束方式：全约束（FC）、对称浮动约束（SFC）。

表2为奇台县季节温度统计数据。由表2可知，QTT台址四季环境温度在–20℃～30℃范围内变化。以天线俯仰角90°、参考温度20℃为例，分析天线结构在不同约束方式、不同环境温度条件下，反射面精度（RRMS）的变化趋势，计算结果如图9所示。

表2 奇台县四季温度Tab.2 Seasonal temperature in Qitai County

图9 不同工况下的R RMSFig. 9 RRMS under different working conditions

由图9可知，俯仰轴的约束方式对反射面精度影响很大，两种约束方式下的RRMS变化趋势相同，但对称浮动约束方式下的反射面RRMS只有全约束方式下的1/80。所以在天线的设计建造过程中应尽量减小座架对俯仰轴的热应力约束，即允许反射体沿轴承轴向有一定游隙。目前实际天线中已自然存在此效应，这也是夏天与冬天夜晚温度虽相差很大，但天线仍可以正常工作的原因。

由上述分析可知，俯仰轴的对称浮动约束方式更接近实际情况，此约束方式下的天线在均匀温升时，反射面精度几乎不受影响。

2.2.2日照温度场及变形场分析

2011年1月1日—11月1日，QTT台址晴天约占40.3%，为所有天气中占比最高；台址风速≤3.0m/s的占比约为50.3%，风速≤4.0 m/s的占比为69.5%[7]。选取典型工况对天线结构进行全天的温度场、变形场时程分析：2011年6月21日，晴天无云，风速分别取0、1.5、3.0 m/s，采用俯仰轴对称浮动约束；天线俯仰角姿态如图10所示，分别取0°、45°和90°，方位角取90°。

图10 各俯仰角姿态Fig. 10 Attitude of each pitch angle

天线仰天姿态下温度场分布特征及机理探讨如下：图11至图14分别为06：00、08：00、14：00、18：00的天线桁架与主反射面的温度场分布情况（俯仰角90°、风速3.0m/s）。由图11至图14可知，日照热在天线结构上产生的温度场呈现出不均匀、时变的特点：日出时太阳位于天线的东偏北方向，日落时太阳位于天线的西偏北方向；06：00～18：00，太阳在反射面上的直射点由西南向东南逆时针转动。在06：00和18：00，由于太阳侧照天线及反射面遮挡，天线温度场的不均匀性更为明显，桁架最高温度分别为24.0 ℃、37.1℃，均高于同时刻反射面的温度最大值（23.2 ℃、35.7℃）。08：00、14：00的反射面可分为3个比较均匀的温度区域，即高温区、过渡区和低温区；这两个时刻反射面最高温度分别为32.4℃、45.2℃，均高于同时刻桁架的最大温度值27.9℃、40.3℃，反射面与桁架的温差分别为4.5℃、4.9℃。

图11 06：00天线温度场分布Fig. 11 Distribution of antenna temperature field at 06：00

图12 08：00天线温度场分布Fig. 12 Distribution of antenna temperature field at 08：00

图13 14：00天线温度场分布Fig. 13 Distribution of antenna temperature field at 14：00

图14 18：00天线温度场分布Fig. 14 Distribution of antenna temperature field at 18：00

天线仰天姿态下温差特点及机理探讨如下：图15、16给出了背架的温度/温差变化曲线。可见，不同风速下背架的温度/温差变化规律相同：从05：00到19：00，太阳高度角大于0°，背架受到太阳辐射，故其最高温度高于空气温度。不同风速下背架的温度最大值均出现在14：00左右，最高温度可达39.7℃（3.0m/s）、41.7℃（1.5m/s）。从19：00到次日凌晨05：00，背架温差较小，且其温度低于空气温度，这是由于夜晚无太阳直接照射，而远空温度低于空气温度，背架向远空散热。背架的温差双峰分别出现在06：00、18：00，风速3.0m/s时分别为4.8℃、5.1℃，风速1.5m/s时分别为6.5℃、6.9℃。在05：00、19：00时，太阳侧照背架，所以不同风速下背架的温差会急剧变化；从05：00到19：00，不同风速下背架温差分别在4～5℃（3.0m/s）、5～7℃（1.5m/s）范围内变化。

图15 背架温度Fig. 15 Temperature of backup structure

图16 背架温差Fig. 16 Temperature difference of backup structure

风速对温度场的影响及机理探讨如下：由图15、16中不同风速曲线对比可见，随着风速增大，背架的温升、温差减小。这是由于在太阳辐射强度一定的条件下，风速越大，对流换热越强烈，背架高温区相对于空气的温升越小，而不同风速下背架的最低温度本就与空气温度相仿，故总体温差减小。由分析依据可知，上述结论在天线任意俯仰角下都成立。

仰天姿态下天线热变形规律分析如下：图17～20分别为对应于图11～14温度场的反射面热变形形状。在非均匀温度场的作用下，反射面变形可分为两部分：一部分是由结构整体温升产生的变形，另一部分是由结构温度梯度产生的变形。前者占主导，但其对反射面精度几乎没有影响；后者导致反射面出现以扭转、畸形为主的不均匀变形。

图17 06：00主反射面变形云图Fig. 17 Deformation cloud of the main reflector at 06：00

图18 08：00主反射面变形云图Fig. 18 Deformation cloud of the main reflector at 08：00

图19 14：00主反射面变形云图Fig. 19 Deformation cloud of the main reflector at 14：00

图20 18：00主反射面变形云图Fig. 20 Deformation cloud of the main reflector at 18：00

反射面不同时刻的热变形具有不同特点，整体沿俯仰轴呈现出“元宝”状的变形趋势，反射面边缘变形与扭转角度较大，且温差越大其不均匀变形程度越大（对比图19及20）。

2.2.3日照热误差及其补偿分析

如图4所示，定义RRMSd为变形后反射面相对于原设计抛物面的均方根值，其反映了反射面偏离原设计面的程度，即绝对变形情况。变形后反射面点B到原来位置（A）的偏差记为 ∆1，其计算如下：

则RRMSd计算公式为：

图21为俯仰角90°、不同风速下RRMSd曲线。仰天姿态下反射面绝对变形随时间、风速变化规律及其机理分析如下：05：00至19：00，绝对变形量随太阳高度角先增大后减小，最大值出现在12：00左右，原因是此时太阳辐射强度最大；在05：00和19：00，RRMSd会发生“跳跃”式突变，因此时天线经历了太阳入射角度的状态变化。风速方面：RRMSd随风速增大而减小，因为风速越大，天线与环境间的换热速率越快，结构整体的温升越小，由热引起的绝对变形量就越小。

图21 俯仰角90°，不同风速下R RMSdFig. 21 RRMSd at different wind speeds at pitch Angle of 90°

仰天姿态下反射面宏观变形时空特征分析：由最佳吻合抛物面的定义可知其参数变化反映了反射面宏观的位姿变化情况，故由式（21）即可得到副面的位置调整量UX、UY、UZ，亦可用于评价反射面的宏观刚体移动。

图22为俯仰角90°下副反射面补偿量。

图22 不同风速下副反射面位置调整量Fig. 22 Position adjustment of subreflector at different wind speeds

由图22可得出以下结论：

1）副反射面沿X轴（俯仰轴）的位置调整量关于12：00呈“奇函数”对称。其原因为：04：00到07：00，太阳在天线的东偏北方向升起（见图11），侧照X轴正半轴一侧的天线结构，被照射侧结构温度高于另一侧，使反射体受热膨胀沿X轴正向产生偏移，此时天线结构的温差最大，故补偿量达到最大；从07：00至10：00，随着太阳高度角增大，太阳照射区逐渐“翻过”反射面，出现在X轴负半轴一侧的天线结构上（见图2），其变形趋势与04：00至7：00的相反，故副面沿X向的补偿量开始减小；10：00至12：00，太阳照射位置靠近天线中心，导致结构温度的不均匀性减小，所以补偿量幅值较小；中午12：00之后，太阳与天线相对位置的变化情况与上午时段刚好相反，故补偿量的变化规律也与上午时刻相反。

2）05：00至14：00，随天线结构温升的增大，其热膨胀变形量一直增大，故副反射面沿Z轴（焦轴）的补偿量单调增大；14：00至20：00，随天线结构温升的减小，补偿量单调减小。

3）风速的影响：对比图22中不同风速的结果可知，补偿量的变化趋势相同但幅值不同。其原因为：风速增大，天线与环境间的换热速率越快，导致背架温差减小，使得其热变形减小，补偿量自然就减少了。

结合温度场、变形场、绝对变形量以及补偿量数据，QTT日照热变形的规律及机理总结如下：

1）太阳侧面照射天线时，即照射反射面凸面，不均匀变形最剧烈；从侧照到部分照射反射面凹面或者相反的过程（见图2），结构不均匀变形会经历峰值；其发生具体时间取决于天线–太阳所呈现出的相对姿态。

2）温度分布的非均匀性是结构不均匀变形的主要原因，温差越大不均匀变形程度越大。

3）风速不改变温度场的空间分布，只影响其幅值；不同风速下同一姿态反射面的热变形具有相同的变化趋势，热变形随风速增大而减小。

在QTT热误差控制及其补偿中，可利用上述规律，通过连续测量天线的温度场、风速及太阳辐射强度等数据，对未来反射面的热变形变化趋势进行预测，形成补偿方案。

2.2.4补偿后反射面精度分析

图23～25为采用副反射面位置补偿后，不同工况下反射面RRMS的变化曲线。由图23～25可知，RRMS曲线皆形如“马鞍”，且数值远小于RRMSd曲线。RRMS在05：00、19：00会发生急剧变化，因此刻太阳侧照天线；在13：00，RRMS出现谷底。不同风速下RRMS时间变化趋势相同，但数值随风增强而减小，这是由于大风加快了天线与环境间的对流换热，背架整体的温差变小，反射面不均匀热变形随之减小，使得反射面形状精度变高。

俯仰角为45°、90°时的RRMS具有与类似变化趋势，但波峰/波谷出现的时间及持续时间不同，这是由于不同姿态下太阳–天线相对位置不同，从而导致太阳“翻过”反射面的时间不同。

图26、27为同一时刻不同风速下主反射面光程差的分布云图，对比可知：风速对反射面光程差的空间分布趋势没有影响；但随风速增大，背架温差变小，反射面的光程差减小。图27、28为不同时刻同一风速下光程差的分布云图，可见其分布绕天线方位轴转动，与太阳入射角有相关性。对比左右图可知，采用副反射面位置补偿可有效减小天线的热误差。