基于-Matlab-手動編程了解有限元的基本原理

1、基于 Matlab 手動編程了解有限元的基本原理盧鳴飛 復(fù)旦大學(xué)力學(xué)與工程科學(xué)系1背景闡述 在有限元方法飛速發(fā)展的今天，各類有限元軟件也越來越普及，功能越來越強大。并且，針對各個不同行業(yè) 的專業(yè)有限元軟件也猶如雨后春筍般地涌現(xiàn)出來，大大降低了有限元軟件的使用門檻。但其背后的理論基礎(chǔ) 卻是大同小異。 本文試圖通過梳理有限元最基本的原理，讓大家了解如何通過數(shù)值計算來求解一個簡單的實際問題。通過原 理來了解數(shù)值計算本質(zhì)，能夠更有利于大家理解、掌握有限元技術(shù)。2問題描述一電線塔如圖 1 所示，所有桿件間的連接均為鉸接。已知彈性模量 2.0x1011 Pa，密度為 7.8x103 kg/m3。桿 件的截

2、面積 A=2.827433 x10-3 m2。求解：系統(tǒng)的前十階固有頻率。F (t) 1 1 1 1 1 1 1 11111輸出點0.080.11圖 1 電線塔模型3理論分析 由于是桿件系統(tǒng)，所以只需要用到桿件系統(tǒng)的質(zhì)量矩陣與剛度矩陣即可：c2EA cs-c2s2-cs-cs-s2 2010r Al 0201k (e) =m(e) =l 對稱c2cs 6 1020s2 0102在固有頻率計算中，采用響應(yīng)分析的 Houbolt 方法，先計算 x1，然后計算 x-1，再計算 x2，x3，其中 x0 是由初始條件給出的。這樣就可以得到全部的位移矩陣了。由如下框圖可以明確表示 Houbolt 方法的步

3、進過 程：計算得出 x1計算得出 x-1計算得出 x2，x3.f1= f - m x ,1m =m +c,47Dt 26Dtf = f + (2m + Dt c)x + ( 6 m + 2c)x120Dt0+(m +c)x63Dt 2Dt0k =m +c + k211f1Dt 26Dt= (k + m )x1kxi+1 = fi+1,i = 1, 2, n -1f= f+ ( 5 m + 3 c)x - ( 4 m +c)x3i+1i+1Dt 2DtiDt 22Dti-1+ ( 1 m +c)x1Dt 23Dti-2x= Dt2 x + 2x - x-1001然后根據(jù)得到的位移矩陣，分別計算速

4、度矩陣和加速度矩陣：x1i+1=(11x-18x + 9x- 2x),6Dti+1ii-1i-2x1i+1=(2x- 5x + 4x- x).Dt 2i+1ii-1i-24源程序程序見最后。在 Matlab 中運行 Main.m 文件即可計算得到結(jié)果。5計算結(jié)果系統(tǒng)的前十階固有頻率：NOFREQUENCY124.0320281.21983112.62794204.89075225.81946344.76107382.97268441.95049450.569010533.2769Main.mclc clearCOOR = -1 0 1 2 0 1 0 1 0 1 0 1 0 1-1 -1 -2

5、2 2 30 0 0 0 1 1 2 2 3 3 4 4 5 5 4.544 4.5 4 4;NOD = 1 2 3 2 3 4 5 5 5 6 6 7 7 878 99 9 10 10 11 11 11 12 12 13 11 11 1315 15 16 12 12 14 18 18 195 5 5 6 6 6 6 7 8 7 8 8 9 9 10 10 10 11 12 11 12 12 13 14 13 14 14 15 16 1516 17 17 18 19 18 19 20 20;BOU = 1 1 2 2 3 3 4 41 2 1 2 1 2 1 2;PRO = 2.0e11 2.

6、898119e-67.8e32.827433e-3; t = 2;delta = 0.001;c = 0.01;drawplot(COOR,NOD);ID,KA MA = Matrix(COOR,NOD,BOU,PRO); V,D = eig(KA,MA);Output(10,V,D);F = loadf(7,1,delta,ID,t);x,velx,accx = Houbolt(KA,MA,V,c,F,delta); Pout(t,delta,x,velx,accx,ID,20);drawplot.mfunction drawplot(COOR,NOD) n = size(NOD,2);fi

7、gure for i=1:nplot(COOR(1,NOD(1:2,i),COOR(2,NOD(1:2,i),o-,LineWidth,2);hold on; axis square;axis(equal);endMatrix.mfunction ID,KA,MA = Matrix(COOR,NOD,BOU,PRO) NUMNP = size(COOR,2);NUME = size(NOD,2); ID = zeros(2,NUMNP);for n = 1:size(BOU,2) ID(BOU(2,n),BOU(1,n) = 1;endNEQ = 0;for n = 1:NUMNP for i

8、 = 1:2if ID(i,n) 0 for j = 1:4if LM(j) 0endendendendendKA(LM(i),LM(j) = KA(LM(i),LM(j)+Ke(i,j);MA(LM(i),LM(j) = MA(LM(i),LM(j)+Me(i,j);Output.mfunction Output(m,V,D) n = size(V,2);if m nfprintf(Error! The maximum is %d.nPlease check Output(m,KA,MA) to make sure m%dn,n,n);returnendomega = sqrt(diag(D

9、)/2/pi; fprintf(NOFREQUENCYn); for i = 1:mfprintf(%5d%11.4fn,i,omega(i);endfprintf(nMain vibration mode of the top %d :n,m); for j = 1:mfprintf(%9d,j);end fprintf(n); for i = 1:nfor j = 1:mfprintf(%9.4f,V(i,j);endend fprintf(n);loadf.mfunction F = loadf(i,j,delta,ID,t) m = max(max(ID);n = t/delta;F

10、= zeros(m,n+1); m = ID(j,i);F(:,1) = 0;for k = 1:0.2/delta F(m,k+1) = 5000*k*delta;endfor k = 0.2/delta+1:0.3/delta; F(m,k+1) = -20000*k*delta+5000;endfor k = 0.3/delta+1:0.5/delta F(m,k+1) = 10000*k*delta-4000;endfor k = 0.5/delta+1:0.6/delta F(m,k+1) = -10000*k*delta+6000;endfor k = 0.6/delta+1:n

11、F(m,k+1) = 0;endHoubolt.mfunction x,velx,accx = Houbolt(KA,MA,V,c,F,delta) K = V*KA*V;M = V*MA*V; F = V*F;n = size(F,1);m = size(F,2); temp = zeros(n,2); x = zeros(n,m); velx = zeros(n,m); accx = zeros(n,m); x(:,1) = 0;velx(:,1) = 0; for i = 1:naccx(i,1) = (F(i,1)-c*velx(i,1)-K(i,i)*x(i,1)/M(i,i);en

12、dfor i = 1:nx(i,2) =(F(i,2)+(2*M(i,i)+delta/2*c)*accx(i,1)+(6/delta*M(i,i)+2*c)*velx(i,1)+.(6/delta2*M(i,i)+3/delta*c)*x(i,1)/(2/delta2*M(i,i)+11/6/delta*c+K(i,i)+. (4/delta2*M(i,i)+7/6/delta*c);endfor i = 1:ntemp(i,2) = delta2*accx(i,1)+2*x(i,1)-x(i,2);temp(i,1) = 6*delta2*accx(i,1)+6*delta*velx(i,

13、1)+9*x(i,1)-8*x(i,2);endfor i = 1:nx(i,3) = (F(i,3)+(5/delta2*M(i,i)+3/delta*c)*x(i,2)- (4/delta2*M(i,i)+3/2/delta.*c)*x(i,1)+(M(i,i)/delta2+c/3/delta)*temp(i,1)/(2/delta2*M(i,i)+11/6/delta.*c+K(i,i);endfor j = 4:mfor i =1:nx(i,j) = (F(i,j)+(5/delta2*M(i,i)+3/delta*c)*x(i,j-1)- (4/delta2*M(i,i)+3/2/

14、.delta*c)*x(i,j-2)+(M(i,i)/delta2+c/3/delta)*x(i,j- 3)/(2/delta2*M(i,i).+11/6/delta*c+K(i,i);endendx = V*x; for j = 1:m-1for i =1:nif j=1a = temp(i,2);b = temp(i,1); else if j=2a = x(i,1);b = temp(i,2); elsea = x(i,j-1);b = x(i,j-2);endendendendvelx(i,j+1) = (11*x(i,j+1)-18*x(i,j)+9*a-2*b)/6/delta; accx(i,j+1) = (2*x(i,j+1)-5*x(i,j)+4*a-b);velx = V*velx; accx = V*accx; endPout.mfunction Pout(t,delta,x,velx,accx,ID,i) m=(x(ID(1,i),:).2+x(ID(2,i),:).2).0.5;n=velx(ID(1,i),:);p=accx