基于MATLAB有限元法結構分析

黃 輝

(1.湖南省永州市交通局,湖南永州 425006; 2.中南大學土木建筑學院,湖南長沙 410075)

0 引言

目前在工程技術領域內常用的數值模擬方法包括有限單元法、邊界單元法、離散單元法和有限元差分法等。但就實用性和應用的廣泛性而言,主要還是有限單元法。MATLAB是當今國際科技領域內具有較大影響力及活力的軟件。它起源于矩陣運算,并已經發展成一種高度集成的計算機語言。它提供了強大的科學運算、靈活的程序設計流程、高質量的圖形可視化與界面設計、便捷的與其他程序和語言接口的功能。工程人員通過自己編寫的M函數文件或調用MATLAB提供的工具箱,可以高效求解較復雜的工程問題,并對系統進行動態仿真和數值計算結果的圖形輸出。

在一般的結構問題分析中,平面梁元(Beam Element)是使用較多的二維有限元之一。本文以平面梁元為例探討有限元法在MATLAB中的應用。平面梁元的系數有彈性模量E、慣性矩I和長度L。每個梁元有2個節點,并且假定它是水平的。梁元有4個自由度——每個節點有2個自由度(橫位移和轉角)。約定位移向上為正,轉角逆時針方向為正。對一個有n個節點的結構而言,其整體剛度矩陣K的維數為2 n×2 n。忽略軸向形變,單元剛度矩陣k給定如式(1)[1-3]。

1 MATLAB有限元分析

1.1 運用MATLAB進行有限元方法解決實際問題的步驟[2]

運用MATLAB進行有限元方法解決實際問題的6個步驟簡述如下:

1)離散化域:將一個表示結構或連續體的求解域離散為若干個不同的有限大小和形狀的子域(單元),子域和子域之間通過節點進行連接(手動進行離散);

2)寫出單元剛度矩陣:用MATLAB寫出域內每個單元的單元剛度矩陣;

3)集成整體剛度矩陣:用MATLAB使用直接剛度法集成整體剛度矩陣;

4)引入邊界條件:引入支座、外加荷載和位移等邊界條件(手動進行實現);

5)求解方程:分解(手動進行)整體剛度矩陣并用高斯消去法(MATLAB)求解方程組;

6)后處理:根據求得的結構節點位移矢量,求解支反力、單元節點力和單元應力。

1.2 編寫M函數文件

借助于MATLAB強大的編程功能,寫出有關平面梁元結構的5個M函數文件:

1)BeamElementStiffness(E,I,L,)——該函數用于計算彈性模量E、轉動慣量I、長度L的梁元的單元剛度矩陣。它返回4×4的單元剛度矩陣k。

2)BeamElementAssemble(K,k,i,j)——該函數將連接節點i和節點j的梁元的單元剛度矩陣k集成到整體剛度矩陣K。每集成一個單元,該函數都將返回2n×2n的整體剛度矩陣K。

3)BeamElementForces(k,u)——該函數用單元剛度矩陣k和單元節點位移矢量u計算單元節點力矢量。它返回4×1的單元節點力矢量f。

4)BeamElementShearDiagram(f,L)——該函數繪制節點力矢量為f和長度為L的單元剪力圖。

5)BeamElementMomentDiagram(f,L)——該函數繪制節點力矢量為f和長度為L的單元彎矩曲線圖。

在MATLAB軟件集成環境中,使用上述5個M函數文件時,首先根據實際結構給出函數所需要用到的參數(E,I,L,K,k,i,j),然后可任意調用5個M函數文件中的一個進行計算。因本文重點不在講述MATLAB軟件的使用方法,如何在MATLAB集成環境中編寫并使用M函數文件請參閱文獻[4]。

2 實例分析

下面通過一個實例驗證MATLAB的有限元求解全過程。如圖1所示的連續梁[5],已知E=210GPa,I=50×10-6m4。求節點2、節點3和節點4的轉角,并繪出單元的剪力圖和彎矩圖(計算單位取kN和m)。

圖1 三單元四節點連續梁

2.1 模型簡化及計算

將結構劃分為3個單元4個節點,單元①由節點1和節點2組成;單元②由節點2和節點3組成;單元③由節點3和節點4組成。首先,需要用等效節點荷載分別替代單元①和③的分布荷載和非節點集中荷載。然后根據M函數求出單元的剛度矩陣k1,k2,k3;該結構有4個節點,所以整體剛度矩陣是8×8的;由于有3個單元,故3次調用BeamElementAssemble()函數,求得整體剛度矩陣K。引入本題的邊界條件后,分解(手動進行)整體剛度矩陣K并用高斯消去法(MATLAB)求解,最后解得節點2、3、4的轉角分別是0.000 6、-0.001 2、0.002 0 rad(正號表示逆時針方向;負號表示順時針方向)。

2.2 結果后處理

首先用MATLAB命令求出節點1、2、3、4的支反力和內力(剪力和彎矩)。建立結構的節點位移矢量U,然后計算結構的節點荷載矢量F。故節點1、2、3、4的垂直支反力分別是3.994 6、-8.478 3、7.724 7、-3.240 5 kN(正號表示方向向上;負號表示方向向下)。節點1、2、3、4的彎矩分別是3.994 6、7.5、-15、15 kN?m(正號表示逆時針方向;負號表示順時針方向)。很顯然,滿足力的平衡。

然后建立單元位移矢量u1、u2和u3,然后用BeamElementForces()函數求出單元矢量f1、f2和f3。需要注意的是,由于分布荷載和非節點集中荷載,所以需要對單元①和③的力進行修正。例如,修正后單元矢量f1為[18.994 6;11.494 6;11.005 4;0.489 1]。

最后調用BeamElementShearDicgram()和BeamElementMomentDicgram()函數分別繪出單元的剪力圖和彎矩圖。本文僅給出單元①的剪力圖(圖2)和彎矩圖(圖3)。

2.3 MATLAB實現過程(程序代碼)

clear

E=210e6;

I=50e-6;

L1=3;

L2=3;

L3=4;

k1=BeamElementStiffness(E,I,L1);

k2=BeamElementStiffness(E,I,L2);

k3=BeamElementStiffness(E,I,L3);

K=zeros(8,8);

K=BeamAssemble(K,k1,1,2);

K=BeamAssemble(K,k2,2,3);

K=BeamAssemble(K,k3,3,4)

k=[K(4,4)K(4,6)K(4,8);K(6,4)K(6,6)K(6,8);K(8,4)K(8,6)K(8,8)]

f=[7.5;-15;15]

u=kf

U=[0;0;0;u(1);0;u(2);0;u(3)]

F=K*U

u1=[U(1);U(2);U(3);U(4)];

u2=[U(3);U(4);U(5);U(6)];

u3=[U(5);U(6);U(7);U(8)];

f1=BeamElementForces(k1,u1);

f1=f1-[-15;-7.5;-15;7.5]

圖2 單元①的剪力圖(kN)

圖3 單元①的彎矩圖(kN?m)

f2=BeamElementForces(k2,u2)

f3=BeamElementForces(k3,u3);

f3=f3-[-15;-15;-15;15]

BeamElementShearDiagram(f1,L1)

BeamElementShearDiagram(f2,L2)

BeamElementShearDiagram(f3,L3)

BeamElementMomentDiagram(f1,L1)

BeamElementMomentDiagram(f2,L2)

BeamElementMomentDiagram(f3,L3)

3 結語

本文以平面梁元為例,闡述了基于MATLAB有限元法結構分析。經實踐證明,該方法簡便、可行、實用。雖然文中以平面梁元為例介紹其應用,但其中的原理和方法可推廣到其它單元剛度矩陣,甚至三維空間和非線性剛度矩陣之中,對此,將另文做進一步探討。

[1]徐榮橋.結構分析的有限元法與MATLAB程序設計[M].北京:人民交通出版社,2006.

[2](德)P.I.Kattan,韓來彬譯.MATLAB有限元分析與應用[M].北京:清華大學出版社,2004.

[3](美)Chandrupatla,T,R.,Belegundu,A.D.,曾攀譯.工程中的有限元方法(第3版)[M].北京:清華大學出版社,2006.

[4]隋思漣,王 巖.MATLAB語言與工程數據分析[M].北京:清華大學出版社,2009.

[5]包世華.結構力學(第3版,下冊)[M].武漢:武漢理工大學出版社,2008.

[6]周水興,陳山林.MATLAB在有限元剛度矩陣推導中的應用[J].重慶交通學院學報,2007,26(2):29-31.

[7]景振毅,張澤兵,董 霖.MATLAB 7.0實用寶典[M].北京:中國鐵道出版社,2009.

[8]丁 星.應用矩陣變換繪制平面剛架內力圖和位移圖[J].四川建筑科學研究,2010(3).