2. 西北工业大学 无人机特种技术重点实验室, 陕西 西安 710065;
3. 军械工程学院 无人机工程系, 河北 石家庄 050003
太阳能无人机可用于执行大范围侦查、监测和通信中继等任务，是近年来高空长航时无人机主要的研究对象之一。受能源能量密度低的瓶颈限制，为满足长航时续航任务需求，这类飞机一般采用大展弦比机翼设计，并大量使用超轻复合材料以尽可能减轻结构重量；在气动载荷的作用下，翼尖最大弯曲变形可达半翼展的25%[1, 2, 3]。随变形的增加，无人机的结构特性、气动特性以及质量特性均发生较大改变，采用线弹性结构模型不能够精确地分析此类柔性飞机的配平与飞行载荷特性[4, 5, 6, 7]。尽可能高精度地模拟这类大柔性无人机的配平及飞行载荷特性，是这类飞机研究的关键技术之一。
大柔性飞机的配平及飞行载荷分析是其结构设计的主要依据，考虑线性静气动弹性效应的配平方法比较成熟，典型的代表有MSC.Flightloads模块等[8]；当结构弹性变形较大时，该方法难以准确地描述飞机的静气动弹性特性。目前，考虑几何非线性效应且适用于大柔性飞机配平及飞行载荷分析的方法主要有：①Patil及其合作者基于几何精确本征梁模型开发非线性气动弹性及飞行动力学分析程序包NATASHA(nonlinear aeroelastic trim and stability of HALE aircraft[4, 5])；②Cesnik及其合作者基于应变梁模型开发的非线性气动弹性与飞行动力学耦合分析工具包UM/NAST(the university of Michigan′s nonlinear aeroelastic simulation toolbox)[6, 7]。杨超等采用更新的Lagrange(UL)格式求解结构几何大变形，研究了某常规布局飞机缩比模型的配平与飞行载荷特性[9]。本征梁模型与应变梁模型具有较低的阶数，但自变量不直接且通用性差；位移梁模型求解中的TL(total Lagrangian)法和UL法一般采用增加载荷步的途径提高计算精度，从而显著地降低了求解效率。CR位移有限元法把结构的运动分解为刚体运动和纯粹的弹性变形，不像传统的TL法和UL法那样从非线性位移-应变关系出发推导切线刚度矩阵，且具有更高的求解精度和计算效率[11, 12]，在大位移小应变几何非线性问题的求解中愈加受到重视。文献[3]采用CR静力学理论耦合CFD求解技术研究了大柔性太阳能无人机静气动弹性问题；文献[10]采用线性化的CR动力学理论研究了大柔性太阳能无人机非线性气动弹性稳定性等问题，并开发了适用于大柔性飞机非线性气动弹性求解的计算程序。
本文采用CR有限元技术耦合片条气动力模型提出了一种可以考虑几何非线性效应的大柔性飞机配平及飞行载荷分析方法，给出了配平迭代求解流程并编写了相应的Fortran计算程序。然后，以类“太阳神”布局太阳能无人机为例，以定高、恒速作为配平前提，研究了不同有效载荷下无人机的纵向配平及飞行载荷特性，最后讨论了所提出的非线性迭代配平方法的收敛性等。
1 几何非线性结构、气动、重力及分布式推力系统建模大尺度大柔性太阳能无人机结构重量轻，且一般采用分布式螺旋桨推力系统；一个完整的大柔性太阳能无人机配平及飞行载荷求解系统应包含结构、气动、重力及分布式推力系统等子模块。
1.1 几何非线性结构建模采用CR有限元理论计算结构大变形问题时，一般采用增量法求解。
以T表示单元坐标系下位移增量到总体坐标系下位移增量的位移增量转换矩阵：
单元坐标系下的广义节点力fl为:
式中,Kl是线性刚度矩阵,与线性梁单元的刚度矩阵一致,pl为单元坐标系下的节点位移。
总体坐标系下广义节点力向量f为:
KT就是所求的切线刚度矩阵,其中TTKlT为材料刚度矩阵,KTσ(fl)为几何刚度矩阵,是单元坐标系下单元节点力fl的函数且与结构变形相关。
那么增量平衡方程可以表示为:
1.2 片条气动力模型
当地有效迎角由其弦线与当地来流方向的夹角表示。
 |
| 图 1 剖面气动载荷示意图 |
根据单元坐标系下的速度投影表示入射来流速度的大小及有效迎角:
那么,气动载荷投影在单元坐标系下的表达式为:
当翼段布置有舵面并且偏转时,升力、阻力、力矩等可以表示为:
1.3 重力模型
重力描述在惯性坐标系下:
通过坐标转换将在惯性坐标系中定义的重力转换到变形后的无人机单元坐标系中,这时无人机的分布式重力模型为:
由于机翼剖面质心可能不在弹性轴上,导致重力产生力矩,将惯性坐标系下的惯性力投影到单元坐标系下,对各个轴产生的力矩:
式中,
。
螺旋桨推进系统一般安装在机翼前缘之前,如图2所示,yp为推力轴线到弹性轴的距离,ϕp为推力轴线与弦线的夹角(即安装角)。
 |
| 图 2 推进系统 |
把推力及其反扭矩转换到单元坐标系下,其表达式为:
2 大柔性飞机非线性配平及飞行载荷迭代求解流程
建立全机配平方程时,选择地轴系I作为惯性参考系,以原点始终为全机瞬时质心的平均轴系作为机体轴参考系。变形后无人机的刚体平衡方程为:
在考虑几何非线性静气动弹性效应的大柔性飞机配平分析中,机翼结构刚度、气动力的大小和方向、全机质心位置、螺旋桨推进系统作用点等参数均与变形有关,在求解的过程中需要随变形的增加而不断更新。因此,大柔性飞机非线性配平参数不能直接解出,需要迭代求解。构造非平衡力表达式如(16)式所示。
收敛判据函数Jtrim的表达为:
以定高、恒速平飞为例,假设配平参数分别为推力Tp、升降舵偏角δe和全机迎角α,令列向量S=(Tp,δe,α)T;采用Newton-Raphson迭代,构造雅克比切线矩阵Sk,T:
具体求解时,雅克比切线矩阵Sk,T可以根据有限差分法,在当前配平参数附近引入微量增量,通过数值法求解得到。
以Ftrim,k表示第k次配平迭代残差,那么,配平参数增量ΔS可以通过下式求解得到:
则有:
在配平迭代求解的过程中，若令任意迭代步下的弹性变形为零，则可以解得刚性飞机假设下的配平参数；若始终以初始位形状态下的结构刚度矩阵参与迭代求解，则可以得到线弹性结构假设下的配平参数。
综上所述，非线性配平迭代求解流程如图3所示：
 |
| 图 3 非线性配平迭代求解流程 |
以文献[5]等给出的类“太阳神”布局无人机作为算例，其几何尺寸及外形正视图、俯视图如图4所示，弦长为2.44 m，展长为72.8 m，根稍比为1；外侧机翼具有10°的上反角；舵面位于机翼后缘，贯穿整个机翼；5个电机布置在平直段内；3个吊舱A、B、C的位置如图4所示。
 |
| 图 4 无人机外形及几何尺寸 |
机翼剖面气动特性、结构参数以及载荷舱参数如表1、表2所示：
吊舱A、B的有效载荷为22.68 kg，仅改变中间吊舱C中的有效载荷(27.216~254 kg)，有效载荷质心在弹性轴下0.914 m，假设飞行速度为12.19 m/s，飞行高度为海平面高度，并假设舵面偏转角一致，各个推进系统的推力大小相同。
3.2 配平及飞行载荷分析以推进系统推力、舵面偏转角及全机迎角作为配平变量，3个配平参数随有效载荷增加的变化曲线如图5~图7所示。
 |
| 图 5 配平推力随有效载荷的变化 |
 |
| 图 6 配平舵面偏转角随有效载荷的变化 |
 |
| 图 7 配平迎角随有效载荷的变化 |
从图5中可以看出，随有效载荷的增加，无论是刚性飞机还是线性、非线性结构模型配平得到的单个电机的拉力基本一致并且保持不变，这主要是因为在所有的配平范围内，无人机的迎角比较小，并且采用的气动力模型中假设阻力系数与迎角无关。
搭载不同重量的有效载荷时，采用刚性、线性和非线性结构模型求解得到的配平状态下的舵面偏转角，随有效载荷的增加而变化的趋势如图6所示。
有效载荷较小时，采用刚性、线弹性模型及非线性模型假设得到的配平舵偏角基本一致，随着有效载荷的增加，刚性与线弹性模型假设解得的结果偏差逐渐增加。因此，变形较大时，刚性及线弹性结构模型不再适用于柔性无人机的配平求解。
无人机配平所需的迎角随不同的有效载荷的变化趋势如图7所示。
从图7可以看出，随着有效载荷的增加，采用线性结构模型配平所需的迎角始终小于刚性飞机的配平迎角；有效载荷较小时，采用非线性结构模型的配平迎角与线性理论基本一致，然而，当有效载荷大于180 kg时，采用非线性结构模型所需的配平迎角大于刚性飞机的配平迎角；这是因为有效载荷较大时，配平状态下的弯曲变形改变了气动载荷的作用方向，降低了全机的有效升力系数，而线弹性结构模型无法描述几何大变形引起的刚体旋转效应，所计算得到的弯曲位移在展向的分量偏小。
从上述推力、舵偏角、迎角等配平参数随有效载荷的增加而解得的结果可以看出，本文所提出的非线性配平方法与文献值均吻合得比较好，验证了本文所提出的大柔性飞机非线性配平方法的精度及所编写的程序的有效性。
变形及载荷是结构设计的主要评价标准之一，不合理的预估配平状态下的载荷及弹性变形可能会使大柔性飞机结构设计产生严重的偏差。随有效载荷的增加，翼尖位移及扭转角的变形趋势如图8~图10所示。
 |
| 图 8 翼尖弯曲位移随有效载荷的变化 |
 |
| 图 9 翼尖展向位移随有效载荷的变化 |
 |
| 图 10 扭转变形随有效载荷的变化 |
由图8中的计算结果可以看出，当翼尖弯曲变形小于展长的3.5%时，采用线弹性结构模型描述机翼变形进而对柔性飞机配平时仍然具有较好的精度，因此可以把3.5%展长的弯曲变形作为是否需要采用非线性结构求解器的衡量标准之一。
由图9中的计算结果可以看出，随着变形的增加，线性理论及非线性理论配平解得的翼尖展向位移均体现出减小有效展长的趋势;但其本质有所区别，其中线性理论的展向位移是由于机翼外段具有初始上反角引起的，而非线性理论还计及了几何大变形中的刚体旋转效应，因此变形较大时，非线性理论解得的展向位移远大于线性理论值。
由图10中的计算结果可以看出，线性理论与非线性理论解得的配平状态下的翼尖扭转角差距较小，这主要是因为算例中机翼扭转刚度较大，扭转变形小。
计及重力、气动力等因素的节点总载荷并投影在翼根坐标系下，在展向的分布如图11所示。
 |
| 图 11 配平状态下机翼总载荷分布 |
由图11中的计算结果可以看出，当弹性变形较大时，靠近翼尖区域的气动载荷在侧向的分量所占的比重也越大，造成的有效升力的减小是不可忽略的。
3.3 非线性配平迭代求解收敛性讨论由于采用了循环迭代法求解无人机的配平参数，有必要对循环迭代求解的收敛性进行讨论，以评价所提出的方法求解柔性飞机配平参数时收敛的速度、精度及稳定性。在每一轮循环的配平迭代中选择收敛残差Jtrim<10-8作为收敛标准，以满载荷状态为例，研究其配平迭代的收敛性。
从图12可以看出,算例非线性配平一般需要迭代12次以上，才能收敛。而静气动弹性模块循环求解时，一般需要5次迭代才能收敛；因此，一个完整的非线性配平中气动/结构迭代求解的次数将大于5×12次，若采用CFD技术(如有限体积法)求解气动力，计算5×12次气动力的花费将是非常大的。另外，采用非线性结构模型与线性结构模型配平时循环迭代次数基本一致，并没有因为几何非线性因素的引入而显著增加循环迭代次数。
 |
| 图 12 满载配平时配平参数收敛历程 |
太阳能无人机的展向尺度、结构柔度一般比较大，且允许较大的机翼弹性变形，配平时需要考虑大变形引起的几何非线性因素的影响，若仍采用刚性飞机或线弹性结构模型假设，得到的配平结果将严重偏离实际情况，基于此结论研究无人机的飞行品质等特性也是不可靠的。
本文基于CR有限元技术和片条气动力模型所提出的非线性配平与飞行载荷分析方法具有较好的精度及收敛性，满足大柔性太阳能无人机的配平及飞行载荷分析对求解精度的需求，并且具有较好的通用性。
| [1] | Kirk Flittie, Bob Curtin. Pathfinder Solar-Powered Aircraft Flight Performance[R]. AIAA-1998-4446 |
| Click to display the text | |
| [2] | 昌敏, 周洲, 郑志成. 太阳能飞机原理及总体参数敏度分析[J]. 西北工业大学学报, 2010, 28(5):792-796 Chang Min, Zhou Zhou, Zheng Zhicheng. Flight Principles of Solar-Powered Airplane and Sensitivity Analysis of Its Conceptual Parameters[J]. Journal of Northwestern Polytechnical University, 2010, 28(5):792-796 (in Chinese) |
| Cited By in Cnki (12) | Click to display the text | |
| [3] | 王伟, 周洲, 祝小平, 王睿. 考虑几何非线性效应的大柔性太阳能无人机静气动弹性分析[J]. 西北工业大学学报, 2014, 32(4):499-504 Wang Wei, Zhou Zhou, Zhu Xiaoping, Wang Rui. Static Aeroelastic Characteristics Analysis of a Very Flexible Solar Powered UAV with Geometrical Nonlinear Effect Considered[J]. Journal of Northwestern Polytechnical University, 2014, 32(4):499-504 (in Chinese) |
| Cited By in Cnki (1) | Click to display the text | |
| [4] | Patil M J, Hodges D H, Cesnik C E S. Nonlinear Aeroelasticity and Flight Dynamics of High-Altitude Long-Endurance Aircraft[R]. AIAA-1999-1470 |
| Click to display the text | |
| [5] | Patil M J, Hodges D H. Flight Dynamics of Highly Flexible Flying Wings[J]. Journal of Aircraft,2006,43(6):1790-1798 |
| Click to display the text | |
| [6] | Cesnik C E S, Su Weihua. Nonlinear Aeroelastic Modleing and Analysis of Fully Flexible Aircraft[R]. AIAA-2005-2169 |
| Click to display the text | |
| [7] | Su W, Cesnik C E S. Dynamic Response of Highly Flexible Flying Wings[J]. AIAA Journal, 2011, 19(2):324-339 |
| Click to display the text | |
| [8] | Rodden W P, Johnson E H. MSC/Nastran Aeroelastic Analysis User's Guide V68[M]. MSC Software Corporation, 1994 |
| [9] | Yang C, Wang L B, Xie C C, et al. Aeroelastic Trim and Flight Load Analysis of Flexible Aircraft with Large Deformations[J]. Sci China Tch Sci, 2012, 55:2700-2711 |
| Click to display the text | |
| [10] | 王伟, 周洲, 祝小平, 王睿. 几何大变形太阳能无人机非线性气动弹性稳定性研究[J]. 西北工业大学学报, 2015, 33(1):1-8 Wang Wei, Zhou Zhou, Zhu Xiaoping, Wang Rui. Exploring Aeroelastic Stability of Very Flexible Solar Powered UAV with Geometrically Large Deformation[J]. Journal of Northwestern Polytechnical University, 2015, 33(1):1-8 (in Chinese) |
| Cited By in Cnki (0) | |
| [11] | Crisfield MA, Galvanetto U, Jelenic' G. Dynamics of 3-D Co-Rotational Beams[J]. Computational Mechanics, 1997, 20:507-519 |
| Click to display the text | |
| [12] | Le T N, Battini J M, Hjiaj M. Dynamics of 3D Beam Elements in a Corotational Context: A Comparative Study of Established and New Formulation[J]. Finite Elements in Analysis and Design, 2012, 61:97-111 |
| Click to display the text |
2. Science and Technology on UAV Laboratory, Northwestern Polytechnical University, Xi'an 710065, China;
3. Department of UAV, Ordance Engineering College, Shijiazhuang 050003, china