多自由度參激系統(tǒng)穩(wěn)定性分析的數(shù)值解法

(哈爾濱工業(yè)大學(xué) 航天學(xué)院，哈爾濱 150001)

1 引言

工程上許多結(jié)構(gòu)受周期基礎(chǔ)激勵(lì)、周期軸向力和周期變剛度、變阻尼等因素的影響而具有參激振動(dòng)特性，因此，以Mathieu方程為代表的周期系數(shù)參激系統(tǒng)動(dòng)力穩(wěn)定問題是具有重要理論意義和工程價(jià)值的研究課題[1,2]?，F(xiàn)有研究通常采用Floquet方法[3,4]、Bolotin方法[5]或Hsu方法[6]確定系統(tǒng)的穩(wěn)定性。在多自由度參激系統(tǒng)中，Bolotin方法可繞開求解Floquet指數(shù)，直接快速獲得系統(tǒng)的簡單不穩(wěn)定區(qū)域，但無法得到組合不穩(wěn)定區(qū)域[7,8]。Hsu方法可用于確定參激系統(tǒng)的簡單不穩(wěn)定區(qū)域和組合不穩(wěn)定區(qū)域，但該方法需對周期進(jìn)行離散處理，計(jì)算過程較為繁復(fù)，不利于多自由度系統(tǒng)的計(jì)算。Floquet方法非常適合求解周期時(shí)變系數(shù)系統(tǒng)周期解，通過周期分量與指數(shù)特征分量之積的設(shè)解形式導(dǎo)出特征行列式，適用于求解各種不穩(wěn)定區(qū)域。Pei[9]在求解陀螺系統(tǒng)時(shí)指出，Bolotin方法相對于Floquet方法會(huì)擴(kuò)大動(dòng)力失穩(wěn)區(qū)，并以Floquet乘子穿越單位圓形式解釋了成因。Dai等[10]采用哈爾小波和Floquet方法分析受軸向周期激勵(lì)圓錐殼的動(dòng)力穩(wěn)定區(qū)域，討論一次近似和二次近似設(shè)解下主不穩(wěn)定區(qū)域的差異。兩者的求解思想都是根據(jù)Floquet理論將其穩(wěn)定性問題轉(zhuǎn)化為求解行列式特征值問題。Takahashi[11]通過對含F(xiàn)loquet指數(shù)的特征值行列式求解，獲得了多自由度系統(tǒng)的簡單不穩(wěn)定區(qū)域和組合不穩(wěn)定區(qū)域，限于當(dāng)時(shí)計(jì)算能力，僅分析了三自由度參激系統(tǒng)。Lee[12]將其推廣到多自由度周期參激系統(tǒng)。應(yīng)祖光等[13]基于狀態(tài)空間分析、Floquet理論與特征值分析法[14],發(fā)展了多自由度周期參激系統(tǒng)穩(wěn)定性的數(shù)值直接法，并用于研究斜拉索的不穩(wěn)定區(qū)。隨著現(xiàn)代計(jì)算機(jī)計(jì)算能力及效率的大幅提高，使得更高自由度系統(tǒng)的計(jì)算得以實(shí)現(xiàn)，適用于各種非線性周期系數(shù)系統(tǒng)的Floquet方法也得到更為廣泛的運(yùn)用和推廣。

Bolotin法和Floquet法是學(xué)者們研究參激不穩(wěn)定的主要方法，但現(xiàn)有研究的關(guān)注對象主要是主不穩(wěn)定區(qū)域。本文將Floquet法與類似于Bolotin周期數(shù)設(shè)解方式相結(jié)合，得到適用于多自由度參激系統(tǒng)穩(wěn)定性的數(shù)值方法，可求解分析主不穩(wěn)定區(qū)域和組合不穩(wěn)定區(qū)域，以期為參激系統(tǒng)動(dòng)力失穩(wěn)問題提供求解手段。

2 穩(wěn)定性特征值問題及物理模型

一般的周期剛度動(dòng)力學(xué)擾動(dòng)方程為[2]

(1)

式中M,C,G,K和S分別為質(zhì)量陣、阻尼陣、陀螺陣、剛度陣和剛度變化矩陣，皆為對稱正定的常量實(shí)矩陣，α和β為周期時(shí)變剛度矩陣控制參數(shù)，ω為周期時(shí)變頻率。

基于Floquet理論將擾動(dòng)解x(t)表示成指數(shù)特征分量與周期分量之積，并展開成Fourier級(jí)數(shù)的形式，有

(2)

式中λ為Floquet指數(shù)，[b]0,[b]k和[a]k為系數(shù)列向量。

將式(2)代入式(1)，通過諧波平衡可導(dǎo)出無限個(gè)代數(shù)方程，寫為

{λ2E2+λE1+E0}X=0

(3)

式中X={[b]k,[b]0,[a]k}和Ei(i=1,2,3)為可確定的矩陣。

構(gòu)建如下特征值問題。

|D-λI|=0

(4)

(5)

通過求解矩陣D的特征值可直接確定擾動(dòng)解x(t)的穩(wěn)定性，若特征值實(shí)部均為負(fù)值，則方程有有限解，系統(tǒng)是穩(wěn)定的；當(dāng)特征值實(shí)部出現(xiàn)一個(gè)正值時(shí)，方程解發(fā)散，系統(tǒng)不穩(wěn)定。

由Bolotin方法可知，在不穩(wěn)定區(qū)域邊界上微分方程組具有周期1和周期2的周期解，其對應(yīng)設(shè)解形式分別為[10]

(6)

(7)

本文采用類似Bolotin方法的設(shè)解形式，將設(shè)解分為類周期1和類周期2兩類設(shè)解[15]。

x=eλ t[b]0+

(8)

(9)

進(jìn)而得到響應(yīng)的D矩陣，具體形式見附錄。

3 數(shù)值算例

3.1 物理模型

兩端簡支旋轉(zhuǎn)軸受周期性軸向力f(t)=fS+fDcos(ωt)作用時(shí)，其轉(zhuǎn)速為Ω，如圖1所示。其中，固定坐標(biāo)系(XYZ)和旋轉(zhuǎn)坐標(biāo)系(xyz)繞軸線轉(zhuǎn)動(dòng)，x軸與X軸重合。根據(jù)Euler梁理論，忽略慣性矩和剪切效應(yīng)，采用假設(shè)模態(tài)離散化，并進(jìn)行無量綱化處理，可得到方程(1)的形式[16]，其中α=0，Ω=0.5。

圖1 受周期載荷的旋轉(zhuǎn)軸系

Fig.1 Rotating shaft under periodic loads

3.2 設(shè)解近似項(xiàng)數(shù)的影響

圖2為不同設(shè)解項(xiàng)數(shù)下得到的無阻尼參激系統(tǒng)動(dòng)力失穩(wěn)結(jié)果，在β-ω平面域內(nèi)劃分為穩(wěn)定和不穩(wěn)定兩個(gè)區(qū)域，且隨著設(shè)解近似項(xiàng)數(shù)的增加，主不穩(wěn)定區(qū)域幾乎不變，但出現(xiàn)了多個(gè)次不穩(wěn)定區(qū)域，主要為頻率ω1+ω2,(ω1+ω2)/2,(ω1+ω2)/3的不穩(wěn)定區(qū)域。在參數(shù)域上取點(diǎn)，采用 Runge -Kutta 法進(jìn)行數(shù)值驗(yàn)證，如圖3所示。在主不穩(wěn)定區(qū)域的邊界附近，點(diǎn)A處于穩(wěn)定區(qū)域，其響應(yīng)是穩(wěn)定的，而取跨過邊界到達(dá)不穩(wěn)定區(qū)域的點(diǎn)B，其響應(yīng)是發(fā)散的，說明該主不穩(wěn)定區(qū)域的邊界精度較高。取二階不穩(wěn)定區(qū)域邊界上點(diǎn)C和點(diǎn)D分析發(fā)現(xiàn)，當(dāng)p=2時(shí)，C和D兩點(diǎn)均穩(wěn)定；當(dāng)p=3時(shí)，點(diǎn)C穩(wěn)定，點(diǎn)D不穩(wěn)定；數(shù)值驗(yàn)證結(jié)果符合p=3時(shí)情況。通過更高階不穩(wěn)定區(qū)域的分析發(fā)現(xiàn)，精確確定n階不穩(wěn)定區(qū)域要求設(shè)解近似項(xiàng)數(shù)p≥n。

圖2 不同設(shè)解項(xiàng)數(shù)下Floquet法獲得的動(dòng)力失穩(wěn)

Fig.2 Dynamic instability obtained by Floquet method with different solution terms

3.3 設(shè)解周期數(shù)的影響

為了比較本文方法與Bolotin方法，對方程參數(shù)做如下修改。

圖4所示為Bolotin法周期1和周期2設(shè)解(3次近似)以及Floquet法類周期1和類周期2設(shè)解(3次近似)下得到的動(dòng)力不穩(wěn)定區(qū)域?？梢?，Bolotin方法確定的不穩(wěn)定區(qū)域是在不同周期數(shù)(1和2)設(shè)解下交替得到的，而本文方法可以得到1階和2階簡單不穩(wěn)定區(qū)域和組合不穩(wěn)定區(qū)域，這兩種類周期數(shù)設(shè)解獲得的不穩(wěn)定區(qū)域幾乎完全重合，并且隨著近似次數(shù)p的提高，前n(n

圖3 時(shí)間歷程

Fig.3 Time history of the system

圖4 Bolotin法與本文Floquet法在不同周期數(shù)設(shè)解下的動(dòng)力失穩(wěn)對比

Fig.4 Comparison of dynamic instability obtained by Floquet method and Bolotin method with different periodic terms

3.4 阻尼對系統(tǒng)的影響

圖5給出了系統(tǒng)在不同阻尼比下的動(dòng)力失穩(wěn)結(jié)果，其中曲線上方的區(qū)域?yàn)椴环€(wěn)定區(qū)域。

可以看出，不同阻尼比下獲得的主不穩(wěn)定區(qū)域與文獻(xiàn)[16]的數(shù)值結(jié)果吻合。隨著阻尼比的增大，穩(wěn)定區(qū)域整體呈擴(kuò)大趨勢，不穩(wěn)定區(qū)域邊界的變化趨勢像一條彎曲的繩子緩慢拉直，各個(gè)次不穩(wěn)定區(qū)域縮小，邊界更加平滑。這也是在大部分有阻尼參激振動(dòng)系統(tǒng)中，學(xué)者們對次不穩(wěn)定區(qū)域研究少或略去不研究的原因之一。但是，本文的研究表明，在小阻尼情況下，次不穩(wěn)定區(qū)域的研究還是非常有必要的。

圖5 不同阻尼下的動(dòng)力失穩(wěn)

Fig.5 Dynamic instability with different damping parameters

4 結(jié) 論

本文借助特征值分析法和Floquet方法，通過擾動(dòng)解的存在條件得到特征值方程，進(jìn)而得到適用于多自由度參激系統(tǒng)穩(wěn)定性的數(shù)值求解方法。通過對一個(gè)受周期常激勵(lì)的兩自由度旋轉(zhuǎn)軸系統(tǒng)進(jìn)行動(dòng)力穩(wěn)定性分析，并與Bolotin方法得到的動(dòng)力失穩(wěn)圖對比，得出如下結(jié)論。

(1) 高階動(dòng)力不穩(wěn)定區(qū)域的確定依賴于設(shè)解近似項(xiàng)數(shù)，設(shè)解項(xiàng)數(shù)越多，對次不穩(wěn)定區(qū)域影響越大，精確確定高階不穩(wěn)定區(qū)域要求設(shè)解近似項(xiàng)數(shù)大于該階不穩(wěn)定區(qū)域的階數(shù)。

(2) 類比于Bolotin法，通過改進(jìn)質(zhì)量矩陣的取值來改進(jìn)Floquet法，按類周期1和類周期2設(shè)解均能得到完整的不穩(wěn)定區(qū)域，并且能夠極大地減少計(jì)算量，加快求解速度。

(3) 阻尼將在整體上縮小不穩(wěn)定區(qū)域，使動(dòng)力失穩(wěn)圖中的不穩(wěn)定區(qū)域邊界趨于平滑，并在整體上呈現(xiàn)出將曲折的不穩(wěn)定區(qū)域邊界拉直的效果。

本文的研究為一般多自由度周期參激阻尼系統(tǒng)的參激穩(wěn)定性分析提供了簡潔的直接數(shù)值解法。

附 錄：

類周期1設(shè)解，如式(8)，可以得到p次近似下(k=2p)，Ei(i=0,1,2)的具體矩陣形式:

類周期2設(shè)解時(shí)，如式(9)，可以得到p次近似下(k=2p-1)，Ei(i=0,1,2)的具體矩陣形式