shangzai.rar_finite element_三角单元_三角形单元_三角形有限元_三角有限元

PROGRAM SIMPLE DIMENSION ALOAD(150,2),NFIX(120,2),SL(300),DISPL(300) DIMENSION DISP(6),STRESS(260,3),REACT(2) INTEGER HBW COMMON /BB/LNODS(260,3),COORDS(150,2) COMMON /WORK/D(3,3),B(3,6),DB(3,6),ELSTIF(6,6) COMMON /AA/SS(7300,50) OPEN (7,FILE='IN.DAT',STATUS='OLD') OPEN (8,FILE='OUT.DAT',STATUS='UNKNOWN') OPEN (9,FILE='DRAW.DAT',STATUS='UNKNOWN') HBW=0 C READ IN DATA READ(7,*) MAXNEL,MAXNOD,MAXLOD,MAXFIX,NCNTRL IF(NCNTRL.EQ.1) THEN WRITE(8,'(/5X,''This is problem of plane stress''//)') ELSE WRITE(8,'(/5X,''This is problem of plane strain''//)') ENDIF WRITE(8,1002) MAXNEL,MAXNOD,MAXLOD,MAXFIX C INITIALIZE VECTORS AND MATRICES MAXVAR=MAXNOD*2 DO I=1,MAXNOD DO J=1,2 ALOAD(I,J)=0.0 NFIX(I,J)=0.0 ENDDO ENDDO DO I=1,MAXVAR SL(I)=0.0 DISPL(I)=0.0 DO J=1,50 SS(I,J)=0.0 ENDDO ENDDO DO I=1,MAXNEL DO J=1,3 STRESS(I,J)=0.0 ENDDO ENDDO C READ IN ELEMENT DEFINITIONS WRITE(8,1003) DO NNEL=1,MAXNEL READ(7,*) NEL,(LNODS(NEL,I),I=1,3) WRITE(8,1005) NEL,(LNODS(NEL,I),I=1,3) NIC1=LNODS(NEL,1) NIC2=LNODS(NEL,2) NIC3=LNODS(NEL,3) IF(IABS(NIC1-NIC2).GT.HBW) HBW=IABS(NIC1-NIC2) IF(IABS(NIC2-NIC3).GT.HBW) HBW=IABS(NIC2-NIC3) IF(IABS(NIC3-NIC1).GT.HBW) HBW=IABS(NIC3-NIC1) ENDDO HBW=(HBW+1)*2 WRITE(8,'(//1X,4HHBW=,I5)') HBW C READ IN NODAL CO-ORDINATES WRITE(8,1006) DO NOD=1,MAXNOD READ(7,*) NIC,(COORDS(NIC,I),I=1,2) WRITE(8,1008) NIC,(COORDS(NIC,I),I=1,2) ENDDO C READ IN APPLIED LOADS AND FIXITIES DO J=1,MAXLOD READ(7,*) NIC,(ALOAD(NIC,I),I=1,2) ENDDO DO J=1,MAXFIX READ(7,*) NIC,(NFIX(NIC,I),I=1,2) ENDDO WRITE(8,1009) DO I=1,MAXNOD WRITE(8,1012) I,(ALOAD(I,J),J=1,2),(NFIX(I,J),J=1,2) ENDDO C READ IN ELASTIC PROPERTIES READ(7,*) YM,PR,T WRITE(8,1014) YM,PR,T IF(NCNTRL.EQ.2) THEN YM=YM/(1.-PR*PR) PR=PR/(1.-PR) T=1.0 ENDIF KK=0 DO NEL=1,MAXNEL NIC1=LNODS(NEL,1) NIC2=LNODS(NEL,2) NIC3=LNODS(NEL,3) X1=COORDS(NIC1,1) Y1=COORDS(NIC1,2) X2=COORDS(NIC2,1) Y2=COORDS(NIC2,2) X3=COORDS(NIC3,1) Y3=COORDS(NIC3,2) AREA=(X2*Y3-X3*Y2-X1*(Y3-Y2)+Y1*(X3-X2))/2.0 IF(AREA.LT.1.0E-15) THEN WRITE(8,'(1X,''ERROR NO.'',4I8)') NEL,NIC1,NIC2,NIC3 WRITE(*,'(1X,''ERROR NO.'',4I8)') NEL,NIC1,NIC2,NIC3 KK=1 ENDIF ENDDO IF(KK.EQ.1) STOP C SET UP STRUCTURAL STIFFNESS MATRIX EC1=YM/(1-PR*PR) EC2=PR EC12=EC1*(1-EC2)/2.0 DO I=1,3 DO J=1,3 D(I,J)=0.0 ENDDO ENDDO D(1,1)=EC1 D(1,2)=EC1*EC2 D(2,1)=EC1*EC2 D(2,2)=EC1 D(3,3)=EC12 DO NEL=1,MAXNEL WRITE(*,'(1X,4HNEL=,I5)') NEL CALL SDB(NEL,AREA) DO I=1,6 DO J=1,I C=0.0 DO K=1,3 C=C+B(K,I)*DB(K,J) ENDDO ELSTIF(I,J)=C*AREA*T ELSTIF(J,I)=ELSTIF(I,J) ENDDO ENDDO C ADD IN ELEMENT STIFFNESS INTO STRUCTURAL MATRIX DO I=1,3 DO II=1,2 ISTRST=(LNODS(NEL,I)-1)*2+II IELEMT=(I-1)*2+II DO J=1,3 DO JJ=1,2 JSTRST=(LNODS(NEL,J)-1)*2+JJ JELEMT=(J-1)*2+JJ IF(ISTRST.LE.JSTRST) THEN JSTRST=JSTRST-ISTRST+1 SS(ISTRST,JSTRST)=SS(ISTRST,JSTRST) & +ELSTIF(IELEMT,JELEMT) ENDIF ENDDO ENDDO ENDDO ENDDO ENDDO C ADD IN APPLIED LOADS DO I=1,MAXNOD DO J=1,2 SL(2*I-2+J)=ALOAD(I,J)+SL(2*I-2+J) ENDDO ENDDO C APPLY FIXITIES I=0 DO NIC=1,MAXNOD DO NVB=1,2 I=I+1 IF(NFIX(NIC,NVB).NE.0) THEN SS(I,1)=SS(I,1)+1.0E+30 ENDIF ENDDO ENDDO C CARRY OUT GAUSSIAN REDUCTION M1=MAXNOD*2 M2=M1-1 DO I=1,M2 I1=I+1 I2=I+HBW-1 IF(I2.GT.M1) I2=M1 DO K=I1,I2 FACT=SS(I,K-I+1)/SS(I,1) NSF1=HBW-K+I DO J=1,NSF1 SS(K,J)=SS(K,J)-FACT*SS(I,J+K-I) ENDDO SL(K)=SL(K)-FACT*SL(I) ENDDO ENDDO C BACK SUBSTITUTE DO I=1,M1 II=M1-I+1 PIVOT=SS(II,1) SS(II,1)=0.0 DO J=1,HBW K=II+J-1 IF(K.LE.M1) THEN SL(II)=SL(II)-SS(II,J)*DISPL(K) ENDIF ENDDO DISPL(II)=SL(II)/PIVOT ENDDO WRITE(8,1018) DO I=1,MAXNOD WRITE(8,1019) I,DISPL(2*I-1),DISPL(2*I) ENDDO WRITE(8,1017) DO I=1,MAXNOD II=0 DO J=1,2 REACT(J)=0.0 IF(NFIX(I,J).NE.0) THEN III=(I-1)*2+J REACT(J)=-DISPL(III)*1.0E30 II=1 ENDIF ENDDO IF(II.EQ.1) WRITE(8,1025) I,REACT ENDDO C CALCULATE STRESSES WRITE(8,1020) DO NEL=1,MAXNEL CALL SDB(NEL,AREA) DO I=1,3 II=2*LNODS(NEL,I)-1 JJ=2*LNODS(NEL,I) DISP(2*I-1)=DISPL(II) DISP(2*I)=DISPL(JJ) ENDDO DO I=1,3 DO J=1,6 STRESS(NEL,I)=STRESS(NEL,I)+DB(I,J)+DISP(J) ENDDO ENDDO WRITE(8,1021) NEL,(STRESS(NEL,J),J=1,3) ENDDO CLOSE (7) CLOSE (8) CLOSE (9) STOP 1002 FORMAT(6X,'TOTAL NUMBER OF ELEMENTS=',I5 & /25X,'NODES=',I5/25X,'LOADS=',I5/22X,'FIXITIES=',I5) 1003 FORMAT(///5X,'ELEMENT DEFINITIONS'/5X,'ELEMENT', & 2X,'NODE1',2X,'NODE2',2X,'NODE3'/) 1005 FORMAT(5X,I5,3(2X,I5)) 1006 FORMAT(///5X,'NODAL CO-ORDINATES'/5X,'NODE',6X,'X',9X,'Y'/) 1008 FORMAT(5X,I4,2F10.3) 1009 FORMAT(///5X,'APPLIED LOADS AND FIXITIES'/5X,'NODE',5X, & 'APPLIED LOADS',8X,'FIXITIES'/) 1012 FORMAT(5X,I4,3X,2E10.3,2I5) 1014 FORMAT(///5X,'MATERIAL PROPERTIES'//5X,'YOUNGS MODULUS FOR', & 'MATERIAL=',E10.3/5X,'POISSONS RATIO FOR MATERIAL=',E10.3/ & 5X,'CONSTANT ELEMENT THICKNESS=',E10.3,//5X,'DATA COMPLETE') 1017 FORMAT(//5X,'NODAL REACTIONS'//6X,4HNODE,6X,6HX COMP,6X,6HY COMP/) 1018 FORMAT(///5X,'NODAL DISPLACEMENTS'///19X,6HX COMP,9X,6HY COMP,//) 1019 FORMAT(5X,I5,5X,F10.5,5X,F10.5) 1020 FORMAT(///5X,'ELEMENT STRESSES'///21X,'SIGMA X-X',11X, & 'SIGMA Y-Y',13X,'TAU X-Y'/) 1021 FORMAT(5X,I5,3(5X,F15.5)) 1025 FORMAT(I10,2(2X,E10.3)) END SUBROUTINE SDB(NEL,AREA) COMMON /BB/LNODS(260,3),COORDS(150,2) COMMON /WORK/D(3,3),B(3,6),DB(3,6),ELSTIF(6,6) DO I=1,3 DO J=1,6 B(I,J)=0.0 ENDDO ENDDO NIC1=LNODS(NEL,1) NIC2=LNODS(NEL,2) NIC3=LNODS(NEL,3) X1=COORDS(NIC1,1) Y1=COORDS(NIC1,2) X2=COORDS(NIC2,1) Y2=COORDS(NIC2,2) X3=COORDS(NIC3,1) Y3=COORDS(NIC3,2) AREA=(X2*Y3-X3*Y2-X1*(Y3-Y2)+Y1*(X3-X2))/2.0 B(1,1)=(Y2-Y3)/(2.0*AREA) B(1,3)=(Y3-Y1)/(2.0*AREA) B(2,4)=(X1-X3)/(2.0*AREA) B(2,2)=(X3-X2)/(2.0*AREA) B(1,5)=(Y1-Y2)/(2.0*AREA) B(2,6)=(X2-X1)/(2.0*AREA) DO I=1,3 B(3,2*I-1)=B(2,2*I) B(3,2*I)=B(1,2*I-1) ENDDO DO I=1,3 DO J=1,6 DB(I,J)=0.0 DO K=1,3 DB(I,J)=DB(I,J)+D(I,K)*B(K,J) ENDDO ENDDO ENDDO END