用于光滑/非光滑殼的穩定化新型協同轉動4 節點四邊形殼單元

李忠學，胡萬波

(浙江大學土木工程系，杭州 310058)

在單元中引入協同轉動法，可采用線性單元理論求解幾何非線性問題，而且協同轉動框架不依賴于特定的某個單元，因此可以利用現有的、性能較好的線性單元作為協同轉動單元公式的“核”，建立新的單元公式，使一些原本只能用于解決小位移與小轉角問題的單元公式可用于解決大位移、大轉角問題[1?4]。

在進行殼單元有限元分析時，往往會遇到閉鎖現象[5?6]。閉鎖現象的實質是單元的假設位移模式無法準確地描述實際位移模式，導致單元剛度偏硬，位移的計算值偏小[7]。在殼單元有限元計算中，剪切閉鎖和膜閉鎖的影響較為嚴重[8]。對于以受彎為主的單元，當其膜應力很小或接近于零時，采用數值積分將會導致膜閉鎖問題，使單元的膜剛度被夸大而使求解失真。對于薄壁梁元或薄板單元、薄殼單元，在采用數值積分計算單元剛度矩陣時，常會由于剪切閉鎖而導致單元剪切剛度被夸大，從而使得到的結構變形偏小[9]。

解決閉鎖問題常采用2 類方法：一類是通過改變積分方案來消除閉鎖，如縮減積分法；另一類是改變應變能表達式，如各種假定應變法[10?14]。縮減積分法有簡單、精度高、計算效率高等優點，但由于使用的積分點少于高斯積分所需的積分點數，在求解過程中，單元剛度矩陣的秩有時會出現缺陷，導致無法捕捉到某些變形模式，從而導致求解失敗。

自Zienkiewicz 等[15]首次提出縮減積分法以來，為了控制因縮減積分引起的零能模態，使單 元 更 穩 定，Kosloff[16]、Taylor[17]、Flanagan[18]、Liu[19?20]、王麗等[21]相繼提出了各種改進方法。MacNeal[22]基于釋放多余約束的等參原理，發展了QUAD4 單元。Hughes 等[23]在MacNeal[22]啟發下，基于縮減積分法，橫向變形采用9 節點雙二次冪形函數進行插值，轉動變量采用4 節點雙線性形函數插值，構造出一種沒有零能模態的新型4 節點四邊形雙線性等參單元。Belytschko 和Bindeman[24]基于穩定性剛度與不完全積分剛度兩者最大特征值成比例的思路，提出了穩定性參數的取值方法，選用合理的參數來控制零能模態，又不至于引起閉鎖現象。Belytschko 和Tsay[25]通過增廣B 矩陣并設置穩定性參數實現穩定化過程，提出了擾動零能模態控制法。然而，擾動法也無法完全消除潛在的零能模態，且計算結果容易受穩定性參數的影響。Belytschko 和Leviathan[26]基于Hu-Washizu 變分原理，在4 節點四邊形單元中，利用混合張量插值法來消除剪切閉鎖，通過廣義應變來控制零能模態，提出了采用物理穩定化零能模態控制法的QPH 單元。李忠學等[27]采用假定應變法構造出穩定應變來消除零能模態，發展了一種適用于光滑殼的新型協同轉動9 節點四邊形單元。

本文發展了一種新型協同轉動4 節點四邊形殼單元，采用矢量型轉動變量描述旋轉，這些變量在非線性增量求解過程中可以采用簡單的加法進行更新，在整體坐標系和局部坐標系下都能得到對稱的單元切線剛度矩陣，可提高計算效率和節省存儲空間。為消除單元中可能出現的閉鎖現象，本文采用了單點積分方案計算內力矢量和單元切線剛度矩陣，并引入Belytschko 等[26]的物理穩定化零能模態控制法，以避免因切線剛度矩陣缺秩而產生的零能模態。新型單元在4 種典型算例中均表現出良好的性能。

1 單元的新型協同轉動法描述

單元協同轉動框架如圖1 所示。初始構形下，四邊形兩條對角線的方向矢量v130和v240的定義如下：

圖1 協同轉動框架Fig. 1 Description of the co-rotational framework

2 局部坐標系下的單元公式

將應變 ε分解為3 部分，即：膜應變 εm、彎曲應變zlχ、出平面的剪切應變 γ。各應變按下式計算：

3 穩定化的單點積分法

圖2 剪切應變插值關聯點分布圖Fig. 2 Interpolation points for shear strain

通過混合張量插值應變法[28]得到自然坐標系下的單元剪切應變：

4 整體坐標系下的單元公式

5 算例分析

將引入物理穩定化零能模態控制法的4 節點四邊形殼單元簡記為PO4Q。

5.1 徑向受拉的開口圓柱殼

如圖3 所示的開口圓柱殼，兩端自由，幾何尺寸為：長度L=10.35 m ，半徑R=4.953 m，厚度a=0.094 m ，材料的彈性模量E=10.5 MPa，泊松比μ=0.3125。圓柱殼在中部受一對大小相等、方向相反的徑向集中拉力F=40 kN作用。利用結構和荷載的對稱性，取殼的1/8 結構進行分析。

圖3 徑向受拉圓柱殼Fig. 3 Pull-out of an open cylindrical shell

計算荷載作用點的撓度時，分別采用8×16、12×24 的單元網格，得到結構的荷載-位移曲線與Jiang 和 Chernuka[30]采用8×14 單元網格的4 節點ANS 退化殼單元以及Campello 等[31]采用8×16×2單元網格的6 節點三角形曲殼單元的計算結果進行比較，如圖4 所示。可見，當荷載達到 20 kN左右時，殼體出現了輕微的跳躍屈曲現象，該單元成功跟蹤到了這一過程。得到的臨界荷載和荷載-位移曲線與其他文獻[30?31]中的結果吻合較好。

圖4 徑向受拉圓柱殼在加載點處的荷載-位移曲線Fig. 4 Load-displacement curve of cylindrical shell

5.2 懸臂梯形板

圖5 所示的懸臂梯形板，一端固支，另一端受到平面內均布剪切荷載。結構幾何尺寸如下：長L=48 mm ，高H1=16 mm、H2=44 mm。材 料的彈性模量E=1.0 MPa ，泊松比μ=1/3，厚度a=1.0 mm 。均布荷載p=P/H1，其中，P=λPref，Pref=1 N。

圖5 懸臂梯形板Fig. 5 A cantilever trapezoidal plate

圖6 給出了分別采用8×8、16×16 和32×32 個PO4Q 單元和ANSYS[32]有限元軟件采用8×8、16×16 和32×32 個shell181 單元的計算結果。表1給出了荷載參數為1 時采用16×16、32×32 個PO4Q單元和shell181 單元的計算結果。可見，PO4Q 單元略微偏硬，但采用32×32 個PO4Q 單元與shell181單元的計算結果的相對誤差僅5.40%，PO4Q 單元的收斂性和計算精度基本令人滿意。

圖6 懸臂板在點A 處的荷載-位移曲線Fig. 6 Load-displacement curves at tip A of cantilever plate

表1 懸臂板在點A 處的位移(λ=1)Table 1 Displacement at tip A of cantilever plate (λ=1)

5.3 直角形折板

圖7 所示的直角形折板，一端完全約束，另一端自由，并施加豎向均布荷載。結構幾何尺寸為：長L=6 m ，寬B=3 m ，厚度a=0.03 m。材料的彈性模量E=30 MPa ，泊松比μ=0.0。自由端受均布荷載p=P/B，其中，P=λPref，Pref=3 N。

圖7 直角形折板Fig. 7 A right-angle cantilever subjected to distributed force

為檢驗單元的收斂性，分別采用3×6×2 和6×12×2 的PO4Q 單元網格對該折板進行分析，自由端的荷載-位移曲線在圖8 中給出。作為對比，還給出了采用ANSYS[32]軟件中的shell181、shell281以及Li 等[33]采用QUAD9 單元的計算結果。可見，采用PO4Q 單元網格得到的結果與ANSYS[32]和Li 等[33]的計算結果吻合較好。

圖8 直角形折板在自由端處的荷載-位移曲線Fig. 8 Load-displacement curves of right-angle cantilever

5.4 鐮刀形殼

如圖9 所示的鐮刀形殼，結構左端固支，在右側自由端的中點A 處作用一平面內集中荷載。結構的幾何尺寸為：L=10 m，R=5 m ，寬B=1 m，厚度a=0.01 m 。材料的彈性模量E=30 MPa，泊松比μ=0.3 。自由端作用的集中荷載值為：P=λPref，Pref=0.01 N。

圖9 鐮刀形殼Fig. 9 A cantilever sickle shell subjected to a lateral tip force

分別采用(15+15)×3 和(20+20)×4 的單元網格計算該鐮刀形殼自由邊中點的位移。其結果和Chró?cielewski 等[34]采用基于合應力的半混合殼單元SEMe9、基于位移-轉動的殼單元CAMe9 以及標準退化殼單元SELe9 得到的結果吻合較好(見圖10)。

圖10 自由邊中點的荷載-位移曲線Fig. 10 Load-displacement curves of cantilever sickle shell

該鐮刀形殼在不同荷載作用下的變形如圖11所示。

圖11 不同荷載作用下鐮刀形殼變形圖Fig. 11 Deformed shape of cantilever sickle shell under different levels of tip force

6 結論

本文發展了一種適用于解決光滑和非光滑板殼大位移和大轉角彈性變形問題的新型協同轉動4 節點四邊形殼單元。單元中采用矢量型轉動變量描述旋轉，采用單點積分法消除閉鎖現象，引入物理穩定化方法消除單點積分引起的零能模態。通過對2 個光滑殼和2 個非光滑殼的非線性分析，并將結果與參考文獻中的單元或ANSYS 軟件中的單元的計算結果進行對比，表明該單元已有效地減輕了閉鎖現象，并消除了零能模態的不利影響，單元的可靠性、計算效率與計算精度是令人滿意的。