1、平面三角形结点有限元程序1、程序名：FEMT3.FOR ，FEMT3.EXE2、程序功能本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种荷载的作用，在计算自重时 y 轴取垂直向上为正；能处理非零已知位移，如支座沉降的作用。主要输出的内容包括：结点位移、单元应力、主应力、第一主应力与 x 轴的夹角以及约束结点的支座反力。程序采用 Visual Fortran 5.0 编制而成，输入数据全部采用自由格式。3、程序流程及框图图-1 程序流程图输入数据信息形成 K形成 R计算，输出应力解 KR，输出 启 动停 机图-2 程序框图其中，各子程序的功能如下：INPUT输入结点坐
2、标、单元信息和材料参数；MR形成结点自由度序号矩阵；FORMMA形成指标矩阵 MA（N ）并调用其他功能子程序，相当于主控程序；DIV取出单元的 3 个结点号码和该单元的材料号并计算单元的 bi，c i 等；MGK形成整体劲度矩阵并按一维存放在 SK（NH）中；LOAD形成整体结点荷载列阵 F；OUTPUT输出结点位移或结点荷载；TREAT由于有非零已知位移，对 K 和 F 进行处理；DECOMP整体劲度矩阵的分解运算；FOBA前代、回代求出未知结点位移 ；ERFAC计算约束结点的支座反力；KRS计算单元劲度矩阵中的子块 Krs。4、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件
3、的名字。在 运 行 程 序 之 前 ， 必 须 根 据 程 序 中 输 入 要 求 建 立 一 个 存 放 原 始 数 据 的 文 件 ， 这 个文 件 的 名 字 由 少 于 8 个 字 符 或 数 字 组 成 。 数 据 文 件 包 括 如 下 内 容 ： 总控信息，共一条，9 个数据NP，NE，NM，NR ，NI ，NL，NG，ND，NCNP结点总数；NE单元总数；NM材料类型总数；NR约束结点总数；NI问题类型标识， 0 为平面应力问题， 1 为平面应变问题；NL受荷载作用的结点的数目；NG考虑自重作用为 1，不计自重为 0；ND非零已知位移结点的数目；NC要计算支座约束反力的结点数目
4、。 材料信息，共 NM 条，每条依次输入EO，VO ，W ，tEO弹性模量（kN/m 2） ；VO泊松比；W材料容重（kN/m 3） ；t单元厚度（m） 。这些信息都存放在数组 AE（ 4，NM）中。 坐标信息，共 NP 条，每条依次输入IP，X，YIP结点号；X，Y分别为结点的 x 坐标和 y 坐标。坐标信息存放在数组 X(2，NP)中。 单元信息，共 NE 条，每条依次输入JE，L ，Io ，Jo，MoJE单元号；L为该单元的材料类型号。Io，Jo，Mo分别为该单元 i， j， m 的整体编码。单元信息存放在数组 MEO(4,NE)中。 约束信息共 NR 条，每条依次输入一个数 IP xy
5、IP结点号；Ix，I y分别为该结点的约束情况，如果方向受约束时填 0，如果自由则填 1。 荷载信息，共 NL 条，每条依次输入 yIP，F x，F yIP结点号；Fx，F y分别为该结点的 x，y 方向的荷载分量（kN ） 。结点号存放在数组 NF（NL）中，结点荷载分量存放在数组 FV（2，NL）中。 若 ND0，输入非零已知位移信息，共 ND 条，每条依次输入IP，u x，u yIP结点号；ux， uy分别为该结点 x， y 方向已知位移分量（m ） ，若其中某方向为自由，则其相应分量为 0。结点号存放在数 NDI（ND）中，已知位移分量存放在数组 DV(2，ND) 中。 支座反力信息，
6、共 NC 条，每条依次输入IP，M1，M2 ， M3，M4IP 支座结点号；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，；MA(1000)主元素序号指
7、标矩阵；JR(2，500)结点自由度序号矩阵；ME(3)存放单元结点 i， j， m 的整体编码；NN(6)单元结点自由度序号；BI(3)，CI(3)单元劲度矩阵计算公式中的 bi，b j，b m 和 ci，c j，c m；S三角形单元的面积；H11，H 12，H 21，H 22单元劲度矩阵中子块 Krs 的 4 个元素。5、算例一个正方形弹性体，厚度为 1m，四边受单位均布法向力作用，由于对称性，取其 1/4进行计算，其有限元网格如图 2-3 所示，设 ， ，不考虑自27/1042mknE167.0重。该问题的精确解应力为 1 ， 1 ， 0。xkN/myxy图 -3 有限元网格（1）输入文
8、件数据6 4 1 5 0 3 0 0 52000.0 0.0 0.0 1.01 0.0 2.02 0.0 1.03 1.0 1.04 0.0 0.05 1.0 0.06 2.0 0.01 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 05 1 06 1 01 -0.5 -0.53 -1.0 -1.06 -0.5 -0.51 1 0 0 02 1 2 3 04 2 0 0 05 2 3 4 06 4 0 0 0（2）输出文件结果NODAL DISPLACEMENTSNODE X-COMP. Y-COMP.1 0.00000E+00 -
9、0.10000E-022 0.00000E+00 -0.50000E-033 -0.50000E-03 -0.50000E-034 0.00000E+00 0.00000E+005 -0.50000E-03 0.00000E+006 -0.10000E-02 0.00000E+00ELEMENT STRESSESELEMENT X-STRESS Y-STRESS XY-STRESS MAX-STRESS MIN-STRESS ANGLE1 -1.000 -1.000 0.000 -1.000 -1.000 90.0002 -1.000 -1.000 0.000 -1.000 -1.000 90
10、.0003 -1.000 -1.000 0.000 -1.000 -1.000 90.0004 -1.000 -1.000 0.000 -1.000 -1.000 90.000NODE STRESSESNODE X-STRESS Y-STRESS XY-STRESS MAX-STRESS MIN-STRESS ANGLE1 -1.000 -1.000 0.000 -1.000 -1.000 90.0002 -1.000 -1.000 0.000 -1.000 -1.000 90.0003 -1.000 -1.000 0.000 -1.000 -1.000 90.0004 -1.000 -1.0
11、00 0.000 -1.000 -1.000 90.0005 -1.000 -1.000 0.000 -1.000 -1.000 90.0006 -1.000 -1.000 0.000 -1.000 -1.000 90.000NODAL REACTIONSNODE X-COMP Y-COMP1 0.000 0.0002 1.000 0.0004 0.500 0.5005 0.000 1.0006 0.000 0.0006、源程序C FINITE ELEMENT PROGRAM FOR TWO DIMENSIONALC TRIANGLE ELEMENTC DIMENSION K(800000),
12、COOR(2,3000),AE(4,11),* MEL(5,2000),MA(6000)CHARACTER*32 datCOMMON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC WRITE(*,300)300 FORMAT(/ * :*/+PLEASE INPUT FILE NAME OF DATA)READ(*,*)dataOPEN (4,FILE=data,STATUS=OLD)OPEN (7,FILE=OUT,STATUS=UNKNOWN)READ (4,*) NP,NE,NM,NR,NI,NL,NG,ND,NC C WRITE (*,400) NP,NE,NM,NR
13、,NI,NL,NG,ND,NC C WRITE (7,400) NP,NE,NM,NR,NI,NL,NG,ND,NCCALL 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,R) CALL DECOMP (SK,MA) CALL FOBA (SK,MA,R) WRITE(*,650) WRITE(7,650)CALL OUTPUT(JR,R) WRITE(*,700)
14、 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 FORMAT(1X,TOTAL DEGREES OF FREEDOM N=,*I4,1X,TOTAL-STORAGE ,NH=,I5,1X,*MAX-SEMI-BANDWIDTH MX=,I3)550 FORMAT(/
15、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,ELEMENT STRESSES/5X,*ELEMENT,5X,X-STRESS,3X,Y-STRESS,*2X,XY-STRESS,1X,MAX-STRESS,1X,*MIN-STRESS,6X,ANGLE/)STOP EN
16、D 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) 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
17、(2,*),AE(4,*),MEL(3,*)COMMON /CA/ NP,NE,NM,NRCOMMON /CC/ N,MX,NHDO 70 I=1,NPREAD(4,*) IP,X,YCOOR(1,IP)=XCOOR(2,IP)=Y70 CONTINUEDO 11 J=1,NEREAD(4,*)NEE,NME,(MEL(I,NEE),I=1,3)MEL(3,NEE)=NME11 CONTINUEDO 10 I=1,NPDO 10 J=1,210 JR(J,I)=1DO 20 I=1,NRREAD(4,*) IP,IX,IYJR(1,IP)=IXJR(2,IP)=IY20 CONTINUEN=0
18、DO 30 I=1,NPDO 30 J=1,2IF (JR(J,I) 30,30,2525 N=N+1JR(J,I)=N30 CONTINUEDO 55 J=1,NMREAD (4,*)jj,(AE(I,jj),I=1,4)C WRITE(*,910)jj,(AE(I,jj),I=1,4)IF(NI.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) endif55 CONTINUE910 FORMAT (/20X,MATERIAL PROPERTIES/(3X,I5,4(1x,E
19、8.3)RETURNENDC *SUBROUTINE CBAND (MA,JR,MEL)DIMENSION MA(*),JR(2,*),MEL(3,*),NN(6)COMMON /CA/ NP,NE,NM,NRCOMMON /CC/ N,MX,NHDO 65 I=1,N65 MA(I)=0DO 90 IE=1,NEDO 75 K=1,3IEK=MEL(K,IE)DO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 CONTINUE75 CONTINUEL=NDO 80 I=1,6NNI=NN(I)IF(NNI.EQ.0) GO TO 80IF(NNI.LT.L)
20、L=NNI80 CONTINUEDO 85 M=1,6JP=NN(M)IF(JP.EQ.0) GO TO 85JPL=JP-L+1IF(JPL.GT.MA(JP) MA(JP)=JPL85 CONTINUE90 CONTINUEMX=0MA(1)=1DO 10 I=2,NIF(MA(I).GT.MX) MX=MA(I)MA(I)=MA(I)+MA(I-1)10 CONTINUENH=MA(N)WRITE (*,500) N,MX,NHWRITE (7,500) N,MX,NH500 FORMAT (/5X,FREEDOM N=*,I5,3X,SEMI-BANDWI. MX=,I5,3X,* S
21、TORAGE NH=,I7)RETURNEND C*SUBROUTINE SK0(SK,R,COOR,MEL,MA,JR,AE)DIMENSION AE(4,*),COOR(2,*),MEL(3,*),JR(2,*),R(*),*MA(*),SK(*),SKE(6,6),NN(6) COMMON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NCCOMMON /CB/ EO,VO,W,T,A,H11,H12,H21,H22,*ME(3),BI(3),CI(3) COMMON /CC/ N,NHDO 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
