基于POD降阶方法的复合材料曲壁板颤振响应特性研究

周 建 ， 杨智春

(1.西安航天动力研究所，西安 710100； 2. 西北工业大学 ，西安 710072)

基于POD降阶方法的复合材料曲壁板颤振响应特性研究

周 建1， 杨智春2

(1.西安航天动力研究所，西安 710100； 2. 西北工业大学 ，西安 710072)

建立了三维复合材料曲壁板的气动弹性有限元方程，将本征正交分解方法(POD)应用于三维复合材料曲壁板的非线性颤振响应降阶分析中，通过POD方法构造三维复合材料曲壁板颤振响应的POD模态，然后将系统的运动方程变换到POD模态坐标下，通过数值积分方法计算三维复合材料曲壁板的颤振响应，与传统的模态缩减法计算结果相比，结果很好的吻合，且大大节省了计算时间。

曲壁板；壁板颤振；本征正交分解

自20世纪50年代飞行器速度达到超音速以来，发生了一些因壁板颤振而引起的飞行事故[1-2]。壁板颤振是飞行器表面蒙皮结构由于空气动力、惯性力和弹性力的相互耦合作用而产生的一种气动弹性不稳定现象，是动力学系统的一种自激振动。由于几何非线性的影响，在壁板发生颤振时并不会瞬间发生破坏，而是表现为固定幅值的极限环振动。

在壁板的非线性颤振响应研究中，由于壁板颤振系统的维数较大，很多研究者基于模态叠加的思路来减缩壁板颤振系统的维数[3-4]。ZHOU等将此种模型技术应用在受热壁板的非线性颤振响应分析中，并与采用迦辽金方法的分析结果进行了对比，说明了该方法的正确性。研究表明，为了得到具有足够精度的壁板颤振响应的近似解，对于三维壁板，至少需要保留36阶固有模态。用如此多的模态阶数计算壁板的非线性颤振响应，不仅时间开销巨大，而且对于颤振抑制器的设计带来了非常大的困难[5-6]，为了解决这些问题，则需要寻找新的降阶方法，以尽可能低阶的模型来描述壁板非线性颤振系统。在2002年，GUO和MEI[7-8]提出了气动弹性模态的概念，并利用壁板的气动弹性模态进行系统降阶，分析表明，这种基于气动弹性模态的模型降阶技术能大大提高壁板非线性颤振响应的分析效率，但由于气动弹性模态不具有正交性，在气动弹性模态坐标系下，各阶模态之间存在质量耦合和刚度耦合；另外为了保证气动弹性模态为实模态和提高求解精度，计算气动弹性模态时的气流动压应该接近颤振临界动压。近些年来，本征正交分解(Proper Orthogonal Decomposition，POD)方法在动力学系统降维中得到了广泛的应用。POD技术主要依赖于已有的数据结果[9-10]，因此很难对其进行物理解释，理解的关键在于由POD技术降阶得到的模态(POMs)是相互正交的，而固有模态只是对质量和刚度正交。应用POD方法对系统进行降维有很多优点：①POMs是从实验数据或数值计算数据中计算提取，不需要计算结构矩阵；②可以直接利用特征值来判断并提取所需要的模态数，简单方便；③POMs是最优的正交基，用其降维效率高和效果好。

在航空领域，许多研究人员采用POD降阶方法来计算翼型表面流场[11-15]，减少了计算时间，提高了计算效率。也有一些研究者将本征正交分解方法应用于壁板非线性颤振响应的降阶分析中[16-18]，但他们对于此降阶方法的应用仅仅停留在二维平壁板非线性颤振的基础上，结构模型相对简单，对于工程应用还存在一定的差距。

本文建立了三维复合材料曲壁板的非线性颤振运动方程，将POD降阶建模方法应用于三维复合材料曲壁板的非线性颤振响应特性研究中，并与传统的模态缩减法计算结果进行了对比。

1 三维复合材料曲壁板颤振系统有限元模型

采用TESSLER和HUGHES[19]提出的MIN3三角形三结点Mindlin板单元来建立三维复合材料曲壁板颤振有限元模型，MIN3单元每个结点有5个自由度，包含三个位移自由度(面内位移ui和vi，横向位移wi)和两个转动自由度(绕x轴θxi和绕y轴θyi)。

则MIN3单元的结点位移向量可以表示为：

{w}={wbθwm}T

(1)

其中：

{wb}=[w1w2w3]T

{θ}=[θx1θx2θx3θy1θy2θy3]T

{wm}=[u1u2u3v1v2v3]T

单元的位移场可以通过结点位移插值求得：

(2)

图1示出了三维曲壁板的示意图，a为x方向曲壁板的跨度，b为y方向曲壁板的长度，H为曲壁板的最大拱高，h为曲壁板的厚度，V∞为来流速度。图2示出了曲壁板上任意一点位移的定义，其中，u0、v0和w0为曲壁板上任意一点的初始几何位置坐标；u、v和w为曲壁板上任意一点沿三个位移自由度方向的位移改变量。

图1 三维曲壁板示意图Fig.1 The curved panel geometry

图2 曲壁板上任意一点位移的定义Fig.2 Displacement Definition of curved panel

三维曲壁板上任意一点处的总应变-位移关系可表示为：

(3)

曲壁板的横向剪切应变可以表示为：

(4)

壁板的应力-应变关系可以表示为：

(5)

式中:{σ}为壁板面内力，{τ}为壁板的面内剪切应力，[Q]和[Qs]分别为弹性模量矩阵和剪切模量矩阵。

(6)

式中：

从虚功原理出发，单元的内力虚功表达式：

(7)

式中:S为单元的面积，αs为横向剪切校正因子[19]。

带有曲率修正的一阶活塞理论气动力可以表示为[20]：

(8)

定义无量纲动压和无量纲气动阻尼：

(9)

则

(10)

忽略壁板面内方向和挠度方向之间的惯性耦合，则单元的外力虚功为：

(11)

由虚功平衡方程：

δWint-δWext=0

(12)

组装单元矩阵，得到三维复合材料曲壁板的颤振运动方程：

(13)

式中:[M]为质量矩阵，[C]为阻尼矩阵，[K0]为小挠度变形下的弹性刚度矩阵，λ[Ax]为气动刚度矩阵， [K]s为剪切刚度矩阵，[K]θ0为初始几何曲率引起的刚度修正矩阵，[N1]0，[N1]θ0，[N1]Nb，[N1]Nm，[N1]Nθ0为一阶非线性矩阵，与结点位移的一次项有关，[N2]为二阶非线性矩阵，与结点位移的二次项有关，λ{Pw0,x}为气动静载荷。

令：

[Aa]=λ[Ax]

[KL]=[K0]+αs[K]s+[K]θ0

[N1]=[N1]0+[N1]θ0+[N1]Nb+

[N1]Nm+[N1]Nθ0

{P}=-λcosΛ{Pw0,x}

则式(13)可以表示为：

2 POD降阶建模方法

POD降阶方法的基本思想是利用n维空间的数据样本(称为“快照”)提取m(m≪n)维最优基，构建一个最优子空间，形成降阶系统，使降阶系统与全阶系统的误差达到最小。在POD方法中首先要得到“快照”矩阵S，“快照”矩阵的数据一般采用实验测试或数值模拟的方法获得，本文的快照矩阵是由36阶固有模态截断法计算得到的曲壁板非线性颤振响应数据(通常认为是精确解)构造而成：

(15)

式中:W(xn,yn,tm)是三维复合材料曲壁板有限元中编号为n的结点在tm时刻的响应，m为选取时间点的个数。

通过“快照”矩阵构造关联矩阵：

P=STS

(16)

求解得到关联矩阵P的特征值λi和特征向量vi(i=1,2,…,m)，λi是对应vi包含能量大小的量度，用于维数估计。将特征值进行降序排列λ1≥λ2≥…≥λm，并将对应的特征向量进行重新排列。通过以下标准对POD模态进行截断：

(17)

式中:ε为根据精度要求设置的参数标准；l为符合标准的POD模态的数目。将前l维特征向量vi组成矩阵V=[v1,v2,…,vl]，将这些向量进行线性叠加Ψ=SV，即可提取POD模态Ψi(i=1,2,…,l)。对三维复合材料曲壁板非线性颤振系统在物理空间中的方程进行POD模态坐标变换，即有：

W(x,y,t)=Ψ·a(t)

(18)

式中:a=[a1a2…al]T为POD降阶模态下的坐标向量。

将式(18)代入三维复合材料曲壁板颤振运动方程(14)，引入下列POD模态坐标变换式：

其中，[N1]i是三维复合材料曲壁板对应于第i阶POD模态Ψi的一阶非线性刚度矩阵；[N2]i,j是三维复合材料曲壁板对应于第i阶POD模态Ψi和第j阶POD模态Ψj组合时的二阶非线性刚度矩阵。

于是，POD模态坐标系下三维复合材料曲壁板的颤振运动方程为：

(19)

式中:

通过POD模态坐标变换，系统的自由度数目由全部结点自由度数目之和缩减为POD模态截断后保留模态数之和，于是三维复合材料曲壁板颤振运动方程的维数大大降低，颤振运动方程(19)的状态空间方程形式可写为：

(20)

式中:I为l×l阶单位矩阵。

下文的颤振响应分析中，采用四阶龙格-库塔方法对式(20)进行数值积分，积分时间步长为Δτ=0.000 01，且均是提取复合材料曲壁板中心位置结点的位移响应来表征曲壁板的非线性颤振响应特性。

3 算例研究

本文采用两个三维复合材料曲壁板模型来说明POD降阶方法的有效性，一种是短玻璃纤维曲壁板，可视作准各向同性材料曲壁板，边界条件为四边固支，其几何尺寸为0.38 m×0.305 m×0.002 m，最大拱高为H=0.002 m；另一种是石墨/环氧铺层曲壁板，铺层方式为[-40/40/-40]，边界条件为四边简支，其几何尺寸为0.381 m×0.305 m×0.001 22 m，最大拱高为H=0.001 22 m。石墨/环氧材料及玻璃纤维材料的力学性能参数分别见表1和表2。

表1 石墨/环氧材料的力学性能参数

表2 玻璃纤维材料的力学性能参数

图3示出了动压λ=1 100时，采用36阶固有模态降阶模型计算得到的准各向同性曲壁板中点处非线性颤振响应的位移时间历程和相图，颤振响应的位移峰值和谷值分别为-0.203 6、-0.674 2，这是因为曲壁板的曲率引起的静气动载荷使得曲壁板发生静气动弹性变形，颤振响应以静气动弹性变形位置为平衡点的振动。图4示出了采用4阶固有模态降阶模型计算得到的准各向同性曲壁板中点处非线性颤振响应的位移时间历程和相图，尽管曲壁板的非线性颤振响应仍然表现为极限环运动，但是其颤振响应的峰、谷值为0.450 8、-1.152，与采用36阶固有模态降阶模型计算得到的结果相差很大；而能够接近采用36阶固有模态降阶模型计算得到的非线性颤振响应，至少需要截取20阶固有模态，图5示出了采用20阶固有模态降阶模型计算得到的准各向同性曲壁板中点处非线性颤振响应的位移时间历程和相图，其峰值与谷值分别为-0.202 8、-0.669 7，计算时间为12 549 s。从利用36阶固有模态降阶模型计算的非线性颤振响应中提取“快照”样本，通过POD方法计算曲壁板的POD降阶模态，其对应的前8阶POD降阶模型特征值如表3所示，由表中可以看出，低阶特征值占得能量比例较大，前4阶POD模态占总能量的99.99%，因此，选取前4阶POD模态计算准各向同性曲壁板中点处的非线性颤振响应，计算时间仅为419.316 s，如图6所示，从图中可以看出，通过4阶POD模态降阶模型计算得到的曲壁板极限环颤振响应与36阶固有模态降阶模型计算得到的结果相比较，无论在定性方面还是定量方面，都吻合的很好，其颤振响应的峰、谷值为-0.203 9、-0.671 2；与至少选取20阶固有模态降阶模型相比，显然通过POD降阶的方法大大降低了计算曲壁板颤振响应的时间。

图3 36阶固有模态降阶模型计算得到的位移时间历程和相图Fig.3 Time history and phase plot using 36 normal modes

图4 4阶固有模态降阶模型计算得到的位移时间历程和相图Fig.4 Time history and phase plot using 4 normal modes

图5 20阶固有模态降阶模型计算得到的位移时间历程和相图Fig.5 Time history and phase plot using 20 normal modes

λ1λ2λ3λ4λ5λ6λ7λ811.14.80.150.010.0012.6E-0041.8E-0049.2E-005

图6 4阶POD模态降阶模型计算得到的位移时间历程和相图Fig.6 Time history and phase plot using 4 POD modes

图7 36阶固有模态降阶模型计算得到的位移时间历程和相图Fig.7 Time history and phase plot using 36 normal modes

图7示出了动压为λ=540情况下，采用36阶固有模态降阶模型计算得到三维复合材料曲壁板中点处非线性颤振响应的位移时间历程和相图，颤振响应的位移峰、谷值为-0.473 6，-0.669 5；图8示出了用6阶固有模态降阶模型计算得到的复合材料曲壁板中点处非线性颤振响应的位移时间历程和相图，通过比较可以发现，在采用6阶固有模态降阶模型计算时，复合材料曲壁板的颤振系统响应收敛到静态平衡点，并不是极限环运动，因此用6阶固有模态降阶模型计算复合材料曲壁板的非线性颤振特性是不可行的，而能够接近采用36阶固有模态降阶模型计算得到的非线性颤振响应，至少要截取20阶固有模态，图9示出了采用20阶固有模态降阶模型计算得到的复合材料曲壁板中点处非线性颤振响应的位移时间历程和相图，其峰值与谷值分别为-0.470 5、-0.670 4，计算时间为11 319 s。以36阶固有模态降阶模型得到的非线性颤振响应的时间历程作为“快照”样本，通过POD方法计算曲壁板的POD降阶模态，其对应的前8阶POD模态降阶模型特征值如表4所示，可见其前6阶POD模态占总能量的99.99%，故选取前6阶POD模态对曲壁板颤振模型进行降阶，计算复合材料曲壁板中点处的非线性颤振响应的位移时间历程及相图，计算时间为735.931s，如图10所示，其位移响应的峰、谷值为-0.473，-0.671 3，可以看出，用6阶POD模态降阶模型很好的模拟复合材料曲壁板的颤振极限环响应特性，与至少选取20阶固有模态降阶模型相比大大降低了计算曲壁板颤振响应的时间。

图8 6阶固有模态降阶模型计算得到的位移时间历程和相图Fig.8 Time history and phase plot using 6 normal modes

图9 20阶固有模态降阶模型计算得到的位移时间历程和相图Fig.9 Time history and phase plot using 20 normal modes

λ1λ2λ3λ4λ5λ6λ7λ89.90.60.028.6E-0043.5E-0041.2E-0042.1E-0051.4E-005

图10 6阶POD模态降阶模型计算得到的位移时间历程和相图Fig.10 Time history and phase plot using 6 POD modes

现在的问题是：由某一动压下三维复合材料曲壁板颤振响应得到的POD降阶模型是否适用于其它动压下曲壁板颤振响应分析。图11示出了分别POD降阶模型与36阶固有模态降阶模型计算曲壁板颤振响应计算结果比较图。图11(a)是采用动压为λ=1040时准各向同性曲壁板非线性颤振响应数据计算得到的POD模态进行降阶，得到POD模态降阶模型，再用该模型计算准各向同性曲壁板在其它动压下的非线性颤振响应幅值，与用36阶固有模态降阶模型计算的非线性颤振响应幅值进行比较；图11(b)是采用动压为λ=540时复合材料曲壁板非线性颤振响应数据计算得到的POD模态降阶模型，再用该模型计算复合材料曲壁板在其它动压下的非线性颤振响应幅值，与用36阶固有模态降阶模型计算的非线性颤振响应幅值进行比较。从比较的结果可以看出，通过POD模态降阶模型计算的非线性颤振响应与采用36阶固有模态降阶模型计算的结果基本一致。说明，由某一动压下三维复合材料曲壁板颤振响应得到的POD降阶模型适用于其它动压下曲壁板颤振响应分析。

图11 用POD降阶模型与36阶固有模态降阶模型计算的响应比较Fig.11 Comparison using 36 normal modes and POD modes

4 结 论

本文通过有限元方法建立了三维复合材料曲壁板的非线性颤振运动方程，并将本征正交分解(POD)降阶方法应用到了三维复合材料曲壁板非线性颤振响应分析中，通过与传统的模态缩减法计算结果相比较，指出无论定性还是定量方面，其结果都很好的一致，且采用POD降阶建模方法可以有效的降低了系统的维数并大大节省了计算时间，提高了计算分析效率。

[ 1 ] GARRICK I E, REED Ⅲ W H. Historical development of aircraft flutter[J]. Journal of Aircraft, 1981,18(11): 897-912.

[ 2 ] ONG C C. Flutter of a heated shield panel [R]. NASA-CR-122855, BELLCOMM Inc. TM-71-1013-6, 1971.

[ 3 ] ZHOU R C, XUE D Y, MEI C. Finite element time domain-modal formulation for nonlinear flutter of composite panels[J]. AIAA Journal, 1994,32(10): 2044-2052.

[ 4 ] ABDEL-MOTAGALY K, CHEN R, MEI C. Nonlinear flutter of composite panels under yawed supersonic flow using finite elements[J]. AIAA Journal, 1999,37(9):1025-1032.

[ 5 ] ABDEL-MOTAGALY K, DUAN B, MEI C. Active control of nonlinear panel flutter under yawed supersonic flow[J]. AIAA Journal, 2005,43(3):671-680.

[ 6 ] ZHOU R C, MEI C, HUANG J K. Suppression of nonlinear panel flutter at supersonic speeds and elevated temperatures[J]. AIAA Journal, 1996,34(2):347-354.

[ 7 ] GUO X Y, MEI C. Application of aeroelastic modes on nonlinear supersonic panel flutter at elevated temperatures[J]. Computers and Structures, 2006(84): 1619-1628.

[ 8 ] GUO X Y, MEI C. Using aeroelastic modes for nonlinear panel flutter at arbitrary supersonic yawed angle[J]. AIAA Journal, 2003,41(2):272-279.

[ 9 ] FEENY B F, KAPPAGANTU R. On the physical interpretation of proper orthogonal modes in vibrations[J]. Journal of Sound and Vibration, 1998,211(4):607-616.

[10] 于海, 陈予恕. 高维非线性动力学系统降维方法的若干进展[J]. 力学进展, 2009,39(2): 154-164. YU Hai, CHEN Yushu. Recent developments in dimension reduction methods for High-Dimension dynamical systems[J].Advances in Mechanics, 2009, 39(2): 154-164.

[11] 赵松原, 黄明恪. 模拟退火算法和POD降阶模态计算在翼型反设计中的应用[J]. 空气动力学报, 2007,25(2): 236-240. ZHAO Songyuan, HUANG Mingke. Application of simulated annealing method and reduced order models based on POD to airfoil inverse design problems[J]. Acta Aerodynamic Sinica,2007,25(2):236-240.

[12] 赵松原, 黄明恪. POD降阶算法中对基模态表达的改进[J]. 南京航空航天大学学报, 2006,38(2): 131-135. ZHAO Songyuan, HUANG Mingke. Modification to basic modes of reduced order model for computing airfoil flow field[J]. Journal of Nanjing Univesity of Aeronautics & Astronautics, 2006,38(2):131-135.

[13] DAN X, MIN X, DOWELL E H. Proper orthogonal decomposition reduced-order model for nonlinear aeroelastic oscillations[J]. AIAA Journal, 2014,52(2):229-241.

[14] ATTAT P J, DOWELL E H, WHITE J R. Reduced order nonlinear system identification methodology[J]. AIAA Journal, 2006,44(8):1895-1904.

[15] LEGRESLEY P A, ALONSO J J. Airfoil design optimization using reduced order models based on proper orthogonal decomposition[R]. AIAA-2000-2545,2000.

[16] MORTARA S A, SLATER J, BERAN P. Analysis of nonlinear aeroelastic panel response using proper orthogonal decomposition[J]. Journal of Vibration and Acoustics, 2004,126(3):416-421.

[17] EPUREANU B I, TANG L S, PAIDOUSSIS M P. Coherent structures and their influence on the dynamics of aeroelastic panels[J]. International Journal of Non-Linear Mechanics, 2004(39):977-991.

[18] LUCIA D J, BERAN P S, KING P I. Reduced order modeling of an elastic panel in transonic flow[R]. AIAA-2002-1594,2002.

[19] TESSLER A, HUGHES T J R. A three-node mindlin plate element with improved transverse shear[J]. Computer Methods in Applied Mechanics and Engineering, 1985,50(1):71-101.

[20] AZZOUZ M S, PRZEKOP A, GUO X Y, et al. Nonlinear flutter of shallow shell under yawed supersonic flow using FEM[R]. AIAA-2003-1516,2003.

Flutter response characteristics of composite curved panels based on POD method

ZHOU Jian1, YANG Zhichun2

(1. Xi’an Aerospace Propulsion Institute, Xi’an 710100, China； 2. Northwestern Polytechnical University, Xi’an 710072, China)

The equations of motion for nonlinear flutter of curved composite panels were developed with the finite element method. The reduced order modes constructed with the proper orthogonal decomposition (POD) method were used in reducing the order of these equations, and the equations of motion were transformed into a reduced nonlinear system under the POD modal coordinates, then the reduced equations were solved in time domain by using the numerical integration method. Compared with the results calculated using the traditional modal reduction method, the results using POD method based on reduced order models agreed well with the former, and also saved the computation time greatly.

curved panel; panel flutter; proper orthogonal decomposition

国家自然科学基金(11072198)

2015-09-16 修改稿收到日期：2016-12-27

周建 男，博士，工程师，1987年3月生

杨智春 男，教授，博士生导师，1964年2月生 E-mail：yangzc@nwpu.edu.cn

V215.3

A

10.13465/j.cnki.jvs.2017.01.006