韓勤鍇 秦朝燁 褚福磊

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

引言

薄壁殼體結(jié)構(gòu)在土木工程、機(jī)械工程和航空航天工程中有著廣泛的應(yīng)用.在載荷參數(shù)值和橫向振動固有頻率的特定組合下，殼體結(jié)構(gòu)在面內(nèi)周期力作用下可能發(fā)生不穩(wěn)定的橫向振動，導(dǎo)致參數(shù)失穩(wěn)，危及結(jié)構(gòu)安全.半個(gè)世紀(jì)以來，具有不同幾何、邊界條件和荷載類型的殼體結(jié)構(gòu)參數(shù)失穩(wěn)問題得到了廣泛研究[1,2].

從研究對象看，現(xiàn)有研究[3-7]多數(shù)集中于圓柱殼，而對圓錐殼的參數(shù)失穩(wěn)研究相對較少.基于Marguerre型動力方程，Ye[8]分析了扁錐薄殼在周期性橫向和面內(nèi)載荷作用下的非線性振動和動力失穩(wěn).Ng等[9]利用廣義微分求積(GDQ)方法研究了邊界條件對周期邊緣載荷作用下截錐殼參數(shù)失穩(wěn)的影響.Sofiyev[10-12]對復(fù)合材料層合錐殼的熱致動力失穩(wěn)和屈曲行為進(jìn)行了一系列研究.近期，不少學(xué)者探索了由碳納米管增強(qiáng)復(fù)合材料[13]、功能梯度材料[14]以及金屬泡沫材料[15]制成的截錐殼非線性動力失穩(wěn)問題.

上述研究并未考慮旋轉(zhuǎn)效應(yīng).對于周期性軸向載荷作用下的旋轉(zhuǎn)圓柱殼，Ng等[16]和Liew等[17]報(bào)道了旋轉(zhuǎn)引起的科里奧利力和離心力對失穩(wěn)區(qū)域有顯著影響.Lam和Hua[18]指出，將研究從旋轉(zhuǎn)圓柱殼擴(kuò)展到旋轉(zhuǎn)圓錐殼所需的分析量相當(dāng)大.此外，圓錐殼是更為一般的殼體結(jié)構(gòu)形式.因此，研究旋轉(zhuǎn)對周期軸向載荷作用下旋轉(zhuǎn)錐殼參數(shù)失穩(wěn)的影響具有重要意義.

GDQ方法最早由Shu和Richards[19]提出，用于直接求解工程問題的控制方程.目前，GDQ方法已廣泛應(yīng)用于殼體結(jié)構(gòu)的振動分析[20-22].本文將基于該方法建立旋轉(zhuǎn)錐殼的振動模型.考慮旋轉(zhuǎn)效應(yīng)后，周期軸向載荷作用下薄壁錐殼結(jié)構(gòu)的動力學(xué)問題屬于一類具有參數(shù)激勵的陀螺動力系統(tǒng).由于Floquet乘數(shù)假設(shè)不滿足陀螺系統(tǒng)，經(jīng)典Bolotin方法不能用于此類系統(tǒng)的參數(shù)穩(wěn)定性分析[23,24].需借助更為普適的Hill方法[25].

因此，本文擬開展周期軸向力作用下旋轉(zhuǎn)圓錐薄殼動力穩(wěn)定性研究.基于Donnell薄殼理論[26]推導(dǎo)旋轉(zhuǎn)錐殼的動力學(xué)方程，采用GDQ法和Hill法分析系統(tǒng)在周期軸向載荷作用下的參數(shù)不穩(wěn)定性.通過與文獻(xiàn)結(jié)果的對比，驗(yàn)證分析模型的準(zhǔn)確性.研究不同轉(zhuǎn)速和邊界條件時(shí)系統(tǒng)主參數(shù)穩(wěn)定區(qū)的變化，討論軸向力相對變化幅值以及幾何設(shè)計(jì)參數(shù)(錐角、長徑比和厚徑比)對不穩(wěn)定區(qū)的影響規(guī)律.

1 分析模型

1.1 周期軸向力

本文考慮各向同性的截?cái)嘈蛨A錐薄殼，其彈性模量為E，密度為ρ，泊松比為υ.圖1給出了以角速度Ω繞其對稱軸旋轉(zhuǎn)的圓錐薄殼的幾何結(jié)構(gòu)和坐標(biāo)系.其中，α為錐角，L為軸向長度，h為厚度，a和b是錐殼兩端的半徑.錐殼的基準(zhǔn)面取其厚度中面，其中正交坐標(biāo)系(x-θ-z)是固定的，r=r(x)是任意坐標(biāo)點(diǎn)(x,θ,z)的半徑.旋轉(zhuǎn)錐殼在經(jīng)向x、周向θ和法向z方向上的變形分別由u,v,w定義.

圖1 旋轉(zhuǎn)圓錐薄殼的幾何結(jié)構(gòu)和坐標(biāo)系

旋轉(zhuǎn)錐殼在圖1所示的子午方向上承受周期性載荷Na(t)，其形式可假定為恒定載荷上疊加一個(gè)小正弦擾動量，即表示為如下形式

Na(t)=ηaNcr(1+εcosωat)

(1)

式中，ηa和ε分別表示恒定軸向力的幅值和相對變化幅值.Ncr表示當(dāng)錐角α=0時(shí)的圓柱薄殼的屈曲載荷，有Ncr=Eh/[R(2(1-υ2)]1/2，其中R表示柱殼的半徑.本研究中，令R=a.

1.2 旋轉(zhuǎn)薄壁錐殼振動模型

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

(L+LT+LI+LN)u=0

(2)

式中，u=[u,v,w]T.L表示非旋轉(zhuǎn)錐殼的微分算子矩陣，可參考文獻(xiàn)[18].LT表示初始環(huán)向張力引起的微分算子矩陣，可表示為

(3)

式中

(4a)

(4b)

(4c)

(4d)

(4e)

(4f)

式中，r=a+xsinα.LI表示慣性力微分算子矩陣，可表示為

(5)

式中

(6a)

(6b)

LI13=LI31=ρhΩ2sinαcosα

(6c)

(6d)

(6e)

(6f)

LN表示軸向力微分算子矩陣，可表示為

(7)

其中

(8a)

(8b)

后續(xù)分析中，將考慮四種邊界條件：小端固支和大端固支(Cs-Cl)、小端簡支和大端固支(Ss-Cl)、小端固支和大端簡支(Cs-Sl)以及小端簡支和大端簡支(Ss-Sl).

1.3 基于GDQ方法的離散模型

在應(yīng)用GDQ方法之前，需確定圓錐薄殼自由振動的位移場為

(9)

L*u*=0

(10)

其中，u*=[U(x),V(x),W(x)]T是未知的振型空間函數(shù)向量，L*是維度為3×3的微分算子矩陣.GDQ方法的思路[19]：一個(gè)足夠光滑的函數(shù)在一個(gè)離散點(diǎn)上相對于坐標(biāo)方向的導(dǎo)數(shù)可以近似地表示為所有離散點(diǎn)上函數(shù)值的加權(quán)線性和.以函數(shù)U(x)為例，GDQ方法的基本做法如下:

(11)

其中，N是x方向上離散網(wǎng)格點(diǎn)的總數(shù).Cijm是與第m階導(dǎo)數(shù)相關(guān)的加權(quán)系數(shù).選取x方向上余弦函數(shù)點(diǎn)作為離散網(wǎng)格點(diǎn)，即xi={1-cos[(i-1)/(N-1)π]}/2L，(i=1,2,…,N).顯然，越靠近錐殼的端點(diǎn)，離散網(wǎng)格點(diǎn)的分布就越密集.值得注意的是，錐殼兩個(gè)端點(diǎn)的坐標(biāo)為x1=0和xN=L.V(x)和W(x)導(dǎo)數(shù)的離散表達(dá)式可與式(11)中U(x)的表達(dá)式類似.

將這些GDQ方法的表達(dá)式應(yīng)用于常微分控制方程組[式(10)]，得到一組線性代數(shù)方程組:

(12)

(13)

式中U(m)(xi)，V(m)(xi)和W(m)(xi)可分別由式(11)確定，且m=1,2,3,4.

(14)

式中H0,H1,H2,Ha分別表示維度為系數(shù)矩陣J×J的系數(shù)矩陣(J=3N-8).需要說明的是，Ha僅由軸向力作用而引起的.d表示第J階模態(tài)振型列向量，有

(15)

通過求解方程(14)的多項(xiàng)式特征值問題，可以得到給定轉(zhuǎn)速下旋轉(zhuǎn)錐殼的前后行波振動頻率.在運(yùn)動方程的推導(dǎo)中，錐殼的位移場可以視作

(16)

與多項(xiàng)式特征值問題的推導(dǎo)類似，旋轉(zhuǎn)錐殼在周期性軸向載荷作用下的運(yùn)動方程可表示如下:

(17)

式中f=[q1(t),q2(t), …,qJ(t),p1(t),p2(t), …,pJ(t)]T表示維度為2J×1廣義自由度列向量，而質(zhì)量、陀螺和剛度系數(shù)矩陣可表示為

(18)

顯然，在考慮周期性軸向載荷的情況下，圓錐薄殼的控制方程具有Mathieu-Hill型的周期性變化系數(shù).這類系統(tǒng)又稱為參數(shù)激勵系統(tǒng)，系統(tǒng)的參數(shù)穩(wěn)定性是人們關(guān)注的主要問題.

2 參數(shù)穩(wěn)定性分析方法

公式(17)的解的穩(wěn)定性將通過Hill方法[25]進(jìn)行研究.該方法基于Floquet理論，其主要思路是式(17)的解可以寫成指數(shù)部分和周期部分的乘積.用復(fù)Fourier級數(shù)展開表示周期部分，這個(gè)解可以寫成

(19)

f=eλωatΦA(chǔ)n

(20)

(21)

(22)

式中，J1=diag([…kI…])和J2=diag([…k2I…]).通過三角函數(shù)運(yùn)算，cosωatΦ可重新表示為

(23)

式中

(24)

將式(20)～式(24)代入式(17)，應(yīng)用諧波平衡條件，得到如下代數(shù)方程:

(25)

式中，IM=diag([…M…])，IG=diag([…G…])，IK=diag([…K0…])和IK=diag([…Ka…]).為了使式(17)具有式(19)形式的非零解，式(25)的系數(shù)矩陣的行列式需為零(ωa≠0):

(26)

式(26)可用于分析給定軸向力的變化頻率ωa和相對振幅ε的系統(tǒng)穩(wěn)定性.如果系統(tǒng)是穩(wěn)定的，所有特征值λ的實(shí)部為負(fù)，指數(shù)部分隨著時(shí)間的推移而減小.另一方面，如果至少有一個(gè)特征值λ具有正實(shí)部，則系統(tǒng)是不穩(wěn)定的.為了得到穩(wěn)定分析的近似數(shù)值特征值，只需少量的Nk就可以滿足精度要求.

3 分析實(shí)例

3.1 模型驗(yàn)證

在開展穩(wěn)定性研究之前，需要對所提出的薄壁錐殼模型和GDQ方法的結(jié)果進(jìn)行驗(yàn)證.采用ωs表示非旋轉(zhuǎn)錐殼的固有頻率.考慮旋轉(zhuǎn)后，特定模態(tài)(m,n)的前/后行波頻率用ωb和ωf表示，其中，n=1, 2, …表示周向波的數(shù)量，而m=1, 2, …則表示相應(yīng)駐波模式中軸向半波的數(shù)量.以[ρa(bǔ)2(1-υ2)/E]1/2為轉(zhuǎn)速、頻率的量綱，則無量綱化后的轉(zhuǎn)速和行波頻率分別表示為Ω*,ωa*,ωs*,ωb*,ωf*.

表1給出了具有Cs-Cl和Ss-Sl邊界條件的非旋轉(zhuǎn)各向同性圓錐殼的固有頻率ωs*的比較.以Irie等[27]的數(shù)值積分方法得到的結(jié)果為基準(zhǔn)值，所分析的圓錐殼參數(shù)為m=1,υ=0.3,α=30°,h/a=0.01.結(jié)果表明，在增加離散網(wǎng)格點(diǎn)的情況下GDQ方法的結(jié)果很快得到收斂，且與文獻(xiàn)結(jié)果具有較好的一致性.將旋轉(zhuǎn)各向同性圓柱殼在不同轉(zhuǎn)速下的前后行波頻率與Sun等[28]的結(jié)果進(jìn)行比較，如表2所示.結(jié)構(gòu)參數(shù)為υ=0.3,α=0°,h=a/500,L=5a，邊界條件為Ss-Sl.從中也可發(fā)現(xiàn)類似的一致性，表明本研究所建立的分析模型是正確的，用GDQ方法得到的固有頻率結(jié)果是可信的.為了精度要求，在后續(xù)分析中將總網(wǎng)格點(diǎn)數(shù)設(shè)為16.

表1 不同邊界條件時(shí)非旋轉(zhuǎn)圓錐殼的固有頻率比較

表2 旋轉(zhuǎn)圓柱殼在不同轉(zhuǎn)速下的前后行波頻率比較

3.2 不同轉(zhuǎn)速下的參數(shù)不穩(wěn)定區(qū)

分析的圓錐薄殼參數(shù)為m=1,n=2,υ=0.3,α=20°,h/a=0.008,L/a=10，邊界條件為Cs-Cl.計(jì)算了不同轉(zhuǎn)速時(shí)周期軸向力作用下主參數(shù)不穩(wěn)定區(qū)，如圖2所示.考慮拉伸軸向力(ηa=0.2)，且相對變化幅值ε由0.05增加至0.5.可以看出，隨著ε的增加，不穩(wěn)定區(qū)域?qū)ΨQ地?cái)U(kuò)大.與對稱軸相對應(yīng)的激勵頻率即為前后行波頻率之和，即ωa*=ωb*+ωf*.顯然，對于不旋轉(zhuǎn)的錐殼，有ωa*=2ωs*.考慮旋轉(zhuǎn)后，不穩(wěn)定區(qū)域向高頻率范圍整體移動，尤其是當(dāng)轉(zhuǎn)速大于0.002時(shí)，如圖4(c)所示.但在一定的轉(zhuǎn)速范圍內(nèi)，不穩(wěn)定區(qū)域的位置受轉(zhuǎn)速的影響不大.當(dāng)邊界條件改為Ss-Sl、Cs-Sl或Ss-Cl時(shí)，也有類似現(xiàn)象.由于篇幅的限制，這里不給出結(jié)果.在后續(xù)分析中，只考慮Cs-Cl邊界條件，并討論系統(tǒng)參數(shù)對不穩(wěn)定區(qū)域的影響.

圖2 不同轉(zhuǎn)速時(shí)周期軸向力作用下旋轉(zhuǎn)錐殼的主參數(shù)不穩(wěn)定區(qū)

3.3 參數(shù)影響分析

3.3.1 定常軸向力幅值

圖3繪制了具有不同定常軸向力幅值時(shí)對應(yīng)不同模態(tài)的不穩(wěn)定區(qū)域變化.考慮了三種軸向力幅值(ηa=0.1,0.2,0.3)，相對幅值固定為ε=0.5，轉(zhuǎn)速為Ω*=0.005，其他參數(shù)與上一節(jié)模型參數(shù)相同.隨著ηa的增加，不穩(wěn)定區(qū)向高頻區(qū)整體移動.這是由于增加恒定軸向載荷會增加固有頻率值，進(jìn)而使失穩(wěn)區(qū)域的起始點(diǎn)向更高的頻率范圍移動.對于給定的ε，增加恒定軸向載荷ηa也會增加周期軸向載荷的變化幅值，剛度參數(shù)激勵得到加強(qiáng)，使得不穩(wěn)定區(qū)的范圍有所擴(kuò)大，如圖3(a)所示.隨著周向波數(shù)的增加[模態(tài)(1,3)和模態(tài)(1,4)，如圖3(b)和圖3(c)所示]，不穩(wěn)定區(qū)的變化規(guī)律與模態(tài)(1,1)類似.

圖3 定常軸向力幅值對不穩(wěn)定區(qū)的影響

3.3.2 錐角

本節(jié)討論錐角對不穩(wěn)定區(qū)的影響.分析的圓錐薄殼參數(shù)為υ=0.3,h/a=0.008,L/a=10, Ω*=0.005,ηa=0.2,ε=0.5.考慮三種錐角值(α=20°, 40°, 60°)，結(jié)果如圖4所示.對于模態(tài)(1,1)，增大α值，失穩(wěn)區(qū)整體向低頻段移動，且不穩(wěn)定區(qū)范圍有所擴(kuò)大.對于模態(tài)(1,3)和模態(tài)(1,4)，首先使不穩(wěn)定區(qū)域向高頻段移動而后向低頻段移動，如圖4(b)和圖4(c)所示.不穩(wěn)定區(qū)范圍也有所增加.與模態(tài)(1,1)的情況相比，模態(tài)(1,3)和模態(tài)(1,4)的不穩(wěn)定區(qū)的移動和范圍擴(kuò)大量相對較小.這表明，隨著周向波的增加，錐角對失穩(wěn)區(qū)的影響呈減弱趨勢.這一現(xiàn)象也與Ng等[9]的發(fā)現(xiàn)一致.

圖4 錐角對不穩(wěn)定區(qū)的影響

3.3.3 厚徑比

本節(jié)研究厚徑比對不同模態(tài)的不穩(wěn)定區(qū)影響，如圖5所示.分析的圓錐薄殼參數(shù)為υ=0.3,α=20°,L/a=10, Ω*=0.005,ηa=0.2,ε=0.5.分析了三種厚徑比值，即h/a=0.005, 0.008, 0.012.結(jié)果表明，增大h/a值不僅會使不穩(wěn)定區(qū)向高頻段移動，而且會使不穩(wěn)定區(qū)明顯擴(kuò)大.對于周向波較大的情況，如圖5(c)所示的模態(tài)(1,4)，不穩(wěn)定區(qū)隨h/a的移動和擴(kuò)大并未減弱.

圖5 厚徑比對不穩(wěn)定區(qū)的影響

3.3.4 長徑比

討論長徑比對不穩(wěn)定區(qū)的影響.分析的圓錐薄殼參數(shù)為υ=0.3,α=20°,h/a=0.008,Ω*=0.005,ηa=0.2,ε=0.5.當(dāng)L/a=10, 15和20時(shí)，不穩(wěn)定區(qū)的結(jié)果如圖6所示.可以發(fā)現(xiàn)，隨著L/a的增加，不穩(wěn)定區(qū)主要出現(xiàn)在低頻段.對于周向波數(shù)較少的失穩(wěn)模態(tài)，如模態(tài)(1,1)，除失穩(wěn)區(qū)的整體移動外，失穩(wěn)區(qū)寬度的變化不大.但是，對于周向波數(shù)較多的失穩(wěn)模態(tài)，增大L/a值也會減小失穩(wěn)區(qū)寬度，如圖6(c)所示.

圖6 長徑比對不穩(wěn)定區(qū)的影響

4 結(jié)論

本文首先推導(dǎo)了旋轉(zhuǎn)錐殼的動力學(xué)方程，采用GDQ法和Hill法分析了系統(tǒng)在周期軸向載荷作用下的參數(shù)不穩(wěn)定性.計(jì)算并討論了多個(gè)不穩(wěn)定區(qū)隨工況和幾何參數(shù)的變化規(guī)律.結(jié)果表明：提高轉(zhuǎn)速會導(dǎo)致不穩(wěn)定區(qū)沿頻率軸移動，但對不穩(wěn)定寬度影響不大.增加恒定拉伸軸向載荷，不僅會顯著增加失穩(wěn)寬度，而且會導(dǎo)致失穩(wěn)區(qū)域向更高的頻率范圍移動.錐角、厚徑比或長徑比的變化都會導(dǎo)致不穩(wěn)定區(qū)沿頻率軸移動.錐角和厚徑比會增大失穩(wěn)寬度(長徑比會減小).隨著周向波個(gè)數(shù)的增加，錐角對失穩(wěn)區(qū)的影響逐漸減弱，而厚徑比的影響則基本保持不變.

上述研究基于線性和各向同性薄壁殼體假設(shè).若殼體是非線性的且各向異性，例如考慮大變形、功能梯度或復(fù)合材料層合殼材料.進(jìn)一步的工作應(yīng)重點(diǎn)探討非線性因素的引入對系統(tǒng)動力穩(wěn)定性的影響.