中厚橢球殼自由振動動力剛度法分析

陳旭東， 葉康生

(1.江蘇科技大學 船舶與建筑工程學院，江蘇 張家港　215600；2. 清華大學 土木工程系 土木工程安全與耐久教育部重點實驗室，北京　100084)

?

中厚橢球殼自由振動動力剛度法分析

陳旭東1， 葉康生2

(1.江蘇科技大學 船舶與建筑工程學院，江蘇 張家港215600；2. 清華大學 土木工程系 土木工程安全與耐久教育部重點實驗室，北京100084)

摘要：介紹精確動力剛度法分析中厚橢球殼自由振動具體實施方法，據環向波數不同將中厚橢球殼自由振動分解為一系列確定環向波數的一維振動；利用控制方程Hamilton形式建立動力剛度關系，用常微分方程求解器COLSYS求解控制方程獲得單元動力剛度，用Wittrick-Williams算法求得該環向波數下橢球殼自振頻率。數值算例給出中厚圓球殼及橢球殼不同邊界條件的自振頻率，驗證動力剛度法高效、可靠、精確。

關鍵詞：橢球殼；自由振動；動力剛度法；Wittrick-Williams算法；Hamilton形式

橢球殼廣泛用于土木工程、化工、機械、航空航天等重要行業。準確求解橢球殼的自由振動，可為進一步動力分析提供基本參數，具有重要理論意義及工程價值。Lamb[1]最早分析橢球殼自由振動，DiMaggio等[2-5]對扁平及長橢球殼自由振動進行研究。以上分析均基于薄殼理論。對較厚殼體，應考慮剪切變形及轉動慣量影響。Sai Ram等[6]基于一階剪切變形理論采用有限元法計算復合材料球殼的自振頻率；Hosseini-Hashemi等[7]用輔助勢函數及分離變量法獲得球殼面板自由振動解析解；Su等[8]用傅里葉級數及Rayleigh-Ritz法分析中厚分段球殼的自由振動；Kang等[9-10]用三維理論結合Ritz法求解變厚度球殼的自由振動。

本文視結構為無限自由度分布參數體系，將控制方程降階后用常微分方程求解器COLSYS[11-12]獲得精確動力剛度，用Wittrick-Williams算法[13-14](W-W算法)對頻率計數，求得結構自振頻率。

動力剛度法因分析精確，多用于桿系結構自由振動問題[15-17]。在板殼領域，El-Kaabazi等[18]將動力剛度法用于軸對稱薄圓柱殼自由振動分析。本文從考慮橫向剪切變形及轉動慣量的中厚殼理論出發，推導中厚橢球殼自由振動控制方程與動力剛度關系，該方法可直接用于雙曲殼、拋物殼等旋轉殼的自由振動求解。

1基本方程

圖1　橢球殼子午線及坐標系示意圖Fig.1 Diagram of an elliptical shell meridian and coordinates

橢球殼見圖1，設u，v，w分別為橢球殼中面某點沿子午線方向、環向及法向的線位移；φ為該點法線與對稱軸z的夾角；θ為緯圓環向；E為橢球殼彈性模量；ν為泊松比；ρ為密度；h為厚度，沿環向封閉。在OXZ平面內，子午線方程為

(1)

據中厚殼理論，橢球殼中面應變與位移關系為

{ε}=[L]{Δ}

(2)

式中：{ε}={εφ,εθ,εφθ,χφ,χθ,χφθ,γφ,γθ}T為應變向量；{Δ}={u,v,w,ψφ,ψθ}T為位移向量；ψφ，ψθ分別為中面沿φ、θ方向角位移；[L]為微分算子矩陣，即

類似，橢球殼內力、中面應變關系為

{N}=[D]{ε}

(4)

式中： {N}={Nφ,Nθ,Nφθ,Mφ,Mθ,Mφθ,Qφ,Qθ}T為內力向量；[D]為剛度矩陣，即

(5)

圖2　橢球殼橫斷面及頂點開口r0Fig.2 Transverse section of an ellipsoidal shell and opening r0 at apex

由于微分算子L中出現1/r項，對完全橢球殼，在子午線頂點因r=0會造成L矩陣奇異。本文在頂點開小圓孔，即在子午線頂點開半徑為r0圓孔，見圖2，r0非常小時可認為橢球殼在頂點接近封閉。

2自由振動對偶系

在橢球殼上截取一微元，據幾何關系，此微元面積為dA=rRdθdφ。將微元應變能在整個中面S上積分，得橢球殼自由振動應變能為

(6)

橢球殼自由振動動能為

(7)

由Hamilton原理

(8)

可推出殼體自由振動運動方程。取動位移函數為

(9)

式中：n為環向波數；ω為自振圓頻率；un(φ)，vn(φ)，wn(φ)，ψφn(φ)，ψθn(φ)為子午向振型函數。

將該動位移函數代入運動方程，可得中厚殼自由振動控制方程，寫成Hamilton形式為

[J]{z′}=[S]{z}

(10)

式中：( )′為對φ求導。

(11)

式中：I為5階單位矩陣；{z}為狀態向量，由環向波數n下振型對偶量(位移、內力)組成，即

{z}={q1q2q3q4q5p1p2p3p4p5}T

(12)

式中：內力量為

(13)

位移量為

q1=un,q2=vn,q3=wn,q4=ψφn,q5=ψθn

(14)

式(10)中[S]為10階對稱矩陣，記sij為元素，則其中上三角非零元素分別為

(15)

3動力剛度

圖3　橢球殼沿子午線方向殼段劃分Fig.3 The mesh division along the meridian

對式(10)中控制方程，據位移量及內力量，類似桿件理論建立橢球殼自由振動動力剛度關系。由于方程的復雜性，式(10)用常微分方程求解器COLSYS進行數值求解。取橢球殼曲率半徑與z軸夾角φ為基本坐標，并按φ將橢球殼沿子午線劃分成ne個殼段單元，見圖3。

以單元(e)為例，設兩端φ坐標為φa、φb，單元端位移向量g0gggggge及力向量{F}e與位移、內力函數間關系為

g0gggggge={q1(φa),q2(φa),q3(φa),q4(φa),q5(φa),

q1(φb),q2(φb),q3(φb),q4(φb),q5(φb)}T

{F}e={-p1(φa),-p2(φa),

-p3(φa),-p4(φa),-p5(φa),

p1(φb),p2(φb),p3(φb),p4(φb),p5(φb)}T

(16)

據剛度定義，求解式(10)時對單元端部節點位移向量依次取10個不同單位向量，即對式(10)依次施加邊界條件為

g0gggggge={ej},(j=1,…,10)

(17)

式中：{ej}為第j元素單位值單位向量。

解得此時單元端部節點力向量，依次排列則得橢球殼殼段單元動力剛度矩陣，即

(18)

將所有殼段單元動力剛度陣按常規有限元集成，即可獲得橢球殼自由振動的整體動力剛度矩陣K。

4Wittrick-Williams算法

W-W算法認為結構所有頻率中低于給定值ω*的頻率個數為

J=J0+s{K(ω*)}

(19)

式中：J0為低于ω*的單元固端頻率總數；s{ }為負數計數；s{K(ω*)}為用高斯消去法將整體動力剛度矩陣K(ω*)消成上三角矩陣后(不換行)主對角線上負元素個數。

(20)

(21)

圖4　COLSYS對殼段單元(e)內部子區間劃分Fig.4 The sub-mesh division of COLSYS on element (e)

(22)

(23)

式中：

(24)

(25)

5數值算例

通過3個算例驗證本文方法的正確性。

5.1算例1

一端固定一端自由球殼段自由振動。令子午線方程(1)中a=b，則橢球殼退化為標準球殼。Gautham等[19]用有限元方法分析頂點有30°、60°截斷開口的中厚半球殼自由振動，邊界條件為下端固定、上端自由，見圖5。

圖5　截斷開口球殼示意圖Fig.5 Diagram of hemi-spherical caps with cutout

5.2算例2

圖6　兩端自由橢球殼示意圖Fig.6 Diagram of a free-free ellipsoidal shell

表1　30°開口半圓球殼的頻率參量Ω

表2　60°開口半圓球殼的振動頻率參量Ω

表3　b/a=1/2時兩端自由橢球殼自振頻率Ω

表4　b/a=2時兩端自由橢球殼自振頻率Ω

5.3算例3

頂點封閉半橢球殼自由振動。文獻[20]用三維彈性理論結合Ritz法對上頂點封閉、下端自由的半橢球殼進行分析，殼體參數為b/a=3，h/a=1/10，泊松比ν=0.3。為避免頂點幾何算子中1/r的奇異性，開小口，開口大小為r0/a=10-5。本文亦用開小口方法避免頂點處奇異性，開口大小同文獻[20]。用本文方法與文獻[20]的三維彈性Ritz法所求半橢球殼在不同環向波n下的最小非零頻率參量Ω(表達式同例2)見表5，可見吻合很好。

表5　半橢球殼的自振頻率Ω

6結論

闡述應用動力剛度法分析中厚橢球殼自由振動過程，據結構的軸對稱性，將橢球殼自由振動問題降為一維自由振動問題，由一維問題的內力量、位移量建立單元動力剛度關系，用COLSYS求解動力剛度關系滿足的常微分方程邊值問題，獲得單元動力剛度。并用W-W算法獲得結構頻率。

本文方法通用性強，可直接推廣至其它軸對稱旋轉殼(如雙曲殼、拋物殼等)的自由振動分析。

參 考 文 獻

[1] Lamb H. On the vibration of a spherical shell[J]. Proceedings of the London Mathematical Society, 1883, 14: 50-56.

[2] DiMaggio F L, Silbiger A. Free extensional torsional vibrations of a prolate spheroidal shell[J]. Journal of Acoustical Society of America, 1961, 33(1): 56-58.

[3] Penzes L E, Burgin G. Free vibration of thin isotropic oblate-spheroidal shells[J]. Journal of Acoustical Society of America, 1966, 39(1): 8-13.

[4] Niordson F I. Free vibrations of thin elastic spherical shells[J].International Journal of Solids and Structures,1984,20(7):667-687.

[5] Al-Jumaily A M, Najim F M. An approximation to the vibrations of oblate spheroidal shells[J]. Journal of Sound and Vibration, 1997, 204(4): 561-574.

[6] Sai Ram K S, Sreedhar Babu T. Free vibration of composite spherical shell cap with and without a cutout [J]. Computers and Structures, 2002, 80(23): 1749-1756.

[7] Hosseini-Hashemi S, Fadaee M. On the free vibration of moderately thick spherical shell panela new exact closed-form procedure[J]. Journal of Sound and Vibration,2011,330(17): 4352-4367.

[8] Su Z, Jin G, Ye T. Free vibration analysis of moderately thick functionally graded open shells with general boundary conditions [J]. Composite Structures,2014,117:169-186.

[9] Kang J H, Leissa A W. Three-dimensional vibrations of thick spherical shell segments with variable thickness[J]. International Journal of Solids and Structures,2000,37: 4811-4823.

[10] Kang J H, Leissa A W. Three-dimensional vibration analysis of solids and hollow hemispheres having varying thicknesses with and without axial conical holes[J]. Journal of Vibration and Control, 2004, 10: 199-214.

[11] Ascher U, Christiansen J, Russell R D. Collocation software for boundary value ODEs[J]. ACM Transactions on Mathematical Software, 1981, 7(2): 209-222.

[12] Ascher U, Christiansen J, Russell R D. Algorithm 569, COLSYS: collocation software for boundary value ODEs [D2] [J]. ACM Transactions on Mathematical Software,1981,7(2):223-229.

[13] Williams F W, Wittrick W H. An automatic computational procedure for calculating natural frequencies of skeletal structures[J]. International Journal of Mechanical Sciences, 1970, 12(9): 781-791.

[14] Wittrick W H, Williams F W. A general algorithm for computing natural frequencies of elastic structures[J]. Quarterly Journal of Mechanics and Applied Mathematics, 1971, 24(3): 263-284.

[15] Yuan S, Ye K, Xiao C, et al. Exact dynamic stiffness method for non-uniform Timoshenko beam vibrations and Berboulli-Euler column buckling[J]. Journal of Sound and Vibration, 2007, 303(3/4/5): 526-537.

[16] 葉康生,趙雪健. 動力剛度法求解平面曲梁面外自由振動問題[J]. 工程力學,2012,29：1-8.

YE Kang-sheng, ZHAO Xue-jian. Dynamic stiffness method for out-of-plane free vibration analysis of planar curved beams [J]. Engineering Mechanics, 2012, 29: 1-8.

[17] Su H, Banerjee J R, Cheung C W. Dynamic stiffness formulation and free vibration analysis of functionally graded beams [J]. Composite Structures, 2013, 106: 854-862.

[18] El-Kaabazi N, Kennedy D. Calculation of natural frequencies and vibration modes of variable thickness cylindrical shells using the Wittrick-Williams algorithm[J]. Computers and Structures, 2012, 104/105: 4-12.

[19] Gautham B P, Ganesan N. Free vibration analysis of thick spherical shells[J]. Computers and Structures,1992,45(2): 307-313.

[20] Shim H J, Kang J H. Free vibration of solid and hollow hemi-ellipsoids of revolution from a three-dimensional theory [J]. International Journal of Engineering Science,2004,42:1793-1815.

Free vibration analysis of moderately thick elliptical shells using the dynamic stiffness method

CHENXu-dong1，YEKang-sheng2

(1. School of Naval Architecture and Civil Engineering, Jiangsu University of Science and Technology, Zhangjiagang 215600,China；2. Key Laboratory of Civil Engineering Safety and Durability of China Education Ministry, Department of Civil Engineering, Tsinghua University, Beijing 100084, China)

Abstract:The application of exact dynamic stiffness method to the free vibration analysis of moderately thick elliptical shells was introduced. The free vibration of moderately thick elliptical shells was decomposed into a series of one-dimensional vibration problems corresponding to structural vibration modes with different circumferential wave numbers. For each one-dimensional vibration problem, the governing equation was written in Hamilton form, from which the dynamic stiffness expression of the one-dimensional problem was derived. The governing equations were solved by using the ordinary differential equations solver COLSYS and the dynamic stiffnesses of elements were obtained. By applying the Wittrick-Williams algorithm, the natural frequencies under the vibration mode with a specific circumferential wave number were found. Numerical examples of moderately thick spherical and elliptical shells with different boundary conditions were given, showing that the dynamic stiffness method is robust, reliable and accurate.

Keywords:elliptical shells; free vibration; dynamic stiffness method; Wittrick-Williams algorithm; Hamilton form

中圖分類號:TU311.3

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2016.06.015

通信作者葉康生 男，博士，副教授，1972年生

收稿日期：2015-07-09修改稿收到日期：2015-10-14

基金項目：國家自然科學基金項目(51078198)；清華大學自主科研計劃(2011THZ03);江蘇省雙創博士

第一作者 陳旭東 男，博士，講師，1984年生