平面三角形单元常应变单元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 THICKglobal FORCE FIXEDglobal BMATX DMATX SMATX AREAglobal ASTIF ASLOD ASDISPglobal 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否则为0BMATX—单元应变矩阵(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 调用子程生成单刚，组成总刚并加入约束信息ASSEMBLE()%----------------------------------------------------------------------------------------------------- 6 调用子程生成荷载向量FORMLOAD()%----------------------------------------------------------------------------------------------------- 7 计算结点位移向量ASDISP=ASTIF\ASLOD'%----------------------------------------------------------------------------------------------------- 8调用子程计算单元应力WRITESTRESS()%******************************************************************* 9 关闭输出数据文件fclose(FP2);%*******************************************************************读取ASSEMBLE子程%*******************************************************************function ASSEMBLE()% 所引用的全局变量：global NPION NELEM NVFIX LNODS ASTIF THICKglobal BMATX SMATX AREA FIXED%------------------------------------------------------------------------------------------------ ---- % 计算单刚并生成总刚ASTIF(1:2*NPION,1:2*NPION)=0; %张成特定大小总刚矩阵并置0 ●说明：建立单元刚度矩阵ASTIF，该矩阵的行列数均为2*NPION ，NPION表示结点数，每个结点有两个方向的力和位移。

for i=1:NELEMFORMSMATX(i) %调用应力子程序ESTIF=BMATX'*SMATX*THICK*AREA; %求解单元刚度矩阵a=LNODS(i,:); %临时向量,用来记录当前单元的节点编号for j=1:3for k=1:3ASTIF((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);%跟据节点编号对应关系将单元刚度分块叠加到总刚矩阵中endendend●说明：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即将单元中的位移编码对应到整体的位移编码。