二维功能梯度壁板热颤振本征问题的精确解

代林桐，邢誉峰

(北京航空航天大学 固体力学所，北京 100083）

颤振是指弹性结构在气流中由于受到气动力、弹性力和惯性力的耦合作用而发生的振幅不衰减的自激振动，对结构的使用寿命影响很大，甚至危及结构的安全［1］。

随着飞行马赫数的不断提高，飞行器面临的气动环境越来越严峻。热气动弹性是近年来最热门的方向之一，研究成果不断涌现。作为有着优良的力学性能和耐热特性的功能梯度材料(FGM）板，研究其气弹及热气动弹性问题具有重要的理论和应用价值。这也给气动弹性问题的研究带来了新的挑战，那就是在分析中必须考虑高马赫数下气动加热所引起的热效应对飞行器结构气动弹性的影响［2］。

Praveen和Reddy［3］对功能梯度壁板进行了非线性瞬态热弹性分析，并讨论了温度场对壁板响应的影响。苑凯华和邱志平［4］利用有限元法建立了壁板在热效应下的运动微分方程，基于非线性动力学模型研究了对壁板施加控制后对颤振极限环幅值的影响。Hosseini等［5］对超声速气流中功能梯度平板的热气弹问题进行了深入分析，建立了基于冯·卡门薄板大变形理论的结构方程并考虑了温度变化对壁板材料性能的影响，通过一阶活塞理论模拟气动力，并采用Galerkin方法对壁板的控制方程进行求解。李丽丽和赵永辉［6］研究了热环境下四边固支壁板结构频率特性的变化，进而利用p-k法进行了颤振分析，研究表明，热效应对结构的频率特性有很大的影响，并影响颤振边界。Shahverdi和Khalafi［7］利用广义微分求积的方法研究了功能梯度曲板的热气动弹性行为。黄小林等［8］基于复合材料薄板理论和有气流偏角的气动压力的一阶活塞模型，用Galerkin方法分析了气流偏角、热环境等因素对FGM板固有频率和颤振临界速压的影响。

壁板颤振分析最终归结为利用不同的结构和气动力模型求解颤振微分方程的问题，其解法包括数值解法和解析解法。

已有文献主要采用Galerkin方法［9］及Rayleigh-Ritz方法［10］分析壁板的非线性颤振特性。虽然这2种方法精度比较高，但仅能从数值解的角度得到颤振特性。而解析解法，如直接求解法［11］、半逆法［12］，这几种方法都能直接从求解本征方程入手来研究颤振特性。

基于活塞理论，Li和Song［11］采用直接求解法求得了Kirchhoff和Mindlin板在不同边界下的颤振边界，但没有考虑热效应的影响。Sun和Xing［13］根据经典板理论和一阶活塞理论建立了二维层合板在简支、固支和自由边界条件下的气动弹性模型，得到了壁板颤振本征解的统一显式形式，并与Galerkin方法进行了比较，但也未考虑超声速下气动热的影响。

目前，鲜见到公开发表的关于功能梯度壁板热颤振本征问题精确解的工作。针对该问题，本文基于经典薄板理论和一阶活塞理论求得了不同边界下壁板热颤振本征问题的精确解并分析了其颤振特性。

1 颤振微分方程

考虑如图1所示的二维FGM壁板模型，其弦长为a，展长为无限长，厚为h，NTx相当于热应力的作用，T为温度。坐标系建立在板的中面，原点在板的中心处。壁板的上表面为陶瓷层，下表面为金属层。

图1 两端简支壁板几何模型Fig.1 Geometric model of simply supported panel at both ends

FGM板的弹性模量E、热膨胀系数α、热传导系数K及密度ρ等物理属性量皆按照式(1）进行估计［14］：

式中：P(z）为板内任意一点的材料参数(弹性模量、泊松比、剪切模量等）；Pt和Pb分别为上表面(用t表示）和下表面(用b表示）对应的材料参数；n为材料梯度指数。

功能梯度材料通常工作在一些高温环境，其不同属性对温度的敏感程度是不一样的。设弹性模量E和热膨胀系数α随温度的变化规律为［15］

式中：Pi(i＝－1，0，1，2，3）为材料某一物理属性的温度相关系数，其值是唯一的。

由式(2）可知，要获得有效物理属性，除了知道其对应的温度系数以外，还需要知道结构温度场分布。当材料热传导系数为常数时，由一维热传导稳态方程所得到的是一个线性的温度场。对于功能梯度材料组成的板而言，当热传导系数k(i）是厚度坐标的函数而不再是一个常数时，这时得到一个非线性的温度场。该非线性温度场满足如下一维热传导方程：

相关文献给出了其近似理论解［16］。

假定壁板上表面受气动力作用，另一侧为静止空气。气体流速为V，密度为ρa，马赫数为Ma。对于本文后续选择的Si3N4和SUS304作为组成成分的FGM板，泊松比υ沿着厚度方向的变化对结果几乎没有影响，因此本文在后续计算中取泊松比为2种材料的泊松比的均值υ＝0.28［17］，称为等效泊松比。

下面基于经典薄板理论和Ham ilton变分原理，推导二维壁板颤振控制微分方程。在薄板理论中，位移分量可以表示为

式中：u0为FGM 矩形薄板几何中面x方向的位移；w为FGM薄板的挠度。

在二维壁板颤振问题分析中，展长假设为无限长，因此不考虑y方向位移的作用。在下面公式推导中，采用的是单位宽度。

FGM薄板的应变可以表示为

应力为

式中：εx(T）为温度应变，即

热环境中FGM矩形薄板的总应变能为

FGM薄板自由振动的动能为

对于功能梯度板来说，由于材料不均匀，材料中性面与几何中面并不重合，设二者之间距离为z0。在中性面上位移u＝0，从式(4）可得

由式(6）可得

轴力为

本文等效泊松比为常数，因此利用轴力等于0得到

实际上，即使泊松比不是常数，中性面也是存在的，这不同于三维壁板情况［18］。

根据一阶活塞理论［19］，气动力的形式为

式中：q为动压。

气动力所做虚功δW的积分为

根据Ham ilton原理的广义形式为

可以得到二维功能梯度壁板在超声速流下的运动微分方程，如下：

其中：Deq为等效刚度；为热应力。

本文用到的简支和固支边界条件如表1所示。

表1 简支和固支边界条件Table 1 Boundary conditions for sim p le support and clam p

下面将根据边界条件求方程(17）的本征精确解，并分析梯度指数、温度场对颤振边界的影响。

2 颤振精确本征解和颤振机理

2.1 精确本征解

设挠度w有如下分离变量的形式：

式中：颤振的稳定性取决于本征值Ω，其实部β代表幅值，而虚部ω代表振动频率。若β＝0，壁板将发生颤振。

将式(19）代入控制微分方程(17）可得

式(20）即是二维壁板颤振的本征微分方程。根据边界条件求解方程(20）可以得到颤振问题的模态函数φ和频率方程。引入无量纲量ξ＝x／a。为了求得本征微分方程的通解，φ可写为

式中：A为待定系数；λ为x方向的空间本征根。

将式(21）代入方程(20）可得

或写成

这即是二维壁板颤振本征代数方程，其中R、p和k都是无量纲系数，具体形式为

求解方程(23）可以得到二维壁板的空间本征根，为一对实数根和一对复数根：

则本征函数φ的通解可表示为

将本征根(25）代入本征代数方程(23），得到如下重要关系式：

当给定壁板和气流参数时，便可确定R与p，因此α1、β1和ϑ中只有一个独立变量，这里选ϑ为独立变量。同时k也取决于ϑ，根据式(24）可知Ω取决于k。将式(19）与式(26）代入表1所示的边界条件，可求得频率方程与本征函数的系数(见表2），进而求出时间本征值Ω，其反应了壁板的振动特性。

表2 二维壁板颤振频率方程和本征函数Table 2 Eigensolutions of two-dim ensional panel flutters

2.2 颤振机理

颤振的发生必然是因为壁板振动过程中有气动力的作用。下面通过分析颤振方程和本征根Ω的性质来判断壁板为何会发生颤振。

3 算例与影响因素参数分析

首先，针对固有频率和颤振参数，把精确解与有限元分析结果和Galerkin方法结果进行比较，以验证精确解的正确性，再对影响颤振特性的参数进行分析。

3.1 典型算例的计算结果比较

选择如表3中所示陶瓷和金属材料组成的梯度材料。只考虑弹性模量E、热膨胀系数α与温度相关，热传导系数为Kc＝9.19 W／(m·K）和Km＝12.04 W／(m·K），其中c表示陶瓷材料，m表示金属材料。考虑非线性温度场，上表面温度350 K，下表面温度为300 K，梯度指数n为5，弦长a为0.5 m，跨厚比a／h＝250。

表3 功能梯度材料弹性常数Table 3 Elastic constants of FGM

1）与有限元结果的比较

现有有限元商业软件如ANSYS、ABAQUS等都未提供成熟的功能梯度材料分析模块。功能梯度板最主要的结构特点在于其材料特性在板厚度方向不断变化，因而在有限元模型中，将功能梯度板认为成层合板。层合板的每层均为各向同性材料，每层的材料参数由该层所处的位置确定。考虑了总层数分别为5层和10层2种情况。表4及表5给出了根据式(1）计算得到的各分层密度、弹性模量。

表4 分层为5层时密度和弹性模量Table 4 Density and m odulus of elasticity for 5 layers

表5 分层为10层时密度和弹性模量Table 5 Density and m odulus of elasticity for 10 layers

表6将所得功能梯度板的前两阶固有频率与ABAQUS软件计算所得结果进行了对比。有限元网格的划分为10×2 500，且采用平面应力单元。可以看出，精确解与ABAQUS结果存在差异，当功能梯度板的分层增加时，差异在减小。当分层数为10层时，SC边界情况的第2阶频率之间的差异已经小于0.413%。由此可以推断，建立的FGM有限元模型是可行的，所得精确固有频率是正确的。

表6 不同边界条件下FGM 板频率Table 6 Frequency of FGM p late under different boundary conditions

2）与Galerkin方法结果的比较

在颤振理论分析领域，Galerkin方法的应用是非常广泛的。下面把精确解与Galerkin结果进行比较，以验证所得颤振结果的正确性。考虑一均匀壁板，其刚度为D＝148.73 N·m2，密度为ρ＝1 600 kg／m3，弦长a为0.3 m，跨厚比a／h＝200。不考虑温度场作用，边界条件为两端简支。Galerkin方法选用前三阶简支模态。表7给出了一、二阶固有频率ω1、ω2，以及颤振时对应的颤振频率ωf及临界马赫数Maf。可以看出，二者结果吻合较好，验证了精确颤振频率和马赫数的正确性。如果Galerkin方法选用更多阶的模态，则二者结果将更加接近。

表7 本文解与Galerkin方法结果的对比Table 7 Com parison of result between present method and Galerkin method

3.2 参数分析

下面讨论梯度指数、温度场、边界条件等对FGM板的临界动压和频率的影响。

3.2.1 梯度指数n对临界动压的影响

功能梯度材料属于非均匀材料，其组分是随坐标变化的。由式(1）可知，梯度指数n决定了2种材料成分各占多少。对于上表面是陶瓷，下表面是金属的壁板，当n为0时，壁板为各向同性的纯陶瓷壁板。当n为无穷大时表示纯金属壁板。采用表3所示材料，考虑两边简支边界条件，a／h＝250，在不考虑温度场的情况下研究无量纲临界动压λ*随梯度指数n的变化关系，λ*＝ρaV2a3／(MaD0），D0＝Emh3／［12(1－υ2）］，Em表示金属材料在常温T＝300 K下的弹性模量。如图2所示，临界动压值随梯度指数的增加出现先下降较快，而后下降变缓的现象，当指数较大时，临界动压值趋于稳定。

图2 临界动压值随梯度指数变化曲线Fig.2 Critical dynamic pressure versus gradient index

3.2.2 温度对临界颤振频率的影响

图3 临界颤振频率随温度的变化曲线Fig.3 Critical flutter frequency versus temperature

当ΔT＝100 K、n＝5、a／h＝250时，在两端边界为简支的条件下，图4给出FGM板在均匀温度场及非线性温度场下的一阶频率与二阶频率随着马赫数的变化情况。可以得出，2种温度场下均发生频率耦合型颤振，在均匀温度场下，壁板的临界颤振频率为408.24 Hz，临界马赫数为2.90。而在非线性温度场下，壁板的临界颤振频率为423.83 Hz，较均匀温度场时提高了3.82%，临界马赫数为3.02，较均匀温度场时提高了4.14%。

图4 两种热环境下频率随马赫数的变化曲线Fig.4 Frequency versus Mach number in two thermal environments

3.2.3 边界条件对临界动压和颤振频率的影响

表8为简支、固支及其组合边界条件下的功能梯度板的临界颤振频率和临界动压的比较。此处分析了均匀温度场及非线性温度场下，3种边界条件下的无量纲临界颤振频率及临界动压的大小，在相同条件下，三者的λ*与ω*的关系为：两端固支大于一端简支一端固支，后者又大于两端简支。

表8 不同边界下FGM 板的临界颤振频率和临界动压的比较Table 8 Com parison of flutter frequency and critical dynam icpressure of FGM plate under different boundary conditions

4 结论

基于经典薄板理论和一阶活塞理论建立了气动力作用下二维功能梯度壁板的热颤振模型，采用分离变量法得到了颤振问题的精确解。结论如下：

1）系统刚度由结构弹性刚度和气动刚度组成。从数学上而言，颤振现象的发生是由于气动刚度致使系统本征根变成了复数，其实部决定系统振动是衰减、等幅还是发散。

2）功能梯度材料能够有效提高热环境下壁板的颤振边界。

3）在相同条件下，3种边界条件下的临界颤振频率及临界动压大小关系为：两端固支大于一端简支一端固支，后者又大于两端简支。