四邊形八節(jié)點等參元matlab程序剖析

1、四邊形八節(jié)點等參元matlab程序 懸臂鋼梁，尺寸如圖一所示；v=0.3。h=1，E=2.1e11. 圖一 懸臂鋼梁 圖二 單元劃分與結(jié)點編號 1 四邊形八節(jié)點等參元matlab程序 Matlab 輸出結(jié)果 2 四邊形八節(jié)點等參元matlab程序 附錄： 有限元ANSYS分析結(jié)果 采用PLANE183單元（四邊形八節(jié)點）單元得出的結(jié)構(gòu)Y向最大位移為-0.216E-04。約等于matlab平面四邊形八節(jié)點等參元結(jié)點Y向最大位移-2.4024E-5。 3 四邊形八節(jié)點等參元matlab程序 附錄： %-四邊形八節(jié)點等參元 matlab計算程序- % 主 程 序 %*%* % 2012年 % 本程序

2、只能處理集中荷載作用下的情況 % 只輸出了節(jié)點位移、單元中心點的應(yīng)力 %*%* % 變量說明 % E v h % 彈性模量 泊松比 厚度 % NPOIN NELEM NVFIX NNODE NFPOIN % 總結(jié)點數(shù) , 單元數(shù), 約束結(jié)點個數(shù), 單元節(jié)點數(shù) ,受力結(jié)點數(shù) % COORD LNODS % 結(jié)構(gòu)節(jié)點整體坐標數(shù)組, 單元定義數(shù)組, % FPOIN FORCE FIXED % 結(jié)點力數(shù)組， 總體荷載向量, 約束信息數(shù)組 % HK DISP % 總體剛度矩陣,結(jié)點位移向量 %* clear all format short e FP1=fopen(bjd.txt,rt); %打開數(shù)據(jù)文

3、件 %讀入控制數(shù)據(jù) E=fscanf(FP1,%f,1); %彈性模量 v=fscanf(FP1,%f,1); % 泊松比 h=fscanf(FP1,%f,1); %厚度 NELEM=fscanf(FP1,%d,1); %單元數(shù) NPOIN=fscanf(FP1,%d,1); % 總結(jié)點數(shù) NNODE=fscanf(FP1,%d,1); %單元節(jié)點數(shù) NFPOIN=fscanf(FP1,%d,1); %受力結(jié)點數(shù) NVFIX=fscanf(FP1,%d,1); %約束結(jié)點個數(shù) LNODS=fscanf(FP1,%f,NNODE,NELEM); % 單元定義： 單元結(jié)點號（逆時針） COORD=

4、fscanf(FP1,%f,2,NPOIN); % 結(jié)點號 x,y坐標(整體坐標下) FPOIN=fscanf(FP1,%f,3,NFPOIN); 4 四邊形八節(jié)點等參元matlab程序 % 節(jié)點力：結(jié)點號、X方向力(向右正)，Y方向力(向上正) FIXED=fscanf(FP1,%d,3,NVFIX); %約束信息數(shù)組(n,3) n:受約束節(jié)點數(shù)目, (n,1):約束點號 %(n,2)與(n,3)分別為約束點x方向和y方向的約束情況,受約束為1否則為0 %* %* %=平面應(yīng)力問題的求解= % %* %* % % 剛度矩陣的生成 %計算剛度矩陣，并對約束條件進行處理 Ke=zeros(2*N

5、NODE,2*NNODE); % 單元剛度矩陣并清零 HK=zeros(2*NPOIN,2*NPOIN); % 張成總剛矩陣并清零 %調(diào)用子程序 生成單元剛度矩陣 for m=1:NELEM %m為單元號 Ke=K(E,v,h,. COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),. COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),. COORD(LNODS(m,5),1),COORD(LNODS(m,5),2),. COORD(LNODS(m,7),1),COORD(LNODS(m,7),2); %調(diào)用單元剛度矩陣 a=LNODS

6、(m,:); %臨時向量,用來記錄當前單元的節(jié)點編號 %對總剛度矩陣的處理 for j=1:8 for k=1:8 HK(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+. Ke(j*2-1:j*2,k*2-1:k*2); end end end % %對荷載向量進行處理 FORCE=zeros(2*NPOIN,1); % 張成總荷載向量并清零 for i=1:NFPOIN b1=FPOIN(i,1)*2-1;b2=FPOIN(i,1)*2; %FPION(i,1)為作用點 5 四邊形八節(jié)點等

7、參元matlab程序 FORCE(b1)=FPOIN(i,2); %FPION(i,2)為x方向的節(jié)點力 FORCE(b2)=FPOIN(i,3); %FPION(i,3)為y方向的節(jié)點力 end % %將約束信息加入總剛,總荷載 for i=1:NVFIX if FIXED(i,2)=1 c1=2*FIXED(i,1)-1; HK(c1,:)=0; %將一約束序號處的總剛列向量清0 HK(:,c1)=0; %將一約束序號處的總剛行向量清0 HK(c1,c1)=1; %將行列交叉處的元素置為1 FORCE(c1)=0; end if FIXED(i,3)=1 c2=2*FIXED(i,1);

8、HK(c2,:)=0; HK(:,c2)=0; HK(c2,c2)=1; FORCE(c2)=0; end end % %= %= DISP=HKFORCE %計算節(jié)點位移向量 %= %= %求解單元應(yīng)力 stress=zeros(3,NELEM); for m=1:NELEM u(1:16)=0; d=LNODS(m,:); %臨時向量,用來記錄當前單元的節(jié)點編號 for i=1:NNODE u(i*2-1:i*2)=DISP(d(i)*2-1:d(i)*2); %從總位移向量中取出當前單元的節(jié)點位移 end D=(E/(1-v*v)*1 v 0;v 1 0;0 0 (1-v)/2;%彈性矩

9、陣 %形成應(yīng)變矩陣BM 6 四邊形八節(jié)點等參元matlab程序 BM=zeros(3,16); for i=1:NNODE J=Jacobi(COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),. COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),. COORD(LNODS(m,5),1),COORD(LNODS(m,5),2),. COORD(LNODS(m,7),1),COORD(LNODS(m,7),2),0,0); N_s,N_t=DHS(0,0); B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i); B2i=-J

10、(2,1)*N_s(i)+J(1,1)*N_t(i); BM(1:3,2*i-1:2*i)=B1i 0;0 B2i;B2i B1i/det(J); end stressm=D*BM*u; stress(:,m)=stressm; end stress %輸出應(yīng)力 function Ke=K(E,v,h,x1,y1,x3,y3,x5,y5,x7,y7) %=單元剛度矩陣= % E 彈性模量 % v 泊松比 % h 厚度 % x1,y1,x3,y3,x5,y5,x7,y7 為4個角結(jié)點的坐標 %矩陣尺寸為16 x 16 Ke=zeros(16,16); D=(E/(1-v*v)*1 v 0;v 1

11、 0;0 0 (1-v)/2;%彈性矩陣 %高斯積分 采用 3 x 3 個積分點 書74頁 W1=5/9;W2=8/9;W3=5/9; %加權(quán)系數(shù) W=W1 W2 W3; r=15(1/2)/5; x=-r 0 r;%積分點 for i=1:3 for j=1:3 B=eleB(x1,y1,x3,y3,x5,y5,x7,y7,x(i),x(j); J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,x(i),x(j); Ke=Ke+W(i)*W(j)*B*D*B*det(J)*h; end end function B=eleB(x1,y1,x3,y3,x5,y5,x7,y7,s

12、,t) 7 四邊形八節(jié)點等參元matlab程序 %調(diào)用導函數(shù) N_s,N_t=DHS(s,t); %求Jacobi矩陣 J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t); %求應(yīng)變矩陣B B=zeros(3,16); for i=1:8 B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i); B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i); B(1:3,2*i-1:2*i)=B1i 0;0 B2i;B2i B1i; end B=B/det(J); function J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t)

13、%-Jacobi- %單元坐標 %2,4,6,8點的坐標 x2=(x1+x3)/2;y2=(y1+y3)/2; x4=(x3+x5)/2;y4=(y3+y5)/2; x6=(x5+x7)/2;y6=(y5+y7)/2; x8=(x7+x1)/2;y8=(y7+y1)/2; x=x1 x2 x3 x4 x5 x6 x7 x8; y=y1 y2 y3 y4 y5 y6 y7 y8; %調(diào)用形函數(shù)對局部坐標的導數(shù) N_s,N_t=DHS(s,t); %求Jacobi矩陣的行列式的值 x_s=0;y_s=0; x_t=0;y_t=0; for i=1:8 x_s=x_s+N_s(i)*x(i);y_s

14、=y_s+N_s(i)*y(i); x_t=x_t+N_t(i)*x(i);y_t=y_t+N_t(i)*y(i); end J=x_s y_s;x_t y_t; function N=shape(s,t) 8 四邊形八節(jié)點等參元matlab程序 %, N(1) = (1-s)*(1-t)*(-s-t-1)/4; N(3) = (1+s)*(1-t)*(s-t-1)/4; N(5) = (1+s)*(1+t)*(s+t-1)/4; N(7) = (1-s)*(1+t)*(-s+t-1)/4; N(2) = (1-t)*(1+s)*(1-s)/2; N(4) = (1+s)*(1+t)*(1-t

15、)/2; N(6) = (1+t)*(1+s)*(1-s)/2; N(8) = (1-s)*(1+t)*(1-t)/2; function N_s,N_t=DHS(s,t) %形函數(shù)求導 %, N_s(1)=-1/4*(1-t)*(-s-t-1)-1/4*(1-s)*(1-t); N_s(3)=1/4*(1-t)*(s-t-1)+1/4*(1+s)*(1-t); N_s(5)=1/4*(1+t)*(s+t-1)+1/4*(1+s)*(1+t); N_s(7)=-1/4*(1+t)*(-s+t-1)-1/4*(1-s)*(1+t); N_s(2)=1/2*(1-s)*(1-t)-1/2*(1+s

16、)*(1-t); N_s(4)=1/2*(1+t)*(1-t); N_s(6)=1/2*(1-s)*(1+t)-1/2*(1+s)*(1+t); N_s(8)=-1/2*(1+t)*(1-t); N_t(1)=-1/4*(1-s)*(-s-t-1)-1/4*(1-s)*(1-t); N_t(3)=-1/4*(1+s)*(s-t-1)-1/4*(1+s)*(1-t); N_t(5)=1/4*(1+s)*(s+t-1)+1/4*(1+s)*(1+t); N_t(7)=1/4*(1-s)*(-s+t-1)+1/4*(1-s)*(1+t); N_t(2)=-1/2*(1+s)*(1-s); N_t(4)=1/2*(1+s)*(1-t)-1/2*(1+s)*(1+t); N_t(6)=1/2*(1+s)*(1-s); N_t(8)=1/2*(1-s)*(1-t)-1/2*(1-s)*(1+t); bjd.txt 文件數(shù)據(jù) 9 四邊形八節(jié)點等參元matlab程序 2.1E11 0.3 1 5 28 8 1 3 1 2 3 13 20 19 18 12 3 4 5 14 22 21 20 13 5 6 7 15 24 23 22 14 7 8 9 16 26 25 24 15 9 10 11 17 28 27