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

第11章程序设计345 IF(N2,GE. LDGOTO 60 N2=NI-1 DO 80 I=LlN H=Xl(+1)-X1(D) AI(1,4)=Y1()!计算插值区间内插值节点处的x1,K2,K3,K (1,3)=A()H (1,2)=3.·(Y1(1+1)Y1())(2*A(I)H+A(+1)*H I(,1)=2·(Yl(I}Yl(+1)+(A(+A(+1)H 100 CONTINUE 0o200l=1,NSl-1 LI=KQ2) NI=KQ2(1+1) STEP=1/NI-Ll) Sx=X1(H+1)Xl(1) DO 200 J=LLNI 2(=X1()+S1·SX!计算滑裂面与土条侧面边界线交点的X,y值 Y2(J)=AI(,4)+S1*(AI(,3)HS1°(AI(,2H+SlAI(,1) 1=S1+STEP NW=KQ2(NSI WRITE(6, N=,NW WRITE(6, 105 M, X2(M), Y2(M), M=I, NW) 105 FORMAT(1X,1s,2F10.3) STOP END 3.光滑曲线滑裂面例题 (1)数据文件 Spline-ldat: 1-51522-35.000 6-46.314-28.235 计算成果 Spline-1 151 2-50480-33.466 -47.356-29.317 6-46.314-28.235 7-45.15727.280 8-44.000-26.545 9-42.842-25.975 10-41.685-25.511 l1-40.528-25.10 13-38.213-24.495 16-34.718-24.243 17-33.539-24.352 19-31.180-24.754 030.000-25.000

第 11 章 程序设计 345 N2=N2-1 IF(N2.GE.L1)GOTO 60 N2=N1-1 DO 80 I=L1,N2 H=X1(I+1)-X1(I) AI(I,4)=Y1(I) ! 计算插值区间内插值节点处的κ1,κ2,κ3,κ4 AI(I,3)=A(I)*H AI(I,2)=3.*(Y1(I+1)-Y1(I))-(2.*A(I)+A(I+1))*H 80 AI(I,1)=2.*(Y1(I)-Y1(I+1))+(A(I)+A(I+1))*H 100 CONTINUE DO 200 I=1,NS1-1 L1=KQ2(I) N1=KQ2(I+1) STEP=1./(N1-L1) SX=X1(I+1)-X1(I) S1=0. DO 200 J=L1,N1 X2(J)=X1(I)+S1*SX ! 计算滑裂面与土条侧面边界线交点的 X Y 值 Y2(J)=AI(I,4)+S1*(AI(I,3)+S1*(AI(I,2)+S1*AI(I,1))) 200 S1=S1+STEP NW=KQ2(NS1) WRITE(6,*)'N= ',NW WRITE(6,105)(M,X2(M),Y2(M),M=1,NW) 105 FORMAT(1X,I5,2F10.3) STOP END 3. 光滑曲线滑裂面例题 (1) 数据文件 Spline-1.dat 5 1 -51.522 -35.000 6 -46.314 -28.235 10 -41.685 -25.511 15 -35.898 -24.224 20 -30.000 -25.000 0 计算成果 Spline-1.c N= 20 1 -51.522 -35.000 2 -50.480 -33.466 3 -49.439 -31.978 4 -48.397 -30.580 5 -47.356 -29.317 6 -46.314 -28.235 7 -45.157 -27.280 8 -44.000 -26.545 9 -42.842 -25.975 10 -41.685 -25.511 11 -40.528 -25.107 12 -39.370 -24.765 13 -38.213 -24.495 14 -37.055 -24.311 15 -35.898 -24.224 16 -34.718 -24.243 17 -33.539 -24.352 18 -32.359 -24.530 19 -31.180 -24.754 20 -30.000 -25.000

346土质边坡德定分析一原理 4.曲线-直线组合滑裂面例题 (1)数据文件 Spline-2dat 51.522-35.000 6-46.314-28.235 10-41.68525.511 1535.898 4.224 20-30.000-25.000 (2)计算成果 Spline-2c 151.522-35.00 2-50.480-33.459 6-46.314-28.235 7-45.157-27.300 8-44000-26.583 11-40.528-25511 12-39370-25.511 13-38.213-25.511 14-37.055-25.511 15-35.898-25.511 16-34.718-25 17-33.539-25.307 1832.359-25204 19-31.180-25.102 20-30.000-25.000 11.3.4计算土条重量和物理力学参数的程序 PROFILE FOR 1.说明 本子程序根据112.1介绍的方法计算每个土条的重量、条底的强度指标和孔隙水压力a 2.源程序 DIMENSION SH(20), M4(20), SHNE(20), M4NE(20) DIMENSION XU(80), HSL(80), X3(80), Y3(80) (10), LINES(80), YSF(80)B3/KD COMMON/AI/IN, INI, PF(20), PC(20), PDW(20), PDS(20), IPH, ICPH(40, 3) COMMON/A2/NN, NJ(20), NHY, DHRP, WL(80), BDA(80), CA(80), FAA(80) COMMON/A3/RW, RUL, UWL/SWEAK/TWEAK, KQ3(50), LWK, IC3(40, 3) COMMON/A6/X(80,Y(80,XN(99),YN(9),X2(80,Y2(80),IC(80,3),N COMMON/A7/(80), RU(80), C(80), F(80)/A8/MYK, KXSX, MTE, KL COMMON/A1S/UWLI, NCOUN, MLA(80), PF1(20), PCI(20), MID(20), PNOR(20) NHT, MRD INTEGER4 OPTION(6),OP1(6) COMMON/OPP/OPTION OPI OPEN(5, FILE= STATUS=UNKNOWN) PEN(6, FILE= STATUS=UNKNOWN) C IwR5=0,1分别相应顶部有、无拉力缝; YTENSION=拉力缝底Y坐标 DS.CX.CY=滑弧深度,圆心坐标

346 土质边坡稳定分析 原理 ⋅ 方法 ⋅ 程序 4. 曲线-直线组合滑裂面例题 (1) 数据文件 Spline-2.dat 5 1 -51.522 -35.000 6 -46.314 -28.235 10 -41.685 -25.511 15 -35.898 -24.224 20 -30.000 -25.000 1 3 (2) 计算成果 Spline-2.c N= 20 1 -51.522 -35.000 2 -50.480 -33.459 3 -49.439 -31.965 4 -48.397 -30.565 5 -47.356 -29.306 6 -46.314 -28.235 7 -45.157 -27.300 8 -44.000 -26.583 9 -42.842 -26.011 10 -41.685 -25.511 11 -40.528 -25.511 12 -39.370 -25.511 13 -38.213 -25.511 14 -37.055 -25.511 15 -35.898 -25.511 16 -34.718 -25.409 17 -33.539 -25.307 18 -32.359 -25.204 19 -31.180 -25.102 20 -30.000 -25.000 11. 3. 4 计算土条重量和物理力学参数的程序 PROFILE.FOR 1. 说明 本子程序根据 11.2.1 介绍的方法计算每个土条的重量 条底的强度指标和孔隙水压力 2. 源程序 DIMENSION SH(20),M4(20),SHNE(20),M4NE(20) DIMENSION XU(80),HSL(80),X3(80),Y3(80) COMMON/CRU/RUS(10),LINES(80),YSF(80)/B3/IKD COMMON/A1/IN,IN1,PF(20),PC(20),PDW(20),PDS(20),IPH,ICPH(40,3) COMMON/A2/NN,NJ(20),NHY,DHRP,WL(80),BDA(80),CA(80),FAA(80) #,LRU(20) COMMON/A3/RW,RU1,UWL/SWEAK/IWEAK,KQ3(50),LWK,IC3(40,3) COMMON/A6/X(80),Y(80),XN(99),YN(99),X2(80),Y2(80),IC(80,3),N COMMON/A7/W(80),RU(80),C(80),F(80)/A8/MYK,KXSX,MTE,KL COMMON/A15/UWL1,NCOUN,MLA(80),PF1(20),PC1(20),MID(20),PNOR(20) # ,NHT,MRD INTEGER*4 OPTION(6),OP1(6) COMMON/OPP/OPTION,OP1 OPEN(5,FILE=' ',STATUS='UNKNOWN') OPEN(6,FILE=' ',STATUS='UNKNOWN') C IWR5=0,1 分别相应顶部有 无拉力缝 YTENSION=拉力缝底 Y 坐标 C DS,CX,CY=滑弧深度 圆心坐标

第11章程序设计347 READ(S, N !土条总数 READ(5,)RW,UWL!水容重,坡外水位 READ5,*)IN!线段总数 DO 3021=L IN READ(5,*)(C(D)J=1,3)!线段信息 302 CONTINUE !点总数 READ,*)(L,XN(D),YN①)=1,NN)!点坐标 DOI=LN READ(,)I,(D,Y ENDDO IF(PH EQ.O)GOTO DO I=LIPH READ5,*)(ICPH(,D),J=1,3)!线段信息 18 READ(S, )IN, ( L, PF(), PC), PF1(), PCI(), PDW(I), PDS) %LRU IFl, IN) CONTINUE 对土条信息重新处理,保证IC,1)的x坐标比IC(1,2)的小 DO 37=l.INI M2=C(I,2) IF(XN(MI). LTXN(M2)GO TO 37 IC(, 2FMI 37 CONTINUE F(PH EQ0)GO TO 672 C REARRANGE ICPH(I )IN THE SAME WAY AS IC(, J) DO 38 I=LIPH MI=ICPH(I, 1) M2=ICPH(, 2) IF(XN(MD).LTXN(M2))GO TO 38 PH(, IFICPH(, 2) H(,2)=MI 38 CONTINUE 672 CONTINUE GTENSION=O FIND THE LINES WHICH INTERCEPT THE CONSIDERED SLICE MX:IC(. YYY=YN(MX) DO 813=lINI DO813Il=1,2 Y=IC(, D IF(YNGY)GT.YYYYYY=YNY 813 CONTINUE YYY=YYY+100 X3(1=X2(1) Y3(1)=Y2(1) DO211K=1 QN=0. W(k=0 DHREO WL(KFYYY

第 11 章 程序设计 347 READ(5,*) N ! 土条总数 READ(5,*) RW,UWL ! 水容重 坡外水位 READ(5,*) IN1 ! 线段总数 DO 302 I=1,IN1 READ(5,*) (IC(I,J),J=1,3) !线段信息 302 CONTINUE READ(5,*)NN ! 点总数 READ(5,*) (L,XN(I),YN(I),I=1,NN) ! 点坐标 DO I=1,N READ(5,*) I,X(I),Y(I) ENDDO READ(5,*)IPH IF(IPH.EQ.0)GOTO 18 DO I=1,IPH READ(5,*) (ICPH(I,J),J=1,3) !线段信息 ENDDO 18 READ(5,*)IN,(L,PF(I),PC(I),PF1(I),PC1(I),PDW(I),PDS(I), %LRU(I),I=1,IN) CONTINUE C 对土条信息重新处理 保证 IC(I,1)的 x 坐标比 IC(I,2)的小 DO 37 I=1,IN1 M1=IC(I,1) M2=IC(I,2) IF(XN(M1).LT.XN(M2))GO TO 37 IC(I,1)=IC(I,2) IC(I,2)=M1 37 CONTINUE IF(IPH.EQ.0)GO TO 672 C REARRANGE ICPH(I,J) IN THE SAME WAY AS IC(I,J) DO 38 I=1,IPH M1=ICPH(I,1) M2=ICPH(I,2) IF(XN(M1).LT.XN(M2))GO TO 38 ICPH(I,1)=ICPH(I,2) ICPH(I,2)=M1 38 CONTINUE 672 CONTINUE GTENSION=0 C FIND THE LINES WHICH INTERCEPT THE CONSIDERED SLICE MX=IC(1,1) YYY=YN(MX) DO 813 I=1,IN1 DO 813 II=1,2 IY=IC(I,II) IF(YN(IY).GT.YYY)YYY=YN(IY) 813 CONTINUE YYY=YYY+100. X3(1)=X2(1) Y3(1)=Y2(1) DO 211 K=1,N MF=0 QN=0. W(K)=0. DHR=0 WL(K)=YYY

348土质边坡德定分析一原理·方法·程序 LPHS=IC(IN1, 3) IF(IPHEQ.O)GO TO 7 DO 73511=1.IPH MI=ICPH(Il, 1) M2=ICPH(1, 2) IF(X(K).LTXN(MD)GO TO 73: IF(X(K).GE XN(M2))GO TO 735 L(K=YN(MI)+X(K)XN(MI*(YNM2)YN(MI)(XN(M2)XN(MI) LPHS=ICPH(1, 3 73 CONTINUE CONTINUE PP=0 FT=0 l=0 DO 20711=LINI MI=IC(l, 1) M2=lC(,2) M3=I1 IF(X(K).GEXN(M2)GO TO 207 SH2=YNMl+(X(K)XN(MI))*(YN(M2YN(MIXN(M2}XNMI) 152 FORMAT(1X,**15,2F15.6) DSH2=ABS((Y(K)-SH2)) DSHI=SH2-Y(K) ASSURE THAT THE STRENGTH PARAMETERS OF THE UNDERLYING SOIL IF(DSH2. GT..00005AND DSHIGT.0.001GO TO 207 SH(KIFSH2 M4(K1FM3 207 CONTINUE IF(KI. NE IGO TO 707 GO TO 708 707LX=M4(1) HN=SH(I) DO45L1=2.K1 IF(SH(LI).GT HN)GO TO 45 H(L1) HN=SHLl) 45 CONTINUE 708Ml=C(LX,1) M2=IC(LX, 2) DELXXN(M2)-XN(MD IF(ABS(DELX)GT0.001)BBK=YN(M2)-YN(MI)(XN(M2)XN(MI X3(K)=x2(K) Y3(KFYNMI-+(X3(K)-XN(MD))BBK 139 FORMAT(IX, KI, LX, M2, MI, XM2, XMI, 415, 2F10.3) LINES(KFLX !96 12.29 47 DO 46 LI=LKI LSB=M4(LI M4(LIFIC(LSB, 3) CONTINU IF(WL(K). GT.Y(K)GO TO 734 DO 530 ISP=l K SDIF=ABS(WL(K)-SH(ISP)) IF(SDIF.LT0.00001GO TO 734

348 土质边坡稳定分析 原理 ⋅ 方法 ⋅ 程序 LPHS=IC(IN1,3) IF(IPH.EQ.0)GO TO 739 DO 735 I1=1,IPH M1=ICPH(I1,1) M2=ICPH(I1,2) IF(X(K).LT.XN(M1))GO TO 735 IF(X(K).GE.XN(M2))GO TO 735 IF(ABS(XN(M2)-XN(M1)).LT.0.001)GOTO 735 WL(K)=YN(M1)+(X(K)-XN(M1))*(YN(M2)-YN(M1))/(XN(M2)-XN(M1)) LPHS=ICPH(I1,3) 735 CONTINUE 739 CONTINUE PP=0. CT=0. FT=0. K1=0 DO 207 I1=1,IN1 M1=IC(I1,1) M2=IC(I1,2) M3=I1 IF(X(K).LT.XN(M1))GO TO 207 IF(X(K).GE.XN(M2))GO TO 207 SH2=YN(M1)+(X(K)-XN(M1))*(YN(M2)-YN(M1))/(XN(M2)-XN(M1)) 152 FORMAT(1X,'**',I5,2F15.6) DSH2=ABS((Y(K)-SH2)) DSH1=SH2-Y(K) C ASSURE THAT THE STRENGTH PARAMETERS OF THE UNDERLYING SOIL IF(DSH2.GT.0.00005.AND.DSH1.GT.0.001)GO TO 207 K1=K1+1 SH(K1)=SH2 M4(K1)=M3 207 CONTINUE IF(K1.NE.1)GO TO 707 LX=M4(1) GO TO 708 707 LX=M4(1) HN=SH(1) DO 45 L1=2,K1 IF(SH(L1).GT.HN)GO TO 45 LX=M4(L1) HN=SH(L1) 45 CONTINUE 708 M1=IC(LX,1) M2=IC(LX,2) DELX=XN(M2)-XN(M1) IF(ABS(DELX).GT.0.001)BBK=(YN(M2)-YN(M1))/(XN(M2)-XN(M1)) X3(K)=X2(K) Y3(K)=YN(M1)+(X3(K)-XN(M1))*BBK 139 FORMAT(1X,'K1,LX,M2,M1,XM2,XM1',4I5,2F10.3) LINES(K)=LX !96.12.29 47 DO 46 L1=1,K1 LSB=M4(L1) M4(L1)=IC(LSB,3) 46 CONTINUE IF(WL(K).GT.Y(K))GO TO 734 DO 530 ISP=1,K1 SDIF=ABS(WL(K)-SH(ISP)) IF(SDIF.LT.0.00001)GO TO 734

第11章程序设计349 530 CONTINUE KI=KI+I SH(KIFWL(K) M4(KI=LPHS 734 MSP=l DO209l=1K1 IF(SHO).LTSHD)GO TO 209 MSP=MSP+ 209 CONTINUE SHNECMSPFSHO) M4NE(MSPFM40J IFO.GE KI)GO TO 229 GO TO 27 9 Kl=KI+ SHNE(KlFY(K) YSF(KFSHNE(I) M4NE(KI FM4NE(KI-1) SH3=SHNE(K1)-SHNE(D) F((KI-1).LEO)GO TO 641 IF(SHNE(D).GTUWL)W(KFRW"(SHNE(1)-UWL) DO 2081=2.K M5=M4NE(I-I Ll=WL(K)-0.0001 IF(SHNE(I-1). LT. WLDPDNO=PDW(M5) IF(SHNE(I-1).GE.WLIPDNO=PDS(M5) DWK=DH·PDNO W(K户=W(K)+DWK 42 CT=CT+DH"PC(M5) IF(PC(M5).LT0.00001.AND PF(M5).LT0.00001)GO TO 589 FT=FT+DH"PF(M5) 589 DDH=DDH+DH DHR=DHR+DH 08 CONTINUE DHRPESH3-DHR M7=M4NE(KD) IF(ME.NEO)M7=IC3 (MF, 3) MLA(KFM PP=0 XU(KSHNE(1) HSL(K)SH3 IF(IPH EQ. O)GO TO 736 (Y(K).LTWL(K)GOTO 590 WLM=WL(K) F(KEQ1OR K EQ NGO TO 590 FORU(M7).EQ0 AND SHNE(1).LEUWL)PP=(Y(K)-WLM)*RW FRU(M7).EQ0 AND SHNE(1).GTUWL)PP=(Y(K)-UWL*RW 590 CONTINUE WTl=Y(K)-UWL IF(WTI. LTOWTI=O W(KFW(K)-WTI*RW PP=PPWT1·RW 736 IF(ABS(SH3).LT0.001.OR ABS(W(K).LT0.001GO TO 641

第 11 章 程序设计 349 530 CONTINUE K1=K1+1 SH(K1)=WL(K) M4(K1)=LPHS 734 J=1 27 MSP=1 DO 209 I=1,K1 IF(J.EQ.I)GO TO 209 IF(SH(J).LT.SH(I))GO TO 209 MSP=MSP+1 209 CONTINUE SHNE(MSP)=SH(J) M4NE(MSP)=M4(J) IF(J.GE.K1)GO TO 229 J=J+1 GO TO 27 229 K1=K1+1 SHNE(K1)=Y(K) YSF(K)=SHNE(1) M4NE(K1)=M4NE(K1-1) SH3=SHNE(K1)-SHNE(1) IF((K1-1).LE.0)GO TO 641 DDH=0. IF(SHNE(1).GT.UWL)W(K)=RW*(SHNE(1)-UWL) DO 208 I=2,K1 M5=M4NE(I-1) WL1=WL(K)-0.0001 IF(SHNE(I-1).LT.WL1)PDNO=PDW(M5) IF(SHNE(I-1).GE.WL1)PDNO=PDS(M5) DH=SHNE(I)-SHNE(I-1) DWK=DH*PDNO W(K)=W(K)+DWK 42 CT=CT+DH*PC(M5) IF(PC(M5).LT.0.00001.AND.PF(M5).LT.0.00001)GO TO 589 FT=FT+DH*PF(M5) GO TO 208 589 DDH=DDH+DH DHR=DHR+DH 208 CONTINUE DHRP=SH3-DHR M7=M4NE(K1) IF(MF.NE.0)M7=IC3(MF,3) MLA(K)=M7 PP=0 XU(K)=SHNE(1) HSL(K)=SH3 IF(IPH.EQ.0)GO TO 736 IF(Y(K).LT.WL(K))GOTO 590 WLM=WL(K) IF(K.EQ.1.OR.K.EQ.N)GO TO 590 IF(LRU(M7).EQ.0.AND.SHNE(1).LE.UWL)PP=(Y(K)-WLM)*RW IF(LRU(M7).EQ.0.AND.SHNE(1).GT.UWL)PP=(Y(K)-UWL)*RW 590 CONTINUE WT1=Y(K)-UWL IF(WT1.LT.0)WT1=0 W(K)=W(K)-WT1*RW PP=PP-WT1*RW 736 IF(ABS(SH3).LT.0.001.OR.ABS(W(K)).LT.0.001)GO TO 641