黃 輝

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

0 引言

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

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

1 MATLAB有限元分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 實例分析

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

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

2.1 模型簡化及計算

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

2.2 結(jié)果后處理

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

然后建立單元位移矢量u1、u2和u3,然后用BeamElementForces()函數(shù)求出單元矢量f1、f2和f3。需要注意的是,由于分布荷載和非節(jié)點集中荷載,所以需要對單元①和③的力進(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實現(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)實踐證明,該方法簡便、可行、實用。雖然文中以平面梁元為例介紹其應(yīng)用,但其中的原理和方法可推廣到其它單元剛度矩陣,甚至三維空間和非線性剛度矩陣之中,對此,將另文做進(jìn)一步探討。

[1]徐榮橋.結(jié)構(gòu)分析的有限元法與MATLAB程序設(shè)計[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在有限元剛度矩陣推導(dǎo)中的應(yīng)用[J].重慶交通學(xué)院學(xué)報,2007,26(2):29-31.

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

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