李玉韋,王 博,程耿東

(大連理工大學工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室工程力學系,大連 116024)

基于新型降階模型的網(wǎng)格加筋圓柱殼頻率分析

李玉韋,王 博,程耿東

(大連理工大學工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室工程力學系,大連 116024)

針對網(wǎng)格加筋圓柱殼頻率分析計算量大的問題，提出一種基于多項式及梁單元形函數(shù)的模型降階方法，即運用多項式及梁單元形函數(shù)將復雜結(jié)構(gòu)有限元的節(jié)點位移轉(zhuǎn)化為主節(jié)點位移，實現(xiàn)模型降階。并通過控制多項式階數(shù)及梁單元個數(shù)來調(diào)整降階模型的分析精度，而降階模型的計算效率較精細模型有顯著提高。以網(wǎng)格加筋圓柱殼為例，對比本文提出的模型降階方法與3-D實體模型及其他模型降階方法的頻率分析結(jié)果，結(jié)果表明，提出的模型降階方法不僅能捕捉到結(jié)構(gòu)的整體模態(tài)，還能反映出結(jié)構(gòu)部分局部模態(tài)，適用范圍廣，能夠為工程結(jié)構(gòu)設計提供簡單有效的計算模型。

多項式函數(shù)；梁單元形函數(shù)；模型降階；網(wǎng)格加筋圓柱殼；頻率分析

0 引言

網(wǎng)格加筋圓柱殼由于較高的比剛度和比強度廣泛應用于航空航天結(jié)構(gòu)，如運載火箭的燃料貯箱和級間段等[1-2]。運載火箭結(jié)構(gòu)的彎曲、扭轉(zhuǎn)、縱向等整體模態(tài)及部分局部模態(tài)是載荷計算、姿態(tài)穩(wěn)定性分析、速率陀螺儀安放位置確定等工作的基礎(chǔ)。詳細描述結(jié)構(gòu)細節(jié)的有限元模型往往具有幾十萬甚至上千萬個自由度，幾乎不可能采用這樣的模型進行結(jié)構(gòu)動力特性分析，而且基于此類有限元模型得到的頻率分析結(jié)果中含有大量的局部模態(tài)，反而干擾了對結(jié)構(gòu)總體性能的認識。因此，在結(jié)構(gòu)設計初期，選取合適的計算模型是開展相關(guān)分析的前提。

為提高網(wǎng)格加筋圓柱殼頻率分析的計算效率，邢譽峰等[3]、潘忠文等[4]研究了加筋圓柱殼彎曲剛度的面積等效方法和慣性矩等效方法，分別分析了面積等效時截面慣性矩的精度和慣性矩等效時面積的等效精度，提出了用于加筋圓柱殼頻率分析的等效厚度模型，并指出等效厚度模型可直接用于計算結(jié)構(gòu)的彎曲和縱向頻率，但扭轉(zhuǎn)頻率計算精度較差，需進行修正。Cheng等[5]、Cai等[6-7]通過改進漸近均勻化方法(NIAH)預測三維周期性結(jié)構(gòu)剛度特性，將復雜結(jié)構(gòu)等效為具有類似剛度特性的簡單結(jié)構(gòu)。這種方法不需要復雜的數(shù)學推導，可以利用商業(yè)軟件中提供的各種單元和建模技術(shù)，使單胞模型保持較小的規(guī)模。鄭淑飛等[8-9]提出一種基于變形修正的動力模型簡化方法，該方法以加筋圓柱殼精細有限元模型為基礎(chǔ)，根據(jù)運動同步性假設將結(jié)構(gòu)分成若干同步性區(qū)域，建立結(jié)構(gòu)整體位移與局部位移模式之間的轉(zhuǎn)換關(guān)系，實現(xiàn)模型降階，但位移轉(zhuǎn)換陣中未考慮原結(jié)構(gòu)的轉(zhuǎn)角自由度，對含有梁、殼等元組件結(jié)構(gòu)的分析精度較差。在此基礎(chǔ)上，王文勝等[10，12]、Cheng等[11]將復雜結(jié)構(gòu)模型劃分為不同的梁段，通過梁平截面假設及位移插值函數(shù)將每一梁段降階為一個梁超單元，按照一定的次序組裝成超梁模型，并引入橫向剪切系數(shù)修正梁超單元的剛度陣。針對具有大開口的結(jié)構(gòu)，本文通過數(shù)值方法引入了表示截面翹曲變形的翹曲基向量，進一步提高了扭轉(zhuǎn)頻率的計算精度，本文稱該方法為平截面梁模型方法。文獻[12]通過網(wǎng)格加筋結(jié)構(gòu)、含開口光筒結(jié)構(gòu)驗證了平截面梁模型方法的有效性，但該方法僅適用于計算結(jié)構(gòu)的整體頻率，無法獲得結(jié)構(gòu)的局部模態(tài)。

本文在基于多項式函數(shù)的動力模型降階方法的基礎(chǔ)上，提出一種基于多項式及梁單元形函數(shù)混合的模型降階方法，該方法將結(jié)構(gòu)沿軸向劃分若干個梁單元，利用多項式及梁單元形函數(shù)建立所有節(jié)點位移與梁單元節(jié)點位移之間的轉(zhuǎn)換關(guān)系，實現(xiàn)結(jié)構(gòu)降階。建立的降階模型不僅可用于計算結(jié)構(gòu)的整體頻率，還可通過增加多項式項數(shù)及梁單元個數(shù)獲得結(jié)構(gòu)的部分局部模態(tài)。文中首先以網(wǎng)格加筋圓柱殼為例，對比本文提出的模型降階方法與3-D實體模型及其他模型降階方法的頻率分析結(jié)果。然后，針對含開口網(wǎng)格加筋圓柱殼結(jié)構(gòu)，文中通過數(shù)值方法引入翹曲變形改善降階模型扭轉(zhuǎn)頻率的計算精度。數(shù)值結(jié)果表明，本文提出的模型降階方法計算精度較高。

1 模型降階方法

模型降階是利用原結(jié)構(gòu)有限元模型的數(shù)學或力學特征，選擇一組合適的減縮基向量，將復雜結(jié)構(gòu)動力學問題降階。不考慮阻尼影響，原結(jié)構(gòu)動力分析的有限元控制方程可表示為：

(1)

引入位移轉(zhuǎn)換矩陣T(r×k)(r?k)，復雜結(jié)構(gòu)的位移可表示為：

U=TQ

(2)

其中，Q(k×1)是降階模型的位移。將式(2)帶入式(1)，并左乘TT，得到降階模型的動力學方程：

(3)

其中，MR=TTMT,KR=TTKT是降階模型的質(zhì)量陣和剛度陣，F(xiàn)R=TTF是降階模型的載荷向量。通過求解特征方程(4)獲得降階模型的固有頻率。

(KR-λkkMR)φkk=0

(4)

其中，λkk為結(jié)構(gòu)第kk階頻率，φkk為與其對應的特征向量。

不同模型降階方法的主要區(qū)別在于減縮基向量T的構(gòu)造，不同的減縮基構(gòu)造方法可能導致不同的計算精度及計算效率，文獻[13-16]給出了不同的模型降階方法，此處不再贅述。

2 基于多項式及梁單元形函數(shù)的模型降階方法

考慮如圖1所示的復雜結(jié)構(gòu)，OXYZ為全局坐標系，將結(jié)構(gòu)沿軸向劃分為若干個梁單元，取每個梁單元的節(jié)點為主節(jié)點。首先利用多項式建立任意一點與其梁軸線上的投影點之間的位移轉(zhuǎn)換關(guān)系，再通過梁單元形函數(shù)建立投影點與主節(jié)點之間的位移轉(zhuǎn)換關(guān)系，通過兩次位移轉(zhuǎn)換可以建立梁單元內(nèi)任意節(jié)點位移與其主節(jié)點位移之間的轉(zhuǎn)換關(guān)系，進而得到用于模型降階的減縮基向量。具體步驟如下：

結(jié)構(gòu)上任意一點與其投影點之間的位移轉(zhuǎn)換關(guān)系如式(5)所示：

uj=Rjqi

(5)

展開形式為：

(6)

ujx=qix1+yqix2+zqix3+y2qix4+yzqix5+z2qix6=Fτjqixτ

ujy=qiy1+yqiy2+zqiy3+y2qiy4+yzqiy5+z2qiy6=Fτjqiyτ

ujz=qiz1+yqiz2+zqiz3+y2qiz4+yzqiz5+z2qiz6=Fτjqizτ

(7)

位移轉(zhuǎn)換矩陣Rj為：

(8)

投影點位移與其所在梁單元主節(jié)點位移的轉(zhuǎn)換關(guān)系為：

(9)

表2 梁單元形函數(shù)

由式(5)及式(9)，可知結(jié)構(gòu)上任意一點與其所在梁單元主節(jié)點之間的位移轉(zhuǎn)換陣為：

(10)

梁單元i內(nèi)所有節(jié)點的位移與主節(jié)點的位移轉(zhuǎn)換關(guān)系為：

例如：在學習乙酸乙酯的制備時，通過動畫模擬的技術(shù)可以清晰直觀的把乙酸和乙醇反應時的斷鍵、成鍵部位展示出來，學生會對這類多原子的化合物的反應機理產(chǎn)生更深入的了解，提高了化學實驗教學的教學效率

(11)

(12)

3 數(shù)值算例

通過兩個具體算例的頻率分析，比較本文提出的模型降階方法與其他模型降階方法的分析結(jié)果，證明本方法的有效性。

3.1 網(wǎng)格加筋圓柱殼算例

算例1為圖2所示的正置正交網(wǎng)格加筋圓柱

殼，其中，橫向筋條11根，軸向筋條4根。筒殼長10020mm，外半徑為1000mm，蒙皮厚10mm，筋條高度為50mm，筋條厚20mm。采用鋁合金材料，彈性模量為73GPa，泊松比為0.3，密度為2.7×10-3g/mm3。采用3-D實體單元(solid45)建立有限元模型，并沿筋條高度及厚度方向局部加密，整體網(wǎng)格密度為100mm，共有17360個節(jié)點，52080個自由度。

為對比不同梁單元個數(shù)對降階模型計算精度的影響，本文利用B4梁單元建立投影點位移與其所在梁單元主節(jié)點之間的位移轉(zhuǎn)換關(guān)系，采用3階多項式(N=3)建立截面任意一點位移與投影點位移的轉(zhuǎn)化關(guān)系。表3給出了固支-自由邊界條件下的不同梁單元個數(shù)(Nm)對網(wǎng)格加筋圓柱殼結(jié)構(gòu)整體頻率的影響。

表3 固支-自由狀態(tài)下不同梁單元個數(shù)對整體頻率的影響

同樣以圖2所示網(wǎng)格加筋柱殼為例，表4討論了多項式階數(shù)(N)對降階模型頻率計算精度的影響，假定降階模型沿軸向方向劃分20個B4梁單元。為說明本文提出方法的有效性，采用模態(tài)置信因子對降階模型和原結(jié)構(gòu)的模態(tài)作比較。模態(tài)置信因子MAC[17]定義為：

(13)

其中，φA代表原模型計算得到的模態(tài)，φR代表降階模型的模態(tài)。一般情況下，當MAC值大于0.7時，認為二者相似性較好；當MAC值小于0.5時，相似性較差。

表4 固支-自由狀態(tài)下多項式階數(shù)對降階模型頻率的影響(單位：Hz)

由表4可知，降階模型的計算精度能滿足工程設計的需要，得到的整體模態(tài)與原模型的模態(tài)具有極高的相似性，降階模型合理地保留了原模型的特征，計算效率有很大提高。此外，本文提出的模型降階方法不僅能捕捉到原模型的整體模態(tài)，還能獲得原結(jié)構(gòu)某些局部模態(tài)，如當多項式階數(shù)N=7時，降階模型得到的局部頻率計算誤差為7.05%，計算時間為原精細模型的43%。進一步提高多項式階數(shù)至N=8時，局部頻率的計算誤差為5.21%，計算時間為精細模型的63%。圖3～圖6給出部分模態(tài)對比圖。因此，本文提出的模型降階方法進一步拓寬了降階模型的適用范圍，為設計人員提供了一個簡單有效的計算模型。

3.2 含開口網(wǎng)格加筋圓柱殼算例

為滿足適用性和維修性要求，在運載火箭結(jié)構(gòu)上要設計各種開口，含開口網(wǎng)格加筋圓柱殼發(fā)生扭轉(zhuǎn)變形時，截面上各點沿軸線方向有不同的位移，產(chǎn)生了翹曲變形，Vlasov[18]提出了薄壁桿件約束扭轉(zhuǎn)時截面變形規(guī)律，認為橫截面上各點的翹曲軸向變形沿桿件橫截面按各點的扇形坐標分布，而對于復雜結(jié)構(gòu)，計算扇形坐標是非常困難的。本文通過數(shù)值方法[12]實現(xiàn)截面變形的精確描述，該方法通過對模型施加單位扭矩，利用靜力分析計算得到各個截面上的節(jié)點位移，進而通過Schmidt正交化引入翹曲基向量。具體步驟如下：

梁單元內(nèi)所有節(jié)點的變形可表示為：

(14)

(15)

梁單元內(nèi)各主節(jié)點的位移分量可表示為：

(16)

(17)

(18)

將式(18)代入式(12)得到考慮翹曲變形的位移轉(zhuǎn)換陣T。

算例2為圖7所示的含開口網(wǎng)格加筋圓柱殼，結(jié)構(gòu)尺寸及筋條分布情況和算例1相同，矩形開口位于軸向1/2位置處，開口長度L=1980mm，寬度H=692.8mm，采用鋁合金材料，彈性模量為73GPa，泊松比為0.3，密度為2.7×10-3g/mm3。有限元模型共有節(jié)點16912個，自由度為50736個。

表5給出了引入翹曲基向量前后的降階模型頻率計算結(jié)果與ANSYS計算結(jié)果的對比，其中第1、2、4、5階為彎曲模態(tài)，第3階為扭轉(zhuǎn)模態(tài)，降階模型軸向仍采用20個B4梁單元。

表5 固支-自由狀態(tài)下含開口網(wǎng)格加筋圓柱殼降階模型頻率計算結(jié)果(單位：Hz)

由表5可看出，僅通過增加多項式階數(shù)可提高降階模型的計算精度，當N=7時，彎曲頻率的計算誤差最高為3.07%，但扭轉(zhuǎn)頻率的計算誤差仍高達18.21%，這是由于降階模型不能精確描述截面的扭轉(zhuǎn)變形，導致扭轉(zhuǎn)剛度增大。引入截面翹曲基向量后，扭轉(zhuǎn)頻率的計算誤差降為0.56%，降階模型的計算精度大幅提高，計算時間為精細模型的34.1%。圖8給出了考慮翹曲變形前后降階模型得到的扭轉(zhuǎn)模態(tài)，對比圖中紅色區(qū)域的變形可看出，引入翹曲基向量的降階模型能更好地描述具有大開口結(jié)構(gòu)的實際變形。

4 結(jié)論

本文針對網(wǎng)格加筋圓柱殼結(jié)構(gòu)提出一種新的模型降階方法，該方法將復雜結(jié)構(gòu)劃分為若干個梁單元，將每個梁單元的節(jié)點作為主節(jié)點，利用多項式函數(shù)建立結(jié)構(gòu)任意節(jié)點與其所在軸線上投影點之間的位移轉(zhuǎn)換關(guān)系，再利用梁單元形函數(shù)建立投影點與主節(jié)點之間的位移轉(zhuǎn)換關(guān)系，進而實現(xiàn)模型降階。針對含開口網(wǎng)格加筋結(jié)構(gòu)，文中通過數(shù)值方法引入翹曲變形，提高降階模型扭轉(zhuǎn)頻率的計算精度。

本文提出的模型降階方法不需要大量的矩陣運算，計算效率高，可利用商業(yè)軟件作為一個黑箱來實現(xiàn)，并通過調(diào)整多項式階數(shù)及梁單元的個數(shù)控制降階模型的計算精度。文中首先以網(wǎng)格加筋柱殼為例驗證該降階方法的有效性，計算結(jié)果表明該降階模型合理地保留了原結(jié)構(gòu)的特征，不僅能用于計算結(jié)構(gòu)的整體頻率，還可捕捉到結(jié)構(gòu)的某些局部模態(tài)。其次，文中將該方法應用于含開口網(wǎng)格加筋柱殼的頻率分析中，通過數(shù)值方法引入翹曲變形提高降階模型扭轉(zhuǎn)頻率的計算精度。因此，本文提出的模型降階方法擴展了降階模型的應用范圍，為設計人員提供了一個簡單有效的計算模型。

本文提出的模型降階方法可以很方便建立不同約束條件下降階模型，在此基礎(chǔ)上，研究多工況、多邊界條件下結(jié)構(gòu)的動力特性分析及動力優(yōu)化設計是下一步開展的工作。

[1] Hao P, Wang B, Tian K, et al. Efficient optimization of cylindrical stiffened shells with reinforced cutouts by curvilinear stiffeners[J]. AIAA Journal, 2016, 54 (4):1-14.

[2] Wang B, Tian K, Hao P, et al. Numerical-based smeared stiffener method for global buckling analysis of grid-stiffened composite cylindrical shells[J]. Composite Structures, 2016, 15(2): 807-815.

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

[4] 潘忠文,王小軍,馬興瑞，等. 基于梁模型的蒙皮加筋結(jié)構(gòu)縱橫扭一體化建模研究[J].中國科學: 技術(shù)科學, 2014, 44 (5) :517-524.

[5] Cheng G D ,Cai Y W, Xu L. Novel implementation of homogenization method to predict effective properties of periodic materials[J]. Acta Mechanica Sinica, 2013, 29 (4):550-556.

[6] Cai Y W, Xu L, Cheng G D. Novel numerical implementation of asymptotic homogenization method for periodic plate structures[J]. International Journal of Solids and Structures, 2014, 51 (1) :284-292.

[7] 蔡園武. 周期性板結(jié)構(gòu)的漸近均勻化方法及微結(jié)構(gòu)優(yōu)化 [D]. 大連：大連理工大學. 2014.

[8] Zheng S F, Ding H. A model reduction method for dynamic analysis based on deformation modification and local rigid body mode[J]. Mechanics in Engineering, 2008, 30 (9):31-35.

[9] 鄭淑飛,丁樺. 基于變形修正的局部剛體化動力模型簡化方法[J].力學與實踐, 2008, 30 (6) :31-34.

[10] 王文勝,程耿東,郝鵬. 基于超梁降階模型的蒙皮加筋圓柱殼頻率分析[J].固體火箭技術(shù), 2015, 38 (3) :401-406.

[11] Cheng G D, Wang W S. Fast dynamic analysis of complicated beam-type structure based on reduced super beam model[J]. AIAA Journal, 2014, 52 (5) :952-963.

[12] Wang W S, Cheng G D, Li Q H. Fast dynamic performance optimization of complicated beam-type structures based on two new reduced physical models[J]. Engineering Optimization, 2013, 45 (7):835-850.

[13] Guyan R J. Reduction of stiffness and mass matrices[J]. AIAA Journal, 1965, 3 (2) :380.

[14] Irons B. Structural eigenvalue problems-elimination of unwanted variables[J]. AIAA Journal, 1965, 3 (5):961-962.

[15] Carrera E, Miglioretti F, Petrolo M. Computations and evaluations of higher-order theories for free vibration analysis of beams[J]. Journal of Sound and Vibration, 2012, 331 (19) :4269-4284.

[16] Carrera E, Zappino E, Filippi M. Free vibration analysis of thin-walled cylinders reinforced with longitudinal and transversal stiffeners[J]. Journal of Vibration and Acoustics, 2013, 135(1):011019.

[17] Allemang R J, Brown, D L. A correlation coefficient for modal vector analysis[M]. Proceedings of the Proceeding of the 1st international modal analysis conference, Orlando, FL. Union College, Schenectady, NY, 1982:110-116.

[18] Vlasov V Z. Thin-Walled elastic beams(2nd edn)[M].Moscow, 1959:57-162.

Frequency Analysis of Stiffened Cylinder Based on NovelReduced Order Model

LI Yu-wei, WANG Bo, CHENG Geng-dong

(State Key Laboratory of Structural Analysis for Industrial Equipment, Department of Engineering Mechanics,Dalian University of Technology, Dalian 116024, China)

A novel reduced order model is proposed to determine the natural frequencies of the stiffened cylinder in this study, which is established by using a reduction basis along with the polynomial interpolation function and shape function of beams. The basic idea is to convert the displacements of FEM nodes to the master nodes with a few generalized DOFs. The proposed ROM can gain a significant reduction of computational cost without sacrifice of accuracy. Typical stiffened cylinder are studied and the outcomes are validated by benchmark studies. The results show that the reduced order model has the ability to identify shell lobe-type modes, which can be served as an efficient alternative model during initial design.

Polynomial interpolation function; Shape function of beams; Model reduction; Stiffened cylinder; Frequency analysis

2017-03-27；

2017-04-27基金項目：國家“973”計劃(2014CB049000)；國家自然科學基金(11372062, 11402049)

李玉韋(1989-)，男，在讀博士，主要研究方向：結(jié)構(gòu)動力模型降階及優(yōu)化。E-mail:liyuwei@mail.dlut.edu.cn通訊作者：王博(1978-)，男，教授，主要從事結(jié)構(gòu)多學科優(yōu)化、航天先進材料與結(jié)構(gòu)設計和輕質(zhì)結(jié)構(gòu)抗沖擊分析與設計。E-mail:wangbo@dlut.edu.cn

V415

A

2096-4080(2017)01-0041-08