变速旋转圆柱薄壳动力稳定性研究*

韩勤锴 秦朝烨 褚福磊

(清华大学 机械系，北京 100084)

引言

旋转薄壁圆柱壳体结构在航空发动机等高速旋转机械中得到大量应用，其动态特性长期以来受到国内外学者的重视[1,2].早在上世纪60年代，DiTaranto和Lessen[3]，Srinivasan和Lauterbach[4]首次发现了旋转壳的行波振动现象和科氏力的影响.Huang和Soedel[5，6]研究了承受谐波移动载荷的旋转柱壳的强迫共振现象，并着重讨论了科氏加速度的影响.基于Love壳理论和无网格方法，Hua和Lam[7]，Liew等[8,9]研究了边界条件对薄壁旋转柱壳固有频率的影响规律.以航空发动机旋转薄壁鼓筒结构为对象，不少国内学者，如曹航等[10]、韩清凯等[11]、孙树鹏等[12]、李文达等[13]，分别开展了较为系统的理论分析与试验测试研究.上述研究均考虑转速是恒定不变的.

实际中，由于外界的干扰以及发动机输出转速的波动，导致旋转壳体的转速是随时间变化的.时变转速将导致系统刚度和陀螺系数也是时变的.这类系统通常又被称之为参数振动系统[14].由于时变刚度和陀螺系数引起的参数激励将使得系统动态特性发生重要变化，在某些工况下，系统可能失稳并出现剧烈的振动[15].然而，具有周期转速的旋转壳体动力稳定性问题尚未得到重视.

事实上，不少类型的非稳态旋转结构，例如旋转梁[16]、旋转盘[17]和叶盘结构[18，19]等，已有不少学者深入探讨了变转速诱发的动力失稳问题.众所周知，在周期轴向力作用下，壳体结构也可能出现剧烈的横向振动甚至失稳现象.文献[20]详细评述了这方面研究进展.尽管Ng等[21，22]和Liew等[23]已经报道了旋转效应对于周期轴向力作用下圆柱壳的参数稳定性的影响，但是在他们的研究中转速是恒定不变的.若同时考虑变转速和周期轴向力，系统将同时存在两个参数激励，且可能具有不同的幅值和相位.作者前期的理论和实验研究已经表明[24]: 在一定范围内调节两个参数激励源的相位，将使得系统不稳定区减小甚至消失.因此，可以预见的是同时考虑两种时变因素将使得系统动力稳定性发生显著变化.

因此，本文提出了变速旋转圆柱薄壳动力稳定性研究.基于Donnell薄壳理论[25]建立同时考虑变转速和周期轴向力的圆柱壳体振动微分方程；采用多尺度方法，推导系统主参数不稳定区和组合不稳定区边界的解析表达式；分别探讨仅考虑周期轴向力、仅考虑变转速以及同时考虑两种时变因素时，系统主参数不稳定区和组合不稳定区的变化规律；通过与文献结果以及数值结果的对比，验证解析结果的准确性.

1 分析模型

图1给出了承受周期轴向力的变速旋转圆柱壳体结构，其长度为L，厚度为h，半径为R，杨氏模量为E，密度为ρ，泊松比为υ.柱壳的轴向、周向和径向振动分别由u,v,w表示.周期轴向力ηa(t)和变转速Ω(t)分别假定为简谐波形式，频率均为ω但相位不同(分别为φa和φr)，表达式表示如下

图1 承受周期轴向力的变速旋转圆柱壳体Fig.1 Thin cylindrical shell under variable speed rotation subjected to periodic axial force

ηa(t)=η0ηcr(1+εacos(ωt+φa))

(1a)

Ω(t)=Ω0(1+εrcos(ωt+φr))

(1b)

式中ηcr表示轴向压应力失稳阈值，η0表示无量纲轴向力幅值，εa表示轴向力变化的相对幅值，Ω0表示转速最大值，εr表示变转速的相对幅值.

根据Donnell薄壳理论[25]，图1所示的旋转圆柱壳体结构的动力学方程可表示为

(2)

式中，u=[u,v,w]T，γ=ρR2(1υ2)/E.L，LR，LP分别表示非旋转壳、转速和周向力引起的微分算子矩阵，可表示为

(3)

式中

(4a)

(4b)

(4c)

(4d)

(4e)

(4f)

(5)

式中

(6a)

(6b)

(6c)

(6d)

(6f)

(7)

其中

LP11=LP12=LP13=LP21=0

(8a)

LP22=LP23=LP31=LP32=0

(8b)

(8c)

考虑到柱壳的边界条件为两端简支，则其振动位移可表示为

(9)

其中qmnj(t),pmnj(t)表示一般坐标.引入WI=cosnθsinλmx，VI=sinnθsinλmx和模态振型比αI=Amnj/Cmnj，βI=Bmnj/Cmnj，系统模态振型可表示为

(10)

将式(9)代入式(2)中，并左乘和右乘模态振型，可得到系统两耦合的二阶振动微分方程如下所示

(11a)

(11b)

其中i=1,2,…,I，各个参数可表示为

Mi=γ(πL/2)(1-ν2)(1+βIβI+αIαI)

(12a)

Gi=-γ(πL/2)(1-ν2)(βI+βI)

(12b)

K0i=(πL/2)K*

(12c)

(12d)

Ksi=γ(πL/4)(βI+βI)

(12e)

Kpi=(πL/4)(R2λmλm)

(12f)

其中

(13)

(14)

考虑转速后，薄壁柱壳结构成为陀螺系统.原系统的固有模态将分为两个：前向和后向行波模态.因陀螺力的存在使得两个模态存在耦合现象，如式(11)所示.不仅如此，由于变转速和周期轴向力的作用，薄壁柱壳的控制方程中出现了时变陀螺和刚度系数.这种系统是一种典型的参数激励陀螺系统，其动力稳定性是现有研究的关键问题.

2 多尺度方法求解系统稳定性

本节将采用多尺度方法求解系统的动力稳定性问题.显然，存在两个参数激励源，分别由周期轴向力和时变转速所引起.它们的相对幅值由εa和εr所决定.在本研究中，假定二者具有相同的阶次，即εr=ε，εr=με，其中μ=Ο(1).通过引入参数：λ0=2GiΩ0/Mi，α0=(K0i+KciΩ02+2Kpiη0ηcr+ε2KciΩ02/2)/Mi，fv1=(KciΩ02+μKpiη0ηcrcosφ)/Mi，fv2=KsiΩ0/Mi，fv3=μKpiη0ηcrsinφ/Mi，fv4=KciΩ02/(2Mi)，式(11)可重新表示为

(15a)

(15b)

采用多尺度方法，系统稳态解可表示为

q=q0(T0,T1)+εq0(T0,T1)+…

(16a)

p=p0(T0,T1)+εp0(T0,T1)+…

(16b)

将式(16)代入式(15)中，并按照ε的阶次整理可得到

(17a)

(17b)

2(ωfv2+fv3)q0sinωT0

(18a)

(18b)

式(17)的一般解可表示为

q0=A1(T1)eiω1T0+A2(T1)eiω2T0+cc

(19a)

(19b)

其中ω12和ω22是如下方程的解

(20)

这里假定ω12和ω22是正值，且ω2>ω1.将式(19)代入式(18)可得到

(21)

(22)

2.1 组合不稳定区边界

这里首先考虑组合不稳定区边界的推导.令ω=ω1+ω2+εσ，其中σ为待定的调节参数.因此，我们可以得到这样的关系：ei(ω-ω1)T0=ei(ω2T0+σT1)和ei(ω-ω2)T0=ei(ω1T0+σT1).由于式(21)和式(22)均是线性方程，因此特解可分别得到.为了确定参数不稳定条件，我们寻找具有eiω1T0，eiω2T0形式的特解

q1=P1(T1)eiω1T0+Q1(T1)eiω2T0

(23a)

p=P2(T1)eiω1T0+Q2(T1)eiω2T0

(23b)

将式(23)代入式(21)和式(22)中，可得到

(24a)

(24b)

和

(25a)

(25b)

其中

(26a)

(26b)

(27a)

(27b)

根据式(20)，可知式(24)和式(25)的系数矩阵奇异.因此，式(24)和式(25)的解不存在，除了

(28a)

(28b)

因此，可以得到

(29a)

(29b)

考虑到式(20)所给出的关系，可以得到

(30a)

(30b)

式(29)的非平凡解为

A1=a1eiκT1,A2=a2ei(κ+σ)T1

(31)

其中

(32a)

Λ=4Γ1Γ2

(32b)

从式(32)可知，当σ2≥Λ时，A1,A2是有界的，反之，A1,A2是无界的.从式(30)可以看出，当fv3≠0时，Γ1,Γ2是复数，σ2<Λ总是满足的.因此，在这种情况下，系统是不稳定的.通过检查fv3的表达式，我们可以进一步得到系统保持不稳定的条件是：φ≠0,π, 2π, ….当φ=0,π, 2π, …时，可得到fv3=0和

(33)

因此，不稳定区的边界可表示为

(34)

2.2 主参数不稳定区边界

按照上一节中相似的推导流程，可得到主参数不稳定区的边界表达式为

(35)

(36)

其中式(35)表示的是前向行波模态，而式(36)则是后向行波模态.

3 数值算例分析

3.1 模型验证

在开展动力稳定性研究前，需验证本文所建立的旋转壳体动力学模型.以模态(1,1)和(1,3)为例，计算了前向和后向行波模态频率随转速的变化，如表1和表2所示.分析的圆柱薄壳参数为：h=R/500，L=5R，E=206 GPa，ρ=7850 kg/m3，υ=0.3.作为对比，表中也给出了文献[26]的计算结果，且进行了无量纲化处理，量纲为(ρR2(1υ2)/E)1/2.可以看出，本文结果与文献结果符合很好，从而验证了旋转壳体模型的准确性.

表1 模态(1,1)的前向和后向行波频率随转速变化结果Table 1 Results of forward and backward traveling wave frequencies of mode (1,1) rotating speed

表2 模态(1,3)的前向和后向行波频率随转速变化结果Table 2 Results of forward and backward traveling wave frequencies of mode (1,3) rotating speed

3.2 仅考虑周期轴向力作用

本节只考虑周期轴向力作用，而转速则被认为是恒定不变的，典型结果如图2和图3所示.圆柱薄壳参数为：h=R/100，L=2R，E=206 GPa，ρ=7850 kg/m3，υ=0.3.固有频率和转速均做了无量纲化处理.为了验证解析结果，本文引入文献[15]所介绍的数值方法，在所关心的参数域内，通过求解离散状态转移矩阵的特征值问题，逐点判断系统稳定性.可以看出，解析解与数值结果符合很好，说明本文推导的不稳定边界表达式可信.考虑旋转后，固有模态被分为前向和后向两个行波模态.通常情况下，将存在三个不稳定区，其中两个是主参数不稳定模态，而另一个则是组合不稳定模态.不稳定区的起点可表示为：2ω1*, 2ω2*和ω1*+ω2*，同时也标注在图中.但是，对于周期轴向力作用下的旋转壳，仅发现了组合不稳定区，而两个主参数不稳定区并未出现.此外，从图2和图3可以看出，本文结果与文献结果存在明显差异.随着转速的增加，文献[27]所给出的不稳定区位置几乎保持不变.这是不合理的，因为转速的增加必然导致固有频率的变化，进而使得不稳定区位置发生整体移动.不仅如此，对于模态(1,1)，文献[27]的结果表明，考虑转速后，不稳定区的起点为一区域而不是一点，如图2(b)和图2(c)所示.这同样是不合理的，因为当系统不存在参数激励时，不稳定区应为一点，而不是一个区域.

图2 仅考虑周期轴向力时关于模态(1,1)的不稳定区：(a)Ω*=0；(b)Ω*=0.1；(c)Ω*=0.2Fig.2 Instability region of mode (1,1) when only periodic axial force is considered: (a)Ω*=0；(b)Ω*=0.1；(c)Ω*=0.2

图3 仅考虑周期轴向力时关于模态(1,2)的不稳定区：(a)Ω*=0；(b)Ω*=0.1；(c)Ω*=0.2Fig.2 Instability region of mode (1,2) when only periodic axial force is considered: (a)Ω*=0；(b)Ω*=0.1；(c)Ω*=0.2

3.3 仅考虑周期转速

本节仅考虑周期转速的影响.求得对应模态(1,1)，(1,2)，(1,3)和(1,4)的系统不稳定区，分别如图4和图5所示.图中“P1”和“P2”分别表示对应于前向和后向行波模态的主参数不稳定区，而“C”则表示组合不稳定区.从图中可以看出，理论结果与数值结果符合很好，表明解析结果可信.

图4 仅考虑变转速时不稳定区：(a)模态(1,1)；(b)模态(1,2)Fig.4 Instability regions for only variable rotations: (a) mode (1,1); (b) mode (1,2)

图5 仅考虑变转速时不稳定区： (a)模态(1,3)；(b)模态(1,4)Fig.5 Instability regions for only variable rotations: (a) mode (1,3); (b) mode (1,3)

3.4 同时考虑变转速和周期轴向力

图6 考虑轴向压缩力时对应模态(1,3)的不稳定区间随相位和相对幅值的变化：(a) φ=0；(b) φ=180Fig.6 Instability regions of modes (1,3) varies with phase and relative amplitude for axial compression force: (a) φ=0; (b) φ=180

在图7(b)中，可以发现εa的值需要大于A/B×εr=4.0243×εr≈0.8，以使得组合不稳定区出现.显然地在0.05范围内，不存在组合不稳定区.对于图6(a)，图8(a)和图9(b)，其反映的是同一类情况：即增加εa的值，不稳定区将逐渐减小.若εa的值大于或等于特定值，不稳定区将消失.对于图6(a)，可以得到该临界值为C/B×εr≈0.018.而对于图8(a)和图9(b)，该值可计算得到C/B×εr≈0.235和C/B×εr≈0.217.发现用于模态(1,4)不稳定区消失的临界εa值要远远大于模态(1,3)的情况，说明模态(1,4)的不稳定区更难以抑制.图6(b)，图8(b)和图9(a)所反映的是一类情况，而图7(a)则属于另一类情况.我们可以从图7(a)发现，εa的值要大于C/B×εr≈0.004以使得不稳定区出现.随着εa的增加，这些图表明组合不稳定区不断放大，特别是对模态(1,4).因此，可以总结对于图6(b)，图8(b)和图9(a)所示的情况，施加周期轴向力将使组合不稳定区放大.

图7 考虑轴向拉伸力时对应模态(1,3)的不稳定区间随相位和相对幅值的变化：(a) φ=0；(b) φ=180Fig.7 Instability regions of modes (1,3) varies with phase and relative amplitude for axial tensile force: (a) φ=0; (b) φ=180

图8 考虑轴向压缩力时对应模态(1,4)的不稳定区间随相位和相对幅值的变化：(a) φ=0；(b) φ=180Fig.8 Instability regions of modes (1,4) varies with phase and relative amplitude for axial compression force: (a) φ=0; (b) φ=180

图9 考虑轴向拉伸力时对应模态(1,4)的不稳定区间随相位和相对幅值的变化：(a) φ=0；(b) φ=180Fig.9 Instability regions of modes (1,4) varies with phase and relative amplitude for axial tensile force: (a) φ=0; (b) φ=180

4 结论

针对变速旋转圆柱壳体动力稳定性问题，采用多尺度方法研究了周期轴向力和变转速的影响规律，并通过与文献结果以及数值结果的对比，验证解析结果的准确性.主要结论如下：

(1)考虑旋转后，周期轴向力作用下的圆柱壳体仅存在组合不稳定区.增加转速将使得不稳定区的整体平移，但其对不稳定区宽度影响不大.增加轴向力将不仅增加不稳定区宽度，而且同时导致不稳定区的整体移动.

(2)仅考虑变转速时，系统总存在主参数不稳定区，其宽度将随转速的增加而增大.对于特定模态，组合不稳定区可能不存在，其存在条件可通过解析方法进行判断.增加转速值将增加不稳定区的宽度.而且，也可能导致组合不稳定区的出现.

(3)当变转速和周期轴向力同时作用时，主要影响组合不稳定区，而主参数不稳定区则几乎不受影响.某些条件下，变转速引起的组合不稳定区将被周期轴向力削弱(甚至消失).但是，在另一些条件下，该不稳定区只会被不断放大.上述削弱或放大的条件也已通过多尺度方法解析给出，并通过数值分析进行了验证.

分析均是基于线性和各向同性圆柱壳体模型.若壳体是非线性的且各向异性，例如考虑大变形、功能梯度或复合材料层合壳材料.非线性因素的引入将使得系统动力稳定性产生变化，这方面的工作尚待深入开展.