八节点等参单元平面有限元Word文件下载.docx
《八节点等参单元平面有限元Word文件下载.docx》由会员分享，可在线阅读，更多相关《八节点等参单元平面有限元Word文件下载.docx（40页珍藏版）》请在冰豆网上搜索。
单元刚度矩阵为
其中积分采用三点高斯积分，
（高斯积分点的总数）,
和
或
是加权系数，
是单元内的坐标。
。
对于三点高斯积分,高斯积分点的位置：
，
.
单元等效节点荷载
结构刚度矩阵
结构结点荷载列阵
注意，对于式和式中
的理解不是简单的叠加而是集成。
总刚平衡方程
从式求出
将
回代入式和式，得到
2.1.
2.2.有限元分析的模块组织
一个典型的有限元分析过程主要包括以下几个步骤：
1）读输入数据,定义节点及单元数组。
2）由边界条件计算方程个数，赋值荷载列阵。
3）读入在带状存储的总刚度矩阵中单元和载荷信息.
4）定义总刚度阵数组。
5）组装总刚度阵。
6）解方程组。
其流程图可见图22.
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.问题说明
图51所示一简支梁，高3m,长18m,承受均布荷载
，取
m，作为平面应力问题.
图52
由于结构和荷载对称,只对右边一半进行有限单元计算,如图52所示，而在y轴各节点布置水平连杆支座。
5.2.ANSYS分析结果
ANSYS分析中,采用plane82单元，在板单元上边施加均布荷载
，在y轴上的各结点约束UX方向的自由度，在板单元右下角的结点约束UX自由度,结点布置、网格划分方案如图53所示。
图53
图54位移图和图55各单元
图分别给出了结构的位移图和
应力云图。
图54位移图
图55各单元
图
从图54位移图和图55各单元
图可以看到,跨中的位移和正应力最大,与简支梁承受均布荷载，跨中挠度和正应力最大的力学常识相符合，可以初步认为模型是正确的。
表格1给出了
的截面上的正应力
和表格2给出了
处各横截面
方向位移，其中
的单位为
表格1
考查点的y/m
1.5
1。
0.5
—0。
5
-1.0
—1。
正应力
-270.20
—178。
25
-88.57
0.00
88.57
178。
270.21
表格2
考查点的x/m
3。
4。
6。
7。
9.0
方向位移（10—6）
-0。
3485
3380
3069
—0.2571
—0.1914
-0.1133
5.3.自编程序分析结果
用自编程序进行分析时，采用了整个结构模型进行计算,其他条件不变，因编程水平有限,在后处理阶段没有给出节点处的位移与应力，而只能得到高斯积分点处的位移与应力，得到积分点处的应力后，根据数值分析相关知识采用了插值进行处理，从而得到
和
方向位移，分别列出表格如下：
表格3
1.0
-0.5
—1.0
—297。
—196.06
-93.25
93
196。
08
297.23
表格4
3.0
7.5
9。
3481
3376
3065
-0.2568
-0.1910
1125
5.4.结果分析比较
为了更直观的比较自编程序和ANSYS的计算结果，将表格1和表格3的数据绘入图56，将表格2和表格4的数据绘入图57。
图56应力比较
图57位移比较
自编程序所得结果与ANSYS分析结果进行比较发现，应力偏大而位移偏小.分析其中的原因，一方面是编程水平有限,自编程序有很多不完善之处,有很多因素没有考虑（温度、变形、非线性等），所得结果可信度不是很高，只能得到应力以及位移大概的分布情况；
另一方面是在后处理阶段,在对高斯积分点结果进行处理时,误差难以避免，还有一方面的原因是在进行求解时保留数据精度不够，造成误差较大，并且进行数值运算时可能会造成大数吃小数,从而引起结果的差异.
参考文献
[1]（美）史密斯（Smith，I.m.）著;
王崧等译。
有限单元法编程（第三版）［M］.北京:
电子工业出版社,2003
[2]王勖成.有限单元法[M］.北京:
清华大学出版社，2003
[3]宋建辉,涂志刚。
MATLAB语言及其在有限元编程中的应用［J]。
湛江师范学院学报。
2003。
6（24）：
101-105
[4]郑大玉,连宏玉，周跃发。
有限元编程与C语言[J].黑龙江商学院学报。
1997.3（13）:
23-28
[5]王伟,刘德富.有限元编程中应用面向对象编程技术的探讨［J].三峡大学学报.2001.2（23）：
124-128
[6]汪文立，吕士俊。
二级C语言程序设计教程[M］。
北京：
中国水利水电出版社，2006
[7]赵翔龙。
FORTRAN90学习教程[M］。
北京:
北京大学出版社，2002
附录程序源代码
＃include〈stdio.h>
＃include〈math.h>
/＊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;
i〈=3；
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（）；
Thewideis％0.3f，thehightis%0。
3f”,wide,hight）；
meshing（wide,hight,mesh）;
printf（”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；
for（i=1;
mesh［1］;
npoin+=3＊mesh[0］+2；
printf（”%d"
，npoin）；
i〈=npoin;
node［i］.nodenu=i；
i〈=nelem;
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；
i〈j;
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++）
=nnode；
scanf（”％d"
lnods[i］［ielem］）;
/＊Inputrestrictionmessage＊/
printf（"
double-digitissymboloftherestriction，1representatesstable，2representatesgiveddisplacement\n"
i=1;
while（i）
{scanf（"
%d"
i）;
if（i！
｛scanf（"
％d”，＆j）;
%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（0.6）;
posgp[2]=0;
posgp[3］=sqrt（0.6）;
/*weightcoefficient*/
weigp［1］=5.0/9.0；
weigp[2]=8.0/9。
0；
weigp[3]=5.0/9。
/*numberofstress＊/
nstre=3;
/*formelasticmatrix＊/
for（ielem=1;
{printf（”Thenumberofelementis%d\n"
，ielem）;
for（i=1；
i〈=8；
{lnode=lnods［i]［ielem]；
for（j=1；
j〈=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；
i〈=16；
j〈=i；
{kevab++；
estif[kevab]=0.0；
/＊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]*/
for（ievab=1；
ievab〈=16;
ievab++）
｛for（istre=1;
istre〈=nstre;
istre++）
｛dbmat［istre］=0.0；
for（jstre=1;
jstre<
=nstre;
jstre++）
dbmat［istre］=dbmat[istre］+bmatx[jstre］［ievab］＊dmatx[jstre][istre］;
for（jevab=1;
jevab<
=ievab；
{kevab=kevab+1；
btdbm=0。
0;
for（istre=1；
istre<
=nstre；
btdbm=btdbm+dbmat［istre］＊bmatx[istre]［jevab］;
estif[kevab］=estif［kevab］+btdbm*dvolu;
｝}
/*floatsfr2（floats,floatt）
{
floats2，t2，ss，tt,st,stt，sst,st2;
/*Shapefunction*/
/＊thesevariablesaretosimplifythefor
升级会员 