MATLAB在結構力學分析中的應用

吳小希 周小波 李彬

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

關鍵詞: MATLAB結構有限元彎矩圖

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

引言

結構力學[3]中，常利用傳統的力法與位移法求解超靜定結構，力法是幾何問題，位移法把復雜的幾何圖乘轉化為代數運算，但它們基本未知量很多時，系數構成的矩陣計算巨大，兩者都不能滿足科研工作者的需要。應用MATLAB軟件豐富可靠的矩陣運算、數據處理、圖形繪制等便利工具，可使得計算和圖象一體化。對于結構力學計算是十分有利的工具。

1基本方法

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

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

圖1

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

第一步,利用MATLAB函數

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

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

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;％分解法和高斯消去法,得到結點位移

u = [0.00100.0006 ；0.0011-0.0003]

％結點2、3的結點位移

U=[0;0;u;0;0]; ％結構各節點位移矢量

第五步，M a t lab函數

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

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

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

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

第三步，計算等效節點載荷。按照結構力學的方法可以求得

M = [-9.333 9.333]

第四步，引入邊界條件，節點2，3 的轉角為a1 a2；其余為0；得出邊界條件下結構剛度矩陣k

第五步，MATLAB函數

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．結束語

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

參考文獻：

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

[2]馬曉光,于國清．MATLAB在結構力學中的應用．白城師范學院學報[J]2006，

20（4）99-102．

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

作者簡介：

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

注：文章內所有公式及圖表請用PDF形式查看。