平面三角形单元常应变单元matlab程序的编制.doc
1三角形常应变单元程序的编制与使用有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法(或变分里兹法)而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。有限元分析的基本步骤可归纳为三大步结构离散、单元分析和整体分析。对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好的适应性，而矩形或四边形单元较三节点三角形有更高的计算精度。MATLAB语言是进行矩阵运算的强大工具，因此，用MATLAB语言编写有限元中平面问题的程序有优越性。本章将详细介绍如何利用MATLAB语言编制三角形常应变单元的计算程序，程序流程图见图1。有限元法中三节点三角形分析结构的步骤如下1)整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。2)单元分析，建立单元刚度矩阵。3)整体分析，建立总刚矩阵。4)建立整体结构的等效节点荷载和总荷载矩阵5)边界条件处理。6)解方程，求出节点位移。7)求出各单元的单元应力。8)计算结果整理。计算结果整理包括位移和应力两个方面；位移计算结果一般不需要特别的处理，利用计算出的节点位移分量，就可画出结构任意方向的位移云图；而应力解的误差表现在单元内部不满足平衡方程，单图1程序流程图开始输入初始数据生成单刚集成总刚施加约束信息生成荷载向量边界条件处理计算结点位移计算单元应力计算结果整理结束2元与单元边界处应力一般不连续，在边界上应力解一般与力的边界条件不相符合。11程序说明三角形常应变单元求解结构主程序功能运用有限元法中三角形常应变单元解平面问题的计算主程序。基本思想单元结点按右手法则顺序编号。荷载类型可计算结点荷载。说明主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。1程序准备ATSHORTE设定输出类型CLEARALL清除所有已定义变量CLC清屏说明ATSHORTE－设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法；CLEARALL－清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况；CLC－清屏，使屏幕在本程序运行开始时2全局变量定义GLOBALNNODENPIONNELEMNVFIXNFORCECOORDLNODSYOUNGPOISSTHICKGLOBALFORCEFIXEDGLOBALBMATXDMATXSMATXAREAGLOBALASTIFASLODASDISPGLOBALFP1说明3NNODE单元结点数，NPION总结点数，NELEM单元数，NVFIX受约束边界点数，NFORCE结点力数，COORD结构结点坐标数组，LNODS单元定义数组，YOUNG弹性模量，POISS泊松比，THICK厚度FORCE节点力数组N,3N受力节点数目,N,1作用点,N,2X方向,N,3Y方向；FIXED约束信息数组N,3N受约束节点数目,N,1约束点N,2与N,3分别为约束点X方向和Y方向的约束情况,受约束为1否则为0BMATX单元应变矩阵36，DMATX单元弹性矩阵33，SMATX单元应力矩阵36，AREA单元面积ASTIF总体刚度矩阵，ASLOD总体荷载向量，ASDISP结点位移向量FP1数据文件指针3打开文件FP1 TXT , RT 打开输入数据文件存放初始数据说明FP1 TXT , RT －打开已存在的输入数据文件TXT，且设置其为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1来执行FP2 OUTPUTTXT , WT －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。该文件设置为可写格式，可在程序执行过程中向输出文件写入数据。4读入程序控制信息NPIONFSCANFFP1, D ,1结点个数(结点编码总数)NELEMFSCANFFP1, D ,1单元个数(单元编码总数)NFORCEFSCANFFP1, D ,1％结点荷载个数NVFIXFSCANFFP1, D ,1％受约束边界点数YOUNGFSCANFFP1, E ,1％弹性模量POISSFSCANFFP1, F ,1％泊松比THICKFSCANFFP1, D ,1厚度LNODSFSCANFFP1, D ,3,NELEM 单元定义数组(单元结点号)说明4建立LNODS矩阵，该矩阵指出了每一单元的连接信息。矩阵的每一行针对每一单元，共计NELEM；每一列相应为单元结点号(编码)、按逆时针顺序输入。命令中，3,NELEM’表示读取NELEM行3列数据赋值给LNODS矩阵。显然，LNODSI,13依次表示I单元的I，J，K结点号。COORDFSCANFFP1, F ,2,NPION 结点坐标数组说明建立COORD矩阵，该矩阵用来存储各结点X，Y方向的坐标值。从FP1文件中读取全部结点个数NPOIN的坐标值，从1开始按顺序读取。COORDI，12表示第I个结点的X，Y坐标。FORCEFSCANFFP1, F ,3,NFORCE 结点力数组说明N,3N受力结点数目,N,1作用点,N,2X方向,N,3Y方向FIXEDFSCANFFP1, D ,3,INF 约束信息数组说明N,3N受约束节点数目,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调用子程生成单刚，组成总刚并加入约束信息ASSEMBLE6调用子程生成荷载向量5LOAD7计算结点位移向量ASDISPASTIFASLOD 8调用子程计算单元应力WRITESTRESS9关闭输出数据文件FCLOSEFP2读取ASSEMBLE子程FUNCTIONASSEMBLE所引用的全局变量GLOBALNPIONNELEMNVFIXLNODSASTIFTHICKGLOBALBMATXSMATXAREAFIXED计算单刚并生成总刚ASTIF12NPION,12NPION0张成特定大小总刚矩阵并置0说明建立单元刚度矩阵ASTIF，该矩阵的行列数均为2NPION，NPION表示结点数，每个结点有两个方向的力和位移。FORI1NELEMSMATXI调用应力子程序ESTIFBMATX SMATXTHICKAREA求解单元刚度矩阵ALNODSI,临时向量,用来记录当前单元的节点编号FORJ13FORK13ASTIFAJ21AJ2,AK21AK2A