有限元单元刚度矩阵

实现功能：输入欲求的单元号,得到该单元的面积和该单元刚度矩阵；输入完所有的单元号,得到所有的单元刚度矩阵后,直接得到半带宽存储的数组。（所有的单元刚度矩阵和半带宽数组分别存在所有单元的刚度矩阵.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 2

main()

{int nj,ne,i,j,k,a,n,*lnd;

int JDS[3];xj,yj,xo,yo;

float B[3][6],b1[3][6];D[3][3],E,u,t;S[3][6],s1[6][3];K[6][6];Kk[2*jds][2*jds]={{0}};

float SK[2*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;i

for(j=0;j<2;j++)

fscanf(fp1,"%f",&xy[i*2+j]);

for (i=0;i

for (j=0;j<3;j++)

fscanf(fp2,"%d",&lnd[i*3+j]);

for(n=0;n

printf("请输入所求的单元号:");

scanf("%d",&a);

i=lnd[(a-1)*3+0],j=lnd[(a-1)*3+1],k=lnd[(a-1)*3+2];

JDS[0]=i;JDS[1]=j;JDS[2]=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.4f\n",a,ae);

fprintf(fp4,"\n第%d单元面积为：%7.4f ",a,ae);

b1[0][0]=y2-y3,b1[0][2]=y3-y1;b1[0][4]=y1-y2,b1[1][1]=x3-x2;b1[1][3]=x1-x3,b1[1][5]=x2-x1;

b1[2][0]=b1[1][1],b1[2][1]=b1[0][0];b1[2][2]=b1[1][3],b1[2][3]=b1[0][2];

b1[2][4]=b1[1][5],b1[2][5]=b1[0][4];b1[1][0]=0,b1[0][1]=0, b1[1][2]=0,b1[0][3]=0,b1[1][4]=0,b1[0][5]=0; for(i=0;i<3;i++)

for(j=0;j<6;j++)

B[i][j]=(0.5/ae)*b1[i][j];

fscanf(fp3,"%f %f %f",&E,&u,&t);

printf("弹性模量、泊松比和厚度分别为:%7.4f %7.4f %7.4f\n",E,u,t);

fprintf(fp4,"\n弹性模量、泊松比和厚度分别为:%7.4f %7.4f %7.4f",E,u,t);

D[0][0]=E/(1-u*u), D[0][1]=E*u/(1-u*u), D[0][2]=0;

D[1][0]=E*u/(1-u*u),D[1][1]=E/(1-u*u), D[1][2]=0;

D[2][0]=0, D[2][1]=0, D[2][2]=E*(1-u)/(2*(1-u*u));

for(i=0;i<3;i++)

for(j=0;j<6;j++)

for(k=0,S[i][j]=0;k<3;k++)

S[i][j]=S[i][j]+D[i][k]*B[k][j];

printf("第%d单元的应力矩阵为:\n",a);

fprintf(fp4,"\n单元%d单元的应力矩阵\n",a);

for(i=0;i<3;i++)

{for(j=0;j<6;j++)

{printf("%7.4f ",S[i][j]);

fprintf(fp4,"%7.4f ",S[i][j]);}

fprintf(fp4,"\n");printf("\n");}

for(i=0;i<3;i++)

for(j=0;j<6;j++)

s1[j][i]=S[i][j];

for(i=0;i<6;i++)

{for(j=0;j<6;j++)

for(k=0, K[i][j]=0;k<3;k++)