DQFEM求解正交各向異性薄板單元靜力問題

鄧小環(huán)，許華南，陳水梅，黃清云

（龍巖學院 資源工程學院，福建 龍巖 364000）

1 理論基礎(chǔ)

通過采用DQFEM（微分求積有限元法）[1]，分析正交各向異性薄板單元[2]，并將所求得與有限元軟件Abaqus所得結(jié)果進行比較，進一步驗證微分求積有限元法的精確性。

本文采用薄板理論[3]，即不考慮橫向剪切變形。在彈性力學中，將兩個平行的面和垂直于這兩個平行面的棱柱或柱體所圍成的物體稱為板，如圖1所示即為矩形薄板示意圖。

圖1 矩形薄板示意圖Fig.1 Schematic diagram of a rectangular plate

利用薄板理論在分析平板彎曲時假設[4]，將平板彎曲問題簡化為二維問題，且全部應力和應變可以用板中面的撓度w表示，即：

2 DQFEM求解正交各向異性材料薄板問題

2.1 二維函數(shù)的微分求解法則

對于二維函數(shù)f（x,y），通常可用兩個一維函數(shù)的積表示[5]，即：

經(jīng)典的微分求積方法通常適合用于矩形區(qū)域[6]。根據(jù)微分求積法則，一維函數(shù)p（x）在結(jié)點xi處對坐標x的r階導數(shù)和q（y）在結(jié)點yi處對坐標y的s階導數(shù)可由結(jié)點函數(shù)值分別表示為：

為方便推導二維微分求積有限單元矩陣，引入以下矩陣和向量：

由下面的遞推關(guān)系可得到權(quán)系數(shù)為：

2.2 DQFEM求解薄板單元

對于線彈性問題，總勢能泛函包括應變能和外力勢[7]，即：

式中：ε和D分別是應變向量和材料矩陣；u是位移向量；q是荷載向量。

結(jié)構(gòu)的動能泛函為：

對于薄板彎曲問題，撓度函數(shù)也可用拉格朗日函數(shù)表示[8]為：

薄板的本構(gòu)關(guān)系：

因為εz=0，γzx=0，γzy=0，可得：

其應變-位移關(guān)系[9]為：

定義如下向量：

應用微分求積法則，可得結(jié)點應變向量為：

將式（25）和式（22）帶入式（17）和式（18），離散薄板勢能泛函為：

C=diag（Ck），其中分別為x和y方向的高斯-洛巴托積分系數(shù)，

q為結(jié)點荷載列向量，其形式同式（24）。

為了滿足薄板單元之間的C1連續(xù)條件，施加固支邊界條件需要對位移向量做一些修正。邊界結(jié)點參數(shù)應該既包括結(jié)點撓度，也包括轉(zhuǎn)角，于是有：

式中：wkx=（?w/?x）k，wky=（?w/?y）k，wkxy=（?2w/?x?y）k。

為方便書寫轉(zhuǎn)換矩陣T，此處設板x方向結(jié)點數(shù)M=2，y方向結(jié)點數(shù)N=2，則：

則轉(zhuǎn)換矩陣T為：

將式（30）帶入式（28），并進行變分，得到與結(jié)點位移向量w相對應的微分求積有限薄板單元的質(zhì)量矩陣、剛度矩陣和荷載列向量，即：

綜上可得：

2.3 算例

薄板如圖1所示，板長寬均為1，板厚h=0.01，板密度為1，彈性模量E=1，泊松比v=0.3，x方向結(jié)點數(shù)為M，y方向結(jié)點數(shù)為N。薄板對側(cè)兩邊界約束條件為簡支約束，在結(jié)點上施加均勻力為1N的力，計算板固定結(jié)點的位移值，并與有限元軟件的計算結(jié)果進行對比。如下分別為網(wǎng)格密度為10×10，20×20，30×30，40×40，50×50，100×100薄板示意圖：

圖2 薄板有限元軟件網(wǎng)格密度示意圖Fig.2 Grid density diagram of thin plate finite element software

通過DQFEM對上述4種網(wǎng)格密度進行編程分析，取薄板相同位置的結(jié)點位移值，得到如下表1至表3所表示的位移值。

表1 有限元軟件與DQFEM計算薄板結(jié)點位移值及誤差對比Tab.1 Comparison of finite element software and DQFEM for calculating plate node displacement values and errors

表2 DQFEM計算薄板結(jié)點位移值及誤差對比Tab.2 The calculation of displacement value and error comparison of thin plate node by DQFEM

表3 有限元軟件計算薄板結(jié)點位移值及誤差對比Tab.3 The displacement value and error comparison of thin plate node calculated by finite element software

由表2可得采用DQFEM計算得到的薄板結(jié)點位移值，網(wǎng)格大小對于位移值取值影響不大，即該方法可采用較少的結(jié)點數(shù)解決問題。由表3可知采用有限元軟件得到的位移值隨網(wǎng)格大小的不同而相差較大，表1則說明有限元軟件需要更密的網(wǎng)格才能接近采用DQFEM計算得到的位移值，進一步說明DQFEM方法的高效性及精確性。

3 小結(jié)

對于正交各向異性薄板單元，采用微分求積有限元法離散最小勢能原理的泛函，通過相應的Matlab程序得到所求問題的剛度矩陣，并最終得到正交各向異性薄板單元在受壓工況下的位移值，并采用有限元軟件ABAQUS中各向異性殼單元模擬了相應的工況，數(shù)值模擬結(jié)果和上述所求計算結(jié)果進行系統(tǒng)對比分析，進一步驗證了微分求積有限元法的精確性。