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

吳小希 周小波 李彬

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

關(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)，力法是幾何問題，位移法把復(fù)雜的幾何圖乘轉(zhuǎn)化為代數(shù)運算，但它們基本未知量很多時，系數(shù)構(gòu)成的矩陣計算巨大，兩者都不能滿足科研工作者的需要。應(yīng)用MATLAB軟件豐富可靠的矩陣運算、數(shù)據(jù)處理、圖形繪制等便利工具，可使得計算和圖象一體化。對于結(jié)構(gòu)力學(xué)計算是十分有利的工具。

1基本方法

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

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

圖1

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

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

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

L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); ％ 局部坐標中桿件長度

第二步, 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];% 總體坐標中建立各單元的剛度 矩陣

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

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é)點位移

u = [0.00100.0006 ；0.0011-0.0003]

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

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

第五步，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

解：對連續(xù)結(jié)構(gòu)單元進行編號十分重要，梁單元是既有局部坐標又有總體坐標的二維有限元，用線性函數(shù)表示，主程序根據(jù)交互輸入的原始數(shù)據(jù)形成單元剛度矩陣，再根據(jù)整體剛度矩陣集成規(guī)則，將單元剛度矩陣形成整體剛度矩陣。通過引人支承條件，然后分解和高斯消去法解方程，得到結(jié)點位移，進而求出各單元桿端彎矩。

第一步，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]; %單元剛度

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

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

M = [-9.333 9.333]

第四步，引入邊界條件，節(jié)點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) 自動生成彎矩圖3:

2．結(jié)束語

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

參考文獻：

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

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

20（4）99-102．

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

作者簡介：

吳小希(1987-)，男（漢族），湖南新化,湖南科技大學(xué)研究生，主要研究領(lǐng)域為橋梁的振動控制。單位：湖南科技大學(xué)土木工程學(xué)院橋梁研究所

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