有限元单元刚度矩阵
实现功能：输入欲求的单元号,得到该单元的面积和该单元刚度矩阵；输入完所有的单元号,得到所有的单元刚度矩阵后,直接得到半带宽存储的数组。（所有的单元刚度矩阵和半带宽数组分别存在所有单元的刚度矩阵.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
for (i=0;i
for (j=0;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
for(j=0;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
for(j=0;j
for(k=0,S[i][j]=0;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
{for(j=0;j
{printf("%7.4f ",S[i][j]);
fprintf(fp4,"%7.4f ",S[i][j]);}
fprintf(fp4,"\n");printf("\n");}
for(i=0;i
for(j=0;j
s1[j][i]=S[i][j];
for(i=0;i
{for(j=0;j
for(k=0, K[i][j]=0;k
K[i][j]=K[i][j]+s1[i][k]*B[k][j]*ae*t;}
printf("第%d单元的单元刚度阵为:\n",a);
fprintf(fp4,"\n第%d单元的单元刚度阵为\n",a);
for(i=0;i
{for(j=0;j
{ printf("%7.4f ",K[i][j]);
fprintf(fp4,"%7.4f ",K[i][j]);}
fprintf(fp4,"\n");printf("\n");}
for(i=0;i
{
}
}
for(i=0;i
}
fclose(fp1);fclose(fp2);fclose(fp3);fclose(fp4);fclose(fp5);fclose(fp7);//fclose(fp6);
}
for(j=0;j=xo&&(yj-xo)=xj&&(yj-xj)=xj&&(yo-xj)=xo&&(yo-xo)