基于新型降階模型的網格加筋圓柱殼頻率分析

李玉韋,王 博,程耿東

(大連理工大學工業裝備結構分析國家重點實驗室工程力學系,大連 116024)

基于新型降階模型的網格加筋圓柱殼頻率分析

李玉韋,王 博,程耿東

(大連理工大學工業裝備結構分析國家重點實驗室工程力學系,大連 116024)

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

多項式函數；梁單元形函數；模型降階；網格加筋圓柱殼；頻率分析

0 引言

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

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

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

1 模型降階方法

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

(1)

引入位移轉換矩陣T(r×k)(r?k)，復雜結構的位移可表示為：

U=TQ

(2)

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

(3)

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

(KR-λkkMR)φkk=0

(4)

其中，λkk為結構第kk階頻率，φkk為與其對應的特征向量。

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

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

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

結構上任意一點與其投影點之間的位移轉換關系如式(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)

位移轉換矩陣Rj為：

(8)

投影點位移與其所在梁單元主節點位移的轉換關系為：

(9)

表2 梁單元形函數

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

(10)

梁單元i內所有節點的位移與主節點的位移轉換關系為：

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

(11)

(12)

3 數值算例

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

3.1 網格加筋圓柱殼算例

算例1為圖2所示的正置正交網格加筋圓柱

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

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

表3 固支-自由狀態下不同梁單元個數對整體頻率的影響

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

(13)

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

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

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

3.2 含開口網格加筋圓柱殼算例

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

梁單元內所有節點的變形可表示為：

(14)

(15)

梁單元內各主節點的位移分量可表示為：

(16)

(17)

(18)

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

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

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

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

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

4 結論

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

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

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

[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] 潘忠文,王小軍,馬興瑞，等. 基于梁模型的蒙皮加筋結構縱橫扭一體化建模研究[J].中國科學: 技術科學, 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] 蔡園武. 周期性板結構的漸近均勻化方法及微結構優化 [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].固體火箭技術, 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-)，男，在讀博士，主要研究方向：結構動力模型降階及優化。E-mail:liyuwei@mail.dlut.edu.cn通訊作者：王博(1978-)，男，教授，主要從事結構多學科優化、航天先進材料與結構設計和輕質結構抗沖擊分析與設計。E-mail:wangbo@dlut.edu.cn

V415

A

2096-4080(2017)01-0041-08