平面四邊形四節(jié)點(diǎn)等參單元Fortran源程序

1、C *C * FINITE ELEMENT PROGRAM *C * FOR Two DIMENSIONAL ELASticity PROBLEM *C * WITH 4 NODE *C * PROGRAM ELASTICITY character*32 dat,cch DIMENSION SK(80000),COOR(2,300),AE(4,11),MEL(5,200), & WG(4),JR(2,300),MA(600),R(600),iew(30),STRE(3,200) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH COMMON

2、 /CMN3/ RF(8),SKE(8,8),NN(8) WRITE(*,*)'PLEASE ENTER INPUT FILE NAME' READ(*,'(A)')DAT OPEN(4,FILE=dat,STATUS='OLD') OPEN(7,FILE='OUT',STATUS='UNKNOWN') READ(4,*)NP,NE,NM,NR WRITE(7,'(A,I6)')'NUMBER OF NODE-NP=',np WRITE(7,'(A,I6)')'

3、;NUMBER OF ELEMENT-NE=',ne WRITE(7,'(A,I6)')'NUMBER OF MATERIAL-NM=',nm WRITE(7,'(A,I6)')'NUMBER OF surporting-NC=',Nr CALL INPUT (JR,COOR,AE,MEL) CALL CBAND (MA,JR,MEL) DO I=1,NH SK(I)=0.0 enddo CALL SK0(SK,MEL,COOR,JR,MA,AE) do I=1,N R(I)=0.0 enddopause 'aaa

4、'stop READ(4,*)NCP,NBE,iz WRITE(*,'(5i8)')NCP,NBE,iz WRITE(7,'(5i8)')NCP,NBE,iz IF(NCP.GT.0)CALL CONCR(NCP,R,JR) IF(NBE.GT.0) CALL BODYR(NBE,R,MEL,COOR,JR,AE) IF(iz.GT.0)then do jj=1,iz READ (4,*)Js,nse,(WG(I),I=1,4) read(4,*)(iew(m),m=1,nse) CALL FACER(iew,NSE,R,MEL,COOR,JR,WG)

5、enddo endif CALL DECOP (SK,MA) CALL FOBA (SK,MA,R) CALL OUTDISP(NP,R,JR) CALL STRESS (COOR,MEL,JR,AE,R,STRE) WRITE(7,'(A)')' PROGRAM SAFF HAS BEEN ENDED' WRITE(*,'(A)')' PROGRAM SAFF HAS BEEN ENDED' STOPc RETURN ENDC * SUBROUTINE INPUT (JR,COOR,AE,MEL) DIMENSION JR(2,

6、*),COOR(2,*),AE(4,*),MEL(5,*) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ 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,4) MEL(5,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)=

7、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) WRITE(*,910) JJ,(AE(I,JJ),I=1,4)55 CONTINUE910 FORMAT (/20X,'MATERIAL PROPERTIES'/(3X,I5,4(1x,E8.3) RETURN ENDC * SUBROUTINE CBAND (MA,JR,MEL

8、) DIMENSION MA(*),JR(2,*),MEL(5,*),NN(8) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH DO 65 I=1,N65 MA(I)=0 DO 90 IE=1,NE DO 75 K=1,4 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,2*4 NNI=NN(I) IF(NNI.EQ.0) GO TO 80 IF(NNI.LT.L) L=NNI 80 CONTINUE DO 8

9、5 M=1,2*4 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=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(7,'(A,I8)')'TOTAL DEGREES OF FREEDOM-N= ',N WRITE(7,'(A,I8)')'MAX-SEMI-B

10、ANDWIDTH-MX=',MX WRITE(7,'(A,I8)')'TOTAL-STORAGE-NH=',NH 500 FORMAT (/5X,'FREEDOM N=' *,I5,3X,'SEMI-BANDWI. MX=',I5,3X, * 'STORAGE NH=',I7) RETURN ENDC * SUBROUTINE SK0(SK,MEL,COOR,JR,MA,AE) DIMENSION SK(*),MEL(5,*),COOR(2,*),JR(2,*),MA(*), * AE(4,*),XYZ(2

11、,4),iven(4) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH COMMON /CMN3/ RF(8),SKE(8,8),NN(8) COMMON /CMN4/ NEE,NME COMMON /GAUSS/ RSTG(3),H(3) H(1)=0.5555555555555560 H(2)=0.8888888888888890 H(3)=H(1) RSTG(1)=-0.7745966692414830 RSTG(2)=0.00 RSTG(3)=-RSTG(1) DO 10 IE=1,NE NEE=IE NME=MEL(5,IE) DO 7

12、5 K=1,4 IEK=MEL(K,IE) iven(k)=IEK DO 95 M=1,2 JJ=2*(K-1)+M NN(JJ)=JR(M,IEK)95 XYZ(M,K)=COOR(M,IEK)75 CONTINUE CALL STIF(XYZ,AE,iven) DO 60 I=1,8 DO 60 J=1,8 II=NN(I) JJ=NN(J) IF (JJ.EQ.0).OR.(II.LT.JJ) GO TO 60 JN=MA(II)-(II-JJ) SK(JN)=SK(JN)+SKE(I,J) 60 CONTINUE 70 CONTINUE write(7,1111) (ske(i,j),

13、j=1,8),i=1,8)1111format(2x,8f12.2) 10 CONTINUE RETURN ENDC * SUBROUTINE STIF(XYZ,AE,iven) DIMENSION AE(4,*),DNX(2,4),XYZ(2,*),iven(*), * RJAC(2,2) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH COMMON /CMN3/ RF(8),SKE(8,8),NN(8) COMMON /CMN4/ NEE,NME COMMON /GAUSS/ RSTG(3),H(3) DO 40 I=1,8 RF(I)=0.

14、00 DO 30 J=1,8 SKE(I,J)=0.00 30 CONTINUE 40 CONTINUE E=AE(1,NME) U=AE(2,NME) GAMA=AE(3,NME) D1=E*(1.00-U)/(1.00+U)*(1.00-2.00*U) D2=E*U/(1.00+U)*(1.00-2.00*U) D3=E*0.50/(1.00+U) DO 120 I=1,4 II=2*(I-1) I1=II+1 I2=II+2 DO 115 J=1,4 JJ=2*(J-1) J1=JJ+1 J2=JJ+2 DXX=0 DXY=0 DYX=0 DYY=0 DO 99 IS=1,3 S=RST

15、G(IS) SH=H(IS) DO 98 IR=1,3 R=RSTG(IR) RH=H(IR) CALL FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE) DNIX=DNX(1,I) DNIY=DNX(2,I) DNJX=DNX(1,J) DNJY=DNX(2,J) DXX=DXX+DNIX*DNJX*DET*RH*SH DXY=DXY+DNIX*DNJY*DET*RH*SH DYX=DYX+DNIY*DNJX*DET*RH*SH DYY=DYY+DNIY*DNJY*DET*RH*SH 98 CONTINUE 99 CONTINUE SKE(I1,J1)=DXX*D1+

16、DYY*D3 SKE(I2,J2)=DYY*D1+DXX*D3 SKE(I1,J2)=DXY*D2+DYX*D3 SKE(I2,J1)=DYX*D2+DXY*D3115 CONTINUE120 CONTINUE RETURN ENDC * SUBROUTINE CONCR(NCP,R,JR) DIMENSION R(*),JR(2,*),XYZ(2) DO 100 I=1,NCP READ (4,*) IP,PX,PY XYZ(1)=PX XYZ(2)=PY DO 95 J=1,2 L=JR(J,IP) IF(L.EQ.0) GO TO 95 R(L)=R(L)+XYZ(J) 95 CONTI

17、NUE100 CONTINUE RETURN ENDC * SUBROUTINE BODYR(NBE,R,MEL,COOR,JR,AE) DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*), & AE(4,*),XYZ(2,4),iven(4) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH COMMON /CMN3/ RF(8),SKE(8,8),NN(8) COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2) COMMON /GAUSS/ RSTG(3),H(3) H(1)=1.

18、0 H(2)=1.0 RSTG(1)=-0.5773502691896260 RSTG(2)=-RSTG(1) DO 10 IE=1,NBE DO I=1,8 RF(I)=0.00 ENDDOc READ(4,*)NEE NEE=ie NME=MEL(5,NEE) GAMA=AE(3,NME) DO 75 K=1,4 IEK=MEL(K,NEE) iven(k)=iek DO 95 M=1,2 JJ=2*(K-1)+M NN(JJ)=JR(M,IEK)95 XYZ(M,K)=COOR(M,IEK)75 CONTINUE DO 99 IS=1,2 S=RSTG(IS) SH=H(IS) DO 9

19、8 IR=1,2 RR=RSTG(IR) RH=H(IR) CALL FUN8 (XYZ,RR,S,DET) DO 30 I=1,4 J=2*I RF(J)=RF(J)-FUN(I)*RH*SH*DET*GAMA 30 CONTINUE 98 CONTINUE 99 CONTINUE CALL ASLOAD (R) 10 CONTINUE RETURN ENDC * SUBROUTINE FACER(iew,NSE,R,MEL,COOR,JR,WG) DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),wg(*) * ,XYZ(2,4),iew(*),PR(2)

20、 COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH COMMON /CMN3/ RF(8),SKE(8,8),NN(8) COMMON /CMN4/ NEE,NME COMMON /GAUSS/ RSTG(3),H(3) H(1)=1.0 H(2)=1.0 RSTG(1)=-0.5773502691896260 RSTG(2)=-RSTG(1) nwf=0 nnf=0 ir=wg(1)+0.1 if(ir.eq.1)nwf=1 if(ir.eq.2)nnf=1 DO 510 IE=1,NSE DO I=1,8 RF(I)=0.00 ENDDO ne

21、e=iew(ie) DO 575 K=1,4 IEK=MEL(K,NEE) DO 595 M=1,2 JJ=2*(K-1)+M NN(JJ)=JR(M,IEK)595 XYZ(M,K)=COOR(M,IEK)575 CONTINUE IF(NWF.EQ.1) then GAMA=WG(2) Z0=WG(3) NSU=WG(4)+0.1 CALL SURLOD (NSU,XYZ,PR,Z0,GAMA,1) endif IF(NNF.EQ.1) then q=WG(2) NSU=WG(4)+0.1 do j=1,2 PR(J)=q enddo CALL SURLOD (NSU,XYZ,PR,Z0,

22、GAMA,2) endif CALL ASLOAD (R) 510 CONTINUE RETURN ENDC * SUBROUTINE SURLOD (NSU,XYZ,PR,Z0,GAMA,NSI) DIMENSION XYZ(2,*),RST(3),PR(2),KCRD(4),KFACE(2,4), & FVAL(4),NODES(2),FACT(4) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH COMMON /CMN3/ RF(8),SKE(8,8),NN(8) COMMON /CMN4/ NEE,NME COMMON /CMN5

23、/ FUN(4),PN(2,4),XJAC(2,2) COMMON /GAUSS/ RSTG(3),H(3) DATA KCRD/1,1,2,2/ DATA KFACE/1, 4, * 2, 3, * 1, 2, * 4, 3/ DATA FVAL/-1.00,1.00,-1.00,1.00/ FACT(1)=1.0 FACT(2)=-1.0 FACT(3)=-1.0 FACT(4)=1.0 FACTNUS=FACT(NSU) DO I=1,2 J=KFACE(I,NSU) NODES(I)=J ENDDO IF (NSI.EQ.1) THEN DO I=1,2 J=NODES(I) Z=Z0

24、-XYZ(2,J) PR(I)=0.00 IF (Z.GT.0.00) PR(I)=Z*GAMA ENDDO ENDIF ML=KCRD(NSU) IF(ML.EQ.1)MM=2 IF(ML.EQ.2)MM=1 RST(ML)=FVAL(NSU) DO 70 LX=1,2 RST(MM)=RSTG(LX) CALL FUN8 (XYZ,RST(1),RST(2),DET) PXYZ=0.00 DO 25 I=1,2 J=NODES(I) PXYZ=PXYZ+FUN(J)*PR(I) 25 CONTINUE A1=XJAC(MM,2) A2=-XJAC(MM,1) 30 DO 60 I=1,2

25、J=NODES(I) K2=2*J K1=K2-1 Q=PXYZ*FUN(J)*H(LX)*FACTNUS RF(K1)=RF(K1)+Q*A1 RF(K2)=RF(K2)+Q*A2 60 CONTINUE 70 CONTINUE RETURN ENDC * SUBROUTINE ASLOAD (R) DIMENSION R(*) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN3/ RF(8),SKE(8,8),NN(8) DO 20 I=1,8 L=NN(I) IF (L.EQ.0) GO TO 20 R(L)=R(L)+RF(I) 20 CONTINUE RET

26、URN ENDC * SUBROUTINE DECOP (SK,MA) DIMENSION SK(*),MA(*) COMMON /CMN2/ N,MX,NH DO 50 I=2,N L=I-MA(I)+MA(I-1)+1 K=I-1 L1=L+1 IF (L1.GT.K) GO TO 30 DO 20 J=L1,K IJ=MA(I)-I+J M=J-MA(J)+MA(J-1)+1 IF (L.GT.M) M=L MP=J-1 IF (M.GT.MP) GO TO 20 DO 10 LP=M,MP IP=MA(I)-I+LP JP=MA(J)-J+LP SK(IJ)=SK(IJ)-SK(IP)

27、*SK(JP) 10 CONTINUE 20 CONTINUE 30 IF (L.GT.K) GO TO 50 DO 40 LP=L,K IP=MA(I)-I+LP LPP=MA(LP) SK(IP)=SK(IP)/SK(LPP) II=MA(I) SK(II)=SK(II)-SK(IP)*SK(IP)*SK(LPP) 40 CONTINUE 50 CONTINUE RETURN ENDC * SUBROUTINE FOBA (SK,MA,R) DIMENSION SK(*),MA(*),R(*) COMMON /CMN2/ N,MX,NH DO 10 I=2,N L=I-MA(I)+MA(I

28、-1)+1 K=I-1 IF (L.GT.K) GO TO 10 DO 5 LP=L,K IP=MA(I)-I+LP R(I)=R(I)-SK(IP)*R(LP) 5 CONTINUE 10 CONTINUE DO 20 I=1,N II=MA(I)45 R(I)=R(I)/SK(II)20 CONTINUE DO 30 J1=2,N I=2+N-J1 L=I-MA(I)+MA(I-1)+1 K=I-1 IF (L.GT.K) GO TO 30 DO 25 J=L,K IJ=MA(I)-I+J55 R(J)=R(J)-SK(IJ)*R(I)25 CONTINUE30 CONTINUE RETU

29、RN ENDC * SUBROUTINE STRESS(COOR,MEL,JR,AE,R,STRE) DIMENSION XYZ(2,4),DNX(2,4),AE(4,*),STRE(3,*), & COOR(2,*),MEL(5,*),JR(2,*),RJAC(2,2),SIG(3), & B(3,8),R(*),iven(4) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN3/ RF(8),SKE(8,8),NN(8) COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2) DO 106 IE=1,NE NME=MEL(5,

30、IE) DO 300 K=1,4 IEK=MEL(K,IE) DO 310 M=1,2310 XYZ(M,K)=COOR(M,IEK) DO 320 M=1,2 JRR=2*(K-1)+M320 NN(JRR)=JR(M,IEK)300 CONTINUE E=AE(1,NME) U=AE(2,NME) D1=E*(1.00-U)/(1.00+U)*(1.00-2.00*U) D2=E*U/(1.00+U)*(1.00-2.00*U) D3=0.50*E/(1.00+U) SS=0.0 RR=0.0 CALL FDNX (XYZ,DNX,DET,RR,SS,RJAC,iven,IE) DO 30

31、 I=1,4 II=2*(I-1) J1=II+1 J2=II+2 BI=DNX(1,I) CI=DNX(2,I) B(1,J1)=BI B(2,J1)=0. B(3,J1)=CI B(1,J2)=0. B(2,J2)=CI B(3,J2)=BI 30 CONTINUE DO 55 II=1,3 SIG(II)=0.00 55 CONTINUE DO 70 K=1,8 NA=NN(K) IF (NA.EQ.0) GO TO 70 DO 60 L=1,3 SIG(L)=SIG(L)+B(L,K)*R(NA) 60 CONTINUE 70 CONTINUE SX=D1*SIG(1)+D2*SIG(

32、2) SY=D2*SIG(1)+D1*SIG(2) SXY=D3*SIG(3) STRE(1,IE)=SX STRE(2,IE)=SY STRE(3,IE)=SXY106 CONTINUE CALL OUTSTRE(NE,STRE) RETURN ENDC * SUBROUTINE FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE) DIMENSION XYZ(2,*),DNX(2,*),RJAC(2,2),iven(*) COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2) CALL FUN8 (XYZ,R,S,DET) IF (DET.LT.1

33、.0E-5)THEN WRITE(7,600) NEE,R,S,det WRITE(7,*) (iven(m),m=1,4) STOP ENDIF REC=1.00/DET RJAC(1,1)=REC*XJAC(2,2) RJAC(2,2)=REC*XJAC(1,1) RJAC(2,1)=-REC*XJAC(2,1) RJAC(1,2)=-REC*XJAC(1,2) DO 30 K=1,4 DO 20 I=1,2 DNX(I,K)=0. DO 25 M=1,2 DNX(I,K)=DNX(I,K)+RJAC(I,M)*PN(M,K) 25 CONTINUE 20 CONTINUE 30 CONT

34、INUE600 FORMAT (1X,'ERR0R* NEGTIVE OR ZERO ' * 'JACOBIAN DETERMINANT FOR ' * 'ELEMENT'/'ELE.=',I5,' R=',F10.5,6X,'S=',F10.5, * 'det=',f12.5) RETURN ENDC * SUBROUTINE FUN8 (XYZ,R,S,DET) DIMENSION XYZ(2,*),XI(4),ETA(4) COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2) DATA XI/-1.0,1.0,1.0,-1.0/ DATA ETA/-1.0,-1.0,1.0,1.0/ DO 10 I=1,4 G1=(1.0+XI(I)*R) G2=(1.0+ETA(I)*S) FUN(I)=0.25*G1*G2 PN(1,I)=0.25*XI(I)*G2 PN(2,I)=0.25*ETA(I)*G1 10 CONTINUE DO 80 I=1,2 DO 75 J=1,2 DET=0.00 DO 70 K=1,4 DET=DET+PN(I,K)*XYZ(