王文勝,程耿東,郝 鵬

(1.河南科技大學 土木工程學院 工程力學系，洛陽 471023；2.大連理工大學 工程力學系 工業(yè)裝備結構分析國家重點實驗室，大連 116023)

?

基于超梁降階模型的蒙皮加筋圓柱殼頻率分析

王文勝1,程耿東2,郝 鵬2

(1.河南科技大學 土木工程學院 工程力學系，洛陽 471023；
2.大連理工大學 工程力學系 工業(yè)裝備結構分析國家重點實驗室，大連 116023)

在基于梁平截面假設的復雜結構動力模型降階方法基礎上，文章提出了建立梁超單元和超梁降階模型的方法。通過對比考慮剪切變形的Timoshenko梁單元剛度陣，給出了梁超單元剛度矩陣的修正方法，提高了超梁降階模型的計算精度。以蒙皮加筋圓柱殼為例，對比了超梁降階模型與幾種等效梁模型的頻率分析結果。結果表明，超梁降階模型建立方法靈活簡單，具有較高的計算精度，能夠為火箭結構動特性分析提供參考。

梁超單元；超梁降階模型；剪切變形；蒙皮加筋圓柱殼；頻率分析

0 引言

蒙皮加筋圓柱殼具有較高的剛度、較輕的重量等特性，被廣泛應用于運載火箭等結構。由于運載火箭結構的長細比較大，因此在工程研制中，經(jīng)常采用梁模型來模擬復雜的火箭結構[1-4]。針對彎曲、扭轉(zhuǎn)和縱向的中低階振動模態(tài)的火箭結構頻率動特性分析，在火箭結構總體設計中具有十分重要的地位，它是載荷計算、姿態(tài)穩(wěn)定性分析、POGO振動分析等工作的基礎，而選取合適的計算模型則是開展相關分析的前提。

近年來，隨著有限元分析軟件和計算機求解能力的快速發(fā)展，已經(jīng)可進行全箭結構的三維動力學建模研究工作。文獻[5]采用了三維動力學模型，分析了日本H-II運載火箭的俯仰、偏航和扭轉(zhuǎn)振動之間的耦合效應。采用這種模型分析的優(yōu)點是便于研究速率陀螺等敏感元件安裝位置的局部振型和振型斜率等問題；缺點是計算耗時，分析結果中大量的局部模態(tài)掩蓋了少量設計關心的整體結構縱、橫和扭轉(zhuǎn)模態(tài)。而降階梁模型具有較少的自由度，能夠快速地獲取結構的整體動力學響應特性，為火箭結構設計提供重要設計參數(shù)。

結構力學[6]中雖然給出了這類蒙皮加筋結構桁條面積等效和剛度等效的方法，但對桁條作用的物理原理考慮不夠，沒有明確給出這種等效方法的適用范圍。由于實際工程中對火箭的低階頻率要求較高，直接用教科書的方法難以滿足航天工程的要求。邢譽峰等[7-9]提出了基于桁條等效厚度方法的梁模型，分析了蒙皮加筋圓柱殼彎曲剛度等效方法的精度，研究了等效方法中的扭轉(zhuǎn)剛度計算問題，并指出等效厚度梁模型可直接用于計算彎曲和縱向固有振動頻率。該方法對于火箭整體結構的相關分析具有實際應用和參考價值，但扭轉(zhuǎn)頻率的計算結果偏高，需要進行修正，本文稱這種梁模型為等效厚度梁模型。欒宇[10]給出了根據(jù)結構靜力響應計算這類結構等效剛度及截面參數(shù)的梁模型降階方法，本文稱為靜力等效梁模型。本文在基于梁平截面假設的動力模型降階方法[11-12]的基礎上，提出了建立梁超單元的方法，且通過與Timoshenko梁單元剛度矩陣對比，對梁超單元剛度陣引入了考慮剪切變形的修正，得到了計算精度較高的超梁降階模型。通過典型蒙皮加筋圓柱殼算例，與等效厚度梁模型及靜力等效梁模型進行了對比，數(shù)值結果表明了超梁降階模型的有效性。

1 超梁降階模型的建立

針對細長梁式結構，文獻[11-12]中提出了基于梁平截面假設的動力模型降階方法，此處稍作敘述。經(jīng)典的Euler梁彎曲理論以及考慮剪切變形的Timoshenko梁彎曲理論中，都假定變形前垂直于中心線的截面變形后仍保持為平面[13]。如圖1所示，根據(jù)細長梁式結構的特點，在發(fā)生非局部變形時，截面的變形可近似為剛體位移，可采用平截面假定來描寫。根據(jù)這一假定，任意一個截面i上任一點j和截面形心處的6個基本準剛體模態(tài)之間可建立如下位移映射關系：

uj=Rjqi

(1)

其展開形式為

(2)

式中uj表示第j點的位移向量；qi表示截面i形心處的準剛體模態(tài)；Rj表示j點的位移轉(zhuǎn)換矩陣；xcj,ycj,zcj表示第j點坐標與截面i形心坐標的差值。

圖1 模型降階示意圖

對于結構上第i個截面，截面內(nèi)的m個節(jié)點均凝聚到該截面形心處，則有：

Ui=Tiqi

(3)

式中Ui=(ui1,…,uim)T表示這一截面內(nèi)m個節(jié)點的位移向量；Ti=(Ri1,…,Rim)T表示這一截面的位移轉(zhuǎn)換矩陣。

在此基礎上，本文取任意2個相鄰截面(截面i和截面i+1)之間的結構部分為一梁段，梁段上所有節(jié)點的位移U可表示為

(4)

假設梁段的剛度陣和質(zhì)量陣分別為k、m，利用式(4)可得到：

kR=TTkT,mR=TTmT

(5)

式(4)中，Q是12×1的列向量，因此式(5)中凝聚得到kR、mR均為12×12的矩陣，分別表示凝聚后的梁超單元的剛度陣和質(zhì)量陣。式(4)中的T表示建立梁超單元的降階基向量。kR、mR可按照序列組裝成一個超梁降階模型的總剛度陣和總質(zhì)量陣，如果模型均勻，只需要建立一次梁超單元；反之，則需要建立若干個梁超單元。

2 超梁降階模型的修正

Timoshenko梁單元引入了橫向剪切變形來修改梁單元剛度陣中的彎曲剛度項，具體推導過程可在文獻[14-16]中查到?？紤]一根主軸與坐標x重合，梁截面關于x-y平面對稱的空間均勻梁，其Timoshenko梁單元剛度矩陣可表示為

(6)

(7)

(8)

如果Iy=Iz，則φy=φz，可統(tǒng)一表示成φ。

由式(7)可得

(9)

又因為

(10)

將式(9)及式(10)代入到式(8)中，φ值可求得：

(11)

kT=kb1+kb2+kb3+kat

(12)

式中kb1只包含(Y1,Y2,Z1,Z2)項；kb2只包含(Y3,Z3)項；kb3只包含(Y4,Z4)項；kat只包含與軸向和扭轉(zhuǎn)剛度有關的(X,S)項。

上述通過降階得到的梁超單元剛度陣中，本文認為其彎曲剛度項與Timoshenko梁單元具有相同的物理意義，同樣可分解為與式(12)類似的4部分：

kR=kRb1+kRb2+kRb3+kRat

(13)

(14)

將整個結構沿長度方向剖分為p個這樣的梁段，則修正過彎曲剛度項后的超梁降階模型的總剛度陣可表示為

(15)

以ANSYS軟件為例，在建立實際殼體結構的復雜有限元模型后，本文采用APDL(ANSYSParametricDesignLanguage)二次開發(fā)語言實現(xiàn)梁段剛度矩陣k、質(zhì)量矩陣m的數(shù)據(jù)輸出。在ANSYS的輸出數(shù)據(jù)中沒有模型降階減縮基T矩陣，本文同時利用APDL語言程序?qū)崿F(xiàn)了T矩陣的建立和輸出。利用式(5)得到梁超單元的剛度陣kR和質(zhì)量陣mR，并將kR分解為4個部分，提取剛度陣中的(X,S)值計算剪切修正系數(shù)進行修正；然后，組裝超梁降階模型的總剛度陣和質(zhì)量陣，調(diào)用特征值分析程序?qū)崿F(xiàn)原結構的快速頻率分析。

3 等效厚度梁模型和靜力等效梁模型概述

為了對比方便，本節(jié)簡述一下等效厚度梁模型及靜力等效梁模型。文獻[7-9]中給出了基于桁條等效厚度方法建立的等效梁模型，在蒙皮和桁條材料彈性常數(shù)接近時，等效光筒圓柱結構的平均厚度為

(16)

在此基礎上，對于蒙皮加筋結構的等效梁模型，其拉伸、彎曲和扭轉(zhuǎn)剛度由材料彈性常數(shù)和橫截面面積、慣性矩、極慣性矩確定。其中橫截面面積、慣性矩、極慣性矩按式(17)計算。

(17)

(18)

利用式(16)～式(18)計算得到這些結構參數(shù)后，可很方便地在有限元軟件中建立相應的等效梁模型，用來模擬原蒙皮加筋結構。本文利用ANSYS軟件實現(xiàn)這一過程。

4 數(shù)值算例

本章通過2個具體算例的頻率分析，以比較本文提出的超梁降階模型與其它模型。

4.1 均勻蒙皮加筋圓柱殼算例

算例1為圖2所示自由-自由的均勻ISOGIRD網(wǎng)格蒙皮加筋圓柱殼。模型長9 976 mm，外徑2 000 mm，蒙皮厚度為10 mm，網(wǎng)格筋條高30 mm，筋條厚25 mm。采用鋁合金，彈性模量70 GPa，泊松比0.3，密度2.7×10-3g/mm3。采用ANSYS殼單元(SHELL181)建立有限元模型，3 472個節(jié)點，20 832個自由度。

圖2 均勻網(wǎng)格蒙皮加筋圓柱殼有限元模型及局部放大圖

利用式(16)可求出光筒殼等效厚度為16.2 mm，并根據(jù)式(17)計算其橫截面面積、慣性矩、極慣性矩。

δN=0.014 877 mm,θp=0.158×10-3rad

φT=0.38×10-7rad

進而可按照式(18)計算靜力等效梁模型參數(shù)如下：

分別按照上面得到的數(shù)據(jù)，利用BEAM188梁單元，建立與超梁降階模型自由度數(shù)相等的等效梁模型，用于模擬原結構的動力學分析。

基于超梁降階模型方法，將原結構模型沿軸向劃分為55個梁段。由于模型是均勻結構，因此只需要建立一個梁超單元，便可按照次序組裝成具有55個梁超單元，336個自由度的超梁降階模型。對于這類結構，k取值為2.0[17]，提取梁超單元剛度陣中的(X、S)值，并計算得到φ=0.340 8，則按照式(14)計算并取近似值β=0.73，進而可按照式(15)修正超梁降階模型。

在此基礎上，分別利用上述3種降階梁模型模擬了原結構模型在自由-自由、一端固支及兩端固支邊界條件下的特征值分析，計算結果分別如表1～表3所示。

表1 自由-自由狀態(tài)下幾種降階梁模型頻率計算結果對比

表2 一端固支狀態(tài)下幾種降階梁模型頻率計算結果對比

表3 兩端固支狀態(tài)下幾種降階梁模型頻率計算結果對比

從表1～表3中的頻率計算結果可知，靜力等效梁模型計算精度最高，但建立一次降階梁模型至少需要3次原結構靜力分析，如果原結構參數(shù)改變(如長度變化)，則需要重新求解式(18)中的參數(shù)。等效厚度梁模型計算精度較差，尤其是扭轉(zhuǎn)頻率的計算結果，說明等效厚度梁模型不能直接用于計算扭轉(zhuǎn)頻率，文獻[8]中給出了具體的修正方法，此處不再敘述。

超梁降階模型的各階頻率計算結果與相應階次ANSYS軟件的計算結果很接近。當對實際工程結構低階動力特性分析精度要求不很高時，為設計人員提供了一個有效簡單的計算模型，且建立模型工作量很小，只要利用已經(jīng)建立的原結構有限元模型的數(shù)據(jù)建立一個梁超單元，便可組裝成不同長度的超梁降階模型，剪切修正系數(shù)的計算過程也較為簡單。需要指出的是超梁降階模型一端固支一端自由及兩端固支頻率結果，是在自由-自由狀態(tài)的降階模型上施加相應的邊界約束條件計算得到。

為了進一步修改和使用降階梁模型，工程師們更希望獲得降階梁模型的截面參數(shù)，即得到降階梁模型的軸向剛度、抗彎剛度、抗扭剛度和質(zhì)量分布等。表4分別給出了利用3種梁模型計算得到的圖1所示的ISOGIRD網(wǎng)格蒙皮加筋圓柱殼的軸向剛度EA及截面抗彎剛度EI。

表4 3種梁模型的軸向及抗彎剛度

利用原結構截面幾何形狀計算截面參數(shù)時，不同截面的軸向剛度和抗彎剛度均不同。因此，這里沒有把它們簡單列入表中。對于本算例，按截取位置不同(軸向加筋處和環(huán)向加筋處)，分別計算得到相應的軸向剛度EA及截面抗彎剛度EI為：軸向加筋處，軸向剛度7.104 2×109N，截面抗彎剛度34.630×1014N·mm2；環(huán)向加筋處相應值為17.052×109N及82.828×1014N·mm2。從這些結果中可看到，等效厚度梁模型截面參數(shù)與原結構軸向加筋處的截面參數(shù)基本一致，這是因為建立等效厚度梁模型時，只考慮了沿軸向筋條的等效。超梁降階模型的截面參數(shù)介于原結構軸向加筋與環(huán)向加筋兩處的截面之間，這是因為超梁降階模型考慮了環(huán)向加肋的影響，更真實地反映了原結構的截面性質(zhì)。此外，由于本文在超梁降階模型上考慮了剪切變形修正，能更好地用于近似實際結構。

4.2 非均勻蒙皮加筋圓柱殼算例

算例2為圖3所示沿長度方向加筋非均勻的自由-自由圓柱殼結構，模型由蒙皮、軸向及環(huán)向桁條組成。模型長20 m，蒙皮厚度0.01 m，桁條沿徑向高度均為0.1 m。將結構分為3部分：第一部分長度10 m，直徑1 m，軸向桁條厚度為0.05 m，環(huán)向桁條厚度為0.03 m；第二部分長度2 m，軸向桁條厚度為0.04 m，環(huán)向桁條厚度為0.05 m；第三部分長8 m，直徑0.8 m，軸向桁條厚度為0.03 m，環(huán)向桁條厚度為0.04 m。所用材料與算例1一致。模型共有2 064個節(jié)點，2 520個單元，12 384個自由度。

圖3 非均勻網(wǎng)格加筋圓柱殼模型及局部放大圖

根據(jù)結構形式的不同，在結構第一、三部分分別建立1種梁超單元，第二部分建立6種梁超單元；然后，將超單元進行組裝，將原結構模型降階為具有42個梁超單元，258個自由度的超梁降階模型。分別提取各種梁超單元剛度陣中的(X、S)值，并計算其修正系數(shù)，本文取這些系數(shù)的平均值β=0.95作為該超梁降階模型的修正系數(shù)。同時，利用式(16)～式(18)，分別計算等效厚度及靜力等效梁模型的參數(shù)，利用ANSYS軟件建立與超梁降階模型自由度數(shù)相等的等效梁模型。表5給出了自由-自由邊界條件下算例2采用3種梁模型的頻率計算結果。

表5所示計算結果表明，對于非均勻蒙皮加筋結構，超梁降階模型也能用于原結構的快速近似動力分析，且計算精度也高。在此基礎上，在自由-自由狀態(tài)的降階模型上施加相應的邊界約束條件，即可計算相應約束條件下原結構的頻率，只需非常小的額外工作量。

表5 非均勻網(wǎng)格加筋圓柱殼降階梁模型頻率計算結果對比

5 結論

(1)在基于梁平截面假設的結構動力模型降階方法的基礎上，給出了建立梁超單元和超梁降階模型的方法。

(2)通過對比Timoshenko梁單元剛度矩陣，考慮了剪切變形，給出梁超單元剛度矩陣的修正方法，建立了計算精度較高的超梁降階模型。

(3)典型蒙皮加筋圓柱殼模型算例的頻率分析結果表明，超梁降階模型具有較高計算精度，且可處理不同情況和約束條件下的結構模型，對火箭結構頻率動特性分析具有實際應用和參考價值。

[1] Zipfel P H. Modeling and simulation of aerospace vehicle dynamics[M]. 2nded.. American Institute of Aeronautics and Astronautics, Inc., Virginia., USA, 2007.

[2] 王龍生,張德文. 火箭結構有限元分析的若干問題[J]. 強度與環(huán)境, 1988(3): 45-53.

[3] 王毅,朱禮文,王明宇,等. 大型運載火箭動力學關鍵技術及其進展綜述[J]. 導彈與航天運載技術, 2000, 243(1): 29-37.

[4] 邱吉寶,王建民. 運載火箭模態(tài)試驗仿真技術研究新進展[J]. 宇航學報, 2007, 28(3): 515-521.

[5] Yoshiki Morino. Vibration test of 1/5 H-II launch vehicle[R]. AIAA 87-D783.

[6] 龔堯南. 結構力學 [M]. 北京: 北京航空航天大學出版社, 2001.

[7] 潘忠文,王旭,邢譽峰,等. 基于梁模型的火箭縱橫扭一體化建模技術[J]. 宇航學報, 2010, 31(5): 1310-1316.

[8] 潘忠文, 邢譽峰, 楊陽. 蒙皮加筋圓柱殼扭轉(zhuǎn)頻率的三種計算模型[J]. 北京航空航天大學學報, 2011, 37(9): 1156-1159.

[9] 邢譽峰,潘忠文,楊陽. 蒙皮加筋圓柱殼彎曲頻率的三種計算模型[J]. 北京航空航天大學學報, 2012, 38(4): 1-6.

[10] 欒宇. 航天器結構中螺栓法蘭連接的動力學建模方法研究 [D]. 大連: 大連理工大學, 2011.

[11] 王文勝,程耿東,李取浩. 基于梁平截面假設的復雜細長結構動力模型降階方法[J]. 計算力學學報, 2012, 29(3): 295-332.

[12] Wang Wen-sheng, Cheng Geng-dong, Li Qu-hao. Fast dynamic performance optimization of complicated beam type structure based on two new reduced physical models[J]. Engineering Optimization, 2013, 45(7):835-850.

[13] Timoshenko S P, Gere J M. Mechanics of materials[M]. 1sted.. Van Nostrand Reinhold Company, New York, 1972.

[14] Cowper G R. The shear coefficient in Timoshenko's beam theory[J]. Journal of Applied Mechanics, 1966, 33(2): 335-340.

[15] Narayanaswami R, Adelman H M. Inclusion of transverse shear deformation in finite element displacement formulations[J]. AIAA Journal, 1974, 12:1613-1614.

[16] Severn R T. Inclusion of shear deflection in the stiffness matrix for a beam[J]. The Journal of Strain Analysis for Engineering Design, 1970, 5(4):239-241.

[17] Robert D C, David S M, Michael E P. Concepts and applications of Finite Element Analysis[M]. 4thed.. John Wiley & Sons Inc, New York, 2001.

(編輯：薛永利)

Frequency analysis of stiffened cylindrical shell based on reduced super beam model

WANG Wen-sheng1, CHENG Geng-dong2, HAO Peng2

(1.Department of Engineering Mechanics, School of Civil Engineering,Henan University of Science and Technology, Luoyang 471023, China；2.Department of Engineering Mechanics, State Key Laboratory of Structural Analysis for Industrial Equipment,Dalian University of Technology, Dalian 116023, China)

A type of super beam element and reduced super beam model were presented by means of the model reduction method, which is based on the plane cross-section assumption of beam. The stiffness matrix of the super beam element was modified by considering the effect of shear deformation in an elegant way, resulting in a reduced super beam model with high accuracy. Typical stiffened cylindrical shell model, served as the illustrative example for frequency analysis, were employed to compare the reduced super beam model with several other equivalent beam models. The results show that the reduced super beam model is easy to construct with high accuracy, which can provide references for the dynamic analysis of rocket structures.

super beam element；reduced super beam model；shear deformation；stiffened cylindrical shell；frequency analysis

2014-07-03；

2014-08-11。

國家自然科學基金(11402077,11402049)；河南科技大學科學研究基金(13350046)；中國博士后科學基金(2014M551070)。

王文勝(1983—)，講師，研究方向為復雜結構動力模型降階及結構優(yōu)化設計。E-mail:wswang@live.cn

V414

A

1006-2793(2015)03-0401-06

10.7673/j.issn.1006-2793.2015.03.020