近十几年来，可用于执行侦查、监测和通信中继等任务的高空长航时无人机引起了研究者广泛的研究兴趣。其中，太阳能无人机作为新一代高空超长航时无人机，其展弦比一般大于30，机翼结构面密度小于3 kg/m²，在气动载荷的作用下，翼尖最大弯曲变形可达半翼展的25%^{[1, 2, 3]}。随变形的增加，无人机的结构动力学特性及气动特性均发生较大改变，线性理论不能够精确地分析此类变形对结构、气动特性及气动弹性稳定性的影响^{[4, 5, 6, 7, 8]}。尽可能高精度模拟几何大变形对这类无人机气动弹性稳定性的影响，是目前研究的热点之一。

太阳能无人机在气动载荷的作用下，产生较大的弯曲变形，但机翼结构的局部变形仍然满足线性位移-应变关系，属于典型的大位移小应变几何非线性问题。目前，大展弦比大柔性机翼一般等效成几何非线性梁模型，借助有限元技术求解，采用的单元类型主要有基于位移的梁单元、基于应变的梁单元和基于速度和内力的本征梁单元^[8]。采用位移有限元法描述几何非线性问题按照参考构形的不同，一般分为3类：①TL(total lagrangian)法；②UL(updated lagrangian)法；③CR法。为了提高计算精度，使用TL法和UL法一般需要增加载荷步，从而显著地增加计算时间，引起了广大学者基于CR列式法对有限旋转问题的研究^{[9, 10, 11, 12]}。

几何大变形柔性机翼的气动弹性问题吸引了诸多学者的研究兴趣，Patil和Hodges采用本征梁模型，首次把几何非线性问题引入到大柔性机翼气动弹性问题的研究当中^{[4, 5, 6]}；谢长川^[13]采用线化UL有限元法也对这类问题开展了一定的研究。非线性气动弹性问题的求解一般比较繁琐，实际应用中需要舍弃一些非线性项而保留最主要的非线性项以简化计算；如何建立有效的简化计算方法，是这类非线性气动弹性问题研究的热点之一。由于机翼大振幅振动中回转惯性力矩的存在，考虑几何非线性效应的气动弹性运动方程一般要在时域内求解，但在非线性静平衡状态附近对气动弹性运动方程线性化，并引入准模态小扰动振动假设，在频域内研究大柔性机翼的气动弹性问题仍然是比较有效的工程分析方法之一；采用这种方法解得的临界发散速度称为非线性颤振速度，相应的临界发散频率称为非线性颤振频率。本文基于CR有限元理论，推导了大柔性机翼的切线刚度矩阵和质量矩阵，建立了考虑几何非线性效应的大柔性无人机结构动力学模型；引入准模态小扰动振动假设，对动力学模型进行线性化；采用建立在局部气流坐标系下的Therdorson片条非定常气动力，通过P-k法对无人机的非线性颤振速度和频率进行了求解。首先以文献算例为例对所发展的气动弹性稳定性分析方法及编写的Fortran代码进行了验证；然后，研究了几何大变形对类"太阳神"布局太阳能无人机的气动弹性稳定性的影响；最后对弯曲刚度、扭转刚度、弹性轴位置及剖面质心位置对几何大变形无人机的气动弹性稳定性的影响进行了深入的研究。

CR有限元法充分利用了几何非线性问题的小应变、大位移特征，弹性变形在单元坐标系内描述，位移-应变关系仍然是线性的。按照文献^[3]中介绍的单元坐标系、节点坐标系，这里直接引出单元坐标系下的节点位移表达式：

点击浏览公式

对(1)式进行微分得到单元坐标系下位移增量δp_l到总体坐标系下位移增量δp的增量转换矩阵T：

点击浏览公式

单元坐标系下的广义节点力f_l为：

点击浏览公式

总体坐标系下广义节点力向量f为：

点击浏览公式

点击浏览公式

点击浏览公式

K_T就是所求的切线刚度矩阵，其中T^TK_lT为材料刚度矩阵，K_Tσ(f_l)为几何刚度矩阵，是单元坐标系下单元节点力f_l的函数，更加详细的推导可参考文献^[3]等。

一般情况下，梁单元截面质心与刚心的距离较小，经典梁单元集中质量矩阵和一致质量矩阵均不考虑惯性力引起的弯扭耦合效应^{[11, 12]}。把大展弦比柔性机翼等效成非线性梁模型时，质心与刚心的距离一般不可忽略，需要对现有的质量矩阵进行适当的改进。本文在一致质量矩阵的基础上考虑弯扭耦合效应，重新构造了惯性力矩阵，得到切线质量矩阵。

大柔性机翼在变形的过程中可以认为其局部翼剖面没有发生翘曲^[12]，则剖面内所有点具有相同的角速度，截面线动量与角动量为：

点击浏览公式

截面惯性力和惯性力矩可以表示为:

点击浏览公式

准模态假设的主要思想是大柔性结构在外载荷的作用下发生几何非线性弹性变形,结构振动假设为在此弹性变形后的结构位形下的小扰动振动^{[4, 14]}。那么,在定常载荷作用下的静平衡状态附近,结构的小扰动振动方程为:

点击浏览公式

点击浏览公式

可以看出,切线质量矩阵M_T是结构位形的函数。

方程(9)的解可以假设为:

点击浏览公式

点击浏览公式

(11)式具有非零解的条件为:|K_T-ω²M_T|=0,称为系统的特征方程。求解方程(11)可得到准模态假设下系统的各个自由振动频率ω_i及系统的各个主振型X_i,各个主振型按质量矩阵标准化,并引入实模态矩阵A:

点击浏览公式

薄翼振动理论的真实绕流特点是机翼后缘不断产生自由涡,并顺气流向后流动形成尾涡区,与翼面附着涡一起诱导出翼面下洗。根据准模态假设及局部有效攻角小于失速极限假设,直接引用机翼作弯扭耦合振动时的Theodorson非定常气动力和气动力矩的表达式:

点击浏览公式

在局部气流坐标系中,h向下为正;气动力L向上为正;气动力距M使机翼前缘抬头为正; k=bω／V称为减缩频率或折合频率,是一个无量纲量;b是半弦长,C(k)=F(k)+iG(k)为Therdorson函数,是减缩频率k的复函数。

采用P-k法求解颤振方程时,认为在颤振临界附近,振动虽然不是简谐,但接近简谐,气动力计算时仍使用机翼作简谐振动的非定常气动力表达式。机翼的运动可以表述为:

点击浏览公式

由方程(11)得到每个片条的节点振动状态,经过坐标转换可得到局部气流坐标系下的振动:

点击浏览公式

局部气流坐标系由结构的静平衡状态确定,绕来流方向对总体坐标系旋转θ₃即可得到:

点击浏览公式

经过插值得到每个片条的振动:

点击浏览公式

根据方程(13)计算得到局部气流坐标系下的片条非定常气动力,然后根据功互等原理,将片条在局部气流坐标轴系下的非定常气动载荷,按照等效节点力的形式插值到相邻节点上,计算公式为:

点击浏览公式

在静平衡位置附近引入小扰动振动假设,根据方程(9),忽略阻尼后,机翼的气动弹性运动方程为:

点击浏览公式

Q(t)为机翼在静平衡位置附近振动时产生的非定常气动力。

引入广义位移q及实模态变换:

点击浏览公式

点击浏览公式

化简后,得到模态坐标系下的机翼气动弹性运动方程:

点击浏览公式

假设机翼分为N个片条,每个片条在其局部气流坐标系下具有沉浮和扭转2个自由度,对每个片条的非定常气动载荷分析后,组装总体气动力列向量,前乘以A^T后得到广义气动载荷矩阵,其它的求解与线性气动弹性系统P-k法的求解一致。

首先以文献标准算例为例,对大柔性机翼的气动弹性稳定性进行求解,以说明所发展的非线性气动弹性稳定性分析方法的有效性。待求解大柔性机翼的几何参数、刚度参数和质量参数如表 1所示:

在翼尖集中力的作用下,机翼分别产生了3.125％、6.25％、9.375％、12.5％、15.625％半展长的翼尖弯曲变形,如图 1所示:

图 1 非线性翼尖弯曲位移

图选项

定义由各个静平衡位置附近解得的准模态振动频率及振型为非线性振动频率和非线性振动振型;大柔性机翼的静变形改变了其惯性特性及刚度特性,非线性振动频率随变形的增加而变化的趋势如图 2所示。

图 2 不同变形基准下的准模态频率

图选项

大柔性机翼的几何大变形严重改变了其振动特性,因此采用线弹性理论求解这类机翼的动气动弹性问题时将会产生严重的偏差。采用P-k法对方程(20)进行特征值跟踪,得到各个静变形基准下的V-g和V-f曲线,进而得到机翼不同静平衡状态下的非线性颤振速度和频率,与文献值的对比见表 2。

从上表中可以看出,采用本文所发展的非线性气动弹性稳定性分析方法与文献^[14]的计算结果吻合的比较好,说明所发展的气动弹性稳定性分析方法具有较好的精度。

以类“太阳神”布局大柔性太阳能无人机为例,研究几何大变形对大柔性太阳能无人机气动弹性稳定性的影响。其几何、结构参数及飞行条件如图 3及表 3所示:

图 3 太阳能无人机几何模型

图选项

按照文献^[3]介绍的变形求解方法,解得无人机在气动载荷作用下的静变形,以翼尖变形分别为展长的5％、10％、15％、20％为例研究几何大变形对大柔性太阳能无人机气动弹性稳定性的影响。无人机变形后的位形如图 4所示。

图 4 太阳能无人机变形后位形图

图选项

当翼尖弯曲变形大于展长的5％时,靠近翼尖的结构变形中刚体旋转效应变的越来越明显,这是几何大变形问题的主要特征之一。另外,局部气动载荷的作用方向也发生较大变化。

在静平衡位置附近,引入准模态假设解得的大柔性太阳能无人机的非线性振动频率如图 5所示。初位形时结构的非线性振动频率与振型与线性解是一致的。

图 5 准模态频率

图选项

图 6 非线性颤振分支V-g曲线

图选项

图 7 非线性颤振分支V-f曲线

图选项

随变形的增加,扭转为主的非线性振动频率降低的比较快,可见几何大变形影响此太阳能无人机的气动弹性稳定性的主要途径之一就是改变了其扭转振动特性。

根据初位形解得的非线性颤振速度与频率分别为73.17 m／s和1.626 Hz;翼尖挠度为展长的20％时,非线性颤振速度降低到66.86／ s,非线性颤振频率降低到1.527 Hz。显然,采用线性理论解此类非线性气动弹性稳定性问题将引起较大偏差。太阳能无人机一般具有较大尺度的翼展,并且允许较大的弹性变形,气动弹性设计时需要着重考虑几何大变形的影响。

将气动弹性系统中的弦向弯曲刚度、扭转刚度、弹性轴位置、剖面质心位置等参数作为变量,以第3.2节中的变形作为基准,即静平衡状态下的变形量作为已知条件,反求已知变形下不同刚度参数的切线刚度矩阵;采用第2节中所发展的气动弹性稳定性分析方法,研究上述参数变化对无人机气动弹性稳定性的影响。

定义弦向弯曲刚度与垂直弯曲刚度的比值为弯曲刚度比:

点击浏览公式

根据定义,3.2节中模型的η_w值为200;在3.2节的模型基础上,改变弦向弯曲刚度进而改变弯曲刚度比,得到的非线性颤振速度变化趋势如图 8所示。

图 8 非线性颤振速度随弯曲刚度比的变化

图选项

较小的弯曲刚度比时,随着水平弯曲刚度的增加,气动弹性稳定性对刚度比的敏感度非常大,在设计时需要注意弦向弯曲刚度与垂直弯曲刚度的分配。当50<η_w时,所有变形状态下的气动弹性稳定性对弯曲刚度比的敏感度都比较低,弦向弯曲刚度在很大范围内对气动弹性稳定性的影响都很小,此时弦向弯曲刚度具有较大的设计空间。另外,当50<η_w时,弹性变形的增加降低了非线性颤振速度,与3.2节中得到的结论一致。

按照(24)式定义扭转刚度比,以3.2节中的扭转刚度为基准,研究不同扭转刚度下的太阳能无人机气动弹性稳定性特征。

点击浏览公式

从图 9中可以看出，扭转刚度的变化对太阳能无人机气动弹性稳定性的影响非常显著；在3.2节中结构模型的基础上，增加扭转刚度可以有效改善其气动弹性稳定性；可见，在这类大柔性太阳能无人机的初步设计阶段，需要着重注意扭转刚度的设计。另外，随着扭转刚度的增加，几何大变形对气动弹性稳定性的影响也越来越大，并且在很宽的扭转刚度变化范围内，几何大变形的影响都是非常显著的，因此这类太阳能无人机的扭转刚度设计时还需要考虑机翼几何大变形的影响。

图 9 非线性颤振速度随扭转刚度比的变化

图选项

大柔性太阳能无人机的弹性轴与质心轴在弦向的站位是非常重要的气动弹性设计量。在3.2节研究中，剖面质心与刚心均位于30%弦长处，这里依然假定剖面质心与刚心重合，以其在弦向的站位作为变量，研究剖面刚心与质心重合时的弦向站位对大柔性太阳能无人机的气动弹性稳定性的影响。

从图 10中可以看出，弹性轴的后移严重削弱了大柔性太阳能无人机的气动弹性稳定性特性，例如在50%弦长处的非线性颤振速度相当于25%弦长处的25%左右；弹性轴在24%至32%弦长范围内，非线性颤振速度随弹性轴的后移降低的非常快；在32%至50%弦长范围内，非线性颤振速度降低的比较慢。可见，弹性轴越靠近前缘对气动弹性稳定性越有利。需要注意的是，在研究的弹性轴站位范围内，几何大变形对气动弹性稳定性均产生不利的影响。

图 10 非线性颤振速度随弹性轴站位的变化

图选项

进一步考虑剖面质心在弦向的站位对气动弹性稳定性的影响。弹性轴选择在50%的弦长处并固定不变，以剖面质心在弦向的站位作为变量并考虑几何大变形效应，研究剖面质心对大柔性太阳能无人机气动弹性稳定性的影响。

由图 11中的计算结果可以看出，剖面质心前移15%的弦长时，非线性颤振速度增加大约100%；并且剖面质心越靠近前缘，几何大变形的影响越小。

图 11 剖面质心站位对非线性颤振速度的影响

图选项

本文基于CR有限元动力学模型和片条非定常气动力模型，建立了适用于几何大变形大柔性飞机气动弹性稳定性研究的分析模型；对大柔性太阳能无人机的气动弹性稳定性问题进行了深入的研究，主要得到以下结论：

1)基于CR理论的动力学建模，可以考虑几何大变形对无人机的刚度、惯量特性的影响，基于此动力学模型所发展的非线性气动弹性稳定性分析方法具有较好的精度，可用于大柔性飞行器的气动弹性稳定性研究。

2)太阳能无人机一般具有较大尺度的翼展，并且允许较大的弹性变形，而几何大变形对无人机的气动弹性稳定性产生了不利影响，其气动弹性设计需要着重考虑几何大变形的影响。

3)结构刚度的分配及弹性轴和剖面质心在弦向的站位等因素对气动弹性稳定性的影响非常显著，几何大变形引起的不利影响可以通过上述参数的合理设计得到改善。