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

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/)')'输出计算结果:'

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

边坡稳定计算-简化毕肖甫法相关推荐

线性代数行列式计算之升阶法
线性代数行列式计算之升阶法声明与简介线性代数行列式计算之升阶法是利用行列式展开式的性质(行列式等于某一行或列乘其对应的代数余子式)在原有的行列式上增加1行或列1和0,增加之后方便消除其它行或列,子 ...
定积分的计算（换元法）习题
前置知识:定积分的计算(换元法) 习题1 已知 f ( x ) f(x) f(x),计算 ∫ a b f ′ ( 2 x ) d x \int_a^bf'(2x)dx ∫abf′(2x)dx 解:原 ...
MATLAB入门之rref计算简化矩阵行阶梯形式
rref 计算简化矩阵行阶梯形式例如: 求解线性方程组的解 %3x+7y+2z=3 %6x+8y+10z=5 %5x+6y+14z=7 >> A=[3,7,2;6,8,10;5,6,14 ...
cass高程点内插插件_聊聊CASS土方计算那些事-DTM法
DTM法,也可以叫三角网法.个人认为,如果数据质量过关,计算设计面是平面的土方工程,结果比较准确的一种方法.我常常用来此法来验算方量,和方格网法计算的结果对比.此方法计算过程几乎全部采用实测高程点.只 ...
##（C语言） CSP 201612-2 工资计算（打表法）（100分）
试题编号: 201612-2 试题名称: 工资计算时间限制: 1.0s 内存限制: 256.0MB 问题描述小明的公司每个月给小明发工资,而小明拿到的工资为交完个人所得税之后的工资.假设他一个月的 ...
计算流体力学有限体积法
有限体积法是计算流体力学中的一种数值模拟方法,它通过对流体的体积进行有限的划分,在每一个体素中分别求解流体的物理量,并通过体素间的相互作用得到整个流体的特性.有限体积法能够模拟复杂的流动状态,并且可以 ...
【Python】圆周率 Pi (π) 的计算（蒙特卡罗法+公式法）
引言圆周率(Pi)是圆的周长与直径的比值,一般用希腊字母 π 表示,是数学中最重要和最奇妙的数字之一.本文教你如何使用 Python 编程实现圆周率的简单计算. 计算蒙特卡罗法 1×1 的正方形里 ...
平均风向风速计算（单位矢量法）
单位矢量法计算公式: C# 代码 static void Main(string[] args){int[] Met_Dir_Test = { 259, 259, 259, 259, 259, 259 ...
（旋转体体积的计算）利用元素法简单解答空间几何体问题——高等数学
相信很多人初学的时候和我一样对这种三维空间的几何体计算方面有困难.我也曾百度过关于几何体体积/表面积的求法,但是始终不是很明白百度上的那种方法.这篇文章让你彻底理解这个万能的几何思想:"元素 ...

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

边坡稳定计算-简化毕肖甫法相关推荐

最新文章

热门文章