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

- 1三角形常应变单元程序的编制与使用有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好的适应性，而矩形或四边形单元较三节点三角形有更高的计算精度。Matlab 语言是进行矩阵运算的强大工具，因此，用 Matlab 语言编写有限元中平面问题的程序有优越性。本章将详细介绍如何利用 Matlab 语言编制三角形常应变单元的计算程序，程序流程图见图 1。有限元法中三节点三角形分析结构的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。2）单元分析，建立单元刚度矩阵。3）整体分析，建立总刚矩阵。4）建立整体结构的等效节点荷载和总荷载矩阵5）边界条件处理。6）解方程，求出节点位移。7）求出各单元的单元应力。8）计算结果整理。计算结果整理包括位移和应力两个方面；位移计算结果一般不需要特别的处理，利用计算出的节点位移分量，就可画出结构任意方向的位移云图；而应力解的误差表现在单元内部不满足平衡方程，单 图 1 程序流程图开始输入初始数据生成单刚集成总刚施加约束信息生成荷载向量边界条件处理计算结点位移计算单元应力计算结果整理结束- 2元与单元边界处应力一般不连续，在边界上应力解一般与力的边界条件不相符合。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 FIXED global BMATX DMATX SMATX AREAglobal ASTIF ASLOD ASDISPglobal FP1 说明：- 3NNODE单元结点数，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) %单元定义数组（单元结点号） 说明：- 4建立 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 调用子程生成荷载向量- 5function FORMLOAD()%-7 计算结点位移向量ASDISP=ASTIFASLOD%-8 调用子程计算单元应力function 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);%跟据节点编号对应关系将单元刚度分块叠加到总刚矩阵中- 6endendend 说明：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 即将单元中的位移编码对应到整体的位移编码。%-%将约束信息加入总刚（置 0 置 1 法）NUM=1; %计数器(当前已分析的节点数 )i=0; %计数器(当前已处理的约束数 )tmp(NVFIX)=0; %临时存被约束的序号while iNVFIXfor j=-1:0if FIXED(NUM,j+3)=1 %若发现约束i=i+1; %计数器自增tmp(i)=FIXED(NUM)*2+j; %求约束序号endendNUM=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 计算整体- 7约束序号，将序号放入 tmp 行向量中的 i 列。%-for i=1:NVFIXASTIF(1:2*NPION,tmp(i)=0; %将一约束序号处的总刚列向量清 0ASTIF(tmp(i),1:2*NPION)=0; %.将一约束序号处的总刚行向量清 0ASTIF(tmp(i),tmp(i)=1; %将行列交叉处的元素置为 1end 说明：后处理法中置 0 置 1 法，设 （包括 ） ，则将总刚中的主元素jjC0jK jj 换为 1，j 行和 j 列的其他元素均改为零。%*% 读取 FORMSMATX 子程%* function FORMSMATX(ELEMENT) %计算应力矩阵%引用所需的全局变量global NPION NELEM COORD LNODS YOUNG POISSglobal BMATX DMATX SMATX AREA%-%生成弹性矩阵 Da=YOUNG/(1-POISS2); 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; 说明：平面应力问题的弹性矩阵 ,其中，E 为弹性模量，21,0,2ED- 8为泊松比。%-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 表示求矩阵行列式的值， ，其中 分别表示一个mjiyxA,1,2),(jix三角形单元的三个节点坐标。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) 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); 说明：应变矩阵 ),(,0,21mjilbcABlll rem 表示求余函数，rem（x，y）命令返回的是 x-n.*y，当 y 0 时，n=fix(x./y)，其中 fix 为最近的整数取整。%-SMATX=DMATX*BMATX; %求应力矩阵 S=D*B- 9%*% 读取 FORMLOAD 子程%* function FORMLOAD() %本子程生成荷载向量%-%所需引用的全局变量global ASLOD NPION NFORCE FORCE%-ASLOD(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 说明：FORCE(i,1)为作用点,FORCE(i,2:3)为 x,y 方向的结点力。%-%将有约束处的荷载置为 0NUM=1; i=0; tmp(NVFIX)=0; while iNVFIXfor j=-1:0if FIXED(NUM,j+3)=1 i=i+1; tmp(i)=FIXED(NUM)*2+j; endendNUM=NUM+1;endfor i=1:NVFIXASLOD(tmp(i)=0; end 说明：置0置1法，同上。%*ASDISP=ASTIFASLOD %计算结点位移向量- 10%*% 读取 WRITESTRESS 子程%* function WRITESTRESS()%求解内力，即两个方向的正应力与一个剪应力（ )xy,%-%所引用的全局变量global NELEM NPION SMATX ASDISP LNODS%-ELEDISP(1:6)=0; %当前单元的结点位移向量 说明：ELEDIS 为单元的结点位移 mjievuvua%-for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2); endFORMSMATX(i) % %调用子程求应力矩阵iSTRESS=SMATX*ELEDISP %求内力end 说明：通过循环，依次从结构位移列阵中对号，赋值给第 i 个单元的单元位移向量 ELEDISP。%- 111.2 程序应用举例【例 1】利用三角形三节点常应变单元程序计算图 2 所示的结构。设弹性模量 ，泊松比为 0，厚度为 1。 %-1E-输入数据文件 input.txt 为：6 4 1 6 1.0e0 0.0 13 1 25 2 43 2 5 6 3 5 图 2 0.0 2.00.0 1.01.0 1.00.0 0.01.0 0.02.0 0.01 0 -11 1 02 1 04 1 15 0 16 0 1%-说明：第一行：读入程序控制信息3654- 12NPION=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,