八节点等参单元平面有限元(共37页)

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

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

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

4、序采用二维等带宽存储整体刚度矩阵，乘大数法引入约束，等带宽高斯消去法求解位移。在有限元程序中，变量数据需赋值的可分为节点信息，单元信息，载荷信息等。对于一个节点来说，需以下信息：节点编号（整型），节点坐标（实型），节点已知位移（实型），节点载荷（实型），边界条件（实型）和节点温度（实型）等。同样，对于一个单元来说，需以下信息：单元的节点联接信息（整型），单元类型信息（桁架、梁、板、壳等）（整型） ，单元特性信息（厚度、内力矩等）（实型），材料信息（弹性模量，泊松比等）（实型）等。在FORTRAN 程序中，以上这些变量混合在一起，很难辨认，使程序的可读性不好，如需要进行单元网络的自适应划分，节点

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

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

7、投放坐标变换输入原始参数计算总刚规模形成总刚方程向总节点荷载阵投放形成单元荷载阵调整几何、弹性矩阵调整单元位移列阵图 22其流程图可见Error! Reference source not found.。3. 程序变量及函数说明3.3.3.1. 主要变量说明：主要程序变量说明 wide分析模型的宽 hight分析模型的高wsize宽度方向网格划分尺寸hsize高度方向网格划分尺寸npoin节点总数nelem单元总数struct node500节点结构体变量struct element 500单元结构体变量ifpre500节点约束信息posgp3高斯积分点weigp3权系数gpcod29高斯积分

8、点总体坐标bmatx316单元变形矩阵dmatx33单元弹性矩阵xjacm22雅可比矩阵xjaci22雅可比矩阵的逆矩阵djacb雅可比矩阵行列式的值estif136单元刚度矩阵shape8单元形函数deriv28形函数对局部坐标的导数cartd28形函数对整体坐标的导数A30000总体刚度矩阵eload16等效节点荷载gpcod29高斯积分点的总体坐标3.1.3.2. 主要函数说明主要函数说明void meshing（ ）网格划分void coordinate( )节点整体坐标计算void input( )信息输入void stif( )形成单元刚度矩阵void sfr2（）计算形函数对当前

9、积分点及形函数对局部坐标的到数值void jacobi2( )形成雅可比矩阵void modps( )形成弹性矩阵void bmatps( )形成应变矩阵void load( )计算等效节点荷载void asem( )形成整体刚度矩阵和整体荷载列阵void solve( )求解整体方程void stre( )计算单元应力4. 程序流程图任何一个有限元程序处理过程，都可以划分为三个阶段：1) 前处理：读入数据和生成数据。2) 处理器：有限元的矩阵计算、组集和求解。3) 后处理：输出节点位移、单元应变等。为了更清楚地描述程序，我们给出了程序的流程图。 5. 程序应用与ANSYS分析的比较为了验证程

10、序的正确性，我们取了一个简单框架结构，分别用自编程序和ANSYS进行分析和比较。4.5.5.1. 问题说明Error! Reference source not found.所示一简支梁，高3m，长18m，承受均布荷载，取m，作为平面应力问题。图 51图 52由于结构和荷载对称，只对右边一半进行有限单元计算，如Error! Reference source not found.所示，而在y轴各节点布置水平连杆支座。5.2. ANSYS分析结果ANSYS分析中，采用plane82单元，在板单元上边施加均布荷载，在y轴上的各结点约束UX方向的自由度，在板单元右下角的结点约束UX自由度，结点布置、网

11、格划分方案如Error! Reference source not found.所示。图 53Error! Reference source not found.和Error! Reference source not found.分别给出了结构的位移图和应力云图。图 54 位移图图 55 各单元图从Error! Reference source not found.和Error! Reference source not found.可以看到，跨中的位移和正应力最大，与简支梁承受均布荷载，跨中挠度和正应力最大的力学常识相符合，可以初步认为模型是正确的。Error! Reference sou

12、rce not found.给出了的截面上的正应力和Error! Reference source not found.给出了处各横截面方向位移，其中 的单位为，的单位为。表格 1考查点的y/m1.51.00.50-0.5-1.0-1.5正应力-270.20-178.25-88.570.0088.57178.25270.21表格 2考查点的x/m01.53.04.56.07.59.0方向位移(10-6)-0.3485-0.3380-0.3069-0.2571-0.1914-0.113305.3. 自编程序分析结果用自编程序进行分析时，采用了整个结构模型进行计算，其他条件不变，因编程水平有限，在

13、后处理阶段没有给出节点处的位移与应力，而只能得到高斯积分点处的位移与应力，得到积分点处的应力后，根据数值分析相关知识采用了插值进行处理，从而得到的截面上的正应力和 处各横截面方向位移，分别列出表格如下：表格 3考查点的y/m1.51.00.50-0.5-1.0-1.5正应力-297.2-196.06-93.250.0093196.08297.23表格 4考查点的x/m01.53.04.56.07.59.0方向位移(10-6)-0.3481-0.3376-0.3065-0.2568-0.1910-0.112505.4. 结果分析比较为了更直观的比较自编程序和ANSYS的计算结果，将Error!

14、Reference source not found.和Error! Reference source not found.的数据绘入Error! Reference source not found.，将Error! Reference source not found.和Error! Reference source not found.的数据绘入Error! Reference source not found.。图 56 应力比较图 57 位移比较 自编程序所得结果与ANSYS分析结果进行比较发现，应力偏大而位移偏小。分析其中的原因，一方面是编程水平有限，自编程序有很多不完善之处，有很

15、多因素没有考虑（温度、变形、非线性等），所得结果可信度不是很高，只能得到应力以及位移大概的分布情况；另一方面是在后处理阶段，在对高斯积分点结果进行处理时，误差难以避免，还有一方面的原因是在进行求解时保留数据精度不够，造成误差较大，并且进行数值运算时可能会造成大数吃小数，从而引起结果的差异。 参考文献1 (美)史密斯(Smith,I.m.)著;王崧等译.有限单元法编程(第三版)M.北京:电子工业出版社,20032 王勖成.有限单元法M.北京:清华大学出版社,20033 宋建辉,涂志刚.MATLAB语言及其在有限元编程中的应用J.湛江师范学院学报.2003.6(24):101-1054 郑大玉,

16、连宏玉, 周跃发. 有限元编程与C语言J. 黑龙江商学院学报.1997.3(13):23-285 王伟,刘德富.有限元编程中应用面向对象编程技术的探讨J.三峡大学学报.2001.2(23):124-1286 汪文立, 吕士俊.二级C语言程序设计教程M. 北京:中国水利水电出版社,20067 赵翔龙.FORTRAN 90学习教程M.北京:北京大学出版社,2002附录 程序源代码#include<stdio.h>#include<math.h>/*The gemotry of the model*/ float mesh2; float wide,hight; float

17、wsize,hsize,young,poiss,thick; int i,j,k,knelem,npoin,elem500,ielem; float coord21,props31,lnods8500,ifpre1; float posgp3,weigp3,estif136,elcod28; int npoin,nelem,kevab,igaus,jgaus,kgasp,ngash,djacb; float gpcod29,bmatx316,dmatx33; float shape8,deriv28; float xjacm22,xjaci22,elcod28; float cartd28;

18、float bmatx316,dmatx33; float nload1; float press32,pgash2,dgash2;struct node int nodenu;/*The number of the node*/ float cor2;/*The coordinate of the node*/ float displacement2;/*The displacement of the node*/ float load2; /*The load of the node*/ float boundary2; node500; struct int elementnu;/*Th

19、e number of element*/ int localnu8;/*Local number*/ int localcorx8; int localcory8; /*Local coordinate of node*/ float qx,qy,t;/*The stress and strain*/ element500;void meshing(float wide,float hight,float mesh2) printf("Please input the meshing sizen"); scanf("%f%f",&wsize,&

20、amp;hsize); getchar(); mesh0=wide/wsize; mesh1=hight/hsize; /*The global coordinate of node*/ void coordinate() int i,j=1; float x,y; for(i=1;i<=2*mesh1+1;i+) x=0; while(x<=wide) if(i%2!=0) nodej.cor1+=wsize/2; nodej.cor2+=(i-1)*hsize; j+; else nodej.cor1+=wsize; nodej.cor2+=(i-1)*hsize; j+; m

21、ain() int top500,bottom500;/*The number of top and bottom element*/ void input(); void stif(); void jacobi2( ); void load(); void asem(); void solve( ); void stre(); npoin=8; for(i=1;i<=8;i+) lnodsi1=i;for(i=1;i<=3;i+) scanf("%f",&propsi1); printf("The EX is %enThe uo is %8f

22、nThe bt is %8f",props11,props21,props31); getch(); printf("Please input the gemotry of the modeln"); scanf("%f%f",&wide,&hight); getchar(); printf("The wide is %0.3f,the hight is %0.3f",wide,hight); getchar(); meshing(wide,hight,mesh); printf("The wide

23、 size is %f,the hight size is %fn",mesh0,mesh1); input(); getchar(); nelem=mesh0*mesh1; for(i=1;i<=nelem;i+) /*The element number*/ elementi.elementnu=i; for(i=1;i<mesh0;i+) npoin+=5; for(i=1;i<mesh1;i+) npoin+=3*mesh0+2; printf("%d",npoin); for(i=1;i<=npoin;i+) nodei.node

24、nu=i; for(i=1;i<=nelem;i+) for(j=1;j<=8;j+) elementi.localnuj=j; for (i=1;i<=mesh0;i+) bottomi=i; j=1; for(i=mesh0*(mesh1-1)+1;i<=nelem;i+) topj+=i; for(i=1;i<j;i+) printf("the top number is %dn",topi); printf("%dn",element1.localnu8); coordinate(); stif(); jacobi2

25、( ); load(); asem(); solve( ); stre();getchar();/*data input*/void input() int nnode=8; int notreadchar;/*input element node number*/printf("input element node numbern");for(ielem=1;ielem<=nelem;ielem+) for(i=1;i<=nnode;i+) scanf("%d",&lnodsiielem);/*Input restriction m

26、essage*/printf("double-digit is symbol of the restriction,1 representates stable,2 representates gived displacementn"); i=1; while(i) scanf("%d",&i); if(i!=0) scanf("%d",&j); scanf("%d",&ifprej); else break; /*Element stiffness matrix*/void stif()

27、int kevab,nstre,ievab,istre,lnode,jstre,jevab; float kgasp,exisp,etasp, dvolu; float btdbm,dbmat3; void sfr();/*Gauss point*/ posgp1=-sqrt(0.6); posgp2=0; posgp3=sqrt(0.6); /*weight coefficient*/ weigp1=5.0/9.0; weigp2=8.0/9.0; weigp3=5.0/9.0; /*number of stress*/ nstre=3; /*form elastic matrix*/ fo

28、r(ielem=1;ielem<=nelem;ielem+) printf("The number of element is %dn",ielem);for(i=1;i<=8;i+) lnode=lnodsiielem; for(j=1;j<=2;j+) coordjlnode=nodelnode.corj; elcodji=coordjlnode; young=props11; poiss=props21; thick=props31; modps(young,poiss); /*The initial value of k*/ kevab=0; fo

29、r(i=1;i<=16;i+) for(j=1;j<=i;j+) kevab+; estifkevab=0.0; /*The gauss point shape function and deriv*/ kgasp=0; for(igaus=1;igaus<=3;igaus+) for(jgaus=1;jgaus<=3;jgaus+) kgasp=kgasp+1; exisp=posgpigaus; etasp=posgpjgaus; printf("The number of gauss point is %dn",kgasp); sfr2(exi

30、sp,etasp); jacob2(ielem,djacb,kgasp,elcod); dvolu=djacb*weigpigaus*weigpjgaus*thick; bmatps()/*The initial value of elastic matrix D*/ kevab=0; for(ievab=1;ievab<=16;ievab+) for(istre=1;istre<=nstre;istre+) dbmatistre=0.0; for(jstre=1;jstre<=nstre;jstre+) dbmatistre=dbmatistre+bmatxjstreiev

31、ab*dmatxjstreistre; for(jevab=1;jevab<=ievab;ievab+) kevab=kevab+1; btdbm=0.0; for(istre=1;istre<=nstre;istre+) btdbm=btdbm+dbmatistre*bmatxistrejevab; estifkevab=estifkevab+btdbm*dvolu; /*float sfr2(float s,float t) float s2,t2,ss,tt,st,stt,sst,st2;/*Shape function*/*these variables are to si

32、mplify the formula */ 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; /*shape function*/ shape1=(-1.0+st+ss+tt-sst-stt)/4.0; shape2=(1.0-t-ss+sst)/2.0; shape3=(-1.0-st+ss+tt-sst+stt)/4.0; shape4=(1.0+s-tt-stt)/2.0; shape5=(-1.0+st+ss+tt+sst+stt)/4.0; shape6=(1.0+t-ss-s

33、st)/2.0; shape7=(-1.0-st+ss+tt+sst-stt)/4.0; shape8=(1.0-s-tt+stt)/2.0; /* differential coefficient to local coordinate*/ deriv11=(t+s2-st2-tt)/4.0; deriv12=-s+st; deriv13=(-t+s2-st2+tt)/4.0; deriv14=(1.0-tt)/2.0; deriv15=(t+s2+st2+tt)/4.0; deriv16=-s-st; deriv17=(-t+s2+st2-tt)/4.0; deriv18=(-1.0+tt

34、)/2.0; deriv21=(s+t2-ss-st2)/4.0; deriv22=(-1.0+ss)/2.0; deriv23=(-s+t2-ss+st2)/4.0; deriv24=-t-st; deriv25=(s+t2+ss+st2)/4.0; deriv26=(1.0-ss)/2.0; deriv27=(-s+t2+ss-st2)/4.0; deriv28=-t+st; */*Jacobi matrix*/void jacobi2( ) /* xjacm22:jacobi matrix;xjaci22:j-1;elcod28:local coordinate*/ float cart

35、d28,gpcod29; int i,j,k; for(i=1;i<=2;i+) gpcodikgasp=0.0; for(j=1;j<=8;j+) gpcodikgasp=gpcodikgasp+elcodij*shapej; for(i=1;i<=2;i+) for(j=1;j<=2;j+) xjacmij=0.0; for(k=1;k<=8;k+) xjacmij=xjacmij+derivik*elcodjk; /*Caculate the value of Jacobi matrix*/ djacb=xjacm11*xjacm22-xjacm12*xja

36、cm21; printf("The det of jacobi matrix is %fn",djacb); if(djacb<=1e-6) printf("The jacobi det of %d is less or equal 0n",ielem); /*The element of j-1*/ xjaci11=xjacm22/djacb; xjaci22=xjacm11/djacb; xjaci12=-xjacm12/djacb; xjaci21=-xjacm21/djacb; /*The deriv of shape function t

37、o global coordinate*/ for(i=1;i<=2;i+) for(k=1;k<=8;k+) cartdik=0.0; for(j=1;j<=2;j+) cartdik=cartdik+xjaciij*derivjk; /* Form elastic matris */*float modps() float bmatx316,dmatx33; float constant; int i,j,y,p; y=props11; p=props21; for(i=1;i<=3;i+) for(j=1;j<=3;j+) dmatxij=0.0; /*If

38、 Ntype=1,it is plane stress question*/constant=y/(1.0-p*p); dmatx11=constant; dmatx22=constant; dmatx12=constant*p; dmatx21=constant*p; dmatx33=constant*(1.0-p)/2.0; */*Form strain matrix*/*void bmatps() float cartd28,shape8,deriv28,bmatx316,dmatx33; int npoin,nelem,ngash,mgash; float gpcod29;/*The

39、initial value ofB*/ for(i=1;i<=16;i+) for(j=1;j<=3;j+) bmatxji=0.0; /*Form B*/ ngash=0; for(i=1;i<=8;i+) mgash=ngash+1; ngash=mgash+1; bmatx1mgash=cartd1i; bmatx1ngash=0.0; bmatx2mgash=0.0; bmatx2ngash=cartd2i; bmatx3mgash=cartd2i; bmatx3ngash=cartd1i; */ /*equal load of node*/void load()fl

40、oat nload500; float elcod23,eload16,posgp3,weigp3; float press32,pgash2,dgash2,noprs3; int iedge,nedge,kount; float sgmar,tau,s,dvolu,pxcom,pycom,n,m,t; int lnode,number,nloca; /*The number of load side*/ /*element load*/ /*The position of gauss point*/ printf("input the number of load side*/n&

41、quot;); scanf("%d",&number); posgp1=-sqrt(0.6); posgp2=0.0; posgp3=sqrt(0.6); /*The weight of gauss point*/ weigp1=5.0/9.0; weigp2=8.0/9.0; weigp3=5.0/9.0; /*The initial load of element*/ for(i=1;i<=nelem;i+) nloadi=0.0; for(i=1;i<=number;i+) /*The initial value of equal node loa

42、d*/ for(i=1;i<=16;i+) eloadi=0; scanf("%d",&nedge); for(iedge=1;iedge<=nedge;iedge+) for(i=1;i<=3;i+) scanf("%f%f",&sgmar,&tau); /*if sgmar=0,input node load q(1,2,3)and t(1,2,3)*/ if(fabs(sgmar)<1e-8)&&(fabs(tau)<1e-8) for(j=1;j<=3;j+) for(i

43、=1;i<=2;i+) scanf("%f",&pressji); else for(i=1;i<=3;i+) pressi1=sgmar; pressi2=tau; /*The coordinate of the load node*/ for(i=1;i<=3;i+) /*input load node*/ printf("input node load numbern"); scanf("%d",&noprsi); lnode=noprsi; for(j=1;j<=2;j+) coordj

44、lnode=nodelnode.corj; elcodji=coordjlnode; t=-1.0; for(igaus=1;igaus<=3;igaus+) s=posgpigaus; sfr(s,t); for(j=1;j<=2;j+) pgashj=0.0; dgashj=0.0; /*current point q,t and the value of deriv*/ pgashj=pgashj+pressij*shapei; dgashj=dgashj+elcodji*deriv1i; dvolu=weigpigaus; pxcom=dgash1*pgash2-dgash2*pgash1; pycom=dgash1*pgash1+dgash2*pgash2; for(i=1;i<=8;i+) nloca=lnodsiielem; if(nloca=noprs1) 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; /*if the local node numner >8,then it is 1*/ if(k>8)