八节点等参单元平面有限元Word格式.docx
- 文档编号:17517027
- 上传时间:2022-12-06
- 格式:DOCX
- 页数:38
- 大小:237.33KB
八节点等参单元平面有限元Word格式.docx
《八节点等参单元平面有限元Word格式.docx》由会员分享,可在线阅读,更多相关《八节点等参单元平面有限元Word格式.docx(38页珍藏版)》请在冰豆网上搜索。
由于
是
、
的函数,在
中的
要按照复合函数来求导,即
从而有
因此,单元应力矩阵为
单元刚度矩阵为
其中积分采用三点高斯积分,
(高斯积分点的总数),
和
或
是加权系数,
是单元内的坐标.。
对于三点高斯积分,高斯积分点的位置:
,
。
单元等效节点荷载
结构刚度矩阵
结构结点荷载列阵
注意,对于式和式中
的理解不是简单的叠加而是集成。
总刚平衡方程
从式求出
将
回代入式和式,得到
2.1.
2.2.有限元分析的模块组织
一个典型的有限元分析过程主要包括以下几个步骤:
1)读输入数据,定义节点及单元数组。
2)由边界条件计算方程个数,赋值荷载列阵。
3)读入在带状存储的总刚度矩阵中单元和载荷信息。
4)定义总刚度阵数组。
5)组装总刚度阵。
6)解方程组。
7)
图22
其流程图可见Error!
Referencesourcenotfound.。
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]
等效节点荷载
高斯积分点的总体坐标
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
Referencesourcenotfound.和Error!
Referencesourcenotfound.分别给出了结构的位移图和
应力云图。
图54位移图
图55各单元
图
从Error!
Referencesourcenotfound.可以看到,跨中的位移和正应力最大,与简支梁承受均布荷载,跨中挠度和正应力最大的力学常识相符合,可以初步认为模型是正确的。
Referencesourcenotfound.给出了
的截面上的正应力
和Error!
处各横截面
方向位移,其中
的单位为
表格1
考查点的y/m
正应力
表格2
考查点的x/m
方向位移(10-6)
5.3.自编程序分析结果
用自编程序进行分析时,采用了整个结构模型进行计算,其他条件不变,因编程水平有限,在后处理阶段没有给出节点处的位移与应力,而只能得到高斯积分点处的位移与应力,得到积分点处的应力后,根据数值分析相关知识采用了插值进行处理,从而得到
和
方向位移,分别列出表格如下:
表格3
93
表格4
5.4.结果分析比较
为了更直观的比较自编程序和ANSYS的计算结果,将Error!
Referencesourcenotfound.的数据绘入Error!
Referencesourcenotfound.,将Error!
图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<
>
/*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;
=8;
lnods[i][1]=i;
for(i=1;
=3;
scanf("
%f"
props[i][1]);
printf("
TheEXis%e\nTheuois%8f\nThebtis%8f"
props[1][1],props[2][1],props[3][1]);
getch();
Pleaseinputthegemotryofthemodel\n"
wide,&
hight);
getchar();
Thewideis%,thehightis%"
wide,hight);
meshing(wide,hight,mesh);
Thewidesizeis%f,thehightsizeis%f\n"
mesh[0],mesh[1]);
input();
nelem=mesh[0]*mesh[1];
=nelem;
i++)/*Theelementnumber*/
element[i].elementnu=i;
mesh[0];
npoin+=5;
mesh[1];
npoin+=3*mesh[0]+2;
%d"
npoin);
=npoin;
node[i].nodenu=i;
for(j=1;
j<
j++)
element[i].localnu[j]=j;
for(i=1;
=mesh[0];
bottom[i]=i;
j=1;
for(i=mesh[0]*(mesh[1]-1)+1;
top[j++]=i;
j;
thetopnumberis%d\n"
top[i]);
%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<
ielem++)
=nnode;
lnods[i][ielem]);
/*Inputrestrictionmessage*/
double-digitissymboloftherestriction,1representatesstable,2representatesgiveddisplacement\n"
i=1;
while(i)
{scanf("
i);
if(i!
{scanf("
j);
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;
{printf("
Thenumberofelementis%d\n"
ielem);
{lnode=lnods[i][ielem];
=2;
{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;
=16;
=i;
{kevab++;
estif[kevab]=;
/*Thegausspointshapefunctionandderiv*/
kgasp=0;
for(igaus=1;
igaus<
igaus++)
for(jgaus=1;
jgaus<
jgaus++)
{kgasp=kgasp+1;
exisp=posgp[igaus];
etasp=posgp[jgaus];
Thenumberofgausspointis%d\n"
kgasp);
sfr2(exisp,etasp);
jacob2(ielem,djacb,kgasp,elcod);
dvolu=djacb*weigp[igaus]*weigp[jgaus]*thick;
bmatps()
/*Theinitialvalueofelasticmatrix[D]*/
for(ievab=1;
ievab<
ievab++)
{for(istre=1;
istre<
=nstre;
istre++)
{dbmat[istre]=;
for(jstre=1;
jstre<
jstre++)
dbmat[istre]=dbmat[istre]+bmatx[jstre][ievab]*dmatx[jstre][istre];
for(jevab=1;
jevab<
=ievab;
{kevab=kevab+1;
btdbm=;
for(istre=1;
btdbm=btdbm+dbmat[istre]*bmatx[istre][jevab];
estif[kevab]=estif[k
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 八节 单元 平面 有限元