用于光滑/非光滑殼的穩(wěn)定化新型協(xié)同轉(zhuǎn)動(dòng)4 節(jié)點(diǎn)四邊形殼單元

李忠學(xué)，胡萬波

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

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

在進(jìn)行殼單元有限元分析時(shí)，往往會(huì)遇到閉鎖現(xiàn)象[5?6]。閉鎖現(xiàn)象的實(shí)質(zhì)是單元的假設(shè)位移模式無法準(zhǔn)確地描述實(shí)際位移模式，導(dǎo)致單元?jiǎng)偠绕?，位移的?jì)算值偏小[7]。在殼單元有限元計(jì)算中，剪切閉鎖和膜閉鎖的影響較為嚴(yán)重[8]。對于以受彎為主的單元，當(dāng)其膜應(yīng)力很小或接近于零時(shí)，采用數(shù)值積分將會(huì)導(dǎo)致膜閉鎖問題，使單元的膜剛度被夸大而使求解失真。對于薄壁梁元或薄板單元、薄殼單元，在采用數(shù)值積分計(jì)算單元?jiǎng)偠染仃嚂r(shí)，常會(huì)由于剪切閉鎖而導(dǎo)致單元剪切剛度被夸大，從而使得到的結(jié)構(gòu)變形偏小[9]。

解決閉鎖問題常采用2 類方法：一類是通過改變積分方案來消除閉鎖，如縮減積分法；另一類是改變應(yīng)變能表達(dá)式，如各種假定應(yīng)變法[10?14]?？s減積分法有簡單、精度高、計(jì)算效率高等優(yōu)點(diǎn)，但由于使用的積分點(diǎn)少于高斯積分所需的積分點(diǎn)數(shù)，在求解過程中，單元?jiǎng)偠染仃嚨闹扔袝r(shí)會(huì)出現(xiàn)缺陷，導(dǎo)致無法捕捉到某些變形模式，從而導(dǎo)致求解失敗。

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

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

1 單元的新型協(xié)同轉(zhuǎn)動(dòng)法描述

單元協(xié)同轉(zhuǎn)動(dòng)框架如圖1 所示。初始構(gòu)形下，四邊形兩條對角線的方向矢量v130和v240的定義如下：

圖1 協(xié)同轉(zhuǎn)動(dòng)框架Fig. 1 Description of the co-rotational framework

2 局部坐標(biāo)系下的單元公式

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

3 穩(wěn)定化的單點(diǎn)積分法

圖2 剪切應(yīng)變插值關(guān)聯(lián)點(diǎn)分布圖Fig. 2 Interpolation points for shear strain

通過混合張量插值應(yīng)變法[28]得到自然坐標(biāo)系下的單元剪切應(yīng)變：

4 整體坐標(biāo)系下的單元公式

5 算例分析

將引入物理穩(wěn)定化零能模態(tài)控制法的4 節(jié)點(diǎn)四邊形殼單元簡記為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作用。利用結(jié)構(gòu)和荷載的對稱性，取殼的1/8 結(jié)構(gòu)進(jìn)行分析。

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

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

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

5.2 懸臂梯形板

圖5 所示的懸臂梯形板，一端固支，另一端受到平面內(nèi)均布剪切荷載。結(jié)構(gòu)幾何尺寸如下：長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 個(gè)PO4Q 單元和ANSYS[32]有限元軟件采用8×8、16×16 和32×32 個(gè)shell181 單元的計(jì)算結(jié)果。表1給出了荷載參數(shù)為1 時(shí)采用16×16、32×32 個(gè)PO4Q單元和shell181 單元的計(jì)算結(jié)果?？梢?，PO4Q 單元略微偏硬，但采用32×32 個(gè)PO4Q 單元與shell181單元的計(jì)算結(jié)果的相對誤差僅5.40%，PO4Q 單元的收斂性和計(jì)算精度基本令人滿意。

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

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

5.3 直角形折板

圖7 所示的直角形折板，一端完全約束，另一端自由，并施加豎向均布荷載。結(jié)構(gòu)幾何尺寸為：長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

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

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

5.4 鐮刀形殼

如圖9 所示的鐮刀形殼，結(jié)構(gòu)左端固支，在右側(cè)自由端的中點(diǎn)A 處作用一平面內(nèi)集中荷載。結(jié)構(gòu)的幾何尺寸為：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 的單元網(wǎng)格計(jì)算該鐮刀形殼自由邊中點(diǎn)的位移。其結(jié)果和Chró?cielewski 等[34]采用基于合應(yīng)力的半混合殼單元SEMe9、基于位移-轉(zhuǎn)動(dòng)的殼單元CAMe9 以及標(biāo)準(zhǔn)退化殼單元SELe9 得到的結(jié)果吻合較好(見圖10)。

圖10 自由邊中點(diǎn)的荷載-位移曲線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 結(jié)論

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