八節(jié)點(diǎn)等參單元平面有限元

圖57。圖STYLEREF1\s5SEQ圖表\*ARABIC\s16應(yīng)力比較圖STYLEREF1\s5SEQ圖表\*ARABIC\s17位移比較自編程序所得結(jié)果與ANSYS分析結(jié)果進(jìn)行比較發(fā)現(xiàn)，應(yīng)力偏大而位移偏小。分析其中的原因，一方面是編程水平有限，自編程序有很多不完善之處，有很多因素沒(méi)有考慮（溫度、變形、非線性等），所得結(jié)果可信度不是很高，只能得到應(yīng)力以及位移大概的分布情況；另一方面是在后處理階段，在對(duì)高斯積分點(diǎn)結(jié)果進(jìn)行處理時(shí)，誤差難以避免，還有一方面的原因是在進(jìn)行求解時(shí)保留數(shù)據(jù)精度不夠，造成誤差較大，并且進(jìn)行數(shù)值運(yùn)算時(shí)可能會(huì)造成大數(shù)吃小數(shù)，從而引起結(jié)果的差異。

參考文獻(xiàn)(美)史密斯(Smith,I.m.)著;王崧等譯.有限單元法編程(第三版)[M].北京:電子工業(yè)出版社,2003王勖成.有限單元法[M].北京:清華大學(xué)出版社,2003宋建輝,涂志剛.MATLAB語(yǔ)言及其在有限元編程中的應(yīng)用[J].湛江師范學(xué)院學(xué)報(bào).2003.6(24):101-105鄭大玉,連宏玉,周躍發(fā).有限元編程與C語(yǔ)言[J].黑龍江商學(xué)院學(xué)報(bào).1997.3(13):23-28王偉,劉德富.有限元編程中應(yīng)用面向?qū)ο缶幊碳夹g(shù)的探討[J].三峽大學(xué)學(xué)報(bào).2001.2(23):124-128汪文立,呂士俊.二級(jí)C語(yǔ)言程序設(shè)計(jì)教程[M].北京:中國(guó)水利水電出版社,2006趙翔龍.FORTRAN90學(xué)習(xí)教程[M].北京:北京大學(xué)出版社,2002附錄程序源代碼#include<stdio.h>#include<math.h>/*Thegemotryofthemodel*/floatmesh[2];floatwide,hight;floatwsize,hsize,young,poiss,thick;inti,j,k,knelem,npoin,elem[500],ielem;floatcoord[2][1],props[3][1],lnods[8][500],ifpre[1];floatposgp[3],weigp[3],estif[136],elcod[2][8];intnpoin,nelem,kevab,igaus,jgaus,kgasp,ngash,djacb;floatgpcod[2][9],bmatx[3][16],dmatx[3][3];floatshape[8],deriv[2][8];floatxjacm[2][2],xjaci[2][2],elcod[2][8];floatcartd[2][8];floatbmatx[3][16],dmatx[3][3];floatnload[1];floatpress[3][2],pgash[2],dgash[2];structnode{intnodenu;/*Thenumberofthenode*/floatcor[2];/*Thecoordinateofthenode*/floatdisplacement[2];/*Thedisplacementofthenode*/floatload[2];/*Theloadofthenode*/floatboundary[2];}node[500];struct{intelementnu;/*Thenumberofelement*/intlocalnu[8];/*Localnumber*/intlocalcorx[8];intlocalcory[8];/*Localcoordinateofnode*/floatqx,qy,t;/*Thestressandstrain*/}element[500];voidmeshing(floatwide,floathight,floatmesh[2]){printf("Pleaseinputthemeshingsize\n");scanf("%f%f",&wsize,&hsize);getchar();mesh[0]=wide/wsize;mesh[1]=hight/hsize;}/*Theglobalcoordinateofnode*/voidcoordinate(){inti,j=1;floatx,y;for(i=1;i<=2*mesh[1]+1;i++){x=0;while(x<=wide)if(i%2!=0){node[j].cor[1]+=wsize/2;node[j].cor[2]+=(i-1)*hsize;j++;}else{node[j].cor[1]+=wsize;node[j].cor[2]+=(i-1)*hsize;j++;}}}main(){inttop[500],bottom[500];/*Thenumberoftopandbottomelement*/voidinput();voidstif();voidjacobi2();voidload();voidasem();voidsolve();voidstre();npoin=8;for(i=1;i<=8;i++)lnods[i][1]=i;for(i=1;i<=3;i++)scanf("%f",&props[i][1]);printf("TheEXis%e\nTheuois%8f\nThebtis%8f",props[1][1],props[2][1],props[3][1]);getch();printf("Pleaseinputthegemotryofthemodel\n");scanf("%f%f",&wide,&hight);getchar();printf("Thewideis%0.3f,thehightis%0.3f",wide,hight);getchar();meshing(wide,hight,mesh);printf("Thewidesizeis%f,thehightsizeis%f\n",mesh[0],mesh[1]);input();getchar();nelem=mesh[0]*mesh[1];for(i=1;i<=nelem;i++)/*Theelementnumber*/element[i].elementnu=i;for(i=1;i<mesh[0];i++)npoin+=5;for(i=1;i<mesh[1];i++)npoin+=3*mesh[0]+2;printf("%d",npoin);for(i=1;i<=npoin;i++)node[i].nodenu=i;for(i=1;i<=nelem;i++)for(j=1;j<=8;j++)element[i].localnu[j]=j;for(i=1;i<=mesh[0];i++)bottom[i]=i;j=1;for(i=mesh[0]*(mesh[1]-1)+1;i<=nelem;i++)top[j++]=i;for(i=1;i<j;i++)printf("thetopnumberis%d\n",top[i]);printf("%d\n",element[1].localnu[8]);coordinate();stif();jacobi2();load();asem();solve();stre();getchar();}/*datainput*/voidinput(){intnnode=8;intnotreadchar;/*inputelementnodenumber*/printf("inputelementnodenumber\n");for(ielem=1;ielem<=nelem;ielem++)for(i=1;i<=nnode;i++)scanf("%d",&lnods[i][ielem]);/*Inputrestrictionmessage*/printf("double-digitissymboloftherestriction,1representatesstable,2representatesgiveddisplacement\n");i=1;while(i){scanf("%d",&i);if(i!=0){scanf("%d",&j);scanf("%d",&ifpre[j]);}elsebreak;}}/*Elementstiffnessmatrix*/voidstif(){intkevab,nstre,ievab,istre,lnode,jstre,jevab;floatkgasp,exisp,etasp,dvolu;floatbtdbm,dbmat[3];voidsfr();/*Gausspoint*/posgp[1]=-sqrt(0.6);posgp[2]=0;posgp[3]=sqrt(0.6);/*weightcoefficient*/weigp[1]=5.0/9.0;weigp[2]=8.0/9.0;weigp[3]=5.0/9.0;/*numberofstress*/nstre=3;/*formelasticmatrix*/for(ielem=1;ielem<=nelem;ielem++){printf("Thenumberofelementis%d\n",ielem);for(i=1;i<=8;i++){lnode=lnods[i][ielem];for(j=1;j<=2;j++){coord[j][lnode]=node[lnode].cor[j];elcod[j][i]=coord[j][lnode];}}young=props[1][1];poiss=props[2][1];thick=props[3][1];modps(young,poiss);/*Theinitialvalueof[k]*/kevab=0;for(i=1;i<=16;i++)for(j=1;j<=i;j++){kevab++;estif[kevab]=0.0;}/*Thegausspointshapefunctionandderiv*/kgasp=0;for(igaus=1;igaus<=3;igaus++)for(jgaus=1;jgaus<=3;jgaus++){kgasp=kgasp+1;exisp=posgp[igaus];etasp=posgp[jgaus];printf("Thenumberofgausspointis%d\n",kgasp);sfr2(exisp,etasp);jacob2(ielem,djacb,kgasp,elcod);dvolu=djacb*weigp[igaus]*weigp[jgaus]*thick;bmatps()/*Theinitialvalueofelasticmatrix[D]*/kevab=0;for(ievab=1;ievab<=16;ievab++){for(istre=1;istre<=nstre;istre++){dbmat[istre]=0.0;for(jstre=1;jstre<=nstre;jstre++)dbmat[istre]=dbmat[istre]+bmatx[jstre][ievab]*dmatx[jstre][istre];}for(jevab=1;jevab<=ievab;ievab++){kevab=kevab+1;btdbm=0.0;for(istre=1;istre<=nstre;istre++)btdbm=btdbm+dbmat[istre]*bmatx[istre][jevab];estif[kevab]=estif[kevab]+btdbm*dvolu;}}}}}/*floatsfr2(floats,floatt){floats2,t2,ss,tt,st,stt,sst,st2;/*Shapefunction*//*thesevariablesaretosimplifytheformula*/s2=s*2.0;t2=t*2.0;ss=s*s;tt=t*t;st=s*t;stt=s*t*t;sst=s*s*t;st2=s*t*2.0;/*shapefunction*/shape[1]=(-1.0+st+ss+tt-sst-stt)/4.0;shape[2]=(1.0-t-ss+sst)/2.0;shape[3]=(-1.0-st+ss+tt-sst+stt)/4.0;shape[4]=(1.0+s-tt-stt)/2.0;shape[5]=(-1.0+st+ss+tt+sst+stt)/4.0;shape[6]=(1.0+t-ss-sst)/2.0;shape[7]=(-1.0-st+ss+tt+sst-stt)/4.0;shape[8]=(1.0-s-tt+stt)/2.0;/*differentialcoefficienttolocalcoordinate*/deriv[1][1]=(t+s2-st2-tt)/4.0;deriv[1][2]=-s+st;deriv[1][3]=(-t+s2-st2+tt)/4.0;deriv[1][4]=(1.0-tt)/2.0;deriv[1][5]=(t+s2+st2+tt)/4.0;deriv[1][6]=-s-st;deriv[1][7]=(-t+s2+st2-tt)/4.0;deriv[1][8]=(-1.0+tt)/2.0;deriv[2][1]=(s+t2-ss-st2)/4.0;deriv[2][2]=(-1.0+ss)/2.0;deriv[2][3]=(-s+t2-ss+st2)/4.0;deriv[2][4]=-t-st;deriv[2][5]=(s+t2+ss+st2)/4.0;deriv[2][6]=(1.0-ss)/2.0;deriv[2][7]=(-s+t2+ss-st2)/4.0;deriv[2][8]=-t+st;}*//*Jacobimatrix*/voidjacobi2(){/*xjacm[2][2]:jacobimatrix;xjaci[2][2]:[j]-1;elcod[2][8]:localcoordinate*/floatcartd[2][8],gpcod[2][9];inti,j,k;for(i=1;i<=2;i++){gpcod[i][kgasp]=0.0;for(j=1;j<=8;j++)gpcod[i][kgasp]=gpcod[i][kgasp]+elcod[i][j]*shape[j];}for(i=1;i<=2;i++)for(j=1;j<=2;j++){xjacm[i][j]=0.0;for(k=1;k<=8;k++)xjacm[i][j]=xjacm[i][j]+deriv[i][k]*elcod[j][k];}/*CaculatethevalueofJacobimatrix*/djacb=xjacm[1][1]*xjacm[2][2]-xjacm[1][2]*xjacm[2][1];printf("Thedetofjacobimatrixis%f\n",djacb);if(djacb<=1e-6)printf("Thejacobidetof%dislessorequal0\n",ielem);/*Theelementof[j]-1*/xjaci[1][1]=xjacm[2][2]/djacb;xjaci[2][2]=xjacm[1][1]/djacb;xjaci[1][2]=-xjacm[1][2]/djacb;xjaci[2][1]=-xjacm[2][1]/djacb;/*Thederivofshapefunctiontoglobalcoordinate*/for(i=1;i<=2;i++)for(k=1;k<=8;k++){cartd[i][k]=0.0;for(j=1;j<=2;j++)cartd[i][k]=cartd[i][k]+xjaci[i][j]*deriv[j][k];}}/*Formelasticmatris*//*floatmodps(){floatbmatx[3][16],dmatx[3][3];floatconstant;inti,j,y,p;y=props[1][1];p=props[2][1];for(i=1;i<=3;i++)for(j=1;j<=3;j++)dmatx[i][j]=0.0;/*IfNtype=1,itisplanestressquestion*/constant=y/(1.0-p*p);dmatx[1][1]=constant;dmatx[2][2]=constant;dmatx[1][2]=constant*p;dmatx[2][1]=constant*p;dmatx[3][3]=constant*(1.0-p)/2.0;}*//*Formstrainmatrix*//*voidbmatps(){floatcartd[2][8],shape[8],deriv[2][8],bmatx[3][16],dmatx[3][3];intnpoin,nelem,ngash,mgash;floatgpcod[2][9];/*Theinitialvalueof[B]*/for(i=1;i<=16;i++)for(j=1;j<=3;j++)bmatx[j][i]=0.0;/*Form[B]*/ngash=0;for(i=1;i<=8;i++){mgash=ngash+1;ngash=mgash+1;bmatx[1][mgash]=cartd[1][i];bmatx[1][ngash]=0.0;bmatx[2][mgash]=0.0;bmatx[2][ngash]=cartd[2][i];bmatx[3][mgash]=cartd[2][i];bmatx[3][ngash]=cartd[1][i];}}*//*equalloadofnode*/voidload(){floatnload[500];floatelcod[2][3],eload[16],posgp[3],weigp[3];floatpress[3][2],pgash[2],dgash[2],noprs[3];intiedge,nedge,kount;floatsgmar,tau,s,dvolu,pxcom,pycom,n,m,t;intlnode,number,nloca;/*Thenumberofloadside*//*elementload*//*Thepositionofgausspoint*/printf("inputthenumberofloadside*/n");scanf("%d",&number);posgp[1]=-sqrt(0.6);posgp[2]=0.0;posgp[3]=sqrt(0.6);/*Theweightofgausspoint*/weigp[1]=5.0/9.0;weigp[2]=8.0/9.0;weigp[3]=5.0/9.0;/*Theinitialloadofelement*/for(i=1;i<=nelem;i++)nload[i]=0.0;for(i=1;i<=number;i++)/*Theinitialvalueofequalnodeload*/for(i=1;i<=16;i++)eload[i]=0;scanf("%d",&nedge);for(iedge=1;iedge<=nedge;iedge++){for(i=1;i<=3;i++)scanf("%f%f",&sgmar,&tau);/*ifsgmar=0,inputnodeloadq(1,2,3)andt(1,2,3)*/if((fabs(sgmar)<1e-8)&&(fabs(tau)<1e-8))for(j=1;j<=3;j++)for(i=1;i<=2;i++)scanf("%f",&press[j][i]);elsefor(i=1;i<=3;i++){press[i][1]=sgmar;press[i][2]=tau;}/*Thecoordinateoftheloadnode*/for(i=1;i<=3;i++){/*inputloadnode*/printf("inputnodeloadnumber\n");scanf("%d",&noprs[i]);lnode=noprs[i];for(j=1;j<=2;j++)coord[j][lnode]=node[lnode].cor[j];elcod[j][i]=coord[j][lnode];}t=-1.0;for(igaus=1;igaus<=3;igaus++){s=posgp[igaus];sfr(s,t);for(j=1;j<=2;j++){pgash[j]=0.0;dgash[j]=0.0;/*currentpointq,tandthevalueofderiv*/{pgash[j]=pgash[j]+press[i][j]*shape[i];dgash[j]=dgash[j]+elcod[j][i]*deriv[1][i];}}dvolu=weigp[igaus];pxcom=dgash[1]*pgash[2]-dgash[2]*pgash[1];pycom=dgash[1]*pgash[1]+dgash[2]*pgash[2];for(i=1;i<=8;i++){nloca=lnods[i][ielem];if(nloca==noprs[1]){j=i+2;break;}}j=i+2;kount=0;for(k=i;k<=j;k++){kount=kount+1;n=(k-1)*2+1;m=n+1;/*ifthelocalnodenumner>8,thenitis1*/if(k>8){n=1;m=2;}/*sumtogettheequalnodeload*/eload[n]=eload[n]+shape[kount]*pxcom*dvolu;eload[m]=eload[m]+shape[kount]*pycom*dvolu;}}}print("Theelementloadmatrxis/n");for(i=1;i<=16;i++)printf("%f\n",eload[i]);}/*Toformglobalstiffnessmatrixandloadmatrix*/voidasem(){floatv[1],a[1];floatlnods[8][1],ifpre[1],nload[1],maxa[1];floatncodf[2],estif[136],eload[16];floatdlarge,x;intipoin,kmin,lnode,khigh,kdofn,nn,nnm,nwk,iwk,in;intinodi,inode,icodf,idofn,ievab,icoln,jnode,lnodj,jdofn,jevab,jrow,ldis,lnodi;maxa[1]=1;for(ipoin=1;ipoin<=npoin;ipoin++){kmin=npoin+1;for(ielem=1;ielem<=nelem;ielem++)for(i=1;i<=8;i++){lnode=lnods[i][ielem];if(lnode==ipoin)break;}for(i=1;i<=8;i++){lnode=lnods[i][ielem];if(lnode>=kmin)break;}kmin=lnode;khigh=(ipoin-1)*2+j+1;for(j=1;j<=2;j++){kdofn=(ipoin-1)*2+j+1;maxa[kdofn]=maxa[kdofn-1]+khigh+j;}}nn=npoin*2;nnm=nn+1;nwk=maxa[nnm]-1;printf("Thesumis%d\n",nwk);for(iwk=1;iwk<=nwk;iwk++)a[iwk]=0.0;for(in=1;in<=nn;nn++)v[in]=0.0;dlarge=1.0e+15;for(ielem=1;ielem<=nelem;ielem++){printf("Thenumberofelementis%d\n",ielem);for(inode=1;inode<=8;inode++){inodi=lnods[inode][ielem];icodf=ifpre[lnodi];ncodf[1]=icodf/10;ncodf[2]=icodf-ncodf[1]*10;for(idofn=1;idofn<=2;idofn++){ievab=(inode-1)*2+idofn;icoln=(lnodi-1)*2+idofn;for(jnode=1;jnode<=8;jnode++){lnodj=lnods[jnode][ielem];for(jdofn=1;jdofn<=2;jdofn++){jevab=(jnode-1)*2+jdofn;jrow=(lnodj-1)*2+jnode;if(jrow>icoln)break;if(jevab>ievab)kevab=jevab*(jevab-1)/2+ievab;elsekevab=ievab*(ievab-1)/2+jevab;ldis=maxa[icoln]+(icoln-jrow);a[ldis]=a[ldis]+estif[kevab];}}if(nload[ielem]>0)v[icoln]=v[icoln]+eload[ievab];if(ncodf[idofn]==0)break;kevab=ievab*(ievab+1)/2;ldis=maxa[icoln];a[ldis]=a[ldis]+dlarge*estif[kevab];}}}printf("Entergiveddisplacement\n");for(i=1;i<=npoin;i++){icodf=ifpre[i];ncodf[1]=icodf/10;ncodf[2]=icodf-ncodf[1]*10;for(j=1;j<=2;j++){if(ncodf[j]==2)scanf("%f",&x);icoln=(i-1)*2+j;v[icoln]=a[ldis]*x;}}}/*Tosolvetheequation*/voidsolve(){floatv[1],a[1],maxa[1];intnpoin,nelem,ntype,nmats;intl,k,nn,nnm,nwk,n,ku,kl,kh,kn;floatic,b,c,klt,ki,nd,kk,m1,m2,kul;nn=npoin*2;nnm=nn+1;for(n=1;n<=nn;n++){kn=maxa[n];kl=kn+1;ku=maxa[n+1]-1;kh=ku-kl;if(kh<0)if(a[kn]<=0){printf("*****Stoprun!Coefficientmatrixisnotillegal!*****Non+positivepivotforequation%d,pivot=%f\n",n,a[kn]);break;}elseif(kh==0){k=n;b=0.0;for(kk=kl;kk<=ku;kk++){k=k-1;ki=maxa[k];c=a[kk]/a[ki];b=b+c*a[kk];a[kk]=c;a[kn]=a[kn]-b;}if(kh>0){k=n-kh;ic=0;klt=ku;for(j=1;j<=kh;j++){ic=ic+1;klt=klt-1;ki=maxa[k];nd=maxa[k+1]-ki-1;if(nd<=0)k=k+1;c=0;for(l=1;l<=kk;l++){m1=ki+l;m2=klt+l;c=c+a[m1]*a[m2];}}printf("Reduceright--hand--sideloadvector\n");for(n=1;n<=nn;n++){kl=maxa[n]+1;ku=maxa[n+1]-1;kul=ku-kl;if(kul>=0){k=n;c=0.0;for(kk=kl;kk<=ku;kk++){k=k-1;c=c+a[kk]*v[k];v[n]=v[n]-c;}}}printf("Back--substitute\n");for(n=1;n<=nn;n++){k=maxa[n];v[n]=v[n]/a[k];}if(nn!=1){n=nn;for(l=2;l<=nn;l++){kl=maxa[n]+1;ku=maxa[n+1]-1;kul=ku-kl;if(kul>=0){k=n;for(kk=kl;kk<=ku;kk++){k=k-1;v[k]=v[k]-a[kk]*v[n];}}n=n-1;}printf("Outputdisplacement\n");printf("xdirection,ydirection\n");m2=0;for(i=1;i<=npoin;i++){m1=m2+1;m2=m1+1;printf("Thenumberofgausspointis%d,%16.5e%16.5e\n",i,v[m1],v[m2]);}}}}}/*Togetelementstress*//*voidstre()*/{floatv[1],lnods[8][500];floatposgp[3],elcod[2][8],eldis[16],be[3],s[3],sp[2];intievab,mdofn,ldofn,l,lnode;flo