平面三角形单元常应变单元matlab程序的编制.docx
《平面三角形单元常应变单元matlab程序的编制.docx》由会员分享，可在线阅读，更多相关《平面三角形单元常应变单元matlab程序的编制.docx（15页珍藏版）》请在冰豆网上搜索。
平面三角形单元常应变单元matlab程序的编制
平面三角形单元常应变单元
matlab程序的编制
三角形常应变单元程序的编制与使用
有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于
变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算
机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。
有限元分析的基本步骤可归纳为三大步：
结构离散、单元分析和整体分析对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三
个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边
界有更好的适应性，而矩形或四边形单元较三节点三角形有更高的计算精度。
Matlab语言是进行矩阵运算的强大工具，因此，用Matlab语言编写有限兀中平面冋题的程序有优越性。
本章将详细介绍如何利用Matlab语言编制三角形常应变单元的计算程序，程序流程图见图1。
有限元法中三节点三角形分析结构的步骤如
下：
1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2）单元分析，建立单元刚度矩阵。
3）整体分析，建立总刚矩阵。
4）建立整体结构的等效节点荷载和总荷载矩
阵
5）边界条件处理。
6）解方程，求出节点位移。
7）求出各单元的单元应力。
8）计算结果整理。
计算结果整理包括位移和应力两个方面；位移计算结果一般不需要特别的处理，利用计算出的节点位移分量，就可画出结构任意方向的位移云图；而应
力解的
误差表现在单元内部不满足平衡方程,
开始
J
输入初始数据
审
生成单刚集成总刚
*
施加约束信息
生成荷载向量
边界条件处理
计算结点位移
计算单元应力
计算结果整理
结束
单元
与单元边界处应力一般不连续，在边界上应力解一般与力的边界条件不相符合。
1.1程序说明
%%（******************************************************************
%三角形常应变单元求解结构主程序
%%（******************************************************************
功能：
运用有限元法中三角形常应变单元解平面问题的计算主程序。
基本思想：
单元结点按右手法则顺序编号。
荷载类型：
可计算结点荷载。
说明：
主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。
%
1程序准备
formatshorte%设定输出类型
clearall%清除所有已定义变量
clc%清屏
说明：
formatshorte-设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法；
clearall—清除所有已定义变量，目的是在本程序的运行过程中，
不会发生变量名相同等可能使计算出错的情况；
clc—清屏，使屏幕在本程序运行开始时
%
2全局变量定义
globalNNODENPIONNELEMNVFIXNFORCECOORD
LNODSYOUNGPOISSTHICK
global
FORCE
FIXED
global
BMATX
DMATX
SMATX
AREA
global
ASTIF
ASLOD
ASDISP
globalFP1
说明：
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打开文件
FP仁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（）
Q%*******************************************************************
9关闭输出数据文件
fclose（FP2）;
Q%*******************************************************************
读取ASSEMBLE子
壬阜0/*****************************************************************
functionASSEMBLE（）
%所引用的全局变量：
globalNPIONNELEMNVFIXLNODSASTIFTHICKglobalBMATXSMATXAREAFIXED
%
%计算单刚并生成总刚
ASTIF（1:
2*NPION,1:
2*NPION）=O;%张成特定大小总刚矩阵并
置0
说明：
Anayjstiffness
建立单元刚度矩阵ASTIF，该矩阵的行列数均为2*NPION，NPION表示结点数，每个结点有两个方向的力和位移%
fori=1:
NELEM
FORMSMATX（i）%调用应力子程序
ESTIF=BMATX'*SMATX*THICK*AREA;%求解单元刚度矩阵
a=LNODS（i,:
）;%临时向量，用来记录当前单元的
节点编号
forj=1:
3
fork=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即将单元中的位移编码对应到整体的位移编码。
%
NVFIX受约束边界点数
%将约束信息加入总刚（置0置1法）NUM=1;%计数器（当前已分析的节点数）i=0;%计数器（当前已处理的约束数）
tmp（NVFIX）=0;%临时存被约束的序号
whileivNVFIX
%若发现约束
%计数器自增
%求约束序号
forj=-1:
0
ifFIXED（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（
升级会员 