《土质边坡稳定分析》教学资料（原理、方法、程序）第十一章 程序设计

340土质边坡德定分析一原理·方法程序 READ(5,*)LSL外边坡线总数 READ(S,*)LNUM()l=1,LSL)外边坡线编号 READ(5,NN!点总数 READ,*)L,XN(I,YN(D)1=1,NN)点坐标 703 FORMAT(T25,***) WRITE(6, 704) ORMAT(IOX, THEABSCISSAVALUESOFTHEUPPERANDLOWER') CALL SEARCLNUM,C,DS,CXCY,XN,YN,LSL, NLOW NUPP)!寻找上、下交点所在线段 CALL DIVI(NIC,IWRS,DS,CXCY,XN,YN!确定上、下交点坐标,条分 S, YTENSION, NLOW, NUPP X, Y, ALI 711 FORMAT(TS, NO, T17, X,T32,Y, T47,ALF) DO309I=1 WRITE(6,710)L, X),Y(),ALFO 710 FORMAT(1X,15,3F15.6) 309 CONTINUE 12 CALLCLOSEFL END SUBROUTINE CLOSEFL CLOSE(5) CLOSE(6) RETURN SUBROUTINE DIVI(N, IC, TWRS, DS, CX, CY,XN,YN S, YTENSION, NLOW,NUPP,X,Y,ALF) INTEGERNS, N, K, IC(80,3), KI, KTEMP, TWRS,J,J1,JTEMPNLOW, NUPP REAL RO, DS, CY,XN(80), YN(80),SK, T, XL, CX, XUl, YTENSION REAL YUI,SJ, XU, ROO, YU, X2(80), Y2(80), X(80), Y(80), ALF(80) S(E1,F1E2,F2)=(F2-F1)(E2-E1) NS=N-1 RO=DS-CY KI=IC(NLOW, 2) IF(XN(K).GTXN(KDGOTO190 KTEMP=K KI=KTEMP 190 IF(ABS(XN(K)-XN(KI).LT0.00005 ) GOTO23 SK=S(XN(K),YN(K), XN(KD,YN(KD) T=SK*XN(K)-YN(K)+CY XL=(CX+SK*T+SQRT(1+SK*2)RO*2-(SK*CX-T)*2))(1+SK**2) !按式(16)计算下交点的横坐标 GOTO28 23 XLFXNK 28 CONTINUE IF(WRS NE.O)THEN XUI=CX-SQRT(RORO-(YTENSION-CY*2) YUI=CY+SQRT(RO**2-ABS(XUl-CX)**2) J=IC(NUPP, D) IF(XNO).GT.XNJDGOTO191 JTEMP- JI=JTEMP

340 土质边坡稳定分析 原理 ⋅ 方法 ⋅ 程序 READ(5,*)LSL!外边坡线总数 READ(5,*)(LNUM(I),I=1,LSL)!外边坡线编号 READ(5,*)NN!点总数 READ(5,*)(L,XN(I),YN(I),I=1,NN)!点坐标 WRITE(6,703) 703 FORMAT(T25,'******************'//) WRITE(6,704) 704 FORMAT(10X,'THEABSCISSAVALUESOFTHEUPPERANDLOWER') CALL SEAR(LNUM,IC,DS,CX,CY,XN,YN,LSL,NLOW,NUPP)!寻找上 下交点所在线段 CALL DIVI(N,IC,IWR5,DS,CX,CY,XN,YN 确定上 下交点坐标 条分 $,YTENSION,NLOW,NUPP,X,Y,ALF) WRITE(6,711) 711 FORMAT(T5,'NO.',T17,'X',T32,'Y',T47,'ALF') DO 309 I=1,N WRITE(6,710)I,X(I),Y(I),ALF(I) 710 FORMAT(1X,I5,3F15.6) 309 CONTINUE 12 CALLCLOSEFL END SUBROUTINE CLOSEFL CLOSE(5) CLOSE(6) RETURN END SUBROUTINE DIVI(N,IC,IWR5,DS,CX,CY,XN,YN $,YTENSION, NLOW,NUPP,X,Y,ALF) INTEGERN S,N,K,IC(80,3),K1,KTEMP,IWR5,J,J1,JTEMP,NLOW,NUPP REAL RO,DS,CY,XN(80),YN(80),SK,T,XL,CX,XU1,YTENSION REAL YU1,SJ,XU,ROO,YU,X2(80),Y2(80),X(80),Y(80),ALF(80) S(E1,F1,E2,F2)=(F2-F1)/(E2-E1) NS=N-1 RO=DS-CY K=IC(NLOW,1) K1=IC(NLOW,2) IF(XN(K).GT.XN(K1))GOTO190 KTEMP=K K=K1 K1=KTEMP 190 IF(ABS(XN(K)-XN(K1)).LT.0.00005)GOTO23 SK=S(XN(K),YN(K),XN(K1),YN(K1)) T=SK*XN(K)-YN(K)+CY XL=(CX+SK*T+SQRT((1+SK**2)*RO**2-(SK*CX-T)**2))/(1+SK**2) !按式(11.6)计算下交点的横坐标 GOTO28 23 XL=XN(K) 28 CONTINUE IF(IWR5.NE.0)THEN XU1=CX-SQRT(RO*RO-(YTENSION-CY)**2) YU1=CY+SQRT(RO**2-ABS(XU1-CX)**2) ENDIF J=IC(NUPP,1) J1=IC(NUPP,2) IF(XN(J).GT.XN(J1))GOTO191 JTEMP=J J=J1 J1=JTEMP

第11章程序设计341 191 F(ABS(XNO)-XNOI)LT0.00005)GOTO82 J=S(XN),YN),XNI),YNJD) T=SJ"XNOFYNOHCY XU=Cx+SJ*T-SQRI(1+SJ**2)RO**2-(SJCX-)*2))(1+SJ**2) XU上交点的横坐标 ROO=RO**2-ABS(XU-CX)**2 IF( ROO.LT0 RETURN无交点,需修改RO,重新计算 YU=CY+SQRT(RO*#2-ABS(XU-CX) 2) IF( CYGTYN( J)AND CYGT.YN(J) RETURN凸弧否定,重算 IF(WRSNE.0 AND YU. LTYUI)XU=XUI XU=XNO) DO541=1,NS划分土条求各点坐标 X2(=(-)+(XLXU(NS-1)+XU!土条宽度B=(XL-XU(NS-1) Y2(FCY+SQRT(RO*2-ABS(X2()-CX)**2) 41 CONTINUE IF(ABS(Y2(NS-Y2(1)). LT0.OlTHEN WRITE(6, THE SLIP SURFACE INTERCEPTS A HORIZONTAL SURFACE RETURN ENDIF N=NS+1 Y2(NFY2(NS 计算土条中点坐标 DO544=2.NS X(=(X2(1)+X2(1-1)2 Y(=Y2(+Y2(-1)y ALF(I=ATAN(Y2(1)Y2(-1))(X2(1}X2(-1)) 4 CONTINUE ALF(IFALF(2) ALF(NALF(NS) X(NX2(N) X(1)=x2(1) Y(1)=Y2(1) RETURN END SUBROUTINE SEAR(LNUM,IC, DS, CX, CY,XN,YN, LSL, NLOW, NUPP) REAL R, DS, CY, R2, CX, XN(80),YN(80),Rl, D, XXL, SLO INTEGER I, I1, K, LNUM(80),J1, J2, IC(80,3), IN, NLOW,NUPP LSL I=1NLOW=下交点所在直线段的线段号 I1=0NUPP=上交点所在直线段的线段号 12K=LNUMD JI=IC(K, 1) J2=IC(K, 2) R= DS-CY!R=圆弧半 R2=SQRT(ABS(CX-XN(2)**2+ABS(CY-YNO2)* *2) RI=SQRT(ABS(CX-XNO1)**2+ABS(CY-YNJD))**2) R,R2为圆心至该线段两端点的距离 IF(REQ R2AND II.NEO)GOTOll IF(RGTR1 ANDR GTR2)GOTO11!无交点 IF(RLTR1 ANDRLTR2GOTO19有两个或零个交点 GOT020 19 CALL DD(XNJI),YND), XN(2), YN(J2), CX, CY, D, IN) XXL=XN(2)-XNO1

第 11 章 程序设计 341 191 F(ABS(XN(J)-XN(J1)).LT.0.00005)GOTO82 SJ=S(XN(J),YN(J),XN(J1),YN(J1)) T=SJ*XN(J)-YN(J)+CY XU=(CX+SJ*T-SQRT((1+SJ**2)*RO**2-(SJ*CX-T)**2))/(1+SJ**2) !XU 上交点的横坐标 ROO=RO**2-ABS(XU-CX)**2 IF(ROO.LT.0)RETURN!无交点,需修改 RO,重新计算 YU=CY+SQRT(RO**2-ABS(XU-CX)**2) IF(CY.GT.YN(J).AND.CY.GT.YN(J1))RETURN!凸弧,否定,重算 IF(IWR5.NE.0.AND.YU.LT.YU1)XU=XU1 GOTO20 82 XU=XN(J) 20 X2(1)=XU DO541I=1,NS!划分土条,求各点坐标 X2(I)=(I-1)*(XL-XU)/(NS-1)+XU!土条宽度 B=(XL-XU)/(NS-1) Y2(I)=CY+SQRT(RO**2-ABS(X2(I)-CX)**2) 541 CONTINUE IF(ABS(Y2(NS)-Y2(1)).LT.0.01)THEN WRITE(6,*)'THE SLIP SURFACE INTERCEPTS A HORIZONTAL SURFACE' RETURN ENDIF N=NS+1 X2(N)=X2(NS) Y2(N)=Y2(NS) C 计算土条中点坐标 DO544I=2,NS X(I)=(X2(I)+X2(I-1))/2. Y(I)=(Y2(I)+Y2(I-1))/2. ALF(I)=ATAN((Y2(I)-Y2(I-1))/(X2(I)-X2(I-1))) 544 CONTINUE ALF(1)=ALF(2) ALF(N)=ALF(NS) X(N)=X2(N) Y(N)=Y2(N) X(1)=X2(1) Y(1)=Y2(1) RETURN END SUBROUTINE SEAR(LNUM,IC,DS,CX,CY,XN,YN,LSL,NLOW,NUPP) REAL R,DS,CY,R2,CX,XN(80),YN(80),R1,D,XXL,SLO INTEGER I,I1,K,LNUM(80),J1,J2,IC(80,3),IN,NLOW,NUPP,LSL I=1 !NLOW=下交点所在直线段的线段号 I1=0 !NUPP=上交点所在直线段的线段号 12K=LNUM(I) J1=IC(K,1) J2=IC(K,2) R=DS-CY!R=圆弧半径 R2=SQRT(ABS(CX-XN(J2))**2+ABS(CY-YN(J2))**2) R1=SQRT(ABS(CX-XN(J1))**2+ABS(CY-YN(J1))**2) !R1,R2 为圆心至该线段两端点的距离 IF(R.EQ.R2.AND.I1.NE.0)GOTO11 IF(R.GT.R1.AND.R.GT.R2)GOTO11 !无交点 IF(R.LT.R1.AND.R.LT.R2)GOTO19 !有两个或零个交点 GOTO20 19 CALL DD(XN(J1),YN(J1),XN(J2),YN(J2),CX,CY,D,IN) XXL=XN(J2)-XN(J1)

342土质边坡德定分析一原理·方法程序 IF(ABS(XXL).LT0.0001THEN SLO=L0E-6 ELSE SLO=(YN(J2YN(J1)XXL!斜率 ENDIF IF(NEQ 1 OR SLO. LEO)GOTOll C第一条外边坡线为水平同时圆弧与其有两个交点不接受继续搜索 IF( R. LT D ORREQ. D)GOTO1无交点 NLOW=LNUM NUPPENLOW RETURN 0 IF(I NE.OGOTO14 NLOW=LNUM 11=11+I 11 CONTINUE IF(EQ 1 AND R.LTRIAND R GT R2)RETURN l=+1 IF(GT.LSLGOTO14 GOTO12 14 NUPP=LNUM 15 CONTINUE RETURN SUBROUTINE DD(X1, Y1, X2, Y2, XC,YC, D, IN) REAL XI, Y1, X2, Y2, XC,YC, D REAL XXL, AK, AKF X,YXT,YT NTEGER IN XXLEX2-XI IF(ABS(XXL). LT0.001)XXL=1.0E-6 IF(ABS((Y2-Y1)XXL). LT0.OOIRETURN AK=(Y2-YI/XXL AKF=l/AK X=(AK*Xl-Y1+AKF*XC+YC)(AK+AKF) Y=AK*XAK*X1+Y1!从圆心作垂直与线段的线获与线段的交点 X1OX-X2 YT=(YY1)*(Y-Y2) IF(XT.LT0OR YT. LT.OGOTO12 IN=11交点不在线段内 RETURN 12D=SQRT(X-XC)**2+(YYO)*2)!线段与圆心的垂直距离 RETURN END 3.例题 数据文件 CIRCLE DAT 00. 24.235-35.164-42.046 21 2 3

342 土质边坡稳定分析 原理 ⋅ 方法 ⋅ 程序 IF(ABS(XXL).LT.0.0001)THEN SLO=1.0E-6 ELSE SLO=(YN(J2)-YN(J1))/XXL !斜率 ENDIF IF(IN.EQ.1.OR.SLO.LE.0)GOTO11 C 第一条外边坡线为水平,同时圆弧与其有两个交点,不接受,继续搜索 IF(R.LT.D.OR.R.EQ.D)GOTO11 !无交点 NLOW=LNUM(I) NUPP=NLOW RETURN 20 IF(I1.NE.0)GOTO14 NLOW=LNUM(I) I1=I1+1 11 CONTINUE IF(I.EQ.1.AND.R.LT.R1.AND.R.GT.R2)RETURN I=I+1 IF(I.GT.LSL)GOTO14 GOTO12 14 NUPP=LNUM(I) 15 CONTINUE RETURN END SUBROUTINE DD(X1,Y1,X2,Y2,XC,YC,D,IN) REAL X1,Y1,X2,Y2,XC,YC,D REAL XXL,AK,AKF,X,Y,XT,YT INTEGER IN D=9999. IN=0 XXL=X2-X1 IF(ABS(XXL).LT.0.001)XXL=1.0E-6 IF(ABS((Y2-Y1)/XXL).LT.0.001)RETURN AK=(Y2-Y1)/XXL AKF=1./AK X=(AK*X1-Y1+AKF*XC+YC)/(AK+AKF) Y=AK*X-AK*X1+Y1 !从圆心作垂直与线段的线,获与线段的交点 XT=(X-X1)*(X-X2) YT=(Y-Y1)*(Y-Y2) IF(XT.LT.0.OR.YT.LT.0)GOTO12 IN=1 !交点不在线段内 RETURN 12 D=SQRT((X-XC)**2+(Y-YC)**2) !线段与圆心的垂直距离 RETURN END 3. 例题 数据文件 CIRCLE.DAT 0 0. -24.235 -35.164 -42.046 21 9 1 2 3 2 3 1 3 4 1 2 5 3 5 6 2 6 7 2

第11章程序设计343 000-25000 70.000-35.000 -40.000-27.000 54.000-31.000 70.000-31.000 52.000-24.000 10-70.000-24.000 计算成果 CIRCLEC: THE ABSCISSA VALUES OF THE UPPER AND LOWEI NO 51.522050 1.094647 50.955680 3-49822940-31.978910 -48.690200 30491210 863838 -47.557460 29278430 770531 46.424730 28265840 685123 7-45.29199027411070 605365 41.893780 25.566700 387709 11-40.761050 -25.147810 319849 1239628310-24.813490 13-38495570 24.558860 188260 14-37.362830-24.380470 123833 3509736024244130 003743 33.964630-24.284500 067425 1832831890 -131382 4.584810 -.19588 30.566420 4848610 261233 30.000050 5.000020 261233 11.3.3构造任意形状滑裂面的程序 SPLINE FOR 1.说明 本程序(刘德贵,1980)用来计算已知第二种边界条件(曲线两端点处的二阶导数值) 三次样条函数的一阶导数值,并利用样条函数对一元函数进行插值,进而构筑包含直线段的 复合滑裂面。对于滑裂面上控制点区间为直线段时,则直接线性插值,而对于非直线段区间 则运用上述介绍的方法进行三次样条插值。滑裂面被NS1个点分为NS1-1段。从上交点向 下交点编号为1,2,…NS1-1。程序用LS代表此NS1-1段中为直线段的线段总数。在数组 LSA(I)(=1,2,…LS)中存入这些直线段的编号 LS,滑裂面上NS1-1条线段中为直线段的总数 SA(LS),滑裂面上线段为直线段的编号;

第 11 章 程序设计 343 7 8 2 5 9 3 9 10 3 3 1 2 3 10 1 .000 -25.000 2 -30.000 -25.000 3 -50.000 -35.000 4 -70.000 -35.000 5 -40.000 -27.000 6 -50.000 -29.000 7 -54.000 -31.000 8 -70.000 -31.000 9 -52.000 -24.000 10 -70.000 -24.000 计算成果 CIRCLE.C THE ABSCISSA VALUES OF THE UPPER AND LOWER NO. X Y ALF 1 -51.522050 -35.000000 1.094647 2 -50.955680 -33.901810 1.094647 3 -49.822940 -31.978910 .969015 4 -48.690200 -30.491210 .863838 5 -47.557460 -29.278430 .770531 6 -46.424730 -28.265840 .685123 7 -45.291990 -27.411070 .605365 8 -44.159260 -26.687410 .529825 9 -43.026520 -26.076840 .457513 10 -41.893780 -25.566700 .387709 11 -40.761050 -25.147810 .319849 12 -39.628310 -24.813490 .253490 13 -38.495570 -24.558860 .188260 14 -37.362830 -24.380470 .123833 15 -36.230100 -24.275990 .059923 16 -35.097360 -24.244130 -.003743 17 -33.964630 -24.284500 -.067425 18 -32.831890 -24.397580 -.131382 19 -31.699150 -24.584810 -.195887 20 -30.566420 -24.848610 -.261233 21 -30.000050 -25.000020 -.261233 11. 3. 3 构造任意形状滑裂面的程序SPLINE.FOR 1. 说明 本程序 刘德贵 1980 用来计算已知第二种边界条件(曲线两端点处的二阶导数值) 三次样条函数的一阶导数值 并利用样条函数对一元函数进行插值 进而构筑包含直线段的 复合滑裂面 对于滑裂面上控制点区间为直线段时 则直接线性插值 而对于非直线段区间 则运用上述介绍的方法进行三次样条插值 滑裂面被 NS1 个点分为 NS1−1 段 从上交点向 下交点编号为 1, 2, , NS1−1 程序用 LS 代表此 NS1−1 段中为直线段的线段总数 在数组 LSA(I)(I=1, 2, , LS)中存入这些直线段的编号 LS, 滑裂面上 NS1−1 条线段中为直线段的总数 LSA(LS) 滑裂面上线段为直线段的编号

344土质边坡稳定分析一原理·方法·程序 AI(14),滑裂面上控制点处四个与样条插值有关的数,K,K2,K3,K4; A(),滑裂面上控制点处的一阶导数值m B(D),一维数组,工作单元 2.源程序 SUBROUTINE SPI(LS, LSA, X2, Y2) DIMENSION KO2(50),X1(50,Y1(50),AI(50,4)LSA(60),A(50),B(50) DIMENSION X2(80),Y2(80) COMMON /A23/NS1, KQ2, X1, Y1/A13/DD, IWRl, IWR2, IWR3 COMMON/WALL/TWALL, G WALL, HMW,EWALL, ETA INTEGER4 OPTION(6),OP1(6) COMMON/OPP/OPTION. OPI OPEN(S, FILE STATUS=UNKNOWN OPEN(6, FILE STATUS=UNKNOWN READ(S NSI DO=LNSI READO, )KQ2), Xl,YI) ENDDO 如果滑裂面上控制点总数NS1=2,则按直线段处理并计算滑裂面与土条侧面边界线交点的xy值 LSAFO ENDDO IF(LS.NEO)READ(, (LSA I=l, LS) sA(S+IFNS C 滑裂面上第NS个控制点对应的土条侧面边界线的编号 Xl(①),Yl(),滑裂面上第I个控制点的x,y值 N1=0 DO 100IL=LLS+ LI=NI+1 N2=NI-1 IF(LSEQ0GOTO83!滑裂面上没有直线段 IF(NIGE NS1)GOTO 40 AI(Nl,1=0.!计算直线段左端点处K1,K2,K,K4 AI(Nl,3)=Y1(N1+1}Y1(N1) AI(N1, 4FYI(ND) 40IF(N1L1EQ0GOTO100!滑裂面上没有直线段 83B(Ll=1.5°(Y1(L+1Y1(L1)/(X1(L1+1)-Xl(L1) A(L1)=0.5!插值区间左端点处的一阶导数值 CN=3.0*(Y1(Nl)Y1(N2)(X1(Nl}X1(N2) IF(N1-L1EQ.1)GOTO55!插值区间只具有左右两个端点 L2=L1+1 DO 50 J=L2 N2 HI=XIO)(-1) H=XIO+1)-XI) AF=H1/(H+H1)!计算α值 BT=3.*(1.-AF)°Y1(J-Y1(J-1)MHI+AF(Y1(J+1}Y1)/H) 计算β值 (JAF(2.+(1.-AF)*A(J-)!计算插值区间内插值节点处的一阶导数值m 50B(J=(BT-(1.-AF)*B(J-1)/(2+(1.-AF)A(J-1) 5A(NI=(CN-B(N2)/(2.+A(N2) A(N2=A(N2)*A(N2+1)+B(N2)!计算插值区间右端点处的一阶导数值mn

344 土质边坡稳定分析 原理 ⋅ 方法 ⋅ 程序 AI(I,4) 滑裂面上控制点处四个与样条插值有关的数 κ1, κ2, κ3, κ4 A(I) 滑裂面上控制点处的一阶导数值m j B(I) 一维数组 工作单元 2. 源程序 C SUBROUTINE SP1(LS,LSA,X2,Y2) DIMENSION KQ2(50),X1(50),Y1(50),AI(50,4),LSA(60),A(50),B(50) DIMENSION X2(80),Y2(80) COMMON /A23/NS1,KQ2,X1,Y1/A13/IDD,IWR1,IWR2,IWR3 COMMON/WALL/IWALL,GWALL,HMW,EWALL,ETA INTEGER*4 OPTION(6),OP1(6) COMMON/OPP/OPTION,OP1 OPEN(5,FILE=' ',STATUS='UNKNOWN') OPEN(6,FILE=' ',STATUS='UNKNOWN') READ(5,*)NS1 DO I=1,NS1 READ(5,*)KQ2(I),X1(I),Y1(I) ENDDO C 如果滑裂面上控制点总数 NS1=2,则按直线段处理并计算滑裂面与土条侧面边界线交点的 x,y 值; READ(5,*)LS DO I=1,NS1-1 LSA(I)=0 ENDDO IF(LS.NE.0)READ(5,*)(LSA(I),I=1,LS) LSA(LS+1)=NS1 C 滑裂面上第 NS1 个控制点对应的土条侧面边界线的编号 C X1(I) Y1(I) 滑裂面上第 I 个控制点的 x y 值 N1=0 DO 100 IL=1,LS+1 L1=N1+1 N1=LSA(IL) N2=N1-1 IF(LS.EQ.0)GOTO 83 ! 滑裂面上没有直线段 IF(N1.GE.NS1) GOTO 40 AI(N1,1)=0. ! 计算直线段左端点处κ1,κ2,κ3,κ4 AI(N1,2)=0. AI(N1,3)=Y1(N1+1)-Y1(N1) AI(N1,4)=Y1(N1) 40 IF(N1-L1.EQ.0)GOTO 100 ! 滑裂面上没有直线段 83 B(L1)=1.5*(Y1(L1+1)-Y1(L1))/(X1(L1+1)-X1(L1)) A(L1)=-0.5 ! 插值区间左端点处的一阶导数值 CN=3.0*(Y1(N1)-Y1(N2))/(X1(N1)-X1(N2)) IF(N1-L1.EQ.1)GOTO 55 ! 插值区间只具有左右两个端点 L2=L1+1 DO 50 J=L2,N2 H1=X1(J)-X1(J-1) H=X1(J+1)-X1(J) AF=H1/(H+H1) ! 计算α值 BT=3.*((1.-AF)*(Y1(J)-Y1(J-1))/H1+AF*(Y1(J+1)-Y1(J))/H) !计算β值 A(J)=-AF/(2.+(1.-AF)*A(J-1)) ! 计算插值区间内插值节点处的一阶导数值 mj 50 B(J)=(BT-(1.-AF)*B(J-1))/(2.+(1.-AF)*A(J-1)) 55 A(N1)=(CN-B(N2))/(2.+A(N2)) 60 A(N2)=A(N2)*A(N2+1)+B(N2) ! 计算插值区间右端点处的一阶导数值 mn