劉慧善，朱燈林

(河海大學機電工程學院，江蘇 常州 213022)

HUGHES等[1]提出的等幾何分析方法(isogeometric analysis，IGA)以CAD模型的NURBS/B-Spline[2-3]函數(shù)作為基函數(shù)來構造力學分析中的位移場，對復雜模型進行精確建模，可以采用稀疏單元劃分獲得較高的分析精度，提高了運算效率，同時便于CAD/CAE集成，特別適用于復雜結構和板殼結構的分析和優(yōu)化設計。

板殼結構廣泛應用于各種工程實際中，如機械結構、建筑結構等。對板殼結構的力學性能分析可采用三維實體單元式殼單元進行計算。HUGHES[4]提出了利用高階多項式對板殼單元插補的方法，進行板殼結構的數(shù)學建模。但是高階多項式插補存在著效率低、建模能力差等問題。針對復雜板殼結構的等幾何分析[5-6]，DORNISCH[7]推導了等幾何分析過程中坐標系的建立方法以及坐標系的相互轉換。KANG[8]提出了裁剪單元方法，對復雜曲面模型進行了分析。KIENDL[9]提出了等幾何Kirchhoff-love殼單元，并用彎曲片方法處理多片組合殼體問題。詹雙喜[10]根據(jù)Kirchhoff-love殼理論推導了薄殼自由振動分析等幾何方法的列式。BENSOND[11]結合Kirchhoff-love和Reissner-Mindlin兩種殼單元特點發(fā)展了混合單元。

本文對Reissner-Mindlin(RM)板殼結構的等幾何分析建模和位移的插補方法進行了研究，分別采用了兩種位移插補方法：一種是根據(jù)HUGHES[4]提出的高階多項式插補方法推導的RM殼單元IGA位移插補方法，另一種是根據(jù)KANG[8]提出的位移插補方法推導的RM殼單元IGA位移插補方法。運用這兩種方法對相關實例進行了分析計算，比較了兩種插補方法的精確度和計算效率。研究結果可為板殼結構的等幾何分析研究和應用提供參考。

1 基于IGA的殼體幾何模型構造

1.1 殼體的定義

殼體結構通常采用其參考面和厚度方向的參數(shù)對其進行幾何定義，如圖1所示，其數(shù)學模型如下：

(1)

圖1 殼體定義

1.2 RM殼體坐標系

其中

式中x,ξ,x,η表示對參數(shù)坐標系求偏導。

圖2 RM殼局部坐標系(lamina[1]坐標系)

由此可得從全局坐標系到殼局部坐標系的轉換矩陣：

式中：i,j分別代表矩陣的行和列。

在傳統(tǒng)殼結構分析[14-16]中，殼內任一點的位移可由參考面沿總體坐標系x,y,z方向的3個位移分量ua,va,wa及對應的θa1,θa2繞節(jié)點坐標系轉動所確定。由于傳統(tǒng)殼單元節(jié)點位于參考面上，因此各節(jié)點的局部坐標系按文獻[1]中的方法確定即可。但是在IGA中，由于控制節(jié)點通常不在參考平面內，這就給節(jié)點處的坐標系的建立以及后續(xù)位移場的插補帶來了困難。為了在每個控制節(jié)點處建立相應的局部坐標系，本文采用如下的插補方法。

(2)

(3)

式中：np為參考面上控制節(jié)點數(shù)；Ω為參數(shù)空間積分區(qū)域。

圖3 IGA控制節(jié)點坐標系(fiber[1]坐標系)

則有：

(4)

根據(jù)式(4)可得：

(5)

2 殼體位移插補方法

(6)

(7)

(8)

第二種是根據(jù)KANG在文獻[8]中提出的位移插值方法推導得到的RM殼單元IGA位移插補方法。在插值過程中直接采用角度進行插值，然后再將轉角轉換為各個方向的位移。具體表達式如下：

(9)

3 剛度矩陣的求解

采用第一種位移插補方法時，應變位移函數(shù)矩陣定義如下:

B=[B1B2…Ben]

(10)

應變可以寫成如下形式：

(11)

對式(6)求導可得：

應變位移函數(shù)矩陣的形式如下：

(12)

其中

(13)

(14)

其中：

采用第二種位移插補方法時，應變位移函數(shù)矩陣定義與式(10)～式(12)一樣。對式(9)求導可得：

應變位移函數(shù)矩陣中，對式(13)的求解與前面所述相同，式(14)的具體計算過程如下：

4 積分簡化

由于板殼結構在厚度方向上的尺寸與其他方向相差較大，在分析過程中可能會出現(xiàn)剪切鎖死或者零勢能的情況。為避免這些情況的發(fā)生，通常采用簡化積分的方法，即在有限元剛度積分的過程中對部分自由度采用高階高斯積分，其余部分采用低階高斯積分[1]。

對于殼單元而言，在對應變函數(shù)矩陣進行積分時，通常對橫向剪切應力和膜應力積分部分進行簡化[1]，簡化方法如下：

式中：×表示不需要精確的積分;?表示需要精確積分。

5 實例分析

基于上述兩種插補方法，首先在MATLAB中分別編程對殼體結構進行靜力學分析，并將分析結果與ANSYS分析結果進行比較以證明分析結果的正確性。然后對兩種位移插補方法的分析結果進行比較，得出兩種插補方法各自的優(yōu)缺點。

5.1 懸臂薄板

如圖4所示，懸臂薄板尺寸為10m×10m×0.1m，其中一邊固定，在相對的另一邊中點處受F=1kN的集中載荷，材料彈性模量E=1×107MPa，泊松比υ=0.3。

圖4 懸臂薄板模型

圖5為懸臂板模型采用ANSYS分析得到的結果。圖6是采用第一種插補方法通過MATLAB編程分析得到的懸臂薄板位移云圖。圖7為第二種插補方法通過MATLAB編程分析得到的懸臂薄板位移云圖。圖8所示為載荷加載處位移隨著單元數(shù)變化的曲線。圖9所示為分析時間與單元數(shù)目的關系曲線。由圖5～9可知，對于平面薄板而言，兩種插補方法得到的結果基本與ANSYS分析結果一致，第二種插補方法分析效率更高。

圖5 懸臂薄板ANSYS位移云圖

圖6 第一種插補方法位移云圖

圖7 第二種插補方法位移云圖

圖8 兩種插補方法位移收斂曲線圖

圖9 兩種插補方法分析時間圖

5.2 半圓環(huán)

如圖10所示，半圓環(huán)半徑為10m，兩底端固定，圓環(huán)中心處受集中載荷F=1kN的作用，材料彈性模量E=1×107MPa，泊松比υ=0.3。

圖10 半圓環(huán)模型

圖11為半圓環(huán)模型采用ANSYS分析得到的結果。圖12為第一種插補方法通過MATLAB編程分析得到的半圓環(huán)位移云圖。圖13為第二種插補方法通過MATLAB編程分析得到的半圓環(huán)位移云圖。圖14為載荷施加處位移隨單元數(shù)目的變化曲線。圖15為分析時間與單元數(shù)目的關系曲線。由圖11～15可知，兩種插補方法在最終均收斂，且收斂結果與ANSYS接近，同樣，也是第二種插補方法求解效率更高。

圖11 半圓環(huán)ANSYS位移云圖

圖12 第一種插補方法位移云圖

圖13 第二種插補方法位移云圖

圖14 兩種插補方法位移收斂曲線圖

圖15 兩種插補方法分析時間圖

6 結束語

本文討論了RM殼單元的等幾何分析計算建模方法，分析了基于等幾何分析的兩種插補方法在收斂性、效率及在分析精度方面的差異。第一種插補方法具有較高的分析精度，但在分析過程中隨著單元數(shù)目的增多，對控制節(jié)點坐標系統(tǒng)求解所消耗的時間也隨之增多；第二種插補方法在分析過程中省去了求解控制點坐標系統(tǒng)的過程，因而插補模型二求解效率更高。

兩種插補方法在采用單元數(shù)目較少的情況下，結果有一定差異但差異不大，而在單元劃分足夠細時，兩種方法得到的結果基本一樣，但第二種插補方法的計算效率更高。因此在分析板殼結構時，為了提升分析效率可直接采用第二種位移插補方法。

參考文獻：

[1] HUGHES T J R . Isogeometric Analysis:Toward Integration of CAD and FEA[M].New York：Wiley,2009.

[2] PIEGL Les,TILLER Wayne.非均勻有理B樣條[M].趙罡，穆國旺，王拉柱,譯. 2版.北京：清華大學出版社,2010,1-314.

[3] SEVILLA R,MéNDEZ S F,HUERTA A.3D NURBS- enhanced finite element method (NEFEM)[J].Int. J. Numer. Methods Engrg,2008,88:103-125.

[4] HUGHES T J R . The Finite Element Method:Linear Static and Dynamic Finite Element Analysis[M]. Englewood Cliffs,New Jersey:Prentice-Hall Inc, 1987:383-405.

[5] 吳紫俊,黃正東,左兵權,等. 等幾何分析研究概述[J]. 機械工程學報,2015(5):114-129.

[6] 陳濤，莫蓉，萬能. 等幾何分析中Dirichlet邊界條件的配點施加方法[J]. 機械工程學報. 2012(5): 157-164.

[7] DORNISCH W.An efficient and robust rotational formulation for isogeometric Reissner-Mindlin shell elements[J].Computer Methods in Applied Mechanics and Engineering, 2016,303:1-34.

[8] KANG Pilseong.Isogeometric shape optimization of trimmed shell structures[J].Struct Multidisc Optim ,2016,53:825-845.

[9] KIENDL J.Isogeometric shell analysis with Kirchhoff-love elements[J].Computer Method in Applied Mechanics and Engineering,2009,198:3902-3914.

[10] 詹雙喜. 薄殼結構自由振動分析的等幾何方法[C]//第16屆北方七省市區(qū)力學學會學術會議.河南：鄭州大學出版社，2016:237-244.

[11] BENSOND J.Blended isogeometric shells[J].Computer Method in Applied Mechanics and Engineering,2013,255:133-146.

[12] 申文斌.張量分析與彈性力學[M].北京：科學出版社， 2016.

[13] GREEN A E, ZERNA W.Theoretical Elasticity [M]. 2nd ed. New York:Dover Publications,INC.

[14] 王勖成.有限單元法[M].北京：清華大學出版社,2008,334-417.

[15] 曾攀.有限元分析及應用[M].北京：清華大學出版社,2004,144-269.

[16] MOAVENI Saeed.有限元分析——ANSYS理論與應用 [M].王崧,譯.3版.北京：電子工業(yè)出版社,2013.