基于MATLAB的有限元結(jié)構(gòu)分析

浮 濤

（鄭州市交通規(guī)劃勘察設(shè)計(jì)研究院，河南 鄭州 450000）

0 引言

有限元法是當(dāng)今工程分析中應(yīng)用最為廣泛的數(shù)值分析方法。在以往的結(jié)構(gòu)有限元分析編程實(shí)踐中，通常采用BASIC、FORTRAN、C、C++語(yǔ)言等進(jìn)行編程，相比較而言，MATLAB 被視為有效快捷的工程分析工具，其優(yōu)點(diǎn)是可將編程計(jì)算、數(shù)據(jù)分析及數(shù)據(jù)圖形化成功地集成在工作環(huán)境中，被廣泛應(yīng)用于多領(lǐng)域。使用者可編譯M 文件或者運(yùn)用軟件自身提供的強(qiáng)大工具箱來(lái)進(jìn)行動(dòng)態(tài)仿真，可求解出一系列工程問(wèn)題，且求出的數(shù)值結(jié)果具有可視化、圖形功能完善等特點(diǎn)。

在求解普通結(jié)構(gòu)的工程分析中，使用較多的二維有限元是平面梁?jiǎn)卧⑵矫骅旒軉卧推矫鎰偧軉卧?，三種單元都具有兩種坐標(biāo)系，分別為總體坐標(biāo)系和局部坐標(biāo)系。

平面梁?jiǎn)卧c平面剛架單元的參數(shù)有E、I、A、L（E代表彈性模量，I代表截面慣性矩，A代表截面面積，L代表單元長(zhǎng)度），每個(gè)平面梁?jiǎn)卧c平面剛架單元均有2個(gè)節(jié)點(diǎn)，并且相對(duì)于總體坐標(biāo)系的X軸正向逆時(shí)針的傾斜角為θ。其中平面梁?jiǎn)卧? 個(gè)自由度（每個(gè)節(jié)點(diǎn)有2 個(gè)自由度即1 個(gè)位移自由度和1 個(gè)轉(zhuǎn)角自由度）；平面剛架單元有6個(gè)自由度（每個(gè)節(jié)點(diǎn)有3個(gè)自由度即2個(gè)位移自由度和一個(gè)轉(zhuǎn)角自由度）[1-3]。

平面桁架單元的參數(shù)有E、A、L（E代表彈性模量，A代表截面面積，L代表單元長(zhǎng)度），每個(gè)平面桁架單元有2個(gè)節(jié)點(diǎn)，并且相對(duì)于總體坐標(biāo)系的X軸正向逆時(shí)針的傾斜角為θ。平面桁架單元有4 個(gè)自由度（每個(gè)節(jié)點(diǎn)有2個(gè)自由度即2個(gè)位移自由度）。

上述平面單元均約定位移向上為正，轉(zhuǎn)角逆時(shí)針?lè)较驗(yàn)檎?/p>

本文采用MATLAB 進(jìn)行有限元分析的方法，并對(duì)某平面桁架結(jié)構(gòu)、平面剛架結(jié)構(gòu)進(jìn)行實(shí)例分析研究。

1 MATLAB有限元分析

1.1 MATLAB有限元分析的基本步驟

結(jié)構(gòu)的離散化，寫出單元?jiǎng)偠染仃?，集成整體剛度矩陣，引入邊界條件，求解方程組，結(jié)果處理。

1.2 編寫M函數(shù)文件

1.2.1 平面桁架單元的M函數(shù)文件

（1）PTEStiffness（E,A,L,θ）—計(jì)算平面桁架單元（E,A,L,θ）的單剛矩陣。

（2）PTAssemble（K,k,i,j）—將連接節(jié)點(diǎn)i、j的平面桁架單元的單剛矩陣k集成到整剛矩陣K。

（3）PTEForce（E,A,L,θ,u）—計(jì)算單元節(jié)點(diǎn)力矢量。（4）PTEStress（E,L,θ,u）—計(jì)算單元應(yīng)力。

1.2.2 平面剛架單元的M函數(shù)文件

（5）PFEStiffness（E,A,I,L,θ）—計(jì)算平面剛架單元(E,A,I,L,theta)的單剛矩陣。

（6）PFAssemble（K,k,i,j）—將連接節(jié)點(diǎn)i、j的平面剛架單元的單剛矩陣k集成到整剛矩陣K。

（7）PFEForces（E,A,I,L,θ,u）—計(jì)算平面剛架單元的節(jié)點(diǎn)力矢量。

（8）PFEADiagram（f,L）—繪制節(jié)點(diǎn)力矢量為f和長(zhǎng)度為L(zhǎng)的單元軸力圖。

（9）PFESDiagram（f,L）—繪制節(jié)點(diǎn)力矢量為f和長(zhǎng)度為L(zhǎng)的單元剪力圖。

（10）PFEMDiagram（f,L）—繪制節(jié)點(diǎn)力矢量為f和長(zhǎng)度為L(zhǎng)的單元彎矩圖。

在MATLAB軟件集成環(huán)境中，若進(jìn)行平面梁?jiǎn)卧Y(jié)構(gòu)分析，使用（1）～（4）這4個(gè)M函數(shù)文件，根據(jù)實(shí)際結(jié)構(gòu)給出函數(shù)所需要的參數(shù)直接調(diào)用即可；若進(jìn)行平面桁架單元結(jié)構(gòu)分析，使用（5）～（10）這6個(gè)M函數(shù)文件，根據(jù)實(shí)際結(jié)構(gòu)給出函數(shù)所需要的參數(shù)直接調(diào)用即可[4]。

MATLAB 軟件的使用方法，本文不再詳細(xì)敘述，用MATLAB 軟件寫出矩陣方程和在MATLAB 集成環(huán)境中編寫并使用M函數(shù)文件請(qǐng)參閱文獻(xiàn)[5]。

2 實(shí)例分析

2.1 平面桁架結(jié)構(gòu)

如圖1 的桁架，已知各桿E及A均相同，且E=210GPa，A=1E-2m2，L1=400mm，L2=300mm，L3=500mm，L4=400mm，求解各桿件內(nèi)力。

圖1 桁架結(jié)構(gòu)圖

根據(jù)M函數(shù)文件PTEStiffness（）求出單元的剛度矩陣k1、k2、k3、k4，該結(jié)構(gòu)有四個(gè)節(jié)點(diǎn)，所以整體剛度矩陣為8×8。由于該結(jié)構(gòu)有四個(gè)桁架單元，故需要四次調(diào)用函數(shù)文件PTAssemble（），求得整體剛度矩陣。

被約束的自由度有：U1x=U1y=U2y=U4x=U4y=0。則活動(dòng)自由度編號(hào)為3，5，6。活動(dòng)自由度對(duì)應(yīng)的節(jié)點(diǎn)載荷F3=20000N，F(xiàn)5=0，F(xiàn)6=-25000N，采用高斯消去法進(jìn)行求解，最后求得節(jié)點(diǎn)2的水平位移3.81e-3（mm）,節(jié)點(diǎn)3的水平位移7.94e-4（mm），節(jié)點(diǎn)3的豎向位移3.125e-3（mm）。建立結(jié)構(gòu)的節(jié)點(diǎn)位移矢量U,然后計(jì)算結(jié)構(gòu)的節(jié)點(diǎn)力矢量F。由此可求得節(jié)點(diǎn)1的水平支反力-15.833kN（方向向左），垂直支反力3.125kN（方向向上），節(jié)點(diǎn)2的水平支反力20kN（方向向右），垂直支反力21.875kN（方向向上），節(jié)點(diǎn)4的水平支反力4.167kN（方向向左），垂直支反力0，滿足力的平衡條件。

建 立 單 元 位 移 矢 量u1、u2、u3和u4，然 后 調(diào) 用PTEForces（）求解各桿軸力，調(diào)用PTEStress（）求解各桿應(yīng)力，見(jiàn)表1所示。

表1 各單元軸力及應(yīng)力[4]

2.2 平面剛架結(jié)構(gòu)

如圖2的剛架，已知各桿E、I及A均相同，且A=1E-2m2，I=3E-5m4，L1=6m，L2=8m，L3=6m，E=210GPa，求節(jié)點(diǎn)2和3的位移和轉(zhuǎn)角，并采用MATLAB繪出單元2的剪力圖和彎矩圖。

圖2 剛架結(jié)構(gòu)圖

根據(jù)M函數(shù)文件PFEStiffness（）求出單元的剛度矩陣k1、k2、k3，該結(jié)構(gòu)有四個(gè)節(jié)點(diǎn)，所以整體剛度矩陣為12×12。由于該結(jié)構(gòu)有三個(gè)剛架單元，故需要三次調(diào)用函數(shù)文件PFAssemble（），求得整體剛度矩陣。

被約束的自由度有：U1x=U1y=R1z=0；U4x=U4y=R4z=0。則活動(dòng)自由度編號(hào)為4，5，6，7，8，9?；顒?dòng)自由度對(duì)應(yīng)的節(jié)點(diǎn)載荷F4=-30000N，F(xiàn)5=0，F(xiàn)6=0，F(xiàn)7=0，F(xiàn)8=0，F(xiàn)9=-20000kN·m，采用高斯消去法進(jìn)行求解，最后求得節(jié)點(diǎn)2 的水平位移-0.0611m（向左），豎向位移0，轉(zhuǎn)角位移0.0078rad（逆時(shí)針）；節(jié)點(diǎn)3 的水平位移-0.0610m（向左），豎向位移0，轉(zhuǎn)角位移0.0043rad（逆時(shí)針）。

建立結(jié)構(gòu)的節(jié)點(diǎn)位移矢量U,然后計(jì)算結(jié)構(gòu)的節(jié)點(diǎn)力矢量F。由此可求得節(jié)點(diǎn)1 的水平支反力13.187kN（方向向右），垂直支反力7.158kN（方向向上），彎矩為47.753kN·m（順時(shí)針）節(jié)點(diǎn)4 的水平支反力16.813kN（方向向右），垂直支反力7.158kN（方向向下），彎矩為54.983kN·m（順時(shí)針），滿足力的平衡條件。

表2 建立單元位移矢量u1、u2和u3，然后調(diào)用PFE Forces（）函數(shù)求出單元矢量f1、f2和f3。

表2 單元軸力、剪力、彎矩值[4]（單位：kN）

最后調(diào)用PFEADiagram（）、PFESDiagram（）以及PFEMDiagram（）可分別繪出單元的軸力圖、剪力圖和彎矩圖。單元2 的軸力圖、剪力圖和彎矩圖見(jiàn)圖3、圖4、圖5。

圖3 單元2軸力圖（單位：N）

圖4 單元2剪力圖（單位：N）

圖5 單元2彎矩圖（單位：N）

3 結(jié)束語(yǔ)

本文以平面桁架結(jié)構(gòu)和平面剛架結(jié)構(gòu)為例，利用MATLAB 進(jìn)行有限元結(jié)構(gòu)分析，經(jīng)實(shí)踐證明，該方法簡(jiǎn)便、可行、實(shí)用。本文雖僅從平面的角度為例驗(yàn)證其可行性，但其原理和方法可推廣到其他單元，如板殼、實(shí)體單元的剛度矩陣，甚至三維空間和非線性剛度矩陣之中。同時(shí)，該方法可為相關(guān)工程人員繼續(xù)深入研究提供一定參考依據(jù)。