📄 s.txt
字号:
REAL AB1
" INTEGER N,NS,OPTION"
" COMMON/AS/ASP"
" COMMON/CIRCLE/ID_C,CX,CY,DS"
"C OPTION=0, SPENCER法"
"C OPTION=1, BISHOP法"
"C OPTION=2, 瑞典法"
"C OPTION=3, 工程师团法"
"C OPTION=4, 简化法"
"C ID_C=0, 任意形状滑裂面； ID_C不为零，圆弧滑裂面"
" REAL X2(80),Y2(80),W(80),X(80),Y(80),C(80),F(80),WL(80)"
"$,FXO(80),FX(80),BN(80)"
" REAL BDA(80),FQUH(80),FDIS(80),YSF(80),RU(80),ALF(80),RW"
" write(6,*)'输入数据文件名'"
"OPEN(5,FILE=' ',STATUS='UNKNOWN')"
" write(6,*)'输入计算成果文件名'"
"OPEN(6,FILE=' ',STATUS='UNKNOWN')"
" RW=9.8"
" AB1=0.0"
" WRITE(6,*)'边坡稳定分析'"
" READ(5,*)OPTION,ID_C"
" IF(ID_C.NE.0)READ(5,*)CX,CY,DS"
" READ(5,*) NS !土条总数"
" READ(5,*)(X2(I),I=1,NS)"
" READ(5,*)(Y2(I),I=1,NS)"
" N=NS+1"
" DO 302 I=1,N"
" READ(5,*)W(I),RU(I),C(I),F(I),BDA(I),WL(I),FQUH(I),FDIS(I),YSF(I)"
302 CONTINUE
"!读入两端点和土条底中点的数据,其中涉及力的物理量均以单宽土条计"
" !X2,Y2=坐标,W=土条重量,RU=孔压系数,C=凝聚力=,F=磨擦系数,BDA=土容重,"
" !WL=浸润线Y坐标,FQUH,FDIS,YSF=水平地震力的大小,Y坐标,垂直水平地震力的大小"
IF(OPTION.EQ.3)THEN
" READ(5,*)ASP"
" WRITE(6,*)'陆军工程师团法，平均坝坡=',ASP"
" ASP=ASP*3.14159/180."
" ENDIF "
"WRITE(6,709)RW"
"709 FORMAT(1X,'水容重=',F10.3//)"
"WRITE(6,702)"
"702 FORMAT(T25,'DATA FOR ASSUMED SIDE FORCE FUNCTION')"
"WRITE(6,703)"
"703 FORMAT(T25,'************************************'//)"
NS=N-1
" !计算土条底中点的坐标和条底倾角"
"DO 208 I=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)))
208 CONTINUE
ALF(1)=ALF(2)
ALF(N)=0
X(N)=X2(NS)
Y(N)=Y2(NS)
X(1)=X2(1)
Y(1)=Y2(1)
"WRITE(6,716)"
"716 FORMAT(T5,'NO. ',T14,'X',T24,'Y',T34,'W',T44,'RU',T54,'C',"
"#T64,'F',T74,'BDA',T84,'WL',T94,'QH',T104,'YE')"
"DO 901 I=1,N"
" WRITE(6,718) I,X(I),Y(I),W(I),RU(I),C(I),F(I),BDA(I),WL(I)"
" #,FQUH(I),FDIS(I)"
"901 CONTINUE"
"718 FORMAT(1X,I5,11F10.3)"
"WRITE(6,719)"
"719 FORMAT(1X,/'NOTE:')"
"WRITE(6,720)"
"720 FORMAT(5X,'W=土条的单宽重量')"
"WRITE(6,721)"
"721 FORMAT(5X,'RU=孔隙水压力系数')"
"WRITE(6,722)"
"722 FORMAT(5X,'C=凝聚力')"
"WRITE(6,723)"
"723 FORMAT(5X,'F=摩擦系数')"
"WRITE(6,724)"
"724 FORMAT(5X,'BDA=土条的平均容重')"
"WRITE(6,725)"
"725 FORMAT(5X,'WL=浸润线的Y坐标')"
"c WRITE(6,701)"
"c 701 FORMAT(T5,'N0. ',T14,'AB1')"
"c WRITE(6,*)'迭代过程'"
" SELECT CASE (OPTION)"
" CASE(0)"
"WRITE(6,701)"
"701 FORMAT(/10X,'SPENCER METHOD')"
"CALL SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,4)"
" IWALL=0"
"DO 301 I=1,80"
FXO(I)=0
FX(I)=1
"301 CONTINUE !tan(BETA(I))=FXO(I)+ALAM*FX(I)"
KXYX=0
"CALL MP(N,KXYX,KS,AB1,DQ,DM,BN,FX,FXO,"
"&C,F,X2,Y2,X,Y,W,RU,FQUH,WL,FDIS,RW,DB,DF)"
" CASE(1)"
"WRITE(6,704)"
"704 FORMAT(/20X,'BISHOP法'/)"
"CALL SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION)"
" CASE(2)"
"WRITE(6,705)"
"705 FORMAT(/20X,'瑞典法'/)"
"CALL SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION)"
" CASE(3)"
"WRITE(6,731)"
"731 FORMAT(/10X,'工程师团法')"
"CALL SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,4)"
" IWALL=0"
"DO 311 I=1,80"
FXO(I)=0
FX(I)=1
"311 CONTINUE !tan(BETA(I))=FXO(I)+ALAM*FX(I)"
KXYX=1
"CALL MP(N,KXYX,KS,AB1,DQ,DM,BN,FX,FXO,"
"&C,F,X2,Y2,X,Y,W,RU,FQUH,WL,FDIS,RW,DB,DF) "
" CASE(4)"
"WRITE(6,789)"
"789 FORMAT(/20X,'简化法法'/)"
"CALL SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,OPTION)"
" CASE DEFAULT"
" END SELECT"
END
" SUBROUTINE SAFE2(X2,F,FQUH,FDIS,W,RU,N,X,ALF,C,AB1,ID)"
" REAL X2(80),F(80),FQUH(80),W(80),RU(80),ALF(80),X(80),C(80)"
" REAL AB1,FDIS(80)"
"C INTEGER N"
" REAL RO1,DM,RM,BX1,CDN,DG(80),G3,DX1,G,GX"
" REAL O,S1,FRIC,S2,S,DRM,P1,P2,DDM,DBX1,DCDN"
" COMMON/CIRCLE/ID_C,CX,CY,DS"
INTEGER NS
" IWALL=0"
" RO1=1.0"
NS=N-1
KL=2
DM=0.
RM=0.
BX1=0.
CDN=0.
DG(2)=0.
G3=0.
"DO 93 I=2,NS"
DX1=X2(I)-X2(I-1)
"C DO 79 J=1,NS+1"
C IF(I.EQ.NOO(J))THEN
C DG(I+1)=DG(I)-F(I+1)*ALF(I+1)+F(I)*ALF(I)
C GOTO 93
C ENDIF
C 79 CONTINUE
DG(I+1)=DG(I)+(F(I)-F(I+1))*ALF(I)
93 G3=G3+(F(I)*ALF(I)+DG(I))*DX1
" G3=G3/(X(N)-X(1))"
"DO 202 I=2,NS"
IF(ID_C.EQ.0)THEN
" RO1=1"
" ELSE"
" RO1=(FDIS(I)-CY)/(DS-CY)"
" ENDIF"
G=ALF(I)
GX=G*F(I)-G3+DG(I)
DX1=X2(I)-X2(I-1)
O=FQUH(I)*SIN(G)
S1=W(I)*(COS(G)-RU(I)/COS(G))
FRIC=F(I)
S2=(S1-O)*FRIC
" S=C(I)/COS(G)"
DRM=S+S2
DRM=DRM*DX1
P1=W(I)*SIN(G)
P2=FQUH(I)*RO1
DDM=(P1+P2)*DX1
DM=DM+DDM
RM=RM+DRM
DBX1=-DRM-SIN(G)*GX*W(I)*DX1
BX1=BX1+DBX1
DCDN=DRM*GX
CDN=CDN+DCDN
202 CONTINUE
"AB1=-BX1/DM+CDN/BX1 !通过简化法求得AB1的迭代初值"
" IF(ID.EQ.4)WRITE(6,*)'简化法安全系数= ',AB1"
" IF(ID.EQ.2)AB1=RM/DM"
" IF(ID.EQ.2)WRITE(6,*)'瑞典法安全系数= ',AB1"
" IF(ID.NE.1)RETURN"
" WRITE(6,*)'迭代过程'"
LM=1
41 DT=0.
DM=0.
RM=0.
"DO 21 I=2,NS"
G=ALF(I)
DX1=X2(I)-X2(I-1)
S1=W(I)*(1-RU(I))*F(I)
S2=S1+C(I)
S=COS(G)*(1+F(I)*TAN(G)/AB1)
DRM=DX1*S2/S
RM=RM+DRM
P1=W(I)*SIN(G)
P2=FQUH(I)*RO1
DDM=(P1+P2)*DX1
DM=DM+DDM
DT=DT+DRM*TAN(G)*F(I)/(AB1**2+AB1*F(I)*TAN(G))
21 CONTINUE
"AB2=RM/DM !BISHOP方法"
IF(AB2.LT.0.OR.AB2.GE.15) THEN
"WRITE(6,*)'F.S NOT REASONABLE = ',AB2"
AB2=9999.9
RETURN
ENDIF
" WRITE(6,740) LM,AB2"
"740 FORMAT(1X,I5,F10.5)"
IF(ABS(AB2-AB1).LT.0.0001)GO TO 49
IF(LM.LT.10)GO TO 9191
"WRITE(6,9192)"
"9192 FORMAT(1X,'ITERATION TROUBLES IN BISHOP METHOD')"
GO TO 49
"9191 CK=1-DT/DM !牛顿迭代法"
AB1=AB1-(AB1-AB2)/CK
LM=LM+1
GO TO 41
"49 WRITE(6,203)AB1"
"203 FORMAT(/5X,'BISHOP 法安全系数=',F10.4)"
RETURN
END
"SUBROUTINE DETE(AB1,RW,DQ,DM,ALAM,KS,BN,FX,FXO,"
"&C,F,X2,Y2,W,RU,FQUH,YSF,WL,FDIS,N,DB,KXYX,DF)"
"REAL AB1,DQ,DT,DM,ALAM,DF"
"REAL C(80),F(80),X2(80),Y2(80),W(80),RU(80),FXO(80),FX(80)"
"REAL FQUH(80),YSF(80),WL(80),FDIS(80),BN(80)"
"INTEGER KS,N,KXYX"
"REAL X(80),ALF(80),G,Y(80)"
"REAL EXT,EXR,EYR,EMT,EMR,FMT,DMT,DEG,DX1"
"REAL DMG,BS,R1X,TX,QX,PT112,R1XT,POX"
"REAL POXT,POY,DPT112,DPX112,BSG,H3,YT(80),GF(80)"
"REAL GFX(80),GFY(80),BET(80),EFX(80)"
"REAL G1,TA,TA1,OP,FI,FI1,CE,OQ,OQ1,THO,CW"
"REAL SEGB,SEGB1,FI0,AXS11,OTEM,DDQ,PX112,DBS"
"REAL SX,SXG,DTX1,DTX2,DTX,TX1,DDP,DBS1,DBS2,EPX112"
"REAL AAKX1,DR1X1,DR1X2,DR1X,DR1XT"
"REAL DG,DDT,GAVE,BETA,DGX,OX,OY,HSL"
"REAL TAA,OQS2,OQS1,DG1,BT(80),EXR1,EXT1,EYR1"
"REAL EYT1,XCE,YCE,XRE,YRE,DR2X11,DR2X12,DR2X1"
"REAL GS,DPOY1,DPOY2,OY1,DDM,DDMT,PM,PMT,DDEG,DDMG"
"REAL OQ0,SEGB0,DGG,DMG1,DEG1,OXX"
"REAL TEMP1,TEMP2,GG,DDM2,DM1,EWALL,GWALL,PX"
"REAL DQM,DMM,DTM,DMTM,DEGM,DMGM,S1,HMW,DDM1,ETA"
"REAL GTENSION,HITE,WP,RW"
"INTEGER EYT,IWALL,NS"
" !本程序计算Morgenstern-Price,Spencer和工程师团法"
"IWALL=0. !土压力的作用位置"
" HMW=0."
GTENSION=0.
HITE=0.
" GWALL=0. !土压力的大小"
" !GWALL=1 计算土压力设定条件下的安全系数"
" !GWALL=2 计算主动土压力"
" IF(KS.EQ.0)THEN !KS调用本程序段标志符,=0表示第1次调用"
"WRITE(6,705)"
"705 FORMAT(1X,/'VALUE OF FX0(I) IN ORDER OF I')"
"WRITE(6,706)(FXO(I),I=1,N)"
"706 FORMAT(1X,5F10.5)"
"WRITE(6,707)"
"707 FORMAT(1X,/'VALUE OF FX(I) IN ORDER OF I')"
"WRITE(6,708)(FX(I),I=1,N)"
"708 FORMAT(1X,5F10.5)"
"WRITE(6,726)"
"726 FORMAT(1X,/'迭代过程:'/)"
"WRITE(6,712)"
"712 FORMAT(T10,'NO.',T24,'DQ',T35,'DM',T55,'AB1',T65,'ALAM')"
" ENDIF"
" KS=1"
DM=0.
DT=0.
EXT=gtension
EYT=0
EXR=0.
EYR=0.
EMT=-gtension*hite/3.
EMR=0.
FMT=0.
DQ=0.
DMT=0.
DEG=0.
DMG=0.
BS=0.
R1X=0.
TX=0.
QX=0.
PT112=0.
PX112=0.
R1XT=0.
POX=0.
POXT=0.
POY=0.
DPT112=0.
DPX112=0.
BSG=0.
H3=HITE/3.
YT(1)=Y2(1)-H3
DQ=DQ-GTENSION
GF(1)=-DQ
DM=DM-GTENSION*H3!96.8.31
" NS=N-1"
"DO 303 I=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)))
303 CONTINUE
ALF(1)=ALF(2)
ALF(N)=0.
X(N)=X2(NS)
Y(N)=Y2(NS)
"X(1)=X2(1) !X,Y各土条的中点坐标(在滑裂面上)"
Y(1)=Y2(1)
"C 601 IF(KXYX.NE.0) GOTO 401 !KXYX=0,采用Morgenstern-Price方法"
"DO 304 I=1,N"
" BET(I)=ATAN(FXO(I)+ALAM*FX(I))"
IF(C(I).LT.0.001.AND.F(I).LT.0.001) BET(I)=0
304 CONTINUE
401 CONTINUE
"EWALL=BET(NS) !土压力的作用角度"
GFX(1)=GF(1)*COS(BET(1))
GFY(1)=GF(1)*SIN(BET(1))
"DO 305 I=2,NS"
PX=X(I)-X(I-1)
DX1=X2(I)-X2(I-1)
G=ALF(I)
G1=ALF(1)
TA=BET(I)
TA1=BET(1)
OP=G-TA
FI=ATAN(F(I)/AB1)
FI1=ATAN(F(1)/AB1)
CE=C(I)/AB1
OQ=FI-G
OQ1=FI1-G1
THO=SIN(TA)*X(I)-COS(TA)*Y(I)-(SIN(TA)*X(1)-COS(TA)*Y(1))
CW=1
IF(I.GT.2)CW=(X2(I-1)-X2(I-2))/(X2(I)-X2(I-1))
IF((OQ+TA).GT.1.5.OR.(OQ1+TA1).GT.1.5.OR.TA.GT.
$ 1.5.OR.TA1.GT.1.5) GOTO 317
SEGB=(1./COS(OQ+TA))**2
SEGB1=(1./COS(OQ1+TA1))**2
FMT=FMT+FQUH(I)*(Y(I)-FDIS(I))*DX1
FI0=ATAN(F(I+1)/AB1)
AXS11=ABS((FI-FI0))
OTEM=ABS(OQ+TA)*180./3.14259
IF(OTEM.LT.95.AND.OTEM.GT.85)RETURN
DDQ=(W(I)*SIN(OQ)+CE*COS(FI)/COS(G)-RU(I)*W(I)*SIN(FI)/COS(G)
"& -FQUH(I)*COS(OQ))/COS(OQ+TA) !DDQ=p(x)sec(FI-ALF-BETA)"
IF(I.GT.2) GOTO 307
PX112=-TAN(OQ1+TA1)*COS(TA1)**2.*(FX(1))
DBS1=TAN(OQ1+TA1)
307 DBS2=TAN(OQ+TA)
EPX112=-DBS2*COS(TA)**2.*(FX(I))
DBS=(DBS1*CW+DBS2)/(1+CW)
DBS1=DBS2
BS=BS+DBS*(BET(I)-BET(I-1))
BSG=BS+DBS*(BET(I)-BET(I-1))*(0.5*DX1)/PX
IF((-BS).GT.30.OR.BSG.GT.30)GOTO 317
"SX=2.71828**(-BS) !SX=s(x)*cos(FI-ALF+BETA)"
SXG=2.71828**BSG
IF(I.GT.2) GO TO 309
DTX1=(SIN(TA1)-TAN(G1)*COS(TA1))
309 DTX2=(SIN(TA)-TAN(G)*COS(TA))/SX
DTX=(DTX1*CW+DTX2)/(1+CW)
DTX1=DTX2
TX1=TX+DTX*PX*0.5
TX=TX+DTX*(X(I)-X(I-1))
DDP=(-FQUH(I)*SIN(TA)*SIN(FI)+CE*COS(OP)*COS(FI)/
# COS(G)+W(I)*SIN(FI)*COS(TA)-RU(I)*W(I)*SIN(FI)*COS(OP)/
# COS(G))*COS(FI)/(-AB1*COS(OQ+TA)**2)
AAKX1=DDP/DDQ
IF(I.GT.2) GO TO 310
DR1X1=(-SIN(FI1)*COS(FI1)*SEGB1/AB1)
310 DR1X2=(-SIN(FI )*COS(FI )*SEGB/AB1)
DR1X=(DR1X1*CW+DR1X2)/(1+CW)
DR1X1=DR1X2
DR1XT=DR1X*TX1
R1XT=R1XT+DR1XT*(BET(I)-BET(I-1))
R1X =R1X +DR1X *(BET(I)-BET(I-1))
DG=DDQ*SX*DX1
"DQ=DQ+DG !DQ=DG的积分,DG=p(x)s(x),力的平衡"
DDT=DG*(AAKX1-R1X)
"DT=DT+DDT !DT=DDT的积分,DDT=p(x)s(x)t(x),力矩平衡"
GF(I)=-DQ*SXG
GFX(I)=GF(I)*COS(BET(I))
GFY(I)=GF(I)*SIN(BET(I))
GAVE=(GFY(I)+GFY(I-1))/2.
BETA=(BET(I)+BET(I-1))/2.
DGX=GFX(I)-GFX(I-1)
OX=(GFX(I-1)+GFX(I))/2.
OY=(GFY(I-1)+GFY(I))/2.
IF(WL(I).GE.Y(I))THEN
WP=0.
ELSE
WP=0.5*(ru(i)*w(i))**2/rw
ENDIF
OXX=OX-WP
EFX(I)=OXX
IF(OXX.GT.0.) HSL=Y(I)-YSF(I)
TAA=COS(TA)*COS(TA)*FX(I)
OQS2=(BET(I+1)+BET(I))/2
OQS1=(BET(I-1)+BET(I))/2
DG1=GF(I)*SIN(G-OQS2)-GF(I-1)*SIN(G-OQS1)
BN(I)=W(I)*DX1*COS(G)+DG1-FQUH(I)*DX1*SIN(G)
BT(I)=(BN(I)-RU(I)*W(I)*DX1/COS(G))*TAN(FI)+CE*DX1/COS(G)
EXR1=BN(I)*SIN(G)-BT(I)*COS(G)
EXR=EXR+EXR1
EXT1=FQUH(I)*DX1
EXT=EXT+EXT1
EYR1=-BN(I)*COS(G)-BT(I)*SIN(G)
EYR=EYR+EYR1
EYT1=W(I)*DX1
EYT=EYT+EYT1
XCE=X(1)
YCE=Y(1)
XRE=X(I)-XCE
YRE=Y(I)-YCE
EMT=EMT+EXT1*(FDIS(I)-YCE)-EYT1*XRE !moment of body force
EMR=EMR+EXR1*YRE-EYR1*XRE !moment of internal force
if(abs(dx1).lt.0.001)goto 317
BN(I)=BN(I)*COS(G)/DX1
BT(I)=BT(I)*COS(G)/DX1
IF(KXYX.EQ.0) GOTO 402
" GOTO 305"
402 IF(I.GT.2) GO TO 312
DR2X11=SEGB1*COS(TA1)**2.*(FX(1))
312 DR2X12=SEGB *COS(TA )**2.*(FX(I))
DR2X1=(DR2X11*CW+DR2X12)/(1+CW)
DR2X11=DR2X12
GS=ALF(I-1)
POX=POX-DR2X1*(ALF(I)-GS)
POXT=POXT-DR2X1*(ALF(I)-GS)*TX1
IF(I.GT.2) GOTO 313
DPOY1=(COS(FI1)/(COS(G1)*COS(OQ1+TA1)))*COS(TA1)**2.*(FX(1))
313 DPOY2=(COS(FI)/(COS(G)*COS(OQ+TA)*SX))*COS(TA)**2.*(FX(I))
OY1=(DPOY1*CW+DPOY2)/(1+CW)
DPOY1=DPOY2
POY=POY+OY1*PX
DDM=DG*TX
DDMT=DG*(AAKX1*TX-R1XT)
PX112=PX112+DPX112
PT112=PT112+DPT112
DPX112=0
DPT112=0.
PM=POX-PX112
PMT=POXT+POY-PT112
DDEG=DG*PM
DDMG=DG*PMT
DM=DM+DDM
DMT=DMT+DDMT
DEG=DEG+DDEG
DMG=DMG+DDMG
"IF(AXS11-0.005) 305,305,315"
G=ALF(I+1)
315 OQ0=FI0-G
IF((OQ0+TA).GT.1.5) GOTO 317
SEGB0=(1./COS(OQ0+TA)**2)
DR2X11=SEGB0*COS(TA)**2.*(FX(I))
DPX112=-TAN(FI0+BET(I+1)-G)*COS(BET(I+1))**2*FX(I+1)
# +TAN(OQ+TA)*COS(TA)**2.*FX(I)
DPT112=DPX112*TX
305 CONTINUE
TEMP1=EPX112+PM
TEMP2=PMT-TX*TAN(OQ+TA)*TAA
IF(IWALL.EQ.2)THEN
GG=-DQ
DDM1=HMW*COS(EWALL)+TX*SX
DDM2=HMW*COS(EWALL)/SX+TX
DM1=DM
DM=DM1+GG/SX*DDM1-FMT
DDDM2=-HMW*COS(EWALL)/SX*TEMP1-TX*TEMP1+TEMP2
DGG=-DEG*(HMW*COS(EWALL)/SX+TX)
DMG1=DMG
DEG1=DEG
DMG=DMG1+GG*(-HMW*COS(EWALL)/SX*TEMP1-HMW/SX*SIN(EWALL)
"# *FX(N)*COS(EWALL)*COS(EWALL)-TX*TEMP1+TEMP2)-"
# DEG*(HMW*COS(EWALL)/SX+TX)
DB=-DM/DMG
ETA=1.0+DQ/GWALL/SX
EXT=EXT-GWALL*COS(EWALL)*(1-ETA)
EYT=EYT-GWALL*SIN(EWALL)*(1-ETA)
EMT=EMT+GWALL*SIN(EWALL)*(X(N)-X(1))*(1-ETA)-
% GWALL*COS(EWALL)*(Y(N)-Y(1)-HMW)*(1-ETA)
GOTO 317
ENDIF
IF(IWALL.EQ.1)THEN
DQM=SX*GWALL
DDM1=HMW*COS(EWALL)+TX*SX
DMM=GWALL*DDM1
DTM=-GWALL*SX*R1X
DMTM=-GWALL*SX*R1XT
DEGM=TEMP1*SX*GWALL
DMGM=TEMP2*SX*GWALL
DMGM=DMGM-GWALL*SIN(EWALL)*FX(N)*COS(EWALL)
"$ *COS(EWALL)*HMW"
DQ=DQ+DQM
DM=DM+DMM
"DT=DT+DTM !DT=DQ对AB1的偏导数"
"DMT=DMT+DMTM !DMT=DM对AB1的偏导数"
"DEG=DEG+DEGM !DEG=DQ对ALAM的偏导数"
"DMG=DMG+DMGM !DMG=DM对ALAM的偏导数"
EXT=EXT-GWALL*COS(EWALL)
EYT=EYT-GWALL*SIN(EWALL)
EMT=EMT+GWALL*SIN(EWALL)*(X(N)-X(1))-GWALL*COS(EWALL)
$ *(Y(N)-Y(1)-HMW)
ENDIF
316 CONTINUE
if(abs(dm).gt.1.0e+20)goto 317
if(abs(dmt).gt.1.0e+20)goto 317
"DM=DM-FMT !DQ=Gn !DM=Mn"
S1=DEG*DMT-DT*DMG
" IF(KXYX.EQ.1)THEN"
" DF=-DQ/DT"
" GOTO 317"
" ENDIF"
IF(ABS(S1).LE.0.01)GOTO 317
DF=(DQ*DMG-DM*DEG)/S1
DB=(-DQ*DMT+DM*DT)/S1
317 RETURN
END
C
" SUBROUTINE MP(N,KXYX,KS,AB1,DQ,DM,BN,FX,FXO,"
"&C,F,X2,Y2,X,Y,W,RU,FQUH,WL,FDIS,RW,DB,DF)"
"REAL AB1,DF,DQ,DM,ALAM,DB"
" INTEGER KS,N,KXYX"
" REAL X2(80),Y2(80),W(80),X(80),Y(80),C(80),F(80),WL(80)"
" REAL FQUH(80),FDIS(80),YSF(80),RU(80),RW"
" REAL GAV,PLA,FA,FB,PLA1"
" REAL BN(80),FX(80),FXO(80)"
" COMMON/AS/ASP"
" !KXYX力矩平衡控制符,=0时同时计算力和力矩的平衡"
" FA=0."
" FB=0."
" PLA=0."
" IF(KXYX.EQ.0)THEN"
GAV=(Y(N)-Y(1))/(X(N)-X(1))
"DO 277 I=2,N"
PLA=PLA+(X(I)-X(I-1))*(FX(I)+FX(I-1))/2.
"277 CONTINUE"
PLA1=(FA+FB)*(X(N)-X(1))/2.
" ALAM=(GAV*(X(N)-X(1))-PLA1)/PLA !ALAM的初值,参见Chen&Morgenstern(1983)"
" ELSE"
" ALAM=ASP"
" ENDIF"
" MQ=1"
15 CONTINUE
"CALL DETE(AB1,RW,DQ,DM,ALAM,KS,BN,FX,FXO,"
"&C,F,X2,Y2,W,RU,FQUH,YSF,WL,FDIS,N,DB,KXYX,DF)"
"WRITE(6,*)MQ,DQ,DM,AB1,ALAM"
" AB1=AB1+DF"
IF(KXYX.EQ.0)ALAM=ALAM+DB
" IF(MQ.GE.20)GO TO 44"
" MQ=MQ+1 "
" IF(ABS(DF).GT.0.1E-4)GO TO 15"
IF(ABS(DB).GT.0.1E-3)GO TO 15
44 CONTINUE
IF(MQ.GE.20)THEN
" WRITE(6,*)'THE ITERATION FAILS'"
" RETURN"
" ENDIF"
"WRITE(6,727)"
"727 FORMAT(1X,//'THE RESULT IS:')"
"WRITE(6,710)ALAM"
"710 FORMAT(1X,'LAMBDA=',F12.6)"
"WRITE(6,711)AB1"
"711 FORMAT(1X,'FACTOR OF SAFETY=',F12.6/)"
RETURN
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -