MATLAB在結(jié)構(gòu)力學(xué)分析中的應(yīng)用

吳小希 周小波 李彬

摘 要:傳統(tǒng)的手算方法解超靜定結(jié)構(gòu)工作量繁重，有時(shí)甚至是不可能，運(yùn)用結(jié)構(gòu)有限元編程的一般方法，通過(guò)兩個(gè)實(shí)例的對(duì)照，展示MATLAB在結(jié)構(gòu)力學(xué)分析中的應(yīng)用，MATLAB具有高性能，方法具有普遍的適用性,實(shí)現(xiàn)彎矩圖自動(dòng)繪制。

關(guān)鍵詞: MATLAB結(jié)構(gòu)有限元彎矩圖

Abstract：While using the traditional manual method to resolve complex statically indeterminate structures, it is heavy workloads, sometimes even impossible，using finite element programming of the general method, Based on two examples, This paper introduces a method of application of MATLAB in structure mechanics, MATLAB has the advantages of high performance, it can be applied to many kinds of structures, realization of automatic drawing bending moment diagram.

Key words: MATLAB; Finite element; Bend moment diagram

引言

結(jié)構(gòu)力學(xué)[3]中，常利用傳統(tǒng)的力法與位移法求解超靜定結(jié)構(gòu)，力法是幾何問(wèn)題，位移法把復(fù)雜的幾何圖乘轉(zhuǎn)化為代數(shù)運(yùn)算，但它們基本未知量很多時(shí)，系數(shù)構(gòu)成的矩陣計(jì)算巨大，兩者都不能滿(mǎn)足科研工作者的需要。應(yīng)用MATLAB軟件豐富可靠的矩陣運(yùn)算、數(shù)據(jù)處理、圖形繪制等便利工具，可使得計(jì)算和圖象一體化。對(duì)于結(jié)構(gòu)力學(xué)計(jì)算是十分有利的工具。

1基本方法

MATLAB結(jié)構(gòu)有限元編程的基本思路是先分后合，即將結(jié)構(gòu)分成各個(gè)單元和節(jié)點(diǎn)，桁架與剛架已經(jīng)離散化，對(duì)于連續(xù)系統(tǒng)這一步極其重要，然后進(jìn)行單元分析，集成整體剛度矩陣，引入邊界條件，最后解方程。在求解平面桁架結(jié)構(gòu)，雖然結(jié)構(gòu)簡(jiǎn)單，用手算可得各桿件的軸力，但重復(fù)的過(guò)程太多，現(xiàn)在使用MATLAB語(yǔ)言來(lái)編制有限元位移法的程序時(shí)，則編程的難度明顯降低，對(duì)有限元位移法的概念的理解更加深入，編程所需時(shí)間也大大減少。

圖1為一平面桁架，各桿E=70GPaA=0.004，試用矩陣位移法求解各桿軸力

圖1

解：平面桁架元是既有局部坐標(biāo)又有總體坐標(biāo)的二維有限元；對(duì)各結(jié)點(diǎn)和單元進(jìn)行編號(hào)，建立結(jié)構(gòu)坐標(biāo)系( 圖1 )

第一步,利用MATLAB函數(shù)

y=Plane Truss Element Length(x1, y1, x2, y2)

L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); ％ 局部坐標(biāo)中桿件長(zhǎng)度

第二步, MATLAB函y=Plane Truss Element Stiffness(E ,A ,L ,theta)

x=theta*pi/180; C= cos (x); S=sin(x);

y=E*A/L*[C*C C*S -C*C -C*S; C*S S*S -C*S -S*S;-C*C -C*S C*C C*S;-C*S -S*S C*S S*S];% 總體坐標(biāo)中建立各單元的剛度 矩陣

第三步,建立整體剛度陣。該結(jié)構(gòu)有4個(gè)節(jié)點(diǎn),每個(gè)節(jié)點(diǎn)有兩個(gè)自由度(可考慮支座沉降),為了得到整體剛度陣K，首先利用生成一個(gè)8×8的0矩陣，因?yàn)樵摻Y(jié)構(gòu)有4個(gè)單元，所以4次調(diào)用M a t lab的Plane Truss Assemble函數(shù);其中K為整體剛度陣, k為單元?jiǎng)偠汝? i j為單元兩端在整體節(jié)點(diǎn)上的編號(hào)。

y=Plane Truss Assemble (K, k, i , j)

K (2*i-1, 2*i-1) =K (2*i-1, 2*i-1) +k (1, 1);

K (2*i-1, 2*i) =K (2*i-1, 2*i) + k (1, 2);

K (2*i-1, 2*j-1) = K (2*i-1, 2*j-1)+ k (1,3);

K (2*i-1, 2*j) =K (2*i-1, 2*j) +k (1, 4);

K (2* i , 2*i-1) =K (2* i, 2*i-1) +k (2, 1);

K (2*i, 2*i) =K (2*i, 2*i) +k (2, 2);

K (2*i , 2*j-1)=K(2*i,2*j-1)+k(2,3);

K (2*i, 2*j) =K (2*i, 2*j) + k (2, 4);

K (2*j-1, 2*i-1) =K (2*j, 2*i-1) +k (3, 1);

K (2*j-1, 2*i) =K (2*j-1, 2*i) + k (3, 2);

K (2*j-1, 2*j-1) =K (2*j-1, 2*j-1) +k (3, 3);

K (2*j-1, 2*j) =K (2*j-1, 2*j) + k (3, 4);

K (2*j, 2*i-1) =K (2*j, 2*i-1) + k (4, 1);

K (2*j, 2*i) =K (2*j, 2*i) +k (4, 2);

K (2*j, 2 *j-1) =K (2*j,2*j-1)+k (4,3);

K (2*j, 2*j) = K (2*j, 2*j) +k (4, 4);

y=K;

第四步k=K(3:6,3:6);％邊界條件下剛度矩陣

f=[0;30;30;0];％形成荷載向量

u=kf;％分解法和高斯消去法,得到結(jié)點(diǎn)位移

u = [0.00100.0006 ；0.0011-0.0003]

％結(jié)點(diǎn)2、3的結(jié)點(diǎn)位移

U=[0;0;u;0;0]; ％結(jié)構(gòu)各節(jié)點(diǎn)位移矢量

第五步，M a t lab函數(shù)

Plane Truss Element Force (E, A ,L ,theta ,u)

x=theta*pi/180;C=cos(x);S=sin(x); y=E*A/L*[-C -S C S]*u;

可得：F1 =39.8018F2=0F3 = 28.5646 F4 = -13.8618 F5 = 9.801F6 = -20.1982

%各桿件的軸力

圖2E=210GPa, I=5*10-6 q=7KN/M，繪制彎矩圖。

圖2

解：對(duì)連續(xù)結(jié)構(gòu)單元進(jìn)行編號(hào)十分重要，梁?jiǎn)卧羌扔芯植孔鴺?biāo)又有總體坐標(biāo)的二維有限元，用線(xiàn)性函數(shù)表示，主程序根據(jù)交互輸入的原始數(shù)據(jù)形成單元?jiǎng)偠染仃?，再根?jù)整體剛度矩陣集成規(guī)則，將單元?jiǎng)偠染仃囆纬烧w剛度矩陣。通過(guò)引人支承條件，然后分解和高斯消去法解方程，得到結(jié)點(diǎn)位移，進(jìn)而求出各單元桿端彎矩。

第一步，MATLAB函數(shù)k=Beam Element Stiffness(E,I,L)

y=E*I/(L*L*L)*[12 6*L -12 6*L; 6*L 4*L*L -6*L 2*L*L;-12 -6*L 12 -6*L; 6*L 2*L*L -6*L 4*L*L]; %單元?jiǎng)偠?/p>

第二步，整體剛度的建立，兩者都是二維有限元，程序相同，根據(jù)劃分單元數(shù)，多次調(diào)用函數(shù)。

第三步，計(jì)算等效節(jié)點(diǎn)載荷。按照結(jié)構(gòu)力學(xué)的方法可以求得

M = [-9.333 9.333]

第四步，引入邊界條件，節(jié)點(diǎn)2，3 的轉(zhuǎn)角為a1 a2；其余為0；得出邊界條件下結(jié)構(gòu)剛度矩陣k

第五步，MATLAB函數(shù)

f=Beam Element Forces（k ,u）求得桿端彎矩

M1 (3.111,-6.222)M2= (-6.222, -6.222) M3= (-6.222,3.111);%桿端彎矩

繪制彎矩圖：在命令窗口輸入如下命令：

q=7

l1=linspace (0, 4);

M1=linspace (3.111,-6.2 22);

l= [4, 6, 8];

mid=q*4^2/8-(6.222+6.222)/2

M= [ -6.222 , mid,-6.222]

P= polyfit (l, M, 2)

l2=4:0.1:8

M2=P (1)*l2.^2+P(2)*l2+P(3);

l3=linspace (8, 12);

M3=linspace (-6.222, 3.111);

plot (l1,M1,l2,M2,l3,M3) 自動(dòng)生成彎矩圖3:

2．結(jié)束語(yǔ)

通過(guò)兩個(gè)例子表明, 在結(jié)構(gòu)力學(xué)中引入MATLAB，簡(jiǎn)單的分步編程，即可完成有關(guān)問(wèn)題的過(guò)程分析、大量計(jì)算和繪圖，平面桁架與連續(xù)梁?jiǎn)卧{(diào)用同一函數(shù)，即可求出整體剛度矩陣，大大提高了效率。更多地了解和掌握MATLAB,對(duì)于我們的教學(xué)和科研工作將是十分有益的。

參考文獻(xiàn)：

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

[2]馬曉光,于國(guó)清．MATLAB在結(jié)構(gòu)力學(xué)中的應(yīng)用．白城師范學(xué)院學(xué)報(bào)[J]2006，

20（4）99-102．

[3]龍馭球，包世華．結(jié)構(gòu)力學(xué)I (第二版)［M］．北京：高等教育出版社，2006．

作者簡(jiǎn)介：

吳小希(1987-)，男（漢族），湖南新化,湖南科技大學(xué)研究生，主要研究領(lǐng)域?yàn)闃蛄旱恼駝?dòng)控制。單位：湖南科技大學(xué)土木工程學(xué)院橋梁研究所

注：文章內(nèi)所有公式及圖表請(qǐng)用PDF形式查看。