边坡稳定计算-简化毕肖甫法

c     文件名: BISHOP.FOR

C     边坡稳定计算(一)  (简化毕肖甫法)

      CHARACTER DFALE*12,YN*1

      COMMON /NB/NB,NS,N /C/C(10),T(10),GS(10)

     *      //X(41),YW(41),Y(10,41),Q(41),DX

      COMMON /TOL/TOL1,DS/EQ/EQH/QA/QA(23)

      COMMON /XU/XU(24)/YU/YU(10,24)/WL/WL(24)

      COMMON /XA/XA,XB,YA,YB

      WRITE(*,'(2X,A\)')

     $  '已有原始数据文件(Y--已有  N--没有)吗?[N]:'

      READ(*,'(A1)')YN

      IF(YN.EQ.'y')YN='Y'

      IF(YN.NE.'Y')YN='N'

      IU=0

      IF(YN.EQ.'Y')THEN

        DFALE='BISHOP.TXT'

        WRITE(*,'(2X,3A\)')'输入已有数据文件名[',DFALE,']:'

        READ(*,'(A12)')DFALE

        IF(DFALE.EQ.' ')DFALE='BISHOP.TXT'

        IU=7

        OPEN(IU,FILE=DFALE,STATUS='OLD')

      END IF

      WRITE(*,'(/2X,A)')'输入数据部份:'

      IF(YN.NE.'Y')WRITE(*,'(2X,A\)')

     $  '输入: 土坡表面几何特征点数NB,土层数NS,土条数N='

      READ(IU,*)NB,NS,N

      WRITE(*,'(/3A14/I10,I16,I14)')

     $ '地表特征点数','土层数','土条数',NB,NS,N

      IF(YN.NE.'Y')WRITE(*,'(2X,A\)')

     $'输入: 滑弧最小深度DS,水平地震系数QEH,相对精度TOL1='

      READ(IU,*)DS,EQH,TOL1

      WRITE(*,'(/A10,2A12/F10.2,2F12.2)')

     $ '滑弧深度','地震系数','相对精度',DS,EQH,TOL1

      IF(YN.NE.'Y')WRITE(*,'(2X,A\)')

     $'输入: 圆心可能的左侧坐标XA,右侧坐标XB,垂直坐标下限YA,上限YB='

      READ(IU,*)XA,XB,YA,YB

      WRITE(*,'(/2X,A/A10,3A12/F10.2,3F12.2)')'圆心初始坐标:',

     $ 'XA  ','XB  ','YA  ','YB  ',XA,XB,YA,YB

      IF(YN.NE.'Y')WRITE(*,'(2X,A/2X,A\)')'输入地表特征点水平座标:',

     $  'X1,X2,...Xnb='

      READ(IU,*)(XU(J),J=1,NB)

      IF(YN.NE.'Y')THEN

        WRITE(*,'(2X,A)')'输入地表特征点下各土层界面的垂直座标:'

        DO 40 I=1,NS

          WRITE(*,'(2X,A,I2,A\)')'输入第',I,' 层: Y1,Y2...Ynb='

          READ(IU,*)(YU(I,J),J=1,NB)

40      CONTINUE

      ELSE

        READ(IU,*)((YU(I,J),J=1,NB),I=1,NS)

      END IF

      WRITE(*,'(/2X,A/A10,A10,10(A7,I1,A4))')'输出地表特征点座标:',

     $ '特征点号','X(米)',('Y',J,'(米)',J=1,NS)

      DO 50 I=1,NB

        WRITE(*,'(A,I2,A,F11.2,10F12.2)')'   (',I,' ) ',

     $  XU(I),(YU(J,I),J=1,NS)

50    CONTINUE

      IF(YN.NE.'Y')WRITE(*,'(2X,A/2X,A\)')

     $'输入地表特征点下地下水位或浸润线的垂直座标:','WL1,WL2...WLnb='

      READ(IU,*)(WL(J),J=1,NB)

      IF(YN.NE.'Y')WRITE(*,'(2X,A/2X,A\)')

     $  '输入地表特征点间的荷载:','QA1,QA2...QAnb='

      READ(IU,*)(QA(J),J=1,NB)

      WRITE(*,'(/2X,A/A10,2A16)')'输出地表特征点间的地下水位与荷载:',

     $ '特征点号','地下水位X(米)','荷载(KN/m**2)'

      WRITE(*,'(5X,A,I2,A,F16.2,F16.2)')

     $  ('(',I,' )',WL(I),QA(I),I=1,NB)

      IF(YN.NE.'Y')THEN

        WRITE(*,'(2X,A)')'输入各土层的粘聚力,容重,内摩擦角:'

        DO 70 I=1,NS

          WRITE(*,'(2X,A,I2,A\)')'输入第',I,' 层: Ci,GSi,Ti='

          READ(IU,*)C(I),GS(I),T(I)

70      CONTINUE

      ELSE

        READ(IU,*)(C(I),GS(I),T(I),I=1,NS)

      END IF

      WRITE(*,'(/2X,A)')'输出各土层的粘聚力,容重,内摩擦角:'

      WRITE(*,'(/A8,3A12)')'土层号','粘聚力','容重','内摩擦角'

      WRITE(*,'(3X,A,I2,A,F10.2,2F12.2)')

     $  ('(',I,' )',C(I),GS(I),T(I),I=1,NS)

      IF(IU.GT.1)CLOSE(IU)

C     ------------------------------

      WRITE(*,'(//2X,A/)')'输出计算结果:'

C

      DX=(XU(NB)-XU(1))/FLOAT(N)

      X(1)=XU(1)

      Y(1,1)=YU(1,1)

      IF(NS.EQ.1) GOTO 26

      DO 20 K=2,NS

      Y(K,1)=YU(K,1)

      IF(Y(K,1).LT.Y(K-1,1)) Y(K,1)=Y(K-1,1)

20    CONTINUE

26    YW(1)=WL(1)

      J=2

      N1=N+1

      DO 21 I=2,N1

      X(I)=XU(1)+DX*FLOAT(I-1)

23    IF(X(I).LE.XU(J)) GOTO 22

      J=J+1

      GOTO 23

22    A=(X(I)-XU(J-1))/(XU(J)-XU(J-1))

      Y(1,I)=YU(1,J-1)+A*(YU(1,J)-YU(1,J-1))

      IF(NS.EQ.1) GOTO 25

      DO 24 K=2,NS

      Y(K,I)=YU(K,J-1)+A*(YU(K,J)-YU(K,J-1))

24    CONTINUE

25    YW(I)=WL(J-1)+(WL(J)-WL(J-1))*A

      IF(YW(I).LT.Y(1,I)) YW(I)=Y(1,I)

21    Q(I-1)=QA(J-1)

      WRITE(*,'(//20X,A)')'COODINATERS OF CIRCLE CENTRE'

      WRITE(*,'(/10X,60A)')('+',K=1,60)

      CALL CENTER

      STOP' '

      END

C     ----------------------------------------

      SUBROUTINE CENTER

      DIMENSION XX(5),YY(5),RA(5),FA(5)

      COMMON /NB/NB,NS,N/C/C(10),T(10),GS(10)//

     *         X(41),YW(41),Y(10,41),Q(41),DX

      COMMON /XA/XA,XB,YA,YB/TOL/TOL1,DS

      WRITE(*,'(6X,4A16)')'XC','YC','R ','Fs'

      DXA=(XB-XA)/4.0

      DYA=(YB-YA)/4.0

      XX(1)=(XA+XB)/2.0

      YY(1)=(YA+YB)/2.0

      CALL RADIUS (XX(1),YY(1),RA(1),FA(1))

      WRITE(*,'(A,4F16.2)')' FA(1)',XX(1),YY(1),RA(1),FA(1)

10    XX(2)=XX(1)-DXA

      XX(3)=XX(2)

      XX(4)=XX(1)+DXA

      XX(5)=XX(4)

      YY(2)=YY(1)-DYA

      YY(3)=YY(1)+DYA

      YY(4)=YY(2)

      YY(5)=YY(3)

      J=1

      FM=FA(1)

      FAMIN=FA(1)

      DO 14 I=2,5

      CALL RADIUS(XX(I),YY(I),RA(I),FA(I))

      WRITE(*,'(A,I1,A,4F16.2)')' FA(',I,')',XX(I),YY(I),RA(I),FA(I)

      FM=FM+FA(I)

      IF(FAMIN.LT.FA(I)) GOTO 14

      J=I

      FAMIN=FA(J)

14    CONTINUE

      FM=FM/5.0

      AA=ABS((FAMIN-FM)/FAMIN)

      IF(AA.LE.TOL1) GO TO 16

      FA(1)=FA(J)

      XX(1)=XX(J)

      YY(1)=YY(J)

      DXA=DXA/2.0

      DYA=DYA/2.0

      GOTO 10

16    WRITE (*,'(/10X,60A)')('+',I=1,60)

      WRITE(*,'(/10X,A,F7.3/46X,A,F7.3/47X,A,F7.2/30X,A,F7.4)')

     $ 'CENTERS COODINATE OF SLIDING RADIUS XC=',XX(J),'YC=',YY(J),

     $ 'R=',RA(J),'FACTOR OF SAFTY Fs=',FA(J)

      WRITE (*,'(/10X,60A)')('+',I=1,60)

      RETURN

      END

C     ------------------------------

      SUBROUTINE RADIUS(XC,YC,RC,FC)

      DIMENSION RB(3),FB(3)

      COMMON/XU/XU(24)/YU/YU(10,24)/NB/NB,NS,N

        COMMON /C/C(10),T(10),GS(10)//X(41),YW(41),

     *        Y(10,41),Q(41),DX/XA/XA,XB,YA,YB

      COMMON /TOL/TOL1,DS

      R1=SQRT((XC-XU(1))**2+(YC-YU(1,1))**2)

      R2=SQRT((XC-XU(NB))**2+(YC-YU(1,NB))**2)

      RMAX=AMIN1(R1,R2)

      RMIN=RMAX

      DO 10 I=2,NB

      R1=SQRT((XC-XU(I))**2+(YC-YU(1,I))**2)

      IF(R1.LT.RMIN) RMIN=R1

10    CONTINUE

      DR=(RMAX-RMIN)/4.0

      RB(1)=(RMAX+RMIN)/2.0

      CALL FACTOR(RB(1),XC,YC,FB(1))

12    RB(2)=RB(1)-DR

      RB(3)=RB(1)+DR

      CALL FACTOR (RB(2),XC,YC,FB(2))

      CALL FACTOR (RB(3),XC,YC,FB(3))

      FM=(FB(1)+FB(2)+FB(3))/3.0

      J=1

      FMIN=FB(1)

      DO 14 I=2,3

      IF(FMIN.LT.FB(I)) GOTO 14

      J=I

      FMIN=FB(J)

14    CONTINUE

      IF(ABS((FMIN-FM)/FMIN).LE.TOL1) GO TO 16

      RB(1)=RB(J)

      FB(1)=FMIN

      DR=DR/2.0

      GOTO 12

16    FC=FMIN

      RC=RB(J)

      RETURN

      END

C     -----------------------------------

      SUBROUTINE FACTOR(R,XC,YC,F)

      DIMENSION DL(60),CA(60),SA(60),YR(61)

      COMMON /NB/NB,NS,N/C/C(10),T(10),GS(10)

      COMMON /EQ/EQH//X(41),YW(41),Y(10,41),

     *    Q(41),DX

      COMMON /XA/XA,XB,YA,YB/TOL/TOL1,DS

      N1=N+1

      DO 10 I=1,N1

      IF(ABS(XC-X(I)).GT.R) GO TO 11

      YR(I)=SQRT(R**2-(X(I)-XC)**2)+YC

      IF(YR(I).GE.Y(1,I)) GO TO 10

11    YR(I)=Y(1,I)

10    CONTINUE

      DO 15 I=1,N

      DL(I)=SQRT((X(I+1)-X(I))**2+(YR(I)-YR(I+1))**2)

      CA(I)=(X(I+1)-X(I))/DL(I)

      SA(I)=(YR(I)-YR(I+1))/DL(I)

15    CONTINUE

      F=1.0

16    WS=0.0

      Z=0.0

      DO 20 I=1,N

      YRM=(YR(I)+YR(I+1))/2.0

      YM=(Y(1,I)+Y(1,I+1))/2.0

      YT=YM

      YEQ=(YT+YRM)/2.

      YWM=(YW(I)+YW(I+1))/2.0

      IF (YT.GE.YWM) U=YRM-YT

      U=YRM-YWM

      IF(U.LT.0.0) U=0.0

C     IF(YRM.LE.YM) GO TO 20

      W=DX*Q(I)

      IF(NS.EQ.1) GO TO 25

      NS1=NS-1

      DO 22 J=1,NS1

      YM=(Y(J+1,I+1)+Y(J+1,I))/2.0

      IF(YRM.GT.YM) GO TO 22

      JA=J

      GO TO 24

22    W=W+GS(J)*DX*(YM-(Y(J,I+1)+Y(J,I))/2.0)

25    JA=NS

24    W=W+GS(JA)*DX*(YRM-(Y(JA,I+1)+Y(JA,I))/2.0)

      WQH=.25*W*EQH

      WQV=WQH/3.

      WS=WS+(W+WQV)*SA(I)+WQH*YEQ/YRM

      SIT=SIN(3.141592*T(JA)/180.0)

      COT=COS(3.141592*T(JA)/180.0)

      ZU=C(JA)*DL(I)*CA(I)+(W+WQV-U*DL(I)*CA(I))*SIT/COT

      ZL=CA(I)+SA(I)*SIT/(COT*F)

      Z=Z+ZU/ZL

20    CONTINUE

      FA=Z/WS

250   FORMAT(1X,'FC=',E15.7)

      IF(ABS(FA-F)/FA.LE.TOL1) GO TO 30

      F=FA

      GOTO 16

30    F=FA

      IF(F.LE.0.0) F=100000.0

      RETURN

      END

 

 BISHOP.TXT

8 2 20
0.00 0.00 0.02
10.00 60.00 50.00 120.00
20.00 50.00 54.00 54.50 57.00 59.50 70.00 80.00
160.30 160.30 155.00 155.00 150.00 150.00 150.00 150.00
160.30 160.30 156.00 156.00 156.00 156.00 156.00 156.00
162.00 162.00 162.00 162.00 162.00 162.00 162.00 162.00
0.00 0.00 0.00 0.00 0.00 0.50 0.00 0.00
2.80 1.99 18.00
2.50 2.00 22.00