有限元上机报告

弹性力学及有限元基础上机实践报告PAGEPAGE1 在以常应变三角形单元计算弹性力学平面问题原理的基础上，编写可计算平面应力问题，亦可计算平面应变问题的通用程序。给出程序设计的流程图，程序源代码和程序应用的工程算例及其结果分析。一、程序流程图此程序采用本人比较熟悉C语言编写。C语言相比于FORTRAN语言拥有更强大的实用编程环境，同时由于提供了指针变量、动态内存分配函数和结构体变量等，使其编写的软件在维护性、可读性和内存利用效率方面具有明显的优势。程序总框图图1程序总框图 整个程序由一个主调主程序(主函数main（）)和七个子程序组成，其中数据子程序DATA()用来接受输入的参数和变量。单刚子程序ELEST()为含有参数的函数，功能控制参数iask可取1、2、3，分别计算相应单元的面积、应力矩阵和单元刚度矩阵，主程序不直接调用它，而通过总刚子程序、载荷子程序和求应力子程序间接调用。总刚子程序TOTSTI()用来合成总刚矩阵，载荷子程序LOAD()用来计算合成载荷，支承子程序SUPPOR()用来引入约束，解方程子程序SOLVEQ()用来求解并输出各个结点的位移，求应力子程序STRESS()用来计算和输出应力、主应力及主平面角。2.主函数及各子程序流程图图2voidmain()和voidDATA()图3voidELEST(intmeo,intiask)图4voidTOTSTI() 图5voidLOAD()图6voidSUPPOR()图7voidSOLVEQ()图8voidSTRESS()二、程序源代码/*有限元计算程序*/#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()/*数据输入函数*/{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)/*单元刚度矩阵*/{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()/*应力计算*/{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();}三、工程算例如图所示有一刚性基础上的三角形坝，受齐顶水压力作用，以平面问题计算取μ＝1/6,将结构分成四个单元，按图中所示的网格划分，用程序计算结点位移及各单元中的应力，设弹性模量E＝1.0。（P1=15kN,P2=45kN）图9工程算例图为了能将程序执行后的结果保存在文本文件中，我们在Windows的cmd工具中，运行编译后的可执行文件yxy.exe，并将生成的数据录入aa.txt的文本文件。先将编译后的可执行文件重命名为yxy.exe并放在D盘根目录下，具体方法如下：启动windows的cmd工具，并cd到D盘根目录下；在D:\>的提示符后输入yxy.exe>>D:\aa.txt，回车执行；在cmd对话框中输入程序执行所需要的数据，注意此时由于程序中所要在屏幕上打印的内容全部都录入了aa.txt的文本文件，在对话框中没有输入提示的语句，在输入时要清除输入的流程；回车运行程序至完毕后会有提示符D:\>，此时关闭对话框，程序运行的结果aa.txt既在D盘根目录下。输入数据：基本参数结点数nj单元数ne支杆数nz半带宽ndd结点载荷数npj问题类型码ind646820其它数据弹性模量eo泊松比un材料密度gama单元厚度ae1.00.1670.01.0单元节点码数组JM矩阵｛｛1，2，3｝，｛2，4，5｝，｛2，5，3｝，｛3，5，6｝｝结点坐标数组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｝｝支承数组NZC矩阵｛｛7｝，｛8｝，｛9｝，｛10｝，｛11｝，｛12｝｝结点载荷数组PJ矩阵｛｛0.0，0.0｝，｛15.0，1.0｝，｛45.0，3.0｝｝运行所得结果如下*****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