王鵬+周琪琛+姚姍姍

摘 要：本文基于平面桁架有限元分析的基本原理，利用MATLAB語言編程對有外荷載作用的平面桁架進行有限元分析，結(jié)果表明，通過MATLAB軟件對平面桁架受力分析的結(jié)果與精確解吻合。本文介紹的方法，在平面桁架有限元中具有普遍的適用性，對復雜的平面桁架結(jié)構(gòu)有限元分析有一定的參考價值。

關(guān)鍵詞：平面桁架；MATLAB；有限元分析

MATLAB是以矩陣為基本的運算單元，可以靈活地進行矩陣運算、圖形繪制、編程開發(fā)等，具有編程效率高、可移植性強、計算速度快等特點；有限元分析法是根據(jù)變分原理求解數(shù)學及物理問題的數(shù)值計算方法，它是隨著近年來計算機技術(shù)的迅速發(fā)展而得到的廣泛應用。本文以解決一個實際的平面桁架問題為例，運用有限元分析法，并利用MATLAB軟件進行編程計算來演示MATLAB軟件在平面桁架中的應用。接下來，首先介紹解決平面桁架問題的有限元分析方法。

1 平面桁架有限元分析的基本原理

在用有限元法對平面桁架受力分析中，平面桁架元是分析的基本單元，它是一種二維有限元，每個平面桁架元有2個節(jié)點和3個參數(shù)，參數(shù)分別為長度L、彈性模量E和橫截面積A，當假設桁架元與正方向總體X軸逆時針傾斜θ角時，并令C=Cos（θ），S=Sin（θ），則單元剛度矩陣可表示為：

在有限元分析中，通過單元分析，建立單元剛度矩陣k，然后再將單元剛度矩陣通過剛度集成規(guī)則集合成結(jié)構(gòu)的整體剛度矩陣K，對于一個有n個節(jié)點的結(jié)構(gòu)而言，其整體的剛度矩陣K為2n×2n的矩陣，在實際MATLAB軟件操作中，并不需要編寫函數(shù)程序，而是直接調(diào)用相應的函數(shù)即可，正是因為這樣，在用MATLAB軟件進行桁架受力分析時，可以大大提高效率，節(jié)省時間。整體剛度矩陣的函數(shù)名稱為PlaneTrussAssemble。一旦得到剛度矩陣K，就可以列出方程：

式中，U代表節(jié)點的位移矢量，F(xiàn)是結(jié)構(gòu)節(jié)點的荷載矢量，這兩個邊界條件需要手動賦值，然后利用高斯消元法便可求解上述方程組，一旦求解出為止的位移和支反力，就可以利用方程：

求解單元的節(jié)點力。式中f代表單元節(jié)點力，u是單元節(jié)點的位移矢量，最后，將單元節(jié)點力除以桁架的橫截面積便可以得到單元應力，之后就可以進行結(jié)構(gòu)校核、結(jié)構(gòu)安全性檢驗等一系列力學性實驗分析。

2 有限元法分析的一般步驟

有限元法是一種高效能、普遍使用的數(shù)值計算方法，隨著近年來計算機技術(shù)的迅速進步，有限元法得到了廣泛的使用和發(fā)展，對于有限元法分析的步驟，不同書籍的介紹不盡相同，但大體上可以分為以下幾步：

離散化域：將結(jié)構(gòu)分解為獨立的單元和節(jié)點，對于桁架和鋼架這類的離散系統(tǒng)，這一步可以省略，對于其他的連續(xù)性結(jié)構(gòu)，如板殼，這一步就顯得尤為重要，離散化的好壞直接影響到最后結(jié)果的準確性。

得到單元剛度矩陣：由上一步離散化處理的結(jié)果，寫出每個單元的剛度矩陣，這一步可以通過調(diào)用MATLAB工具箱來完成，相應函數(shù)的使用方法將在下面的實例中做出介紹。

集成整體剛度矩陣：根據(jù)上一步得到的單元剛度矩陣，通過剛度集成規(guī)則，集成結(jié)構(gòu)整體的剛度矩陣，這一步也可以直接用MATLAB工具箱來完成，

引入邊界條件：所謂的邊界條件是指位移、外加荷載、支座類型等，不同結(jié)構(gòu)的邊界條件不盡相同，所以這一步需要手動賦值，具體的操作步驟會在案例中演示。

解方程：在這一步中，需要對整體剛度矩陣進行分解，然后再用高斯消去法求解方程組，在用高斯消去時，有時候需要手動分解矩陣。

后處理：當需要得到其他信息時，如支反力、單元節(jié)點力、單元應力等，還需要進行后處理工作，這一步需要操作者掌握一定的材料力學方面的知識。

3 實例分析

為了進一步說明MATLAB在平面桁架計算中的作用，接下來通過計算一個實例，來演示有限元法分析平面桁架問題的具體操作步驟。

如圖為一平面桁架結(jié)構(gòu)，為了便于計算，

假定結(jié)構(gòu)的彈性模E=200GPa

橫截面積A=2cm2，結(jié)構(gòu)受力如圖所示：

先對結(jié)構(gòu)進行離散化處理，由于桁架結(jié)構(gòu)已經(jīng)是離散化結(jié)構(gòu)，所以我們只需要將結(jié)構(gòu)單元編號即可：

接下來多次調(diào)用MATLAB的PlaneTrussElementStiffness函數(shù)，分別生成單元剛度矩陣k1、k2、k3。然后集成整體剛度矩陣，由于該結(jié)構(gòu)有3個節(jié)點，所以整體剛度矩陣為6×6矩陣，在生成整體剛度矩陣前，要先建立一個6×6的零矩陣，零矩陣可以手動設置，也可以直接調(diào)用MATLAB的zeros（）函數(shù)直接生成相應的零矩陣，得到零矩陣后，再反復調(diào)用planeTrussAssemble函數(shù)生成整體剛度矩陣K，由于本題結(jié)構(gòu)只有3個單元，所以只需調(diào)用3次該函數(shù)即可；得到整體剛度矩陣后，就可以建立該結(jié)構(gòu)的矩陣方程[K]{U}={F}，再輸入邊界條件，本題的邊界條件為U1x=U1y=U2y=0，F(xiàn)2x=0，F(xiàn)3x=5，F(xiàn)3y=-10 。接著手動分解方程并通過MATLAB軟件用高斯消法求解方程，便可得到節(jié)點的位移，最后一步就是后處理，這一步需要看題目要求解什么未知數(shù)，如求外力，用F=K*U等式便可求解，然后調(diào)用PlaneTrussElementstress函數(shù)即可求出。下面是該實例在用MATLAB軟件求解的程序：

% 輸入?yún)?shù) E=200e6；A=2e-4；

% 計算各桿的長度

L1=4；L2=PlaneTrussElementLength（0，0，2，3）；L3=PlaneTrussElementLength（4，0，2，3）；

% 計算單元剛度矩陣，

k1=PlaneTrussElementStiffness（E，A，L1，0）；theta2=atan（3/2）*180/pi；theta3=180-theta2；k2=PlaneTrussElementStiffness（E，A，L2，theta2）；k3=PlaneTrussElementStiffness（E，A，L3，theta3）；

% 建立零矩陣 K=zeros（6，6）；

% 整體的剛度矩陣

K=PlaneTrussAssemble（K，k1，1，2）；K=PlaneTrussAssemble（K，k2，1，3）；K=PlaneTrussAssemble（K，k3，2，3）；

% 引入邊界條件，手動賦值 k=[K（3，3），K（3，5：6）；K（5：6，3），K（5：6，5：6）]；

% 外力 f=[0；-5；-10]；

% 利用高斯消去法解方程組 u=k＼f； U=[0；0；u（1）；0；u（2：3）]；

% 計算節(jié)點力矢量 F=K*U；

% 調(diào)用節(jié)點位移矢量

u1=[U（1）；U（2）；U（3）；U（4）]；u2=[U（1）；U（2）；U（5）；U（6）]；u3=[U（3）；U（4）；U（5）；U（6）]；

% 計算單元應力

sigmal=PlaneTrussElementStress（E，L1，0，u1）；sigma2=PlaneTrussElementStress（E，L2，theta2，u2）；sigma3=PlaneTrussElementStress（E，L3，theta3，u3）；

運行后的部分結(jié)果：

4 結(jié)論

實踐證明，應用MATLAB軟件解平面桁架問題簡便、易行、通用，本文只是列舉了一個簡單的例子，對于求解復雜的平面桁架問題，本文介紹的方法同樣適用。因此，我們倡議將MATLAB軟件運用到教學里，讓更多的學生了解并會用MATLAB軟件求解力學問題，以便提高學生學習效率，激發(fā)學生學習興趣。

參考文獻

[1] P.I Kattan著，韓來彬譯.MATLAB有限元分析與應用[M].北京：清華大學出版社，2004.

[2] 龍奴球，包世華. 結(jié)構(gòu)力學I（第二版）[M].北京：高等教育出版社，2006.