亚毫米波天线有限元模型参数修正复合优化方法∗

张明珠 王海仁 左营喜

(1 中国科学院紫金山天文台 南京 210034)

(2 中国科学技术大学天文与空间科学学院 合肥 230026)

(3 中国科学院射电天文重点实验室 南京 210034)

1 引言

随着射电望远镜不断向更高工作频率、更大口径发展,对天线的指向精度和等效反射面精度(含各反射面及其准直精度)要求越来越高,目前亚毫米波望远镜指向精度要求已趋近1′′,等效反射面精度要求已趋近10µm.在重力、温度、风等载荷作用下,天线的支撑结构、反射面面形以及主、副反射面间的对准关系会发生变化,从而显著影响天线的指向精度和等效反射面精度,因而掌握天线随工况的变形规律并实时修正变形影响十分重要.实时测量天线变形是实现实时修正的最好手段,但限于目前的测量方法和条件,高精度实时测量还难以实现; 天线的变形往往只能在有限工况下测量(例如热致变形),甚至难以测量(例如风致变形).结构有限元分析是目前预测天线变形规律普遍采用的方法.然而由于复合材料(如碳纤维复合材料)和复合结构(如三明治夹心结构)的应用越来越广泛,其机械和物理特性参数随制作工艺或批次变化较大,难以准确预知,而且天线有限元模型在结构上有诸多简化,种种因素导致天线变形的有限元分析结果与实测结果之间存在偏差[1].因此,初期建立的天线有限元模型必须通过优化修正[2−3],使之与实际天线高度相符,才能实际用于天线变形高精度预测.依据修正对象分类,有限元模型修正主要有两类: 矩阵型修正和参数型修正[4].矩阵型修正方法主要针对模型刚度矩阵和质量矩阵进行修正,使模型矩阵特征值和特征向量逼近试验结果[5−6].参数型修正是通过优化调整结构几何参数、物理材料参数等使有限元模型分析结果与实际结构工作性能趋于一致[7−8].关于射电望远镜天线有限元模型修正方法少有研究,文献[9]研究了采用矩阵型修正方法的大型射电望远镜天线的有限元模型修正.本文基于零阶优化法和一阶优化法,提出了亚毫米波天线有限元模型参数修正复合优化方法.为验证这种优化方法的可行性,本文以1.2 m亚毫米波天线有限元模型修正为例,通过假想实验分别对多个材料参数在不同的载荷工况下进行复合优化.将初始建立的有限元模型与假想实际天线作对比,初始建立的有限元模型的物理材料参数值为A,假想实际天线的物理材料参数为A′(未知量),对初始模型进行参数复合优化使其分析结果逼近假想实际天线结构的结果,并且初始模型物理材料参数A能够收敛到假想实际天线结构物理材料参数值A′.这里A和A′既可以分别指单一材料参数值,也可以分别指一组材料参数值.最后对优化修正后的模型进行分析,得到重力载荷下天线主反射面变形与俯仰角变化的关系以及主副面相对位置随俯仰角变化的规律.

2 有限元模型参数修正复合优化方法原理

本文提出的优化方法包括3步: 第1步,建立亚毫米波天线参数化有限元模型,并进行参数灵敏度分析; 第2步,以第1步中得到的高灵敏度参数为优化变量建立优化函数,以初始模型与假想实际天线变形特性的差值最小为目标进行优化; 第3步,在前两步的基础上,基于零阶优化法和一阶优化法对亚毫米波天线有限元模型进行复合优化,使天线有限元模型逐渐逼近实际的天线结构.

2.1 天线有限元模型

本文对一个1.2 m亚毫米波望远镜天线进行优化方法可行性试验.望远镜为卡塞格林式,是斜轴式望远镜,其有限元模型如图1所示,该模型包括机架、背架和反射面.实际工作的望远镜在背架上一角放置有测量仪器,有限元模型中用质量单元模拟该仪器(即图1中标注的位置M).机架整体采用钢结构,背架和副镜支撑架为碳纤维复合材料(CFRP),主反射面(以下简称主面)采用CFRP-铝蜂窝-CFRP三明治复合结构.CFRP一般是各向异性材料,该天线所使用的碳纤维复合材料面板采用准各向同性铺层,在平面上两个方向近似各向同性,在厚度方向上的弹性模量为平面方向上弹性模量的1/10.有限元模型中CFRP和钢的材料属性如表1.

图1 口径1.2 m亚毫米波天线有限元模型Fig.1 The finite element model of 1.2 m submillimeter-wave antenna

表1 CFRP和钢的材料属性Table 1 The material properties of CFRP and steel

在对模型进行优化实验之前,已对模型进行过重力变形、热变形以及模态分析,分析结果表明有限元模型的正确性(这里正确性是指天线特性随载荷变化的规律合理,即分析结果在合理范围之内,与实际模型相比不一定具有准确性),可以在此初始模型上进行优化实验.本文通过对该天线的有限元模型分别施加重力载荷和温度载荷进行变形分析,对天线主面的变形特性进行有限元模型参数复合优化仿真.

2.2 参数灵敏度分析

灵敏度分析的原理是将建立模型的某些输入参数描述成服从某种概率分布的不确定性变量,经过大量的采样点计算统计分析出这些输入参数对输出参数的影响程度以及输入参数和输出参数之间的关系.由于天线有限元模型中所使用材料的属性较多,不同属性的变化对天线变形的影响大小不一,灵敏度分析有助于确定影响大的材料属性,从而使得研究更有效.

灵敏度分析涉及随机输入变量和随机输出参数,其中随机输入变量是指直接影响分析结果的输入数据,随机输出参数是指有限元分析的结果.本文使用的天线有限元模型主要选用材料为钢和CFRP,因此选取这两种材料的属性参数为随机输入变量,将天线主面变形误差作为随机输出参数来进行灵敏度分析.重力和温度载荷下对应的灵敏度分析变量参数设置如表2.

表2 重力和温度载荷下天线主面变形误差灵敏度分析设置Table 2 The sensitivity analysis and settings of the main-reflector deformation errors under gravity and temperature loads

2.3 优化目标设计

结构的优化设计一般是通过建立优化模型,在约束条件下运用优化算法迭代运算得到目标函数的极小值.优化过程涉及设计变量、状态变量和目标函数,一般优化问题表示为[8,10−11]:

设计变量x: 通常是结构的几何参数、材料属性等,是优化的对象,在整个优化过程中不断随着目标函数变化;分别为设计变量的上下限,用于约束设计变量的取值范围;n为设计变量的个数.

状态变量sj1(x)、tj2(x)、wj3(x): 是设计变量x的函数,也是优化设计的约束条件,其中为状态变量的上限,为状态变量的下限;m1、m2、m3分别为状态变量sj1(x)、tj2(x)、wj3(x)的个数.

目标函数f: 是设计变量x的函数,一般求取目标函数的最小值来作为优化设计的目标.

实际天线的变形是通过两种状态下的测量结果相减获得,为了与实际一致,对两种载荷状态下天线有限元模型(可以是同一俯仰姿态施加不同载荷,也可以是在不同俯仰姿态施加相同载荷)进行两次分析,将得到的两个分析结果相减后作为实际的变形结果.在材料属性参数为A的初始有限元模型中,天线主面有N个节点,两种载荷状态下天线主面某一节点i相对于理想抛物面的位移量分别为那么该节点相对变形位移量为

同样在两种载荷状态下,分别对材料属性参数为A′的假想实际天线进行分析,得到的天线主面节点相对于理想抛物面的位移量分别为主面节点相对变形位移量为

优化目标为通过将材料属性参数A设为设计变量,来使得有限元模型与假想实际天线的主面节点变形位移量的差距最小,即令

最小,因此目标函数为:

由于

故选取状态变量为:

因此针对天线主面变形的优化问题表示为:

2.4 零阶优化法和一阶优化法

零阶优化法和一阶优化法都是将约束的优化问题转化为无约束问题来求解,区别在于零阶法不利用一阶信息,能够得到全局最优解但精度低; 而一阶法运用因变量的偏导数,精度高但是收敛速度慢,有可能得到局部最优解.这两种方法各有优缺点,因此将两种优化算法相结合进行复合优化,先用零阶优化法得到精度较低的全局最优解,基于这个解再用一阶优化法提高精度,最后得到精度高的全局最优解,即进行两次复合优化[10−11].

对于(1)式表示的优化问题,零阶法采用最小二乘拟合以及加罚函数的方法进行优化.首先拟合设计变量、状态变量和目标函数的响应函数,拟合后的目标函数为:

其中,a0、ak和bkj均为常系数.其次对设计变量和状态变量加罚函数,将有约束的问题转化为无约束的问题:

其中f0为目标函数参考值,p为响应面参数,X(xk)、S(sj1)、T(tj2)、W(wj3)分别为设计变量x和状态变量sj1(x)、tj2(x)、wj3(x)的罚函数,当设计变量或者状态变量的值接近其上限或者下限时,相应的罚函数的值会急剧增加.

一阶优化同样对设计变量和状态变量加罚函数,由目标函数和罚函数组成无约束的目标函数:

其中f0为目标函数参考值,q为控制约束的参数,Px(xk)以及Ps(sj1)、Pt(tj2)、Pw(wj3)分别为设计变量x和状态变量sj1(x)、tj2(x)、wj3(x)的罚函数.

一阶优化由设计变量的一阶偏导数来确定优化搜索方向.对于第l次优化迭代,其优化搜索方向d(l),设计变量值为x(l),那么下一次迭代的变量值为:

其中x(l+1)为第(l+1)次迭代的设计变量值,hl为线性搜索参数,对应于搜索方向d(l)上的最小步进值.根据(11)式搜索方向可以分为两部分:

优化的迭代循环在达到收敛精度时终止,收敛条件为:

其中τ为目标函数的公差,f(l)为当前迭代步的目标函数值,f(l−1)为前一迭代步的目标函数值,f(b)为最佳设置处的目标函数值.

对于(8)式表示的优化问题,通过零阶优化法和一阶优化法优化得到最优解即最小的f时,所对应的设计变量A若能收敛到假想实际天线的材料属性A′,则表明本文提出的优化方法可行.

3 天线主面系统有限元模型优化仿真

3.1 灵敏度分析

根据2.2节的叙述,对1.2 m亚毫米波天线有限元模型分别施加重力载荷和温度载荷分析输入参数对天线主面变形误差的灵敏度,按照表2中对这两种载荷下灵敏度分析的设置,分析结果分别如图2和图3.其中柱状图左边的刻度表示随机输入变量与天线主面变形误差的相关系数,其最大范围是从−1到1.从图2可以看出,重力载荷下对主面变形误差RMSG影响较大的参数是EX-CARB和DENS-CARB,即碳纤维复合材料(CFRP)的弹性模量和密度,相关系数分别为−0.751和0.586,材料的泊松比对主面变形误差没有影响; 由图3看出,温度载荷下对主面变形误差RMST影响较大的参数为ALPX-CARB、ALPX-ST,分别表示CFRP的热膨胀系数和钢的热膨胀系数,相关系数分别为−0.695和0.595,材料的密度对主面变形误差没有影响.因此接下来分别选取CFRP的弹性模量和热膨胀系数作为设计变量来进行优化仿真.

图2 重力载荷下天线主面变形误差灵敏度分析结果Fig.2 The results of sensitivity analysis of the antenna main-reflector deformation error under gravity load

图3 温度载荷下天线主面变形误差灵敏度分析结果Fig.3 The results of sensitivity analysis of the antenna main-reflector deformation error under temperature load

3.2 重力载荷下主面系统优化仿真

由于工作中的真实天线总会受到重力等载荷的影响,因此不存在无重力影响的理想情况,也就无法测量得到变形后主面节点相对理想抛物面的位移量,实际测量时只能通过两个俯仰角的面形相减来得到重力引起的主面变形误差.根据2.3节的叙述,选取的两个载荷状态分别是在天线俯仰角为0◦和90◦时对有限元模型施加重力载荷,分析得到的天线位移变形云图如图4,背架上放置的仪器位于图4中标注的位置M.

图4 重力载荷下天线变形:(左)天线仰角为0◦;(右)天线仰角为90◦.图中底部的颜色条标识天线的变形位移,单位为mmFig.4 The deformation of the antenna under gravity load:(Left)the elevation angle is 0◦;(Right)the elevation angle is 90◦.The color bar on the bottom of the panels identifies the deformation displacement of the antenna,and the unit is millimeter

根据3.1节重力载荷下的参数灵敏度分析结果,选取EX-CARB即CFRP的弹性模量为初始有限元模型的设计变量,初始值为E0=68000 N·mm−2,假想实际天线中CFRP的弹性模量为E′=50000 N·mm−2(对于低弹性模量类型的CFRP,其弹性模量低于100000 N·mm−2).由(2)–(8)式,记设计变量为E,状态变量为SG,施加重力载荷后主面变形的优化设计为:

零阶优化法和一阶优化法的迭代次数均选为15,先用零阶法对模型进行优化得到全局最优解,在此基础上再用一阶法进行优化提高精度,目标函数fG与设计变量E的迭代过程如图5所示,在第7次迭代取得最优解.经过零阶法和一阶法两次优化达到最优解的迭代次数、目标函数值以及设计变量值如表3.当目标函数达到最小值时,设计变量为50134 N·mm−2,与E′的相对误差为134 N·mm−2,绝对误差为0.27%.

表3 重力载荷下主面变形优化设计结果Table 3 The results of the optimization of the main-reflector deformation under gravity load

图5 目标函数fG随E的变化以及二者的迭代过程Fig.5 The dependence of the objective function fG upon E and the iteration processes

重力载荷下,当天线仰角变化90◦时,假想实际天线结构即CFRP的弹性模量为E′=50000 N·mm−2时分析得到的天线主面变形分布如图6,天线主面变形误差RMSG为7.4529µm; 初始模型(对应于E0=68000 N·mm−2)的天线主面重力变形分布如图7左图,RMSG为6.9771µm,与假想实际天线结构主面变形的偏差分布如图7右图,偏差errord为0.7398µm,优化结果模型(对应于E=50134 N·mm−2)的主面变形分布如图8左图,RMSG为7.4486µm,与假想实际天线结构主面变形的偏差分布如图8右图,errord极小,为0.0492µm.由表3以及图6、7、8可以看出,重力载荷下优化后天线特性(图8左图)与假想实际天线结构的特性(图6)极为相符,并且优化得到的参数E收敛到假想实际天线结构的材料参数E′,初步表明了所提出优化方法的可行性.

图6 重力载荷下假想实际天线结构的主面变形分布(对应于E′=50000 N·mm−2)Fig.6 The main-reflector deformation distribution of the imaginary real structure of the antenna under gravity load(corresponding to E′=50000 N·mm−2)

图7 左: 重力载荷下初始模型的主面变形分布情况(对应于E0=68000 N·mm−2); 右: 初始模型与假想实际天线结构的形变偏差Fig.7 Left: the main-reflector deformation distribution of the initial model under gravity load(corresponding to E0=68000 N·mm−2); Right: the deformation deviation between the initial model and the imaginary real structure of antenna

图8 左: 重力载荷下优化模型的主面变形分布情况(对应于E=50134 N·mm−2); 右: 优化模型与假想实际天线结构的形变偏差Fig.8 Left: the main-reflector deformation distribution of the optimization model under gravity load(corresponding to E=50134 N·mm−2); Right: the deformation deviation between the optimization model and the imaginary real structure of antenna

3.3 温度载荷下天线主面变形优化仿真

选取天线俯仰角为指向水平方向,参考南极Dome A存在的极端温度条件以及望远镜一面朝阳一面背阴情况,对有限元模型分别施加两种温度载荷进行分析.第1种是平均温度为5◦C的温度载荷,载荷分布如图9左图; 第2种是在平均温度为−55◦C的温度载荷下再加上沿天线指向的1◦C·m−1温度梯度载荷,温度载荷分布如图9右图所示.

图9 天线有限元模型温度载荷分布:(左)5◦C温度载荷;(右) −55◦C温度载荷加沿天线指向的1◦C·m−1温度梯度载荷.TMIN表示最低温度载荷,TMAX表示最高温度载荷,右图的颜色条标识温度,单位为◦CFig.9 The temperature load distribution of the finite element model of the antenna:(Left)5◦C mean temperature;(Right)superimposing 1◦C·m−1 temperature gradient on the −55◦C mean temperature.TMIN is the minimum temperature,and TMAX is the maximum.The color bar in the right panel identifies the temperature.The unit is ◦C

根据文献[12]中对CFRP热膨胀系数的叙述,设假想实际天线所使用的CFRP热膨胀系数为α′=1.5×10−6K−1,将初始有限元模型的CFRP热膨胀系数设为设计变量,其初始值为α0=3.0×10−6K−1.由(2)–(8)式,记设计变量为α,状态变量为ST,对望远镜施加温度载荷后主面变形的优化设计表示为:

零阶优化法和一阶优化法的迭代次数均选为15,整个优化迭代过程如图10所示,在第9次迭代取得最优解.依次经过零阶法和一阶法两个阶段的优化后,得到最优解时的迭代次数、相应的目标函数值以及设计变量值在表4中列出,达到最优解时对应的设计变量为1.4935×10−6K−1,与α′的相对误差为0.65×10−8K−1,绝对误差为0.43%,因此α收敛到α′.

表4 温度载荷下主面变形优化设计结果Table 4 The results of the optimization of the main-reflector deformation under temperature load

图10 目标函数fT随α变化以及二者的迭代过程Fig.10 The dependence of objective function fT upon α and the iteration processes

在温度载荷从均匀温度5◦C变化到极端的不均匀温度即在−55◦C的温度载荷上叠加1◦C·m−1的温度梯度载荷,假想实际天线结构即CFRP热膨胀系数为α′=1.5×10−6K−1时分析得到的天线主面热变形分布如图11,RMST为14.1649µm;初始模型(对应于α0=3.0×10−6K−1)的RMST为4.1598µm,变形误差分布如图12左图,与假想实际天线结构的偏差分布如图12右图,errord为10.4215µm; 优化结果模型(对应于α=1.4935×10−6K−1)的主面变形分布如图13左图,RMST为14.2085µm,与假想实际天线结构的偏差分布如图13右图,errord为0.0594µm,说明优化结果模型极其接近假想实际模型.表4以及图11、12、13再次证明了这种优化方法的可行性.

图11 温度载荷下假想实际天线结构的主面变形分布(对应于α′=1.5× 10−6 K−1)Fig.11 The main-reflector deformation distribution of the imaginary real structure of the antenna under temperature load(corresponding to α′=1.5× 10−6 K−1)

图12 左: 温度载荷下初始模型的主面变形分布情况(对应于α0=3.0× 10−6 K−1); 右: 初始模型与假想实际天线结构的形变偏差Fig.12 Left: the main-reflector deformation distribution of the initial model under temperature load(corresponding to α0=3.0× 10−6 K−1); Right: the deformation deviation between the initial model and the imaginary real structure of antenna

图13 左: 温度载荷下优化模型的主面变形分布情况(对应于α=1.4935× 10−6 K−1); 右: 优化模型与假想实际天线结构的形变偏差Fig.13 Left: the main-reflector deformation distribution of the optimization model under temperature load(corresponding to α=1.4935× 10−6 K−1); Right: the deformation deviation between the optimization model and the imaginary real structure of antenna

4 重力载荷下天线变形规律分析

天线在不同俯仰角受重力影响的程度大小不一,主面面型的变形和主副面的相对位置都会随俯仰角的变化而变化.下面将修正优化后的天线有限元模型作为实际天线,分析重力载荷下,随俯仰角逐渐从水平指向变化到天顶指向,天线主面面型的变化规律以及主副面之间的相对位置关系.

天线主面面型的变化通过主面变形误差均方根(Root Mean Square,RMS)来表示,重力作用下变形误差RMS与俯仰角EL的关系如图14所示.一般天线主面变形误差随俯仰角的变化关系曲线是正弦曲线或余弦曲线,但是由于试验对象1.2 m天线采用斜轴式结构,不同于常用的方位俯仰系统,随着天线俯仰角度的变化,整个反射面系统都会绕自身轴线旋转,副面的3根非圆柱型支撑杆以及背架上仪器设备的质量均会使得沿重力方向主副面的支撑状态不断变化,可能使得天线主面变形误差与俯仰角的关系曲线偏离正弦或余弦曲线.

图14 重力载荷下天线主面变形误差RMS随俯仰角变化Fig.14 The dependence of the main-reflector deformation errors RMS upon elevation angle under gravity load

主副面的相对位置关系包含不同工况下副面的位置相对于理想抛物面在空间直角坐标系X、Y、Z方向上的位移以及变形后副面轴线相对于理想抛物面倾角的X、Y分量.其中Z方向即为理想天线的指向方向.

图15中的ux、uy、uz分别表示相对于理想抛物面,重力作用下副面位置在X、Y、Z方向上的位移.随着天线俯仰角从0◦变化到90◦,副面位置在X方向上的位移ux逐渐减小,且都是偏向X轴正方向; 在Y方向上的位移uy则是一直偏向Y轴的负方向,位移量先增大后减小,在35◦时位移量最大; 副面位置在Z方向上的位移uz先偏向正方向逐渐减小,后偏向负方向逐渐增大,在俯仰角为55◦时位移量最小; 图16给出了天线变形前后副面位置偏移的绝对距离随俯仰角的变化,对其进行拟合发现绝对距离与俯仰角度的关系近似是三角函数.一般采用方位俯仰结构的望远镜副面在X方向上的位移不随俯仰角度的变化而变化,但是对于斜轴式结构,在不同俯仰角度,副面在X方向上的位移会不相同.

副面轴线倾角在X、Y方向的分量随俯仰角度变化的关系曲线如图17.在重力载荷下,副面轴线倾角在X方向上的分量随俯仰角的增大而减小;Y方向上的分量在俯仰角约为0◦–75◦之间先增大后减小,在70◦–90◦之间又逐渐增大.同样,对于采用方位俯仰结构的望远镜,变形后副面轴线倾角在X方向上的分量在各个俯仰角度是相同的; 但对于斜轴式结构,副面轴线倾角的X分量会随俯仰角度变化.

图15 重力作用下副面位置在X、Y、Z方向上的位移随俯仰角变化Fig.15 The dependence of the displacements in X, Y,and Z directions of sub-reflector upon elevation angle under gravity load

图16 重力作用下副面位置偏移的绝对距离随俯仰角变化Fig.16 The dependence of the absolute distance value of the sub-reflector position deviation upon elevation angle under gravity load

图17 重力作用下副面轴线倾角的X、Y 方向分量随俯仰角变化关系Fig.17 The dependence of the inclination angles of sub-reflector axis relative to X-direction and Y-direction upon elevation angle under gravity load

5 总结

针对亚毫米波天线有限元模型与实际天线变形特性的差异,本文提出了有限元模型参数修正复合优化方法,并以1.2 m亚毫米波天线在重力载荷和温度载荷下的变形特性为例进行复合优化.首先对1.2 m亚毫米波天线的有限元模型进行参数灵敏度分析,选择对天线主面变形灵敏度高的材料属性参数作为优化变量,将初始模型与假想实际天线主面变形的差作为优化目标,采用零阶优化法与一阶优化法相结合的方法对天线的有限元模型进行优化.在重力载荷情况下,将初始有限元模型中CFRP的弹性模量作为设计变量,初始值为68000 N·mm−2,当达到最优目标即差值最小时对应的设计变量为50134 N·mm−2,与假想实际天线的CFRP的弹性模量(50000 N·mm−2)的绝对误差为0.27%; 在温度载荷下,将初始有限元模型中CFRP的热膨胀系数作为设计变量,初始值为3.0×10−6K−1,当达到最优目标时所对应的设计变量为1.4935×10−6K−1,与假想实际天线的CFRP的热膨胀系数(1.5×10−6K−1)的绝对误差为0.43%; 以上两种情况都表明采用的优化方法可行.将经过修正优化的模型作为实际天线进行重力变形分析,得到天线主面变形误差以及主副面相对位置随俯仰角变化关系,未来可根据实验数据来检验这种关系的正确性.