田 楊，章橋新，余金桂，黃 濤

(1.武漢理工大學(xué) 機(jī)電工程學(xué)院，湖北 武漢 430070；2.湖北工業(yè)大學(xué) 土木建筑與環(huán)境學(xué)院，湖北 武漢 430068)

0 引 言

圓柱殼結(jié)構(gòu)被廣泛運(yùn)用于建筑、航空航天和船舶等領(lǐng)域，加肋圓柱殼更是潛艇的主要結(jié)構(gòu)形式。李學(xué)斌[1]利用波動法討論圓柱殼在不同的邊界條件下固有頻率和模態(tài)振型的計(jì)算精度，結(jié)果表明，當(dāng)波動法運(yùn)用在兩端簡支的細(xì)長圓柱殼時，精度較高。江晨半等[2]基于精細(xì)積分和傳遞矩陣法對變厚度圓柱殼進(jìn)行振動分析，與有限元仿真結(jié)果進(jìn)行比較驗(yàn)證準(zhǔn)確性，并討論邊界條件和結(jié)構(gòu)尺寸對振動特性的影響。何偉東等[3]基于虛擬質(zhì)量法和邊界元法研究水下有限長圓柱殼的聲輻射，分析不同深度的輻射聲壓值。何春雨等[4]利用有限元法通過數(shù)值仿真研究海水管路系統(tǒng)的振動特性。

傳遞矩陣法是一種分析鏈?zhǔn)浇Y(jié)構(gòu)的常用方法，其優(yōu)點(diǎn)是計(jì)算過程清晰明了，方法簡便，結(jié)果精度高，只需將結(jié)構(gòu)寫成一階微分方程形式即可使用傳遞矩陣法。國外較早使用傳遞矩陣法計(jì)算圓柱殼振動特性的是日本學(xué)者IRIE等[5]，除研究雙層圓柱殼外，還研究變厚度圓錐殼[6]等，但在求解結(jié)果時采用冪級數(shù)展開取逼近，導(dǎo)致結(jié)果精度不高。張娟等[7]基于傳遞矩陣法研究漁船推進(jìn)軸系的振動特性。向宇[8]運(yùn)用微分方程和矩陣?yán)碚撏茖?dǎo)傳遞矩陣的精確形式。王祖華等[9]使用傳遞矩陣的精確形式求解短粗環(huán)肋圓柱殼的振動特性，但在考慮環(huán)肋與殼體的相互作用時，忽略面外彎曲和扭轉(zhuǎn)振動。

基于傳遞矩陣法計(jì)算縱肋圓柱殼的振動特性，按照傳遞方向?qū)鬟f矩陣法分為周向傳遞矩陣法和軸向傳遞矩陣法。在考慮環(huán)肋與殼體相互作用時將面內(nèi)和面外作用力全部加上，推導(dǎo)縱肋圓柱殼的場傳遞矩陣和點(diǎn)傳遞矩陣，采用精細(xì)時程積分法求解，保證計(jì)算過程的精度，并討論兩種傳遞矩陣法適用條件及范圍。

1 基本理論公式

基于Flugge殼體理論，殼體的應(yīng)力和應(yīng)變的等式按照位移的形式，內(nèi)力和內(nèi)力矩可表達(dá)為

(1)

式中：Nx和Nθ分別為面內(nèi)的軸向力和周向力；Nxθ和Nθx分別為不同方向的剪切力；Mx和Mθ分別為繞周向和軸向的彎矩；Mxθ和Mθx分別為不同方向的扭矩；K為彎曲剛度，K=Eh/(1-μ2)，其中，E為彈性模量，h為厚度，μ為泊松比；D為薄膜剛度，D=Eh3/12(1-μ2)；u為軸向位移；x為軸向長度；v為周向位移；R為半徑；θ為周向角；w為徑向位移；Qx為橫向剪力。

徑向位移與轉(zhuǎn)角ψ之間的關(guān)系為

(2)

柱坐標(biāo)系的圓柱殼段在自由振動時由慣性力引起的分布載荷為

(3)

式中：px、pθ和pr分別為軸向、周向和徑向的慣性力，其中x、θ和r分別為圓柱殼的軸向、周向和徑向；ρ為材料密度；t為時間；ω為自由振動頻率；時間因子省略。

圓柱殼的振動微分控制方程表示為

(4)

2 縱肋圓柱殼自由振動特性分析

圖1為縱肋圓柱殼結(jié)構(gòu)示例。圖1中，殼體長為L，半徑為R，且該圓柱殼滿足殼體理論中的基本假設(shè)。

圖1 縱肋圓柱殼結(jié)構(gòu)示例

2.1 基于周向傳遞的場傳遞矩陣

周向傳遞矩陣法與軸向傳遞矩陣法最大的區(qū)別在于狀態(tài)矢量元素傳遞的方向不同，當(dāng)需要狀態(tài)矢量元素沿周向傳遞時，殼體位移分量和內(nèi)力展開成三角函數(shù)形式與軸向時有所不同，對于圓柱殼邊界條件為兩端簡支、殼體自由振動時，殼體位移分量和內(nèi)力可以寫成如下形式：

(5)

式中：m為軸向半波數(shù)；標(biāo)有上劃線的量為無量綱變量；Qθ為周向剪力；Sθ為Kelvin-Kirchhoff剪力。

為了方便計(jì)算，同樣引入如下無量綱參數(shù)：

(6)

式中：l和H分別為長度和厚度對應(yīng)的無量綱參數(shù)；λ為頻率因子。

在16個狀態(tài)矢量元素中只保留8個，消去Qx、Q、N、Nθ、Mx和Mxθ，柱殼振動方程可寫成矩陣微分方程：

(7)

式中：SYM表示對稱矩陣。

(8)

解該微分方程的精度決定最終解的精度，為了得到高精度的結(jié)果，采用精細(xì)時程積分法計(jì)算。

2.2 縱肋處點(diǎn)傳遞矩陣

縱肋的存在會改變肋骨與圓柱殼連接處的狀態(tài)向量。肋骨對圓柱殼的作用力主要表現(xiàn)為法向力、切向力、縱向力及縱向力繞軸線的彎矩。將這些作用力充分考慮后得到的縱肋運(yùn)動微分方程[10]為

(9)

利用圓柱殼在縱肋骨處的連續(xù)條件和平衡條件，并將位移和內(nèi)力沿軸向模態(tài)展開，整理后可寫成如下矩陣方程：

(10)

3 數(shù)值算例及討論

按傳遞方向?qū)鰝鬟f矩陣和點(diǎn)傳遞矩陣依次相乘，結(jié)合邊界條件得到頻率方程，采用精細(xì)時程積分法計(jì)算，通過MATLAB軟件編程計(jì)算得到圓柱殼的固有頻率，將固有頻率代入傳遞矩陣，利用各截面狀態(tài)向量非零元素的相對比值得到振動模態(tài)。通過算例1驗(yàn)證方法的正確性，將周向傳遞矩陣法與ANSYS仿真值進(jìn)行對比分析。算例2討論周向傳遞矩陣法和軸向傳遞矩陣法在計(jì)算圓柱殼自由振動特性的使用范圍。

(1)算例1。兩端簡支加外縱肋圓柱殼，圓柱殼體總長L為1.000 m，半徑R為0.500 m，殼體厚度h為0.010 m，肋骨截面尺寸為0.020 m×0.020 m，在圓柱殼上一共添加4根縱肋，結(jié)構(gòu)材料密度ρ為7 850 kg/m3，泊松比μ為0.3，彈性模量E=2.1×1011Pa。固有頻率計(jì)算結(jié)果如表1所示，表1括號內(nèi)的數(shù)字分別為該模態(tài)的周向波數(shù)n和軸向半波數(shù)m。模態(tài)振型如圖2～圖5所示。

由表1可知：推薦方法與ANSYS仿真結(jié)果對比，最大誤差為1.68%，最小誤差為0.01%，能滿足實(shí)際工程要求。圖2和圖3比較2階固有頻率時的模態(tài)振型。圖4和圖5比較6階固有頻率時的模態(tài)振型。推薦方法計(jì)算的模態(tài)振型與ANSYS仿真結(jié)果的模態(tài)振型基本一致，驗(yàn)證方法的正確性。

圖2 2階(4，1)模態(tài)振型(傳遞矩陣法)

圖3 2階(4，1)模態(tài)振型(ANSYS)

圖4 6階(6，1)模態(tài)振型(傳遞矩陣法)

圖5 6階(6，1)模態(tài)振型(ANSYS)

表1 傳遞矩陣法與ANSYS頻率計(jì)算結(jié)果對比 Hz

(2)算例2。兩端簡支不加肋圓柱殼。圓柱殼體總長L為0.200 m，半徑R為0.100 m，殼體厚度h為0.002 m，結(jié)構(gòu)材料密度ρ為7 850 kg/m3，泊松比μ為0.3，彈性模量E=2.1×1011Pa。按周向傳遞矩陣法計(jì)算得到的前10階對稱模態(tài)固有頻率如表2所示。為作對比，表2給出ANSYS計(jì)算結(jié)果及該算例采用軸向傳遞矩陣法計(jì)算的結(jié)果。前10階反對稱模態(tài)固有頻率與前10階對稱模態(tài)固有頻率完全相同，不再單獨(dú)列表顯示。

表2 3種不同方法前10階對稱模態(tài)固有頻率對比結(jié)果 Hz

由表2可知：在m取值較小(低頻部分)時周向傳遞矩陣法計(jì)算結(jié)果與ANSYS吻合較好，而軸向傳遞矩陣法則在n取值較小(高頻部分)時與ANSYS計(jì)算結(jié)果接近。2種方法在計(jì)算圓柱殼的自由振動特性時可以互為補(bǔ)充，即低頻部分采用周向傳遞矩陣法，高頻部分采用軸向傳遞矩陣法。

4 結(jié) 語

基于傳遞矩陣法，結(jié)合精細(xì)時程積分法，建立一種快速計(jì)算縱肋圓柱殼自由振動特性的周向傳遞矩陣法。在驗(yàn)證該方法的有效性上，討論周向傳遞矩陣法和軸向傳遞矩陣法的適用范圍和條件。結(jié)果表明，該方法與ANSYS解吻合較好，計(jì)算精度較高，且軸向傳遞矩陣法適用于圓柱殼高頻部分的計(jì)算，周向傳遞矩陣法適用于圓柱殼低頻部分的計(jì)算，由于只有兩端簡支邊界圓柱殼在軸向位移函數(shù)才可以設(shè)為三角函數(shù)形式，因此周向傳遞矩陣法僅適用于簡支邊界。