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

1、四邊形八節(jié)點等參元matlab程序懸臂鋼梁，尺寸如圖一所示；v=0.3。h=1，E=2.1e11. 圖一 懸臂鋼梁圖二 單元劃分與結(jié)點編號Matlab 輸出結(jié)果 附錄：有限元ANSYS分析結(jié)果采用PLANE183單元（四邊形八節(jié)點）單元得出的結(jié)構(gòu)Y向最大位移為-0.216E-04。約等于matlab平面四邊形八節(jié)點等參元結(jié)點Y向最大位移-2.4024E-5。附錄：%-四邊形八節(jié)點等參元 matlab計算程序-% 主 程 序 %*%*% 2012年 % 本程序只能處理集中荷載作用下的情況% 只輸出了節(jié)點位移、單元中心點的應力%*%*% 變量說明% E v h% 彈性模量 泊松比 厚度% NPOI

2、N 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 allformat short e FP1=fopen('bjd.txt','rt'); %打開數(shù)據(jù)文件%讀入控制數(shù)據(jù)E=fscanf(FP1,'%f',1); %彈性模量 v=fscanf(FP1,'

3、%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)' %

4、 單元定義： 單元結(jié)點號（逆時針）COORD=fscanf(FP1,'%f',2,NPOIN)' % 結(jié)點號 x,y坐標(整體坐標下)FPOIN=fscanf(FP1,'%f',3,NFPOIN)' % 節(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%*%*%=平面應力問題的求解=%*%*% % 剛度矩陣

5、的生成%計算剛度矩陣，并對約束條件進行處理Ke=zeros(2*NNODE,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(L

6、NODS(m,7),2); %調(diào)用單元剛度矩陣a=LNODS(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 endend%對荷載向量進行處理FORCE=zeros(2*NPOIN,1); % 張成總荷載向量并清零for i=1:NFPOIN b1=FPOIN(i,1)*2-1;b2=FPOIN(i,1)*2; %F

7、PION(i,1)為作用點 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); H

8、K(c2,:)=0; HK(:,c2)=0; HK(c2,c2)=1; FORCE(c2)=0; end end%=%=DISP=HKFORCE %計算節(jié)點位移向量%=%=%求解單元應力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;%彈性矩陣 %形成應變矩陣BM

9、 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(2,1)*N_s(i)+J(1,1)*N_t(i); BM(

10、1:3,2*i-1:2*i)=B1i 0;0 B2i;B2i B1i/det(J); end stressm=D*BM*u' stress(:,m)=stressm;endstress %輸出應力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 16Ke=zeros(16,16);D=(E/(1-v*v)*1 v 0;v 1 0;0 0 (1-v)/2;%彈性矩陣%高斯積分 采用 3 x 3

11、個積分點 書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; endendfunction B=eleB(x1,y1,x3,y3,x5,y5,x7,y7,s,t)%調(diào)用導函數(shù)N_s,N_t=DHS(s,t);%求Jacobi矩陣J=

12、Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t);%求應變矩陣BB=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;endB=B/det(J);function J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t)%-Jacobi-%單元坐標%2,4,6,8點的坐標 x2=(x1+x3)/2;y2=(y1+y3)/2;x4=(x3+x5)/2;y4=

13、(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=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);endJ=x_s y_s;x_t y_t;fun

14、ction N=shape(s,t)%,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)/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)*(-

15、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)*(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)-

16、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ù)2.1E11 0.3 1 5 28 8 1 31 2 3 13 20 19 18 123 4 5 14 22 21 20 135 6 7 15 24 23 22 147 8 9 16 26 25