有限元弹性力学上机报告

弹性力学及有限元基础上机实践报告指导老师：班级：机械（茅）101学号： 1004010510姓名：一．程序功能（1）本程序只用于三节点三角形单元，它可计算平面应力问题，亦可计算平面应变问题。

这两类问题用类型码IND来区别：IND输入0的数据——为解平面应力问题。

IND输入1的数据——为解平面应变问题。

（2）程序中的载荷包括结点载荷和自重载荷两种，如有其他非节点载荷，则应事先换算成等效结点载荷。

（3）支撑方式，可以是在任一结点的水平和数值支撑，支撑数量不限。

但平面弹性体至少应具有保证弹性体几何不变性的三个独立的支撑。

（4）弹性体由单一的材料组成，即只有一组E,G,u弹性常数。

（5）弹性体计算简图的厚度是常量。

（6）直要计算机内存允许，本程序对对单元总数，结点总数等没有限制。

（7）本程序输出全部结点位移，对每个但愿输出应力分量，两个主应力及其主平面角。

二．程序中的记号说明nj结点个数ne单元个数nz约束个数ndd半带宽ind问题类型码ind=0平面应力问题ind=1nj2位移分量个数eo弹性模量un泊松比gama材料容重比γte单元厚度ae单元面积jm(ne,3)单元结点码数组nzc(nz)支承数组cjz(nj,2)结点坐标数组pj(npj,2)结点载荷数组b(3,6)几何矩阵d(3,3)弹性矩阵s(3,6)应力矩阵tkz(nj2,ndd)半带储存的整刚矩阵eke(6,6)单元刚度矩阵p(nj2)载荷向量，后来存放结点位移npj1结点载荷个数加1。

是载荷数组的最大行代码meo任一单元码iask子程序elest中的形参，它是计算信息码iask=1求单元面积aeiask=2求应力矩阵[s]iask=3求单元刚度[k]ie,je,me单元meo的三个角点码cm,bm,cj,bj计算形函数中的常cm，bm，cj，bjI,j,k,ii,jj循环参数lh单元行码ldh半带存储的整体刚度矩阵的行码l单元列码lz整体刚度矩阵的列码ld半带存储的整体刚度矩阵的列码i1工作单元pe自重的等效结点荷载mz支杆相应的位移分量码jo最大列码j1工作单元im最大行码c系数比值ld1m工作单元wy(6)单元结点位移向量yl(3)应力向量sigxsigytoxy应力分量“d”x dy 套xypyl平均应力ryl应力圆半径sig工作单元sig1sig2最大、最小主应力geta1工作单元geta主平面角三．子程序名data输入elest（meo，iask）iask=1iask=2iask=3totsti整体刚度矩阵load荷载向量suppor支撑条件solveq位移stress应力四．程序流程及框图此程序采用本人比较熟悉C语言编写。

C语言相比于FORTRAN语言拥有更强大的实用编程环境，同时由于提供了指针变量、动态内存分配函数和结构体变量等，使其编写的软件在维护性、可读性和内存利用效率方面具有明显的优势。

1.程序总框图3图1 程序总框图整个程序由一个主调主程序(主函数main（）)和七个子程序组成，其中数据子程序DATA( )用来接受输入的参数和变量。

单刚子程序ELEST( )为含有参数的函数，功能控制参数iask可取1、2、3，分别计算相应单元的面积、应力矩阵和单元刚度矩阵，主程序不直接调用它，而通过总刚子程序、载荷子程序和求应力子程序间接调用。

总刚子程序TOTSTI( )用来合成总刚矩阵，载荷子程序LOAD( )用来计算合成载荷，支承子程序SUPPOR( )用来引入约束，解方程子程序SOLVEQ( )用来求解并输出各个结点的位移，求应力子程序STRESS( )用来计算和输出应力、主应力及主平面角。

2.主函数及各子程序流程图图2 void main( )和void DA TA( )图4 void TOTSTI( ) 图5 void LOAD( )图6 void SUPPOR( )图7 void SOLVEQ( )五．有限元计算源程序#include "stdio.h"#include "math.h"int nj,ne,nz,ndd,ind,nj2;int jm[100][3],nzc[200],npj1,npj;float eo,un,gama,te,ae;float cjz[100][2],pj[100][2];float b[3][6],d[3][3],s[3][6],eke[6][6],tkz[200][20],p[200];void DATA()/*数据输入函数*/{int i,j;printf("please input 6 base data\n");scanf("%d,%d,%d,%d,%d,%d",&nj,&ne,&nz,&ndd,&npj,&ind);nj2=nj*2;npj1=npj+1;getchar();printf("please input 4 parameters\n");scanf("%f,%f,%f,%f",&eo,&un,&gama,&te);printf("nj=%d,ne=%d,nz=%d,ndd=%d,npj=%d,ind=,%d\n",nj,ne,nz,ndd,npj,ind);printf("eo=%f,un=%f,gama=%f,te=%f\n",eo,un,gama,te);getchar();printf("please input JM matrix\n");for(i=0;i<ne;i++){for(j=0;j<3;j++)scanf("%d",&jm[i][j]);getchar();}printf("please input CJZ matrix\n");for(i=0;i<nj;i++){for(j=0;j<2;j++)scanf("%f",&cjz[i][j]);getchar();}printf("please input NZC matrix\n");for(i=0;i<nz;i++)scanf("%d",&nzc[i]);getchar();printf("please input NIZ matrix\n");for(i=0;i<npj+1;i++){for(j=0;j<2;j++)scanf("%f",&pj[i][j]);getchar();}printf("CJZ matrix as follows\n");for(i=0;i<nj;i++){for(j=0;j<2;j++)printf(" %f\t",cjz[i][j]);printf("\n");}}void ELEST(int meo,int iask)/*单元刚度矩阵*/ {int ie,je,me,i,j,k;float cm,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;}}}}}void TOTSTI()/*总体刚度矩阵*/{int i,j,ii,jj,meo,lii,l,lz,ldh,ld;for(i=0;i<nj;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];}}}}}}void LOAD()/*添加载荷*/{int i,j,meo,ie,je,me;float p0;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;}}}void SUPPOR()/*添加约束*/{int i,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;}}void SOLVEQ()/*求解过程*/{int k,i,j,im,l,m,ii,jo,lh;float c;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("Output Result\n");printf(" u v\n");for(i=0;i<nj2;i++){printf(" %16f",p[i]);if(i%2==1)printf("\n");}printf("\n");}void STRESS()/*应力计算*/{double ceta,cetal,sigx,sigy,toxy,pyl,ryl,sig1,sig2;int meo,i,j,lh,ldh;float wy[6],yl[3];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]-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];pyl=(sigx+sigy)/2;ryl=sqrt((sigx-sigy)*(sigx-sigy)/4+toxy*toxy);sig1=pyl+ryl;sig2=pyl-ryl;if(sigy==sig2)ceta=0;else{cetal=toxy/(sigy-sig2);ceta=90-57.29578*atan(cetal);}printf("meo=%d\n sigx=%f sigy=%f toxy=%f sig1=%f sig2=%f ceta=%f\n", meo,sigx,sigy,toxy,sig1,sig2,ceta);}}void main()/*主程序*/{printf("******************Plane Limited Programes******************\n");printf("\n");DA TA();if(ind!=0){eo=eo/(1-un*un);un=un/(1-un);}TOTSTI();LOAD();SUPPOR();SOLVEQ();STRESS();getch();printf("\n");printf("\n");printf(" 机械（茅)101 1004010510 王洋\n");printf(" 2013.12.16\n");六．工程算例一个结构网格，其中9个单元10各节点，4个枝干，在节点10作用竖向载荷10KN.按平面应力问题计算。