三角形常应变单元matlab程序的编制与使用

三角形常应变单元程序的编制与使用

有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。

有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好的适应性，而矩形或四边形单元较三节点三角

形有更高的计算精度。

Matlab语言是进行矩阵运算的强大工具，因

此，用Matlab语言编写有限元中平面问题的程序

有优越性。本章将详细介绍如何利用Matlab语言

编制三角形常应变单元的计算程序，程序流程图见

图1。

有限元法中三节点三角形分析结构的步骤如

下：

1）整理原始数据，如材料性质、荷载条件、约

束条件等，离散结构并进行单元编码、结点

编码、结点位移编码、选取坐标系。

2）单元分析，建立单元刚度矩阵。

3）整体分析，建立总刚矩阵。

4）建立整体结构的等效节点荷载和总荷载矩

阵

5）边界条件处理。

6）解方程，求出节点位移。

7）求出各单元的单元应力。

8）计算结果整理。计算结果整理包括位移和应

力两个方面；位移计算结果一般不需要特别

的处理，利用计算出的节点位移分量，就可

画出结构任意方向的位移云图；而应力解的

误差表现在单元内部不满足平衡方程，单元与单元边界处应力一般不连续，在边界上应力解一般与力的边界条件不相符合。图1 程序流程图

1.1 程序说明

%******************************************************************* % 三角形常应变单元求解结构主程序

%******************************************************************* ●功能：运用有限元法中三角形常应变单元解平面问题的计算主程序。

●基本思想：单元结点按右手法则顺序编号。

●荷载类型：可计算结点荷载。

●说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算

法的全过程，即实现程序流程图的程序表达。

%----------------------------------------------------------------------------------------------------- 1 程序准备

format short e %设定输出类型

clear all %清除所有已定义变量

clc %清屏

●说明：

format short e －设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法；

clear all －清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况；

clc －清屏，使屏幕在本程序运行开始时

%----------------------------------------------------------------------------------------------------- 2 全局变量定义

global NNODE NPION NELEM NVFIX NFORCE COORD LNODS YOUNG POISS THICK

global FORCE FIXED

global BMATX DMATX SMATX AREA

global ASTIF ASLOD ASDISP

global FP1

●说明：

NNODE—单元结点数，NPION—总结点数，NELEM—单元数，NVFIX—受约束边界点数，NFORCE—结点力数，COORD—结构结点坐标数组，LNODS —单元定义数组，YOUNG—弹性模量，POISS—泊松比，THICK—厚度

FORCE —节点力数组(n,3) n:受力节点数目,(n,1):作用点,(n,2):x方向,(n,3):y 方向；FIXED—约束信息数组(n,3) n:受约束节点数目, (n,1):约束点(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0

BMATX—单元应变矩阵(3*6)，DMATX—单元弹性矩阵(3*3)，SMATX—单元应力矩阵(3*6)，AREA—单元面积

ASTIF—总体刚度矩阵，ASLOD—总体荷载向量，ASDISP—结点位移向量FP1—数据文件指针

3 打开文件

FP1=fopen('input.txt','rt'); %打开输入数据文件存放初始数据

●说明：

FP1=fopen('input.txt','rt'); －打开已存在的输入数据文件input.txt，且设置其为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1来执行

FP2=fopen('output.txt','wt'); －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。该文件设置为可写格式，可在程序执行过程中向输出文件写入数据。

%----------------------------------------------------------------------------------------------------- 4 读入程序控制信息

NPION=fscanf(FP1,'%d',1) %结点个数（结点编码总数）

NELEM=fscanf(FP1,'%d',1) %单元个数（单元编码总数）

NFORCE=fscanf(FP1,'%d',1) ％结点荷载个数

NVFIX=fscanf(FP1,'%d',1) ％受约束边界点数

YOUNG=fscanf(FP1,'%e',1) ％弹性模量

POISS=fscanf(FP1,'%f',1) ％泊松比

THICK=fscanf(FP1,'%d',1) %厚度

LNODS=fscanf(FP1,'%d',[3,NELEM])' %单元定义数组（单元结点号）

●说明：

建立LNODS矩阵，该矩阵指出了每一单元的连接信息。

矩阵的每一行针对每一单元，共计NELEM；每一列相应为单元结点号（编码）、按逆时针顺序输入。

命令中，[3,NELEM]’表示读取NELEM行3列数据赋值给LNODS矩阵。

显然，LNODS(i,1:3)依次表示i单元的i，j，k结点号。

COORD=fscanf(FP1,'%f',[2,NPION])' %结点坐标数组

●说明：

建立COORD矩阵，该矩阵用来存储各结点x，y方向的坐标值。

从FP1文件中读取全部结点个数NPOIN的坐标值，从1开始按顺序读取。

COORD(i，1：2)表示第i个结点的x，y坐标。

FORCE=fscanf(FP1,'%f',[3,NFORCE])' %结点力数组

●说明：

(n,3) n:受力结点数目,(n,1):作用点,(n,2):x方向,(n,3):y方向

FIXED=fscanf(FP1,'%d',[3,inf])' %约束信息数组

●说明：

(n,3) n:受约束节点数目, (n,1):约束点(n,2)与(n,3)分别为约束点x方向和y 方向的约束情况,受约束为1否则为0

●总体说明：

从输入文件FP1中读入结点个数，单元个数，结点荷载个数，受约束边界点数，弹性模量，泊松比，厚度，单元定义数组，结点坐标数组，结点力数组，约束信息数组；

程序中弹性模量仅输入了一个值，表明本程序仅能求解一种材料构成的结构，如：钢筋混凝土结构、钢结构，不能求解钢筋混凝土－钢组合结构。

采用了命令fscanf，其中：’%d’表示读入整数格式，’%f'’表示读入浮点数；1表示读取1个数，[A,B]形式表示读A行B列数组，[A,B]’表示将[A,B]转置，inf 表示正无穷。

%----------------------------------------------------------------------------------------------------- 5 调用子程生成单刚，组成总刚并加入约束信息

function ASSEMBLE()

%----------------------------------------------------------------------------------------------------- 6 调用子程生成荷载向量

function FORMLOAD()

%----------------------------------------------------------------------------------------------------- 7 计算结点位移向量

ASDISP=ASTIF\ASLOD'

%----------------------------------------------------------------------------------------------------- 8调用子程计算单元应力

function WRITESTRESS()

%******************************************************************* 9 关闭输出数据文件

fclose(FP2);

%*******************************************************************

读取ASSEMBLE子

程%*****************************************************************

**

function ASSEMBLE()

% 所引用的全局变量：global NPION NELEM NVFIX LNODS ASTIF THICK

global BMATX SMATX AREA FIXED

%------------------------------------------------------------------------------------------------ ---- % 计算单刚并生成总刚

ASTIF(1:2*NPION,1:2*NPION)=0; %张成特定大小总刚矩阵并置0 ●说明：

建立单元刚度矩阵ASTIF，该矩阵的行列数均为2*NPION ，NPION表示结点数，每个结点有两个方向的力和位移。

%----------------------------------------------------------------------------------------------------- for i=1:NELEM

FORMSMATX(i) %% %调用应力子程序

ESTIF=BMATX'*SMATX*THICK*AREA; %求解单元刚度矩阵

a=LNODS(i,:); %临时向量,用来记录当前单元的节点编号

for j=1:3

for k=1:3

ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=ASTIF((a(j)*2-1):a(j)*2,(

a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2);

%跟据节点编号对应关系将单元刚度分块叠加到总刚矩阵中

end

end

end

●说明：

FORMSMATX(i)调用应力子程序，提取i单元的应力矩阵SMATX；

a=LNODS(i,:)记录i单元的三个结点编号；

for…end循环语句表示行从1到3循环，列从1到3循环，将单刚中的元素叠加至总刚中：

ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)表示总刚中第a(j)*2-1到:a(j)*2行，

第a(k)*2-1到a(k)*2列的元素由单刚中第j*2-1到j*2行，第k*2-1到k*2列的元素叠加而得，a(j)*2即将单元中的位移编码对应到整体的位移编码。

%---------------------------------------------------------------------------------------------------- %将约束信息加入总刚（置0置1法）

NUM=1; %计数器(当前已分析的节点数)

i=0; %计数器(当前已处理的约束数)

tmp(NVFIX)=0; %临时存被约束的序号

while i

for j=-1:0

if FIXED(NUM,j+3)==1 %若发现约束

i=i+1; %计数器自增

tmp(i)=FIXED(NUM)*2+j; %求约束序号

end

end

NUM=NUM+1;

end

说明：

tmp(NVFIX)=0，形成一个元素值均为0的一行NVFIX列的行向量，

执行while语句，首先判断i是否小于控制数据NVFIX，若小于则往下进行，若不小于则退出。

执行for语句，FIXED(NUM,j+3)表示约束信息数组中第NUM行第j+3列的元素，j从-1到0，即j+3表示2到3列，即约束信息数组中描述结点x和y方向受约束的情况，判断FIXED(NUM,j+3)若等于1，则约束数自增，若不等于1，跳出。

FIXED(NUM)表示FIXED(NUM，1)，tmp(i)=FIXED(NUM)*2+j计算整体约束序号，将序号放入tmp行向量中的i列。

%----------------------------------------------------------------------------------------------------- for i=1:NVFIX

ASTIF(1:2*NPION,tmp(i))=0; %将一约束序号处的总刚列向量清0

ASTIF(tmp(i),1:2*NPION)=0; %.将一约束序号处的总刚行向量清0

ASTIF(tmp(i),tmp(i))=1; %将行列交叉处的元素置为1 end

● 说明：

后处理法中置0置1法，设j j C =δ（包括0=j C ），则将总刚中的主元素 K jj 换为1，j 行和j 列的其他元素均改为零。

%******************************************************************* % 读取FORMSMATX 子程

%******************************************************************* function FORMSMATX(ELEMENT) %计算应力矩阵 %引用所需的全局变量 global NPION NELEM COORD LNODS YOUNG POISS global BMATX DMATX SMATX AREA %----------------------------------------------------------------------------------------------------- %生成弹性矩阵D a=YOUNG/(1-POISS^2);

DMATX(1,1)=1*a; DMATX(1,2)=POISS*a; DMATX(2,1)=POISS*a; DMATX(2,2)=1*a;

DMATX(3,3)=(1-POISS)*a/2; ● 说明：

平面应力问题的弹性矩阵?

????

??

????

???--=21,0,00,1,0,,1122μμμμE D ,其中，E 为弹性模量，μ为泊松比。

%-----------------------------------------------------------------------------------------------------

i=ELEMENT; %i 为当前所计算的单元号 %计算当前单元的面积

AREA=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);... 1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);... 1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2);])/2; ● 说明：

det 表示求矩阵行列式的值，m

m j j i

i y x y x y x A ,,1,,1,12

1

,=，其中),,(m j i x i 分别表示一个三角形单元的三个节点坐标。

MATLAB 中若一行中无法写下一个完整的命令，则可以在行尾加入3个连续的英文句号，表示命令余下的部分在下一行出现。

%----------------------------------------------------------------------------------------------------- %生成应变矩阵B for j=0:2 b(j+1)=

、

COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2); c(j+1)=

-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1); end

BMATX=[b(1) 0 b(2) 0 b(3) 0;... 0 c(1) 0 c(2) 0 c(3);... c(1) b(1) c(2) b(2) c(3) b(3)]/(2*AREA); 说明：

应变矩阵),,(,,00,21m j i l b c c b A B l l l l l =????

??????= rem 表示求余函数，rem （x ，y ）命令返回的是x-n.*y ，当y ≠0时，n=fix(x./y)，其中fix 为最近的整数取整。

%----------------------------------------------------------------------------------------------------- SMATX=DMATX*BMATX; %求应力矩阵S=D*B

%******************************************************************* % 读取FORMLOAD 子程

%******************************************************************* function FORMLOAD() %本子程生成荷载向量 %---------------------------------------------------------------------------------------------------- %所需引用的全局变量

global ASLOD NPION NFORCE FORCE

%----------------------------------------------------------------------------------------------------

ASLOD(1:2*NPION)=0; %张成特定大小的0向量 ● 说明：

ASLOD 为总体荷载向量，每个结点有x ，y 两个方向的结点力。

%---------------------------------------------------------------------------------------------------- for i=1:NFORCE

ASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);

end ● 说明：

FORCE(i,1)为作用点,FORCE(i,2:3)为x,y 方向的结点力。

%---------------------------------------------------------------------------------------------------- %将有约束处的荷载置为0 NUM=1;

i=0;

tmp(NVFIX)=0; while i

if FIXED(NUM,j+3)==1 i=i+1; tmp(i)=FIXED(NUM)*2+j; end end

NUM=NUM+1; end

for i=1:NVFIX

ASLOD(tmp(i))=0; end

● 说明：

置0置1法，同上。

%******************************************************************* ASDISP=ASTIF\ASLOD' %计算结点位移向量 %******************************************************************* % 读取WRITESTRESS 子程

%******************************************************************* function WRITESTRESS()

%求解内力，即两个方向的正应力与一个剪应力（xy y x τεε,,)

%---------------------------------------------------------------------------------------------------- %所引用的全局变量

global NELEM NPION SMATX ASDISP LNODS

%----------------------------------------------------------------------------------------------------

ELEDISP(1:6)=0; %当前单元的结点位移向量 ● 说明：

ELEDIS 为单元的结点位移?????

?????????????????=m m j j i i e

v u v u v u a

%---------------------------------------------------------------------------------------------------- for i=1:NELEM for j=1:3

ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2); end

FORMSMATX(i) %% %调用子程求应力矩阵 i

STRESS=SMA TX*ELEDISP' %求内力 end ● 说明：

通过循环，依次从结构位移列阵中对号，赋值给第i 个单元的单元位移向量ELEDISP 。

%----------------------------------------------------------------------------------------------------

1.2 程序应用举例

【例1】利用三角形三节点常应变单元程序计算图2所示的结构。

设弹性模量1 E ，泊松比为0，厚度为1。 %-----------------------------------------------------------------

输入数据文件input.txt 为： 6 4 1 6 1.0e0 0.0 1 3 1 2 5 2 4

3 2 5 6 3 5 图2 0.0 2.0 0.0 1.0 1.0 1.0 0.0 0.0 1.0 0.0 2.0 0.0 1 0 -1 1 1 0 2 1 0

4 1 1

5 0 1

6 0 1

%----------------------------------------------------------------------------------------------- 说明：

第一行：读入程序控制信息

NPION=fscanf(FP1,'%d',1) %结点个数（结点编码总数） NELEM=fscanf(FP1,'%d',1) %单元个数（单元编码总数） NFORCE=fscanf(FP1,'%d',1) ％结点荷载个数 NVFIX=fscanf(FP1,'%d',1) ％受约束边界点数 YOUNG=fscanf(FP1,'%e',1) ％弹性模量

136

5

4

2

POISS=fscanf(FP1,'%f',1) ％泊松比

THICK=fscanf(FP1,'%d',1) %厚度

第二、三、四、五行：读入单元连接信息：

LNODS=fscanf(FP1,'%d',[3,NELEM])'; %单元定义数组，单元结点号，逆时针输入

第六、七、八、九、十、十一行：读入结点坐标

COORD=fscanf(FP1,'%f',[2，NPOIN])'; %结点坐标值，第1列为x坐标值，第2列为y坐标值

第十一行：读入结点荷载信息

FORCE=fscanf(FP1,'%f',[3,NFORCE])'; %结点号，X向结点荷载数值，Y 向结点荷载数值（与坐标轴方向一致为正）

第十二、十六行：读入零位移信息

FIXED=fscanf(FP1,'%d',[3,inf])'; %结点号，X向约束，Y向约束%-----------------------------------------------------------------------------------------------------

1.3 程序的改进要点

上述三角形三节点程序反映了有限元的基本思路，可以计算简单的平面应力或平面应变问题。在熟练掌握了程序的编制与使用后，可在以下几方面对程序进行改进，以加深对矩阵位移法及MATLAB语言编程的理解：

1、本程序的弹性模量仅能输入一个数值，意味着程序仅能计算由同种材料构成的结构。考虑如何改进使程序可计算由不同材料构成的组合结构。

2、本程序仅能计算结点集中荷载类型，考虑如何编制体积力、表面力、跨中集中力等类型的程序。

3、考虑如何编制有支座位移的程序。

4、本程序最后的结果没有生成输出文件，可以编制输出文件。

5、本程序计算的应力没有进行结果处理，可以编制最后结果处理的程序。

6、可以在此程序基础上编制三角形六节点、四边形四节点等程序。

综上所述，本章的三角形三节点常应变程序体现了如何将有限元法的计算方法和过程用MATLAB程序语言表达出来，重点放在程序架构和流程的建立以及算法实现方面，主要依赖手工操作－手工输入初始数据（前处理）、手工绘制计

算结果（后处理），目的是使学生能够清晰、明确地掌握有限元法的基本理论和概念，熟练掌握应用MATLAB语言编制、修改和调试简单程序的能力。