有限元上機(jī)報(bào)告

彈性力學(xué)及有限元基礎(chǔ)上機(jī)實(shí)踐報(bào)告PAGEPAGE1 在以常應(yīng)變?nèi)切螁卧?jì)算彈性力學(xué)平面問題原理的基礎(chǔ)上，編寫可計(jì)算平面應(yīng)力問題，亦可計(jì)算平面應(yīng)變問題的通用程序。給出程序設(shè)計(jì)的流程圖，程序源代碼和程序應(yīng)用的工程算例及其結(jié)果分析。一、程序流程圖此程序采用本人比較熟悉C語言編寫。C語言相比于FORTRAN語言擁有更強(qiáng)大的實(shí)用編程環(huán)境，同時(shí)由于提供了指針變量、動(dòng)態(tài)內(nèi)存分配函數(shù)和結(jié)構(gòu)體變量等，使其編寫的軟件在維護(hù)性、可讀性和內(nèi)存利用效率方面具有明顯的優(yōu)勢(shì)。程序總框圖圖1程序總框圖 整個(gè)程序由一個(gè)主調(diào)主程序(主函數(shù)main（）)和七個(gè)子程序組成，其中數(shù)據(jù)子程序DATA()用來接受輸入的參數(shù)和變量。單剛子程序ELEST()為含有參數(shù)的函數(shù)，功能控制參數(shù)iask可取1、2、3，分別計(jì)算相應(yīng)單元的面積、應(yīng)力矩陣和單元?jiǎng)偠染仃?，主程序不直接調(diào)用它，而通過總剛子程序、載荷子程序和求應(yīng)力子程序間接調(diào)用?？倓傋映绦騎OTSTI()用來合成總剛矩陣，載荷子程序LOAD()用來計(jì)算合成載荷，支承子程序SUPPOR()用來引入約束，解方程子程序SOLVEQ()用來求解并輸出各個(gè)結(jié)點(diǎn)的位移，求應(yīng)力子程序STRESS()用來計(jì)算和輸出應(yīng)力、主應(yīng)力及主平面角。2.主函數(shù)及各子程序流程圖圖2voidmain()和voidDATA()圖3voidELEST(intmeo,intiask)圖4voidTOTSTI() 圖5voidLOAD()圖6voidSUPPOR()圖7voidSOLVEQ()圖8voidSTRESS()二、程序源代碼/*有限元計(jì)算程序*/#include<stdio.h>#include<math.h>intnj,ne,nz,ndd,ind,nj2;intjm[100][3],nzc[200],npj1,npj;floateo,un,gama,te,ae;floatcjz[100][2],pj[100][2];floatb[3][6],d[3][3],s[3][6],eke[6][6],tkz[200][20],p[200];voidDATA()/*數(shù)據(jù)輸入函數(shù)*/{inti,j;printf("Qingshuru6gecanshu\n");scanf("%d,%d,%d,%d,%d,%d",&nj,&ne,&nz,&ndd,&npj,&ind);nj2=nj*2;npj1=npj+1;getchar();printf("Qingshuru4gejisuancanshu\n");scanf("%f,%f,%f,%f",&eo,&un,&gama,&te);getchar();printf("nj=%dne=%dnz=%dndd=%dnpj=%dind=%d\n",nj,ne,nz,ndd,npj,ind);printf("eo=%fun=%fgama=%fte=%f\n",eo,un,gama,te);getchar();printf("QingshuruJMjuzheng\n");for(i=0;i<ne;i++){for(j=0;j<3;j++)scanf("%d",&jm[i][j]);getchar();}printf("QingshuruCJZjuzheng\n");for(i=0;i<nj;i++){for(j=0;j<2;j++)scanf("%f",&cjz[i][j]);getchar();}printf("QingshuruNZCjuzheng\n");for(i=0;i<nz;i++)scanf("%d",&nzc[i]);getchar();printf("QingshuruPJjuzheng\n");for(i=0;i<npj1;i++){for(j=0;j<2;j++)scanf("%f",&pj[i][j]);getchar();}getchar();printf("CJZjuzhengruxia\n");for(i=0;i<nj;i++){for(j=0;j<2;j++){printf("%f",cjz[i][j]);if(j%2==1)printf("\n");}}getchar();}voidELEST(intmeo,intiask)/*單元?jiǎng)偠染仃?/{intie,je,me,i,j,k;floatcm,bm,cj,bj;ie=jm[meo-1][0];je=jm[meo-1][1];me=jm[meo-1][2];cm=cjz[je-1][0]-cjz[ie-1][0];bm=cjz[ie-1][1]-cjz[je-1][1];cj=cjz[ie-1][0]-cjz[me-1][0];bj=cjz[me-1][1]-cjz[ie-1][1];ae=(bj*cm-bm*cj)/2;if(iask>1){for(i=0;i<3;i++){for(j=0;j<6;j++)b[i][j]=0;}b[0][0]=-bj-bm;b[0][2]=bj;b[0][4]=bm;b[1][1]=-cj-cm;b[1][3]=cj;b[1][5]=cm;b[2][0]=-cj-cm;b[2][1]=-bj-bm;b[2][2]=cj;b[2][3]=bj;b[2][4]=cm;b[2][5]=bm;for(i=0;i<3;i++){for(j=0;j<6;j++)b[i][j]=b[i][j]/(2*ae);}d[0][0]=eo/(1-un*un);d[0][1]=eo*un/(1-un*un);d[0][2]=0;d[1][0]=eo*un/(1-un*un);d[1][1]=eo/(1-un*un);d[1][2]=0;d[2][0]=0;d[2][1]=0;d[2][2]=eo/(2*(1+un));for(i=0;i<3;i++){for(j=0;j<6;j++){s[i][j]=0;for(k=0;k<3;k++)s[i][j]=s[i][j]+d[i][k]*b[k][j];}}if(iask>2){for(i=0;i<6;i++){for(j=0;j<6;j++){eke[i][j]=0;for(k=0;k<3;k++)eke[i][j]=eke[i][j]+s[k][i]*b[k][j]*ae*te;}}}}}voidTOTSTI(void)/*總的剛度矩陣*/{inti,j,ii,jj,meo,lii,l,lz,ldh,ld;for(i=0;i<nj2;i++){for(j=0;j<ndd;j++)tkz[i][j]=0;}for(meo=1;meo<=ne;meo++){ELEST(meo,3);for(i=1;i<4;i++){for(ii=1;ii<3;ii++){lii=2*(i-1)+ii;ldh=2*(jm[meo-1][i-1]-1)+ii;for(j=1;j<4;j++){for(jj=1;jj<3;jj++){l=2*(j-1)+jj;lz=2*(jm[meo-1][j-1]-1)+jj;ld=lz-ldh+1;if(ld>0)tkz[ldh-1][ld-1]=tkz[ldh-1][ld-1]+eke[lii-1][l-1];}}}}}}voidLOAD()/*添加載荷*/{inti,j,meo,ie,je,me;floatp0;for(i=0;i<nj2;i++)p[i]=0;if(npj>0){for(i=1;i<npj;i++){j=(int)(pj[i][1]);p[j-1]=pj[i][0];}}if(gama>0){for(meo=1;meo<=ne;meo++){ELEST(meo,1);p0=-gama*ae*te/3;ie=jm[meo-1][0];je=jm[meo-1][1];me=jm[meo-1][2];p[2*ie-1]=p[2*ie-1]+p0;p[2*je-1]=p[2*je-1]+p0;p[2*me-1]=p[2*me-1]+p0;}}}voidSUPPOR()/*添加約束*/{inti,j,mz,j0;for(i=1;i<=nz;i++){mz=nzc[i-1];tkz[mz-1][0]=1;for(j=1;j<ndd;j++)tkz[mz-1][j]=0;if(mz<ndd)j0=mz;elsej0=ndd;for(j=2;j<=j0;j++)tkz[mz-j][j-1]=0;p[mz-1]=0;}}voidSOLVEQ()/*求解過程*/{intk,i,j,im,l,m,ii,jo,lh;floatc;for(k=1;k<=nj2-1;k++){if(nj2>=k+ndd-1)im=nj2;elseim=k+ndd-1;for(i=k+1;i<=im;i++){l=i-k+1;c=tkz[k-1][l-1]/tkz[k-1][0];for(j=1;j<=ndd-l+1;j++){m=j+i-k;tkz[i-1][j-1]=tkz[i-1][j-1]-c*tkz[k-1][m-1];}p[i-1]=p[i-1]-c*p[k-1];}}p[nj2-1]=p[nj2-1]/tkz[nj2-1][0];for(ii=1;ii<=nj2-1;ii++){i=nj2-ii;if(ndd<=nj2-i+1)jo=ndd;elsejo=nj2-i+1;for(j=2;j<=jo;j++){lh=j+i-1;p[i-1]=p[i-1]-tkz[i-1][j-1]*p[lh-1];}p[i-1]=p[i-1]/tkz[i-1][0];}printf("\n");printf("Shuchujiedianweiyi\n");printf("uv\n");for(i=0;i<nj2;i++){printf("%16f",p[i]);if(i%2==1)printf("\n");}printf("\n");}voidSTRESS()/*應(yīng)力計(jì)算*/{doubleceta,ceta1,sigx,sigy,toxy,py1,ry1,sig1,sig2;intmeo,i,j,lh,ldh;floatwy[6],yl[3];printf("Shuchuyinglizhuyinglizhupingmianjiao\n");for(meo=1;meo<=ne;meo++){ELEST(meo,2);for(i=1;i<=3;i++){for(j=1;j<=2;j++){lh=2*(i-1)+j;ldh=2*(jm[meo-1][i-1])+j;wy[lh-1]=p[ldh-1];}}for(i=1;i<=3;i++){yl[i-1]=0;for(j=1;j<=6;j++)yl[i-1]=yl[i-1]+s[i-1][j-1]*wy[j-1];}sigx=yl[0];sigy=yl[1];toxy=yl[2];py1=(sigx+sigy)/2;ry1=sqrt((sigx-sigy)*(sigx-sigy)/4+toxy*toxy);sig1=py1+ry1;sig2=py1-ry1;if(sigy==sig2)ceta=0;else{ceta1=toxy/(sigy-sig2);ceta=90-57.29578*atan(ceta1);}printf("Danyuan%d\n",meo);printf("sigx=%12fsigy=%12ftoxy=%12fsig1=%12fsig2=%}}voidmain()/*主程序*/{printf("*****Pingmianyouxianyuanjisuanchengxu*****\n");printf("\n");DATA();if(ind!=0){eo=eo/(1-un*un);un=un/(1-un);}TOTSTI();LOAD();SUPPOR();SOLVEQ();STRESS();getchar();printf("\n");printf("R-jixie-0630618030303Zhangchunqing\n");printf("2010/11/20\n");getchar();}三、工程算例如圖所示有一剛性基礎(chǔ)上的三角形壩，受齊頂水壓力作用，以平面問題計(jì)算取μ＝1/6,將結(jié)構(gòu)分成四個(gè)單元，按圖中所示的網(wǎng)格劃分，用程序計(jì)算結(jié)點(diǎn)位移及各單元中的應(yīng)力，設(shè)彈性模量E＝1.0。（P1=15kN,P2=45kN）圖9工程算例圖為了能將程序執(zhí)行后的結(jié)果保存在文本文件中，我們?cè)赪indows的cmd工具中，運(yùn)行編譯后的可執(zhí)行文件yxy.exe，并將生成的數(shù)據(jù)錄入aa.txt的文本文件。先將編譯后的可執(zhí)行文件重命名為yxy.exe并放在D盤根目錄下，具體方法如下：?jiǎn)?dòng)windows的cmd工具，并cd到D盤根目錄下；在D:\>的提示符后輸入yxy.exe>>D:\aa.txt，回車執(zhí)行；在cmd對(duì)話框中輸入程序執(zhí)行所需要的數(shù)據(jù)，注意此時(shí)由于程序中所要在屏幕上打印的內(nèi)容全部都錄入了aa.txt的文本文件，在對(duì)話框中沒有輸入提示的語句，在輸入時(shí)要清除輸入的流程；回車運(yùn)行程序至完畢后會(huì)有提示符D:\>，此時(shí)關(guān)閉對(duì)話框，程序運(yùn)行的結(jié)果aa.txt既在D盤根目錄下。輸入數(shù)據(jù)：基本參數(shù)結(jié)點(diǎn)數(shù)nj單元數(shù)ne支桿數(shù)nz半帶寬ndd結(jié)點(diǎn)載荷數(shù)npj問題類型碼ind646820其它數(shù)據(jù)彈性模量eo泊松比un材料密度gama單元厚度ae1.00.1670.01.0單元節(jié)點(diǎn)碼數(shù)組JM矩陣｛｛1，2，3｝，｛2，4，5｝，｛2，5，3｝，｛3，5，6｝｝結(jié)點(diǎn)坐標(biāo)數(shù)組CJZ矩陣｛｛0.0，6.0｝，｛0.0，3.0｝，｛2.0，3.0｝，｛0.0，0.0｝，｛2.0，0.0｝，｛4.0，0.0｝｝支承數(shù)組NZC矩陣｛｛7｝，｛8｝，｛9｝，｛10｝，｛11｝，｛12｝｝結(jié)點(diǎn)載荷數(shù)組PJ矩陣｛｛0.0，0.0｝，｛15.0，1.0｝，｛45.0，3.0｝｝運(yùn)行所得結(jié)果如下*****Pingmianyouxianyuanjisuanchengxu*****Qingshuru6gecanshuQingshuru4gejisuancanshunj=6ne=4nz=6ndd=8npj=2ind=0eo=1.000000un=0.167000gama=0.000000te=1.000000QingshuruJMjuzhengQingshuruCJZjuzhengQingshuruNZCjuzhengQingshuruPJjuzhengCJZjuzhengruxia0.0000006.0000000.0000003.0000002.0000003.0000000.0000000.0000002.0000000.0000004.0000000.000000Shuchujiedianweiyiuv275.97079545.53725470.00800345.0285766