八节点等参单元平面有限元
八节点等参单元平面有限元
分析程序
土木工程学院
目录
1.概述 (1)
2.编程思想 (2)
1 (2)
2 (2)
2.1.八节点矩形单元介绍 (2)
1 (5)
2 (5)
2.1 (5)
2.2.有限元分析的模块组织 (5)
3.程序变量及函数说明 (6)
3. (6)
3 (6)
3.1.主要变量说明： (6)
3.1. (7)
3.2.主要函数说明 (7)
4.程序流程图 (8)
5.程序应用与ANSYS分析的比较 (9)
4 (9)
5 (9)
5.1.问题说明 (9)
5.2.ANSYS分析结果 (10)
5.3.自编程序分析结果 (12)
5.4.结果分析比较 (12)
参考文献 (15)
附录程序源代码 (16)
《计算力学》课程大作业
1.概述
通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。因此，必须寻找新的编程语言。
随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。
C语言最初是为操作系统、编译器以及文字处理等编程而发明的。随着不断完善，它已应用到其它领域，包括工程应用软件的编程。近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。用C 编写的程序，可移植性极好，几乎不用作多少修改，就可在任何一台装有ANSI、C编译器的计算机上运行。C既是高级语言，也是低级语言，也就是说，可用它作数值计算，也可用它对计算机存储进行操作。
2. 编程思想
本程序采用C 语言编程，编制平面四边形八节点等参元程序，用以求解平面结构问题。程序采用二维等带宽存储整体刚度矩阵，乘大数法引入约束，等带宽高斯消去法求解位移。
在有限元程序中，变量数据需赋值的可分为节点信息，单元信息，载荷信息等。对于一个节点来说，需以下信息：节点编号（整型），节点坐标（实型），节点已知位移（实型），节点载荷（实型），边界条件（实型）和节点温度（实型）等。同样，对于一个单元来说，需以下信息：单元的节点联接信息（整型），单元类型信息（桁架、梁、板、壳等）（整型） ，单元特性信息（厚度、内力矩等）（实型），材料信息（弹性模量，泊松比等）（实型）等。
在FORTRAN 程序中，以上这些变量混合在一起，很难辨认，使程序的可读性不好，如需要进行单元网络的自适应划分，节点及单元的修改将非常困难。在进行C 语言编译过程中，采用结构struct 使每个节点信息存储在一个结构体数组中，提高程序的可读性，使数据结构更趋于合理。
2.1. 八节点矩形单元介绍
八节点矩形单元编号如错误!未找到引用源。所示
八节点矩形单元的位移函数为：
222212345678u x y x xy y x y xy αααααααα=+++++++ (2.1)
2222910111213141516v x y x xy y x y xy αααααααα=+++++++ (2.2)
其形函数为
1122334455667788u N u N u N u N u N u N u N u N u =+++++++ (2.3) 1122334455667788v N v N v N v N v N v N v N v N v =+++++++ (2.4)
式(2.3)和式(2.4)中(,)i i N N εη=，并且采用等参变换，则单元的坐标变换式可取为
1122334455667788
1122334455667788X N x N x N x N x N x N x N x N x Y N y N y N y N y N y N y N y N y =+++++++??
=+++++++?
(2.5)
图 2-1
单元应变矩阵为
{}x y xy u x u y u v y x εεεγ??????????
?????==?????????????
??+??????
(2.6)
式(2.6)一般简写为
{}[]{}B εδ= (2.7)
其中[]B 的子块矩阵为
[]i i i i i N x N B y N N y x ?????????
?
?=?
??????
????????
(2.8) 由于i N 是ε、η的函数，在[]i B 中的x 、y 要按照复合函数来求导，即
[]i i i i i i N N N x y x x J N N N x y y y εεεηηη?????????
??????
????????????
???????==????????????????????
??????????????
(2.9) 从而有
[]1
i i i i N N x J N N y εη-??????
??????????=????????????????
(2.10) 因此，单元应力矩阵为
{}[][]{}D B σδ= (2.11) 单元刚度矩阵为
[][][][]T
e
A
K B D B hdxdy =?? (2.12)
其中积分采用三点高斯积分，
33
11
,11
11
1
(,)()
(,)nip
i
j
i j i
i
i j i f d d f W f ξηξη
ωω
ξηξη--===∑∑∑??
(2.13)
2nip n =（高斯积分点的总数），i ω和j ω或i W 是加权系数，i ξ和j η是单元内的坐标.。对
于三点高斯积分，高斯积分点的位置：
11 5.0ξω==，220.0,8.0ξω==，
33 5.0ξω==。
单元等效节点荷载
{}[]
{}T
e
S P N P ds =? (2.14)
结构刚度矩阵
e
e
K K
=∑ (2.15)
结构结点荷载列阵
e e
P P =∑ (2.16)
注意，对于式(2.15)和式(2.16)中e
∑
的理解不是简单的叠加而是集成。
总刚平衡方程
[]{}{}K P δ= (2.17)
从式(2.17)求出
{}[]
{}1
K P δ-= (2.18)
将{}δ回代入式(2.7)和式(2.11)，得到{}ε和{}σ。
2.2.有限元分析的模块组织
一个典型的有限元分析过程主要包括以下几个步骤：
1)读输入数据，定义节点及单元数组。
2)由边界条件计算方程个数，赋值荷载列阵。
3)读入在带状存储的总刚度矩阵中单元和载荷信息。
4)定义总刚度阵数组。
5)组装总刚度阵。
6)解方程组。
其流程图可见错误!未找到引用源。。
图2-2
3.程序变量及函数说明3.1.主要变量说明：
3.2.主要函数说明
4.程序流程图
任何一个有限元程序处理过程，都可以划分为三个阶段：
1)前处理：读入数据和生成数据。
2)处理器：有限元的矩阵计算、组集和求解。
3)后处理：输出节点位移、单元应变等。
为了更清楚地描述程序，我们给出了程序的流程图。
5.程序应用与ANSYS分析的比较
为了验证程序的正确性，我们取了一个简单框架结构，分别用自编程序和ANSYS进行分析和比较。
5.1.问题说明
错误!未找到引用源。所示一简支梁，高3m，长18m，承受均布荷载2
10/
N m，10
210
E Pa
=?，0.167
μ=，取1
t=m，作为平面应力问题。
由于结构和荷载对称，只对右边一半进行有限单元计算，如错误!未找到引用源。所示，而在y轴各节点布置水平连杆支座。
图5-1
图5-2
5.2. ANSYS 分析结果
ANSYS 分析中，采用plane82单元，在板单元上边施加均布荷载210/N m ，在y 轴上的各结点约束UX 方向的自由度，在板单元右下角的结点约束UX 自由度，结点布置、网格划分方案如错误!未找到引用源。所示。
图 5-3
错误!未找到引用源。和错误!未找到引用源。分别给出了结构的位移图和X 应力云图。
图 5-5 各单元X 图
图 5-4 位移图
从错误!未找到引用源。和错误!未找到引用源。可以看到，跨中的位移和正应力最大，与简支梁承受均布荷载，跨中挠度和正应力最大的力学常识相符合，可以初步认为模型是正确的。错误!未找到引用源。给出了0.75x m =的截面上的正应力X σ和错误!未找到引用源。给出了 1.5y m =处各横截面Y 方向位移，其中
X σ的单位为a P ，y δ的单位为m 。
表格 2
5.3. 自编程序分析结果
用自编程序进行分析时，采用了整个结构模型进行计算，其他条件不变，因编程水平有限，在后处理阶段没有给出节点处的位移与应力，而只能得到高斯积分点处的位移与应力，得到积分点处的应力后，根据数值分析相关知识采用了插值进行处理，从而得到0.75x m =的截面上的正应力X σ和 1.5y m =-处各横截面Y 方向位移，分别列出表格如下：
表格 3
表格 4
5.4. 结果分析比较
为了更直观的比较自编程序和ANSYS 的计算结果，将错误!未找到引用源。和错误!未找到引用源。的数据绘入错误!未找到引用源。，将错误!未找到引用源。和错误!未找到引用源。的数据绘入错误!未找到引用源。。
图5-6 应力比较
图5-7 位移比较
自编程序所得结果与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]赵翔龙.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 %,the hight is %",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;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;