平面刚架静力分析有限元程序

有限单元法

平面刚架静力分析程序

#include

#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();

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;

/************主函数**************/

void 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

for(i=0;i

for(i=0;i

for(i=0;i

}

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 loads:%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 type\n");

for(i=0;i

fprintf(outfile,"%5d%10d%10d%10d\n",i+1,jel[i][0],jel[i][1],nae[i]);

fprintf(outfile,"coordinates x and y of nodes\n");

fprintf(outfile,"node x y\n");

for(i=0;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 aera moment of interia elastic modulus\n");

for(i=0;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 r estr.-disp.-no restr.-disp.\n");

for(i=0;i

fprintf(outfile,"%10d%19.3f%19.3f\n",i+1,res[i][0],res[i][1]);

hbw();

for(i=0;i

if(npj==0) goto aa;

fprintf(outfile,"nodal loads\n");

fprintf(outfile,"node px py zm\n");

for(i=0;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

for(i=0;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

{

fscanf(infile,"%f%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

{

nel=floor(pf[np][2]+0.1);

sncs(nel-1);

fix(np);

trmat();

fis(nel-1);

fpj();

}

bb:for(i=0;i

for(j=0;j

for(nel=0;nel

{

sncs(nel);

trmat();

stiff();

matmul();

fis(nel);

addsm();

}

restr();

soleq();

outdis();

force();

}

/*===================求最大半带宽========================*/ void hbw()

{

int nel;

mbw=0;

for(nel=0;nel

{

int ma=abs(jel[nel][0]-jel[nel][1]);

if(mbw

}

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*c2;

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[2];

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

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

{ 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

{ 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

{ m=i-k+1;

c1=tsm[k-1][m-1]/tsm[k-1][0];

for(j=1;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

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

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

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

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

fprintf(outfile,"%10d%15.4f%15.4f%15.4f\n",i+1,r[3*i],r[3*i+1],r[3*i+2]);

}