薄殼結構的向量式有限元屈曲行為分析

王震，趙陽，楊學林

?

薄殼結構的向量式有限元屈曲行為分析

王震1, 2，趙陽1，楊學林2

(1. 浙江大學空間結構研究中心，浙江杭州，310058；2. 浙江省建筑設計研究院，浙江杭州，310006)

基于向量式有限元三角形薄殼單元的基本理論，推導力控制法和位移控制法的基本原理及在向量式有限元中的處理方式，以實現薄殼結構的屈曲和后屈曲行為分析。在此基礎上編制薄殼結構的屈曲計算分析程序，并通過算例驗證。研究結果表明：所編制的向量式有限元薄殼單元程序可以用于薄殼結構的靜、動力屈曲和后屈曲分析，驗證了理論推導和所編制程序的有效性和正確性。無需進行特殊處理，采用位移控制法可有效越過屈曲極值點，跟蹤獲得薄殼結構大位移大轉動、屈曲及后屈曲運動變形的全過程。

薄殼結構；向量式有限元；三角形薄殼單元；屈曲和后屈曲；位移控制；力控制

屈曲和后屈曲是薄殼結構分析中的一類典型問題，即研究薄殼結構在外荷載作用下產生變形以至喪失原有平衡和承載能力的過程，常伴隨有極強的幾何非線性效應。目前對于規則形狀薄殼結構的屈曲問題已有較多的理論研究成果[1?2]。根據結構失穩時平衡狀態的變化特征可分為分枝點失穩和極值點失穩，根據加載特性又可分為靜力失穩和動力失穩；傳統有限元中用于薄殼結構失穩全過程跟蹤分析的主要方法有力控制法、位移控制法和弧長法。ALAMATIAN[3]基于位移增量控制方式并結合動力松弛方法，可獲得精確的結構屈曲荷載和后屈曲行為，數值分析表明其具有較好的收斂速度和計算精度。劉正興等[4]提出交替使用載荷、位移參數確定屈曲路徑、臨界荷載的數值分析方法，避免了一般方法中存在的極值點處剛度矩陣奇異性的影響。LEE等[5]在結構屈曲路徑的預測和校正2個步驟同時引入隱式和顯式的Newton?Raphson算法，提出了改進的隱?顯和顯?隱混合弧長法，并跟蹤獲得半剛性彈塑性空間框架的后屈曲平衡路徑。ZHOU等[6]通過將不平衡荷載向量分解成2個正交向量提出了修正的弧長法，并推導新的約束方程來求解獲得當前荷載步長因子，解決了結構屈曲和材料軟化的非線性分析。梁珂等[7]基于Koiter的初始后屈曲理論和Newton法的增量迭代技術，通過在降階模型中引入攝動載荷和主路徑上的變形位移，提出能夠在非線性屈曲平衡路徑上的任意位置展開攝動求解從而實現自動跟蹤的降階方法。向量式有限元(vector form intrinsic finite element，VFIFE)[8?10]以點值描述和向量力學理論為基礎，通過牛頓第二定律描述質點運動過程來獲得結構的變形和應力。該理論本身不存在幾何線性和非線性的區分，克服了傳統有限元中容易產生的剛度矩陣奇異和非線性方程迭代不收斂問題，適用于結構大位移大轉動、機構運動及復雜結構不連續行為預測分析。對于結構屈曲問題，利用中央差分數值求解質點牛頓運動方程時，無需進行特殊處理即可越過極值點，通過結構構形的步步更新來跟蹤獲得其屈曲和后屈曲行為。本文作者基于向量式有限元三角形薄殼單元基本理論，分別采用力控制法和位移控制法來跟蹤獲得薄殼結構屈曲和后屈曲變形的全過程；在此基礎上編制薄殼結構屈曲計算分析程序，通過靜力和動力屈曲算例分析驗證理論推導和所編制程序的有效性和正確性，并對這2種控制方法的特點及其應用范圍進行探討。

1 向量式有限元薄殼單元基本理論

向量式有限元中結構質量集中于質點上，質點間則為無質量的單元連接，本文所用三角形薄殼單元為三角形CST薄膜單元和三角形DKT薄板單元的線性疊加組成。質點運動滿足平動微分方程和轉動微分 方程：

式中：和分別為質點的質量矩陣和慣量矩陣；和分別為線位移和角位移向量；為阻尼參數；和分別為質點合力和合力矩向量，下標ext和int分別表示外力(矩)和內 力(矩)。

通過中央差分數值求解質點運動式(1)時，需將其轉化為差分形式。

連續時，差分形式為

不連續時，差分形式為

式中：Δ為時間步長，?1=1/，，。

以下給出向量式有限元三角形薄殼單元的主要求解過程，詳細推導可參考文獻[11]。

1.1 單元節點純變形線(角)位移計算

采用逆向運動方法扣除單元節點全位移中的剛體平移和剛體轉動，從而獲得節點純變形線(角)位移，如式(4)所示。剛體平移可取一參考節點的總位移；剛體轉動包括面外和面內轉動，可通過單元面外(內)剛體轉動向量和對應面外(內)逆向轉動矩陣來換算獲得。對于三角形薄殼單元，質點位移的合并和分離為：變形坐標系下單元面內的節點線位移分量作為膜元部分的位移分量，面外節點線位移和角位移分量作為板元部分的位移分量。

1.2 單元節點內力(矩)計算

基于薄殼單元的變形虛功方程，由單元節點純變形線(角)位移求解單元節點內力(矩)，如式(5)所示。引入單元變形坐標系和傳統有限元單元形函數，在單元平面上求解單元節點內力(矩)向量；通過坐標轉換和面外(內)正向轉動矩陣，進而轉換為整體坐標系下內力(矩)向量；最后反向作用于質點上并進行集成即得到單元傳給質點的內力(矩)向量。對于三角形薄殼單元，質點內力的合并和分離為：節點線位移和用于膜元部分的節點內力計算，節點角位移和用于板元部分的節點內力計算。

1.3 應力應變轉換計算

僅當對結構應力、應變進行輸出時才需進行應力應變轉換計算，如式(6)所示。對于每個時間步，應力增量均是基于各自不同的變形坐標系下獲得，因而上一步末變形坐標系下應力需通過坐標轉換矩陣和面外(內)正向轉動矩陣，轉換到整體坐標系下，再轉換到下一步初變形坐標系下，以實現應力的逐步循環計算。對于三角形薄殼單元，需根據膜元部分和板元部分的應力應變性質進行應力應變的合并和分離。

2 力控制法和位移控制法的處理

力控制法和位移控制法是跟蹤獲得薄殼結構屈曲和后屈曲全過程的2種基本方法，以下給出其基本原理及在向量式有限元中的實現方式。

2.1 力控制法

力控制法以力的變化作為自變量，遵循已知結構外力求解位移反應的一般思路。在施加已知變化的外力后，通過求解質點運動式(1)可直接獲得外力作用點或其他位置點的位移變化情況，中央差分位移式(2)和(3)即為力控制方程：

2.2 位移控制法

位移控制法以位移的變化作為自變量，與力控制法相反，其遵循已知結構位移求解外力的逆向思路。在結構位移控制點處假設存在支座，施加已知變化的支座位移后，需另外通過式(2)和(3)的逆形式來獲得支座反力(即對應外力)，不連續時該逆形式即為位移控制方程：

在每個循環分析步驟中，通過控制位移控制點處支座位移()的逐步增大，求出結構在該分析步的外力()。初始位移增量可取小于單元厚度t的量級(如t/10)以滿足單元小變形大變位假定，從而可避免單元純變形過大而導致的變形計算失真。

3 程序實現及算例分析

在傳統有限元中，結構屈曲和后屈曲過程一般會出現大位移、大轉動或大應變，求解平衡方程時為考慮變形的影響應以變形后形狀為基本構形，即作為幾何非線性問題來處理。而向量式有限元方法并不存在幾何線性和非線性的區分，在每一時間子步內均為小變形大變位過程，因而均以變形前形狀為基本構形，對于結構屈曲問題無需進行特殊處理，采用顯式時間積分進行步步更新便可獲得結構的應力和變形。

基于向量式有限元三角形薄殼單元基本理論和位移(力)控制處理方法，本文采用Matlab編制了薄殼結構的屈曲計算分析程序，并分別通過靜力和動力屈曲算例分析驗證理論推導和所編制程序在薄殼結構靜力和動力屈曲分析領域的有效性和正確性，程序實現流程如圖1所示。

圖1 薄殼結構屈曲分析流程圖

3.1 柱面殼的靜力屈曲分析

兩邊簡支的柱面殼結構是分析躍越失穩問題的經典算例[7, 12]，采用位移控制法和力控制法跟蹤薄殼結構屈曲和后屈曲變形的全過程。柱面殼的半徑= 2.54 m，兩縱向邊邊長=0.508 m，厚度t=12.7 mm，曲邊弧度角為=0.2 rad；初始兩縱向邊簡支、兩圓弧曲邊自由，中心受豎向集中荷載作用；殼體的彈性模量=3.102 75 GPa，泊松比=0.3，密度=2 500 kg/m3。采用三角形薄殼單元網格，共計200個單元和121個節點，結構模型如圖2所示。分析時間步長取=2.0×10?5s，阻尼參數取=36(由文獻[13]方法計算獲得)。

(a) 三維圖；(b) 側視圖

圖3所示為外荷載隨頂部中心節點豎向位移的變化曲線，并與文獻[12]結果進行比較。由圖3可見：位移控制法結果與文獻結果符合較好，成功獲得柱面殼大位移大轉動、屈曲及后屈曲的全過程變化曲線，包括屈曲前上升段、屈曲后下降段及后屈曲承載上升段；本文方法所獲得上、下臨界極限荷載分別為2.23 kN和0.53 kN，與文獻[12]結果(2.21 kN和0.56 kN)的相對誤差僅為0.90%和5.36%。力控制法的屈曲和后屈曲上升段均與文獻結果符合，而屈曲后下降段由于荷載的逐步遞增施加表現為快速遞增跳過(無法獲得精確的臨界極限荷載)；且跳過下降段后由于該段位移變化過快，在后屈曲上升初始還存在小幅振蕩效應。

1—位移控制；2—力控制；3—文獻[12]結果。

圖4所示為幾個典型時刻柱面殼屈曲和后屈曲的變形和Mises應力分布云圖(Matlab后處理中自行編制程序繪制獲得)，以更加直觀地觀察其應力變化及變形發展情況。可見：隨著控制位移的遞增施加，柱面殼中心的變形也逐漸增大，位移為10.8 mm時達到極值點而出現屈曲；隨后開始發生快速翻轉(對應荷載?位移的下降段)，位移為19.4 mm到達下降段的極小值點；之后，柱面殼變形再次緩慢增大(對應后屈曲的上升段)。

中心節點豎向位移/mm：(a) 10.8(屈曲)；(b) 19.4(屈曲后下臨界極限)；(c) 25；(d) 30

3.2 圓柱殼的動力屈曲分析

本算例采用位移控制法分析圓柱薄殼結構在軸壓荷載作用下的動力屈曲和后屈曲問題[14?15]。圓柱殼的半徑=0.1 m，軸向長度=0.5 m，厚度t=2.0 mm，初始底端固定，頂端僅有軸向自由度并按強制性位移控制形式進行加載，使該端截面以40 m/s的速度沿軸向發生壓縮位移；殼體的彈性模量=201 GPa，泊松比=0.3，密度=7 850 kg/m3；采用三角形薄殼單元網格，共計2 560個單元和1 320個節點，結構模型如圖5所示。分析時間步長取=1.0×10?6s(對應位移增量步長為4.0×10?5 m)，阻尼參數取=0(即不考慮阻尼效應)。

圖5 圓柱薄殼模型及其網格劃分

圖6所示為頂端的總軸壓荷載隨軸向位移的變化曲線，并與ABAQUS計算結果進行比較。由圖6可見：位移控制法的前3次屈曲情況均與ABAQU符合較好；第3次屈曲之后開始出現差異，但總的變化趨勢仍基本一致。本文方法獲得的前3次屈曲荷載分別為2 286.9，6 125.0和6 888.7 kN，與ABAQUS結果(分別為2 252.7，5 565.0和6 870.1 kN)的相對誤差僅為1.50%，9.14%和0.27%；驗證了本文方法在薄殼結構動力屈曲分析中的有效性和準確性。

軸向位移/mm：(a) 0~80；(b) 0~12(前3次屈曲)1—位移；2—ABAQUS。

圖7所示為幾個典型時刻圓柱薄殼屈曲和后屈曲的變形和Mises應力分布云圖。由圖7可見：隨著控制位移的遞增施加，首先在柱殼頂端出現局部屈曲(第1次屈曲模態)，接著在柱殼底端出現局部屈曲(第2次屈曲模態)，之后的高階屈曲則開始逐步擴展至柱殼整體(整體屈曲模態)。

軸向位移/mm：(a) 0.28(第1次屈曲)；(b) 8.04(第2次屈曲)；(c) 20；(d) 40；(e) 60；(f) 80

由于本文向量式有限元方法基于質點動力分析的理論基礎，在薄殼結構動力屈曲問題中更為適用；屈曲和后屈曲過程通常伴隨有大變形大轉動，該法克服了傳統有限元中容易產生的剛度矩奇異和非線性方程迭代不收斂問題，通過結構構形的步步更新可輕易跟蹤獲得其屈曲和后屈曲行為；相對傳統有限元計算更為準確和迅速。

4 結論

1) 基于向量式有限元三角形薄殼單元的基本理論，推導了力控制法和位移控制法的基本原理及在向量式有限元中的處理方式，以實現薄殼結構的屈曲和后屈曲行為分析。

2) 編制了薄殼結構屈曲計算分析程序，并進行了算例驗證。算例分析表明所編制程序可很好地完成薄殼結構的靜、動力屈曲和后屈曲分析，驗證了理論推導和所編制程序的有效性和正確性。

3) 采用位移控制法可有效跟蹤薄殼結構運動變形的全過程，并獲得薄殼結構失穩后的下降段和精確的上、下臨界極限荷載；而力控制法無法獲得位移的下降段且在后屈曲上升初始會存在小幅度的振蕩效應。

4) 當作用位置點同時存在荷載和位移的下降段時，需引入更為高級的弧長法才可獲得準確的荷載?位移曲線，本文將進行進一步研究。

[1] TENG J G, ROTTER J M. Buckling of thin metal shells[M]. London: Spon Press, 2004: 1?41.

[2] 沙風煥. 圓柱殼的動力屈曲理論及其應用[M]. 北京: 中國建材工業出版社, 2008: 62?89. SHA Fenghuan. Theory and application of dynamic bucking of cylindrical shells[M]. Beijing: China Building Materials Press, 2008: 62?89.

[3] ALAMATIAN J. Displacement-based methods for calculating the buckling load and tracing the post-buckling regions with dynamic relaxation method[J]. Computer & Structures, 2013, 114: 84?97.

[4] 劉正興, 葉榕. 薄殼后屈曲分析的載荷?位移交替控制法[J]. 上海交通大學學報, 1990, 24(3): 38?44. LIU Zhengxing, YE Rong. Alternatively controlling technique of load and displacement for post-buckling analysis of thin shells[J]. Journal of Shanghai Jiaotong University, 1990, 24(3): 38?44.

[5] LEE K, HAN S E, HONG J W. Post-buckling analysis of space frames using concept of hybrid arc-length methods[J]. International Journal of Non-linear Mechanics, 2014, 58: 76?88.

[6] ZHOU L Y, LI Q, LI T M, et al. Improved arc-length method for solving buckling problem[J]. Journal of Southwest Jiaotong University, 2011, 46(6): 922?925.

[7] 梁珂, 孫秦, GURDAL Z. 結構非線性屈曲分析的有限元降階方法[J]. 華南理工大學學報(自然科學版), 2013, 41(2): 105?110. LIANG Ke, SUN Qin, GURDAL Z. Finite element-based order reduction method for nonlinear buckling analysis of structures[J]. Journal of South China University of Technology (Natural Science Edition), 2013, 41(2): 105?110.

[8] TING E C, SHIH C, WANG Y K. Fundamentals of a vector form intrinsic finite element: part I. Basic procedure and a plane frame element[J]. Journal of Mechanics, 2004, 20(2): 113?122.

[9] TING E C, SHIH C, WANG Y K. Fundamentals of a vector form intrinsic finite element: part II. Plane solid elements[J]. Journal of Mechanics, 2004, 20(2): 123?132.

[10] WU T Y, WANG C Y, CHUANG C C, et al. Motion analysis of 3D membrane structures by a vector form intrinsic finite element[J]. Journal of the Chinese Institute of Engineers, 2007, 30(6): 961?976.

[11] 王震, 趙陽, 胡可. 基于向量式有限元的三角形薄殼單元研究[J]. 建筑結構學報, 2014, 35(4): 64?70. WANG Zhen, ZHAO Yang, HU Ke. Triangular thin-shell element based on vector form intrinsic finite element[J]. Journal of Building Structures, 2014, 35(4): 64?70.

[12] HORRIGMOE G, BERGAN P G. Nonlinear analysis of free-form shells by flat finite element[J]. Computer Methods in Applied Mechanics and Engineering, 1978, 16(1): 11?35.

[13] 王震, 趙陽, 胡可. 基于向量式有限元的三角形薄板單元研究[J]. 工程力學, 2014, 31(1): 37?45. WANG Zhen, ZHAO Yang, HU Ke. Triangular thin-plate element based on vector form intrinsic finite element[J]. Engineering Mechanics, 2014, 31(1): 37?45.

[14] 邢小東, 侯飛. 軸向沖擊載荷下圓柱殼動力屈曲的計算機仿真[J]. 機械管理開發, 2008, 23(1): 77?78. XING Xiaodong, HOU Fei. The computer simulation on dynamic buckling of cylindrical shells under axial impulsive load[J]. Mechanical Management and Development, 2008, 23(1): 77?78.

[15] 李帥. 沖擊載荷下圓柱殼非線性動力屈曲的數值研究[D]. 大連: 大連理工大學化工機械與安全學院, 2005: 46?72. LI Shuai. Numerical research for nonlinear dynamic buckling of shell under axial impacted load[D]. Dalian: Dalian University of Technology. School of Chemical Machinery and Safety Engineering, 2005: 46?72.

(編輯 趙俊)

Vector form intrinsic finite element method for buckling analysis of thin-shell structures

WANG Zhen1, 2, ZHAO Yang1, YANG Xuelin2

(1. Space Structure Research Center, Zhejiang University, Hangzhou 310058, China;2. Zhejiang Province Institute of Architectural Design and Research, Hangzhou, Zhejiang 310006, China)

Based on the basic theory of triangular thin-shell element of VFIFE, the basic principles of force-controlled method and displacement-controlled method were derived, and the corresponding processes of VIFIFE were presented. Then buckling and post buckling behavior for thin-shell structures were analyzed. On this basis, a computer program of the buckling analysis for triangular thin-shell element was developed and numerical examples were provided. The results show that the static and dynamic buckling and post buckling analysis can be well performed for thin-shell structures by the developed program, verifying the validity and the correctness of the theoretical derivation and the computer program. Without special processing, the buckling extreme point can be acrossed effectively by adopting displacement-controlled method, and the whole deformation of large deformation and large rotation, buckling and post buckling of thin-shell structures can be tracked.

thin-shell structures; vector form intrinsic finite element; triangle thin-shell element; buckling and post buckling; displacement-controlled method; force-controlled method

10.11817/j.issn.1672-7207.2016.06.033

TU33

A

1672?7207(2016)06?2058?07

2015?06?15；

2015?09?03

國家自然科學基金資助項目(51378459)(Project(51378459) supported by the National Natural Science Foundation of China)

趙陽，教授，博士生導師，從事薄殼結構、鋼結構和向量式有限元研究；E-mail：ceyzhao@zju.edu.cn