中厚橢球殼自由振動(dòng)動(dòng)力剛度法分析

陳旭東， 葉康生

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

?

中厚橢球殼自由振動(dòng)動(dòng)力剛度法分析

陳旭東1， 葉康生2

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

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

關(guān)鍵詞：橢球殼；自由振動(dòng)；動(dòng)力剛度法；Wittrick-Williams算法；Hamilton形式

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

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

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

1基本方程

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

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

(1)

據(jù)中厚殼理論，橢球殼中面應(yīng)變與位移關(guān)系為

{ε}=[L]{Δ}

(2)

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

類似，橢球殼內(nèi)力、中面應(yīng)變關(guān)系為

{N}=[D]{ε}

(4)

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

(5)

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

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

2自由振動(dòng)對偶系

在橢球殼上截取一微元，據(jù)幾何關(guān)系，此微元面積為dA=rRdθdφ。將微元應(yīng)變能在整個(gè)中面S上積分，得橢球殼自由振動(dòng)應(yīng)變能為

(6)

橢球殼自由振動(dòng)動(dòng)能為

(7)

由Hamilton原理

(8)

可推出殼體自由振動(dòng)運(yùn)動(dòng)方程。取動(dòng)位移函數(shù)為

(9)

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

將該動(dòng)位移函數(shù)代入運(yùn)動(dòng)方程，可得中厚殼自由振動(dòng)控制方程，寫成Hamilton形式為

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

(10)

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

(11)

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

{z}={q1q2q3q4q5p1p2p3p4p5}T

(12)

式中：內(nèi)力量為

(13)

位移量為

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

(14)

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

(15)

3動(dòng)力剛度

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

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

以單元(e)為例，設(shè)兩端φ坐標(biāo)為φa、φb，單元端位移向量j5i0abt0be及力向量{F}e與位移、內(nèi)力函數(shù)間關(guān)系為

j5i0abt0be={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)

據(jù)剛度定義，求解式(10)時(shí)對單元端部節(jié)點(diǎn)位移向量依次取10個(gè)不同單位向量，即對式(10)依次施加邊界條件為

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

(17)

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

解得此時(shí)單元端部節(jié)點(diǎn)力向量，依次排列則得橢球殼殼段單元?jiǎng)恿偠染仃?，?/p>

(18)

將所有殼段單元?jiǎng)恿偠汝嚢闯Ｒ?guī)有限元集成，即可獲得橢球殼自由振動(dòng)的整體動(dòng)力剛度矩陣K。

4Wittrick-Williams算法

W-W算法認(rèn)為結(jié)構(gòu)所有頻率中低于給定值ω*的頻率個(gè)數(shù)為

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

(19)

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

(20)

(21)

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

(22)

(23)

式中：

(24)

(25)

5數(shù)值算例

通過3個(gè)算例驗(yàn)證本文方法的正確性。

5.1算例1

一端固定一端自由球殼段自由振動(dòng)。令子午線方程(1)中a=b，則橢球殼退化為標(biāo)準(zhǔn)球殼。Gautham等[19]用有限元方法分析頂點(diǎn)有30°、60°截?cái)嚅_口的中厚半球殼自由振動(dòng)，邊界條件為下端固定、上端自由，見圖5。

圖5　截?cái)嚅_口球殼示意圖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°開口半圓球殼的振動(dòng)頻率參量Ω

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

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

5.3算例3

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

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

6結(jié)論

闡述應(yīng)用動(dòng)力剛度法分析中厚橢球殼自由振動(dòng)過程，據(jù)結(jié)構(gòu)的軸對稱性，將橢球殼自由振動(dòng)問題降為一維自由振動(dòng)問題，由一維問題的內(nèi)力量、位移量建立單元?jiǎng)恿偠汝P(guān)系，用COLSYS求解動(dòng)力剛度關(guān)系滿足的常微分方程邊值問題，獲得單元?jiǎng)恿偠取２⒂肳-W算法獲得結(jié)構(gòu)頻率。

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

參 考 文 獻(xiàn)

[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] 葉康生,趙雪健. 動(dòng)力剛度法求解平面曲梁面外自由振動(dòng)問題[J]. 工程力學(xué),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

文獻(xiàn)標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2016.06.015

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

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

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

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