《平面三角形3节点有限元程序.pdf》由会员分享，可在线阅读，更多相关《平面三角形3节点有限元程序.pdf（10页珍藏版）》请在三一文库上搜索。
1、平面三角形结点有限元程序 1、程序名： FEMT3.FOR ，FEMT3.EXE 2、程序功能 本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种 荷载的作用， 在计算自重时y轴取垂直向上为正； 能处理非零已知位移，如支座沉降的作用。 主要输出的内容包括：结点位移、 单元应力、 主应力、 第一主应力与x轴的夹角以及约束结 点的支座反力。 程序采用 Visual Fortran 5 .0 编制而成，输入数据全部采用自由格式。 3、程序流程及框图 图 -1 程序流程图 输入数据信息 形成 K 形成 R 计算，输出应力 解K R，输出 启 动 停 机 图 -2 程序框图 其
2、中，各子程序的功能如下： INPUT 输入结点坐标、单元信息和材料参数； MR 形成结点自由度序号矩阵； FORMMA 形成指标矩阵MA （N）并调用其他功能子程序，相当于主控程序； DIV 取出单元的3 个结点号码和该单元的材料号并计算单元的bi，ci等； MGK 形成整体劲度矩阵并按一维存放在SK（NH ）中； LOAD 形成整体结点荷载列阵F； OUTPUT 输出结点位移或结点荷载； TREAT 由于有非零已知位移，对K和F进行处理； DECOMP 整体劲度矩阵的分解运算； FOBA 前代、回代求出未知结点位移； ERFAC 计算约束结点的支座反力； KRS 计算单元劲度矩阵中的子块Kr
3、s。 4、输入数据及变量说明 当程序开始运行时，按屏幕提示，键入数据文件的名字。 在运行程序之前，必须根据程序中输入要求建立一个存放原始数据的文件，这个 文件的名字由少于8 个字符或数字组成。数据文件包括如下内容： 总控信息，共一条，9 个数据 NP，NE，NM ，NR，NI ，NL，NG，ND，NC NP 结点总数； NE 单元总数； NM 材料类型总数； NR 约束结点总数； NI 问题类型标识，0 为平面应力问题，1 为平面应变问题； NL 受荷载作用的结点的数目； NG 考虑自重作用为1，不计自重为0； ND 非零已知位移结点的数目； NC 要计算支座约束反力的结点数目。 材料信息，共
4、NM 条，每条依次输入 EO，VO ，W，t EO 弹性模量（kN/m 2） ； VO 泊松比； W 材料容重（kN/m 3） ； t 单元厚度（ m） 。 这些信息都存放在数组AE（4，NM ）中。 坐标信息，共NP 条，每条依次输入 IP，X，Y IP 结点号； X，Y 分别为结点的x坐标和y坐标。 坐标信息存放在数组X(2 ，NP)中。 单元信息，共NE 条，每条依次输入 JE， L ，Io，Jo，Mo JE 单元号； L 为该单元的材料类型号。 Io，Jo，Mo 分别为该单元i，j，m 的整体编码。 单元信息存放在数组MEO(4,NE) 中。 约束信息共NR 条，每条依次输入一个数 I
5、P xy IP 结点号； Ix，Iy 分别为该结点的约束情况，如果方向受约束时填0，如果自由则填1。 荷载信息，共NL 条，每条依次输入y IP，Fx，Fy IP 结点号； Fx，Fy 分别为该结点的x，y方向的荷载分量（kN） 。 结点号存放在数组NF（NL ）中，结点荷载分量存放在数组FV（2，NL）中。 若 ND 0，输入非零已知位移信息，共ND 条，每条依次输入 IP，ux，uy IP 结点号； ux，uy 分别为该结点x，y方向已知位移分量（m） ，若其中某方向为自由，则其相应 分量为 0。 结点号存放在数NDI （ND）中，已知位移分量存放在数组DV(2 ，ND) 中。 支座反力信
6、息，共NC 条，每条依次输入 IP，M1 ，M2 ，M3 ，M4 IP 支座结点号； M1，M2 ，M3，M4 为与该支座结点相关的单元号，若不足4 个，则用 0 补充。支座结 点号存放在数组NCI(NC) 中，相关单元号存放在数组NCE(4，NC)中。 以上数据须按如上顺序存放在数据文件中。除此之外， 程序中还用到其他一些主要变量 和数组，说明如下： N 结构自由度总数； NH 按一维存贮的整体劲度矩阵的总容量； NX 最大半带宽； SK(10000) 维存贮的劲度矩阵； R(1000) 开始存放等效结点荷载，求解方程以后，用来存放结点位移； B(6) 存放单元应力 x，y，xy，1，2，
7、； MA(1000) 主元素序号指标矩阵； JR(2， 500) 结点自由度序号矩阵； ME(3) 存放单元结点i，j，m 的整体编码； NN(6) 单元结点自由度序号； BI(3) ，CI(3) 单元劲度矩阵计算公式中的bi，bj，bm和ci，cj，cm； S 三角形单元的面积； H11，H12，H21， H22 单元劲度矩阵中子块Krs的 4 个元素。 5、算例 一个正方形弹性体，厚度为1m，四边受单位均布法向力作用，由于对称性，取其1/4 进行计算，其有限元网格如图2-3 所示，设 27 /1042mknE，167.0，不考虑自 重。该问题的精确解应力为 x1 2 kN / m， y 1
8、 2 kN / m， xy 0。 图 -3 有限元网格 （1）输入文件数据 6 4 1 5 0 3 0 0 5 2000.0 0.0 0.0 1.0 1 0.0 2.0 2 0.0 1.0 3 1.0 1.0 4 0.0 0.0 5 1.0 0.0 6 2.0 0.0 1 1 3 1 2 2 1 2 4 5 3 1 3 2 5 4 1 5 6 3 1 0 1 2 0 1 4 0 0 5 1 0 6 1 0 1 -0.5 -0.5 3 -1.0 -1.0 6 -0.5 -0.5 1 1 0 0 0 2 1 2 3 0 4 2 0 0 0 5 2 3 4 0 6 4 0 0 0 （2）输出文件结果
9、 NODAL DISPLACEMENTS NODE X-COMP. Y-COMP. 1 0.00000E+00 -0.10000E-02 2 0.00000E+00 -0.50000E-03 3 -0.50000E-03 -0.50000E-03 4 0.00000E+00 0.00000E+00 5 -0.50000E-03 0.00000E+00 6 -0.10000E-02 0.00000E+00 ELEMENT STRESSES ELEMENT X-STRESS Y-STRESS XY-STRESS MAX-STRESS MIN-STRESS ANGLE 1 -1.000 -1.000
10、 0.000 -1.000 -1.000 90.000 2 -1.000 -1.000 0.000 -1.000 -1.000 90.000 3 -1.000 -1.000 0.000 -1.000 -1.000 90.000 4 -1.000 -1.000 0.000 -1.000 -1.000 90.000 NODE STRESSES NODE X-STRESS Y-STRESS XY-STRESS MAX-STRESS MIN-STRESS ANGLE 1 -1.000 -1.000 0.000 -1.000 -1.000 90.000 2 -1.000 -1.000 0.000 -1.
11、000 -1.000 90.000 3 -1.000 -1.000 0.000 -1.000 -1.000 90.000 4 -1.000 -1.000 0.000 -1.000 -1.000 90.000 5 -1.000 -1.000 0.000 -1.000 -1.000 90.000 6 -1.000 -1.000 0.000 -1.000 -1.000 90.000 NODAL REACTIONS NODE X-COMP Y-COMP 1 0.000 0.000 2 1.000 0.000 4 0.500 0.500 5 0.000 1.000 6 0.000 0.000 6、源程序
12、 C FINITE ELEMENT PROGRAM FOR TWO DIMENSIONAL C TRIANGLE ELEMENT C DIMENSION K(800000),COOR(2,3000),AE(4,11), * MEL(5,2000),MA(6000) CHARACTER*32 dat COMMON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC WRITE(*,300) 300 FORMAT(/ * :*/+PLEASE INPUT FILE NAME OF DATA) READ(*,*)data OPEN (4,FILE=data,STATUS=OLD) OPE
13、N (7,FILE=OUT,STATUS=UNKNOWN) READ (4,*) NP,NE,NM,NR,NI,NL,NG,ND,NC C WRITE (*,400) NP,NE,NM,NR,NI,NL,NG,ND,NC C WRITE (7,400) NP,NE,NM,NR,NI,NL,NG,ND,NC CALL INPUT (JR,COOR,MEL,AE) CALL CBAND (MA,JR,MEL) CALL SK0(SK,R,COOR,MEL,MA,JR,AE) CALL LOAD (COOR,MEL,R,JR,AE) IF (ND.GT.0) CALL TREAT (SK,MA,JR
14、,R) CALL DECOMP (SK,MA) CALL FOBA (SK,MA,R) WRITE(*,650) WRITE(7,650) CALL OUTPUT(JR,R) WRITE(*,700) WRITE(7,700) CALL CES (COOR,MEL,JR,R,AE) IF(NC.GT.0) call ERFAC (COOR,MEL,JR,R,AE) 400 FORMAT (/2X,NP=,I3,2X,NE=,I3,2X,NM= *,I3,2X,NR=,I3,2X,NI=I3,2X,NL=,I3,2X, *NG=,I3,2X,ND=,I3,2X,NC=,I3) 500 FORMA
15、T(1X,TOTAL DEGREES OF FREEDOM N=, *I4,1X,TOTAL-STORAGE ,NH=,I5,1X, *MAX-SEMI-BANDWIDTH MX=,I3) 550 FORMAT(/20X,TOTAL STORAGE IS *GREATER THAN 50000) 600 FORMAT(30X,NODAL FORCES/8X,NODE, *11X,X-COMP.,14X,Y-COMP.) 650 FORMAT(/30X,NODAL DISPLACEMENTS/8X, *NODE,13X,X-COMP.,12X,Y-COMP.) 700 FORMAT(/30X,E
16、LEMENT STRESSES/5X, *ELEMENT,5X,X-STRESS,3X,Y-STRESS, *2X,XY-STRESS,1X,MAX-STRESS,1X, *MIN-STRESS,6X,ANGLE/) STOP END C * SUBROUTINE KRS (BR,BS,CR,CS) COMMON /CB/ EO,VO,W,T,A,H11,H12,H21,H22 *,ME(3),BI(3),CI(3) ET=EO*T/(1.0-VO*VO)/A/4.0 V=(1.0-VO)/2.0 H11=ET*(BR*BS+V*CR*CS) H12=ET*(VO*BR*CS+V*BS*CR)
17、 H21=ET*(VO*CR*BS+V*BR*CS) H22=ET*(CR*CS+V*BR*BS) RETURN END C * SUBROUTINE INPUT (JR,COOR,MEL,AE) DIMENSION JR(2,*),COOR(2,*),AE(4,*),MEL(3,*) COMMON /CA/ NP,NE,NM,NR COMMON /CC/ N,MX,NH DO 70 I=1,NP READ(4,*) IP,X,Y COOR(1,IP)=X COOR(2,IP)=Y 70 CONTINUE DO 11 J=1,NE READ(4,*)NEE,NME,(MEL(I,NEE),I=
18、1,3) MEL(3,NEE)=NME 11 CONTINUE DO 10 I=1,NP DO 10 J=1,2 10 JR(J,I)=1 DO 20 I=1,NR READ(4,*) IP,IX,IY JR(1,IP)=IX JR(2,IP)=IY 20 CONTINUE N=0 DO 30 I=1,NP DO 30 J=1,2 IF (JR(J,I) 30,30,25 25 N=N+1 JR(J,I)=N 30 CONTINUE DO 55 J=1,NM READ (4,*)jj,(AE(I,jj),I=1,4) C WRITE(*,910)jj,(AE(I,jj),I=1,4) IF(N
19、I.eq.1)then AE(1,jj)=AE(1,jj)/(1.0-AE(2,jj)*AE(2,jj) AE(2,jj)=AE(2,jj)/(1.0-AE(2,jj) endif 55 CONTINUE 910 FORMAT (/20X,MATERIAL PROPERTIES/(3X,I5,4(1x,E8.3) RETURN END C * SUBROUTINE CBAND (MA,JR,MEL) DIMENSION MA(*),JR(2,*),MEL(3,*),NN(6) COMMON /CA/ NP,NE,NM,NR COMMON /CC/ N,MX,NH DO 65 I=1,N 65
20、MA(I)=0 DO 90 IE=1,NE DO 75 K=1,3 IEK=MEL(K,IE) DO 95 M=1,2 JJ=2*(K-1)+M NN(JJ)=JR(M,IEK) 95 CONTINUE 75 CONTINUE L=N DO 80 I=1,6 NNI=NN(I) IF(NNI.EQ.0) GO TO 80 IF(NNI.LT.L) L=NNI 80 CONTINUE DO 85 M=1,6 JP=NN(M) IF(JP.EQ.0) GO TO 85 JPL=JP-L+1 IF(JPL.GT.MA(JP) MA(JP)=JPL 85 CONTINUE 90 CONTINUE MX
21、=0 MA(1)=1 DO 10 I=2,N IF(MA(I).GT.MX) MX=MA(I) MA(I)=MA(I)+MA(I-1) 10 CONTINUE NH=MA(N) WRITE (*,500) N,MX,NH WRITE (7,500) N,MX,NH 500 FORMAT (/5X,FREEDOM N= *,I5,3X,SEMI-BANDWI. MX=,I5,3X, * STORAGE NH=,I7) RETURN END C* SUBROUTINE SK0(SK,R,COOR,MEL,MA,JR,AE) DIMENSION AE(4,*),COOR(2,*),MEL(3,*),
22、JR(2,*),R(*), *MA(*),SK(*),SKE(6,6),NN(6) COMMON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC COMMON /CB/ EO,VO,W,T,A,H11,H12,H21,H22, *ME(3),BI(3),CI(3) COMMON /CC/ N,NH DO 10 I=1,NH 10 SK(I)=0.0 DO 70 IE=1,NE CALL DIV (IE,COOR,MEL,AE) DO 30 I=1,3 DO 30 J=1,3 CALL KRS (BI(I),BI(J),CI(I),CI(J) SKE(2*I-1,2*J-1)=H11 SKE(2*I-1,2*J)=H12 SKE(2*I,2*J-1)=H21 SKE(2*I,2*J)=H22 30 CONTINUE DO 40 I=1,3