八节点等参单元平面有限元.docx
- 文档编号:4582086
- 上传时间:2022-12-06
- 格式:DOCX
- 页数:38
- 大小:237.33KB
八节点等参单元平面有限元.docx
《八节点等参单元平面有限元.docx》由会员分享,可在线阅读,更多相关《八节点等参单元平面有限元.docx(38页珍藏版)》请在冰豆网上搜索。
八节点等参单元平面有限元
八节点等参单元平面有限元
分析程序
土木工程学院
《计算力学》课程大作业
1.概述
通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、SAP等)进行前处理和后处理,而中间的计算部分一般采用自己编制的程序来运算。
具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。
随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中,但在实际应用中也暴露出一些问题。
有时网格离散化的区域较大,而又限于研究精度的要求,使得划分的网格数目极其庞大,结点数可多达数万个,从而造成计算中要运算的数据量巨大,程序运行的时间较长的弊端,这就延长了问题解决的时间,使得求解效率降低。
因为运行周期长,不利于程序的调试,特别是对于要计算多种运行工况时的情况;同时大数据量处理对计算机的内存和CPU提出了更高的要求,而在实际应用中,单靠计算机硬件水平的提高来解决问题的能力是有限的。
因此,必须寻找新的编程语言。
随着有限元前后处理的不断发展和完善,以及大型工程分析软件对有限元接口的要求,有限元分析程序不应只满足解题功能,它还应满足软件工程所要求的结构化程序设计条件,能够对存储进行动态分配,以充分利用计算机资源,它还应很容易地与其它软件如CAD的实体造型,优化设计等接口。
现在可编写工程应用软件的计算机语言较多,其中C语言是一个较为优秀的语言,很容易满足现在有限元分析程序编程的要求。
C语言最初是为操作系统、编译器以及文字处理等编程而发明的。
随着不断完善,它已应用到其它领域,包括工程应用软件的编程。
近年来,C语言已经成为计算机领域最普及的一个编程语言,几乎世界上所有的计算机都装有C的编译器,从PC机到巨型机到超巨型的并行机,C与所有的硬件和操作系统联系在一起。
用C编写的程序,可移植性极好,几乎不用作多少修改,就可在任何一台装有ANSI、C编译器的计算机上运行。
C既是高级语言,也是低级语言,也就是说,可用它作数值计算,也可用它对计算机存储进行操作。
2.编程思想
本程序采用C语言编程,编制平面四边形八节点等参元程序,用以求解平面结构问题。
程序采用二维等带宽存储整体刚度矩阵,乘大数法引入约束,等带宽高斯消去法求解位移。
在有限元程序中,变量数据需赋值的可分为节点信息,单元信息,载荷信息等。
对于一个节点来说,需以下信息:
节点编号(整型),节点坐标(实型),节点已知位移(实型),节点载荷(实型),边界条件(实型)和节点温度(实型)等。
同样,对于一个单元来说,需以下信息:
单元的节点联接信息(整型),单元类型信息(桁架、梁、板、壳等)(整型),单元特性信息(厚度、内力矩等)(实型),材料信息(弹性模量,泊松比等)(实型)等。
在FORTRAN程序中,以上这些变量混合在一起,很难辨认,使程序的可读性不好,如需要进行单元网络的自适应划分,节点及单元的修改将非常困难。
在进行C语言编译过程中,采用结构struct使每个节点信息存储在一个结构体数组中,提高程序的可读性,使数据结构更趋于合理。
1.
2.
2.1.八节点矩形单元介绍
八节点矩形单元编号如Error!
Referencesourcenotfound.所示
图21
八节点矩形单元的位移函数为:
其形函数为
式和式中
,并且采用等参变换,则单元的坐标变换式可取为
单元应变矩阵为
式一般简写为
其中
的子块矩阵为
由于
是
、
的函数,在
中的
、
要按照复合函数来求导,即
从而有
因此,单元应力矩阵为
单元刚度矩阵为
其中积分采用三点高斯积分,
(高斯积分点的总数),
和
或
是加权系数,
和
是单元内的坐标.。
对于三点高斯积分,高斯积分点的位置:
,
,
。
单元等效节点荷载
结构刚度矩阵
结构结点荷载列阵
注意,对于式和式中
的理解不是简单的叠加而是集成。
总刚平衡方程
从式求出
将
回代入式和式,得到
和
。
1.
2.
2.1.
2.2.有限元分析的模块组织
一个典型的有限元分析过程主要包括以下几个步骤:
1)读输入数据,定义节点及单元数组。
2)由边界条件计算方程个数,赋值荷载列阵。
3)读入在带状存储的总刚度矩阵中单元和载荷信息。
4)定义总刚度阵数组。
5)组装总刚度阵。
6)解方程组。
7)
图22
其流程图可见Error!
Referencesourcenotfound.。
3.程序变量及函数说明
3.
3.
3.1.主要变量说明:
主要程序变量说明
wide
分析模型的宽
hight
分析模型的高
wsize
宽度方向网格划分尺寸
hsize
高度方向网格划分尺寸
npoin
节点总数
nelem
单元总数
structnode[500]
节点结构体变量
structelement[500]
单元结构体变量
ifpre[500]
节点约束信息
posgp[3]
高斯积分点
weigp[3]
权系数
gpcod[2][9]
高斯积分点总体坐标
bmatx[3][16]
单元变形矩阵
dmatx[3][3]
单元弹性矩阵
xjacm[2][2]
雅可比矩阵
xjaci[2][2]
雅可比矩阵的逆矩阵
djacb
雅可比矩阵行列式的值
estif[136]
单元刚度矩阵
shape[8]
单元形函数
deriv[2][8]
形函数对局部坐标的导数
cartd[2][8]
形函数对整体坐标的导数
A[30000]
总体刚度矩阵
eload[16]
等效节点荷载
gpcod[2][9]
高斯积分点的总体坐标
3.1.
3.2.主要函数说明
主要函数说明
voidmeshing()
网格划分
voidcoordinate()
节点整体坐标计算
voidinput()
信息输入
voidstif()
形成单元刚度矩阵
voidsfr2()
计算形函数对当前积分点及形函数对局部坐标的到数值
voidjacobi2()
形成雅可比矩阵
voidmodps()
形成弹性矩阵
voidbmatps()
形成应变矩阵
voidload()
计算等效节点荷载
voidasem()
形成整体刚度矩阵和整体荷载列阵
voidsolve()
求解整体方程
voidstre()
计算单元应力
4.程序流程图
任何一个有限元程序处理过程,都可以划分为三个阶段:
1)前处理:
读入数据和生成数据。
2)处理器:
有限元的矩阵计算、组集和求解。
3)后处理:
输出节点位移、单元应变等。
为了更清楚地描述程序,我们给出了程序的流程图。
5.程序应用与ANSYS分析的比较
为了验证程序的正确性,我们取了一个简单框架结构,分别用自编程序和ANSYS进行分析和比较。
4.
5.
5.1.问题说明
Error!
Referencesourcenotfound.所示一简支梁,高3m,长18m,承受均布荷载
,
,
,取
m,作为平面应力问题。
图52
由于结构和荷载对称,只对右边一半进行有限单元计算,如Error!
Referencesourcenotfound.所示,而在y轴各节点布置水平连杆支座。
5.2.ANSYS分析结果
ANSYS分析中,采用plane82单元,在板单元上边施加均布荷载
,在y轴上的各结点约束UX方向的自由度,在板单元右下角的结点约束UX自由度,结点布置、网格划分方案如Error!
Referencesourcenotfound.所示。
图53
Error!
Referencesourcenotfound.和Error!
Referencesourcenotfound.分别给出了结构的位移图和
应力云图。
图54位移图
图55各单元
图
从Error!
Referencesourcenotfound.和Error!
Referencesourcenotfound.可以看到,跨中的位移和正应力最大,与简支梁承受均布荷载,跨中挠度和正应力最大的力学常识相符合,可以初步认为模型是正确的。
Error!
Referencesourcenotfound.给出了
的截面上的正应力
和Error!
Referencesourcenotfound.给出了
处各横截面
方向位移,其中
的单位为
,
的单位为
。
表格1
考查点的y/m
0
正应力
表格2
考查点的x/m
0
方向位移(10-6)
0
5.3.自编程序分析结果
用自编程序进行分析时,采用了整个结构模型进行计算,其他条件不变,因编程水平有限,在后处理阶段没有给出节点处的位移与应力,而只能得到高斯积分点处的位移与应力,得到积分点处的应力后,根据数值分析相关知识采用了插值进行处理,从而得到
的截面上的正应力
和
处各横截面
方向位移,分别列出表格如下:
表格3
考查点的y/m
0
正应力
93
表格4
考查点的x/m
0
方向位移(10-6)
0
5.4.结果分析比较
为了更直观的比较自编程序和ANSYS的计算结果,将Error!
Referencesourcenotfound.和Error!
Referencesourcenotfound.的数据绘入Error!
Referencesourcenotfound.,将Error!
Referencesourcenotfound.和Error!
Referencesourcenotfound.的数据绘入Error!
Referencesourcenotfound.。
图56应力比较
图57位移比较
自编程序所得结果与ANSYS分析结果进行比较发现,应力偏大而位移偏小。
分析其中的原因,一方面是编程水平有限,自编程序有很多不完善之处,有很多因素没有考虑(温度、变形、非线性等),所得结果可信度不是很高,只能得到应力以及位移大概的分布情况;另一方面是在后处理阶段,在对高斯积分点结果进行处理时,误差难以避免,还有一方面的原因是在进行求解时保留数据精度不够,造成误差较大,并且进行数值运算时可能会造成大数吃小数,从而引起结果的差异。
参考文献
[1](美)史密斯(Smith,.)著;王崧等译.有限单元法编程(第三版)[M].北京:
电子工业出版社,2003
[2]王勖成.有限单元法[M].北京:
清华大学出版社,2003
[3]宋建辉,涂志刚.MATLAB语言及其在有限元编程中的应用[J].湛江师范学院学报.(24):
101-105
[4]郑大玉,连宏玉,周跃发.有限元编程与C语言[J].黑龙江商学院学报.(13):
23-28
[5]王伟,刘德富.有限元编程中应用面向对象编程技术的探讨[J].三峡大学学报.(23):
124-128
[6]汪文立,吕士俊.二级C语言程序设计教程[M].北京:
中国水利水电出版社,2006
[7]赵翔龙.FORTRAN90学习教程[M].北京:
北京大学出版社,2002
附录程序源代码
#include<>
#include<>
/*Thegemotryofthemodel*/
floatmesh[2];
floatwide,hight;
floatwsize,hsize,young,poiss,thick;
inti,j,k,knelem,npoin,elem[500],ielem;
floatcoord[2][1],props[3][1],lnods[8][500],ifpre[1];
floatposgp[3],weigp[3],estif[136],elcod[2][8];
intnpoin,nelem,kevab,igaus,jgaus,kgasp,ngash,djacb;
floatgpcod[2][9],bmatx[3][16],dmatx[3][3];
floatshape[8],deriv[2][8];
floatxjacm[2][2],xjaci[2][2],elcod[2][8];
floatcartd[2][8];
floatbmatx[3][16],dmatx[3][3];
floatnload[1];
floatpress[3][2],pgash[2],dgash[2];
structnode
{intnodenu;/*Thenumberofthenode*/
floatcor[2];/*Thecoordinateofthenode*/
floatdisplacement[2];/*Thedisplacementofthenode*/
floatload[2];/*Theloadofthenode*/
floatboundary[2];
}node[500];
struct
{intelementnu;/*Thenumberofelement*/
intlocalnu[8];/*Localnumber*/
intlocalcorx[8];
intlocalcory[8];/*Localcoordinateofnode*/
floatqx,qy,t;/*Thestressandstrain*/
}element[500];
voidmeshing(floatwide,floathight,floatmesh[2])
{printf("Pleaseinputthemeshingsize\n");
scanf("%f%f",&wsize,&hsize);
getchar();
mesh[0]=wide/wsize;
mesh[1]=hight/hsize;
}
/*Theglobalcoordinateofnode*/
voidcoordinate()
{
inti,j=1;
floatx,y;
for(i=1;i<=2*mesh[1]+1;i++)
{x=0;
while(x<=wide)
if(i%2!
=0)
{node[j].cor[1]+=wsize/2;
node[j].cor[2]+=(i-1)*hsize;
j++;}
else
{node[j].cor[1]+=wsize;
node[j].cor[2]+=(i-1)*hsize;
j++;}
}
}
main()
{inttop[500],bottom[500];/*Thenumberoftopandbottomelement*/
voidinput();
voidstif();
voidjacobi2();
voidload();
voidasem();
voidsolve();
voidstre();
npoin=8;
for(i=1;i<=8;i++)
lnods[i][1]=i;
for(i=1;i<=3;i++)
scanf("%f",&props[i][1]);
printf("TheEXis%e\nTheuois%8f\nThebtis%8f",props[1][1],props[2][1],props[3][1]);
getch();
printf("Pleaseinputthegemotryofthemodel\n");
scanf("%f%f",&wide,&hight);
getchar();
printf("Thewideis%,thehightis%",wide,hight);
getchar();
meshing(wide,hight,mesh);
printf("Thewidesizeis%f,thehightsizeis%f\n",mesh[0],mesh[1]);
input();
getchar();
nelem=mesh[0]*mesh[1];
for(i=1;i<=nelem;i++)/*Theelementnumber*/
element[i].elementnu=i;
for(i=1;i npoin+=5; for(i=1;i npoin+=3*mesh[0]+2; printf("%d",npoin); for(i=1;i<=npoin;i++) node[i].nodenu=i; for(i=1;i<=nelem;i++) for(j=1;j<=8;j++) element[i].localnu[j]=j; for(i=1;i<=mesh[0];i++) bottom[i]=i; j=1; for(i=mesh[0]*(mesh[1]-1)+1;i<=nelem;i++) top[j++]=i; for(i=1;i printf("thetopnumberis%d\n",top[i]); printf("%d\n",element[1].localnu[8]); coordinate(); stif(); jacobi2(); load(); asem(); solve(); stre(); getchar(); } /*datainput*/ voidinput() {intnnode=8; intnotreadchar; /*inputelementnodenumber*/ printf("inputelementnodenumber\n"); for(ielem=1;ielem<=nelem;ielem++) for(i=1;i<=nnode;i++) scanf("%d",&lnods[i][ielem]); /*Inputrestrictionmessage*/ printf("double-digitissymboloftherestriction,1representatesstable,2representatesgiveddisplacement\n"); i=1; while(i) {scanf("%d",&i); if(i! =0) {scanf("%d",&j); scanf("%d",&ifpre[j]);} elsebreak; } } /*Elementstiffnessmatrix*/ voidstif() {intkevab,nstre,ievab,istre,lnode,jstre,jevab; floatkgasp,exisp,etasp,dvolu; floatbtdbm,dbmat[3]; voidsfr(); /*Gausspoint*/ posgp[1]=-sqrt; posgp[2]=0; posgp[3]=sqrt; /*weightcoefficient*/ weigp[1]=; weigp[2]=; weigp[3]=; /*numberofstress*/ nstre=3; /*formelasticmatrix*/ for(ielem=1;ielem<=nelem;ielem++) {printf("Thenumberofelementis%d\n",ielem); for(i=1;i<=8;i++) {lnode=lnods[i][ielem]; for(j=1;j<=2;j++) {coord[j][lnode]=node[lnode].cor[j]; elcod[j][i]=coord[j][lnode]; } } young=props[1][1]; poiss=props[2][1]; thick=props[3][1]; modps(young,poiss); /*Theinitialvalueof[k]*/ kevab=0; for(i=1;i<=16;i++) for(j=1;j<=i;j++) {kevab++; estif[kevab]=; } /*Thegausspointshapefunctionandderiv*/ kgasp=0; for(igaus=1;igaus<=3;igaus++) for(jgaus=1;jgaus<=3;jgaus++) {kgasp=kgasp+1; exisp=posgp[igaus]; etasp=posgp[jgaus]; printf("Thenumberofgausspointis%d\n",kgasp); sfr2(exisp,etasp); jacob2(ielem,djacb,kgasp,elcod); dvolu=djacb*weigp[igaus]*weigp[jgaus]*thick; bmatps() /*Theinitialvalueofelasticmatrix[D]*/ kevab=0; for(ievab=1;ievab<=16;ievab++) {for(istre=1;istre<=nstre;istre++) {dbmat[istre]=; for(jstre=1;jstre<=nstre;jstre++) dbmat[istre]=dbmat[istre]+bmatx[jstre][ievab]*dmatx[jstre][istre]; } for(jevab=1;jevab<=ievab;ievab++) {kevab=kevab+1; btdbm=; for(istre=1;istre<=nstre;istre++) btdbm=btdbm+dbmat[istre]*bmatx[istre][jevab]; estif[kevab]=estif[k
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 八节 单元 平面 有限元