《计算力学》课程大作业
八节点等参单元平面有限元
程序
土木
学院
2011.2
目录
11.
概述
22.
编程思想
22.1.
八节点矩形单元介绍
52.2.
有限元分析的模块组织
63.
程序变量及函数说明
63.1.
主要变量说明:
73.2.
主要函数说明
84.
程序流程图
95.
程序应用与ANSYS分析的比较
95.1.
问题说明
105.2.
ANSYS分析结果
125.3.
自编程序分析结果
125.4.
结果比较分析
14参考文献
15附录 程序源代码
《计算力学》课程大作业
1. 概述
通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、SAP等)进行前处理和后处理,而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中,但在实际应用中也暴露出一些问题。有时网格离散化的区域较大,而又限于研究精度的要求,使得划分的网格数目极其庞大,结点数可多达数万个,从而造成计算中要运算的数据量巨大,程序运行的时间较长的弊端,这就延长了问题解决的时间,使得求解效率降低。因为运行周期长,不利于程序的调试,特别是对于要计算多种运行工况时的情况;同时大数据量处理对计算机的内存和CPU 提出了更高的要求,而在实际应用中,单靠计算机硬件水平的提高来解决问题的能力是有限的。因此,必须寻找新的编程语言。
随着有限元前后处理的不断发展和完善,以及大型工程分析软件对有限元接口的要求,有限元分析程序不应只满足解题功能,它还应满足软件工程所要求的结构化程序
条件,能够对存储进行动态分配,以充分利用计算机资源,它还应很容易地与其它软件如CAD 的实体造型,优化设计等接口。现在可编写工程应用软件的计算机语言较多,其中C语言是一个较为优秀的语言,很容易满足现在有限元分析程序编程的要求。
C语言最初是为操作系统、编译器以及文字处理等编程而发明的。随着不断完善,它已应用到其它领域,包括工程应用软件的编程。近年来,C语言已经成为计算机领域最普及的一个编程语言,几乎世界上所有的计算机都装有C的编译器,从PC机到巨型机到超巨型的并行机,C与所有的硬件和操作系统联系在一起。用C 编写的程序,可移植性极好,几乎不用作多少修改,就可在任何一台装有ANSI、C编译器的计算机上运行。C既是高级语言,也是低级语言,也就是说,可用它作数值计算,也可用它对计算机存储进行操作。
2. 编程思想
本程序采用C语言编程,编制平面四边形八节点等参元程序,用以求解平面结构问题。程序采用二维等带宽存储整体刚度矩阵,乘大数法引入约束,等带宽高斯消去法求解位移。
在有限元程序中,变量数据需赋值的可分为节点信息,单元信息,载荷信息等。对于一个节点来说,需以下信息:节点编号(整型),节点坐标(实型),节点已知位移(实型),节点载荷(实型),边界条件(实型)和节点温度(实型)等。同样,对于一个单元来说,需以下信息:单元的节点联接信息(整型),单元类型信息(桁架、梁、板、壳等)(整型) ,单元特性信息(厚度、内力矩等)(实型),
信息(弹性模量,泊松比等)(实型)等。
在FORTRAN 程序中,以上这些变量混合在一起,很难辨认,使程序的可读性不好,如需要进行单元网络的自适应划分,节点及单元的修改将非常困难。在进行C语言编译过程中,采用结构struct 使每个节点信息存储在一个结构体数组中,提高程序的可读性,使数据结构更趋于合理。
1.1. 八节点矩形单元介绍
八节点矩形单元编号如图 2‑1所示
八节点矩形单元的位移函数为:
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.1)
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.2)
其形函数为
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.3)
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.4)
式(2.4)
中(2.3)
和式 GOTOBUTTON ZEqnNum727261 \* MERGEFORMAT ,并且采用等参变换,则单元的坐标变换式可取为
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.5)
单元应变矩阵为
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.6)
式(2.6)
一般简写为
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.7)
其中
的子块矩阵为
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.8)
由于
是
、
的函数,在
中的
、
要按照复合函数来求导,即
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.9)
从而有
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.10)
因此,单元应力矩阵为
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.11)
单元刚度矩阵为
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.12)
其中积分采用三点高斯积分,
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.13)
(高斯积分点的总数),
和
或
是加权系数,
和
是单元内的坐标.。对于三点高斯积分,高斯积分点的位置:
,
,
。
单元等效节点荷载
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.14)
结构刚度矩阵
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.15)
结构结点荷载列阵
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.16)
注意,对于式(2.16)
中(2.15)
和式 GOTOBUTTON ZEqnNum449883 \* MERGEFORMAT 的理解不是简单的叠加而是集成。
总刚平衡方程
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.17)
从式(2.17)
求出
MACROBUTTON MTPlaceRef \* MERGEFORMAT (2.18)
将
回代入式(2.11)
,得到(2.7)
和式 GOTOBUTTON ZEqnNum216712 \* MERGEFORMAT 和
。
1.2. 有限元分析的模块组织
一个典型的有限元分析过程主要包括以下几个步骤:
1) 读输入数据,定义节点及单元数组。
2) 由边界条件计算方程个数,赋值荷载列阵。
3) 读入在带状存储的总刚度矩阵中单元和载荷信息。
4) 定义总刚度阵数组。
5) 组装总刚度阵。
6) 解方程组。
其流程图可见图 2‑2。
3. 程序变量及函数说明
1.3. 主要变量说明:
主要程序变量说明
wide
分析模型的宽
hight
分析模型的高
wsize
宽度方向网格划分尺寸
hsize
高度方向网格划分尺寸
npoin
节点总数
nelem
单元总数
struct node[500]
节点结构体变量
struct element [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]
高斯积分点的总体坐标
1.4. 主要函数说明
主要函数说明
void meshing( )
网格划分
void coordinate( )
节点整体坐标计算
void input( )
信息输入
void stif( )
形成单元刚度矩阵
void sfr2()
计算形函数对当前积分点及形函数对局部坐标的到数值
void jacobi2( )
形成雅可比矩阵
void modps( )
形成弹性矩阵
void bmatps( )
形成应变矩阵
void load( )
计算等效节点荷载
void asem( )
形成整体刚度矩阵和整体荷载列阵
void solve( )
求解整体方程
void stre( )
计算单元应力
4. 程序流程图
任何一个有限元程序处理过程,都可以划分为三个阶段:
1) 前处理:读入数据和生成数据。
2) 处理器:有限元的矩阵计算、组集和求解。
3) 后处理:输出节点位移、单元应变等。
为了更清楚地描述程序,我们给出了程序的流程图。
5. 程序应用与ANSYS分析的比较
为了验证程序的正确性,我们取了一个简单框架结构,分别用自编程序和ANSYS进行分析和比较。
1.5. 问题说明
图 5‑1所示一简支梁,高3m,长18m,承受均布荷载
,
,
,取
m,作为平面应力问题。
由于结构和荷载对称,只对右边一半进行有限单元计算,如图 5‑2所示,而在y轴各节点布置水平连杆支座。
1.6. ANSYS分析结果
ANSYS分析中,采用plane82单元,在板单元上边施加均布荷载
,在y轴上的各结点约束UX方向的自由度,在板单元右下角的结点约束UX自由度,结点布置、网格划分
如图 5‑3所示。
图 5‑3
图 5‑4 位移图和图 5‑5 各单元图分别给出了结构的位移图和
应力云图。
图 5‑5 各单元
图
从图 5‑4 位移图和图 5‑5 各单元图可以看到,跨中的位移和正应力最大,与简支梁承受均布荷载,跨中挠度和正应力最大的力学常识相符合,可以初步认为模型是正确的。#表格# 1给出了
的截面上的正应力
和表格 2给出了
处各横截面
方向位移,其中
的单位为
,
的单位为
。
表格 1
考查点的y/m
1.5
1.0
0.5
0
-0.5
-1.0
-1.5
正应力
-270.20
-178.25
-88.57
0.00
88.57
178.25
270.21
表格 2
考查点的x/m
0
1.5
3.0
4.5
6.0
7.5
9.0
方向位移(10-6)
-0.3485
-0.3380
-0.3069
-0.2571
-0.1914
-0.1133
0
1.7. 自编程序分析结果
用自编程序进行分析时,采用了整个结构模型进行计算,其他条件不变,因编程水平有限,在后处理阶段没有给出节点处的位移与应力,而只能得到高斯积分点处的位移与应力,得到积分点处的应力后,根据数值分析相关知识采用了插值进行处理,从而得到
的截面上的正应力
和
处各横截面
方向位移,分别列出表格如下:
表格 3
考查点的y/m
1.5
1.0
0.5
0
-0.5
-1.0
-1.5
正应力
-297.2
-196.06
-93.25
0.00
93
196.08
297.23
表格 4
考查点的x/m
0
1.5
3.0
4.5
6.0
7.5
9.0
方向位移(10-6)
-0.3481
-0.3376
-0.3065
-0.2568
-0.1910
-0.1125
0
1.8. 结果分析比较
为了更直观的比较自编程序和ANSYS的计算结果,将表格 1和表格 3的数据绘入图 5‑6,将表格 2和表格 4的数据绘入图 5‑7。
图 5‑6 应力比较
图 5‑7 位移比较
自编程序所得结果与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] 赵翔龙.FORTRAN 90学习教程[M].北京:北京大学出版社,2002
附录 程序源代码
#include
#include
/*The gemotry of the model*/
float mesh[2];
float wide,hight;
float wsize,hsize,young,poiss,thick;
int i,j,k,knelem,npoin,elem[500],ielem;
float coord[2][1],props[3][1],lnods[8][500],ifpre[1];
float posgp[3],weigp[3],estif[136],elcod[2][8];
int npoin,nelem,kevab,igaus,jgaus,kgasp,ngash,djacb;
float gpcod[2][9],bmatx[3][16],dmatx[3][3];
float shape[8],deriv[2][8];
float xjacm[2][2],xjaci[2][2],elcod[2][8];
float cartd[2][8];
float bmatx[3][16],dmatx[3][3];
float nload[1];
float press[3][2],pgash[2],dgash[2];
struct node
{ int nodenu;/*The number of the node*/
float cor[2];/*The coordinate of the node*/
float displacement[2];/*The displacement of the node*/
float load[2]; /*The load of the node*/
float boundary[2];
} node[500];
struct
{ int elementnu;/*The number of element*/
int localnu[8];/*Local number*/
int localcorx[8];
int localcory[8]; /*Local coordinate of node*/
float qx,qy,t;/*The stress and strain*/
}element[500];
void meshing(float wide,float hight,float mesh[2])
{ printf("Please input the meshing size\n");
scanf("%f%f",&wsize,&hsize);
getchar();
mesh[0]=wide/wsize;
mesh[1]=hight/hsize;
}
/*The global coordinate of node*/
void coordinate()
{
int i,j=1;
float x,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()
{ int top[500],bottom[500];/*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++)
lnods[i][1]=i;
for(i=1;i<=3;i++)
scanf("%f",&props[i][1]);
printf("The EX is %e\nThe uo is %8f\nThe bt is %8f",props[1][1],props[2][1],props[3][1]);
getch();
printf("Please input the gemotry of the model\n");
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 size is %f,the hight size is %f\n",mesh[0],mesh[1]);
input();
getchar();
nelem=mesh[0]*mesh[1];
for(i=1;i<=nelem;i++) /*The element number*/
element[i].elementnu=i;
for(i=1;i8,then it is 1*/
if(k>8)
{ n=1;
m=2;
}
/*sum to get the equal node load*/
eload[n]=eload[n]+shape[kount]*pxcom*dvolu;
eload[m]=eload[m]+shape[kount]*pycom*dvolu;
}
}
}
print("The element load matrx is/n");
for(i=1;i<=16;i++)
printf("%f\n",eload[i]);
}
/*To form global stiffness matrix and load matrix*/
void asem()
{float v[1],a[1];
float lnods[8][1],ifpre[1],nload[1],maxa[1];
float ncodf[2],estif[136],eload[16];
float dlarge,x;
int ipoin,kmin,lnode,khigh,kdofn,nn,nnm,nwk,iwk,in;
int inodi,inode,icodf,idofn,ievab,icoln,jnode,lnodj,jdofn,jevab,jrow,ldis,lnodi;
maxa[1]=1;
for(ipoin=1;ipoin<=npoin;ipoin++)
{ kmin=npoin+1;
for(ielem=1;ielem<=nelem;ielem++)
for(i=1;i<=8;i++)
{ lnode=lnods[i][ielem];
if(lnode==ipoin)
break;
}
for(i=1;i<=8;i++)
{ lnode=lnods[i][ielem];
if(lnode>=kmin)
break;
}
kmin=lnode;
khigh=(ipoin-1)*2+j+1;
for(j=1;j<=2;j++)
{ kdofn=(ipoin-1)*2+j+1;
maxa[kdofn]=maxa[kdofn-1]+khigh+j;
}
}
nn=npoin*2;
nnm=nn+1;
nwk=maxa[nnm]-1;
printf("The sum is %d\n",nwk);
for(iwk=1;iwk<=nwk;iwk++)
a[iwk]=0.0;
for(in=1;in<=nn;nn++)
v[in]=0.0;
dlarge=1.0e+15;
for(ielem=1;ielem<=nelem;ielem++)
{ printf("The number of element is %d\n",ielem);
for(inode=1;inode<=8;inode++)
{ inodi=lnods[inode][ielem];
icodf=ifpre[lnodi];
ncodf[1]=icodf/10;
ncodf[2]=icodf-ncodf[1]*10;
for(idofn=1;idofn<=2;idofn++)
{ ievab=(inode-1)*2+idofn;
icoln=(lnodi-1)*2+idofn;
for(jnode=1;jnode<=8;jnode++)
{ lnodj=lnods[jnode][ielem];
for(jdofn=1;jdofn<=2;jdofn++)