有限元单元刚度矩阵.doc

实现功能：输入欲求的单元号,得到该单元的面积和该单元刚度矩阵；输入完所有的单元号,得到所有的单元刚度矩阵后,直接得到半带宽存储的数组。（所有的单元刚度矩阵和半带宽数组分别存在所有单元的刚度矩阵.txt和SK矩阵.txt）第1单元面积为： 0.5000 弹性模量、泊松比和厚度分别为:100.0000 0.3000 0.1000单元1单元的应力矩阵-109.8901 -32.9670 109.8901 0.0000 0.0000 32.9670 -32.9670 -109.8901 32.9670 0.0000 0.0000 109.8901 -38.4615 -38.4615 0.0000 38.4615 38.4615 0.0000 第1单元的单元刚度阵为 7.4176 3.5714 -5.4945 -1.9231 -1.9231 -1.6484 3.5714 7.4176 -1.6484 -1.9231 -1.9231 -5.4945 -5.4945 -1.6484 5.4945 0.0000 0.0000 1.6484 -1.9231 -1.9231 0.0000 1.9231 1.9231 0.0000 -1.9231 -1.9231 0.0000 1.9231 1.9231 0.0000 -1.6484 -5.4945 1.6484 0.0000 0.0000 5.4945 半带宽存储的数组 7.4176 3.5714 -1.9231 -1.6484 -5.4945 -1.9231 7.4176 -1.9231 -5.4945 -1.6484 -1.9231 0.0000 9.8901 -2.3810 -4.9451 4.7619 -3.0220 -0.4579 9.8901 4.7619 -4.9451 -0.7326 0.5495 0.0000 20.8791 -2.3810 0.5495 -2.3810 -10.989 1.6484 13.7363 -2.3810 -3.0220 1.9231 -3.8462 0.0000 3.4341 1.1905 -0.9615 1.9231 0.0000 0.0000 5.2198 1.6484 -2.7473 0.0000 0.0000 0.0000 11.9506 -3.5714 0.0000 0.0000 0.0000 0.0000 6.5934 0.0000 0.0000 0.0000 0.0000 0.0000 #include #define jds 5#define d 2main()int nj,ne,i,j,k,a,n,*lnd; int JDS3;xj,yj,xo,yo; float B36,b136;D33,E,u,t;S36,s163;K66;Kk2*jds2*jds=0; float SK2*jds(d+1)*2=0; float x1,x2,x3,y1,y2,y3,ae,*xy; FILE *fp1,*fp2,*fp3,*fp4,*fp5,*fp7; fp1=fopen(节点坐标.txt,r); fp2=fopen(单元节点编号.txt,r);fp3=fopen(材料参数.txt,r); fp4=fopen(输出结果.txt,a); fp5=fopen(节点总数和单元总数.txt,r);fp7=fopen(SK矩阵.txt,w);fscanf(fp5,%d,&nj);fscanf(fp5,%d,&ne);xy=(float*)malloc(nj*2*sizeof(float); lnd=(int*)malloc(ne*3*sizeof(int); for(i=0;inj;i+) for(j=0;j2;j+) fscanf(fp1,%f,&xyi*2+j); for (i=0;ine;i+) for (j=0;j3;j+)fscanf(fp2,%d,&lndi*3+j);for(n=0;nne;n+)printf(请输入所求的单元号:);scanf(%d,&a);i=lnd(a-1)*3+0,j=lnd(a-1)*3+1,k=lnd(a-1)*3+2;JDS0=i;JDS1=j;JDS2=k;x1=xy(i-1)*2+0;y1=xy(i-1)*2+1;x2=xy(j-1)*2+0;y2=xy(j-1)*2+1;x3=xy(k-1)*2+0;y3=xy(k-1)*2+1;ae=0.5*(x2-x1)*(y3-y1)-(x3-x1)*(y2-y1);printf(第%d单元的面积为：%7.4fn,a,ae);fprintf(fp4,n第%d单元面积为：%7.4f ,a,ae);b100=y2-y3,b102=y3-y1;b104=y1-y2,b111=x3-x2;b113=x1-x3,b115=x2-x1;b120=b111,b121=b100;b122=b113,b123=b102;b124=b115,b125=b104;b110=0,b101=0, b112=0,b103=0,b114=0,b105=0;for(i=0;i3;i+) for(j=0;j6;j+) Bij=(0.5/ae)*b1ij;fscanf(fp3,%f %f %f,&E,&u,&t);printf(弹性模量、泊松比和厚度分别为:%7.4f %7.4f %7.4fn,E,u,t);fprintf(fp4,n弹性模量、泊松比和厚度分别为:%7.4f %7.4f %7.4f,E,u,t);D00=E/(1-u*u), D01=E*u/(1-u*u), D02=0;D10=E*u/(1-u*u),D11=E/(1-u*u), D12=0;D20=0, D21=0, D22=E*(1-u)/(2*(1-u*u);for(i=0;i3;i+) for(j=0;j6;j+) for(k=0,Sij=0;k3;k+) Sij=Sij+Dik*Bkj;printf(第%d单元的应力矩阵为:n,a);fprintf(fp4,n单元%d单元的应力矩阵n,a);for(i=0;i3;i+) for(j=0;j6;j+)printf(%7.4f ,Sij);fprintf(fp4,%7.4f ,Sij); fprintf(fp4,n);printf(n);for(i=0;i3;i+) for(j=0;j6;j+) s1ji=Sij;for(i=0;i6;i+)for(j=0;j6;j+) for(k=0, Kij=0;k3;k+) Kij=Kij+s1ik*Bkj*ae*t;printf(第%d单元的单元刚度阵为:n,a);fprintf(fp4,n第%d单元的单元刚度阵为n,a);for(i=0;i6;i+)for(j=0;j6;j+) printf(%7.4f ,Kij);fprintf(fp4,%7.4f ,Kij); fprintf(fp4,n);printf(n);for(i=0;i6;i+)k=i/2;xj=JDSk*2-1;xo=JDSk*2;if(i%2=0)for(j=0;j=xj&(yj-xj)=xj&(yo-xj)(d+1)*2)SKxj-1yo-xj=SKxj-1yo-xj+Kij;elsefor(j=0;j=xo&(yj-xo)=xo&(yo-xo)(d+1)*2)SKxo-1yo-xo=SKxo-1yo-xo+Kij;f