彈性力學(xué)平面問題有限元程序

有限單元法彈性力學(xué)平面問題有限元程序#defineNODE_NUM300#defineELE_NUM200#defineNODE2NODE_NUM*2#include<stdlib.h>#include<stdio.h>#include<math.h>intnj,ne,nm,nz,npj,jm[ELE_NUM][4],nzc[60],nj2,bw,lxm;floatea,e,um,oul,te,xy[NODE_NUM][4],EA[5],p[NODE2],pj,F[ELE_NUM],etm[5][5];floatek[7][7];double**K,*f;intBAND2(intn,intm,double**K,double*f);double**TwoArrayAlloc(int,int);voidTwoArrayFree(double**);FILE*outf;/*-------------------------------------------*/double**TwoArrayAlloc(intr,intc){double*x,**y;intn;x=(double*)calloc(r*c,sizeof(double));y=(double**)calloc(r,sizeof(double*));for(n=0;n<=r-1;++n)y[n]=&x[c*n];return(y);}/*-------------------------------------------*/voidTwoArrayFree(double**x){free(x[0]);free(x);}/*------------------------------------------*/intBAND2(intn,intm,double**K,double*f){inti,j,t,ij,ji,it,jt,tm,m1;doubles,w;m1=m+1;for(i=1;i<=n;i++){if(K[i-1][m1-1]<=0)return(0);w=0.0;if(i>m1)tm=i-m;elsetm=1;for(j=tm;j<=i;j++){s=0.0;ij=j-i+m1;for(t=tm;t<=j-1;t++){it=t-i+m1;jt=t-j+m1;s=s+K[i-1][it-1]*K[j-1][jt-1]/K[t-1][m1-1];}K[i-1][ij-1]=K[i-1][ij-1]-s;if(j==i)f[i-1]=f[i-1]-w;elsew=w+K[i-1][ij-1]*f[j-1]/K[j-1][m1-1];}}for(i=n;i>=1;i--){s=0.0;if(i>n-m1)tm=n;elsetm=i+m;for(j=i+1;j<=tm;j++){ji=i-j+m1;s=s+K[j-1][ji-1]*f[j-1];}f[i-1]=(f[i-1]-s)/K[i-1][m1-1];}return1;}/*-----------------------------------------*/voidinput(){intjj,j,i,nh,nl;floatdx,dy;FILE*infile;charname1[30],name2[30];R:printf("pleaseenterdata-filename\n");scanf("%s",name1);if((infile=fopen(name1,"r"))==NULL){printf("thedata-filenotexit!");gotoR;}printf("Pleaseenterresult-filename\n");scanf("%s",name2);outf=fopen(name2,"w");fscanf(infile,"%d%d%d%d%d%d5d",&nj,&ne,&nm,&nz,&npj,&lxm);for(i=1;i<=nj;i++){for(j=1;j<3;j++)fscanf(infile,"%f",&xy[i][j]);}for(i=1;i<=ne;i++){for(j=0;j<4;j++)fscanf(infile,"%d",&jm[i][j]);}for(i=1;i<=nm;i++){for(j=1;j<=4;j++)fscanf(infile,"%f",&etm[i][j]);}for(i=1;i<=nz;i++)fscanf(infile,"%d",&nzc[i]);nj2=nj*2;for(i=0;i<nj2;i++)p[i]=0.0;for(i=0;i<npj;i++){fscanf(infile,"%d%f",&jj,&pj);p[jj]=pj;}fprintf(outf,"TheNum.OfNodes:%3d\n",nj);fprintf(outf,"TheNum.OfElements.:%3d\n",ne);fprintf(outf,"TheNum.OfTypeOfSectionCharacteristic:%3d\n",nm);fprintf(outf,"TheNum.OfRestriction:%3d\n",nz);fprintf(outf,"TheNum.OfNodalLoads:%3d\n",npj);fprintf(outf,"LXM=%d\n",lxm);fprintf(outf,"CoordinatesxandyOfNodes:\n");fprintf(outf,"Nodexy\n");for(i=1;i<=nj;i++){fprintf(outf,"%10d%10.2f%10.2f\n",i,xy[i][1],xy[i][2]);}fprintf(outf,"TheNodesNum.OfMem:\n");fprintf(outf,"Mem.TypeLeftRight\n");for(i=1;i<=ne;i++){fprintf(outf,"%10d%10d%10d%10d%10d\n",i,jm[i][0],jm[i][1],jm[i][2],jm[i][3]);}for(i=1;i<=nm;i++){fprintf(outf,"nm=%dE=%f,UM=%f,OUL=%f,T=%f\n",i,etm[i][1],etm[i][2],etm[i][3],etm[i][4]);}fprintf(outf,"ThePositionOfRetriction:");for(i=1;i<=nz;i++){fprintf(outf,"%d",nzc[i]);}fprintf(outf,"\n===========loadsofnodes==============\n");for(i=1;i<=nj2;i++){fprintf(outf,"%12d%12.2f\n",i,p[i]);}if(lxm>0){for(i=1;i<=nm;i++) { etm[i][1]=etm[i][1]/(1-etm[i][2]*etm[i][2]); }}fclose(infile);}/*------------計(jì)算半帶寬函數(shù)----------------*/intbwidth(){intie,i,j,min,max,jj,ib,lk,m,m1,m2,t,nn[8];ib=0;for(ie=1;ie<=ne;ie++){ m=abs(jm[ie][1]-jm[ie][2]); m1=abs(jm[ie][1]-jm[ie][3]); m2=abs(jm[ie][2]-jm[ie][3]); if(m>ib)ib=m; if(m1>ib)ib=m1;if(m2>ib)ib=m2;}bw=2*ib+2;fprintf(outf,"bw=%d\n",bw);returnbw;}/*--------------形成單元?jiǎng)偠染仃嚭瘮?shù)-----------*/voidstiff(intie){inti,j,k,m,n,i1,j1,i2,j2,ii,jj;floatcc,bb,bc,cb,a1,u1,c[4],b[4];i=jm[ie][1];j=jm[ie][2];m=jm[ie][3];/*取出單元節(jié)點(diǎn)號*/n=jm[ie][0];e=etm[n][1];um=etm[n][2];oul=etm[n][3];te=etm[n][4];/*取出材料參數(shù)*/c[3]=xy[j][1]-xy[i][1];b[3]=xy[i][2]-xy[j][2];c[2]=xy[i][1]-xy[m][1];b[2]=xy[m][2]-xy[i][2];c[1]=-c[2]-c[3];b[1]=-b[2]-b[3];/*計(jì)算[B]陣中的bi,ci,bj,cj,bm,cm*/ea=(b[2]*c[3]-b[3]*c[2])/2.0;/*計(jì)算單元面積*/a1=e*te/4.0/(ea-um*um*ea);u1=(1.0-um)/2.0;for(ii=1;ii<=3;ii++)/*按2子塊循環(huán)求出單剛計(jì)算單元面積*/{for(jj=1;jj<=3;jj++){bb=b[ii]*b[jj];cc=c[ii]*c[jj];cb=c[ii]*b[jj];bc=b[ii]*c[jj];i2=2*ii;j2=2*jj;i1=i2-1;j1=j2-1;ek[i1][j1]=a1*(bb+u1*cc);ek[i1][j2]=a1*(um*bc+u1*cb);ek[i2][j1]=a1*(um*cb+u1*bc);ek[i2][j2]=a1*(cc+u1*bb);}}}/*-------集成總剛函數(shù)，存在剛下三角--------------------------*/intekzk(intie){inti1,j1,i,j,i2,j2,ii,jj,ji;for(i1=1;i1<=3;i1++){for(i2=1;i2<3;i2++){i=2*(i1-1)+i2;ii=2*(jm[ie][i1]-1)+i2;for(j1=1;j1<=3;j1++){for(j2=1;j2<3;j2++){j=2*(j1-1)+j2;jj=2*(jm[ie][j1]-1)+j2;ji=bw+jj-ii;if(ji<=bw)K[ii-1][ji-1]=K[ii-1][ji-1]+ek[i][j];}}}}}/*--------用0，1置換法進(jìn)行約束處理-----------------------------*/voidrestrict(){inti,j,i1,j1,ii;for(i=1;i<=nz;i++){j=nzc[i]-1;for(i1=0;i1<bw-1;i1++){K[j][i1]=00.00;ii=j+i1+1;if(ii<nj2)K[ii][bw-i1-2]=00.00;}K[j][bw-1]=1.00;f[j]=0;}}/*----------計(jì)算單元應(yīng)力函數(shù)-----------------*/voidforce(){inti,j,ie,m,n,j1,j2,ii,jj;floataa,pyl,ryl,w[7],c[4],b[4],s1[4][7],yl[4],yl1[4];fprintf(outf,"==========stressesofelements==========\n");fprintf(outf,"NOSXSYSTS1S2CETA\n\n");for(ie=1;ie<=ne;ie++){i=jm[ie][1];j=jm[ie][2];m=jm[ie][3];n=jm[ie][0];e=etm[n][1];um=etm[n][2];oul=etm[n][3];te=etm[n][4];c[3]=xy[j][1]-xy[i][1];b[3]=xy[i][2]-xy[j][2];c[2]=xy[i][1]-xy[m][1];b[2]=xy[m][2]-xy[i][2];c[1]=-c[2]-c[3];b[1]=-b[2]-b[3];ea=(b[2]*c[3]-b[3]*c[2])/2.0;aa=e/(2.0*ea*(1-um*um));w[1]=f[2*i-2];w[2]=f[2*i-1];w[3]=f[2*j-2];w[4]=f[2*j-1];w[5]=f[2*m-2];w[6]=f[2*m-1];for(ii=1;ii<=3;ii++){j1=2*ii-1;j2=j1+1;s1[1][j1]=b[ii]*aa;s1[2][j1]=s1[1][j1]*um;s1[3][j2]=s1[1][j1]*(1-um)/2.0;s1[2][j2]=c[ii]*aa;s1[1][j2]=s1[2][j2]*um;s1[3][j1]=s1[2][j2]*(1-um)/2.0;}for(ii=1;ii<=3;ii++){yl[ii]=0.0;for(jj=1;jj<=6;jj++){yl[ii]=yl[ii]+s1[ii][jj]*w[jj];}}pyl=(yl[1]+yl[2])/2.0;ryl=sqrt(((yl[1]-yl[2])/2.0)*((yl[1]-yl[2])/2.0)+yl[3]*yl[3]);yl1[1]=pyl+ryl;yl1[2]=pyl-ryl;aa=abs(yl[2]-yl1[2]);if(abs(yl[2]-yl1[2])<6.1e-2){yl1[3]=0.0;}elseyl1[3]=90.0-57.29578*atan(yl[3]/(yl[2]-yl1[2]));fprintf(outf,"%2d%10.2f%10.2f%10.2f%10.2f%10.2f%10.2f\n",ie,yl[1],yl[2],yl[3],yl1[1],yl1[2],yl1[3]);}}/*----------------------------------------------