#include<io.h>
#include "string.h"
#include"stdlib.h"
#include"stdio.h"
#include"math.h"
void hbw();
void sncs(int nel);
void fix(int np);
void trmat();
void fis(int nel);
void fpj();
void force();
void stiff(int nel);
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],typ[100],is[6];
FILE *infile,*outfile;
/******主函数*****/
int main()
{ char name1[30],name2[30];
int i,j,nel,np;
printf("please enter data-filename\n");
scanf("%s",name1);
printf("please enter result-filename\n");
scanf("%s",name2);
if((infile=fopen(name1,"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%d",&jel[i][0],&jel[i][1],&nae[i],&typ[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("the data-file not exit!");
exit(1);}
nn=3*nj;
outfile=fopen(name2,"w");
if(outfile==NULL)
{ printf("the result-file not exist!"); exit(1);}
fprintf(outfile," statis analysis of plane frame\n");
fprintf(outfile,"input data\n");
fprintf(outfile,"*******\n");
fprintf(outfile,"control data\n");
fprintf(outfile,"the num. of node:%3d\n",nj);
fprintf(outfile,"the num. of mem:%3d\n",ne);
fprintf(outfile,"the num. of type of section characteristic:%3d\n",nt);
fprintf(outfile,"the num. of restricted degrees of freedom:%3d\n",nr);
fprintf(outfile,"the num. of nodal load:%3d\n",npj);
fprintf(outfile,"the num. of non-nodal loads:%3d\n",npf);
fprintf(outfile,"the num. of nodal degrees of freedom:%3d\n",nn);
fprintf(outfile,"information of mem.\n");
fprintf(outfile," mem. start node end node section type\n");
for(i=0;i<ne;i++)
fprintf(outfile,"%5d%10d%10d%10d%10d\n",i+1,jel[i][0],jel[i][1],nae[i],typ[i]);
fprintf(outfile,"coordinates x and y of nodes\n");
fprintf(outfile," node 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,"information of cross section each mem. \n");
fprintf(outfile," type area moment of interia elastic modulus\n");
for(i=0;i<nt;i++)
fprintf(outfile,"%8d%15.5f%15.5f%20.2f\n",i+1,ae[i][0],ae[i][1],ae[i][2]);
fprintf(outfile,"information restriction \n");
fprintf(outfile," rester.-no restr.-disp.-no restr.-disp.\n");
for(i=0;i<nr;i++)
fprintf(outfile,"%10d%19.3f%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,"nodal loads\n");
fprintf(outfile," node px py zm\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,"%10d%10.2f%10.2f%10.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,"non-nodal loads\n");
fprintf(outfile," loads leng.supp. to load mem type\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.2f%15.2f\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(nel);
matmul();
fis(nel);
addsm(); }
restr();
soleq();
outdis();
force();
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,"Half _Bandwidth Mbw:%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 nel)
{ int i,j,it;
float s1,s2;
it=floor(typ[nel]+0.1);
if(it==1)
{
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];
}
}
else if(it==2)
{
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=3.*eil/(sl*sl); s2=3.*eil/sl;
sm[1][1]=s1; sm[4][1]=-s1;
sm[4][4]=s1; sm[2][1]=s2;
sm[4][2]=-s2; sm[2][2]=3*eil;
}
}
for(i=0;i<6;i++)
{
for(j=0;j<i;j++)
sm[j][i]=sm[i][j];
}
}
else if(it==3)
{
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=3.*eil/(sl*sl); s2=3.*eil/sl;
sm[1][1]=s1; sm[4][1]=-s1;
sm[4][4]=s1;sm[5][1]=s2;
sm[5][4]=-s2; sm[5][5]=3*eil;
}
}
for(i=0;i<6;i++)
{
for(j=0;j<i;j++)
sm[j][i]=sm[i][j];
}
}
else if(it==4)
{
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;
}
}
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.0e30;
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,"*****singularity in row %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,"*****singularity in row %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,"output results\n");
fprintf(outfile,"******\n");
fprintf(outfile,"nodal displacements\n");
fprintf(outfile," node u v ceta\n");
for(i=0;i<nj;i++) fprintf(outfile,"%10d%15.6f%15.6f%15.6f\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," Mem. N1 Q1 M1 N2 Q2 M2 \n");
for(nel=0;nel<ne;nel++){
sncs(nel);
trmat();
stiff(nel);
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,"Nodal Reaction\n");
fprintf(outfile," Node RX RY RM\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]);
}
和平面刚架程序不同点在于，形成单元刚度矩阵的函数stiff()做了改变（第220行），杆系单元单元刚度矩阵可参考结构力学书籍。共分析了四种单元如下图所示。输入每个单元时还要输入单元类型typ[i]。
案例：
跨度5m，层高3m
荷载：集中力100kN，弯矩200kN，均布荷载20kN/m
截面：
输入文件：
输出文件：
经检验和结构力学求解器所得结果相近。
京公网安备 11010502036488号