一种基于CR理论的大柔性机翼非线性气动弹性求解方法

一种基于CR理论的大柔性机翼非线性气动弹性求解方法

王伟1，周洲1，祝小平2，王睿1

(1.西北工业大学航空学院，西安710072；2. 西北工业大学无人机研究所，西安710065)

摘要：大展弦比大柔性机翼在气动载荷的作用下，产生较大的弹性变形，其惯性特性、刚度特性、动气动弹性特性等亦发生较大改变，常规的线性气动弹性分析方法不再适用。基于Co-rotational(CR)理论，推导了机翼变形后的切线刚度矩阵和质量矩阵，建立了考虑几何非线性效应的大柔性机翼结构动力学模型；耦合改进的ONERA非线性非定常气动力模型，提出了一种适用于大柔性机翼的非线性气动弹性求解方法。采用Newmark直接数值积分法及松耦合技术在时域内对气动弹性运动方程进行求解，对所提出的非线性气动弹性求解方法的正确性和精度进行了验证，并研究了大柔性机翼的极限环颤振特性。研究表明：适用于大柔性机翼完整的非线性气动弹性建模需要考虑机翼结构大变形和非定常气动力动态失速等非线性因素；弯曲变形可降低临界极限环颤振速度的15%以上，而前移弹性轴能够有效的提高临界极限环颤振速度；所提出的非线性气动弹性求解方法具有较好的精度和效率，满足大柔性机翼非线性气动弹性的求解需求。

关键词：非线性气动弹性；极限环颤振；CR理论；非定常气动力；动态失速；Newmark积分法

中图分类号:V211.47

文献标志码：A

DOI:10.13465/j.cnki.jvs.2015.19.010

Abstract:Very flexible wings under aerodynamic loads tend to produce larger deformation, it results in significant changes in inertial and stiffness characteristics, and dynamic aeroelastic features, the linear aeroelastic analysis method is no longer applicable. Here, based on the co-rotational(CR) theory, the tangent stiffness matrix and mass matrix of a wing after deformation were derived, the structural dynamic model of very flexible wings considering geometric nonlinearity was then established. Coupled with ONERA dynamic stall model, an efficient method to solve nonlinear aeroelasticity of very flexible wings was proposed. Using Newmark direct integration method and loose coupled algorithms, a numerical procedure for solving nonlinear aeroelastic dynamic equations was presented, and the efficiency and precision of the method were verified through tests. The results showed that structural and aerodynamic nonlinearities should be considered for complete nonlinear dynamic aeroelastic simulations of very flexible wings; the wing’s critical limit cycle oscillation speed decreases 15% or more due to its bending deformation, but it increases through shifting forward the wing’s elastic axis; the proposed method has a good precision and efficiency, and meets requirements of nonlinear aeroelastic analysis of very flexible wings.

基金项目：国家自然科学基金(51475387)

收稿日期：2014-06-12修改稿收到日期：2014-09-30

A CR theory-based approach for solving nonlinear aeroelasticity of very flexible wings

WANGWei1,ZHOUZhou1,ZHUXiao-ping2,WANGRui1(1. College of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China；2. UAV Research Institute, Northwestern Polytechnical University, Xi’an 710065, China)

Key words:nonlinear aeroelasticity; limit cycle oscillation; CR theory; unsteady aerodynamics loads; dynamic stall; Newmark integration method

高空超长航时太阳能无人机在飞行高度上能够填补低轨道卫星和常规能源飞机的空白，执行环境监测、通信中继等任务，具有覆盖区域广、持续时间长和花费低等特点，由于高空长航时的性能要求，这类飞机一般具有较小的结构重量和超大展弦比机翼，在气动载荷的作用下，机翼弯曲变形可达半翼展的50%，机翼的结构动力学特性和气动弹性特性随变形的增加而产生较大改变。建立高精度的大柔性机翼结构动力学模型进而尽可能高精度、高效率的预测大柔性大变形机翼的非线性气动弹性特征是这类飞机的主要研究内容之一。由于面内、面外弯曲变形和扭转变形与结构当前的几何位形有关，并且机翼局部剖面的有效攻角可能很大，这类机翼的气动弹性建模需要考虑结构几何非线性和大攻角时非定常气动力动态失速的影响。

目前考虑几何非线性效应的非线性气动弹性运动方程的求解方法主要有直接数值积分法和线性化后的频域法[1-6]。直接数值积分法可以考虑结构大振幅运动以及气动力的动态失速效应等非线性因素；线性化的频域法一般假设结构在静平衡位置附近作小振幅振动，然后耦合线性频域非定常气动力，采用经典的P-K法或V-g法求解而无法考虑气动非线性的影响[5]。Patil等[4]基于经典ONERA动态失速模型建立了可以描述二维翼段准静态失速的非线性非定常气动力模型，并耦合混合变分形式的弯-弯-扭耦合本征梁模型建立了非线性气动弹性运动方程，开展了针对大柔性机翼的非线性气动弹性的研究。Shams等[6]采用Duhamel积分形式的Wagner函数描述线性时域非定常气动力并耦合非线性Euler-Bernouli理论梁模型建立了描述大柔性机翼的非线性气动弹性运动方程，没有考虑动态失速效应引起的非线性效应。非线性理论梁模型固然具有高精度、高效率等优点，但物理变量不直接，通用性受到一定的限制。随着非线性有限元技术的发展，谢长川等[5]基于线性化的UL有限元理论对复杂结构机翼的非线性颤振问题开展了大量研究，气动弹性运动方程在频域内求解，而无法预测临界速度前后的结构运动行为；在时域内采用TL法或UL法求解结构动力学方程时为了获得较好的精度一般选择非常小的时间步长，从而显著地增加求解时间；随着CR有限元技术的发展，其静力学求解允许适度大的载荷步、动力学方程求解同样允许适度大的时间步，从而引起了许多学者基于CR列式法对有限旋转问题的研究[7-13]。作者在大柔性太阳能无人机的静气动弹性研究中，引入了CR非线性有限元梁模型描述机翼的几何大变形，减少了流固耦合迭代次数，取得了较好的效果[14]。

本文在所编写的大柔性结构几何非线性静力学求解程序的基础上，进一步推导了空间CR梁单元的结构动力学方程，首先以悬臂梁模型为例，研究了大柔性结构在外载荷作用下的强迫振动，验证了基于CR有限元理论所建立的结构动力学模型的精度、效率及所编写程序的正确性。其次，以简谐振荡的NACA0012翼型为例，采用改进的ONERA动态失速模型对典型翼段的非线性非定常气动力进行求解，验证了非定常非线性气动力模型对动态失速效应模拟的精度。最后，耦合上述结构动力学模型和非线性非定常气动力模型，基于松耦合求解技术，提出了一种适用于几何大变形大柔性机翼的非线性气动弹性求解方法，并研究了大柔性机翼的极限环颤振特性以及弹性轴站位对极限环颤振的影响等。

1大柔性机翼结构动力学建模

大柔性机翼在气动载荷的作用下，变形前后的几何位形具有明显的差异，常规的线弹性小变形假设不再成立，而局部结构的材料应力应变关系仍处于线弹性范围内。CR有限元法的主要思想是将结构的变形分为刚体的旋转及平移和单元坐标系内的线弹性变形，平衡方程在单元坐标系内建立，通过坐标转换矩阵及其微分形式构造总体坐标系下的切线刚度矩阵，因而其核心是坐标转换矩阵的建立及其微分形式的推导，其计算精度主要取决于刚度矩阵和坐标系转换矩阵对几何非线性描述的精确性。

1.1CR梁单元静力学平衡方程

CR有限元理论的平衡方程建立在单元坐标系内，通过总体坐标系Ixyz、节点坐标系ui和单元平均构形坐标系um建立单元坐标系ue，见图1。

结构位形确定后，单元节点的三个欧拉角和现时坐标随之确定，按照式(1)旋转总体坐标系得到节点坐标系：

ui(x)=Ti(x)Ixyz

(1)

那么参照文献[9,14]，CR空间梁单元的单元坐标系定义为：

ue1=(d2-d1)/ln

ue2=um2-((uTm2ue1)/2)(ue1+um1)

ue3=um3-((uTm3ue1)/2)(ue1+um1)

(2)

式中，d1与d2分别是节点i,j的现时坐标ln=[(d2-d1)T(d2-d1)]1/2。

图1　柔性结构变形示意图 Fig.1 Sketch of a deformed very flexible structure

可见节点坐标系和单元坐标系都是结构位形的函数，在求解过程中需要随变形的增加而不断更新。

CR有限元法充分利用了几何非线性问题的小应变、大位移特征，弹性变形在单元坐标系内描述，位移-应变关系仍然是线性的。利用方程2建立的单元坐标系与方程1建立的节点坐标系，得到单元坐标系下的节点位移：

ul=ln-lo

2θl1,1=-uTi3ue2+uTi2ue3

2θl1,2=-uTi1ue3+uTi3ue1

2θl1,3=-uTi2ue1+uTi1ue2

2θl2,1=-uTi3ue2+uTi2ue3

2θl2,2=-uTi1ue3+uTi3ue1

2θl2,3=-uTi2ue1+uTi1ue2

(3)

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

δpl=Tδp

(4)

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

f=TTfl=TTKlpl

(5)

δf=δTTKlpl+TTKlδpl=

[KTσ(fl)+TTKlT]δp

(6)

增量平衡方程为：

KTΔp=ΔFe

(7)

KT=KTσ(fl)+TTKlT为所求的切线刚度矩阵，其中TTKlT为材料刚度矩阵，KTσ(fl)为几何刚度矩阵，是单元坐标系下单元节点力fl的函数。

当大柔性机翼上布置有发动机等推进系统时，将会引入随动力进而影响机翼的颤振稳定性[17]。假设该发动机挂载在第i个节点下，Fs(t)表示总体坐标系下随动载荷等效到节点i处的外载荷。那么，第i个节点的总外载荷可以描述为：

Fi,e(t)=Fi,a(t)+Fs(t)

(8)

当随动载荷描述在局部坐标系下时，需要把Fs,i(t)转换到总体坐标系下，转换表达式为：

Fs(t)=Ui,IFs,i(t)

(9)

具体求解时需要注意，随结构运动需要实时更新转换矩阵Ui,I。

1.2CR梁单元动力学平衡方程

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

(10)

那么，截面惯性力和惯性力矩可以表示为：

(11)

截面刚心处的线加速度、角速度和角加速度依次通过节点线加速度、角速度和角加速度插值得到：

(12)

同样可以插值得到截面虚位移，那么惯性力在虚位移上所做的虚功为：

(13)

根据虚功原理，等效节点惯性力所做的虚功为：

δW=f Tmass(δuTi,δθTi,δuTj,δθTj)T

(14)

把方程13代入式(14)得：

由方程(15)可以看出，线加速度是定义在总体坐标系下的，而角加速度是定义在单元坐标系下的，并且惯性力与时间和结构位形有关。

根据方程(15)所得到的惯性力表达式，忽略阻尼后考虑几何非线性效应的大柔性结构动力学平衡方程可以表示为：

Fe-fmass-fi=0

(16)

式中，Fe表示结构t时刻所受的外载荷，fmass表示结构t时刻的惯性力，fi表示结构t时刻的内力，它们均与时间t有关。

1.3非线性动力学平衡方程的时域推进求解

采用隐式Newmark时域积分格式求解大柔性机翼的结构动力学问题时，需要考虑系统在tn+1时的平衡，并采用Newton-Raphson法迭代求解控制平衡方程。Newmark积分法可以认为是对线性加速度法的推广，其插值假设是：

(17)

式中，ΔΨ=UTnΔψ，为随动坐标系下描述节点坐标系Un旋转到Un+1所对应的旋转伪矢量；Δψ为总体坐标系下描述节点坐标系Un旋转到Un+1所对应的旋转伪矢量。

根据上述假设，对方程15进行微分，得：

δfmass=δ([Ue]fmass,in)=MT,massδp

(18)

根据方程(18)得到惯性力的切线表达形式，可以采用预估校正技术求解方程(16)。预估的本质是用tn时刻位移及增量斜率预测tn+1时刻的位移。第一次位移增量Δp0可以通过下式预估求解得到：

(19)

将预估求解得到的位移增量代入方程式(16)，计算tn+1时刻第k次迭代后的非平衡力：

(20)

校正步内采用下式迭代求解，直到满足收敛标准：

gk=(KkT,tn+1+MkT,tn+1)Δpk

(21)

与静力学增量平衡方程的求解基本一致，采用预估校正推进技术求解方程式(16)时，在校正步内一般采用Newton-Raphson或修正的Newton-Raphson迭代求解，只是动力学求解中残差gk随所求解的时间步内位移增量的更新而不断更新。

2非线性非定常气动力模型

采用ONERA动态失速模型描述非线性非定常气动力，在不同文献中，其表达式形式也不尽相同；这里采用改进的ONEAR(EDLin)动态失速模型[18]。非定常升力、力矩作用在四分之一弦长处，表达式为：

(22)

(23)

(24)

把式(23)和(24)写成紧凑矩阵形式：

(25)

(26)

把式(26)写成紧凑矩阵形式：

采用四级四阶Runge-Kutta公式求解改进的ONERA动态失速模型中的微分方程，求解格式为：

(27)

3非线性气动弹性运动方程

第二节中的动力学方程求解是以总体坐标系下的节点线位移、角位移、线速度、线加速度和单元坐标系下角速度、角加速度为未知量的；非线性非定常气动力的求解输入是局部气流坐标系下的片条沉浮线位移、线速度、线加速度、角位移、角速度、角加速度等，因此需要对动力学模型的解向量进行转换和插值。

首先绕Iy旋转90°，再绕现在的x轴旋转90°得到总体气流坐标系，再绕总体气流坐标系的负x轴旋转θ得到局部气流坐标系：

(28)

那么总体坐标系下的线位移、角位移、线加速度、角加速度和单元坐标系下的角速度及角加速度到局部气流坐标系下的转换关系式为：

(29)

经过插值得到每个片条的振动及沉浮、扭转运动：

(30)

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

(31)

采用下式把局部气流坐标系下的节点力转换到总体坐标系下：

(32)

对式(32)得到的节点气动力组装后可以得到总体坐标系下的非定常气动力Qe，代入方程(16)代替外载荷向量Fe，可以得到机翼的气动弹性运动方程：

Qe-fmass-fi=0

(33)

在时域内对式(33)进行直接数值积分求解时，与结构动力学求解中tn+1时刻的外载荷已知不同，tn+1时刻的非定常气动力Qe与tn+1时刻的结构运动状态有关，因而需要预估tn+1时刻的非定常气动力，然后根据解得的tn+1时刻的结构状态更新非定常气动力，根据Newton等距向前插值多项式进行外推：

Qe,tn+1=3Qe,tn-3Qe,tn-1+Qe,tn-2

(34)

松耦合的求解思想是：首先预估下一时刻的非定常气动力，然后带入气动弹性运动方程中求解结构运动状态，以解得的下一时刻结构状态重新计算非定常气动力代替预估值，但不参与迭代求解，而是直接推进到下一时间步的求解；若参与迭代求解，上述计算方法则变成紧耦合推进求解算法。紧耦合迭代求解算法具有精度高的优点但求解易发散，一般很少在非线性气动弹性系统的求解中使用。综上所述，松耦合推进求解方程式(33)的流程见图2。

图2　非线性气动弹性系统松耦合求解流程 Fig.2 Calculation process of the nonlinear aeroelastic system based on loosely coupled technology

4算例及验证

大柔性机翼非线性气动弹性求解所需的几何参数、刚度参数、质量参数及大气参数见表1。

4.1大柔性悬臂梁几何非线性动力学求解

首先把大柔性机翼作为悬臂梁，不考虑气动力的作用，在翼尖施加集中载荷历程见图3(b)，分别采用经典UL法(时间步长为1E-4s)和CR法(时间步长为5E-2s)进行求解，以验证基于CR有限元法采用直接数值积分技术(Newmark法)求解大柔性结构动力学响应的精度等。

表1　大柔性机翼模型参数 [4]

图3　大柔性悬臂梁及其端部加载示意图 Fig.3 Sketch of a very flexible cantilever beam and the tip loading history

图4　不同步长下CR法与UL法求解结果对比 Fig.4 Results comparing between CR method and UL method under variant time steps

由图4中的计算结果可以看出，采用CR法求解大柔性结构的动力学响应问题时允许采用较大的时间步长，且具有较好的求解精度，满足大柔性机翼结构动力学求解对几何非线性因素的描述精度要求。

4.2非线性非定常气动力模型验证

以NACA0012翼型为例，采用ONERA动态失速模型及其改进形式对非定常非线性气动力进行求解，结果见图5。

图5　α=15°+10° cos(0.1τ)的非定常气动力 Fig.5 Unsteady aerodynamics atα=15°+10° cos(0.1τ)

大迎角区域内，经典的ONERA(EDLin)模型求解精度较好，而改进后的ONERA模型虽然略微降低了对大迎角区的模拟精度但是显著的改善了小迎角区的求解精度。

仍然以NACA0012翼型为例，平均迎角分别取6°、12°、17°，h=0.1sin(0.1τ)m，τ=Vb/t；采用改进的ONERA动态失速模型研究完全线性段、部分线性段部分非线性段及完全非线性段内的非定常气动力特性，得到翼段纯沉浮运动下的非定常气动力见图6。

图6　h=0.1sin(0.1τ)m的非定常气动力 Fig.6 Unsteady aerodynamics at h=0.1 sin (0.1τ)m

可见采用改进的ONERA模型无论是线性段还是非线性段均可以较好的模拟二维机翼的非定常气动力，特别是非线性段的动态失速效应。

4.3大柔性机翼非线性气动弹性响应求解

以表1中所给出的大柔性机翼为例，机翼剖面为NACA0012翼型，初始翼尖弯曲位移为4m时，本文与文献4及文献6得到的临界极限环颤振速度均约为28m/s，翼尖位移响应历程与文献值的对比见图7示。

图7　V=28m/s，y 0=4.0m时翼尖弯曲位移响应 Fig.7 Wingtip displacement response at V=28m/s，y 0=4.0m

从图7可以看出，采用所提出的非线性气动弹性求解方法解得的位移响应与文献值6吻合的更好，可能是因为文献4中没有考虑非定常气动力的动态失速效应。翼尖位移响应的相位图见图8。

图8 V=28m/s,y0=4.0m翼尖相位图Fig.8Phase-planeplotsofthewingtipresponseatV=28m/s,y0=4.0m图9 y0=4.0m时不同来流速度下的翼尖弯曲位移响应Fig.9Wingtipdisplacementresponseaty0=4.0mandvariantairspeeds

图10　y 0=4.0m时不同来流速度下的翼尖响应相位图 Fig.10 Phase-plane plots of the wingtip response at y 0=4.0m and variant air speeds

由计算结果可以看出，采用所提出的非线性气动弹性求解方法可以较好的预测大柔性机翼的极限环颤振。进一步增加来流速度，研究来流速度大于临界速度时的大柔性机翼非线性气动弹性响应特性；初始弯曲位移为4.0m。

由计算结果可以看出，当来流速度略大于28m/s时，大柔性机翼的气动弹性翼尖位移响应仍具有周期性，当来流速度增加到32m/s时，翼尖位移响应不再具有周期性。随来流速度的增加，翼尖弯曲位移响应及扭转响应幅值随之增加，且极限环颤振变得越来越复杂。来流速度为28m/s时，极限环颤振每个周期的振动幅值一致，但来流速度增加后，翼尖位移响应需要经过一个幅值较小的周期和一个幅值较大的周期才能够回到原位置。另外，从图10中还可以看出随来流速度的增加，局部动态失速效应越来越明显。

4.4弹性轴站位对大柔性机翼极限环颤振的影响

给定初始条件下，诱发非线性气动弹性极限环颤振的最小速度称为临界极限环颤振速度[4]。结构参数仍以表一中给出的参数作为基准不变，研究弹性轴和剖面质心一起向前缘移动时对临界极限环颤振速度的影响。初始翼尖弯曲位移分别选为4m、3m、2m、1m，得到的临界极限环颤振速度随弹性轴在弦向站位的不同而增减的趋势见图11。

图11　弹性轴站位对极限环颤振的影响 Fig.11 The effects on limit cycle oscillation due to variant elastic axis positions

从图11中可以看出，翼尖弯曲位移的增加将降低机翼的临界极限环颤振速度；弹性轴在50%的弦长处时，初始翼尖弯曲位移为4m时解得的临界极限环颤振速度与初始翼尖弯曲位移为1m时相比降低了15.36%；另外，弹性轴的前移可以有效的提高临界极限环颤振速度；因此，对于这类大柔性机翼，通过合理设计弹性轴的弦向站位，可以有效的改善弯曲变形对机翼气动弹性所引起的不利效应。

4.5积分步长的选择及收敛性

直接数值积分的收敛性和精度可以通过不同时间步长的选取以证明，在上述研究中时间步长均为0.005s，这里进一步补充来流速度为28m/s时，初始弯曲位移为4m，时间步长为0.01s和0.0001s的计算结果。

图12　不同时间步长的积分结果 Fig.12 Comparing of direct integration results between variant time steps

从图12中可以看出，时间步长为0.005s和0.0001s的积分结果基本一致，因此时间步长选为0.005s满足计算精度的要求，并且与时间步长为0.0001s时相比大大减少了求解时间。针对具体的求解问题，对比不同的积分步长对求解精度和稳定性的影响，指导选择合适的积分步长，可以有效地提高求解效率；具体执行时，可选择下一次积分步长为本次积分步长的1/N(N为大于等于2的整数)，两次求解得到的最大响应幅值相差小于幅值的5%时，即可把本次试算积分步长作为采用直接数值积分求解目标算例的较为合理的时间积分步长。

5结论

本文基于CR有限元非线性动力学求解技术和改进的ONERA动态失速非线性非定常气动力模型提出了一种适用于大柔性大变形机翼的非线性气动弹性求解方法，进而对大柔性机翼的非线性气动弹性响应特性进行了较为深入的研究，主要得到以下结论：

(1)本文所提出的非线性气动弹性响应计算方法，以总体坐标系下的位移等为变量，物理意义简单直接，具有较好的求解精度和计算效率，能够满足大柔性机翼非线性气动弹性响应求解的需求，并可以较好的预测大柔性机翼的极限环颤振。

(2)较小的时间步长可以提高求解的精度但计算花费随之增加，在具体的求解中可以预先对比不同积分步长对求解精度的影响，指导选取合适的积分步长，以满足求解的需求。

(3)来流速度大于临界极限环颤振速度时，大柔性机翼的非线性气动弹性响应是非常复杂的，非定常气动力的动态失速效应虽然限制了振幅的增加，但也会带来结构疲劳等问题，在这类机翼的结构设计时需要特别注意。

(4)前移弹性轴可以有效地提高临界极限环颤振速度，增加周期性极限环颤振的速度区间，显著地改善大柔性机翼的非线性气动弹性特性，可以用于指导大柔性机翼弹性轴在弦向的站位设计。进一步的研究工作中可以引入高效的预估求解算子，减少动力学响应求解时的迭代次数，提高求解效率。

参考文献

[1]Patil M J, Hodges D H. On the importance of aerodynamic and structural geometrical nonlinearities in aeroelastic behavior of high-aspect-ration wings[J].Journal of Fluids and Structures,2004,19:905-915.

[2]Patil M J, Hodges D H, Cesnik C E S. Nonlinearaeroelasticity and flight dynamics of high-altitude long-endurance Aircraft [R]. AIAA-99-1470.

[3]Patil M J, Hodges D H. Flightdynamics of highly flexible flying wings [J]. Journal of Aircraft，2006，43(6):1790-1798.

[4]Patil M J, Hodges D J, Cesnik C E S. Limit cycle oscillations in high-aspect-ratio wings[J].Journal of Fluid and Structures,2001,15:107-132.

[5]Xie Chang-chuan, Yang Chao. Linearization methods of nonlinear aeroelastic stability for complete aircraft with high-aspect-ratio wings[J]. Sci China Tech Sci,2011, 54:403-411.

[6]Shams S, Sadr M H, Haddadpour H.An efficient method for nonlinear aeroelasticy of slender wings[J].Nonlinear Dynamics,2012,67:659-681.

[7]周凌远，李乔.基于UL法的CR列式三维梁单元计算方法[J]. 西南交通大学学报，2006，41(6):690-695.

ZHOU Ling-yuan, LI Qiao. Updated lagrangian co-rotational formulation for geometrically nonlinear FE analysis of 3D beam element[J].Journal of Southwest Jiaotong Unoversity,2006,41(6):690-695.

[8]Belytschko T, Schwer L. Large displacement, transient analysis of space frames[J]. International Journal for Numerical Methods in Engineering, 1977, 11:65-84.

[9]Crisfield M A. A consistent Co-rotational formulation for non-linear, three-dimensional, beam element[J]. Computer Methods In Applied Mechanics And Engineering, 1990,81:131-150.

[10]Crisfield M A. Non-linear finite element analysis of solids and structures, Volume 2: Advanced topics[M]. John Wiley & Sons, Chichester, New York, Weinheim, Brisbane, Singapore, Toronto,2000.

[11]Crisfield M A,Galvanetto U, Jelenicé G. Dynamics of 3-D co-rotational beams[J]. Computational Mechanics, 1997,20：507-519.

[12]Ghosh S, Roy D. Consistent quaternion interpolation for objective finite element approximation of geometrically exact beams[J]. Computer Methods In Applied Mechanics And Engineering,2008, 198:555-571.

[13]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.

[14]王伟,周洲,祝小平,等.考虑几何非线性效应的大柔性太阳能无人机静气动弹性分析[J]. 西北工业大学学报, 2014, 32(4):499-504.

WANG Wei, ZHOU Zhou, ZHU Xiao-ping,et al. 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.

[15]Palacios R, Murua J, Cook R. Structural and aerodynamic models in nonlinear flight dynamics of very flexible aircraft [J]. AIAA Journal, 2010,48 (11):2648-2659.

[16]吴国荣,颜桂云，陈福全.柔性梁大柔度动力响应分析的多体系统方法[J].振动与冲击,2007,26(3):87-89.

WU Guo-rong,YAN Gui-yun,CHEN Fu-quan.Large deflection dynamic response analysis of flexible beams bi multibody system method[J].Journal of Vibration and Shock, 2007,26(3):87-89.

[17]张健,向锦武.侧向随动力作用下大展弦比柔性机翼的稳定性[J].航空学报,2010,31(11):2115-2123.

ZHANG Jian, XIANG Jin-wu. Stability of high-aspect-ratio fiexible wings loaded by a lateral follower force[J].Chinese Journal of Aeronautics,2010,31(11):2115-2123.

[18]Laxman V,Venkatesan C.Chaotic response of an airfoil due to aeroelastic coupling and dynamic stall[J].AIAA Jourmal, 2007,45(1):271-280.

第一作者何洪阳男，硕士，1988年生

通信作者陈春俊男，博士，教授，1967年生