平面三角形单元常应变单元matlab程序的编制

1、图1程序流程图三角形常应变单元程序的编制与使用有限元法是求解微分方程边值问题的一种通用数值方法， 该方法是一种基于 变分法或变分里兹法而开展起来的求解微分方程的数值计算方法， 以电脑为 手段，采用分片近似，进而逼近整体的研究思想求解物理问题。有限元分析的根本步骤可归纳为三大步：结构离散、单元分析和整体分析。 对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个 顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有 更好的适应性，而矩形或四边形单元较三节点三角 形有更高的计算精度。Matlab语言是进行矩阵运算的强大工具，因 此，用Matlab语言编写有限兀中

2、平面冋题的程序 有优越性。本章将详细介绍如何利用 Matlab语言 编制三角形常应变单元的计算程序，程序流程图见 图1。有限元法中三节点三角形分析结构的步骤如下：1整理原始数据，如材料性质、荷载条件、约 束条件等，离散结构并进行单元编码、结点 编码、结点位移编码、选取坐标系。2单元分析，建立单元刚度矩阵。3整体分析，建立总刚矩阵。4建立整体结构的等效节点荷载和总荷载矩阵5边界条件处理。6解方程，求出节点位移。7求出各单元的单元应力。8计算结果整理。计算结果整理包括位移和应 力两个方面；位移计算结果一般不需要特别 的处理，利用计算出的节点位移分量，就可 画出结构任意方向的位移云图；而应力解的 误

3、差表现在单元内部不满足平衡方程，单元 与单元边界处应力一般不连续，在边界上应 力解一般与力的边界条件不相符合。1.1 程序说明%*三角形常应变单元求解结构主程序%* 功能：运用有限元法中三角形常应变单元解平面问题的计算主程序。 根本思想：单元结点按右手法那么顺序编号。 荷载类型：可计算结点荷载。说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算 法的全过程，即实现程序流程图的程序表达。%1 程序准备format short e%设定输出类型clear all%去除所有已定义变量clc%清屏说明：format short e 设定计算过程中显示在屏幕上的数字类型为短格式、科 学计

4、数法；clear all 去除所有已定义变量，目的是在本程序的运行过程中，不会 发生变量名相同等可能使计算出错的情况；clc 清屏，使屏幕在本程序运行开始时%2 全局变量定义globalNNODENPIONNELEMNVFIXNFORCE COORDLNODSYOUNGPOISSTHICKglobalFORCEFIXEDglobalBMATXDMATXSMATXAREAglobalASTIFASLODASDISPglobalFP1说明：NNODE 单元结点数，NPION 总结点数， NELEM 单元数，NVFIX 受约束边界点数，NFORCE 结点力数，COORD结构结点坐标数组，LNODS

5、单元定义数组，YOUNG 弹性模量，POISS泊松比，THICK 厚度FORCE 节点力数组(n,3) n:受力节点数目，(n,1):作用点,(n,2):x方向,(n,3):y 方向；FIXED 约束信息数组(n,3) n:受约束节点数目，(n,1):约束点(n,2)与仇3) 分别为约束点 x 方向和 y 方向的约束情况 ,受约束为 1 否那么为 0BMATX 单元应变矩阵 (3*6)， DMATX 单元弹性矩阵 (3*3)，SMATX 单元应力矩阵 (3*6)，AREA 单元面积ASTIF 总体刚度矩阵，ASLOD 总体荷载向量，ASDISP结点位移向量FP1数据文件指针3 翻开文件FP1=

6、fopen(input.txt,rt);%翻开输入数据文件存放初始数据说明：FP1=fopen(input.txt,rt);-翻开已存在的输入数据文件 input.txt，且设置其 为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1 来执行FP2=fopen(output.txt,wt); -翻开输出数据文件， 该文件不存在时， 通过此 命令创立新文件，该文件存在时那么将原有内容全部删除。 该文件设置为可写格式， 可在程序执行过程中向输出文件写入数据。% 4 读入程序控制信息%结点个数结点编码总数%单元个数单元编码总数结点荷载个数弹性模量泊松比%厚度%单元定义数组单元结

7、点号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矩阵 显

8、然，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=fsca nf(FP1,%f,3,NFORCE)% 结点力数组说明：(n,3) n:受力结点数目,(n,1):作用点,(n,2):x方向,(n,3):y方向FIXED=fsca nf(FP1,%d,3,i nf)% 约束信息数组说明：(n,3) n:受约束

9、节点数目，(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 表示正无穷。%-nteJ5调用子程生成单刚，

10、组成总刚并参加约束信息 ASSEMBLE。%6调用子程生成荷载向量FORMLOAD()%7计算结点位移向量ASDISP=ASTIFASLOD%8调用子程计算单元应力WRITESTRESS()Q%*9关闭输出数据文件fclose(FP2);读取ASSEMBLE子程%*fun ctio n ASSEMBLE()% 所引用的全局变量：global NPION NELEM NVFIX LNODS ASTIF THICK global BMATX SMATX AREA FIXED%计算单刚并生成总刚ASTIF(1:2*NPION,1:2*NPION)=O;%张成特定大小总刚矩阵并置 0说明：Array

11、stiffness建立单元刚度矩阵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*

12、2);%跟据节点编号对应关系将单元 刚度分块叠加到总刚矩阵中endendend说明：FORMSMATX(i)调用应力子程序，提取i单元的应力矩阵SMATX ； a=LNODS(i,:)记录i单元的三个结点编号；forend循环语句表示行从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即将单元中的位移编码对应到整体的位移编码。%

13、NVFIX :受约束边界点数%将约束信息参加总刚置0置1法 NUM=1; %计数器(当前已分析的节点数) i=0;%计数器(当前已处理的约束数)tmp(NVFIX)=0; %临时存被约束的序号 while ivNVFIXfor j=-1:0%假设发现约束%计数器自增匚temporary%求约束序号 if FIXED(NUM,j+3)=1i=i+1;tmp(i)=FIXED(NUM)*2+j;endendNUM=NUM+1;end说明： tmp(NVFIX)=0，形成一个元素值均为 0的一行，NVFIX列的行向量， 执行while语句，首先判断i是否小于控制数据NVFIX，假设小于那么往下进 行

14、，假设不小于那么退出。执行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:NVFIXASTIF(1:2*NPION,tmp(i)=0;%将一约束序号处的总刚列向量清0ASTIF(tmp(i),1:2*NPION)=0;

15、%.将一约束序号处的总刚行向量清0ASTIF(tmp(i),tmp(i)=1;%将行列交叉处的元素置为1end说明：0丨，那么将总刚中的主元素后处理法中置0置1法，设j Cj包括CjKjj换为1, j行和j列的其他元素均改为零。%(*%读取FORMSMATX子程%(*fun ctio n FORMSMATX(ELEMENT)% 计算应力矩阵%引用所需的全局变量global NPION NELEM COORD LNODS YOUNG POISS global BMATX DMATX SMATX AREA%生成弹性矩阵Da=YOUNG/(1-POISSA2);DMATX(1,1)=1*a;DMAT

16、X(1,2)=POISS*a;DMATX(2,1)=POISS*a;DMATX(2,2)=1*a;DMATX(3,3)=(1-POISS)*a/2;说明：1, ,0 0,0平面应力问题的弹性矩阵D,1,0，其中，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;说明：1, xi,

17、yidet表示求矩阵行列式的值,1,Xj ,yj其中xi(i, j,m)分别表示一个三1, Xm, ym角形单元的三个节点坐标。MATLAB中假设一行中无法写下一个完整的命令，那么可以在行尾参加3个连续的英文句号，表示命令余下的局部在下一行出现。%生成应变矩阵Bfor j=0:2b(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);endBMATX=b(1

18、) 0b(2)0b(3)0;.0c(1)0c(2)0c(3);.c(1) b(1) c(2)b(2) c(3)b(3)/(2*AREA);说明：bl,0、 1应变矩阵 Bl0, c (I i, j, m)2ACl ,blrem表示求余函数，remx, y命令返回的是x-n.*y，当y 0时，n=fix(x./y)， 其中fix为最近的整数取整。%SMATX=DMATX*BMATX;%求应力矩阵 S=D*BQ%*%读取FORMLOAD子程Q%*fun ctio n FORMLOAD()%本子程生成荷载向量%所需引用的全局变量global ASLOD NPION NFORCE FORCEASLOD

19、(1:2*NPION)=0;%张成特定大小的 0 向量说明：ASLOD 为总体荷载向量，每个结点有 x， y 两个方向的结点力。%for i=1:NFORCEASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end 说明：F0RCE(i,1)为作用点，FORCE(i,2:3)为x,y方向的结点力。%将有约束处的荷载置为 0NUM=1;i=0;tmp(NVFIX)=0;while iNVFIXfor j=-1:0if FIXED(NUM,j+3)=1i=i+1; tmp(i)=FIXED(NUM)*2+j;endend NUM=NUM+1;endf

20、or i=1:NVFIXASLOD(tmp(i)=0;end说明： 置0置1法，同上。%*ASDISP=ASTIFASLOD%计算结点位移向量%*%读取 WRITESTRESS 子程%* function WRITESTRESS()%求解内力，即两个方向的正应力与一个剪应力x, y, xy )%所引用的全局变量global NELEM NPION SMATX ASDISP LNODS%ELEDISP(1:6)=0;%当前单元的结点位移向量说明：ui vi ujELEDIS 为单元的结点位移 ae jvjumvm%for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=

21、ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);endFORMSMATX(i) %调用子程求应力矩阵iSTRESS=SMATX*ELEDISP %求内力 end说明：通过循环， 依次从结构位移列阵中对号， 赋值给第 i 个单元的单元位移向量ELEDISP。%1.2程序应用举例-;1【例1】利用三角形三节点常应变单元程序计算 图2所示的结构。设弹性模量E 1，泊松比为0,厚度为1。%输入数据文件input.txt为:6 4 1 6 1.0e0 0.0 13 1 25 2 43 2 56 3 51 0 -11 1 02 1 04 1 15 0 16 0 1%说明：第一行：读入程序控制信息%结点个数结点编码总数%单元个数单元编码总数%结点荷载个数%受约束边界点数%弹性模量NPION=fsca nf(FP1,%d,1)NELEM=fsca nf(FP1,%d,1)NFORCE=fsca nf(FP1,%d,1)NVFIX=fsca nf(FP1,%d,1)YOUNG=fsca nf(FP1,%e,1)泊松比%厚度POISS=fscanf(FP1,%f,1)THICK=fscanf(FP1,%d,1) 第二、三、四、五行：读入单元连接信息：LNODS=fscanf(FP1,%d,3,NELEM);%单元定义数组，单元结点号，逆时