MATLAB在平面桁架計算中的應用

2017-01-21 16:21:07王鵬周琪琛姚姍姍

卷宗 2016年10期

王鵬+周琪琛+姚姍姍

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

關鍵詞：平面桁架；MATLAB；有限元分析

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

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

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

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

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

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

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

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

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

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

集成整體剛度矩陣：根據上一步得到的單元剛度矩陣，通過剛度集成規則，集成結構整體的剛度矩陣，這一步也可以直接用MATLAB工具箱來完成，

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

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

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

3 實例分析

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

如圖為一平面桁架結構，為了便于計算，

假定結構的彈性模E=200GPa

橫截面積A=2cm2，結構受力如圖所示：

先對結構進行離散化處理，由于桁架結構已經是離散化結構，所以我們只需要將結構單元編號即可：

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

% 輸入參數 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）]；

% 計算節點力矢量 F=K*U；

% 調用節點位移矢量

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）；

運行后的部分結果：