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

1、精選優(yōu)質文檔-傾情為你奉上廣州大學 有限元方法與程序設計 學院： 土木工程學院 專業(yè)： 結構工程 姓名： 曾一凡 學號： * % 平面四邊形八節(jié)點等參元MATLAB程序% 變量說明&（2015級結構工程曾一凡）% YOUNG POISS THICK% 彈性模量 泊松比 厚度% NPOIN NELEM NVFIX NFORCE % 總結點數(shù) 單元數(shù) 約束結點個數(shù) 受力結點數(shù)% COORD LNODS FORCE% 結構節(jié)點整體坐標數(shù)組 單元定義數(shù)組 結點力數(shù)組% ALLFORCE FIXED HK DISP% 總體荷載向量 約束信息數(shù)組 總體剛度矩陣 結點位移向量 % 1本程序計算了節(jié)點

2、位移和單元中心應力并輸出到nonde8out.txt文本里% 2在第四頁舉了一個實例 用MATLAB算出結果再用ANSYS算出的結果對比%=主程序=format short e %設定輸出類型clear %清除內存變量FP1=fopen('nonde8.txt','rt'); %打開初始數(shù)據(jù)文件FP2=fopen('nonde8out.txt','wt'); %打開文件 存放計算結果NPOIN=fscanf(FP1,'%d',1); %結點數(shù)NELEM=fscanf(FP1,'%d',1); %單元

3、數(shù)NFORCE=fscanf(FP1,'%d',1); %作用荷載的結點個數(shù)NVFIX=fscanf(FP1,'%d',1); %約束數(shù)YOUNG=fscanf(FP1,'%e',1); %彈性模量POISS=fscanf(FP1,'%f',1); %泊松比THICK=fscanf(FP1,'%f',1); %厚度LNODS=fscanf(FP1,'%d',8,NELEM)' %單元結點號數(shù)組（逆時針）COORD=fscanf(FP1,'%f',2,NPOIN)' %

4、 結點號 x,y坐標(整體坐標下)FORCE=fscanf(FP1,'%f',3,NFORCE)' %結點力數(shù)組% 節(jié)點力：結點號、X方向力(向右正)，Y方向力(向上正)FIXED=fscanf(FP1,'%d',NVFIX)' % 約束信息：約束對應的位移編碼（共計 NVFIX 組）EK=zeros(2*8,2*8); % 單元剛度矩陣并清零HK=zeros(2*NPOIN,2*NPOIN); % 張成總剛矩陣并清零X=zeros(1,8); %存放單元8個x方向坐標的向量并清零Y=zeros(1,8); %存放單元8個y方向坐標的向量并清零%

5、-求總剛-for i=1:NELEM % 對單元個數(shù)循環(huán) for m=1:8% 對單元結點個數(shù)循環(huán) X(m)=COORD(LNODS(i,m),1); %單元8個x方向坐標的向量 Y(m)=COORD(LNODS(i,m),2); %單元8個y方向坐標的向量 endEK=eKe(X,Y,YOUNG,POISS,THICK); %調用單元剛度矩陣a=LNODS(i,:); %臨時向量,用來記錄當前單元的結點編號 for j=1:8 %對行進行循環(huán)-按結點號循環(huán) for k=1:8 %對列進行循環(huán)-按結點號循環(huán) HK(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK(a

6、(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+. EK(j*2-1:j*2,k*2-1:k*2); % 單剛子塊疊加到總剛中 end endendALLFORCE=FOECEXL(NPOIN,NFORCE,FORCE); % 調用函數(shù)生成荷載向量%-處理約束-for j=1:NVFIX % 對約束個數(shù)進行循環(huán) N1=FIXED(j);HK(1:2*NPOIN,N1)=0; HK(N1,1:2*NPOIN)=0; HK(N1,N1)=1; % 將零位移約束對應的行、列變成零，主元變成1 ALLFORCE(N1)=0;endDISP=HKALLFORCE; % 方程求解，H

7、K先求逆再與力向量左乘得到位移%-求應力-stress=zeros(3,NELEM); % 應力向量并清零for i=1:NELEM % 對單元個數(shù)進行循環(huán) D=(YOUNG/(1-POISS*POISS)*1 POISS 0;POISS 1 0;0 0 (1-POISS)/2; %彈性矩陣 for k=1:8 % 對單元結點個數(shù)進行循環(huán) N2=LNODS(i,k); % 取單元結點號 U(k*2-1:k*2)=DISP(N2*2-1:N2*2); %從總位移向量中取出當前單元的結點位移 B=eBe(X,Y,0,0); %調用單元中心的應變矩陣 end stress(:,i)=D*B*U

8、9;end%=計算單元剛度矩陣函數(shù)=function EK=eKe(X,Y,YOUNG,POISS,THICK)EK=zeros(16,16); % 張成16*16矩陣并清零D=(YOUNG/(1-POISS*POISS)*1 POISS 0;POISS 1 0;0 0 (1-POISS)/2; %彈性矩陣%高斯積分采用3*3個積分點A1=5/9;A2=8/9;A3=5/9; %對應積分的加權系數(shù)A=A1 A2 A3;r=(3/5)(1/2); x=-r 0 r; %積分點for i=1:3 for j=1:3 B=eBe(X,Y,x(i),x(j); %調用應變矩陣B J=Jacobi(X,

9、Y,x(i),x(j); %調用雅可比矩陣J EK=EK+A(i)*A(j)*B'*D*B*det(J)*THICK; endend%=計算雅可比矩陣函數(shù)=function J=Jacobi(X,Y,s,t)N_s,N_t=DHS(s,t);x_s=0;y_s=0;x_t=0;y_t=0;for j=1:8 x_s=x_s+N_s(j)*X(j);y_s=y_s+N_s(j)*Y(j); x_t=x_t+N_t(j)*X(j);y_t=y_t+N_t(j)*Y(j);endJ=x_s y_s;x_t y_t;%=計算應變矩陣函數(shù)=function B=eBe(X,Y,s,t)N_s,N

10、_t=DHS(s,t);J=Jacobi(X,Y,s,t);B=zeros(3,16);for i=1:8 B1=J(2,2)*N_s(i)-J(1,2)*N_t(i); B2=-J(2,1)*N_s(i)+J(1,1)*N_t(i); B(1:3,2*i-1:2*i)=B1 0;0 B2;B2 B1;endB=B/det(J);%=計算形函數(shù)的求導函數(shù)= =function N_s,N_t=DHS(s,t)N_s(1)=1/4*(1-t)*(s-t-1)+1/4*(1+s)*(1-t);N_s(2)=1/2*(1+t)*(1-t);N_s(3)=1/4*(1+t)*(s+t-1)+1/4*(

11、1+s)*(1+t);N_s(4)=1/2*(1-s)*(1+t)-1/2*(1+s)*(1+t);N_s(5)=-1/4*(1+t)*(-s+t-1)-1/4*(1-s)*(1+t);N_s(6)=-1/2*(1+t)*(1-t);N_s(7)=-1/4*(1-t)*(-s-t-1)-1/4*(1-s)*(1-t);N_s(8)=1/2*(1-s)*(1-t)-1/2*(1+s)*(1-t);N_t(1)=-1/4*(1+s)*(s-t-1)-1/4*(1+s)*(1-t);N_t(2)=1/2*(1+s)*(1-t)-1/2*(1+s)*(1+t);N_t(3)=1/4*(1+s)*(s+

12、t-1)+1/4*(1+s)*(1+t);N_t(4)=1/2*(1+s)*(1-s);N_t(5)=1/4*(1-s)*(-s+t-1)+1/4*(1-s)*(1+t);N_t(6)=1/2*(1-s)*(1-t)-1/2*(1-s)*(1+t);N_t(7)=-1/4*(1-s)*(-s-t-1)-1/4*(1-s)*(1-t);N_t(8)=-1/2*(1+s)*(1-s);end%=計算總荷載矩陣函數(shù)=function ALLFORCE=FOECEXL(NPOIN,NFORCE,FORCE) % 本函數(shù)生成荷載向量 ALLFORCE=zeros(2*NPOIN,1); % 張成特定大小

13、的向量，并賦值0 for i=1:NFORCE ALLFORCE(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);%FORCE(i,1)為作用點,FORCE(i,2:3)為x,y方向的結點力 end%-輸出節(jié)點位移和單元中心應力-for i=1:NPOIN fprintf(FP2,'x%d=%dn',i, DISP(2*i-1); %輸出結點x方向的位移 fprintf(FP2,'y%d=%dn',i, DISP(2*i); %輸出結點y方向的位移endfor j=1:NELEM fprintf(FP2,'%d x=

14、%fn',j,stress(1,j); %輸出單元x方向的應力 fprintf(FP2,'%d y=%fn',j,stress(2,j); %輸出單元y方向的應力 fprintf(FP2,'%d xy=%fn',j,stress(3,j); %輸出單元切應力end%-實例計算并用ANSYS進行對比結果-如圖所示一個4m*1m懸臂梁，在3節(jié)點作用1*105N的豎向力，參數(shù)如下：彈性模量2.0*108，泊松比0.3，劃分成四個單元，每個單元八個節(jié)點，單元尺寸是1m*1m。1.用MATLAB進行計算在MATLAB當前工作目錄下存入初始數(shù)據(jù)文件，nonde.tx

15、t文本數(shù)據(jù)內容在表第1列，計算結果（位移和應力）在表第2、3、列，第4列為ANSYS計算的結點豎向位移23 4 1 6 2E08 0.3 11 2 3 4 5 6 7 87 6 5 9 10 11 12 1312 11 10 14 15 16 17 1817 16 15 19 20 21 22 234 04 0.54 13.5 13 13 0.53 03.5 02.5 12 12 0.5 2 02.5 01.5 11 1 1 0.51 01.5 00.5 10 10 0.50 00.5 03 0 -1E0539 40 41 42 43 44x1=-2.e-02y1=-1.e-01x2=-8.e

16、-05y2=-1.e-01x3=2.e-02y3=-1.e-01x4=2.e-02y4=-1.e-01x5=2.e-02y5=-8.e-02x6=1.e-05y6=-8.e-02x7=-2.e-02y7=-8.e-02x8=-2.e-02y8=-1.e-01x9=2.e-02y9=-6.e-02x10=1.e-02y10=-4.e-02x11=-3.e-06y11=-4.e-02x12=-1.e-02y12=-4.e-02x13=-2.e-02y13=-6.e-02x14=1.e-02y14=-2.e-02x15=1.e-02y15=-1.e-02x16=5.e-07y16=-1.e-02x1

17、7=-1.e-02y17=-1.e-02x18=-1.e-02y18=-2.e-02x19=5.e-03y19=-3.e-03x20=0y20=0x21=0y21=0x22=0y22=0x23=-5.e-03y23=-3.e-03-應力-1x=-18069.1y=8572.1xy=-83189.2x=3732.2y=-1485.2xy=-87734.3x=-669.3y=275.3xy=-92666.4x=98.4y=-40.4xy=-67107.ANSYS 豎向位移 1 -0.129512 -0.130633 -0.132164 -0.106345 -0.83270E-016 -0.82963E-017 -0.83027E-018 -0.10684 9 -0.61320E-0110 -0.41589E-0111 -0.41216E-0112 -0.41644E-0113 -0.612