平面四节点等参单元分析程序

变分原理与有限元大作业

平面四节点等参单元分析程序；

姓名：潘清

【

学号：SQ

完成时间：2011-4-26

|

一、概述

!

通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、ABAQUS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。因此，必须寻找新的编程语言。

随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。

C语言最初是为操作系统、编译器以及文字处理等编程而发明的。随着不断完善，它已应用到其它领域，包括工程应用软件的编程。近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。用C 编写的程序，可移植性极好，几乎不用作多少修改，就可在任何一台装有ANSI、C编译器的计算机上运行。C既是高级语言，也是低级语言，也就是说，可用它作数值计算，也可用它对计算机存储进行操作。

二、编程思想

;

本程序采用C语言编程，编制平面四边形四节点等参元程序，用以求解平面结构问题。程序采用二维等带宽存储整体刚度矩阵，乘大数法引入约束，等带宽高斯消去法求解位移，然后求中间高斯点的应力，最后用绕节点平均法讲单元应力等效到节点上，再将结果写到tecplot文件中。

在有限元程序中，变量数据需赋值的可分为节点信息，单元信息，载荷信息等。对于一个节点来说，需以下信息：节点编号（整型），节点坐标（实型），节点已知位移（实型），节点载荷（实型），边界条件（实型）等。同样，对于一个单元来说，需以下信息：单元的

节点联接信息（整型），材料信息（弹性模量，泊松比等）（实型）等。

在FORTRAN

程序中，以上这些变量混合在一起，很难辨认，使程序的可读性不好，如

需要进行单元网络的自适应划分，节点及单元的修改将非常困难。在进行C语言编译过程中，

采用结构struct 使每个节点信息存储在一个结构体数组中，提高程序的可读性，使数据结构更趋于合理。

三、平面四节点等参单元介绍

*

四节点等参单元实际单元与基本单元的映射关系如图3-1所示

(

坐标的映射关系为：

其位移模式和坐标的映射有相同的插值函数，形函数为：

<

单元应变矩阵为：

`

图3-1

{}x y xy u x u y u v y x εεεγ??????????

?????==?????????????

??+??????

上式一般简写为：

{

{}[]{}

B εδ=

其中[]B 的子块矩阵为

[]i i i i i N x N B y N N y x ??????????

?=?

??????

????????

由于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 εεεηηη?????????

??????

????????????

???????==????????????????????

??????????????

！

从而有

[]1i i i i N N x J N N y εη-??????

??????????=????????????????

因此，单元应力矩阵为

{}[][]{}D B σδ=

单元刚度矩阵为

<

[]

[][][]T

e

A

K B D B hdxdy =??

其中积分采用三点高斯积分，

33

11

,11

11

1

(,)()

(,)nip

i

j

i j i

i

i j i f d d f W f ξηξη

ωω

ξηξη--===∑∑∑??

2nip n =（高斯积分点的总数），i ω和j ω或i W 是加权系数，i ξ和j η是单元内的坐标.。对

于三点高斯积分，高斯积分点的位置：

11 5.09.0ξω==，220.0,8.0ξω==，

33 5.0ξω==。

！

结构刚度矩阵

e e

K K =∑

结构结点荷载列阵

e e

P P =∑

注意，对于上两式中

e

∑

的理解不是简单的叠加而是按照对应的自由度集成。

】

总刚平衡方程

[]{}{}K P δ=

从式上式求出：

{}[]{}1

K P δ-=

*

四、有限元分析的模块组织

一个典型的有限元分析过程主要包括以下几个步骤： 1) 读输入数据，定义节点及单元数组。 2) 由边界条件计算方程个数，赋值荷载列阵。 3) 读入在带状存储的总刚度矩阵中单元和载荷信息。 4) (

5)

定义总刚度阵数组。

6) 组装总刚度阵。 7) 解方程组。

其流程图可见下图

五、程序变量及函数说明

、

1、控制信息

np：结构节点总数

ne：结构离散单元总数

nr1，nr2：总的约束的节点数，nr1，x方向；nr2，y方向

nd：每个单元的节点数

[

nf：每个节点的自由度数

ld：集中力载荷的个数

nm：材料的种类

nu1，nu2：非零位移边界条件的节点数，nu1，x方向；nu2，y方向

u1，u2：非零位移的大小，u1，x方向；u2，y方向

%

n=nf*np ：结构的节点位移总数

ndf=nd*nf ：每个单元的节点自由度数

2、输入的原始数据

x(np):节点的x方向坐标

~

y(np):节点的y方向坐标

me（nd，ne）：单元节点的总体编号

nrr(nr1+nr2) ：约束为零的位移所对应的总体位移编号

p（n）：载荷向量

nu（nu1+nu2）：位移载荷

《

mat（6，nm）：材料参数

3、程序中的其他标识符

LD（n）：存放结构刚度阵所以主对角线元素在A（nn）中的序号

IS（ndf）：单元节点位移和节点力在总体位移阵列和载荷阵列中对应的序号~

EK（ndf，ndf）：总体坐标系下的单元刚度矩阵

A（nn）：架构刚度阵下三角变带宽一维压缩存储的数组

nn：数组A的元素个数

RSTG(3):高斯积分点的值

H(3):高斯积分点的加权系数

【

S（6，ne）：各单元的应力分量

XJAC（2,2）：雅阁比矩阵

RJAC（2,2）：雅阁比矩阵的逆

PN（2,4）：4个节点处形函数对局部坐标的导数

DNX（2,4）：4个节点处形函数对整体坐标的导数)

FUN（4）：形函数的值

，

六、计算流程图

%

七、计算结果与abaqus分析结果的比较

1、中间带圆孔平面应力板的分析。

宽40m，长50m，圆孔位于板中心，半径为4m，承受水平方向位移载荷，E=200Pa，

3.0

=

μ，取1

t=m。用abaqus建模离散，并计算。再讲abaqus离散的节点和单元信息拷贝到本程序的输入文件中，用该程序计算，结果输出成tecplot文件，用tecplot可以查看结果。与abaqus的计算结果进行比较，位移和应力云图如下（左边是程序计算结果，右边是abaqus结果，下同）：

：

输出节点位移和应力

2、纤维增强复合材料平面应力板的分析。

;

宽40m ，长50m ，增强纤维位于板中心，纤维半径为8m ，承受水平方向位移载荷，基体材料E=20000Pa ，25.0=μ，纤维材料E=8000Pa ，3.0=μ，取1t =m 。

、

3、半无限大含裂纹板的应力分析

宽40m ，长60m ，裂纹位于板左侧中间位置，裂纹长10m ，承受竖直方向位移载荷，E=200Pa ，

3.0=μ，取1t =m 。

$

结果分析比较：

通过以上三个算例发现该程序可以用于计算单材料、双材料、带孔、含裂纹等各种平面问题，加载条件可以是加集中力和加位移，因此，该程序的适用范围还是比较广的。

以上三个算例的自编程序所得结果与abaqus分析结果进行比较发现，两者的计算结果很接近，而且自编程序对于孔边应力集中和裂尖应力集中都能很好的表达，说明该程序有很好的精度。

第三个算例在裂尖处数值上有些区别，但总的分布规律还是很吻合的，主要是因为本程序是用四节点等参单元，对于应力的奇异性表达效果还不是很好。

#

八、附录

·

附录一：程序代码

#include<>

#include<>

#include<>

~

#include<>

#include<>

float **float_two_array_malloc(int m,int n) ....\n");

}

/

void output() f\t",s[j][i]);

fprintf(fp2,"\n");

}

)

fclose(fp2);

printf("写入完成......\n");

}

(

void tecplot() f\t%.4f\t%.4f\t%.4f\t",x[j]+bl*p[j*2],y[j]+bl*p[j*2+1],p[j*2],p[j*2+1]);

for(i=0;i<6;i++)

{

fprintf(fp3,"%.4f\t",ns[i][j]);

\

}

fprintf(fp3,"\n");

}

for(i=0;i

|

{

for(j=0;j

fprintf(fp3,"%d\t",me[j][i]+1);

fprintf(fp3,"\n");

}

*

fclose(fp3);

printf("写入完成......\n");

}

|

void fld() ....\n");

}

void fis(int ie) ....\n");

@

void cholesky(float*A) ....\n");

}

void stress(int ie)....\n");

}

void main()

{

int i,j,k;

float *A;

input();

fld();

A=float_one_array_malloc(nn);

for(i=0;i

fr(A);

cholesky(A);

for(j=0;j

{

stress(j);

}

printf("求单元应力完成......\n");

stress_average();

output();

tecplot();

}

附录二：输入的原始数据的格式：

9 4 3 1 4 2 3 1 3 0 /依次为np,ne,nr1,nr2,nd,nf,ld,nm,nu1,u1,nu2,u2 -20 15 //各节点的x,y坐标

-20 0

-20 -15

0 15

0 0

0 -15

20 15

20 0

20 -15

1 2 5 4 //各单元节点的整体编号

2 3 6 5

4 5 8 7

5 6 9 8

1 2 3 //x方向约束为零的位移对应的节点编号3 //y方向约束为零的位移对应的节点编号12 10 //集中力载荷

14 40

16 10

1 4 200 //材料信息——起、止单元号，E，MU 7 8 9 //位移非零的节点号