下载文档原格式
三角形常应变单元程序的编制与使用有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。
有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。
对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好的适应性，而矩形或四边形单元较三节点三角形有更高的计算精度。
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'); －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。
三角形常应变单元程序的编制与使用有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。
有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。
对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好Array的适应性，而矩形或四边形单元较三节点三角形有更高的计算精度。
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'); －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。
用常应变三角形单元解弹性力学平面问题的程序******************************************************************** ANALYSIS PROGTAM OF FINITE ELEMENT METHOD ** FOR PLANE STRESS/STRAIN OF TRIANGULAR ELEMENT ** ----- FEMT3.FOR ----- **------------------------------------------------------------- ** Subroutines: 1-SDATA, 2-STE, 3-ATE, 4-DTE, 5-BTE, 6-STIFF ** 7-EQUPE, 8-INSCD, 9-BGSMT, 10-SIGME ********************************************************************DIMENSION LND(50,3),X(100),Y(100),JR(20,3),PJ(20,3),P(200)REAL KS(200,100)OPEN(5,FILE='FEMT3.DAT')OPEN(6,FILE='FEMT3.OUT',STATUS='NEW')READ(5,*) NJ,NE,NS,NPJ,IPS(结点、单元、支承、荷载、类型)WRITE(6,*)' FINITE ELEMENT ANALYSIS IN PLANE PROBLEM'WRITE(6,*)' SOURCE DATA OUTPUT'WRITE(6,20) NJ,NE,NS,NPJ,IPS20 FORMAT(4X,'NJ',3X,'NE',3X,'NS',3X,'NPJ',2X,'IPS'/1X,5I5)IF(IPS.EQ.0) WRITE(6,*)' PLANE STRESS PROBLEM'IF(IPS.EQ.1) WRITE(6,*)' PLANE STRAIN PROBLEM'CALL SDATA(NJ,NE,NS,NW,NPJ,IPS,E,PR,T,V,LND,X,Y,JR,PJ)NJ2=2*NJWRITE(6,50) NJ250 FORMAT(/1X,'DEGREES OF FREEDOM=',I5)WRITE(6,60) NW60 FORMAT(1X,'BAND WIDTH=',I5)CALL STIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS)(总刚6)CALL EQUPE(NJ,NE,NPJ,NJ2,T,V,LND,X,Y,PJ,P)({P}7)CALL INSCD(NS,NW,NJ2,JR,KS,P)(引入支承条件8)CALL BGSMT(NJ,NJ2,NW,KS,P)（解方程9）CALL SIGME(NE,NJ,NJ2,E,PR,LND,X,Y,P)（求应力10）CLOSE(5)CLOSE(6)END*--------------------------------------------------------C SUBPROGRAM-1C INPUT STRUCTURAL DATASUBROUTINE SDATA(NJ,NE,NS,NW,NPJ,IPS,E,PR,* T,V,LND,X,Y,JR,PJ)DIMENSION LND(NE,3),X(NJ),Y(NJ),JR(NS,3),PJ(NPJ,3)READ(5,*) E,PR,T,V(弹性模量、泊松比、厚度、容重)WRITE(6,10) E,PR,T,V10 FORMAT(/6X,'E',10X,'PR',9X,'T',9X,'V'/,4F10.2)READ(5,*)((LND(I,J),J=1,3),I=1,NE)（结点编码）WRITE(6,20)20 FORMAT(/1X,'ELEMENT INFORMATION'/3X,'ELEM',3X,* 'I J K'/)WRITE(6,30)(I,(LND(I,J),J=1,3),I=1,NE)30 FORMAT(1X,4I5)READ(5,*)(X(I),Y(I),I=1,NJ)（结点坐标）WRITE(6,40)40 FORMAT(/1X,'COORDINATES OF NODES'/3X,'NODES',* 8X,'X',13X,'Y')WRITE(6,50)(I,X(I),Y(I),I=1,NJ)50 FORMAT(1X,I5,2E15.6)READ(5,*)((JR(I,J),J=1,3),I=1,NS)（约束信息）WRITE(6,60)60 FORMAT(/1X,'CONSTRAINED NODES'/3X,'NODE',3X,'X',4X,'Y') WRITE(6,70)((JR(I,J),J=1,3),I=1,NS)70 FORMAT(1X,3I5)READ(5,*)((PJ(I,J),J=1,3),I=1,NPJ)（荷载信息）WRITE(6,80)80 FORMAT(/1X,'LOAD CASES'/3X,'NODE',8X,'X',13X,'Y')WRITE(6,90)((PJ(I,J),J=1,3),I=1,NPJ)90 FORMAT(1X,F5.0,2E15.6)100 NW=0（半带宽）DO 110 IE=1,NEDO 110 I=1,3DO 110 J=1,3IW=IABS(LND(IE,I)-LND(IE,J))IF(NW.LT.IW) THENNW=IWENDIF110 CONTINUENW=(NW+1)*2IF(IPS.NE.0) THENE=E/(1.0-PR*PR)PR=PR/(1.0-PR)ENDIFEND*---------------------------------------------------------C SUBPROGRAM-2C CALCULATE ELEMENT STIFFNESS MATRIXSUBROUTINE STE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)DIMENSION LND(NE,3),X(NJ),Y(NJ),B(3,6),D(3,3)REAL KE(6,6)CALL ATE(IE,NJ,NE,LND,X,Y,AE)CALL DTE(E,PR,D)CALL BTE(IE,NJ,NE,LND,X,Y,AE,B)DO 10 I=1,6DO 10 J=1,6KE(I,J)=0.DO 10 K=1,3DO 10 K1=1,310 KE(I,J)=KE(I,J)+B(K,I)*D(K,K1)*B(K1,J)C=AE*TDO 30 I=1,6DO 30 J=1,630 KE(I,J)=KE(I,J)*CEND*------------------------------------------------ C SUBPROGRAM-3C CALCULATE ELEMENT AREASUBROUTINE ATE(IE,NJ,NE,LND,X,Y,AE)DIMENSION LND(NE,3),X(NJ),Y(NJ)I=LND(IE,1)J=LND(IE,2)K=LND(IE,3)XIJ=X(J)-X(I)YIJ=Y(J)-Y(I)XIK=X(K)-X(I)YIK=Y(K)-Y(I)AE=.5*(XIJ*YIK-XIK*YIJ)END*---------------------------------------------- C SUBPROGRAM-4C CALCULATE ELASTICITY MATRIXSUBROUTINE DTE(E,PR,D)DIMENSION D(3,3)DO 10 I=1,3DO 10 J=1,310 D(I,J)=0.D(1,1)=E/(1.-PR*PR)D(1,2)=E*PR/(1.-PR*PR)D(2,1)=D(1,2)D(2,2)=D(1,1)D(3,3)=.5*E/(1.+PR)END*------------------------------------------------ C SUBPROGRAM-5C CALCULATE MATRIX [B]SUBROUTINE BTE(IE,NJ,NE,LND,X,Y,AE,B)DIMENSION LND(NE,3),X(NJ),Y(NJ),B(3,6)I=LND(IE,1)J=LND(IE,2)K=LND(IE,3)DO 10 II=1,3DO 10 JJ=1,610 B(II,JJ)=0.B(1,1)=Y(J)-Y(K)B(1,3)=Y(K)-Y(I)B(1,5)=Y(I)-Y(J)B(2,2)=X(K)-X(J)B(2,4)=X(I)-X(K)B(2,6)=X(J)-X(I)B(3,1)=B(2,2)B(3,2)=B(1,1)B(3,3)=B(2,4)B(3,4)=B(1,3)B(3,5)=B(2,6)B(3,6)=B(1,5)DO 20 I1=1,3DO 20 J1=1,620 B(I1,J1)=.5/AE*B(I1,J1)END*------------------------------------------------------- C SUBPROGRAM-6C CALCULATE GLOBAL STIFFNESS MATRIXSUBROUTINE STIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS) DIMENSION LND(NE,3),X(NJ),Y(NJ)REAL KS(NJ2,NW),KE(6,6)DO 5 I=1,NJ2DO 5 J=1,NW5 KS(I,J)=0.DO 10 IE=1,NECALL STE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)DO 10 I=1,3IZ=LND(IE,I)DO 10 II=1,2IH =2*(I -1)+IIIDH=2*(IZ-1)+IIDO 10 J=1,3JZ=LND(IE,J)DO 10 JJ=1,2L =2*(J -1)+JJIL=2*(JZ-1)+JJIF(IL.GE.IDH) THENIDL=IL-IDH+1KS(IDH,IDL)=KS(IDH,IDL)+KE(IH,L)ENDIF10 CONTINUEEND*-------------------------------------------------------- C SUBPROGRAM-7C CALCULATE NODAL LOAD VECTORSUBROUTINE EQUPE(NJ,NE,NPJ,NJ2,T,V,LND,X,Y,PJ,P) DIMENSION LND(NE,3),X(NJ),Y(NJ),PJ(NPJ,3),P(NJ2) DO 10 I=1,NJ210 P(I)=0.DO 20 I=1,NPJII=PJ(I,1)P(2*II-1)=PJ(I,2)20 P(2*II)=PJ(I,3)30 IF(V.EQ.0.) GOTO 50DO 40 IE=1,NECALL ATE(IE,NJ,NE,LND,X,Y,AE)PE=-V*AE*T/3.DO 40 I=1,3II=LND(IE,I)40 P(2*II)=P(2*II)+PE50 RETURNEND*---------------------------------------------C SUBPROGRAM-8C INTRODUCE BOUNDARY CONDITIONSUBROUTINE INSCD(NS,NW,NJ2,JR,KS,P)DIMENSION P(NJ2),JR(NS,3)REAL KS(NJ2,NW)DO 30 I=1,NSIR=JR(I,1)DO 30 J=2,3IF(JR(I,J).EQ.0) GOTO 30II=2*IR+J-3KS(II,1)=1.DO 10 JJ=2,NW10 KS(II,JJ)=0.IF(II.GT.NW) JO=NWIF(II.LE.NW) JO=IIDO 20 JJ=2,JO20 KS(II-JJ+1,JJ)=0.P(II)=0.30 CONTINUEEND*-------------------------------------------C SUBPROGRAM-9C SOLVE EQUATIONS BY GS METHODSUBROUTINE BGSMT(NJ,NJ2,NW,KS,P)DIMENSION P(NJ2)REAL KS(NJ2,NW)NJ1=NJ2-1DO 20 K=1,NJ1IF(NJ2.GT.K+NW-1) IM=K+NW-1IF(NJ2.LE.K+NW-1) IM=NJ2K1=K+1DO 20 I=K1,IML=I-K+1C=KS(K,L)/KS(K,1)IW=NW-L+1DO 10 J=1,IWM=J+I-K10 KS(I,J)=KS(I,J)-C*KS(K,M)20 P(I)=P(I)-C*P(K)P(NJ2)=P(NJ2)/KS(NJ2,1)DO 40 I1=1,NJ1I=NJ2-I1IF(NW.GT.NJ2-I+1) JO=NJ2-I+1IF(NW.LE.NJ2-I+1) JO=NWDO 30 J=2,JOK=J+I-130 P(I)=P(I)-KS(I,J)*P(K)40 P(I)=P(I)/KS(I,1)WRITE(6,50)50 FORMAT(/1X,'NODAL DISPLACEMENTS'/3X,* 'NODE',5X,'X-DISP.',8X,'Y-DISP.')DO 60 I=1,NJ60 WRITE(6,70) I,P(2*I-1),P(2*I)70 FORMAT(1X,I5,2E15.6)END*--------------------------------------------------- C SUBPROGRAM-10C CALCULATE ELEMENT STRESS MATRIXSUBROUTINE SIGME(NE,NJ,NJ2,E,PR,LND,X,Y,P)DIMENSION LND(NE,3),X(NJ),Y(NJ),D(3,3),B(3,6), * S(3,6),ST(3),P(NJ2),DE(6)WRITE(6,*)WRITE(6,*)' ELEMENT STRESSES'CALL DTE(E,PR,D)DO 50 IE=1,NECALL ATE(IE,NJ,NE,LND,X,Y,AE)CALL BTE(IE,NJ,NE,LND,X,Y,AE,B)DO 10 I=1,3DO 10 J=1,6S(I,J)=0.DO 10 K=1,310 S(I,J)=S(I,J)+D(I,K)*B(K,J)DO 20 I=1,3DO 20 J=1,2IH=2*(I-1)+JIW=2*(LND(IE,I)-1)+J20 DE(IH)=P(IW)DO 30 I=1,3ST(I)=0.DO 30 J=1,630 ST(I)=ST(I)+S(I,J)*DE(J)SGX=ST(1)SGY=ST(2)TXY=ST(3)ASG=(SGX+SGY)*.5RSG=SQRT(.25*(SGX-SGY)**2+TXY*TXY)SGMA=ASG+RSGSGMI=ASG-RSGIF(SGY.EQ.SGMI) CETA=0.IF(SGY.NE.SGMI) CETA=90.-57.29578*ATAN* (TXY/(SGY-SGMI))50 WRITE(6,60) IE,SGX,SGY,TXY,SGMA,SGMI,CETA60 FORMAT(1X,'ELEMENT NO.=',I4/2X,'SIGX=',E10.4, * 2X,'SIGY=',E10.4,2X,'TXY =',E10.4/2X,'SGMA=', * E10.4,2X,'SGMI=',E10.4,2X,'CETA=',E10.4)END。
