基于曲梁模型的大展弦比大柔性机翼颤振分析

段静波, 周洲, 江涛

(1.西北工业大学 航空学院, 陕西 西安 710072; 2.军械工程学院 无人机工程系, 河北 石家庄 050003)



基于曲梁模型的大展弦比大柔性机翼颤振分析

段静波1,2, 周洲1, 江涛2

(1.西北工业大学 航空学院, 陕西 西安 710072; 2.军械工程学院 无人机工程系, 河北 石家庄 050003)

提出一种大展弦比大柔性机翼颤振分析的方法。该方法首先引入准模态假设，将气动载荷作用下发生大静变形的大展弦比大柔性机翼视为一根变曲率曲梁，并将其离散为一系列常曲率曲梁单元，利用机翼静变形结果，通过多项式插值获得各曲梁单元的平均曲率。其次，在曲梁单元内，利用曲梁振动微分方程和Therdorson非定常气动力模型建立曲梁单元的颤振微分方程。然后，运用传递函数法，将曲梁单元的颤振微分方程转换为状态空间形式，并依照有限元组集的思想形成机翼整体平衡方程。最后，通过求解特征值问题获得机翼的颤振速度和颤振频率。通过与已有文献结果的对比，验证了新方法的正确性和有效性。

大展弦比；大柔性机翼；曲梁；非定常气动力；颤振

大展弦比大柔性机翼在气动载荷作用下通常经历较大的结构变形。准确把握这类大变形机翼的气动弹性特性是目前研究的热点之一。国内外学者开展了较多卓有成效的研究工作。Smith Patil等[1-2]较早地开展了基于CFD方法的大展弦比机翼气动弹性特性研究。机翼由Hodges-Dowell梁来模拟，气动力计算则以Euler方程为基础。Garcia[3]应用精确梁理论的有限元方法，并耦合N-S方程的气动力，计算了大展弦比机翼的跨音速静气动弹性特性。Palacois、Beran等[4-5]基于流固耦合方法，分别研究了超大展弦比机翼和翼身融合体飞机的静气动弹性问题。国内方面，谢长川等[6-7]应用“准模态”假设，分析结构几何非线性对大展弦比机翼振动的影响，即认为结构是在大的静变形平衡位置附近作微幅振动，从而可沿用线性系统振动理论进行机翼振动描述，通过算例验证了方法的有效性。周洲等[8-9]基于CR有限元理论，推导了大柔性机翼的切线刚度矩阵和质量矩阵，建立了考虑几何非线性效应的大柔性无人机结构动力学模型，并引入准模态小扰动振动假设，对动力学模型进行线性化，并采用建立在局部气流坐标系下的Therdorson片条非定常气动力，通过p-k法对无人机的非线性颤振速度和频率进行了求解。杨智春和党会学等[10]耦合Euler方程求解器和非线性结构求解器，计算了大展弦比机翼的静变形，并在静变形的基础上，提取结构的剩余刚度进行了非线性颤振特性分析。向锦武等[11]研究了侧向随动力作用下大展弦柔性机翼气动弹性稳定性等问题。

本文在前人研究的基础上，发展了一种大展弦比大柔性机翼颤振分析的方法，通过与已有文献结果的对比，验证了本文方法的正确性和有效性。

1 机翼颤振微分方程建立

1.1 机翼大静变形曲线方程拟合

大展弦比大柔性机翼静变形后可视为一根曲梁。梁横截面的刚心轴线形成一条曲线，在直角坐标系XOY下如图1所示。

图1 变曲率曲梁的示意图

将曲梁划分若干段，每段为一个曲梁单元。在直角坐标系下，每个曲梁单元的起点和终点的坐标以及斜率均可根据机翼静变形给出，从而通过多项式插值可将曲梁单元的曲线方程表示出来，为[12]：

(1)

由曲线的曲率公式可知曲梁单元内任一点的曲率半径为

(2)

每个曲梁单元内的曲率可近似为常数,可取单元两端处的曲率半径平均值作为该单元的平均曲率半径。

1.2 曲梁单元振动微分方程

在曲梁单元内,曲梁的曲率可视为常数R,其自然坐标系如图2所示,x轴表示切向,y轴表示径向,z轴表示竖向。根据文献[13],考虑曲梁质心轴与弹性轴不重合的情形,可获得曲梁六自由度振动方程为

(3)

1.3 非定常气动力模型

在忽略机翼重力影响条件下,机翼颤振时的外力为气动力产生的分布升力以及分布扭矩。本文采用片条理论进行非定常气动力计算。气动力片条在曲梁形状基础位置上定义,且位于曲梁单元中点处。根据Theodorson理论,单位展长的非定常升力与相应的俯仰力矩按下式计算[14]

(4)

式中,k=ωb/V为折合频率,是一个无量纲参数,C(k)为Theodorsen函数,由于k是ω和V的函数,为了后续求解方便,将C(k)写为C(ω,V),其他符号含义同文献[8]。

1.4 机翼单元颤振微分方程

将(4)式代入(3)式后就可得机翼单元颤振时的微分方程

(5)

图2 曲梁单元自然坐标系

2 机翼颤振的传递函数求解

2.1 机翼单元的传递函数

对机翼单元颤振微分方程(6)式进行Fourier变换,并整理可得

定义单元状态向量

(7)

将(6)式改写为状态空间形式的方程

(8)

式中,ge(x,ω)=0,转移矩阵Fe(ω,V)为12×12的方阵,其非零元素为

边界条件为

(9)

式中,l为曲梁单元长度,Mb、Nb分别为单元边界条件选择矩阵,为12×12的方阵,其非零元素为

Nb(10,8)=1,Nb(11,10)=1,Nb(12,12)=1

γe(ω)为由位移或力组成的列向量,其表达式为

方程(8)式的传递函数解为

(10)

式中

(11)

2.2 系统的合成与求解

由于大变形后的机翼划分为若干常曲率曲梁单元来描述整个机翼,其求解方程可借鉴有限元方法的思想进行组集。具体如下:

曲梁单元截面上的内力为

(12)

式中,My、Mz分别为曲梁单元绕x、y、z轴弯矩,Tx分别为曲梁单元绕x轴的扭矩,Nx为沿x轴的轴力,Qy和Qz分别为沿y、z轴的剪力。

将曲梁单元内力写成矩阵形式

(13)

将(10)式代入(13)式后,取x=0、l,可得到曲梁单元两端点处的内力

(14)

式中

上式与有限元法中单元节点力表达式十分相似,fe可视为单元节点内力,Ke(ω,V)可视为单元刚度矩阵,γe(ω)可视为单元节点位移向量。

按照有限元组集方法对各节点统一编号,并进行组集拼接,可得到机翼整体平衡方程为

(15)

式中,K(ω,V)可视为整体刚度矩阵,γ(ω)可视为整体节点位移向量,f为各单元节点内力拼装成的向量。由于本文中将机翼的气动力与机翼作为一个完整的系统来考虑,除此之外,机翼没有受到其它外力作用,因而根据单元节点内力与外载荷平衡,可得出

(16)

根据机翼约束条件,按照有限元方法对整体刚度矩阵K(ω,V)进行边界条件处理[15-16]。

当机翼颤振时,γ(ω)应有非零解,此时须满足条件

(17)

由于K(ω,V)为复矩阵,其行列式值等于零的必要条件为矩阵行列式值的实部与虚部均为零,即

(18)

(18)式2个方程包含2个未知变量V和ω,因而求解(18)式可得到V和ω的解。其中,V即为机翼颤振速度,ω即为机翼颤振频率。通常满足(18)式条件的解可能有多个,其中V最小的一组解即为机翼的颤振特性。

3 结果与分析

3.1 方法验证

为了验证本文方法的正确性,按照文献[8]算例计算大展弦比柔性机翼的气动弹性稳定性。大柔性大展弦比机翼的几何参数、刚度参数和质量参数均与文献[8]一致,见表1所示。

表1 机翼模型参数

通过表1中机翼模型参数,假设机翼为平板翼型可反推出本文方法所需要的参数。

文献[8]给出了机翼翼尖作用集中力下的机翼挠度曲线。本文取翼尖位移分别为机翼半展长的3.125%、6.250%、12.50%的变形情形,分别记为变形一、变形二、变形三。将变形后的机翼视为曲梁,分别拟合出机翼变形后的曲线方程,沿机翼展向划分10个单元后,按(2)式计算每个单元的平均曲率半径,然后利用本文方法计算机翼颤振特性。

表2 机翼颤振结果

图3至图5给出了大展弦比大柔性机翼在不同程度静变形下求解机翼的颤振速度和颤振频率的过程。

图3 Re[det(K)]和Im[det(K)]的等值线图(变形一)

图3为机翼变形一时,(18)式中Re[det(K)]和Im[det(K)]的等值线图。图中Re[det(K)]和Im[det(K)]的零等值线交点即为满足(18)式的解。从图中可以看到,同时满足(18)式且V最小的一组解(V,ω)=(34.5,24.4),即为图3中O点。这表明,机翼变形一时,机翼的颤振速度为34.5 m/s,机翼的颤振频率为24.4 rad/s。

图4 Re[det(K)]和Im[det(K)]的等值线图(变形二)

图5 Re[det(K)]和Im[det(K)]的等值线图(变形三)

同样,图4和图5分别为机翼变形二、变形三时的Re[det(K)]和Im[det(K)]等值线图,由图中O点可以得到机翼的颤振特性。

为了便于比较，将计算得到的机翼3种变形下的颤振速度和颤振频率列入表2中，并与已有文献结果进行对比。从表中可以看出，本文方法与文献的结果基本吻合。两者存在一定差异的主要原因在于：文献中进行气弹分析时，在结构方面进行了动力学降阶，利用了机翼若干低阶模态，模态的选用可能会影响结果的精度；在气动力模型方面，采用的模型考虑了机翼的三维效应。本文方法在结构方面直接采用微分方程，避免了模态降阶，在气动力模型方面采用了片条理论，没有考虑机翼的三维效应。

此外，从表2中还可以看出，随着机翼变形量增大，本文结果与参考文献结果的相对误差逐渐减小。这其中主要的原因是，用曲梁模型模拟机翼时，当机翼变形较小时，曲梁单元的曲率半径误差较大，随着机翼变形增大，曲梁单元的曲率半径误差逐渐减小。

3.2 机翼颤振特性影响分析

采用本文方法，将划分单元数目、弯曲刚度比、扭转刚度、质心弦向位置等作为变量，以上节中的3种变形情形作为基准，研究上述因素对大变形机翼颤振特性的影响。

表3给出了单元数目对机翼颤振速度和颤振频率收敛性的影响。从表中可以看出，当采用10个单元后，结果已收敛。这表明本文方法使用较少单元即可较快收敛，因而非常适合总体设计阶段的机翼颤振特性快速分析要求。

表3 单元数目对机翼颤振结果的影响

图6给出了3种不同变形情形下扭转刚度变化对机翼颤振速度的影响。从图中可以看到，随着扭转刚度增大，不同变形条件下机翼颤振速度都有了提高。但是，机翼变形程度不同时，机翼颤振速度增加是有差别的。在变形一情形下，机翼变形相对较小，机翼颤振速度增幅相对较大；变形三情形下，机翼相对变形较大，机翼颤振速度增幅相对降低。这表明，随着机翼变形越来越大，通过提高扭转刚度改善机翼颤振特性的效果会有所降低。主要原因是：在机翼大变形下，靠近翼尖的结构大变形会显著改变机翼的扭转振动特性，进而影响扭转刚度增加改善机翼颤振特性的效果。

图7给出了3种不同变形情形下机翼弦向弯曲刚度与垂直弯曲刚度比对机翼颤振速度的影响。从图中可以看到，在弦向弯曲刚度与垂直弯曲刚度比较小时，3种变形状态下的机翼颤振速度对弯曲刚度比均比较敏感度，而且机翼颤振速度在机翼变形相对小时更敏感；当弯曲刚度比较大时，机翼颤振速度对其敏感度逐渐降低，此时弦向弯曲刚度具有较大的设计空间。

质心轴在弦向的位置是非常重要的气动弹性设计量，在上节研究中，剖面质心与刚心均位于弦长中点处，这里仍然假定两者重合，以其在弦向的位置作为变量，研究质心弦向位置对大变形机翼颤振特性的影响。图8给出了3种不同变形情形下质心弦向位置对机翼颤振速度的影响。从图中可以看到，在机翼弦向中点前后10%半展长范围内，随着质心沿弦向从前向后移动，机翼颤振速度逐渐降低。这一结论与刚性机翼颤振特性相似。但是，随着机翼变形越来越大，机翼颤振速度降低幅度加快。图中可以看，变形三情形下机翼颤振速度下降快于变形一、变形二。

图6 扭转刚度对机翼颤振速度的影响 图7 弯曲刚度比对机翼颤振速度的影响 图8 质心弦向位置对机翼颤振速度的影响

4 结 论

本文提出了一种大展弦比大柔性机翼颤振的分析方法。该方法将大变形后机翼视为一根变曲率曲梁，从而不仅可以通过曲梁形状考虑大展弦比大柔性机翼变形中的几何非线性因素，而且仍然可以运用曲梁线性振动理论来近似描述大变形机翼的动力学特性。在求解机翼的颤振特性时，采用传递函数法进行求解，方法计算规模小、收敛快，非常适合工程总体设计阶段的快速分析要求，而且该方法还可进一步推广应用于机翼气动弹性响应计算。

[1] Smith M J, Patil M J, Hodges D H. CFD-Based Analysis of Nonlinear Aeroelastic Behavior of High-Aspectratio Wings[R]. AIAA-2001-1582

[2] Patil M J, Hodges D H. Limit-Cycle Oscillations in High-Aspect-Ratio Wings[J]. Journal of Fluid and Structures, 2001,15:07-132

[3] Garcia J A. Numerical Investigation of Nonlinear Aeroelastic Effects on Flexible High-Aspect Ratio Wings[J]. Journal of Aircraft, 2005, 42(4):1025-1036

[4] Palacios R, Cesni C S. Static Nonlinear Aeroelasticity of Flexible Slender Wings in Compressible Flow[R]. AIAA-2005-1945

[5] Beran P S, Hur J Y, Snder R D. Static Nonlinear Aeroelastic Analysis of a Blended wing Body[R]. AIAA-2005-1944

[6] 谢长川,吴志刚,杨超. 大展弦比柔性机翼的气动弹性分析[J]. 北京航空航天大学学报, 2003, 29(12): 1087-1091

Xie Changchuan, Wu Zhigang, Yang Chao. Aeroelastic Analysis of Flexible Large Aspect Ratio Wing[J]. Journal of Beijing University of Aeronautics and Astronautics, 2003, 29(12):1087-1091 (in Chinese)

[7] Xie Changchuan, Leng Jiazhen, Yang Chao. Geometrical Nonlinear Aeroelastic Stability Analysis of a Composite High-Aspect-Ratio Wing[C]∥International Conference on Engineering Dynamics, 2007

[8] 王伟,周洲,祝小平,王睿. 几何大变形太阳能无人机非线性气动弹性稳定性研究[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)

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

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)

[10] 杨智春,党会学,李毅. 大展弦比机翼非线性气动弹性特性的数值模拟研究[C]∥第十一届全国空气弹性学术交流会, 2009

Yang Zhichun, Dang Huixue, Li Yi. Numerical Simulation of Nonlinear Aeroelastic Characteristics for a High-Aspect-Ratio Wing[C]∥11th Conference on Aeroelasticity, 2009 (in Chinese)

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

Zhang Jian, Xiang Jingwu. Stability of High-Aspect-Ratio Flexible Wings Loaded By a Lateral Follower Force[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(11):2115-2124 (in Chinese)

[12] 蒋纯志. 分布传递函数方法在梁杆结构分析中的应用[D]. 长沙：国防科技大学, 2006

Jiang Chunzhi. The Application of Distributed Transfer Function Method to Girders and Beams Analysis[D]. Changsha, National University of Defense Technology, 2006 (in Chinese)

[13] 赵雪健. 平面曲梁自由振动的动力刚度法研究[D]. 北京：清华大学，2010

Zhao Xuejian. Research on Dynamic Stiffness Method for Free Vibration of Planar Curve Beams[D]. Beijing, Tingshua University, 2013 (in Chinese)

[14] 赵永辉. 气动弹性力学与控制[M]. 北京: 科学出版社, 2006

Zhao Yonghun. Aeroelastic Mechanics and Control[M]. Beijing, Science Press, 2006 (in Chinese)

[15] 李海阳. 数学物理问题的数值分布传递函数方法[D]. 长沙:国防科技大学,1999

Li Haiyang. Numerical Transfer Function Method for Mathematic and Physical Problems[D]. Changsha, National University of Defense Technology, 1999 (in Chinese)

[16] 雷勇军. 结构分析的分布参数传递函数方法[D]. 长沙:国防科技大学,1998

Lei Yongjun. Distributed Transfer Function Method for Structural Analysis[D]. Changsha, National University of Defense Technology, 1998 (in Chinese)

A Method for Flutter of the Very Flexible Wing with High-Aspect-Ratio Based on Curve Beam Model

Duan Jingbo1,2, Zhou Zhou1, Jiang Tao2

1.School of Aeronautics, Northwestern Polytechnical University, Xi′an 710072, China 2.Department of UAV Engineering, Ordnance Engineering College, Shijiazhuang 050003, China

A method for flutter of the very flexible wing with a high-aspect-ratio is developed. Firstly, the very flexible wing subjected to aeroelastic load will undergo a large deformation. So the deformed wing can be regarded as a curve beam and divided into a couple of curve beam elements. The mean curvature of every beam element is obtained by mean of polynomial interpolation method. For each element, the flutter differential equations is established by combining the differential equations of the constant curvature beam vibration and the Therdorson' unsteady aerodynamics model. Then, using the distributed transfer function method and the finite element method, the equilibrium equations of the whole wing are obtained. Finally, the wing flutter was carried out by solving a eigenvalue problem. The results are good agreement to the literature solutions and indicate that the present method is accurate and efficient.

high-aspect-ratio; very flexible wing; curve beams; unsteady aerodynamics; flutter

2016-03-02

中国博士后基金(2014M560803)资助

段静波(1982―)，西北工业大学博士后，主要从事飞行器气动弹性研究。

V215.3

A

1000-2758(2016)05-0774-09