輸流圓柱殼固有特性分析的有限元傳遞矩陣法

陳愛志，王浩然，柳貢民

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

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

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

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

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

圖1 含流圓柱殼模型

參考文獻(xiàn)[12] ，將圓柱殼3個(gè)方向的位移進(jìn)行傅里葉余弦展開,具體表達(dá)式如下:

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

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

圖2 兩種殼體模型

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

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

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

其中：

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

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

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

另外，考慮管內(nèi)流體為不可壓縮、非黏性流體時(shí)，極坐標(biāo)系下的拉普拉斯方程可以表示為[13]

對(duì)于管內(nèi)以恒定速度流動(dòng)的流體，其單元質(zhì)量、剛度矩陣和耦合阻尼矩陣可表示為[12]

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

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

2 有限元傳遞矩陣法

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

假設(shè)圓柱殼振動(dòng)為簡諧振動(dòng)形式，則單元位移列向量以及單元力向量可表示為

ω為管道振動(dòng)圓頻率。因此，對(duì)式(12)進(jìn)行變換得

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

其中

分塊以后得

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

其中

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

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

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

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

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

3 計(jì)算精度影響因素分析及算例驗(yàn)證

3.1 計(jì)算精度影響因素分析

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

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

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

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

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

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

可以看出單元數(shù)目太低會(huì)導(dǎo)致計(jì)算結(jié)果準(zhǔn)確性下降，網(wǎng)格數(shù)增加雖能提高計(jì)算精度，卻又會(huì)導(dǎo)致運(yùn)算時(shí)間增加。

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

3.2 算例驗(yàn)證

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

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

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

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

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

圖7 固有頻率對(duì)比

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

4 流體對(duì)圓柱殼固有特性的影響

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

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

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

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

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

對(duì)同一圓柱殼模型，內(nèi)含靜止流體時(shí)固有頻率相對(duì)不含流體時(shí)的固有頻率值會(huì)降低，這是因?yàn)閮?nèi)含的流體相當(dāng)于增加了圓柱殼的質(zhì)量，而圓柱殼固有頻率會(huì)隨著其質(zhì)量的增加而減小。含靜止流體圓柱殼固有頻率計(jì)算結(jié)果與實(shí)驗(yàn)值以及有限元值的對(duì)比見表5和圖9。

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

圖9 含靜水固有頻率對(duì)比

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

圖10 n=1時(shí)固有頻率對(duì)比

這是因?yàn)檫@些固有頻率的計(jì)算結(jié)果之中不僅包含了含流圓柱殼整體振動(dòng)的固有頻率，還包含了內(nèi)含的水柱的固有頻率。

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

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

圖11 無量綱頻率計(jì)算結(jié)果對(duì)比

5 結(jié)語

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

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

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