平面四边形八节点等参元matlab程序

1、泊广州大学有限元方法与程序设计学院：土木工程学院专业：结构工程姓名：曾一凡学号：21115160*%平面四边形八节点等参元MATLABS序%变量说明&（2015级一一结构工程一一曾一凡）%YOUNGPOISSTHICK%弹性模量泊松比厚度%NPOINNELEMNVFIXNFORCE%总结点数单元数约束结点个数受力结点数%COORDLNODSFORCE%结构节点整体坐标数组单元定义数组结点力数组%ALLFORCEFIXEDHKDISP%总体荷载向量约束信息数组总体刚度矩阵结点位移向量%1本程序计算了节点位移和单元中心应力并输出到nonde8out.txt文本里%2在第四页举了一个实例用M

2、ATLABJ出结果再用ANSYSt出的结果对比%=§序formatshorte%设定输出类型clear%青除内存变量FP1=fopen('nonde8.txt','rt');%丁开初始数据文件FP2=fopen('nonde8out.txt','wt');%丁开文件存放计算结果NPOIN=fscanf(FP1,'%d',1);%结点数NELEM=fscanf(FP1,'%d',1);%11元数NFORCE=fscanf(FP1,'%d',1);%作用荷载的结点个数NVFIX

3、=fscanf(FP1,'%d',1);%勺束数YOUNG=fscanf(FP1,'%e',1);%单性模量POISS=fscanf(FP1,'%f,1);柏白松比THICK=fscanf(FP1,'%f,1);%厚度LNODS=fscanf(FP1,'%d',8,NELEM)'%i1元结点号数组(逆时针)COORD=fscanf(FP1,'%f,2,NPOIN)'%结点号x,y坐标(整体坐标下)FORCE=fscanf(FP1,'%f,3,NFORCE)'%结点力数组%节点力：结点号、X方

4、向力(向右正)，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方向坐标的向量并清零%求总刚fori=1:NELEM%对单元个数循环form=1:8%对单元结点个数循环X(m尸COORD(LNODS(i,m),1);%单元8个x方向坐标的向量Y(m)=COORD(L

5、NODS(i,m),2);%单元8个y方向坐标的向量endEK=eKe(X,Y,YOUNG,POISS,THICK);%调用单元刚度矩阵a=LNODS(i,:);%临时向量，用来记录当前单元的结点编号forj=1:8%对行进彳T循环-按结点号循环fork=1:8%对列进彳T循环-按结点号循环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)+.EK(j*2-1:j*2,k*2-1:k*2);%单刚子块叠加到总刚中endendendALLFORCE=FOECEXL(NPOIN,NFORCE,F

6、ORCE);%调用函数生成荷载向量%处理约束forj=1:NVFIX%对约束个数进行循环N1=FIXED(j);HK(1:2*NPOIN,N1)=0;HK(N1,1:2*NPOIN)=0;HK(N1,N1)=1;%将零位移约束对应的行、列变成零，主元变成1ALLFORCE(N1)=0;endDISP=HKALLFORCE;%方程求解，HK先求逆再与力向量左乘得到位移%求应力stress=zeros(3,NELEM);%应力向量并清零fori=1:NELEM%对单元个数进行循环D=(YOUNG/(1-POISS*POISS)*1POISS0;POISS10;00(1-POISS)/2;%弹性矩阵

7、fork=1:8%对单元结点个数进行循环N2=LNODS(i,k);%取单元结点号U(k*2-1:k*2)=DISP(N2*2-1:N2*2);%从总位移向量中取出当前单元的结点位移B=eBe(X,Y,0,0);%调用单元中心的应变矩阵endstress(:,i)=D*B*U'end%="算单元刚度矩阵函数=functionEK=eKe(X,Y,YOUNG,POISS,THICK)EK=zeros(16,16);%张成16*16矩阵并清零D=(YOUNG/(1-POISS*POISS)*1POISS0;POISS10;00(1-POISS)/2;%弹性矩阵%高斯积分采用3*3

8、个积分点A1=5/9;A2=8/9;A3=5/9;%对应积分的加权系数A=A1A2A3;r=(3/5)A(1/2);x=-r0r;%K分点fori=1:3forj=1:3B=eBe(X,Y,x(i),x(j);%调用应变矩阵BJ=Jacobi(X,Y,x(i),x(j);%调用雅可比矩阵JEK=EK+A(i)*A(j)*B*D*B*det(J)*THICK;endend%=H*算雅可比矩阵函数=functionJ=Jacobi(X,Y,s,t)N_s,N_t=DHS(s,t);x_s=0;y_s=0;x_t=0;y_t=0;forj=1:8x_s=x_s+N_s(j)*X(j);y_s=y_s

9、+N_s(j)*Y(j);x_t=x_t+N_t(j)*X(j);y_t=y_t+N_t(D*Y(j);endJ=x_sy_s;x_ty_t;%=!=算应变矩阵函数=functionB=eBe(X,Y,s,t)N_s,N_t=DHS(s,t);J=Jacobi(X,Y,s,t);B=zeros(3,16);fori=1:8B1=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)=B10;0B2;B2B1;endB=B/det(J);%=i="算形函数的求导函数=functionN_s

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

11、(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+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);e

12、nd%="算总荷载矩阵函数=functionALLFORCE=FOECEXL(NPOIN,NFORCE,FORCE)%本函数生成荷载向量ALLFORCE=zeros(2*NPOIN,1);%张成特定大小的向量，并赋值0fori=1:NFORCEALLFORCE(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);%FORCE(i,1历作用点,FORCE(i,2:3)为x,y方向的结点力end%输出节点位移和单元中心应力fori=1:NPOINfprintf(FP2,'x%d=%dn',i,DISP(2*i-1);%俞出结点x方向的位移

13、fprintf(FP2,'y%d=%dn',i,DISP(2*i);%俞出结点y方向的位移endforj=1:NELEMfprintf(FP2,'%dx=%fn',j,stress(1,j);%输出单元x方向的应力fprintf(FP2,'%dy=%fn',j,stress(2,j);%俞出单元y方向的应力fprintf(FP2,'%dxy=%fn',j,stress(3,j);%俞出单元切应力end%实例计算并用ANSYS!行对比结果如图所示一个4m*1m悬臂梁，在3节点作用1*105N的竖向力，参数如下：弹性模量2.0*108

14、,泊松比0.3,划分成四个单元，每个单元八个节点，单元尺寸是1m*1m1.用MATLA避行计算在MATLA由前工作目录下存入初始数据文件，nonde.txt文本数据内容在表第1歹计算结果(位移和应力)在表第2、3、歹I，第4列为ANSYSt算的结点竖向位移234162E080.31x1=-2.335725e-02x16=5.504682e-07ANSYS竖向位移12345678y1=-1.298092e-01y16=-1.135171e-021-0.12951765910111213x2=-8.552411e-05x17=-1.009374e-022-0.1306312111014151617

15、18y2=-1.304313e-01y17=-1.224831e-023-0.132161716151920212223x3=2.414924e-02x18=-1.428159e-024-0.1063440y3=-1.314740e-01y18=-2.498899e-025-0.83270E-0140.5x4=2.339276e-02x19=5.239181e-036-0.82963E-0141y4=-1.063855e-01y19=-3.600960e-037-0.83027E-013.51x5=2.212796e-02x20=08-0.1068431y5=-8.301568e-02y20=

16、09-0.61320E-0130.5x6=1.768304e-05x21=010-0.41589E-0130y6=-8.281498e-02y21=011-0.41216E-013.50x7=-2.217334e-02x22=012-0.41644E-0152.512120.5202.501.511110.5101.500.510100.5000.5030-1E05394041424344y7=-8.298558e-02x8=-2.314206e-02y8=-1.064554e-01x9=2.026594e-02y9=-6.115975e-02x10=1.770574e-02y10=-4.15

17、3420e-02x11=-3.209328e-06y11=-4.112391e-02x12=-1.769400e-02y12=-4.152824e-02x13=-2.028457e-02y13=-6.114672e-02x14=1.428595e-02y14=-2.498661e-02x15=1.009174e-02y15=-1.224757e-02y22=0x23=-5.240129e-03y23=-3.600612e-03应力1x=-18069.7775521y=8572.1772051xy=-83189.4115272x=3732.7435712y=-1485.7694402xy=-87734.9047963x=-669.4351803y=275.0800593xy=-92666.3579854x=98.0150424y=-40.2620194xy=-67107.62081913-0.61222E-0114-0.25027E-0115-0.12282E-0116-0.11371E-0117-0.12285E-0118-0.25067E-0119-0.36203E-02200.0000210.0000220.000023-0.361