1. 引言
圆柱形扁壳作为超声速/高超声速飞行器的一种结构形式，其颤振特性一直是气动弹性研究领域不可忽视的重要问题。对圆柱形扁壳/大开口圆柱壳高超声速气动弹性颤振问题的研究起步较晚，相关文献并不多。最早以圆柱形扁壳为对象的气动弹性研究始于Anderson [1] [2] ，他用Galerkin方法预测了一圆柱壳块的颤振边界。Ganapathi [3] 采用双曲四面体壳单元研究了层合双曲壁板的颤振临界动压及铺层方式的影响。Algazin [4] 用有限差分法求解了任意气流偏角的圆柱形扁壳和扁球壳的颤振临界动压及模态。Wang [5] 用一种样条有限模态法及子结构法对棱柱形复合材料板/壳的颤振进行了研究，并对比他人结果，验证了该方法的计算精度和有效性。Oh [6] 用基于多场分层理论的非线性单元研究了压电圆柱形扁壳的颤振临界动压，考虑了气动力、热和压电效应。Shin [7] 用有限元方法研究了杂交复合材料粘弹性圆柱形扁壳的颤振问题，对比了不同铺层方式对结果的影响。Singha [8] 采用16节点的退化等参壳单元研究了曲率、铺层方式、气流方向以及边界条件对圆柱形扁壳颤振临界动压的影响，并考虑了面内压力及剪应力等因素。
由于圆柱形扁壳结构的复杂性和多样性，其在轴向流中发生颤振的临界动压、周向波数、颤振频率对曲率等参数都非常敏感，目前对这一问题尚缺乏较为系统地研究。为深入了解圆柱形扁壳曲率由小到大变化时，系统的气动弹性特性，在基于活塞理论的气动力作用下，通过建立圆柱形扁壳气动弹性运动方程，采用二维的微分求积法离散，运用特征值方法，求解了不同曲率参数对线性系统的颤振临界动压的影响。
2. 圆柱形扁壳的气动弹性运动方程
基于线性理论，在某一轴向临界流速下，壳体由于流体的动压作用会产生失稳。适用于超声速/高超声速颤振分析的气动力理论主要有 [9] 线性活塞理论(Linear Piston Theory)，准定常理论(Quasisteady Theory)，细长体理论(Slender-body Theory)，势流理论(Exact Potential Theory)。由于活塞理论 [10] 的气动力可以表述为较为简单的解析形式，因此为绝大多数的研究者所采用。一阶线性活塞理论公式
(1)
对于圆柱壳，需要添加曲率修正项，文献 [11] 给出了适用于圆柱壳颤振分析的准定常活塞理论气动力公式
(2)
其中，R为半径。
根据唐纳尔简化理论，大开口圆柱壳或者圆柱形扁壳只包含法向位移的小挠度振动方程，写作如下形式 [12]
(3.a)
(3.b)
其中，，f为Airy应力函数，为弯曲刚度，E为杨氏模量，为泊松比，
h为壳的厚度，R为圆弧半径，为材料密度，x为轴向坐标，为周向坐标，t为时间变量，w为径向位移，f为Airy应力函数。
考虑四边简支的圆柱形扁壳，其边界条件可以写做下列形式
，时，
，时， (4)
由于在振动方程中不显含变量u和v，因此将和写成
(5)
式中，，，，，，，均表示内力分量。
将气动力公式(2)代入方程(3)，并引入下列无量纲量：
(6)
其中，L为轴向长度，为周向角度。无量纲处理后的圆柱形扁壳气动弹性方程如下：
(7.a)
(7.b)
相应边界条件表示为
，时，，即
；；；；
(8)
3. 二维微分求积法离散
依据微分求积法(以下简称DQM)二维问题的离散化方法 [13] ，在，的区域内，分别在和方向上置入Ｎ和Ｍ个网格点，于是在平行于的任一直线上，，函数在网格点处，对的一阶~四阶偏导可近似写为：
(9.a)
同样地，在任一圆弧上，，函数在网格点处，对的一阶~四阶偏导可近似写为
(9.b)
同理，对应力函数的各阶偏导也可以写做类似形式。其中，权系数A和B的确定参见文献 [13] 。
根据式(9.a)和(9.b)，可以将无量纲形式的圆柱形扁壳气动弹性方程(7)写为下列DQM离散后的形式：
(10.a)
(10.b)
边界条件的(8)的DQM形式如下：
，
， (11)
将方程(10)写成DQM的矩阵形式为：
(12.a)
(12.b)
这里为不含边界变量的无量纲位移矩阵，为不含边界变量的无量纲Airy应力函数矩阵，可写做下述形式：
和 (13)
引入边界条件(11)，一般矩阵形式的气动弹性方程可以写做
(14)
此处，[M]为质量阵，[C]为阻尼阵，[Rd]为不含边界变量的右刚度系数阵，[Ld]为不含边界变量的左刚度系数阵，[Rb]为边界变量的右刚度系数阵，[Lb]为边界变量的左刚度系数阵，为包含所有位移变量和应力变量的矩阵，第三项中的求和指标“5”表示方程12中的包含五项系数，第三项中的求和指标“16”表示对应的边界条件方程11中的系数。
利用“Kronecker”转换简化方程(14)的形式。原矩阵形式的位移变量矩阵可以展开为一维的列向量，相应的系数阵可以写成
(15)
这里，[I]为单位阵。由此方程可以简化为下面的形式
(16)
对方程(16)的求解最后转化为标准特征值问题。其中特征值虚部表示系统的无量纲自振频率。随着流速的增大，可以得到对应于动压的一系列特征值，当特征值实部大于零，且虚部不为零时，系统发生颤振失稳，特征值实部为零的动压就称为是颤振临界动压。
4. 曲率对颤振临界动压的影响
取计算参数如文献 [14] ，，，，保持其他参数不变，固定长宽比为1，也就是，反比例变化和。当很大，很小时，可以近似认为曲率为零，圆柱形扁壳方程可以近似作为平板方程。取，，，图1给出了特征值实虚部随无量纲动压变化的情况。与平板气动弹性颤振研究的文献 [14] 比较得知，两者无量纲频率变化规律及无量纲颤振临界动压基本一致，皆为。为研究一带有微小曲率的扁壳，取，，图2给出了特征值实虚部随无量纲动压变化的情况。对于大曲率的扁壳，取，，图3给出了特征值实虚部随无量纲动压变化的情况。
从图1至图3可以看出，在长宽比不变的情况下，随着曲率的增加，越来越多的低阶频率在周向上分布越密集，颤振的频率也随之增大。当时，在相对较大的动压范围内，发生耦合颤振的模态较为单一，随着的减小，越来越多的高阶模态加入到颤振的行列，使情况变得更加复杂，这对结构的颤振设计而言是非常不利的。图4描述了颤振临界动压与的关系曲线，可以看到，随着的减小，呈指数增长。由图可以看到，在保持不变的情况下，当较大时(对应于接近平板的情况)，较大幅度的减小并不能引起颤振临界动压参数的显著增长，这是因为在较大的范围内(图中)，系统仍就处于接近于平板的状态。只有当减小到一定程度时(图中)，曲率的变化才能显著影响颤振临界动压。
(a)(b)
Figure 1. The real parts and imaginary parts of eigenvalues vs. non-dimensional aerodynamic pressure as,
图1.，时，特征值实虚部随无量纲动压变化
Figure 2. The real parts and imaginary parts of eigenvalues vs. non-dimensional aerodynamic pressure as,
图2.，时，特征值实虚部随无量纲动压变化
Figure 3. The real parts and imaginary parts of eigenvalues vs. non-dimensional aerodynamic pressure as,
图3.，时，特征值实虚部随无量纲动压变化
Figure 4. Flutter critical aerodynamic pressure vs.
图4. 颤振临界动压随变化图
5. 结论
本文给出了超声速轴向流中，圆柱形扁壳气动弹性系统的小挠度运动方程，利用二维的微分求积法进行了离散，运用特征值方法求解了不同曲率对线性系统颤振临界动压的影响，结论如下：
1) 曲率很小的圆柱形扁壳结构的颤振临界动压接近于平板的颤振临界动压；
2) 在长宽比不变的情况下，随着曲率的增加，越来越多的低阶频率在周向上分布越密集，发生颤振的频率也随之增大。
3) 随着曲率的增加，系统颤振临界动压呈指数增长。
基金项目
国家自然科学基金(11302181)；高等学校博士点新教师基金(20110184120025)。
*通讯作者。
参考文献