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

四邊形八節(jié)點等參元matlab程序四邊形八節(jié)點等參元matlab程序四邊形八節(jié)點等參元matlab程序四邊形八節(jié)點等參元matlab程序編制僅供參考審核批準生效日期地址：電話：傳真:郵編:懸臂鋼梁，尺寸如圖一所示；v=。h=1，E=.圖一懸臂鋼梁圖二單元劃分與結(jié)點編號Matlab輸出結(jié)果附錄Ⅰ：有限元ANSYS分析結(jié)果采用PLANE183單元（四邊形八節(jié)點）單元得出的結(jié)構(gòu)Y向最大位移為。約等于matlab平面四邊形八節(jié)點等參元結(jié)點Y向最大位移。附錄Ⅱ：%---------------四邊形八節(jié)點等參元matlab計算程序----------------------------%———————————主程序—————————%*******************************************************************%************************************%2012年%本程序只能處理集中荷載作用下的情況%只輸出了節(jié)點位移、單元中心點的應(yīng)力%*******************************************************************%***************%變量說明%Evh%彈性模量泊松比厚度%NPOINNELEMNVFIXNNODENFPOIN%總結(jié)點數(shù),單元數(shù),約束結(jié)點個數(shù),單元節(jié)點數(shù),受力結(jié)點數(shù)%COORDLNODS%結(jié)構(gòu)節(jié)點整體坐標數(shù)組,單元定義數(shù)組,%FPOINFORCEFIXED%結(jié)點力數(shù)組，總體荷載向量,約束信息數(shù)組%HKDISP%總體剛度矩陣,結(jié)點位移向量%******************************clearallformatshorteFP1=fopen('','rt');%打開數(shù)據(jù)文件%%讀入控制數(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=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%*******************************************************************%*******************************************************************%========平面應(yīng)力問題的求解==============%%*******************************************************************%*******************************************************************%—————————————————————%剛度矩陣的生成%計算剛度矩陣，并對約束條件進行處理Ke=zeros(2*NNODE,2*NNODE);%單元剛度矩陣并清零HK=zeros(2*NPOIN,2*NPOIN);%張成總剛矩陣并清零%調(diào)用子程序生成單元剛度矩陣form=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(m,:);%臨時向量,用來記錄當前單元的節(jié)點編號%對總剛度矩陣的處理forj=1:8fork=1:8HK((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);endendend%—————————————————————————————————%對荷載向量進行處理FORCE=zeros(2*NPOIN,1);%張成總荷載向量并清零fori=1:NFPOINb1=FPOIN(i,1)*2-1;b2=FPOIN(i,1)*2;%FPION(i,1)為作用點FORCE(b1)=FPOIN(i,2);%FPION(i,2)為x方向的節(jié)點力FORCE(b2)=FPOIN(i,3);%FPION(i,3)為y方向的節(jié)點力end%—————————————————————————————————%將約束信息加入總剛,總荷載fori=1:NVFIXifFIXED(i,2)==1c1=2*FIXED(i,1)-1;HK(c1,:)=0;%將一約束序號處的總剛列向量清0HK(:,c1)=0;%將一約束序號處的總剛行向量清0HK(c1,c1)=1;%將行列交叉處的元素置為1FORCE(c1)=0;endifFIXED(i,3)==1c2=2*FIXED(i,1);HK(c2,:)=0;HK(:,c2)=0;HK(c2,c2)=1;FORCE(c2)=0;endend%—————————————————————————————————%===========================================================%===========================================================DISP=HK\FORCE%計算節(jié)點位移向量%===========================================================%===========================================================%———————————求解單元應(yīng)力————————————————stress=zeros(3,NELEM);form=1:NELEMu(1:16)=0;d=LNODS(m,:);%臨時向量,用來記錄當前單元的節(jié)點編號fori=1:NNODEu(i*2-1:i*2)=DISP(d(i)*2-1:d(i)*2);%從總位移向量中取出當前單元的節(jié)點位移endD=(E/(1-v*v))*[1v0;v10;00(1-v)/2];%彈性矩陣%形成應(yīng)變矩陣BMBM=zeros(3,16);fori=1:NNODEJ=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(1:3,2*i-1:2*i)=[B1i0;0B2i;B2iB1i]/det(J);endstressm=D*BM*u';stress(:,m)=stressm;endstress%輸出應(yīng)力functionKe=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é)點的坐標%矩陣尺寸為16x16Ke=zeros(16,16);D=(E/(1-v*v))*[1v0;v10;00(1-v)/2];%彈性矩陣%高斯積分采用3x3個積分點書74頁W1=5/9;W2=8/9;W3=5/9;%加權(quán)系數(shù)W=[W1W2W3];r=15^(1/2)/5;x=[-r0r];%積分點fori=1:3forj=1:3B=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;endendfunctionB=eleB(x1,y1,x3,y3,x5,y5,x7,y7,s,t)%調(diào)用導(dǎo)函數(shù)[N_s,N_t]=DHS(s,t);%求Jacobi矩陣J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t);%求應(yīng)變矩陣BB=zeros(3,16);fori=1:8B1i=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)=[B1i0;0B2i;B2iB1i];endB=B/det(J);functionJ=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=(y3+y5)/2;x6=(x5+x7)/2;y6=(y5+y7)/2;x8=(x7+x1)/2;y8=(y7+y1)/2;x=[x1x2x3x4x5x6x7x8];y=[y1y2y3y4y5y6y7y8];%%調(diào)用形函數(shù)對局部坐標的導(dǎo)數(shù)[N_s,N_t]=DHS(s,t);%求Jacobi矩陣的行列式的值x_s=0;y_s=0;x_t=0;y_t=0;fori=1:8x_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_sy_s;x_ty_t];functionN=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ù)求導(dǎo)%ξ,η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);