当前位置：文档之家 › 平面三角形单元有限元程序设计

平面三角形单元有限元程序设计

  • 格式：doc
  • 大小：244.00 KB
  • 文档页数：10

P9 m9 m

一、题目 如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖向载荷。已知：P=150N/m，E=200GPa，=0.25，t=0.1m，忽略自重。试计算薄板的位移及应力分布。

要求： 1. 编写有限元计算机程序，计算节点位移及单元应力。（划分三角形单元，单元数不得少于30个）； 2. 采用有限元软件分析该问题（有限元软件网格与程序设计网格必须一致），详细给出有限元软件每一步的操作过程，并将结果与程序计算结果进行对比（任选取三个点，对比位移值）； 3. 提交程序编写过程的详细报告及计算机程序； 4. 所有同学参加答辩，并演示有限元计算程序。

有限元法中三节点三角形分析结构的步骤如下： 1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。 2）单元分析，建立单元刚度矩阵。 3）整体分析，建立总刚矩阵。 4）建立整体结构的等效节点荷载和总荷载矩阵 5）边界条件处理。 6）解方程，求出节点位移。 7）求出各单元的单元应力。 8）计算结果整理。

一、程序设计 网格划分 如图，将薄板如图划分为6行，并建立坐标系，则

X Y

P X

Y

P

节点编号 单元编号 刚度矩阵的集成 建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。 由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。 通过循环逐个计算：（1）每个单元对应2种单元刚度矩阵中的哪一种； （2）该单元对应总刚度矩阵的那几行哪几列 （3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列 循环又分为3层循环：（1）最外层：逐行计算 （2）中间层：该行逐个计算 （3）最里层：区分为第 奇/偶 数个计算

单元刚度的集成：





''''''215656665656266256561661eZeeeZeZeeeekkkKkkkkkk





边界约束的处理：划0置1法 适用：这种方法适用于边界节点位移分量为已知(含为0)的各种约束。 做法： （1） 将总刚矩阵〔K〕中相应于已知位移行主对角线元素置1，其他元素改为零；同 时将载荷列阵｛R｝中相应元素用已知位移置换。 ◎ 这样，由该方程求得的此位移值一定等于已知量。 （2） 将〔K〕中已知位移相应的列的非主对角成元素也置0，以保持〔K〕的对称性。 ◎ 当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为 保证方程成立，须在方程右端减去已知位移对该方程的贡献——已知位移和相应总刚元素的乘积。 ◎ 若约束为零位移约束时，此步则可省去。 特点： （1） 经以上处理同样可以消除刚性位移(约束足够的前提下)，去掉未知约束反力。 （2） 但这种方法不改变方程阶数，利于存贮。 （3） 不过，若是要求出约束反力，仍要重新计算各个划去的总刚元素。

程序如下： 变量说明 NNODE 单元节点数 NPION 总结点数 NELEM 单元数 NVFIX 受约束边界点数 FIXED 约束信息数组 NFORCE 节点力数 FORCE 节点力数组 COORD 结构节点坐标数组 LNODS 单元定义数组 YOUNG 弹性模量 POISS 泊松比 THICK 厚度 B 单元应变矩阵(3*6) D 单元弹性矩阵(3*3) S 单元应力矩阵(3*6) A 单元面积 ESTIF 单元刚度矩阵 ASTIF 总体刚度矩阵 ASLOD 总体荷载向量 ASDISP 节点位移向量 ELEDISP 单元节点位移向量 STRESS 单元应力

%********************************************************** %初始化 clear format short e %设定输出类型 clear %清除内存变量 NELEM=36 %单元个数（单元编码总数） NPION=28 %结点个数（结点编码总数） NVFIX=2 %受约束边界点数 NFORCE=1 %结点荷载个数 YOUNG=2e11 %弹性模量 POISS=0.25 %泊松比 THICK=0.1 %厚度 LNODS=[1 2 3;2 4 5;2 5 3;3 5 6; 4 7 8;4 8 5;5 8 9;5 9 6; 6 9 10;7 11 12;7 12 8;8 12 13; 8 13 9;9 13 14;9 14 10;10 14 15; 11 16 17;11 17 12; 12 17 18; 12 18 13; 13 18 19; 13 19 14;14 19 20;14 20 15; 15 20 21;16 22 23;16 23 17;17 23 24; 17 24 18;18 24 25;18 25 19;25 19 26; 19 26 20;20 26 27;20 27 21;21 27 28] %单元定义数组（单元结点号） %相应为单元结点号（编码）、按逆时针顺序输入 COORD=[0 0;-0.75 1.5;0.75 1.5;-1.5 3;0 3; 1.5 3;-2.25 4.5;-0.75 4.5;0.75 4.5; 2.25 4.5;-3 6;-1.5 6;0 6;1.5 6;3 6; -3.75 7.5;-2.25 7.5; -0.75 7.5;0.75 7.5; 2.25 7.5;3.75 7.5;-4.5 9;-3 9; -1.5 9;0 9;1.5 9;3 9;4.5 9] %结点坐标数组 %坐标：x，y 坐标（共 NPOIN 组） FORCE=[1 0 -15] %结点力数组（受力结点编号, x 方向,y 方向） FIXED=[22 1 1;28 1 1] %约束信息（约束点，x 约束，y 约束） %有约束为 1，无约束为 0 %********************************************************** %生成单元刚度矩阵并组成总体刚度矩阵 ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置 0 %********************************************************** for i=1:NELEM %生成弹性矩阵 D D= [1 POISS 0; POISS 1 0; 0 0 (1-POISS)/2]*YOUNG/(1-POISS^2) %********************************************************** %计算当前单元的面积 A=-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 %********************************************************** %生成应变矩阵 B for j=0:2 b(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); end B=[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*A); B1( :,:,i)=B; %********************************************************** %求应力矩阵 S=D*B S=D*B; ESTIF=B'*S*THICK*A; %求解单元刚度矩阵 a=LNODS(i,:); %临时向量,用来记录当前单元的节点编号 for j=1:3 for k=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 %********************************************************** %将约束信息加入总体刚度矩阵（对角元素改一法） for i=1:NVFIX if FIXED(i,2)==1 ASTIF(:,(FIXED(i,1)*2-1))=0; %一列为零 ASTIF((FIXED(i,1)*2-1),:)=0; %一行为零 ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1; %对角元素为 1 end

%********************************************************** %生成单元刚度矩阵并组成总体刚度矩阵

合集下载

平面问题的有限元分析及三角形单元的应用

平面问题的有限元分析及三角形单元的应用

第八章 平面问题的有限元分析及三角形单元的应用第一节 概述分析弹性力学平面问题时，最简单的单元式由三个结点组成的三角形单元。

当用以分析平面应力问题时，可将其视为三角板；当用以分析平面应变问题时，则可式为三棱柱。

各单元在结点处为铰结。

图8-1所示位移悬臂梁离散为三角形单元的组合体以矩阵形式列出弹性力学平面问题的基本量和基本方程。

谈形体所受体力分量可表示为[]Tyxy x p p p p p =⎥⎥⎦⎤⎢⎢⎣⎡= (8-1)所受面力分量可表示为[]Tyxy x p p p p p =⎥⎥⎦⎤⎢⎢⎣⎡= （8-2）体内任一点应力分量可表示为[]T xy y x τδδδ= （8-3）任一点的应变分量可表示为[]T xy y x γεεε= （8-4）任一点的位移分量可表示为[]Tv u =δ （8-5）弹性力学平面问题的几何方程的矩阵表达式为⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂+∂∂∂∂∂∂=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=x u y v y v x u xy y x εεεε （8-6） 平面应力问题的物理方程的矩阵表达式为⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡--=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡xy y x xy y x E γεεμμμμτσσ2100010112 （8-7）或简写成为 εσD = （8-8）式中⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡--=2100010112μμμμE D （8-9）称为弹性矩阵。

平面应变问题的物理方程也可写成式（8-8），但须将式（8-9）中的E 换成21μ-E，μ换成21μμ-，因此得出⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡-----+-=)1(22100011011)21)(1()1(22μμμμμμμμμE D （8-10）平衡微分方程及边界条件也可以用矩阵表示，但弹性力学有限元位移法中，通常用虚功方程代替平衡微分方程和应力边界条件。

虚功方程的矩阵表达式为⎰⎰⎰⎰⎰***=+tdxdy tds p f ptdxdy f T T σε （8-11）式中：[]Tv u f ***=，表示虚位移；[]Txyx x ****=γεεε，表示与虚位移相对应的虚应变。

有限元方法 第五章 平面三角形单元

有限元方法 第五章 平面三角形单元

bi 、ci ，aj 、bj 、
第27页，共88页。
cj 和am 、bm 、cm 分别是行列式2的第一行、第二行和第 三行各元素的代数余子式，我们有
⒈ 形函数在各单元节点上的值，具有“本点是1、它点
式中 ui、vi 是节点i在x轴和y轴方向的位移。
第14页，共88页。
在有限单元法中，虽然是用离散化模型来代替原来的连续体， 但每一个单元体仍是一个弹性体，所以在其内部依然是符合弹性
力学基本假设的，弹性力学的基本方程在每个单元内部同样 适用。
从弹性力学平面问题的解析解法中可知，如果弹性体内的 位移分量函数已知，则应变分量和应力分量也就确定了。但是， 如果只知道弹性体中某几个点的位移分量的值，那么就不能直 接求得应变分量和应力分量。因此，在进行有限元分析时，必 须先假定一个位移模式。由于在弹性体内，各点的位移变化情 况非常复杂，很难在整个弹性体内选取一个恰当的位移函数来 表示位移的复杂变化，但是如果将整个区域分割成许多小单元，
题，三角形单元是最简单、也是最常用的单元，在平面应力问题中，
单元为三角形板，而在平面应变问题中，则是三棱柱。
假设采用三角形单元，把弹性体划分为有限个互不重叠的三角形。这 些三角形在其顶点（即节点）处互相连接，组成一个单元集合体， 以替代原来的弹性体。同时，将所有作用在单元上的载荷（包括集 中载荷、表面载荷和体积载荷），都按虚功等效的原则移置到节点 上，成为等效节点载荷。由此便得到了平面问题的有限元计算模型，
所示。
第17页，共88页。
y
vj (Vj )
j
uj (Uj )
um (Um )
vi (Vi )
i ui (Ui )
m
um (Um )
o
x
有限元中三角形单元源程序

有限元中三角形单元源程序

附录：平面问题三角形单元源程序******************************************************************* * 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。

有限元教学程序及使用说明§A-1平面三角形3结点有限元程序

有限元教学程序及使用说明§A-1平面三角形3结点有限元程序

附录有限元教学程序及使用说明§A-1平面三角形３结点有限元程序1、程序名：，2、程序功能本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种荷载的作用，在计算自重时y轴取垂直向上为正；能处理非零已知位移，如支座沉降的作用。

主要输出的内容包括：结点位移、单元应力、主应力、第一主应力与x轴的夹角以及约束结点的支座反力。

程序采用Visual Fortran 编制而成，输入数据全部采用自由格式。

3、程序流程及框图图A-1 程序流程图图A-2 程序框图其中，各子程序的功能如下：INPUT——输入结点坐标、单元信息和材料参数；MR——形成结点自由度序号矩阵；FORMMA——形成指标矩阵MA（N）并调用其他功能子程序，相当于主控程序；DIV——取出单元的3个结点号码和该单元的材料号并计算单元的b i，c i等；MGK——形成整体劲度矩阵并按一维存放在SK（NH）中；LOAD——形成整体结点荷载列阵F；OUTPUT——输出结点位移或结点荷载；TREAT——由于有非零已知位移，对K和F进行处理；DECOMP——整体劲度矩阵的分解运算；FOBA——前代、回代求出未知结点位移 ；ERFAC——计算约束结点的支座反力；KRS——计算单元劲度矩阵中的子块K rs。

4、程序使用说明当程序开始运行时，按屏幕提示，键入数据文件的名字。

在运行程序之前，必须根据程序中输入要求建立一个存放原始数据的文件，这个文件的名字由少于8个字符或数字组成。

数据文件包括如下内容：⑴总控信息，共一条，9个数据NP，NE，NM，NR，NI，NL，NG，ND，NCNP——结点总数；NE——单元总数；NM——材料类型总数；NR ——约束结点总数；NI ——问题类型标识，0为平面应力问题，1为平面应变问题；NL ——受荷载作用的结点的数目；NG ——考虑自重作用为1，不计自重为0；ND ——非零已知位移结点的数目；NC ——要计算支座约束反力的结点数目。