基于帶面內轉角自由度四節點平板殼單元的板殼非線性分析

文穎，戴公連，曾慶元

(中南大學 土木工程學院，湖南 長沙，410075)

薄板、薄殼結構具有較高的面內強度、結構自重輕、便于工場預制以及現場拼裝等優點，在大型土木工程結構，例如大跨度橋梁主梁(板梁、箱梁、組合梁)，大跨度屋蓋結構以及近海鉆井平臺等獲得廣泛應用，受到結構工程師們普遍關注。其設計計算理論已從基于微分方程的經典解析法發展到以有限單元法為代表的數值分析方法[1]。板殼有限元，與桿系有限元一樣，都是有限元理論體系的重要組成部分[2]。平板殼單元作為一類應用較多的板殼單元，歷來是研究的熱點。根據單元變形特點，這些研究成果可集中分成2種類型：一類是反映單元面內變形的膜位移分析，另一類則是單元面外撓曲變形分析。國內外學者在板殼單元膜位移分析領域做了大量工作，龍馭球等[3]總結了這些理論研究成果。學術界普遍認為單純按節點線位移插值得到單元膜位移場一般會引入虛假應變而影響計算結果精度。雖然采用縮減積分，加密有限元網格等手段可以一定程度上彌補上述不足，但也帶來多余零能模式，增加建模工作量等問題。目前公認的解決方案是考慮節點面內轉角自由度的貢獻[2,4],這樣既可以在保證單元具有同等位移精度條件下減少單元自由度數又可解決組裝由共面單元離散的板殼結構剛度矩陣時由于各單元面內轉角自由度的缺失而導致整體剛度矩陣病態問題。為此，Allman[5]提出了平面膜元的二階位移場理論。朱菊芬等[6]基于Allman方法將單元膜位移表示成由節點線位移插值得到的基本位移場和考慮節點旋轉自由度影響的高次位移場的疊加，以避免剪切閉鎖與零能位移模式。Long等[7]基于廣義協調元理論提出一系列高性能膜單元，例如 GT9，GT9M，GR12和GQ12等。馮仲齊等[8]通過引用連續介質力學中旋轉度的概念構造邊界位移按三次變化的 Q20單元，便于與一維梁單元的銜接。夏桂云等[9]將節點角位移作為旋轉自由度用有限條帶構造單元，具有列式簡單、物理概念清晰等特點，可是節點剛性假定增加了單元剛度。此外，也有學者提出單元線位移與轉角變位獨立插值以減輕單元“過剛”現象[4,10]。平板殼單元的面外撓曲分析通常分為2種類型[10]，即不考慮單元橫向剪切變形的 Kirchhoff型單元和計入剪切變形影響的Reissner-Mindlin(R-M)型單元。R-M板殼元雖然只要求C0類連續，但在實際應用時容易產生剪切閉鎖。Kirchhoff型單元因為列式簡單明了，計算效率高而在薄板、殼結構分析方面獲得廣泛應用。國內外學者在構造非協調與協調 Kirchhoff型薄板彎曲單元方面進行了大量研究。本文作者在總結前人工作的基礎上，提出一種面內、面外變形協調的四節點平板殼單元。采用有限條帶描述單元整體位移場。其中，膜位移場由包含節點轉角自由度在內的節點位移參數沿縱、橫向條帶順次插值得到，橫向彎曲變位則表示為雙三次Hermite多項式乘積形式。該單元模型確保位移及應變分量在單元內部和相鄰單元之間連續，物理概念明確，分析精度高。通過選取矩形薄板彈塑性分析和圓柱殼體大變位、大轉動2個經典算例，說明本文提出的板殼單元位移模型可以有效解決一般板殼結構穩定極限承載力分析問題。

1 四節點平板殼單元列式

1.1 基本假定

(1) 單元面外撓曲變形滿足Kirchhoff假定，即單元截面法線在單元變形后仍為直線且保持與單元中面垂直，忽略橫向剪切變形對單元撓度影響。其次，忽略單元纖維沿截面厚度方向擠壓變形。

(2) 基于Von-Karman大撓度、小應變理論，大變位狀態下單元應變只計入撓度w對面內坐標x，y直至二階非線性偏導數的影響，忽略面內位移u，v對坐標x，y非線性偏導數的影響。

1.2 單元空間位移場

首先選取局部坐標系xyz作為描述單元變形的參考系，如圖1所示。坐標系原點O位于單元任意一節點i，x軸和y軸均通過單元的中面且分別與單元邊界ij和邊界 il重合，z軸垂直于單元中面且與x軸和 y軸滿足右手螺旋法則。單元每個節點(角點)有 7個自由度：3個空間線位移自由度u，v，w，3個空間轉動自由度θx，θy，θz和一個扭轉率自由度κ，這樣單元共有28個自由度。本文中行向量用 表示，列向量用｛｝表示，矩陣用[]表示，下標“，”表示對緊跟著的變量求導數，而變量出現次數就是求導的 階數。

根據基本假定 1且忽略單元沿厚度方向擠壓變形，單元內部任意一點的位移可以表示為：

式中：us，vs和ws分別表示單元中面內任意點沿x，y，z坐標線位移；ws,y=θx，-ws,x=θy則分別表示單元截面繞x，y坐標軸旋轉角位移。

為了確定單元中面線位移 us(x，y)，可假想將單元沿y方向離散為連續的橫向條帶。條帶上各點位移

圖1 四節點平板殼單元Fig. 1 A four-node flat shell element

認為沿x方向線性變化，同時假定條帶端點位移沿y方向也滿足線性變化規律，即有

對于單元中面線位移 vs(x，y)，也可以假想沿 x方向單元被劃分成連續分布的條帶。條帶上各點y方向位移沿著條帶方向按3次Hermite多項式變化，而條帶端點y方向位移又可沿x方向線性插值得到。因此，

式中：θzi，θzj，θzk，θzl表示單元節點面內旋轉自由度，這些轉角自由度對單元面內位移的影響通過式(7)表示。這樣既可保證單元邊界位移相互協調，又可避免單元節點“過剛”而影響計算結果精度。

單元中面面外撓曲位移ws(x，y)可以通過沿x和y方向分別3次Hermitian插值得到[11]，即有

需要注意的是，上述單元空間位移模式具有 C1連續性，即在單元邊界上線位移和轉角是協調的；通過引入單元節點面內自旋自由度，克服了一般板殼單元處理單元面內旋轉剛度不便的缺點，因此具有較好的分析精度和較強的適用性。

1.3 單元應變度量

由前述基本假定 1可知，單元應變分量εxz=εzx=εyz=εzy=εzz=0，非零應變分量包括 εxx，εyy，εxy和εyx。又因為剪應力互等εxy=εyx，剩下3個非零獨立應變分量：正應變 εxx，εyy和剪應變 εxy。將這些應變張量分量寫成工程應變形式[12]，當不考慮單元大位移影響時，有

其中：

將式(2)，(7)和(12)代入式(17)，式(17)可縮寫成：

當分析單元大位移、大轉動問題時，需要計及幾何非線性應變的影響[13]，根據基本假定1和2，有

將式(12)代入式(22)后得：

其中：

由式(18)和式(21)，可得一般情況下單元應變度量為：

1.4 單元剛度矩陣

根據連續介質力學理論，單元任意程度變形均可表示成一系列連續有限變形過程。對于一典型有限變形又可由經典Updated-Lagrange列式描述，由此可得線性化的增量平衡方程如下[14]：

其中：D為當前增量變形步內應力-應變本構關系；tσe為單元產生增量變形之前(t時刻)Cauchy應力分量列向量；t+dtR則為單元產生增量變形后(t+dt時刻)節點外荷載列向量。將式(18)和式(21)代入式(23)可得：

α為系數矩陣且有 δe=αde。進一步整理后式(27)可寫為：

其中：

其中：h為單元厚度；Nx，Ny和Nxy為單元膜內力，由截面應力沿厚度方向積分得到；為單元自然變形剛度矩陣；當假定單元為完全彈性時，D為線彈性本構關系矩陣，當考慮材料彈塑性影響時，D則是切線本構關系矩陣。為單元幾何剛度矩陣，反映單元變形對初始節點力的影響。de和Pe分別為節點位移和等效節點外荷載。

2 板殼結構非線性分析

基于經典Updated-Lagrange列式的板殼結構非線性響應分析一般采用增量-迭代求解技術[15]。主要思路是將施加在結構上外荷載(或某些關鍵點的位移)分成離散的增量步以避免處理大位移、大轉動以及考慮材料塑性等問題帶來的不便，隨后根據弧長法等平衡路徑追蹤策略，形成自動增量累加系統。為了求得典型增量步內結構響應(內力和變形)需要引入迭代策略，這是由問題的非線性本質所決定的。結構非線性分析的核心是基于平衡迭代尋求結構非線性響應，故上述增量-迭代方法的關鍵落在典型增量步內平衡迭代上。一般來說，非線性迭代過程可分成3個階段：預測階段、修正階段和平衡校核階段。為了全文敘述完整，下面簡要描述這3個階段：

(1) 預測階段，首先依據初始幾何和物理條件，按式(29)和(30)計算單元切線剛度矩陣，又由“對號入座”法則組拼結構整體切線剛度矩陣和結構整體荷載列陣；其次，求解線性方程組得到節點位移預測值

3 數值算例

算例1 受彎方板彈塑性分析

本算例旨在說明本文提出的四節點平板殼單元可以較好地模擬塑性變形逐步發展的過程。如圖 2(a)所示，承受均布荷載 q作用的四邊簡支方形板的邊長L=1.0 m，板厚h=0.01 m，板件用材的彈性模量E=10.92 kPa，泊松比μ=0.3，屈服極限s=1 600 kPa。Owen和Hinton[16]采用 2×2的 9節點 Largrange單元和Heterosis單元網格劃分1/4板面積，運用截面分層法考慮塑性沿板厚的開展，最終得到板的極限荷載因子25.0。Voyiadjis等[17]計算了板的荷載-中央撓度曲線，如圖2(b)所示。為了便于比較，本文利用對稱性將1/4方板劃分成4×4的網格，計算方板的荷載-變位曲線(如圖2(b)所示)，并得到板的極限荷載因子為24.943。可見，本文計算結果與文獻[16]和[17]的結果吻合。

算例2 圓柱殼屋頂空間大變形分析

如圖 3(a)所示，兩對邊簡支的圓柱殼屋頂承受中央豎向集中力P的作用。柱殼屋頂的厚度為12.7 mm。該算例也是分析結構空間幾何非線性分析的經典考題，不少學者給出各自分析結果[18-19]。屋頂的材料特性參數：E=3.102 75×103MPa，μ=0.3，幾何參數屋面圓弧半徑R=2.54 m，L=0.254 m，θ=0.1 rad。出于對稱性考慮，只取1/4的屋頂進行分析，采用8×8單元網格模式。計算結果如圖3(b)所示，可見用8×8網格分析12.7 mm厚屋頂可以滿足工程計算精度需要。進一步說明采用本文提出的帶面內自旋自由度四節點平板殼單元分析板殼結構非線性問題具有較高的計算精度和求解效率。

圖2 四邊簡支矩形板荷載-中央撓度曲線Fig. 2 Load-deflection curve of simply-supported square plate

圖3 圓柱殼屋頂空間大變位分析Fig. 3 Large deflection analysis of cylindrical roof under a central point load

4 結論

(1) 提出了一種引入面內自旋自由度的四節點平板殼單元。通過給定單元縱、橫向條帶位移模式進而確定單元位移場。此種基于“條帶”的位移模型既能保證相鄰單元邊界線位移與轉角的協調一致又能與梁單元較好耦合(例如模擬帶肋板殼結構)。引入節點面內轉動自由度既能避免殼元共面后結構整體剛度矩陣奇異又可改善殼元面內剛度“過剛”而導致的薄膜閉鎖。

(2) 依據經典有限單元法例程，推導了四節點平板殼單元的小變形(自然變形)剛度矩陣和幾何剛度矩陣的計算表達式。推導思路清晰，過程簡潔，便于工程應用。

(3) 2個算例的數值結果表明：構造四節點平板殼單元不僅物理概念明確而且用于有限元分析計算結果正確、可靠。在分析一般板殼結構時，用較粗略的網格可以獲得滿足工程精度要求的結果。

[1] Hinton E, Owen D R J. Finite element software for plates and shells[M]. Swansea: Pineridge Press, 1984: 157-176.

[2] Cook R D, Malkus D S,Plesha M E, et al. Concepts and applications of finite element analysis[M]. 4th ed. New York:John & Wiley Sons, 2002: 530-580.

[3] 龍馭球, 龍志飛, 岑松. 新型有限元論[M]. 北京: 清華大學出版社, 2005: 5-7.LONG Yuqiu, LONG Zhi-fei, CENG Song. Advanced finite element method in structural engineering[M]. Beijing: Tsinghua University Press, 2005: 5-7.

[4] Hughes T J R, Brezzi F. On drilling degrees of freedom[J].Computer Methods in Applied Mechanics and Engineering, 1989,72(1): 105-121.

[5] Allman D J. A quadrilateral finite element including vertex rotations for plane elasticity analysis[J]. International Journal for Numerical Methods in Engineering, 1988, 26(3): 717-730.

[6] 朱菊芬, 鄭罡. 帶旋轉自由度C0類任意四邊形板(殼)單元[J].計算力學學報, 2000, 17(3): 287-292.ZHU Jufen, ZHENG Gang. A new 4-node C0 quadrilateral plate/shell element with drilling degree of freedom[J]. Chinese Journal of Computational Mechanics, 2000, 17(3): 287-292.

[7] LONG Yuqiu, XU Yin. Generalized conforming triangular membrane element with vertex rigid rotational freedoms[J].Finite Elements in Analysis and Design, 1994, 17(4): 259-271.

[8] 馮仲齊, 梅占馨. 有旋轉自由度的高精度四邊形單元 [J]. 應用力學學報, 2002, 19(1): 119-123.FENG Zhongqi, MEI Zhanxin. A high precision quadrilateral element with rotational degree of freedom[J]. Chinese Journal of Applied Mechanics, 2002, 19(1): 119-123.

[9] 夏桂云, 曾慶元, 李傳習. 用有限條帶構造帶旋轉自由度的矩形膜元[J]. 長沙交通學院學報, 2005, 21(1): 16-20.XIA Guiyun, ZENG Qingyuan, LI Chuanxi. A rectangular membrane element with rotational degree of freedom consisted of finite belts[J]. Journal of Changsha Communications University, 2005, 21(1): 16-20.

[10] 康瀾, 張其林. 帶旋轉自由度的四邊形平板殼單元[J]. 同濟大學學報: 自然科學版, 2009, 37(2): 164-168.KANG Lan, ZHANG Qilin. A quadrilateral flat shell element with drilling degree of freedom[J]. Journal of Tongji University:Natural Science, 2009, 37(2): 164-168.

[11] Przemieniecki J S. Theory of matrix structural analysis[M]. New York: Dover Publication, 1985: 121-122.

[12] 王勖成, 邵敏. 有限單元波基本原理和數值方法[M]. 第2版.北京: 清華大學出版社, 1996: 329-334.WANG Xucheng, SHAO Min. Principles and numerical methods of finitie element methods[M]. 2nd ed. Beijing: Tsinghua University Press, 1996: 329-334.

[13] 王榮輝. 桿、板、殼結構計算理論及應用[M]. 北京: 中國鐵道出版社, 1999: 237-239.WANG Ronghui. Bars, plates and shells: Theory and applications[M]. Beijing: China Railway Publishing House,1999: 237-239.

[14] Bathe K J. Finite element procedure[M]. New York:Prentice-Hall, 1996: 522-528.

[15] Riks E. An incremental approach to the solution of snapping buckling problems[J]. International Journal of Solids and Structures, 1979, 15(7): 529-551.

[16] Owen D R J, Hinton E. Finite elements in plasticity: Theory and practice[M]. Swansea: Pineridge Press, 1980: 370-371.

[17] Voyiadjis G Z, Woelke P. General non-linear finite element analysis of thick plates and shells[J]. International Journal of Solids and Structures, 2006, 43(7/8): 2209-2242.

[18] Yang Y B, Chang J T, Yao J D. A simple nonlinear triangular plate element and strategies of computation for nonlinear analysis[J]. Computer Methods in Applied Mechanics and Engineering, 1999, 178(3/4): 307-321.

[19] Sze K Y, Liu X H, Lo S H. Popular benchmark problems for geometric nonlinear analysis of shells[J]. Finite Elements in Analysis and Design, 2004, 40(11): 1551-1569.