基于傳遞矩陣法的縱肋圓柱殼振動特性計算

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

(1.武漢理工大學 機電工程學院，湖北 武漢 430070；2.湖北工業大學 土木建筑與環境學院，湖北 武漢 430068)

0 引 言

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

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

基于傳遞矩陣法計算縱肋圓柱殼的振動特性，按照傳遞方向將傳遞矩陣法分為周向傳遞矩陣法和軸向傳遞矩陣法。在考慮環肋與殼體相互作用時將面內和面外作用力全部加上，推導縱肋圓柱殼的場傳遞矩陣和點傳遞矩陣，采用精細時程積分法求解，保證計算過程的精度，并討論兩種傳遞矩陣法適用條件及范圍。

1 基本理論公式

基于Flugge殼體理論，殼體的應力和應變的等式按照位移的形式，內力和內力矩可表達為

(1)

式中：Nx和Nθ分別為面內的軸向力和周向力；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為橫向剪力。

徑向位移與轉角ψ之間的關系為

(2)

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

(3)

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

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

(4)

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

圖1為縱肋圓柱殼結構示例。圖1中，殼體長為L，半徑為R，且該圓柱殼滿足殼體理論中的基本假設。

圖1 縱肋圓柱殼結構示例

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

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

(5)

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

為了方便計算，同樣引入如下無量綱參數：

(6)

式中：l和H分別為長度和厚度對應的無量綱參數；λ為頻率因子。

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

(7)

式中：SYM表示對稱矩陣。

(8)

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

2.2 縱肋處點傳遞矩陣

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

(9)

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

(10)

3 數值算例及討論

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

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

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

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

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

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

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

表1 傳遞矩陣法與ANSYS頻率計算結果對比 Hz

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

表2 3種不同方法前10階對稱模態固有頻率對比結果 Hz

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

4 結 語

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