平面四节点等参单元分析程序(共27页).doc

《平面四节点等参单元分析程序(共27页).doc》由会员分享，可在线阅读，更多相关《平面四节点等参单元分析程序(共27页).doc（27页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上变分原理与有限元大作业平面四节点等参单元分析程序姓 名：潘 清学 号：SQ完成时间：2011-4-26一、 概述通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、ABAQUS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时

2、间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。因此，必须寻找新的编程语言。随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。现在可编写工程应用软件的计算机语言较多，其中C语言是一

3、个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。C语言最初是为操作系统、编译器以及文字处理等编程而发明的。随着不断完善，它已应用到其它领域，包括工程应用软件的编程。近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。用C 编写的程序，可移植性极好，几乎不用作多少修改，就可在任何一台装有ANSI、C编译器的计算机上运行。C既是高级语言，也是低级语言，也就是说，可用它作数值计算，也可用它对计算机存储进行操作。二、 编程思想本程序采用C语言编程，编制平面四边形四节点等参元程序，用

4、以求解平面结构问题。程序采用二维等带宽存储整体刚度矩阵，乘大数法引入约束，等带宽高斯消去法求解位移，然后求中间高斯点的应力，最后用绕节点平均法讲单元应力等效到节点上，再将结果写到tecplot文件中。在有限元程序中，变量数据需赋值的可分为节点信息，单元信息，载荷信息等。对于一个节点来说，需以下信息：节点编号（整型），节点坐标（实型），节点已知位移（实型），节点载荷（实型），边界条件（实型）等。同样，对于一个单元来说，需以下信息：单元的节点联接信息（整型），材料信息（弹性模量，泊松比等）（实型）等。在FORTRAN 程序中，以上这些变量混合在一起，很难辨认，使程序的可读性不好，如需要进行单元网络

5、的自适应划分，节点及单元的修改将非常困难。在进行C语言编译过程中，采用结构struct 使每个节点信息存储在一个结构体数组中，提高程序的可读性，使数据结构更趋于合理。三、平面四节点等参单元介绍四节点等参单元实际单元与基本单元的映射关系如Error! Reference source not found.所示图 31坐标的映射关系为：其位移模式和坐标的映射有相同的插值函数，形函数为：单元应变矩阵为：上式一般简写为：其中的子块矩阵为由于是、的函数，在中的、要按照复合函数来求导，即从而有因此，单元应力矩阵为单元刚度矩阵为其中积分采用三点高斯积分，（高斯积分点的总数），和或是加权系数，和是单元内的坐标

6、.。对于三点高斯积分，高斯积分点的位置： ，。结构刚度矩阵结构结点荷载列阵 注意，对于上两式中的理解不是简单的叠加而是按照对应的自由度集成。总刚平衡方程从式上式求出：四、有限元分析的模块组织一个典型的有限元分析过程主要包括以下几个步骤：1) 读输入数据，定义节点及单元数组。2) 由边界条件计算方程个数，赋值荷载列阵。3) 读入在带状存储的总刚度矩阵中单元和载荷信息。4) 定义总刚度阵数组。5) 组装总刚度阵。6) 解方程组。输入边界条件(力、位移)形成各荷载工况的节点荷载阵总刚分解 回代求出位移及输出 计算应变、应力形成单元刚阵单刚向总刚投放坐标变换输入原始参数计算总刚规模形成总刚方程向总节点

7、荷载阵投放形成单元荷载阵调整几何、弹性矩阵调整单元位移列阵其流程图可见下图五、程序变量及函数说明1、控制信息np：结构节点总数ne：结构离散单元总数nr1，nr2：总的约束的节点数，nr1，x方向；nr2，y方向nd：每个单元的节点数nf：每个节点的自由度数ld：集中力载荷的个数nm：材料的种类nu1，nu2：非零位移边界条件的节点数，nu1，x方向；nu2，y方向u1，u2：非零位移的大小，u1，x方向；u2，y方向n=nf*np ：结构的节点位移总数ndf=nd*nf ：每个单元的节点自由度数2、输入的原始数据x(np):节点的x方向坐标y(np):节点的y方向坐标me（nd，ne）：单元

8、节点的总体编号nrr(nr1+nr2) ：约束为零的位移所对应的总体位移编号p（n）：载荷向量nu（nu1+nu2）：位移载荷mat（6，nm）：材料参数3、程序中的其他标识符LD（n）：存放结构刚度阵所以主对角线元素在A（nn）中的序号IS（ndf）：单元节点位移和节点力在总体位移阵列和载荷阵列中对应的序号EK（ndf，ndf）：总体坐标系下的单元刚度矩阵A（nn）：架构刚度阵下三角变带宽一维压缩存储的数组nn：数组A的元素个数RSTG(3):高斯积分点的值H(3):高斯积分点的加权系数S（6，ne）：各单元的应力分量XJAC（2,2）：雅阁比矩阵RJAC（2,2）：雅阁比矩阵的逆PN（2,

9、4）：4个节点处形函数对局部坐标的导数DNX（2,4）：4个节点处形函数对整体坐标的导数FUN（4）：形函数的值六、计算流程图输出节点位移计算局部坐标系下的单元刚度矩阵坐标变换矩阵计算单元的IS(ndf)数组取出单元在总体坐标系下的节点位移单元刚度矩阵向结构刚度矩阵叠加形成单元的IS(ndf)数组计算局部坐标下的单元刚度矩阵坐标变换矩阵总体坐标下的单元刚度矩阵形成LD（n）数组输入原始数据n=nfnpndf=ndnf开始 输入控制信息结束计算局部坐标系下单元节点力应力计算局部坐标系下单元节点位移求解线性方程组求得结构节点位移进行约束处理 输出节点位移和应力七、计算结果与abaqus分析结果的比

10、较1、中间带圆孔平面应力板的分析。宽40m，长50m，圆孔位于板中心，半径为4m，承受水平方向位移载荷，E=200Pa，取m。用abaqus建模离散，并计算。再讲abaqus离散的节点和单元信息拷贝到本程序的输入文件中，用该程序计算，结果输出成tecplot文件，用tecplot可以查看结果。与abaqus的计算结果进行比较，位移和应力云图如下（左边是程序计算结果，右边是abaqus结果，下同）：2、纤维增强复合材料平面应力板的分析。宽40m，长50m，增强纤维位于板中心，纤维半径为8m，承受水平方向位移载荷，基体材料E=20000Pa，纤维材料E=8000Pa，取m。3、半无限大含裂纹板的应

11、力分析宽40m，长60m，裂纹位于板左侧中间位置，裂纹长10m，承受竖直方向位移载荷，E=200Pa，取m。结果分析比较：通过以上三个算例发现该程序可以用于计算单材料、双材料、带孔、含裂纹等各种平面问题，加载条件可以是加集中力和加位移，因此，该程序的适用范围还是比较广的。以上三个算例的自编程序所得结果与abaqus分析结果进行比较发现，两者的计算结果很接近，而且自编程序对于孔边应力集中和裂尖应力集中都能很好的表达，说明该程序有很好的精度。第三个算例在裂尖处数值上有些区别，但总的分布规律还是很吻合的，主要是因为本程序是用四节点等参单元，对于应力的奇异性表达效果还不是很好。八、附录附录一：程序代码

12、#include#include#include#include#includefloat *float_two_array_malloc(int m,int n) /实型二维动态数组分配 float *a; int i,j; a=(float *)malloc(m*sizeof(float *); for (i=0;im;i+) ai=(float *)malloc(n*sizeof(float);for (j=0; jn; j+)aij=0; return a;float *float_one_array_malloc(int n) /实型一维动态数组分配 float *a; int j;

13、 a=(float *)malloc(n*sizeof(float);for (j=0; jn; j+)aj=0; return a;int *int_two_array_malloc(int m,int n) /整型二维动态数组分配 int *a; int i,j; a=(int *)malloc(m*sizeof(int *); for (i=0;im;i+) ai=(int *)malloc(n*sizeof(int);for (j=0; jn; j+)aij=0; return a;int *int_one_array_malloc(int n) /整型一维动态数组分配 int *a;

14、 int j; a=(int *)malloc(n*sizeof(int);for (j=0; jn; j+)aj=0; return a;/全局变量定义int np,ne,nr1,nr2,nd,nf,n,ndf,ld,nm,nu1,nu2;float *x,*y,*p,*pp,u1,u2; /x,y,me,nrr,p为输入信息float coords24,det,fun4,pn24,xjac22,rjac22,dnx24; /求单刚定义的变量int *LD,*me,*nrr,*nu,*IS,nn; /nn为变带宽一维组中元素个数 /IS表示单元节点位移在结构位移阵列中对应的序号float *

15、EK,*s,*A,*ns,*mat; /EK为总体坐标系下单元的刚度矩阵 /A存储结构刚度矩阵的一维变带宽数组 /s存放各单元应力中心点的应力/ns存放各节点的应力void input() /输入函数定义int i,j;FILE *fp1;fp1=fopen(input.dat,r); /要输入的内容放在input.dat中if(fp1=NULL)printf(cannt open the file !n);exit(1);/控制信息读入 fscanf(fp1,%d,&np);fscanf(fp1,%d,&ne); fscanf(fp1,%d,&nr1); fscanf(fp1,%d,&nr2

16、); fscanf(fp1,%d,&nd); fscanf(fp1,%d,&nf);fscanf(fp1,%d,&ld); fscanf(fp1,%d,&nm); fscanf(fp1,%d,&nu1); fscanf(fp1,%f,&u1); fscanf(fp1,%d,&nu2); fscanf(fp1,%f,&u2); n=nf*np; ndf=nf*nd;/printf(%dt%ft%dt%fn,nu1,u1,nu2,u2);/动态数组生成x=float_one_array_malloc(np); y=float_one_array_malloc(np);nrr=int_one_arr

17、ay_malloc(nr1+nr2); /位移约束为零的节点，nr1-x方向约束，nr2-y方向约束nu=int_one_array_malloc(nu1+nu2); /位移约束非零的节点pp=float_two_array_malloc(ld,2); /非零载荷 （集中力）,ppld0对应的自由度序号，ppld1集中力大小p=float_one_array_malloc(n); /载荷矩阵me=int_two_array_malloc(nd,ne); /单元节点的总体编号mat=float_two_array_malloc(4,nm);/单元材料数组，先定义为各向同性材料/原始数据读入 fo

18、r(i=0;inp;i+)fscanf(fp1,%f,&xi);fscanf(fp1,%f,&yi); for(i=0;ine;i+) for(j=0;jnd;j+)fscanf(fp1,%d,&meji); meji=meji-1; /abaqus模型节点编号从1开始，C语言中从0开始for(i=0;inr1;i+) /读入固定位移约束条件fscanf(fp1,%d,&nrri);/x方向约束的节点号nrri=(nrri-1)*2; /转化为自由度序号，abaqus模型节点编号从1开始，C语言中从0开始for(i=0;inr2;i+) fscanf(fp1,%d,&nrrnr1+i);/y方

19、向约束的节点号，abaqus模型节点编号从1开始，C语言中从0开始nrrnr1+i=(nrrnr1+i-1)*2+1; for(i=0;ild;i+) /读入集中力载荷fscanf(fp1,%f%f,&ppi0,&ppi1); for(i=0;in;i+) pi=0.0;for(i=0;ild;i+) j=int(ppi0); pj=ppi1; for(i=0;inm;i+) /读入材料信息 for(j=0;j4;j+)fscanf(fp1,%f,&matji); mat0i=mat0i-1;mat1i=mat1i-1; /abaqus模型节点编号从1开始，C语言中从0开始for(i=0;in

20、u1;i+) /读入非零位移载荷fscanf(fp1,%d,&nui);/x方向约束的节点号nui=(nui-1)*2; /转化为自由度序号，abaqus模型节点编号从1开始，C语言中从0开始for(i=0;inu2;i+) fscanf(fp1,%d,&nunu1+i);/y方向约束的节点号，abaqus模型节点编号从1开始，C语言中从0开始nunu1+i=(nunu1+i-1)*2+1;fclose(fp1);printf(读入信息完成.n);void output() /结果输出函数 int i,j;FILE *fp2;fp2=fopen(output.dat,w); /运算结果放在ou

21、tput.dat中if(fp2=NULL)printf(cannt open the file !n);exit(1); fprintf(fp2,位移为：n); /输出位移for(i=0;in;i+=2) fprintf(fp2,u%d=%f v%d=%fn,i/2+1,pi,i/2+1,pi+1);fprintf(fp2,n各单元应力：Sx,tSy,tSxy,tS1,tS2,tSmisesn); /输出单元应力for(i=0;ine;i+) fprintf(fp2,%dt,i+1);for(j=0;j6;j+)fprintf(fp2,%.4ft,sji);fprintf(fp2,n);fcl

22、ose(fp2);printf(写入output.dat完成.n);void tecplot() /结果生成tecplot文件 int i,j;float bl=1.0; /位移显示比例FILE *fp3;fp3=fopen(tecplot.plt,w); /运算结果放在tecplot.plt中if(fp3=NULL)printf(cannt open the file !n);exit(1); fprintf(fp3, TITLE = test n); fprintf(fp3, VARIABLES=x,y,Ux,Uy,Sx,Sy,Sxy,Smax,Smin,Smisesn); fprintf

23、(fp3, ZONE N=%d,E=%d,F=FEPOINT,ET=quadrilateraln,np,ne); for(j=0;jnp;j+)fprintf(fp3,%.4ft%.4ft%.4ft%.4ft,xj+bl*pj*2,yj+bl*pj*2+1,pj*2,pj*2+1);for(i=0;i6;i+)fprintf(fp3,%.4ft,nsij);fprintf(fp3,n); for(i=0;ine;i+)for(j=0;jnd;j+)fprintf(fp3,%dt,meji+1);fprintf(fp3,n);fclose(fp3);printf(写入tecplot.plt完成.

24、n);void fld() /生成LD数组，存放组对角线元素在A中的位置 int k,i,j,j1,ig,nw; LD=int_one_array_malloc(n); LD0=0; for(k=0;knp;k+) ig=np; for(i=0;ine;i+) for(j=0;jnd;j+) if(meji=k)for(j1=0;j1nd;j1+)if(mej1iig) ig=mej1i; for(i=0;inf;i+)j=nf*k+i;nw=nf*(k-ig)+i+1;if(j!=0) LDj=LDj-1+nw; nn=LDn-1+1; /nn为变带宽一维组中元素个数,LD下标从零开始pri

25、ntf(生成LD数组完成.n);void fis(int ie) /形成单元ie的节点自由度号IS数组 int i,j,ii; /IS=int_one_array_malloc(ndf); for(i=0;ind;i+) for(j=0;jnf;j+) ii=nf*i+j; ISii=nf*meiie+j; void fdnx(int ie,float gs,float gr)/-计算形函数对整体坐标的导数dnx24 int m,u,i,v; float g1,g2,xi4,eta4,rec; xi0=-1;xi1=1;xi2=1;xi3=-1;eta0=-1;eta1=-1;eta2=1;e

26、ta3=1;/-读入单元节点坐标for (i=0;i4;i+)coords0i=xmeiie;coords1i=ymeiie;/-计算形函数和|J|for(m=0;m4;m+)g1=(1+xim*gr);g2=(1+etam*gs);funm=0.25*g1*g2;pn0m=0.25*xim*g2;pn1m=0.25*etam*g1;for(u=0;u2;u+)for(v=0;v2;v+)det=0.0;for(m=0;m4;m+)det=det+pnum*coordsvm;xjacuv=det;det=xjac00*xjac11-xjac01*xjac10;/-雅格比矩阵求逆-rec=1.0

27、/det;rjac00=rec*xjac11;rjac11=rec*xjac00;rjac01=-rec*xjac01;rjac10=-rec*xjac10;/-形函数对整体坐标的导数-for(u=0;u4;u+)for(v=0;v2;v+)dnxvu=0.0;for(m=0;m2;m+)dnxvu=dnxvu+rjacvm*pnmu;void fek(int ie) /形成总体坐标系下四节点平面等参单元ie的刚度矩阵EK（8，8） int i,j,ii,jj,i1,i2,j1,j2,k,l; float D33,E,MU,d1,d2,d3,dxx,dxy,dyx,dyy,gs,gsh,gr,

28、grh,rstg3,h3,dnix,dniy,dnjx,dnjy;rstg0=-0.;rstg1=0;rstg2=-rstg0;h0=0.6;h1=0.9;h2=h0;/-形成D矩阵- for(i=0;inm;i+) /读入材料if(ie=mat0i)E=mat2i;MU=mat3i;D00=E/(1.0-MU*MU);D11=D00;D01=D00*MU;D10=D00*MU;D02=0;D12=0;D20=0;D21=0;D22=E/(2*(1+MU);d1=D00;d2=D01;d3=D22;/-单元刚度矩阵清零- for(i=0;indf;i+)for(j=0;jndf;j+)EKji

29、=0.0;/-for(i=0;ind;i+) /-形成I1，I2，J1，J2-ii=2*i;i1=ii;i2=ii+1;for(j=0;jnd;j+)jj=2*j;j1=jj;j2=jj+1;/-dxx,dxy,dyx,dyy- dxx=0;dxy=0;dyx=0;dyy=0; for(k=0;k3;k+)gs=rstgk; /高斯点的值gsh=hk; /加权系数for(l=0;l3;l+)gr=rstgl; grh=hl; fdnx(ie,gs,gr);dnix=dnx0i;dniy=dnx1i;dnjx=dnx0j;dnjy=dnx1j;dxx=dxx+dnix*dnjx*det*gsh*

30、grh;dxy=dxy+dnix*dnjy*det*gsh*grh;dyx=dyx+dniy*dnjx*det*gsh*grh;dyy=dyy+dniy*dnjy*det*gsh*grh;/-分块形成单元刚度矩阵-EKi1j1=dxx*d1+dyy*d3;EKi2j2=dyy*d1+dxx*d3;EKi1j2=dxy*d2+dyx*d3;EKi2j1=dyx*d2+dxy*d3;void fasum(float*A) /将单元刚度矩阵叠加到结构刚度矩阵中并存到A（nn）中 int i,j,ni,nj,ij;for(i=0;indf;i+)for(j=0;jndf;j+)ni=ISi;nj=IS

31、j; if(nj=ni) ij=LDni-(ni-nj); Aij=Aij+EKij; void fr(float*A)/约束处理子程序int i,j,t;for(i=0;inr1+nr2;i+) /固定位移约束处理j=nrri;t=LDj;At=(float)1.0E30;pj=0.0;for(i=0;inu1;i+) /非零位移约束处理x方向j=nui;t=LDj;At=(float)1.0E30;pj=(float)1.0E30*u1;for(i=nu1;inu1+nu2;i+) /非零位移约束处理y方向j=nui;t=LDj;At=(float)1.0E30;pj=(float)1.0

32、E30*u2;printf(约束处理完成.n);void cholesky(float*A) /用cholesky方法求解线性方程组 int i, j, k,ij,mi; float* z; /增广矩阵 z = float_two_array_malloc(n, n+1); /* 将原始数据导入z中 */ for(i=0; in; i+) for(j=0; j=i; j+) if(i=0) mi=0; else mi=i-(LDi-LDi-1)+1;/第i行第一个非零元素的列号if(jmi) zij=0;else ij=LDi-(i-j); zij=Aij; for(i=0;in;i+)for

33、(j=i+1;jn;j+)zij=zji; for(i=0;in;i+)zin=pi; float temp; /* 分解过程 */ for(k=0; kn; k+) temp = 0; for(i=0; ik; i+) temp = temp + zik * zik; zkk = (float)sqrt(zkk - temp); for(i=k+1; in+1; i+) temp = 0; for(j=0; j=0; k-) temp = 0; for(i=k+1; in; i+) temp = temp + zki * zin; zkn = (zkn - temp) / zkk; /*将位

34、移阵列存入p中*/ for(i=0; in; i+) pi=zin; /*释放z数组*/for (i=0;in;i+) delete zi; delete z;printf(求解线性方程组完成.n); void stress(int ie)/ 求单元中心点的应力 int i,j,ii,j1,j2;float d1,d2,d3,D33,B38,r8,sig3,E,MU,ss,rr,bi,ci,h1,h2,sx,sy,sxy;/-形成D矩阵- for(i=0;inm;i+) /读入材料if(ie=mat0i)E=mat2i;MU=mat3i;D00=E/(1.0-MU*MU);D11=D00;D01=D00*MU;D10=D00*MU;D02=0;D12=0;D20=0;D21=0;D22=E/(2*(1+MU);d1=D00;d2=D01;d3=D22;ss=0.0;rr=0.0;fdnx(ie,ss,rr);/-