陳愛志，王浩然，柳貢民

（1.海裝沈陽局駐哈爾濱汽輪機廠有限責任公司 軍事代表室，哈爾濱 150046；2.耐世特汽車系統(tǒng)有限公司，江蘇 蘇州 215000；3.哈爾濱工程大學 動力與能源工程學院，哈爾濱 150001）

圓柱殼結構常見于船舶、航空航天、水利水電等工程領域。結構振動直接影響圓柱殼結構的可靠性。圓柱殼結構固有特性的計算分析是結構優(yōu)化、安全性設計中的重要環(huán)節(jié)。因而，準確、高效地完成圓柱殼結構的固有特性計算就成為人們關注的重要研究內容。

結構固有特性的分析方法有很多，包括當今應用廣泛的有限元法（FEM）[1]，能高效計算鏈式結構的傳遞矩陣法（TMM）[2]，將二者結合起來的有限元傳遞矩陣法（FETMM）[3]以及常用于瞬態(tài)應力分析的特征線方法（MOC）[4]等。其中，F(xiàn)ETMM是由Dokanish[5]在1972年提出的，該方法將FEM和TMM聯(lián)合起來進行結構的分析，拓寬了傳統(tǒng)TMM的研究范圍并改善了FEM在使用過程中計算效率較低的問題。FETMM處理復雜結構時既具有FEM的優(yōu)點，又兼?zhèn)銽MM的高效率[6]。徐銘陶[7]創(chuàng)建了一種新的始端線和終端線節(jié)點數(shù)目不等的傳遞元件，提出了廣義傳遞矩陣和有限單元組合法，此法便于分析復雜形狀平板的振動。何斌、芮筱亭等[8]基于歐拉梁理論，用FETMM對細長彈箭的振動特性進行了計算。當前，對FETMM的利用局限于對梁和板的分析，較少出現(xiàn)在圓柱殼的結構振動特性的分析中。不過根據(jù)Petyt[9]的研究，可以通過對圓柱殼周向模態(tài)的限制，利用類梁的殼單元來減少自由度，來達到簡化計算的目的。張勇亮[10]利用此理論對內部含有流體的圓柱殼耦合振動波傳播問題進行了研究。

本文利用圓柱殼的軸對稱特性，基于Sanders薄殼理論[11]建立類梁的殼模型，運用MATLAB軟件完成了有限元傳遞矩陣法程序的編寫和圓柱殼固有特性的計算，驗證了有限元傳遞矩陣法與類梁的殼模型結合運用于輸流圓柱殼固有特性研究的準確性和高效性。

1 基于類梁殼模型的圓柱殼流固耦合方程

圖1 含流圓柱殼模型

參考文獻[12] ，將圓柱殼3個方向的位移進行傅里葉余弦展開,具體表達式如下:

其中，n為周向波數(shù)，ux為軸向位移，uθ為徑向位移，ur為周向位移。軸向和徑向位移的模態(tài)假設為cosnθ，周向假設為sinnθ[9]。將位移進行如此變換，則只需進行軸向網格劃分而不必進行周向網格劃分，如此，有限元方程的自由度得到降低，計算量和計算時間大大減少。這樣劃分所得單元稱為類梁的殼單元。

圖2分別為基于殼單元的一般圓柱殼網格劃分模型與基于本文類梁殼單元的圓柱殼網格劃分模型。

圖2 兩種殼體模型

在求解過程中，每個節(jié)點各方向位移的具體形式可以通過相應的形函數(shù)和節(jié)點位移呈現(xiàn)[12]

因此，單元位移表達式可以簡單地寫為

[Ns] 為形函數(shù)，{uˉ}為節(jié)點位移矢量。每個單元有2個節(jié)點，每個節(jié)點包含4個位移，即軸向uˉx、周向uˉr、徑向uˉθ和一個旋轉量φ=?uˉr/?x，節(jié)點位移矢量和形函數(shù)具體形式如下

其中：

ξ=x∕le，le是單元的長度。此處，圓柱殼的徑向位移類似梁的橫向位移。因此，方程式(4)中的形函數(shù)和梁橫向振動的形函數(shù)相同。

參考文獻[12] ，圓柱殼的單元質量、剛度矩陣和力矢量可以分別表示為

其中：ρs為殼體密度，h為殼體厚度，R為殼體中面半徑，[B ] 為應變位移矩陣，[D] 為應力矩陣，pd為流體作用于殼壁的壓力，ps為外部分布力。

另外，考慮管內流體為不可壓縮、非黏性流體時，極坐標系下的拉普拉斯方程可以表示為[13]

對于管內以恒定速度流動的流體，其單元質量、剛度矩陣和耦合阻尼矩陣可表示為[12]

hf表示流體的有效厚度，它取決于周向模數(shù)n和特征值λ。

流體作用于圓柱殼內壁的壓力可以用勢函數(shù)來表達，其中U為流體速度。

2 有限元傳遞矩陣法

將輸流圓柱殼整體沿軸向離散成若干單元，則單元的流固耦合振動微分方程可以表示為[10]

假設圓柱殼振動為簡諧振動形式，則單元位移列向量以及單元力向量可表示為

ω為管道振動圓頻率。因此，對式(12)進行變換得

式(15)可以簡寫成如下形式

其中

分塊以后得

對式(16)進行恒等變換，寫成如下形式

其中

由此得圓柱殼軸向劃分為n個單元時，前后端的狀態(tài)矢量傳遞關系滿足

從輸入端矢量到輸出端矢量的整體傳遞關系可用傳遞矩陣Ut表示

值得注意的是，在運用有限元傳遞矩陣法求解圓柱殼振動問題時，由于矩陣的連乘，會出現(xiàn)計算不穩(wěn)定的情況，導致結果出現(xiàn)較大誤差。本文引入Riccati變換后，可以較好地解決這一問題。圖3給出了引入Riccati變換前后的計算結果對比。

圖3 Riccati傳遞矩陣法與傳統(tǒng)傳遞矩陣法對比

其中，縱坐標Log(miu)為編程求解過程中所設的參數(shù)，其峰值所對應的橫坐標頻率即為計算所得的圓柱殼固有頻率。

3 計算精度影響因素分析及算例驗證

3.1 計算精度影響因素分析

如前所述，將類梁殼模型與Riccati有限元傳遞矩陣法結合進行圓柱殼固有頻率求解。作為對比算例[13]的薄壁圓柱殼與流體具體參數(shù)如表1所示，邊界條件為一端固支一端自由。

表1 輸流圓柱殼與流體參數(shù)

在求解數(shù)學問題時，選擇恰當?shù)姆椒ㄊ琼樌蠼獾年P鍵。本文運用3、4、5次高斯積分分別進行計算，但是在求解過程中發(fā)現(xiàn)運用3次高斯積分所得結果并不理想，具體結果對比如圖4和表2所示。

從圖4和表2的計算結果很容易發(fā)現(xiàn)，相對于3次高斯積分，4次高級積分和5次高斯積分得到的計算結果更接近于實驗值和商業(yè)軟件仿真結果。因此，在接下來的計算分析過程中采用計算結果更精確的5次高斯積分進行求解。

表2 不同高斯積分階數(shù)計算結果對比

簡便和高效是本文方法的突出特點，但在求解過程中發(fā)現(xiàn)軸向劃分的單元數(shù)目會影響算法的準確性和效率。圖5展示了同一模態(tài)下劃分不同軸向單元數(shù)時圓柱殼固有頻率計算結果對比，表3列出不同單元數(shù)下固有頻率的同時，也列出了計算所耗費的時間。

可以看出單元數(shù)目太低會導致計算結果準確性下降，網格數(shù)增加雖能提高計算精度，卻又會導致運算時間增加。

表3 不同單元數(shù)下的固有頻率計算結果與耗時對比

3.2 算例驗證

在3.1所選定的高斯積分階數(shù)(5階)及軸向單元劃分數(shù)(20個)的基礎之上，利用本文方法對圓柱殼的固有特性進行計算。此外，采用商業(yè)軟件ANSYS對同一圓柱殼振動固有頻率進行同步求解，仿真計算中選取的單元類型為SHELL63。

圖4 不同高斯積分階數(shù)計算結果對比

圖5 不同單元數(shù)下的固有頻率計算結果對比

圖6所示為計算得到的不考慮內含流體時不同周向模態(tài)下的圓柱殼振動固有頻率；表4、圖7表現(xiàn)的是實驗值、有限元軟件所得數(shù)據(jù)與本文方法計算結果的對比，從中可以看出本文模型和方法的準確性。

圖6 不同周向模態(tài)下的圓柱殼固有頻率

圖7 固有頻率對比

表4 不含流圓柱殼的固有頻率

4 流體對圓柱殼固有特性的影響

分析圓柱殼固有特性與其內含流體是否有著密切的聯(lián)系。

流體對圓柱殼的作用包括附加質量、阻尼和剛度。其中，作為附加質量的流體質量是變化的，這種變化和周向模態(tài)有關。

圖8表示的是不同周向模態(tài)下流體作為附加質量作用于圓柱殼時的有效厚度值。

圖8 流體在不同周向模態(tài)下的附加質量有效厚度

從圖8中可以發(fā)現(xiàn)當n=1時，有效厚度值為1即流體完全作用于附加質量。隨著n的增加,流體作用于圓柱殼的有效厚度逐漸變小，變化也是由快到慢的。

對同一圓柱殼模型，內含靜止流體時固有頻率相對不含流體時的固有頻率值會降低，這是因為內含的流體相當于增加了圓柱殼的質量，而圓柱殼固有頻率會隨著其質量的增加而減小。含靜止流體圓柱殼固有頻率計算結果與實驗值以及有限元值的對比見表5和圖9。

表5 含靜止流體圓柱殼固有頻率

圖9 含靜水固有頻率對比

圖10顯示了n=1時含流與不含流圓柱殼固有頻率計算結果對比，可以看出含流時同一頻率區(qū)間內固有頻率(Log(miu)的峰值)個數(shù)明顯增加。

圖10 n=1時固有頻率對比

這是因為這些固有頻率的計算結果之中不僅包含了含流圓柱殼整體振動的固有頻率，還包含了內含的水柱的固有頻率。

接下來考慮圓柱殼內含的流動流體對其固有特性的影響。采用的算例來自于Weaver[12]，具體參數(shù)如下：ρf∕ρt=0.128，h∕R=0.01，L∕R=2。無量綱流速Uˉ =U∕U0，無量綱頻率為ωˉ=ω∕ω0。

本文得到的結果和Weaver以及Selmane[13]得到的結果對比如圖11所示?？梢院苊黠@看出圓柱殼固有頻率隨著流速的增加而降低。同時從對比的結果可以發(fā)現(xiàn)，流速較低時，本文方法和Weaver以及Selman得到的結果很相近，但隨著流速的增加，對比周向模態(tài)數(shù)m=2時的結果可以發(fā)現(xiàn)，本文方法與Weaver出現(xiàn)較大不同，但與Selman文中的結果吻合相對較好。這可能是由于Weaver在應用伽遼金法時采用的項較少。總體來看，本文方法在高流速時計算結果的精確度沒有流速較低時好，這可能與本文處理圓柱殼內流體的方式有關。

圖11 無量綱頻率計算結果對比

5 結語

本文結合類梁殼模型及有限元理論，考慮圓柱殼與流體的耦合作用，采用有限元傳遞矩陣法，對輸流圓柱殼的固有振動特性進行了分析。擴展了傳遞矩陣法的研究對象，同時在保持有限元計算精度的前提下，提高了計算速度與效率。通過與實驗和仿真數(shù)據(jù)的對比驗證了本文模型和算法的準確性，同時研究了高斯積分次數(shù)以及劃分的軸向單元個數(shù)對計算精度的影響，以及流體對圓柱殼振動固有頻率的影響

計算結果表明隨著高斯積分次數(shù)和軸向劃分單元數(shù)的增加，計算精度也會隨之增加，但增加的幅度會越來越小，且計算時間也會相應延長，因此在滿足一定計算精度要求的前提，應使2個參數(shù)盡量的小以提高計算效率。

含流與不含流圓柱殼對比，同一頻率區(qū)間內固有頻率個數(shù)明顯增加，同時圓柱殼振動固有頻率隨著流速增加而降低。