平面三角形3节点有限元程序.doc

《平面三角形3节点有限元程序.doc》由会员分享，可在线阅读，更多相关《平面三角形3节点有限元程序.doc（17页珍藏版）》请在三一办公上搜索。

1、平面三角形结点有限元程序1、程序名：FEMT3.FOR，FEMT3.EXE2、程序功能本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种荷载的作用，在计算自重时y轴取垂直向上为正；能处理非零位移，如支座沉降的作用。主要输出的内容包括：结点位移、单元应力、主应力、第一主应力与x轴的夹角以与约束结点的支座反力。程序采用Visual Fortran 5.0编制而成，输入数据全部采用自由格式。3、程序流程与框图输入数据信息形成K形成R计算，输出应力解KR，输出启 动停 机图-1 程序流程图图-2 程序框图其中，各子程序的功能如下：INPUT输入结点坐标、单元信息和材料参数；M

2、R形成结点自由度序号矩阵；FORMMA形成指标矩阵MAN并调用其他功能子程序，相当于主控程序；DIV取出单元的3个结点和该单元的材料号并计算单元的bi，ci等；MGK形成整体劲度矩阵并按一维存放在SKNH中；LOAD形成整体结点荷载列阵F；OUTPUT输出结点位移或结点荷载；TREAT由于有非零位移，对K和F进展处理；DEP整体劲度矩阵的分解运算；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要计算支座约束反力的结点数目。材料信息，共NM条，每条依次输入EO，VO，W，tEO弹性模量kN/m2；VO泊松比；W材料容重kN/m3；t单元厚度m。信息都存放在数组AE4，NM中。坐标信息，共NP条，每条依次输入IP，X，YIP结点号；X，Y分别为结

4、点的x坐标和y坐标。坐标信息存放在数组X(2，NP)中。单元信息，共NE条，每条依次输入JE，L ，Io，Jo，MoJE单元号；L为该单元的材料类型号。Io，Jo，Mo分别为该单元i，j，m的整体编码。单元信息存放在数组MEO(4,NE)中。约束信息共NR条，每条依次输入一个数IP IP结点号；Ix，Iy分别为该结点的约束情况，如果方向受约束时填0，如果自由如此填1。荷载信息，共NL条，每条依次输入yIP，Fx，FyIP结点号；Fx，Fy分别为该结点的x，y方向的荷载分量kN。结点号存放在数组NFNL中，结点荷载分量存放在数组FV2，NL中。假如ND0，输入非零位移信息，共ND条，每条依次输入

5、IP，ux，uyIP结点号；ux，uy分别为该结点x，y方向位移分量m，假如其中某方向为自由，如此其相应分量为0。结点号存放在数NDIND中，位移分量存放在数组DV(2，ND)中。支座反力信息，共NC条，每条依次输入IP，M1，M2，M3，M4IP 支座结点号；M1，M2，M3，M4 为与该支座结点相关的单元号，假如不足4个，如此用0补充。支座结点号存放在数组NCI(NC)中，相关单元号存放在数组NCE(4，NC)中。以上数据须按如上顺序存放在数据文件中。除此之外，程序中还用到其他一些主要变量和数组，说明如下：N结构自由度总数；NH 按一维存贮的整体劲度矩阵的总容量；NX 最大半带宽；SK(1

6、0000)维存贮的劲度矩阵；R(1000)开始存放等效结点荷载，求解方程以后，用来存放结点位移；B(6)存放单元应力x，y，xy，1，2，；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所示，设，不考虑自重。该问

7、题的准确解应力为1，1，0。图-3 有限元网格1输入文件数据6 41 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 13 1 2 2 12 4 5 3 13 2 5 4 15 6 3101 201 4005106101 -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 02输出文件结果 NODAL DISPLACEMENTSNODE X-P. Y-P. 1 0.00000E

8、+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 STRESSESELEMENT 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.

9、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 STRESSESNODE 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.000 -1.000 90.000 3 -1.000 -1.000 0.000 -1.000 -1.000

10、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-P Y-P 1 0.000 0.000 2 1.000 0.000 4 0.500 0.500 5 0.000 1.000 6 0.000 0.0006、源程序C FINITE ELEMENT PROGRAM FOR TWO DIMENSIONALC TRIANGLE ELEMEN

11、TC DIMENSIONK(800000),COOR(2,3000),AE(4,11),* MEL(5,2000),MA(6000) CHARACTER*32 dat MON /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) OPEN (7,FILE=OUT,STATUS=UNKNOWN) READ (4,*) NP,NE,NM,NR,NI,NL,NG,ND,NC C

12、 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) CALL DEP (SK,MA) CALL FOBA (SK,MA,R) WRITE(*,650) WRITE(7,650) CALL OUTPUT(JR,R) WRITE(*,700) WRITE(7,7

13、00) CALL CES (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(/20X,TOTAL STORAGE IS *GREATER THAN 50000)600 FO

14、RMAT(30X,NODAL FORCES/8X,NODE, *11X,X-P.,14X,Y-P.)650 FORMAT(/30X,NODAL DISPLACEMENTS/8X, *NODE,13X,X-P.,12X,Y-P.)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 END C * SUBROUTINE KRS (BR,BS,CR,CS) MON /CB/ EO,VO,W,T

15、,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(2,*),AE(4,*),MEL(3,*) MON /CA/ NP,NE,NM,NR MON /CC

16、/ 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=1,3) MEL(3,NEE)=NME11 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)

17、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) 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,E8.3) RETURN ENDC * SUBROUTINE CBAND (MA,JR,ME

18、L) DIMENSION MA(*),JR(2,*),MEL(3,*),NN(6) MON /CA/ NP,NE,NM,NR MON /CC/ N,MX,NH DO 65 I=1,N65 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 CONTINUE75 CONTINUE L=N DO 80 I=1,6 NNI=NN(I) 80 CONTINUE DO 85 M=1,6 JP=NN(M) JPL=JP-L+1 85 CONTINUE 90 CONTINUE

19、MX=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,NH500 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,*)

20、,JR(2,*),R(*), *MA(*),SK(*),SKE(6,6),NN(6) MON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC MON /CB/ EO,VO,W,T,A,H11,H12,H21,H22, *ME(3),BI(3),CI(3) MON /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

21、*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 J2=ME(I) DO 40 J=1,2 J3=2*(I-1)+J NN(J3)=JR(J,J2) 40 CONTINUE DO 60 I=1,6 DO 60 J=1,6 JJ=NN(I) JK=NN(J) JL=MA(JJ) JM=JJ-JK JN=JL-JM SK(JN)=SK(JN)+SKE(I,J) 60 CONTINUE 70 CONTINUE c WRITE(0,500) (SK(I),I=1,20) 500 FORMAT(/10X,S

22、K=/(6F12.5) RETURN END C * SUBROUTINE LOAD (COOR,MEL,R,JR,AE) DIMENSION AE(4,*),COOR(2,*),MEL(3,*),R(*),JR(2,*), * NF(50),FV(2,50) MON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC MON /CB/ EO,VO,W,T,A,A1(4),ME(3),BB(6) MON /CC/ N,NH DO 10 I=1,N 10 R(I)=0.0 IF(NG) 70,70,30 30 DO 60 IE=1,NE CALL DIV (IE,COOR,MEL,A

23、E) DO 50 I=1,3 J2=ME(I) J3=JR(2,J2) IF(J3) 50,50,40 40 R(J3)=R(J3)-T*W*A/3.0 50 CONTINUE 60 CONTINUE 70 IF(NL) 110,110,80 80 READ(4,*) (NF(I),I=1,NL) READ(4,*) (FV(I,J),I=1,2),J=1,NL) C WRITE(*,500) (NF(I),I=1,NL) C WRITE(7,500) (NF(I),I=1,NL)C WRITE(*,600) (FV(I,J),I=1,2),J=1,NL) C WRITE(7,600) (FV

24、(I,J),I=1,2),J=1,NL) DO 100 I=1,NL JJ=NF(I) J=JR(1,JJ) M=JR(2,JJ) 100 CONTINUE 110 RETURN 500 FORMAT(/20X,NODES OF APPLIED LOAD*NF= */(1X,10I8)600 FORMAT(/30X,LUMPED-LOADS*FV= */(5X,5F15.3) END C * SUBROUTINE TREAT (SK,MA,JR,R) DIMENSION SK(*),MA(*),NDI(75),DV(2,75),JR(2,*),R(*) MON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC MON /CC/ N,NH READ(4,*) (NDI(J),J=1,ND) READ(4,*) (DV(I,J),I=1,2),J=1,ND) C WRITE(*,500) (NDI(J),J=1,ND) C WRITE(7,500) (NDI(J),J=1,ND)C WR