韓勤鍇 秦朝燁 褚福磊

(清華大學(xué) 機(jī)械系，北京 100084)

引言

旋轉(zhuǎn)薄壁圓柱殼體結(jié)構(gòu)在航空發(fā)動(dòng)機(jī)等高速旋轉(zhuǎn)機(jī)械中得到大量應(yīng)用，其動(dòng)態(tài)特性長期以來受到國內(nèi)外學(xué)者的重視[1,2].早在上世紀(jì)60年代，DiTaranto和Lessen[3]，Srinivasan和Lauterbach[4]首次發(fā)現(xiàn)了旋轉(zhuǎn)殼的行波振動(dòng)現(xiàn)象和科氏力的影響.Huang和Soedel[5，6]研究了承受諧波移動(dòng)載荷的旋轉(zhuǎn)柱殼的強(qiáng)迫共振現(xiàn)象，并著重討論了科氏加速度的影響.基于Love殼理論和無網(wǎng)格方法，Hua和Lam[7]，Liew等[8,9]研究了邊界條件對薄壁旋轉(zhuǎn)柱殼固有頻率的影響規(guī)律.以航空發(fā)動(dòng)機(jī)旋轉(zhuǎn)薄壁鼓筒結(jié)構(gòu)為對象，不少國內(nèi)學(xué)者，如曹航等[10]、韓清凱等[11]、孫樹鵬等[12]、李文達(dá)等[13]，分別開展了較為系統(tǒng)的理論分析與試驗(yàn)測試研究.上述研究均考慮轉(zhuǎn)速是恒定不變的.

實(shí)際中，由于外界的干擾以及發(fā)動(dòng)機(jī)輸出轉(zhuǎn)速的波動(dòng)，導(dǎo)致旋轉(zhuǎn)殼體的轉(zhuǎn)速是隨時(shí)間變化的.時(shí)變轉(zhuǎn)速將導(dǎo)致系統(tǒng)剛度和陀螺系數(shù)也是時(shí)變的.這類系統(tǒng)通常又被稱之為參數(shù)振動(dòng)系統(tǒng)[14].由于時(shí)變剛度和陀螺系數(shù)引起的參數(shù)激勵(lì)將使得系統(tǒng)動(dòng)態(tài)特性發(fā)生重要變化，在某些工況下，系統(tǒng)可能失穩(wěn)并出現(xiàn)劇烈的振動(dòng)[15].然而，具有周期轉(zhuǎn)速的旋轉(zhuǎn)殼體動(dòng)力穩(wěn)定性問題尚未得到重視.

事實(shí)上，不少類型的非穩(wěn)態(tài)旋轉(zhuǎn)結(jié)構(gòu)，例如旋轉(zhuǎn)梁[16]、旋轉(zhuǎn)盤[17]和葉盤結(jié)構(gòu)[18，19]等，已有不少學(xué)者深入探討了變轉(zhuǎn)速誘發(fā)的動(dòng)力失穩(wěn)問題.眾所周知，在周期軸向力作用下，殼體結(jié)構(gòu)也可能出現(xiàn)劇烈的橫向振動(dòng)甚至失穩(wěn)現(xiàn)象.文獻(xiàn)[20]詳細(xì)評述了這方面研究進(jìn)展.盡管Ng等[21，22]和Liew等[23]已經(jīng)報(bào)道了旋轉(zhuǎn)效應(yīng)對于周期軸向力作用下圓柱殼的參數(shù)穩(wěn)定性的影響，但是在他們的研究中轉(zhuǎn)速是恒定不變的.若同時(shí)考慮變轉(zhuǎn)速和周期軸向力，系統(tǒng)將同時(shí)存在兩個(gè)參數(shù)激勵(lì)，且可能具有不同的幅值和相位.作者前期的理論和實(shí)驗(yàn)研究已經(jīng)表明[24]: 在一定范圍內(nèi)調(diào)節(jié)兩個(gè)參數(shù)激勵(lì)源的相位，將使得系統(tǒng)不穩(wěn)定區(qū)減小甚至消失.因此，可以預(yù)見的是同時(shí)考慮兩種時(shí)變因素將使得系統(tǒng)動(dòng)力穩(wěn)定性發(fā)生顯著變化.

因此，本文提出了變速旋轉(zhuǎn)圓柱薄殼動(dòng)力穩(wěn)定性研究.基于Donnell薄殼理論[25]建立同時(shí)考慮變轉(zhuǎn)速和周期軸向力的圓柱殼體振動(dòng)微分方程；采用多尺度方法，推導(dǎo)系統(tǒng)主參數(shù)不穩(wěn)定區(qū)和組合不穩(wěn)定區(qū)邊界的解析表達(dá)式；分別探討僅考慮周期軸向力、僅考慮變轉(zhuǎn)速以及同時(shí)考慮兩種時(shí)變因素時(shí)，系統(tǒng)主參數(shù)不穩(wěn)定區(qū)和組合不穩(wěn)定區(qū)的變化規(guī)律；通過與文獻(xiàn)結(jié)果以及數(shù)值結(jié)果的對比，驗(yàn)證解析結(jié)果的準(zhǔn)確性.

1 分析模型

圖1給出了承受周期軸向力的變速旋轉(zhuǎn)圓柱殼體結(jié)構(gòu)，其長度為L，厚度為h，半徑為R，楊氏模量為E，密度為ρ，泊松比為υ.柱殼的軸向、周向和徑向振動(dòng)分別由u,v,w表示.周期軸向力ηa(t)和變轉(zhuǎn)速Ω(t)分別假定為簡諧波形式，頻率均為ω但相位不同(分別為φa和φr)，表達(dá)式表示如下

圖1 承受周期軸向力的變速旋轉(zhuǎn)圓柱殼體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表示軸向壓應(yīng)力失穩(wěn)閾值，η0表示無量綱軸向力幅值，εa表示軸向力變化的相對幅值，Ω0表示轉(zhuǎn)速最大值，εr表示變轉(zhuǎn)速的相對幅值.

根據(jù)Donnell薄殼理論[25]，圖1所示的旋轉(zhuǎn)圓柱殼體結(jié)構(gòu)的動(dòng)力學(xué)方程可表示為

(2)

式中，u=[u,v,w]T，γ=ρR2(1υ2)/E.L，LR，LP分別表示非旋轉(zhuǎn)殼、轉(zhuǎn)速和周向力引起的微分算子矩陣，可表示為

(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)

考慮到柱殼的邊界條件為兩端簡支，則其振動(dòng)位移可表示為

(9)

其中qmnj(t),pmnj(t)表示一般坐標(biāo).引入WI=cosnθsinλmx，VI=sinnθsinλmx和模態(tài)振型比αI=Amnj/Cmnj，βI=Bmnj/Cmnj，系統(tǒng)模態(tài)振型可表示為

(10)

將式(9)代入式(2)中，并左乘和右乘模態(tài)振型，可得到系統(tǒng)兩耦合的二階振動(dòng)微分方程如下所示

(11a)

(11b)

其中i=1,2,…,I，各個(gè)參數(shù)可表示為

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)

考慮轉(zhuǎn)速后，薄壁柱殼結(jié)構(gòu)成為陀螺系統(tǒng).原系統(tǒng)的固有模態(tài)將分為兩個(gè)：前向和后向行波模態(tài).因陀螺力的存在使得兩個(gè)模態(tài)存在耦合現(xiàn)象，如式(11)所示.不僅如此，由于變轉(zhuǎn)速和周期軸向力的作用，薄壁柱殼的控制方程中出現(xiàn)了時(shí)變陀螺和剛度系數(shù).這種系統(tǒng)是一種典型的參數(shù)激勵(lì)陀螺系統(tǒng)，其動(dòng)力穩(wěn)定性是現(xiàn)有研究的關(guān)鍵問題.

2 多尺度方法求解系統(tǒng)穩(wěn)定性

本節(jié)將采用多尺度方法求解系統(tǒng)的動(dòng)力穩(wěn)定性問題.顯然，存在兩個(gè)參數(shù)激勵(lì)源，分別由周期軸向力和時(shí)變轉(zhuǎn)速所引起.它們的相對幅值由εa和εr所決定.在本研究中，假定二者具有相同的階次，即εr=ε，εr=με，其中μ=Ο(1).通過引入?yún)?shù)：λ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)

采用多尺度方法，系統(tǒng)穩(wěn)態(tài)解可表示為

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 組合不穩(wěn)定區(qū)邊界

這里首先考慮組合不穩(wěn)定區(qū)邊界的推導(dǎo).令ω=ω1+ω2+εσ，其中σ為待定的調(diào)節(jié)參數(shù).因此，我們可以得到這樣的關(guān)系：ei(ω-ω1)T0=ei(ω2T0+σT1)和ei(ω-ω2)T0=ei(ω1T0+σT1).由于式(21)和式(22)均是線性方程，因此特解可分別得到.為了確定參數(shù)不穩(wěn)定條件，我們尋找具有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)

根據(jù)式(20)，可知式(24)和式(25)的系數(shù)矩陣奇異.因此，式(24)和式(25)的解不存在，除了

(28a)

(28b)

因此，可以得到

(29a)

(29b)

考慮到式(20)所給出的關(guān)系，可以得到

(30a)

(30b)

式(29)的非平凡解為

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

(31)

其中

(32a)

Λ=4Γ1Γ2

(32b)

從式(32)可知，當(dāng)σ2≥Λ時(shí)，A1,A2是有界的，反之，A1,A2是無界的.從式(30)可以看出，當(dāng)fv3≠0時(shí)，Γ1,Γ2是復(fù)數(shù)，σ2<Λ總是滿足的.因此，在這種情況下，系統(tǒng)是不穩(wěn)定的.通過檢查fv3的表達(dá)式，我們可以進(jìn)一步得到系統(tǒng)保持不穩(wěn)定的條件是：φ≠0,π, 2π, ….當(dāng)φ=0,π, 2π, …時(shí)，可得到fv3=0和

(33)

因此，不穩(wěn)定區(qū)的邊界可表示為

(34)

2.2 主參數(shù)不穩(wěn)定區(qū)邊界

按照上一節(jié)中相似的推導(dǎo)流程，可得到主參數(shù)不穩(wěn)定區(qū)的邊界表達(dá)式為

(35)

(36)

其中式(35)表示的是前向行波模態(tài)，而式(36)則是后向行波模態(tài).

3 數(shù)值算例分析

3.1 模型驗(yàn)證

在開展動(dòng)力穩(wěn)定性研究前，需驗(yàn)證本文所建立的旋轉(zhuǎn)殼體動(dòng)力學(xué)模型.以模態(tài)(1,1)和(1,3)為例，計(jì)算了前向和后向行波模態(tài)頻率隨轉(zhuǎn)速的變化，如表1和表2所示.分析的圓柱薄殼參數(shù)為：h=R/500，L=5R，E=206 GPa，ρ=7850 kg/m3，υ=0.3.作為對比，表中也給出了文獻(xiàn)[26]的計(jì)算結(jié)果，且進(jìn)行了無量綱化處理，量綱為(ρR2(1υ2)/E)1/2.可以看出，本文結(jié)果與文獻(xiàn)結(jié)果符合很好，從而驗(yàn)證了旋轉(zhuǎn)殼體模型的準(zhǔn)確性.

表1 模態(tài)(1,1)的前向和后向行波頻率隨轉(zhuǎn)速變化結(jié)果Table 1 Results of forward and backward traveling wave frequencies of mode (1,1) rotating speed

表2 模態(tài)(1,3)的前向和后向行波頻率隨轉(zhuǎn)速變化結(jié)果Table 2 Results of forward and backward traveling wave frequencies of mode (1,3) rotating speed

3.2 僅考慮周期軸向力作用

本節(jié)只考慮周期軸向力作用，而轉(zhuǎn)速則被認(rèn)為是恒定不變的，典型結(jié)果如圖2和圖3所示.圓柱薄殼參數(shù)為：h=R/100，L=2R，E=206 GPa，ρ=7850 kg/m3，υ=0.3.固有頻率和轉(zhuǎn)速均做了無量綱化處理.為了驗(yàn)證解析結(jié)果，本文引入文獻(xiàn)[15]所介紹的數(shù)值方法，在所關(guān)心的參數(shù)域內(nèi)，通過求解離散狀態(tài)轉(zhuǎn)移矩陣的特征值問題，逐點(diǎn)判斷系統(tǒng)穩(wěn)定性.可以看出，解析解與數(shù)值結(jié)果符合很好，說明本文推導(dǎo)的不穩(wěn)定邊界表達(dá)式可信.考慮旋轉(zhuǎn)后，固有模態(tài)被分為前向和后向兩個(gè)行波模態(tài).通常情況下，將存在三個(gè)不穩(wěn)定區(qū)，其中兩個(gè)是主參數(shù)不穩(wěn)定模態(tài)，而另一個(gè)則是組合不穩(wěn)定模態(tài).不穩(wěn)定區(qū)的起點(diǎn)可表示為：2ω1*, 2ω2*和ω1*+ω2*，同時(shí)也標(biāo)注在圖中.但是，對于周期軸向力作用下的旋轉(zhuǎn)殼，僅發(fā)現(xiàn)了組合不穩(wěn)定區(qū)，而兩個(gè)主參數(shù)不穩(wěn)定區(qū)并未出現(xiàn).此外，從圖2和圖3可以看出，本文結(jié)果與文獻(xiàn)結(jié)果存在明顯差異.隨著轉(zhuǎn)速的增加，文獻(xiàn)[27]所給出的不穩(wěn)定區(qū)位置幾乎保持不變.這是不合理的，因?yàn)檗D(zhuǎn)速的增加必然導(dǎo)致固有頻率的變化，進(jìn)而使得不穩(wěn)定區(qū)位置發(fā)生整體移動(dòng).不僅如此，對于模態(tài)(1,1)，文獻(xiàn)[27]的結(jié)果表明，考慮轉(zhuǎn)速后，不穩(wěn)定區(qū)的起點(diǎn)為一區(qū)域而不是一點(diǎn)，如圖2(b)和圖2(c)所示.這同樣是不合理的，因?yàn)楫?dāng)系統(tǒng)不存在參數(shù)激勵(lì)時(shí)，不穩(wěn)定區(qū)應(yīng)為一點(diǎn)，而不是一個(gè)區(qū)域.

圖2 僅考慮周期軸向力時(shí)關(guān)于模態(tài)(1,1)的不穩(wěn)定區(qū)：(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 僅考慮周期軸向力時(shí)關(guān)于模態(tài)(1,2)的不穩(wěn)定區(qū)：(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 僅考慮周期轉(zhuǎn)速

本節(jié)僅考慮周期轉(zhuǎn)速的影響.求得對應(yīng)模態(tài)(1,1)，(1,2)，(1,3)和(1,4)的系統(tǒng)不穩(wěn)定區(qū)，分別如圖4和圖5所示.圖中“P1”和“P2”分別表示對應(yīng)于前向和后向行波模態(tài)的主參數(shù)不穩(wěn)定區(qū)，而“C”則表示組合不穩(wěn)定區(qū).從圖中可以看出，理論結(jié)果與數(shù)值結(jié)果符合很好，表明解析結(jié)果可信.

圖4 僅考慮變轉(zhuǎn)速時(shí)不穩(wěn)定區(qū)：(a)模態(tài)(1,1)；(b)模態(tài)(1,2)Fig.4 Instability regions for only variable rotations: (a) mode (1,1); (b) mode (1,2)

圖5 僅考慮變轉(zhuǎn)速時(shí)不穩(wěn)定區(qū)： (a)模態(tài)(1,3)；(b)模態(tài)(1,4)Fig.5 Instability regions for only variable rotations: (a) mode (1,3); (b) mode (1,3)

3.4 同時(shí)考慮變轉(zhuǎn)速和周期軸向力

圖6 考慮軸向壓縮力時(shí)對應(yīng)模態(tài)(1,3)的不穩(wěn)定區(qū)間隨相位和相對幅值的變化：(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)中，可以發(fā)現(xiàn)εa的值需要大于A/B×εr=4.0243×εr≈0.8，以使得組合不穩(wěn)定區(qū)出現(xiàn).顯然地在0.05范圍內(nèi)，不存在組合不穩(wěn)定區(qū).對于圖6(a)，圖8(a)和圖9(b)，其反映的是同一類情況：即增加εa的值，不穩(wěn)定區(qū)將逐漸減小.若εa的值大于或等于特定值，不穩(wěn)定區(qū)將消失.對于圖6(a)，可以得到該臨界值為C/B×εr≈0.018.而對于圖8(a)和圖9(b)，該值可計(jì)算得到C/B×εr≈0.235和C/B×εr≈0.217.發(fā)現(xiàn)用于模態(tài)(1,4)不穩(wěn)定區(qū)消失的臨界εa值要遠(yuǎn)遠(yuǎn)大于模態(tài)(1,3)的情況，說明模態(tài)(1,4)的不穩(wěn)定區(qū)更難以抑制.圖6(b)，圖8(b)和圖9(a)所反映的是一類情況，而圖7(a)則屬于另一類情況.我們可以從圖7(a)發(fā)現(xiàn)，εa的值要大于C/B×εr≈0.004以使得不穩(wěn)定區(qū)出現(xiàn).隨著εa的增加，這些圖表明組合不穩(wěn)定區(qū)不斷放大，特別是對模態(tài)(1,4).因此，可以總結(jié)對于圖6(b)，圖8(b)和圖9(a)所示的情況，施加周期軸向力將使組合不穩(wěn)定區(qū)放大.

圖7 考慮軸向拉伸力時(shí)對應(yīng)模態(tài)(1,3)的不穩(wěn)定區(qū)間隨相位和相對幅值的變化：(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 考慮軸向壓縮力時(shí)對應(yīng)模態(tài)(1,4)的不穩(wěn)定區(qū)間隨相位和相對幅值的變化：(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 考慮軸向拉伸力時(shí)對應(yīng)模態(tài)(1,4)的不穩(wěn)定區(qū)間隨相位和相對幅值的變化：(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 結(jié)論

針對變速旋轉(zhuǎn)圓柱殼體動(dòng)力穩(wěn)定性問題，采用多尺度方法研究了周期軸向力和變轉(zhuǎn)速的影響規(guī)律，并通過與文獻(xiàn)結(jié)果以及數(shù)值結(jié)果的對比，驗(yàn)證解析結(jié)果的準(zhǔn)確性.主要結(jié)論如下：

(1)考慮旋轉(zhuǎn)后，周期軸向力作用下的圓柱殼體僅存在組合不穩(wěn)定區(qū).增加轉(zhuǎn)速將使得不穩(wěn)定區(qū)的整體平移，但其對不穩(wěn)定區(qū)寬度影響不大.增加軸向力將不僅增加不穩(wěn)定區(qū)寬度，而且同時(shí)導(dǎo)致不穩(wěn)定區(qū)的整體移動(dòng).

(2)僅考慮變轉(zhuǎn)速時(shí)，系統(tǒng)總存在主參數(shù)不穩(wěn)定區(qū)，其寬度將隨轉(zhuǎn)速的增加而增大.對于特定模態(tài)，組合不穩(wěn)定區(qū)可能不存在，其存在條件可通過解析方法進(jìn)行判斷.增加轉(zhuǎn)速值將增加不穩(wěn)定區(qū)的寬度.而且，也可能導(dǎo)致組合不穩(wěn)定區(qū)的出現(xiàn).

(3)當(dāng)變轉(zhuǎn)速和周期軸向力同時(shí)作用時(shí)，主要影響組合不穩(wěn)定區(qū)，而主參數(shù)不穩(wěn)定區(qū)則幾乎不受影響.某些條件下，變轉(zhuǎn)速引起的組合不穩(wěn)定區(qū)將被周期軸向力削弱(甚至消失).但是，在另一些條件下，該不穩(wěn)定區(qū)只會(huì)被不斷放大.上述削弱或放大的條件也已通過多尺度方法解析給出，并通過數(shù)值分析進(jìn)行了驗(yàn)證.

分析均是基于線性和各向同性圓柱殼體模型.若殼體是非線性的且各向異性，例如考慮大變形、功能梯度或復(fù)合材料層合殼材料.非線性因素的引入將使得系統(tǒng)動(dòng)力穩(wěn)定性產(chǎn)生變化，這方面的工作尚待深入開展.