采用等效刚度有限元模型的复合材料机翼颤振分析

欧阳星， 余雄庆， 王 宇

(南京航空航天大学飞行器先进设计技术国防重点学科实验室， 江苏 南京 210016)

采用等效刚度有限元模型的复合材料机翼颤振分析

欧阳星， 余雄庆， 王 宇

(南京航空航天大学飞行器先进设计技术国防重点学科实验室， 江苏 南京 210016)

采用等效刚度方法，研究了一种适用于复合材料机翼初步设计阶段的动力学和颤振分析的简化结构有限元模型。该方法首先计算不同布局形式的加筋壁板的刚度矩阵，然后将其赋予与加筋壁板平面形状相同的光板(等效板)上，使加筋壁板和等效板具相同的力学性能。以某客机概念方案的机翼为例，建立了反映实际结构详细有限元及其等效刚度有限元模型。通过编程实现了自动化生成等效模型，极大地减少了建模和分析时间。计算结果和对比分析表明，两种模型的固有频率、振动模态和颤振分析结果吻合得很好，从而验证了等效刚度方法在机翼结构动力学和颤振分析方面的准确性。

机翼； 颤振； 复合材料； 有限元； 等效刚度

引 言

机翼主承力结构采用复合材料可有效减轻重量，是现代飞机先进性的重要标志之一。机翼结构必须保证在整个飞行包线不能发生颤振。这一设计要求对机翼结构方案设计有重要影响，因此在机翼的初步方案设计中，就必须对机翼的颤振特性进行分析。

结构动力学分析是颤振分析的基础。适用于复合材料机翼颤振分析的结构动力学分析模型可概括为三类：1)等效梁模型[1-2]；2)等效板模型[3-5]；3)有限元模型[6-8]。等效梁模型是将机翼简化为梁，主要适用于展弦比较大的机翼。等效板模型是将机翼等效至与机翼相同平面形状的一组梯形板上，梯形板的力学性能逼近真实机翼，多用于小展弦比机翼的动力学分析。有限元方法能详细描述机翼结构特征，具有通用和计算精度高等优点。

目前飞机设计的一个重要趋势是在飞机早期设计阶段采用高精度分析方法[9]。这是因为：1)高精度分析模型可用于新概念飞机方案(如翼身融合或联接翼方案)；2)能提高分析结果的可信度，减少设计方案的风险。

有限元模型是结构动力学分析的高精度模型。但详细的有限元模型复杂，其建模过程繁杂，所需计算时间长。由于在飞机初步设计阶段，机翼外形和结构方案还在不断调整之中，需要一种快速的颤振评估方法，而有限元模型的上述缺点阻碍了其在初步设计阶段的应用。解决这个问题的一个措施是采用等效刚度方法对结构有限元模型进行简化。

等效刚度法是将加筋壁板等效到一块相同形状的光板(无加强筋的板)上，并使光板与加筋壁板具有相同的力学性能。基于等效刚度法的有限元模型(以下简称等效刚度有限元模型)能模拟各种加强筋的力学特性，同时避免了对加强筋的有限元建模，大大简化了有限元模型，易于实施参数化机翼结构有限元建模。该方法已应用于结构静强度和稳定性分析。例如，Bradley[10]采用一种等效平板法对翼身融合布局飞机进行有限元静力分析。赵群[11]采用等效刚度法对复合材料加筋板总体失稳性能进行分析。Wenzel[12]采用等效刚度方法开展了以静强度和屈曲作为约束的复合材料垂尾结构重量计算方法的研究。但应用等效刚度有限元模型来分析复材机翼的颤振特性，还鲜见有关资料的报导，动力学特性和颤振特性的预测精度还不甚明了。

为给机翼初步设计阶段提供一种高精度的、快速的结构动力学分析和颤振分析方法，本文应用刚度等效方法简化复合材料机翼结构有限元模型，验证机翼等效刚度有限元模型的计算精度。

1 等效刚度方法

1.1 刚度等效的思路

现代飞机通常采用桁条和蒙皮组成的半硬壳式加筋结构。桁条相互平行或者接近平行，其两端有翼肋支持。桁条有两个主要的作用：一是与蒙皮一起承受轴向拉伸/压缩载荷；二是使加筋壁板的等效中面偏移蒙皮中面，提高了蒙皮局部抗弯刚度，提升了蒙皮局部抗失稳性能。

等效刚度方法的目的就是简化薄壁加筋结构的有限元模型。它的主要思路是：将加强筋离散为一系列的板条元，然后计算加筋壁板的等效中面，根据复合材料经典层合板理论，分别计算蒙皮的刚度系数和各板条元刚度系数，叠加而成等效刚度矩阵；如果是夹芯结构，则是面板刚度系数和夹芯结构刚度系数叠加而成等效刚度矩阵，最后将等效刚度矩阵赋予等效板的属性参数。图1显示了用等效刚度法对不同的加筋结构方案的等效过程。等效刚度方法带来的益处是：当加筋结构方案改变时，只需调整加强筋板条元的外形参数(宽和高)和加强筋布置参数(间距)，无需重新建立有限元模型，从而增强了有限元模型的参数化功能。

图1 各种加筋/夹层壁板等效为光板Fig.1 Modeling approach of different stiffened panels

将刚度矩阵A中的系数定义薄膜材料(Membrane material)和剪切材料(Shear material)，将刚度矩阵B中的系数定义耦合材料(Coupling material)，将刚度矩阵D中的系数定义弯曲材料(Bending material)，然后将这些材料赋予有限元模型中等效板的属性。对于不同布局的加强壁板，计算出不同的刚度矩阵，因此赋予光板不同的力学属性。

1.2 等效刚度参数

(1)

式中Aij，Bij和Dij分别代表层板的面内、耦合以及弯曲刚度系数。

图2 层合板示意图Fig.2 Coordinate locations of plies in a laminate

复合材料蒙皮和加强筋可视为各向异性板，可按照复合材料经典层合板理论计算蒙皮和各板条元的刚度系数，载荷/变形关系如下[13]

(2)

1.3 加筋壁板等效中面计算

计算蒙皮和各板条元的弯曲刚度需要参考中性面的位置。桁条的存在使加强区域的结构中性面偏离蒙皮中面，而桁条之间的非加强区域的中性面仍然与蒙皮中面重合。因此，加筋壁板实际中性面并不是平面而是曲面，如图3所示。

图3 等效中面Fig.3 Equivalent neutral surface

当加筋壁板承受弯矩时，蒙皮和桁条均以中性面产生形变，但由于该实际中性面为曲面，难以确定各板条元以其为参考面的本构关系，所以采用一个等效中性面(图3中xy平面)替代曲面的实际中性面，以便于计算。基于纯弯曲假设，弯曲变形时等效中性面的曲率为кx，整个截面沿筋条方向的合力为零[14]，即

(3)

由于Δhxi=-Δzxc+Δzxi，因此式(3)可写为

(4)

Δzxc表示等效中性面到蒙皮中面的距离，由式(4)可以确定其值为

(5)

1.4 刚度系数平移计算

由于加筋壁板的等效中面偏移蒙皮中面，蒙皮的弯曲刚度得到增强，需要进行平行移面计算。蒙皮和加强筋相对于自身中面刚度系数与等效中面的刚度系数之间有如下关系[15]

(6)

图4 中面平移Fig.4 Translation of neutral surface

1.5 蒙皮刚度矩阵计算

复合材料层合板为各向异性板，设计复合材料机翼蒙皮时，为避免固化时因耦合效应引起翘曲变形，蒙皮采取对称铺层。对于正交各向异性对称层合板，可以忽略拉弯耦合刚度，所以蒙皮的载荷/变形关系可简化为

(7)

对于复合材料加筋板，等效中面偏移出蒙皮中面，这时蒙皮的弯曲刚度由两部分组成，一是蒙皮相对于自身中面弯曲刚度，二是中面偏移增加的弯曲刚度，可由式(6)计算。然而采用对称铺层的蒙皮没有拉弯耦合刚度，故平行移面时它不会对弯曲刚度产生贡献。由于桁条轴向为x向，提高蒙皮在zx平面内的抗弯刚度，蒙皮和加强筋板条元只需在z向进行平移计算。平移后的蒙皮刚度矩阵为

(8)

1.6 桁条刚度矩阵计算

对于复合材料桁条，腹板和上下凸缘可采用不同的铺层，故需要分开处理。以工字型加强筋为例，将其离散为5个板条元，各板条元宽度为wi，厚度为hi，如图5所示。图中蒙皮厚度为tsk，桁条间距为bstr，蒙皮中面到组合截面等效中性面的距离为Δzxc。通过设置板条元的参数，很容易地模拟工程中常见不同截面形状的加强筋。例如，将板条元5的宽度和高度参数设置为0时，加强筋则为J形截面。将板条元4，5的参数同时设为0时，则为T形截面，将板条元2，4的参数设置为0时，则是Z形截面。通过加筋壁板的参数化定义，可以快捷地描述不同加筋板的形状。

图5 开口加筋壁板的典型剖面Fig.5 A typical profile of stiffened panel with open stiffeners

1.6.1x向桁条刚度矩阵

假设桁条的轴向为x向，垂直于桁条向右为y向，垂直蒙皮向下为z向，如图5所示。蒙皮和桁条上的载荷按刚度分配，在单位宽度等效刚度蒙皮上，桁条上的轴向压缩载荷与弯矩表示为各离散板条上载荷的叠加：

(9)

(10)

对于共固化的复合材料蒙皮和桁条，蒙皮接触的板条元可以承受面内剪切载荷，因此需计算该板条元剪切刚度

(11)

式中i=1,2，即只计算板条元1和2的剪切刚度。

同时桁条的扭转刚度也需考虑，对于开口加强筋的扭转刚度计算公式如下[15]

(12)

而对于闭口加强筋的扭转刚度计算公式如下

(13)

式中i=2,3,4,6；Ω1为闭口加筋条壁厚中心线与蒙皮中面所围成的面积；GJ1为闭室Ω1的抗扭刚度。

图6 闭口加筋壁板的典型剖面Fig.6 A typical profile of stiffened panel with closed stiffeners

平行移面计算后，x向桁条的载荷/变形关系为

(14)

1.6.2y向桁条刚度矩阵

同理可得，平移移面计算后的y向桁条的载荷/变形关系为

(15)

而双向加筋布局的桁条刚度矩阵为x向和y向桁条刚度矩阵的叠加。

综上所述，统一的桁条刚度矩阵可写为

(16)

1.7 组装蒙皮桁条刚度矩阵

分别计算蒙皮和离散的加强筋板条元的刚度矩阵后，叠加即可组装加筋壁板的等效刚度矩阵，简化为

(17)

组装后的蒙皮和桁条耦合刚度正负抵消，即等效刚度中的耦合刚度Beq=Bsk+Bst=0。

对于夹芯加强壁板，则是计算上、下面板和夹芯结构对等效中面刚度矩阵，然后叠加而成等效刚度矩阵。

2 算例验证

2.1 机翼模型描述

本文以如图7所示客机概念方案为研究背景，以其机翼作为例子，验证等效刚度有限元模型的结构动力学和颤振计算精度。

图7 某客机概念方案Fig.7 A notional civil jet

该机翼展弦比为10.0，参考面积为95.75 m2，1/4弦线后掠角为24.5°，机翼根部相对厚度0.14。机翼的详细有限元模型和等效刚度有限元模型见图8。两种模型的载荷和边界条件相同。桁条(加强筋)截面为矩形，平行于后梁布置，与复合材料蒙皮共同固化，参见图8(a)所示。详细模型包含了53 697个单元，251 003个自由度。等效模型包含了6 205个单元，26 364个自由度，参见图8(b)所示。可见，等效模型中单元数只有详细模型的单元数量的1/9，模型规模大为降低。

复合材料机翼的蒙皮和桁条铺层方向和比例如表1所示。沿展向划分9个区域，每个区域具有相同的厚度。对于详细有限元模型，由于其包含很多细节信息，有限元建模需要8 h左右，各设计区域复合材料定义和属性的赋予也需要7 h，在Patran的Flightloads模块中定义颤振计算工况需要5 h左右，最后提交Nastran进行颤振计算需要20 min左右。而对于等效有限元模型，本文采用了Patran的PCL编程实现了有限元建模、等效材料计算、属性定义和工况提交、颤振计算等所有步骤自动化进行，总运行时间大概为4 min，颤振计算时间约为0.5 min。

图8 有限元模型对比Fig.8 Comparison of two finite element models

表1 复合材料蒙皮和桁条铺层顺序和比例

2.2 固有频率和振动模态计算

应用MSC.Nastran软件对详细有限元模型和等效刚度有限元模型的固有频率和振动模态进行计算分析，其结果见表2所示。由表2可知，等效刚度模型的各阶振动模态与详细模型的完全一致。两种模型的各阶频率很接近，误差均小于5%。

2.3 颤振分析

采用亚声速偶极子法进行颤振计算，求解方法使用p-k法[16]。在工程分析中，通常取机翼的前5阶模态来做颤振分析。为了验证等效刚度方法的准确性，分别从表2中取详细模型和等效模型的前10阶模态进行颤振计算，计算获得的速度-阻尼(V-g)图和速度-频率(V-ω)图，分别见图9和10所示。从图9和10可以看到，两种模型的速度-阻尼图和速度-频率图的走势一致。

表2 动力学分析结果比较

图9 详细模型颤振分析结果Fig.9 Flutter analysis results of the detailed FE model

图10 等效模型颤振分析结果Fig.10 Flutter analysis results of the equivalent FE model

根据详细模型和等效模型计算获得的颤振速度如表3所示。比较表3中的颤振计算结果，颤振速度误差很小，仅为0.849%，而颤振频率几乎完全吻合。

表3 颤振计算结果比较

3 结 论

为简化复合材料机翼有限元模型，本文研究了一种将加筋壁板结构等效为具有相同力学属性光板的等效刚度方法。通过对比详细有限元模型和等效刚度有限元模型的计算结果，得出以下结论：

1)等效刚度方法大大减小了有限元网格规模和复杂程度，单元数量减至为详细模型的1/9，且有利于机翼结构有限元的参数化建模。

2)两种模型计算出复合材料机翼的固有频率和振动模态一致，最大误差小于5%。

3)两种模型计算出的颤振速度和频率非常一致，颤振速度的误差不到1%。

4)等效刚度有限元模型具有简单快速和准确的特点，可用于复合材料机翼初步设计阶段对颤振特性的评估。

[1] Zhao Yonghui, Hu Haiyan. Structural modeling and aeroelastic analysis of high-aspect ratio composite wings[J]. Chinese Journal of Aeronautics, 2005,18(1):25—30.

[2] Wright J R, Cooper J E. Introduction to Aircraft Aeroelasticity and Loads[M]. Chichester, England: John Wiley & Sons, Ltd, 2007:427—436.

[3] Giles G L. Equivalent plate analysis of aircraft wing box structures with general planform geometry[R]. NASA-TM-87697, 1986.

[4] Livne E, Navarro I. Nonlinear equivalent plate modeling of wing-box structures[J]. Journal of Aircraft, 1999,36(5)：851—865.

[5] 晏飞.复合材料翼面结构综合优化设计技术研究[D].西安:西北工业大学,2001.

Yan Fei. Study on the optimization technologies of composite wing structure[D]. Xi′an: Northwestern Polytechnical University, 2001．

[6] Cesnik C E S, Hodges D H, Patil M J. Aeroelastic analysis of composite wings[R]. AIAA 96-1469-CP. 1996:1 113—1 123.

[7] 李少华，杨智春，谷迎松，等.复合材料连翼的气动弹性剪裁研究[J].西北工业大学学报,2008,26(3):292—296.

Li Shaohua, Yang Zhichun,Gu Yingsong, et al. A different aeroelastic tailoring of a composite joined-wing[J]. Journal of Northwestern Polytechnical University, 2008,26(3):292—296.

[8] 陈全龙,韩景龙，员海玮.考虑发动机推力影响的机翼颤振分析[J].振动工程学报，2012,25(2):110—116.

Chen Quanlong, Han Jinglong, Yuan Haiwei. Flutter analysis of wings subjected to engine thrusts[J]. Journal of Vibration Engineering, 2012,25(2):110—116.

[9] Weck O De, Agte J, Sobieszczanski-Sobieski J, et al. State-of-the art and future trends in multidisciplinary design optimization[R]. AIAA 2007-1905, 2007.

[10]Bradley K R. A sizing methodology for the conceptual design of blended-wing-body transports[D]. George Washington: George Washington University, 2004.

[11]赵群，金海波,丁运亮.加筋板总体失稳分析的等效层合板模型[J].复合材料力学报,2009,26(3):195—201.

Zhao Qun, Jin Haibo, Ding Yunliang. Equivalent laminates model for stiffened panel global buckling analysis[J]. Acta Materiae Compositae Sinica, 2009,26(3):195—201.

[12]Wenzel J, Sinapius M, Gabbert U. Primary structure mass estimation in early phases of aircraft development using the finite element method[J]. CEAS Aeronautical Journal, 2012,3(1):35—44.

[13]Jones R M. Mechanics of Composite Materials[M]. 2nd ed. Philadelphia: Taylor & Francis Inc, 1999:55—118.

[14]Timshenko S, Woinhowsky-Krieger S. Theory of Plates and Shells[M]. 2nd ed. New York: McGraw-Hill Book Company, 1959:36—45.

[15]李顺林.复合材料工作手册[M].北京:航空工业出版社,1988:334—335.

Li Shunlin. An Handbook of Composite Materials[M]. Beijing :Aviation Industry Press, 1988:334—335.

[16]Rodden W P, Harder R L, Bellinger E D. Aeroelastic addition to NASTRAN[R]． NASA CR-3094. Washington: NASA, 1979.

Flutter prediction of composite wing using finite element method with equivalent stiffness

OUYANGXing,YUXiong-qing,WANGYu

(Key Laboratory of Fundamental Science for National Defense-Advanced Design Technology of Flight Vehicle,Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China)

A method using the concept of equivalent stiffness was developed for the modal and flutter analysis of composite wing in preliminary design. The stiffness matrices calculated from different stringer-stiffened panels were applied to FEM shell properties of the same planform panels (i.e., unstiffened panels) that have the same mechanical properties as those of the original stringer-stiffened panels. The advantage of this method is that the complicated definition of composite material and modeling of stiffeners are avoided in the equivalent finite element model (EFEM) and the complicacy of finite element model is reduced, while the feature of structural layout of wing was still retained. A wing of the notional civil jet was used as an example to verify the method. Both the detailed finite element model (DFEM) and the EFEM of wings are created and analyzed for flutter prediction. Especially, the EFEM is automatically created by running scripts, and lots of time is saved from the process of modeling and analysis. The results indicate that the natural frequencies, mode shapes and flutter speed from the two models are in a very good agreement. The accuracy of the EFEM for the modal and flutter analysis of wings is verified. The EFEM provides an attractive method for wing flutter prediction in aircraft preliminary design phase due to its rapid modeling and satisfactory accuracy.

wing; flutter; composite; finite element; equivalent stiffness

2014-03-15；

2014-12-09

V211.41; V215.3+4

A

1004-4523(2015)03-0404-07

10.16385/j.cnki.issn.1004-4523.2015.03.009

欧阳星(1986—)，男，博士研究生。电话：(025)84892102；E-mail：ouyangxing@nuaa.edu.cn

余雄庆(1965—)，男，教授，博士生导师。电话：(025)84892102；E-mail：yxq@nuaa.edu.cn