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

黃 輝

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

0 引言

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

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

1 MATLAB有限元分析

1.1 運(yùn)用MATLAB進(jìn)行有限元方法解決實(shí)際問題的步驟[2]

運(yùn)用MATLAB進(jìn)行有限元方法解決實(shí)際問題的6個(gè)步驟簡述如下:

1)離散化域:將一個(gè)表示結(jié)構(gòu)或連續(xù)體的求解域離散為若干個(gè)不同的有限大小和形狀的子域(單元),子域和子域之間通過節(jié)點(diǎn)進(jìn)行連接(手動(dòng)進(jìn)行離散);

2)寫出單元?jiǎng)偠染仃?用MATLAB寫出域內(nèi)每個(gè)單元的單元?jiǎng)偠染仃?

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

4)引入邊界條件:引入支座、外加荷載和位移等邊界條件(手動(dòng)進(jìn)行實(shí)現(xiàn));

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

6)后處理:根據(jù)求得的結(jié)構(gòu)節(jié)點(diǎn)位移矢量,求解支反力、單元節(jié)點(diǎn)力和單元應(yīng)力。

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

借助于MATLAB強(qiáng)大的編程功能,寫出有關(guān)平面梁元結(jié)構(gòu)的5個(gè)M函數(shù)文件:

1)BeamElementStiffness(E,I,L,)——該函數(shù)用于計(jì)算彈性模量E、轉(zhuǎn)動(dòng)慣量I、長度L的梁元的單元?jiǎng)偠染仃?。它返?×4的單元?jiǎng)偠染仃噆。

2)BeamElementAssemble(K,k,i,j)——該函數(shù)將連接節(jié)點(diǎn)i和節(jié)點(diǎn)j的梁元的單元?jiǎng)偠染仃噆集成到整體剛度矩陣K。每集成一個(gè)單元,該函數(shù)都將返回2n×2n的整體剛度矩陣K。

3)BeamElementForces(k,u)——該函數(shù)用單元?jiǎng)偠染仃噆和單元節(jié)點(diǎn)位移矢量u計(jì)算單元節(jié)點(diǎn)力矢量。它返回4×1的單元節(jié)點(diǎn)力矢量f。

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

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

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

2 實(shí)例分析

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

圖1 三單元四節(jié)點(diǎn)連續(xù)梁

2.1 模型簡化及計(jì)算

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

2.2 結(jié)果后處理

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

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

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

2.3 MATLAB實(shí)現(xiàn)過程(程序代碼)

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 結(jié)語

本文以平面梁元為例,闡述了基于MATLAB有限元法結(jié)構(gòu)分析。經(jīng)實(shí)踐證明,該方法簡便、可行、實(shí)用。雖然文中以平面梁元為例介紹其應(yīng)用,但其中的原理和方法可推廣到其它單元?jiǎng)偠染仃?甚至三維空間和非線性剛度矩陣之中,對此,將另文做進(jìn)一步探討。

[1]徐榮橋.結(jié)構(gòu)分析的有限元法與MATLAB程序設(shè)計(jì)[M].北京:人民交通出版社,2006.

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

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

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

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

[6]周水興,陳山林.MATLAB在有限元?jiǎng)偠染仃囃茖?dǎo)中的應(yīng)用[J].重慶交通學(xué)院學(xué)報(bào),2007,26(2):29-31.

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

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