#include <io.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
void hbw();
void sncs(int nel);
void fix(int np);
void trmat();
void fis(int nel);
void fpj();
void force();
void stiff();
void addsm();
void restr();
void matmul();
void soleq();
void outdis();
float sm[6][6], tsm[300][30],res[60][2],pj[300],t[6][6],d[6][6],r[300],
fo[6],foj[6],pf[200][4],x[100],y[100],ae[10][3],sl,sn,cs,eal,eil;
int nj,ne,nt,nr,npj,npf,nn,mbw,jel[100][2],nae[100],is[6];
FILE *infile,*outfile;
/******主函数*****/
int main()
{ char name1[30],name2[30];
int i,j,nel,np;
printf("程序汉化与输入优化\n");
printf("请输入数据文件名（请使用txt文件格式，并不要输入后缀名）\n");
scanf("%s",name1);
printf("请输出输出文件名，后缀名将会自动补全\n");
scanf("%s",name2);
if((infile=fopen(strcat(name1,".txt"),"r"))!= NULL)
{ fscanf(infile,"%d%d%d%d%d%d",&nj,&ne,&nt,&nr,&npj,&npf);
for(i=0;i<ne;i++) fscanf(infile,"%d%d%d",&jel[i][0],&jel[i][1],&nae[i]);
for(i=0;i<nj;i++) fscanf(infile,"%f%f",&x[i],&y[i]);
for(i=0;i<nt;i++) fscanf(infile,"%f%f%f",&ae[i][0],&ae[i][1],&ae[i][2]);
for(i=0;i<nr;i++) fscanf(infile,"%f%f",&res[i][0],&res[i][1]); }
else
{ printf("未找到%s!",name1);
exit(1);}
nn=3*nj;
outfile=fopen(strcat(name2,".txt"),"w");
if(outfile==NULL)
{ printf("结果文件输出失败!"); exit(1);}
fprintf(outfile,"平面刚架分析结果文件——汉化\n");
fprintf(outfile,"\n读入数据\n");
fprintf(outfile,"*******\n");
fprintf(outfile,"控制数据\n");
fprintf(outfile,"结点总数:%3d\n",nj);
fprintf(outfile,"单元总数:%3d\n",ne);
fprintf(outfile,"单元性质总数:%3d\n",nt);
fprintf(outfile,"约束总数:%3d\n",nr);
fprintf(outfile,"有荷载结点总数:%3d\n",npj);
fprintf(outfile,"非结点荷载数:%3d\n",npf);
fprintf(outfile,"结点自由度总数:%3d\n",nn);
fprintf(outfile,"\n以下为单元结点号与性质：\n");
fprintf(outfile,"\n 单元号 起点坐标 终点坐标 单元性质编号\n");
for(i=0;i<ne;i++)
fprintf(outfile,"%7d%13d%15d%15d\n",i+1,jel[i][0],jel[i][1],nae[i]);
fprintf(outfile,"\n结点的坐标列表：\n");
fprintf(outfile," 结点号 x y\n");
for(i=0;i<nj;i++)
fprintf(outfile,"%10d%10.2f%10.2f\n",i+1,x[i],y[i]);
fprintf(outfile,"\n单元截面类型性质列表：\n");
fprintf(outfile,"\n 种类号 面积 惯性矩 弹性模量\n");
for(i=0;i<nt;i++)
fprintf(outfile,"%8d%15.6f%15.6f%20.6f\n",i+1,ae[i][0],ae[i][1],ae[i][2]);
fprintf(outfile,"\n约束列表\n");
fprintf(outfile,"\n 约束编号 约束位置 约束方向位移/转角\n");
for(i=0;i<nr;i++)
fprintf(outfile,"%10d%25.0f%19.3f\n",i+1,res[i][0],res[i][1]);
hbw();
for(i=0;i<nn;i++) pj[i]=0.0;
if (npj==0) goto aa;
fprintf(outfile,"\n结点荷载列表：\n");
fprintf(outfile,"\n结点号 x向分荷载 y向分荷载 弯矩\n");
for(i=0;i<npj;i++)
{ int no;
float px,py,zm;
fscanf(infile,"%d%f%f%f",&no,&px,&py,&zm);
fprintf(outfile,"%6d %12.2f %12.2f%12.2f\n",no,px,py,zm);
pj[3*no-3]=px; pj[3*no-2]=py; pj[3*no-1]=zm; }
aa: for(i=0;i<nn;i++) r[i]=0.0;
for(i=0;i<nr;i++)
{ int ni;
ni=floor(res[i][0]+0.1)-1;
if(pj[ni]!=0.0) r[ni]=r[ni]-pj[ni]; }
if(npf==0) goto bb;
fprintf(outfile,"\n非结点荷载列表：\n");
fprintf(outfile,"\n 荷载大小 作用点/范围 单元号 荷载类型\n");
for(i=0;i<npf;i++)
{ fscanf(infile,"%f%f%f%f",&pf[i][0],&pf[i][1],&pf[i][2],&pf[i][3]);
fprintf(outfile,"%15.2f%15.2f%15.0f%15.0f\n",pf[i][0],pf[i][1],pf[i][2],pf[i][3]);}
for(np=0;np<npf;np++)
{ nel=floor(pf[np][2]+0.1);
sncs(nel-1);
fix(np);
trmat();
fis(nel-1);
fpj(); }
bb:for(i=0;i<nn;i++)
for(j=0;j<mbw;j++) tsm[i][j]=0.0;
for(nel=0;nel<ne;nel++)
{ sncs(nel);
trmat();
stiff();
matmul();
fis(nel);
addsm(); }
restr();
soleq();
outdis();
force();
fclose(outfile);
system(name2);
return 0;
}
/*求最大半带宽*/
void hbw()
{ int nel;
mbw=0;
for( nel=0;nel<ne;nel++)
{ int ma=abs(jel[nel][0]-jel[nel][1]);
if(mbw<ma) mbw=ma; }
mbw=3*(mbw+1);
fprintf(outfile,"\n最大半带宽:%5d\n",mbw);
}
/*矩阵相乘*/
void matmul()
{ int i,j,k;
for(i=0;i<6;i++)
for(j=0;j<6;j++) d[i][j]=0.0;
for(i=0;i<6;i++)
for(j=0;j<6;j++)
for(k=0;k<6;k++) d[i][j]=d[i][j]+t[k][i]*sm[k][j];
for(i=0;i<6;i++)
for(j=0;j<6;j++) sm[i][j]=0.0;
for(i=0;i<6;i++)
for(j=0;j<6;j++)
for(k=0;k<6;k++)
sm[i][j]=sm[i][j]+d[i][k]*t[k][j];
}
/*求单元常数*/
void sncs (int nel)
{ int ii,jj,k;
float xi,yi,xj,yj;
ii=jel[nel][0];
jj=jel[nel][1];
xi=x[ii-1];
xj=x[jj-1];
yi=y[ii-1];
yj=y[jj-1];
sl=sqrt((xj-xi)*(xj-xi)+(yj-yi)*(yj-yi));
sn=(yj-yi)/sl;
cs=(xj-xi)/sl;
k=nae[nel];
eal=ae[k-1][0]*ae[k-1][2]/sl;
eil=ae[k-1][1]*ae[k-1][2]/sl;
}
/*求单元坐标转换矩阵*/
void trmat()
{ int i,j;
for(i=0;i<6;i++)
for(j=0;j<6;j++)
t[i][j]=0.0; t[0][0]=cs; t[0][1]=sn;
t[1][0]=-sn; t[1][1]=cs; t[2][2]=1.0;
for(i=0;i<3;i++)
for(j=0;j<3;j++) t[i+3][j+3]=t[i][j];
}
/*求单元固端力*/
void fix (int np)
{ float w,c,c1,c2,cc,cc1,cc2;
int i,im;
w=pf[np][0]; c=pf[np][1];
c1=c/sl; c2=c1*c1; cc=sl-c;
cc1=cc/sl; cc2=cc1*cc1;
for(i=0;i<6;i++) fo[i]=0.0;
im=floor(pf[np][3]+0.1);
if(im==1)
{ fo[1]=-w*c*(1.0-c2+c2*c1/2.0);
fo[2]=-w*c*c*(6.0-8.0*c1+3*c2)/12.0;
fo[4]=-w*c-fo[1];
fo[5]=w*c1*c*c*(4.0-3.0*c1)/12.0; }
else if (im==2)
{ fo[1]=-w*(1+2*c1)*cc2;
fo[2]=-w*c*cc2;
fo[4]=-w*(1+2*cc1)*c2;
fo[5]=w*cc*c2; }
else if (im==3)
{ fo[1]=6*w*c1*cc1/sl;
fo[2]=w*cc1*(2-3*cc1);
fo[4]=-fo[2];
fo[5]=w*c1*(2-3*c1); }
else if (im==4)
{ fo[1]=-0.25*w*c*(2.-3.*c2+1.6*c2*c1);
fo[2]=-w*c*c*(2.-3.*c1+1.2*c2)/6.0;
fo[4]=-w*c/2.0-fo[1];
fo[5]=0.25*w*c*c*c1*(1.-0.8*c1); }
else if (im==5)
{ fo[0]=-w*cc1;
fo[3]=-w*c1; }
else if (im==6)
{ fo[0]=-w*c*(1.-0.5*c1);
fo[3]=-0.5*w*c*c1; }
}
/*形成总结点荷载向量*/
void fpj()
{ int i,j,k,ii,jj;
for( k=0;k<6;k++)
foj[k]=0.0;
for(i=0;i<6;i++)
for(j=0;j<6;j++) foj[i]=foj[i]+t[j][i]*fo[j];
for(ii=0;ii<6;ii++) pj[is[ii]]=pj[is[ii]]-foj[ii];
}
/*形成单元刚度矩阵*/
void stiff()
{ int i,j;
float s1,s2;
for(i=0;i<6;i++)
for(j=0;j<6;j++)
sm[i][j]=0.0; sm[0][0]=eal;
sm[3][0]=-eal; sm[3][3]=eal;
s1=12.*eil/(sl*sl); s2=6.*eil/sl;
sm[1][1]=s1; sm[4][1]=-s1 ;
sm[4][4]=s1; sm[2][1]=s2;
sm[5][1]=s2; sm[4][2]=-s2;
sm[5][4]=-s2; sm[2][2]=4*eil;
sm[5][5]=4*eil; sm[5][2]=2*eil;
for(i=0;i<6;i++)
for(j=0;j<i;j++)
sm[j][i]=sm[i][j];
}
/*由单元位移分量码L形成总刚位移分量码IS(L)*/
void fis(int nel)
{ int i,j;
for(i=0;i<2;i++)
for(j=0;j<3;j++) is[3*i+j]=3*(jel[nel][i]-1)+j; }
/*形成结构原始总刚度矩阵*/
void addsm()
{ int i,j;
for( i=0;i<6;i++)
for(j=0;j<6;j++)
{ int kc;
kc=is[j]-is[i];
if(kc>=0) tsm[is[i]][kc]=tsm[is[i]][kc]+sm[i][j]; }
}
/*约束处理*/
void restr()
{ int i;
for(i=0;i<nr;i++)
{ int ni;
ni=floor(res[i][0]+0.1); tsm[ni-1][0]=1.0e25;
pj[ni-1]=tsm[ni-1][0]*res[i][1]; }
}
/*解线性方程组*/
void soleq()
{ float c1;
int k,ni,im,i,m,j,nm,jm;
for(k=1;k<nn;k++)
{ if (fabs(tsm[k-1][0])<=0.000001)
{ fprintf(outfile,"*****%4d行数据异常******\n",k+1);
goto fin; }
ni=k+mbw-1; im=ni;
if(ni>nn) im=nn;
for(i=k+1;i<im+1;i++)
{ m=i-k+1;
c1=tsm[k-1][m-1]/tsm[k-1][0];
for(j=1;j<mbw-m+2;j++)
tsm[i-1][j-1]=tsm[i-1][j-1]-c1*tsm[k-1][j+i-k-1];
pj[i-1]=pj[i-1]-c1*pj[k-1]; }
}
if(fabs(tsm[nn-1][0])<=0.000001)
{ fprintf(outfile,"*****%4d行数据异常******\n",nn);
goto fin; }
pj[nn-1]=pj[nn-1]/tsm[nn-1][0];
for(i=nn-1;i>0;i--)
{ nm=nn-i+1;
jm=nm;
if(nm>mbw) jm=mbw;
for(j=2;j<jm+1;j++)
pj[i-1]=pj[i-1]-tsm[i-1][j-1]*pj[j+i-2];
pj[i-1]=pj[i-1]/tsm[i-1][0]; }
fin: return;
}
/*输出位移向量*/
void outdis()
{ int i;
fprintf(outfile,"\n输出结果\n");
fprintf(outfile,"******\n");
fprintf(outfile,"\n结点位移列表：\n");
fprintf(outfile,"\n 结点号 x向 y向 转角\n");
for(i=0;i<nj;i++) fprintf(outfile,"%10d%15.4f%15.4f%15.4f\n",i+1,pj[3*i],pj[3*i+1],pj[3*i+2]);
}
/*求单元杆端力、支座反力或结点合力*/
void force ()
{ float dj[6],f[6],fj[6],dd[6];
int nel,np,i,j,ip;
fprintf(outfile,"\n单元号 轴力1 剪力1 弯矩1 轴力2 剪力2 弯矩2 \n");
for(nel=0;nel<ne;nel++){
sncs(nel);
trmat();
stiff();
fis(nel);
for(i=0;i<6;i++) dj[i]=pj[is[i]];
for(i=0;i<6;i++) {dd[i]=0.0;f[i]=0.0;}
for(i=0;i<6;i++)
for(j=0;j<6;j++) dd[i]=dd[i]+t[i][j]*dj[j];
for(i=0;i<6;i++)
for(j=0;j<6;j++) f[i]=f[i]+sm[i][j]*dd[j];
if(npf==0) goto a;
for(np=0;np<npf;np++){
ip=floor(pf[np][2]+0.1)-1;
if(ip==nel){ fix(np);
for(j=0;j<6;j++) f[j]=f[j]+fo[j];}
}
a: fprintf(outfile,"%5d%11.2f%11.2f%11.2f%11.2f%11.2f%11.2f\n",
nel+1,f[0], f[1],f[2],f[3],f[4],f[5]);
for(i=0;i<6;i++) fj[i]=0.0;
for(i=0;i<6;i++) fj[i]=0.0;
for(i=0;i<6;i++)
for(j=0;j<6;j++) fj[i]=fj[i]+t[j][i]*f[j];
for(i=0;i<6;i++) r[is[i]]=r[is[i]]+fj[i]; }
fprintf(outfile,"\n结点合力列表：\n");
fprintf(outfile,"\n 结点号 X向合力 Y向合力 合力矩\n");
for(i=0;i<nj;i++)
fprintf(outfile,"%10d%15.4f%15.4f%15.4f\n",i+1,r[3*i],r[3*i+1],r[3*i+2]);
}
输入文件示例：
3 2 1 4 1 2
1 2 1 2 3 1
0 0 0 5 5 5
0.008 0.002 2.5E7
1 0 2 0 7 0 9 0
2 0 0 12000
-6000 5 2 1
-8000 2.5 1 2
输出文件：
平面刚架分析结果文件——汉化
读入数据
控制数据
结点总数: 3
单元总数: 2
单元性质总数: 1
约束总数: 4
有荷载结点总数: 1
非结点荷载数: 2
结点自由度总数: 9
以下为单元结点号与性质：
单元号 起点坐标 终点坐标 单元性质编号
1 1 2 1
2 2 3 1
结点的坐标列表：
结点号 x y
1 0.00 0.00
2 0.00 5.00
3 5.00 5.00
单元截面类型性质列表：
种类号 面积 惯性矩 弹性模量
1 0.008000 0.002000 25000000.000000
约束列表
约束编号 约束位置 约束方向位移/转角
1 1 0.000
2 2 0.000
3 7 0.000
4 9 0.000
最大半带宽: 6
结点荷载列表：
结点号 x向分荷载 y向分荷载 弯矩
2 0.00 0.00 12000.00
非结点荷载列表：
荷载大小 作用点/范围 单元号 荷载类型
-6000.00 5.00 2 1
-8000.00 2.50 1 2
输出结果
结点位移列表：
结点号 x向 y向 转角
1 -0.0000 -0.0000 0.2000
2 0.2500 -0.7500 -0.8000
3 0.0000 -5.8750 -0.0000
单元号 轴力1 剪力1 弯矩1 轴力2 剪力2 弯矩2
1 30000.00 -2000.00 0.00 -30000.00 10000.00 -30000.01
2 10000.00 30000.00 42000.00 -10000.00 0.00 33000.00
结点合力列表：
结点号 X向合力 Y向合力 合力矩
1 2000.0010 30000.0000 0.0000
2 0.0010 0.0000 11999.9941
3 -10000.0020 0.0000 33000.0039
程序解读：
以上代码和图片来自大连理工大学陈廷国教授结构力学专题课件
京公网安备 11010502036488号