基于MATLAB的有限元結構分析

浮 濤

（鄭州市交通規劃勘察設計研究院，河南 鄭州 450000）

0 引言

有限元法是當今工程分析中應用最為廣泛的數值分析方法。在以往的結構有限元分析編程實踐中，通常采用BASIC、FORTRAN、C、C++語言等進行編程，相比較而言，MATLAB 被視為有效快捷的工程分析工具，其優點是可將編程計算、數據分析及數據圖形化成功地集成在工作環境中，被廣泛應用于多領域。使用者可編譯M 文件或者運用軟件自身提供的強大工具箱來進行動態仿真，可求解出一系列工程問題，且求出的數值結果具有可視化、圖形功能完善等特點。

在求解普通結構的工程分析中，使用較多的二維有限元是平面梁單元、平面桁架單元和平面剛架單元，三種單元都具有兩種坐標系，分別為總體坐標系和局部坐標系。

平面梁單元與平面剛架單元的參數有E、I、A、L（E代表彈性模量，I代表截面慣性矩，A代表截面面積，L代表單元長度），每個平面梁單元與平面剛架單元均有2個節點，并且相對于總體坐標系的X軸正向逆時針的傾斜角為θ。其中平面梁單元有4 個自由度（每個節點有2 個自由度即1 個位移自由度和1 個轉角自由度）；平面剛架單元有6個自由度（每個節點有3個自由度即2個位移自由度和一個轉角自由度）[1-3]。

平面桁架單元的參數有E、A、L（E代表彈性模量，A代表截面面積，L代表單元長度），每個平面桁架單元有2個節點，并且相對于總體坐標系的X軸正向逆時針的傾斜角為θ。平面桁架單元有4 個自由度（每個節點有2個自由度即2個位移自由度）。

上述平面單元均約定位移向上為正，轉角逆時針方向為正。

本文采用MATLAB 進行有限元分析的方法，并對某平面桁架結構、平面剛架結構進行實例分析研究。

1 MATLAB有限元分析

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

結構的離散化，寫出單元剛度矩陣，集成整體剛度矩陣，引入邊界條件，求解方程組，結果處理。

1.2 編寫M函數文件

1.2.1 平面桁架單元的M函數文件

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

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

（3）PTEForce（E,A,L,θ,u）—計算單元節點力矢量。（4）PTEStress（E,L,θ,u）—計算單元應力。

1.2.2 平面剛架單元的M函數文件

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

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

（7）PFEForces（E,A,I,L,θ,u）—計算平面剛架單元的節點力矢量。

（8）PFEADiagram（f,L）—繪制節點力矢量為f和長度為L的單元軸力圖。

（9）PFESDiagram（f,L）—繪制節點力矢量為f和長度為L的單元剪力圖。

（10）PFEMDiagram（f,L）—繪制節點力矢量為f和長度為L的單元彎矩圖。

在MATLAB軟件集成環境中，若進行平面梁單元結構分析，使用（1）～（4）這4個M函數文件，根據實際結構給出函數所需要的參數直接調用即可；若進行平面桁架單元結構分析，使用（5）～（10）這6個M函數文件，根據實際結構給出函數所需要的參數直接調用即可[4]。

MATLAB 軟件的使用方法，本文不再詳細敘述，用MATLAB 軟件寫出矩陣方程和在MATLAB 集成環境中編寫并使用M函數文件請參閱文獻[5]。

2 實例分析

2.1 平面桁架結構

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

圖1 桁架結構圖

根據M函數文件PTEStiffness（）求出單元的剛度矩陣k1、k2、k3、k4，該結構有四個節點，所以整體剛度矩陣為8×8。由于該結構有四個桁架單元，故需要四次調用函數文件PTAssemble（），求得整體剛度矩陣。

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

建 立 單 元 位 移 矢 量u1、u2、u3和u4，然 后 調 用PTEForces（）求解各桿軸力，調用PTEStress（）求解各桿應力，見表1所示。

表1 各單元軸力及應力[4]

2.2 平面剛架結構

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

圖2 剛架結構圖

根據M函數文件PFEStiffness（）求出單元的剛度矩陣k1、k2、k3，該結構有四個節點，所以整體剛度矩陣為12×12。由于該結構有三個剛架單元，故需要三次調用函數文件PFAssemble（），求得整體剛度矩陣。

被約束的自由度有：U1x=U1y=R1z=0；U4x=U4y=R4z=0。則活動自由度編號為4，5，6，7，8，9。活動自由度對應的節點載荷F4=-30000N，F5=0，F6=0，F7=0，F8=0，F9=-20000kN·m，采用高斯消去法進行求解，最后求得節點2 的水平位移-0.0611m（向左），豎向位移0，轉角位移0.0078rad（逆時針）；節點3 的水平位移-0.0610m（向左），豎向位移0，轉角位移0.0043rad（逆時針）。

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

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

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

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

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

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

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

3 結束語

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