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

有限单元法平面刚架静力分析程序#include#includestring.h#includestdlib.h#includestdio.h#includemath.hvoid 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 sm66,tsm30030,res602,pj300,t66,d66,r300, fo6,foj6,pf2004,x100,y100,ae103,sl,sn,cs,eal,eil;int nj,ne,nt,nr,npj,npf,nn,mbw,jel1002,nae100,is6;FILE *infile,*outfile;/*主函数*/void main()char name130,name230;int i,j,nel,np;printf(please enter data-filenamen);scanf(%s,name1);printf(please enter result-filenamen);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;ine;i+) fscanf(infile,%d%d%d,&jeli0,&jeli1,&naei);for(i=0;inj;i+) fscanf(infile,%f%f,&xi,&yi);for(i=0;int;i+) fscanf(infile,%f%f%f,&aei0,&aei1,&aei2);for(i=0;inr;i+) fscanf(infile,%f%f,&resi0,&resi1);elseprintf(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 framen);fprintf(outfile,input datan);fprintf(outfile,*n); fprintf(outfile,control datan); fprintf(outfile,the num.of node:%3dn,nj); fprintf(outfile,the num.of mem:%3dn,ne); fprintf(outfile,the num.of type of section characteristic:%3dn,nt); fprintf(outfile,the num.of restricted degrees of freedom:%3dn,nr); fprintf(outfile,the num.of nodal loads:%3dn,npj); fprintf(outfile,the num.of non-nodal loads:%3dn,npf);fprintf(outfile,the num.of nodal degrees of freedom:%3dn,nn);fprintf(outfile,information of mem.n);fprintf(outfile,mem. start node end node typen); for(i=0;ine;i+) fprintf(outfile,%5d%10d%10d%10dn,i+1,jeli0,jeli1,naei); fprintf(outfile,coordinates x and y of nodesn); fprintf(outfile,node x yn); for(i=0;inj;i+)/*第二次输入*/ fprintf(outfile,%10d%10.2f%10.2fn,i+1,xi,yi); fprintf(outfile,information of cross section each mem.n); fprintf(outfile,typeaeramoment of interia elastic modulusn); for(i=0;int;i+) fprintf(outfile,%8d%15.5f%15.5f%20.2fn,i+1,aei0,aei1,aei2); fprintf(outfile,information restrictionn); fprintf(outfile,rester.-norestr.-disp.-norestr.-disp.n);for(i=0;inr;i+) fprintf(outfile,%10d%19.3f%19.3fn,i+1,resi0,resi1); hbw();for(i=0;inn;i+) pji=0.0; if(npj=0) goto aa; fprintf(outfile,nodal loadsn); fprintf(outfile,nodepxpyzmn);for(i=0;inpj;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.2fn,no,px,py,zm); pj3*no-3=px; pj3*no-2=py; pj3*no-1=zm; aa:for(i=0;inn;i+) ri=0.0; for(i=0;inr;i+) int ni; ni=floor(resi0+0.1)-1; if(pjni!=0.0)rni=rni-pjni; if(npf=0) goto bb; fprintf(outfile,non-nodal loadsn); fprintf(outfile,loads leng.supp.to load mem typen); for(i=0;inpf;i+) fscanf(infile,%f%f%f%f%f,&pfi0,&pfi1,&pfi2,&pfi3); fprintf(outfile,%15.2f%15.2f%15.2f%15.2fn,pfi0,pfi1,pfi2,pfi3); for(np=0;npnpf;np+) nel=floor(pfnp2+0.1); sncs(nel-1); fix(np); trmat(); fis(nel-1); fpj(); bb:for(i=0;inn;i+) for(j=0;jmbw;j+) tsmij=0.0; for(nel=0;nelne;nel+) sncs(nel); trmat(); stiff(); matmul(); fis(nel); addsm(); restr(); soleq(); outdis(); force(); /*=求最大半带宽=*/void hbw()int nel;mbw=0;for(nel=0;nelne;nel+)int ma=abs(jelnel0-jelnel1);if(mbwma) mbw=ma;mbw=3*(mbw+1);fprintf(outfile,Half_Bandwidth Mbw:%5dn,mbw);/*矩阵相乘*/void matmul()int i,j,k;for(i=0;i6;i+)for(j=0;j6;j+) dij=0.0;for(i=0;i6;i+)for(j=0;j6;j+)for(k=0;k6;k+) dij=dij+tki*smkj;for(i=0;i6;i+)for(j=0;j6;j+) smij=0.0;for(i=0;i6;i+)for(j=0;j6;j+)for(k=0;k6;k+)smij=smij+dik*tkj;/*求单元常数*/void sncs(int nel)int ii,jj,k; float xi,yi,xj,yj;ii=jelnel0;jj=jelnel1;xi=xii-1;xj=xjj-1;yi=yii-1;yj=yjj-1;sl=sqrt(xj-xi)*(xj-xi)+(yj-yi)*(yj-yi);sn=(yj-yi)/sl;cs=(xj-xi)/sl;k=naenel;eal=aek-10*aek-12/sl;eil=aek-11*aek-12/sl;/*求单元坐标转换矩阵*/void trmat()int i,j;for(i=0;i6;i+)for(j=0;j6;j+)tij=0.0;t00=cs;t01=sn; t10=-sn;t11=cs;t22=1.0;for(i=0;i3;i+)for(j=0;j3;j+) ti+3j+3=tij;/*求单元固端力*/void fix(int np)float w,c,c1,c2,cc,cc1,cc2;int i,im;w=pfnp0;c=pfnp1;c1=c/sl; c2=c1*c1; cc=sl-c;cc1=cc/sl; cc2=cc1*cc1;for(i=0;i6;i+) foi=0.0;im=floor(pfnp3+0.1); if(im=1) fo1=-w*c*(1.0-c2+c2*c1/2.0);fo2=-w*c*c*(6.0-8.0*c1+3*c2)/12.0;fo4=-w*c-fo1;fo5=w*c1*c*c*(4.0-3.0*c1)/12.0; else if(im=2)fo1=-w*(1+2*c1)*cc2;fo2=-w*c*c2;fo4=-w*(1+2*cc1)*c2;fo5=w*cc*c2;else if(im=3)fo1=6*w*c1*cc1/sl;fo2=w*cc1*(2-3*cc1);fo4=-fo2;fo5=w*c1*(2-3*c1);else if(im=4)fo1=0.25*w*c*(2.-3.*c2+1.6*c2*c1);fo2=-w*c*c*(2.-3.*c1+1.2*c2)/6.0;fo4=-w*c/2.0-fo2;fo5=0.25*w*c*c*c1*(1.-0.8*c1);else if(im=5)fo0=-w*cc1;fo3=-w*c1;else if(im=6)fo0=-w*c*(1.-0.5*c1);fo3=-0.5*w*c*c1;/*形成总结点荷载向量*/void fpj() int i,j,k,ii,jj; for(k=0;k6;k+) fojk=0.0; for(i=0;i6;i+) for(j=0;j6;j+) foji=foji+tji*foj; for(ii=0;ii6;ii+) pjisii=pjisii-fojii; /*形成单元刚度矩阵*/void stiff() int i,j; float s1,s2; for(i=0;i6;i+) for(j=0;j6;j+) smij=0.0; sm00=eal; sm30=-eal; sm33=eal; s1=12.*eil/(sl*sl);s2=6.*eil/sl; sm11=s1; sm41=-s1; sm44=s1; sm21=s2; sm51=s2; sm42=-s2; sm54=-s2; sm22=4*eil; sm55=4*eil; sm52=2*eil; for(i=0;i6;i+) for(j=0;ji;j+) smji=smij;/*由单元位移分量码L形成总刚位移分量码IS（L）*/void fis(int nel) int i,j; for(i=0;i2;i+) for(j=0;j3;j+) is3*i+j=3*(jelneli-1)+j; /*形成结构原始总刚度矩阵*/void addsm() int i,j; for(i=0;i6;i+) for(j=0;j=0) tsmisikc=tsmisikc+smij;/*约束处理*/void restr() int i; for(i=0;inr;i+) int ni; ni=floor(resi0+0.1); tsmni-10=1.0e25; pjni-1=tsmni-10*resi1; /*解线性方程组*/void soleq() float c1; int k,ni,im,i,m,j,nm,jm; for(k=1;knn;k+) if(fabs(tsmk-10)nn) im=nn; for(i=k+1;iim+1;i+) m=i-k+1; c1=tsmk-1m-1/tsmk-10; for(j=1;jmbw-m+2;j+) tsmi-1j-1=tsmi-1j-1-c1*tsmk-1j+i-k-1; pji-1=pji-1-c1*pjk-1; if(fabs(tsmnn-10)0;i-) nm=nn-i+1; jm=nm; if(nmmbw)jm=mbw; for(j=2;jjm+1;j+) pji-1=pji-1-tsmi-1j-1*pjj+i-2; pji-1=pji-1/tsmi-10; fin:return; /*输出位移向量*/ void outdis() int i; fprintf(outfile,output resultsn); fprintf(outfile,*n); fprintf(outfile,nodal displacementsn); fprintf(outfile, node u v cetan); for(i=0;inj;i+) fprintf(outfile,%10d%15.6f%15.6f%15.6fn,i+1,pj3*i,pj3*i+1,pj3*i+2);/*求单元杆端力、支座反力或结合点力*/void force () float dj6,f6,fj6,dd6; int nel,np,i,j,ip; fprintf(outfile,Mem.N1 Q1 M1 N2 Q2 M2 n); for(n