***************************************************************
PROGRAM
PROG4
********************************************************************************
* PROGRAM FOR PUMP RUNDOWN FROM POWER
FAILURE
* AND SUBSEQUENT RESTART. PROVISION IS
MADE FOR
* COLUMN SEPARATION, AIR-VACUUM
VALVES, AND
* ONE-WAY VACUUM VALVES. A BYPASS LINE
IS
* AROUND THE PUMPS.
*
********************************************************************************
*
DIMENSION
X[ALLOCATABLE](:,:),VL[ALLOCATABLE](:,:),
$VR[ALLOCATABLE](:,:),H[ALLOCATABLE](:,:),HEAD[ALLOCATABLE](:,:),
$HHIGH[ALLOCATABLE](:,:),HLOW[ALLOCATABLE](:,:),
$VNEWL[ALLOCATABLE](:,:),VNEWR[ALLOCATABLE](:,:),
$HNEW[ALLOCATABLE](:,:),HSTEAD[ALLOCATABLE](:,:),
$THIGH[ALLOCATABLE](:,:),TLOW[ALLOCATABLE](:,:),
$CAV[ALLOCATABLE](:,:),HCAV[ALLOCATABLE](:,:),
$CAVMAX[ALLOCATABLE](:,:),TCAV[ALLOCATABLE](:,:)
DIMENSION
L[ALLOCATABLE](:),D[ALLOCATABLE](:),
$A[ALLOCATABLE](:),PIPEZ[ALLOCATABLE](:),
$F[ALLOCATABLE](:),VZERO[ALLOCATABLE](:),
$NPAR[ALLOCATABLE](:),C[ALLOCATABLE](:),AK[ALLOCATABLE](:),
$SINE[ALLOCATABLE](:),AREA[ALLOCATABLE](:),
$DELTT[ALLOCATABLE](:),RATIO[ALLOCATABLE](:),
$AIRVAC[ALLOCATABLE](:),CTZERO[ALLOCATABLE](:),
$BUBBLE[ALLOCATABLE](:)
DIMENSION
TSTART(10),TRAMP(10),RPMSTRT(10),RPMEND(10),RPMS(10),
$QP(10)
DIMENSION
QN(6),HNSQ(6),TNSQ(6),IOU(6),TIOU(6)
$,QNN(6),HNN(6),C7(6),C8(6),IZ(6)
LOGICAL
FAIL,BYPASS,SAVE,PPLOT,HVPRNT,RERUN,RSTART
LOGICAL BUBBLE,GRAPH,PFILE,NOGOOD
INTEGER PIPE,AIRVAC
REAL LA,L,NEXP,KLB
CHARACTER
TITLE*80,TITLE2*80,VACYES*3,VACNO*3,AV*3
CHARACTER CH,CH12,CH78,CH79,CHN,FNAME*12
DATA
SAVE,PPLOT,HVPRNT,RERUN,PFILE,GRAPH/6*.FALSE./
*
NAMELIST /SPECS/
NPIPES,NPARTS,HRES,HSUMP,ZEND,HATM,QTRY,QACC,
$
TMAX,DTNEW,HVPRNT,PPLOT,PFILE,GRAPH,RERUN,
$ DB,KLB
NAMELIST /PUMPS/ NPUMPS,NSTAGE,RPM,RPMZ,WRSQ,NSTART,QN,HNSQ,TNSQ
NAMELIST /RESTART/
TSTART,TRAMP,RPMSTRT,RPMEND
NAMELIST /OUTCTRL/ IOU,TIOU
*
*------------------------------------------------------------------------------
WRITE(*,791)
READ(*,790) FNAME
OPEN(5,FILE=FNAME)
READ(5,100) TITLE
READ(5,100) TITLE2
READ(5,SPECS)
NP=NPIPES
ALLOCATE
(L(NP),D(NP),A(NP),PIPEZ(NP),F(NP),VZERO(NP),NPAR(NP),
$C(NP),AK(NP),SINE(NP),AREA(NP),DELTT(NP),RATIO(NP),AIRVAC(NP),
$BUBBLE(NP),CTZERO(NP))
READ(5,*)
(PIPE,D(I),L(I),F(I),A(I),PIPEZ(I),AIRVAC(I),I=1,NPIPES)
READ(5,PUMPS)
IF(NSTART.GT.0) READ(5,RESTART)
READ(5,OUTCTRL)
WRITE(*,792)
READ(*,790) FNAME
OPEN(6,FILE=FNAME,STATUS='NEW')
*------------------------------------------------------------------------------
*
III=2
ICOUNT=0
IOUT=IOU(1)
CH12=CHAR(12)
CH=CHAR(27)
CH78=CHAR(78)
CH79=CHAR(79)
CHN=CHAR(3)
DELTN=100000.
RPMZ=RPM
VACYES='YES'
VACNO=' NO'
RSTART=.FALSE.
FAIL=.FALSE.
BYPASS=.FALSE.
IF(PPLOT.OR.HVPRNT.OR.GRAPH) SAVE=.TRUE.
TTTT=10000.
DO 303 I=1,NSTART
RPMS(I)=0.
QP(I)=0.
IF(TSTART(I).LT.TTTT) TTTT=TSTART(I)
303 CONTINUE
PINERT=WRSQ/32.2
PI=3.141592
COEF=RPMZ/RPM
DO 302 I=1,6
QNN(I)=QN(I)*COEF
302
HNN(I)=HNSQ(I)*COEF*COEF
NEXP=1.0
IF(F(1).GT.10.) NEXP=0.85
*
* COMPUTE STEADY
STATE DISCHARGE
IF(QTRY.LT.0.001) QTRY=QN(4)
COEF=12.*144.*144.*16.*NPUMPS*NPUMPS/(64.4*449.*449.*PI*PI)
IF(F(1).GT.10.)
COEF=3.03*(12.**4.87)*(NPUMPS/(449.*.7854))**1.85
EX=2.0
IF(F(1).GT.10.) EX=1.85
DO 320 NZ=1,20
DO 300 I=1,5
IF(QTRY.GT.QN(I).AND.QTRY.LE.QN(I+1)) GO
TO 310
300 CONTINUE
WRITE(6,250)
STOP
310
HN=HNSQ(I)+(QTRY-QN(I))*(HNSQ(I+1)-HNSQ(I))/(QN(I+1)-QN(I))
SUM=0.
DO 301 J=1,NPIPES
IF(F(1).GT.10.)
SUM=SUM+L(J)*QTRY**1.85/((F(J)**1.85)*
$(D(J)**4.87))
IF(F(1).LT.10.)
SUM=SUM+F(J)*L(J)*QTRY*QTRY/D(J)**5
301 CONTINUE
FUNCT=HSUMP-HRES+NSTAGE*HN-COEF*SUM
FPRIME=NSTAGE*(HNSQ(I+1)-HNSQ(I))/(QN(I+1)-QN(I))-2.*COEF*SUM/QTRY
QNEXT=QTRY-FUNCT/FPRIME
IF(ABS(QNEXT-QTRY).LT.QACC) GO TO 350
QTRY=QNEXT
320 CONTINUE
WRITE(6,251)
STOP
350
QLINE=QNEXT*NPUMPS
Q=QNEXT
HPUMP=NSTAGE*HN+HSUMP
COEF=4.*144./(449.*PI)
DO 2 I=1,NPIPES
VZERO(I)=QLINE*COEF/D(I)**2
AREA(I)=0.25*PI*D(I)*D(I)/144.
C(I)=32.2/A(I)
DELTT(I)=L(I)/(NPARTS*(2.*VZERO(I)+A(I)))
IF(I.EQ.NPIPES)
SINE(I)=(ZEND-PIPEZ(I))/L(I)
IF(I.EQ.NPIPES) GO TO 2
SINE(I)=(PIPEZ(I+1)-PIPEZ(I))/L(I)
2 CONTINUE
*
* ** COMPUTE
PUMP BYPASS PARAMETERS **
*
AB=.25*PI*DB**2/144.
CCA=(KLB*AREA(1)**2)/(64.4*AB**2)
CCB=1.0/C(1)
*
* ** COMPUTE
MINIMUM DELTA T **
*
DELT=DELTT(1)
KMIN=1
DO 12 I=2,NPIPES
IF(DELTT(I).GT.DELT) GO TO 12
DELT=DELTT(I)
KMIN=I
12 CONTINUE
IF(RERUN) DELT=DTNEW
DO 13 I=1,NPIPES
IF(I.EQ.KMIN) GO TO 14
ANPAR=L(I)/(DELT*(VZERO(I)+A(I)))
NPAR(I)=ANPAR
TEST=NPAR(I)
IF(ABS(TEST+1-ANPAR).LT..0001)
NPAR(I)=NPAR(I)+1
GO TO 166
14 NPAR(I)=NPARTS
166
IF(NPAR(I).GT.NPARR) NPARR=NPAR(I)
13 CONTINUE
IX=NPARR+1
IZZ=IX+1
ALLOCATE(X(NP,IX),H(NP,IX),HLOW(NP,IX),HHIGH(NP,IX),HEAD(NP,IZZ),
$HNEW(NP,IX),HSTEAD(NP,IX),VL(NP,IX),VR(NP,IX),VNEWL(NP,IX),
$VNEWR(NP,IX),CAV(NP,IX),HCAV(NP,IX),CAVMAX(NP,IX),THIGH(NP,IX),
$TLOW(NP,IX),TCAV(NP,IX))
DO 177 I=1,NPIPES
K=NPAR(I)+1
177
HEAD(I,K+1)=0.
*
WRITE(6,1011) CH,CH78,CHN
WRITE(6,201) TITLE
WRITE(6,2011) TITLE2
WRITE(6,202)
WRITE(6,101)
NPIPES,NPARTS,HRES,HSUMP,HATM,ZEND,TMAX,QACC,DB,KLB
WRITE(6,1010)
(IOU(I),I=1,6),(TIOU(I),I=1,6)
WRITE(6,104)
WRITE(6,105) NPUMPS,NSTAGE,RPMZ,WRSQ,RPM
WRITE(6,106)
WRITE(6,107)
(QN(I),HNSQ(I),TNSQ(I),I=1,6)
WRITE(6,108)
(I,TSTART(I),RPMSTRT(I),RPMEND(I),TRAMP(I),
$I=1,NSTART)
HPUMPP=HPUMP-HSUMP
WRITE(6,1088) QLINE,HPUMPP
WRITE(6,1011) CH12
IF(F(1).LT.10.) WRITE(6,102)
IF(F(1).GT.10.) WRITE(6,1022)
DO 3 I=1,NPIPES
IF(F(1).LT.10.) WRITE(6,103)
I,D(I),L(I),A(I),PIPEZ(I),F(I),
$VZERO(I)
IF(F(1).GT.10.) WRITE(6,1033)
I,D(I),L(I),A(I),PIPEZ(I),F(I),
$VZERO(I)
3 CONTINUE
WRITE(6,1020)
DO 4 I=1,NPIPES
AV=VACNO
IF(AIRVAC(I).GE.1) AV=VACYES
LA=L(I)/A(I)
AINT=1.0-DELT*NPAR(I)/LA
4 WRITE(6,1030)
I,DELTT(I),NPAR(I),SINE(I),LA,AINT,AV
WRITE(6,1040) DELT
*
* CONVERT TO Q/N, H/NSQ, T/NSQ
COEF=60.*550./(2.*PI)
DO 1 I=1,6
QN(I)=QN(I)/(449.*RPM)
TNSQ(I)=NSTAGE*COEF*TNSQ(I)/RPM**3
1 HNSQ(I)=HNSQ(I)/RPM**2
*
*
* ** SET UP
CONDITIONS FOR STEADY STATE (T=0) **
*
H(1,1)=HPUMP
VL(1,1)=100000.
VNEWL(1,1)=VL(1,1)
VR(1,1)=VZERO(1)
DO 20 I=1,NPIPES
DELL=L(I)/NPAR(I)
IF(F(I).LT.10.) AK(I)=12.*F(I)*DELT/(2.*D(I))
IF(F(I).GT.10.)
AK(I)=12.*DELT*195./(2.*D(I)*(F(I)**1.85)*
$(D(I)/12.)**.17)
IF(F(I).GT.10.)
F(I)=195./((F(I)**1.85)*(VZERO(I)**.15)*
$((D(I)/12.)**.17))
DELHF=12.*F(I)*DELL*VZERO(I)**2/(64.4*D(I))
K=NPAR(I)+1
X(I,1)=0.
HEAD(I,1)=H(I,1)-PIPEZ(I)
HCAV(I,1)=PIPEZ(I)-HATM
RATIO(I)=DELT/DELL
CAV(I,1)=0.
CAVMAX(I,1)=0.
TLOW(I,1)=0.
THIGH(I,1)=0.
TCAV(I,1)=0.
BUBBLE(I)=.FALSE.
DO 21 J=2,K
X(I,J)=(J-1)*DELL/L(I)
H(I,J)=H(I,1)-DELHF*(J-1)
HEAD(I,J)=H(I,J)-(PIPEZ(I)+X(I,J)*L(I)*SINE(I))
HCAV(I,J)=PIPEZ(I)+X(I,J)*L(I)*SINE(I)-HATM
CAV(I,J)=0.
CAVMAX(I,J)=0.
TLOW(I,J)=0.
THIGH(I,J)=0.
TCAV(I,J)=0.
VL(I,J)=VZERO(I)
VR(I,J)=VZERO(I)
IF(J.EQ.K) VR(I,J)=100000.
IF(J.EQ.K) VNEWR(I,J)=VR(I,J)
21 CONTINUE
IF(I.EQ.NPIPES) GO TO 20
H(I+1,1)=H(I,K)
VR(I+1,1)=VZERO(I+1)
VL(I+1,1)=100000.
VNEWL(I+1,1)=VL(I+1,1)
20 CONTINUE
DO 23 I=1,NPIPES
IF(AIRVAC(I).GE.1) HCAV(I,1)=PIPEZ(I)
IF(I.GT.1) KK=NPAR(I-1)+1
IF(I.GT.1.AND.AIRVAC(I).GE.1)
HCAV(I-1,KK)=PIPEZ(I)
K=NPAR(I)+1
DO 23 J=1,K
HSTEAD(I,J)=H(I,J)
HLOW(I,J)=H(I,J)
23 HHIGH(I,J)=H(I,J)
PHMAX=-100.
PHMIN=100000.
DO 11 I=1,NPIPES
K=NPAR(I)+1
DO 11 J=1,K
IF(HEAD(I,J).LT.PHMAX) GO TO 16
IPMAX=I
XMAX=X(I,J)
PHMAX=HEAD(I,J)
TTMAX=T
GO TO 17
16 IF(HEAD(I,J).GT.PHMIN) GO TO 17
IPMIN=I
XMIN=X(I,J)
PHMIN=HEAD(I,J)
TTMIN=T
17 CONTINUE
11 CONTINUE
*
* ** WRITE OUT
STEADY STATE CONDITIONS **
*
T=0.0
IENTRY=0
INDEX=TMAX/DELT+1
WRITE(6,1011) CH12
WRITE(6,204)
WRITE(6,205) T
DO 22 I=1,NPIPES
K=NPAR(I)+1
WRITE(6,206)
I,(X(I,J),HEAD(I,J),H(I,J),VL(I,J),
$VR(I,J),CAV(I,J),J=1,K)
22 CONTINUE
II=1
IF(SAVE) CALL
GRAFCAV(II,NP,IX,INDEX,DELT,VL,VR,HEAD,CAV,
$NPAR,IENTRY)
IENTRY=1
HPUMP=HPUMP-HSUMP
WRITE(6,212) RPM,Q,HPUMP
Q=Q/449.
*
*-------------------------------------------------------------------------------
* BEGIN TRANSIENT ANALYSIS
*-------------------------------------------------------------------------------
*
DO 99 II=1,INDEX
T=T+DELT
*
* ** COMPUTE H
AND V AT INTERIOR NODES **
*
DO 30 I=1,NPIPES
K=NPAR(I)
DO 30 J=2,K
RA=RATIO(I)*A(I)
VLEFT=VL(I,J)-RA*(VL(I,J)-VR(I,J-1))
VRITE=VR(I,J)-RA*(VR(I,J)-VL(I,J+1))
HLEFT=H(I,J)-RA*(H(I,J)-H(I,J-1))
HRITE=H(I,J)-RA*(H(I,J)-H(I,J+1))
IF (CAV(I,J).GT.0.) GO TO 31
VNEWL(I,J)=0.5*(VLEFT+VRITE+C(I)*(HLEFT-HRITE)+C(I)*DELT*SINE(I)
$*(VLEFT-VRITE)-AK(I)*(VLEFT*ABS(VLEFT)**NEXP+VRITE*
$ABS(VRITE)**NEXP))
HNEW(I,J)=0.5*(HLEFT+HRITE+(VLEFT-VRITE)/C(I)+DELT*SINE(I)*
$(VLEFT+VRITE)-(AK(I)/C(I))*(VLEFT*ABS(VLEFT)**NEXP-VRITE*
$ABS(VRITE)**NEXP))
VNEWR(I,J)=VNEWL(I,J)
IF (HNEW(I,J).LT.HCAV(I,J)) GO TO 31
GO TO 30
31
C1=VLEFT+C(I)*HLEFT-AK(I)*VLEFT*ABS(VLEFT)**NEXP+C(I)*DELT*VLEFT*
$SINE(I)
C3=VRITE-C(I)*HRITE-AK(I)*VRITE*ABS(VRITE)**NEXP-C(I)*DELT*VRITE*
$SINE(I)
VNEWL(I,J)=C1-C(I)*HCAV(I,J)
VNEWR(I,J)=C3+C(I)*HCAV(I,J)
HNEW(I,J)=HCAV(I,J)
CAV(I,J)=(VNEWR(I,J)-VNEWL(I,J))*DELT*AREA(I)+CAV(I,J)
IF (CAV(I,J).GT. 0.) GO TO 30
DELTAH=ABS(0.5*A(I)*(VNEWR(I,J)-VNEWL(I,J))/32.2)
HNEW(I,J)=HCAV(I,J)+DELTAH
VNEWL(I,J)=0.5*(VNEWL(I,J)+VNEWR(I,J))
VNEWR(I,J)=VNEWL(I,J)
CAV(I,J)=0.
30 CONTINUE
*
* ** COMPUTE H
AND V AT INTERIOR JUNCTIONS **
*
KK=NPIPES-1
DO 40 I=1,KK
K=NPAR(I)+1
VLEFT=VL(I,K)-RATIO(I)*A(I)*(VL(I,K)-VR(I,K-1))
HLEFT=H(I,K)-RATIO(I)*A(I)*(H(I,K)-H(I,K-1))
CCC=VLEFT+C(I)*HLEFT+C(I)*VLEFT*DELT*SINE(I)-AK(I)*VLEFT*
$ABS(VLEFT)**NEXP
VRITE=VR(I+1,1)-RATIO(I+1)*A(I+1)*(VR(I+1,1)-VL(I+1,2))
HRITE=H(I+1,1)-RATIO(I+1)*A(I+1)*(H(I+1,1)-H(I+1,2))
CC=VRITE-C(I+1)*HRITE-C(I+1)*DELT*VRITE*SINE(I+1)-AK(I+1)*VRITE*
$ABS(VRITE)**NEXP
IF (CAV(I,K).GT.0.) GO TO 41
HNEW(I,K)=(-CC*AREA(I+1)+CCC*AREA(I))/(C(I+1)*AREA(I+1)+C(I)*
$AREA(I))
IF(HNEW(I,K).LT.HCAV(I,K)) GO TO 41
HNEW(I+1,1)=HNEW(I,K)
VNEWL(I,K)=CCC-C(I)*HNEW(I,K)
VNEWR(I+1,1)=CC+C(I+1)*HNEW(I+1,1)
GO TO 40
41 VNEWL(I,K)=CCC-C(I)*HCAV(I,K)
VNEWR(I+1,1)=CC+C(I+1)*HCAV(I+1,1)
HNEW(I,K)=HCAV(I,K)
HNEW(I+1,1)=HCAV(I+1,1)
DELCAV=DELT*(VNEWR(I+1,1)*AREA(I+1)-VNEWL(I,K)*AREA(I))
IF(DELCAV.GT.0.) GO TO 43
IF(AIRVAC(I+1).EQ.2.AND.HCAV(I,K).LE.PIPEZ(I+1).AND.
$DELCAV.LT.0.) CTZERO(I+1)=CAV(I,K)
IF(CTZERO(I+1).LT.(5.*AREA(I+1))) GO TO
43
IF(AIRVAC(I+1).EQ.2.AND.HCAV(I,K).LE.PIPEZ(I+1).AND.
$DELCAV.LT.0.) BUBBLE(I+1)=.TRUE.
43 CAV(I,K)=CAV(I,K)+DELCAV
CAV(I+1,1)=CAV(I,K)
IF(.NOT. BUBBLE(I+1)) GO TO 42
HCAV(I,K)=PIPEZ(I+1)-HATM+HATM*(CTZERO(I+1)/CAV(I,K))**1.2
IF(HCAV(I,K).LT.PIPEZ(I+1))
BUBBLE(I+1)=.FALSE.
IF(HCAV(I,K).LT.PIPEZ(I+1)) HCAV(I,K)=PIPEZ(I+1)
HCAV(I+1,1)=HCAV(I,K)
IF(BUBBLE(I+1)) GO TO 40
42 IF(CAV(I,K).GT.0.) GO TO 40
DELH2=ABS((0.5*A(I+1)/32.2)*((AREA(I)/AREA(I+1))*VNEWL(I,K)-
$VNEWR(I+1,1)))
HNEW(I,K)=HCAV(I,K)+DELH2*2.0*AREA(I+1)*A(I)/(A(I+1)*AREA(I)
$+AREA(I+1)*A(I))
HNEW(I+1,1)=HNEW(I,K)
VNEWL(I,K)=VNEWL(I,K)-(32.2/A(I))*(HNEW(I,K)-HCAV(I,K))
VNEWR(I+1,1)=VNEWL(I,K)*AREA(I)/AREA(I+1)
CAV(I,K)=0.
CAV(I+1,1)=0.
40 CONTINUE
*
* ** COMPUTE H
AND V AT DOWNSTREAM END **
* --THIS BOUNDARY CONDITION IS FOR A
CONSTANT HEAD RESERVOIR--
*
NV=NPAR(NPIPES)+1
VLEFT=VL(NPIPES,NV)-RATIO(NPIPES)*A(NPIPES)*(VL(NPIPES,NV)-
$VR(NPIPES,NV-1))
HLEFT=H(NPIPES,NV)-RATIO(NPIPES)*A(NPIPES)*(H(NPIPES,NV)-
$H(NPIPES,NV-1))
CC=VLEFT+C(NPIPES)*HLEFT+C(NPIPES)*VLEFT*DELT*SINE(NPIPES)-
$AK(NPIPES)*VLEFT*ABS(VLEFT)**NEXP
HNEW(NPIPES,NV)=HRES
VNEWL(NPIPES,NV)=CC-C(NPIPES)*HNEW(NPIPES,NV)
*
* ** COMPUTE H AND V AT UPSTREAM END **
*
* BOUNDARY CONDITIONS AT THE PUMP
*
IF(T.GT.TTTT) RSTART=.TRUE.
IF(RSTART) BYPASS=.FALSE.
IF(RSTART) GO TO 554
IF(BYPASS) GO TO 554
556 QQ=Q/RPM
IF(QQ.LE.0.) GO TO 552
DO 550 I=1,5
IF(QQ.GT.QN(I).AND.QQ.LE.QN(I+1)) GO TO
551
550 CONTINUE
ICOUNT=ICOUNT+1
IF(ICOUNT.LT.5) WRITE(6,260)
I=5
TN2=TNSQ(6)
GO TO 553
551 IF(RSTART)
GO TO 554
TN2=TNSQ(I)+(QQ-QN(I))*(TNSQ(I+1)-TNSQ(I))/(QN(I+1)-QN(I))
GO
TO 553
552 TN2=TNSQ(1)
I=1
553
RPM=RPM-30.*TN2*DELT/(PI*PINERT)*RPM*RPM
IF(RPM.LT.0.) BYPASS=.TRUE.
554
VRITE=VR(1,1)-RATIO(1)*A(1)*(VR(1,1)-VL(1,2))
HRITE=H(1,1)-RATIO(1)*A(1)*(H(1,1)-H(1,2))
C1=VRITE-C(1)*HRITE-C(1)*DELT*VRITE*SINE(1)-AK(1)*VRITE*
$ABS(VRITE)**NEXP
*
* IF PUMP RESTART HAS BEGUN, EXIT THE POWER
FAILURE MODE OF
* OPERATION
*
IF(RSTART) GO TO 650
IF(.NOT.BYPASS) GO TO 555
CCC=-HSUMP-C1/C(1)
VPP=(-CCB+SQRT(CCB**2-4.*CCA*CCC))/(2.*CCA)
IF(VPP.GT.0.) GO TO 557
VPP=0.
HPP=-C1/C(1)
GO TO 530
557
HPP=(VPP-C1)/C(1)
GO TO 530
555 DO 500
NZ=1,20
HA=HNSQ(I)
HB=HNSQ(I+1)
QA=QN(I)
QB=QN(I+1)
C55=NSTAGE*RPM*RPM
C77=(HA-HB)/(QA-QB)
C88=HB-QB*C77
C66=HSUMP+C55*C88+C55*C77*C1*AREA(1)/(NPUMPS*RPM)
C99=C55*C77*C(1)*AREA(1)/(NPUMPS*RPM)
HPP=C66/(1.-C99)
IF(HPP.GT.HSUMP) GO TO 558
CCC=-HSUMP-C1/C(1)
VPP=(-CCB+SQRT(CCB**2-4.*CCA*CCC))/(2.*CCA)
IF(VPP.GT.0.) GO TO 559
VPP=0.
HPP=-C1/C(1)
GO TO 530
559
HPP=(VPP-C1)/C(1)
BYPASS=.TRUE.
GO TO 530
558
VPP=C1+HPP*C(1)
IF(VPP.LE.0.) VPP=0.
IF(VPP.LE.0.) HPP=-C1/C(1)
QPP=VPP*AREA(1)/(NPUMPS*RPM)
DO 510 J=1,5
IF(QPP.GT.QN(J).AND.QPP.LE.QN(J+1)) GO TO
520
510 CONTINUE
IF(QPP.GT.0..AND.I.EQ.5) GO TO 530
IF(QPP.GT.0..AND.I.LT.5) I=I+1
IF(QPP.LE.0..AND.I.EQ.1) GO TO 530
IF(QPP.LE.0..AND.I.GT.1) I=I-1
GO TO 500
520 IF(J.EQ.I)
GO TO 530
I=J
500 CONTINUE
WRITE(6,270)
STOP
530
HNEW(1,1)=HPP
VNEWR(1,1)=VPP
IF(BYPASS) GO TO 531
IF(VPP.GT.0.) HPUMP=HPP-HSUMP
IF(HPP.LT.HSUMP) HPUMP=0.
IF(HPP.LT.HSUMP) BYPASS=.TRUE.
IF(VPP.LE.0.)
HPUMP=HNSQ(1)*RPM*RPM*NSTAGE
Q=QPP*RPM
GO TO 531
*......................................................................
* PUMPS RESTARTING
*......................................................................
650 CONTINUE
ICHEK=0
*
* COMPUTE SPEED
OF EACH STARTUP PUMP
*
DO 651 I=1,NSTART
IF(T.LT.TSTART(I)) GO TO 651
NACT=I
IF(T.GT.(TSTART(I)+TRAMP(I))) GO TO 652
RPMS(I)=RPMSTRT(I)+(T-TSTART(I))*(RPMEND(I)-RPMSTRT(I))/TRAMP(I)
GO TO 651
652 RPMS(I)=RPMEND(I)
651 CONTINUE
*
* COMPUTE
LOCATION ON PUMP CURVES FOR FIRST TRY
*
DO 655 J=1,NSTART
IF(RPMS(J).LE.0.) GO TO 655
QP(J)=QP(J)/RPMS(J)
DO 656 I=1,5
IF(QP(J).GE.QN(I).AND.QP(J).LE.QN(I+1))
GO TO 657
656 CONTINUE
I=5
IF(QP(J).LE.0.) I=1
657 IZ(J)=I
HA=HNSQ(I)
HB=HNSQ(I+1)
QA=QN(I)
QB=QN(I+1)
C7(J)=(HA-HB)/(QA-QB)
C8(J)=HB-QB*C7(J)
655 CONTINUE
*
* COMPUTE TERMS
REQUIRED TO GENERATE SOLUTIONS FOR PUMP STARTUP
*
658 SUM7=0.
SUM8=0.
DO 660 I=1,NACT
IF(RPMS(I).LE.0.) GO TO 660
SUM7=SUM7+1.0/(C7(I)*RPMS(I))
SUM8=SUM8+RPMS(I)*C8(I)/C7(I)
660 CONTINUE
*
* FIND SOLUTION
FOR PUMP STARTUP
*
C15=C(1)*AREA(1)-SUM7/NSTAGE
C16=-C1*AREA(1)-HSUMP*SUM7/NSTAGE-SUM8
HPP=C16/C15
VPP=C1+C(1)*HPP
NACT2=NACT
DO 665 I=1,NACT
IF(RPMS(I).LE.0.) GO TO 665
QP(I)=(RPMS(I)/C7(I))*((HPP-HSUMP)/(NSTAGE*RPMS(I)**2)
$-C8(I))
IF(QP(I).LE.0.) NACT2=NACT2-1
IF(QP(I).LE.0.) QP(I)=0.
665 CONTINUE
IF(NACT2.EQ.0) GO TO 654
IF(NACT2.EQ.NACT) GO TO 654
NACT=NACT2
GO TO 658
654
IF(VPP.LE.0.) VPP=0.
IF(VPP.LE.0.) HPP=-C1/C(1)
*
* CHECK TO SEE
IF Q/N ON PROPER SEGMENT
*
689 CONTINUE
NOGOOD=.FALSE.
DO 685 IZZ=1,NSTART
IF(RPMS(IZZ).LE.0.) GO TO 685
QPP=QP(IZZ)/RPMS(IZZ)
I=IZ(IZZ)
DO 686 J=1,5
IF(QPP.GT.QN(J).AND.QPP.LE.QN(J+1)) GO TO
688
686 CONTINUE
IF(QPP.GT.0..AND.I.EQ.5) GO TO 685
IF(QPP.GT.0..AND.I.LT.5) IZ(IZZ)=I+1
IF(QPP.LE.0..AND.I.EQ.1) GO TO 685
IF(QPP.LE.0..AND.I.GT.1) IZ(IZZ)=I-1
GO TO 684
688 IF(J.EQ.I)
GO TO 685
IZ(IZZ)=J
684
NOGOOD=.TRUE.
685 CONTINUE
ICHEK=ICHEK+1
IF(ICHEK.GT.20) WRITE(6,270)
IF(ICHEK.GT.20) STOP
IF(NOGOOD) GO TO 691
IF(HPP.GT.HSUMP) GO TO 699
CCC=-HSUMP-C1/C(1)
VPP=(-CCB+SQRT(CCB**2-4.*CCA*CCC))/(2.*CCA)
HPP=(VPP-C1)/C(1)
BYPASS=.TRUE.
GO TO 699
*
* RECOMPUTE
SOLUTION FOR REVISED POINTS ON CURVE
*
691 DO 690
J=1,NSTART
IF(RPMS(J).LE.0.) GO TO 690
QPP=QP(J)/RPMS(J)
I=IZ(J)
HA=HNSQ(I)
HB=HNSQ(I+1)
QA=QN(I)
QB=QN(I+1)
C7(J)=(HA-HB)/(QA-QB)
C8(J)=HB-QB*C7(J)
690 CONTINUE
GO TO 658
*
* SOLUTION IS
SUCCESSFUL. COMPUTE FINAL VALUES.
*
699
HNEW(1,1)=HPP
VNEWR(1,1)=VPP
IF(BYPASS) GO TO 696
IF(VPP.GT.0.) HPUMP=HPP-HSUMP
*
* ** CHECK
CONTINUITY AT PUMPS **
*
QPUMPS=0.
DO 695 I=1,NSTART
IF(RPMS(I).LE.0.) GO TO 695
QPUMPS=QPUMPS+QP(I)
695 CONTINUE
QLINE=VPP*AREA(1)
IF(ABS(QLINE).LE..0001.AND.ABS(QPUMPS).LE..0001) GO TO 696
IF(ABS((QPUMPS-QLINE)/QLINE).LT.0.01) GO
TO 696
WRITE(*,271) T
STOP
696 CONTINUE
*
* ** LOCATE MAXIMUM
AND MINIMUM HEADS AND H-VALUES **
*
531 DO 61
I=1,NPIPES
K=NPAR(I)+1
DO 61 J=1,K
IF(J.EQ.1) GO TO 532
DELTX=L(I)/(NPAR(I)*(ABS(VNEWL(I,J))+A(I)))
IF(DELTX.LT.DELTN) DELTN=DELTX
532 CONTINUE
IF(J.EQ.K) GO TO 533
DELTX=L(I)/(NPAR(I)*(ABS(VNEWR(I,J))+A(I)))
IF(DELTX.LT.DELTN) DELTN=DELTX
533 CONTINUE
IF(HNEW(I,J).LT.HLOW(I,J)) TLOW(I,J)=T
IF(HNEW(I,J).LT.HLOW(I,J))
HLOW(I,J)=HNEW(I,J)
IF(HNEW(I,J).GT.HHIGH(I,J)) THIGH(I,J)=T
IF(HNEW(I,J).GT.HHIGH(I,J))
HHIGH(I,J)=HNEW(I,J)
IF(CAV(I,J).GT.CAVMAX(I,J)) TCAV(I,J)=T
IF(CAV(I,J).GT.CAVMAX(I,J))
CAVMAX(I,J)=CAV(I,J)
HEAD(I,J)=HNEW(I,J)-(PIPEZ(I)+X(I,J)*L(I)*SINE(I))
IF(HEAD(I,J).LT.PHMAX) GO TO 66
IPMAX=I
XMAX=X(I,J)
PHMAX=HEAD(I,J)
TTMAX=T
GO TO 67
66 IF(HEAD(I,J).GT.PHMIN) GO TO 67
IPMIN=I
XMIN=X(I,J)
PHMIN=HEAD(I,J)
TTMIN=T
67 CONTINUE
61 CONTINUE
*
* ** WRITE OUT
H,V,AND HEAD VALUES AND TEST FOR TMAX **
*
IF(FAIL) GO TO 71
IF(T.LT.TIOU(III)) GO TO 70
IOUT=IOU(III)
III=III+1
70 CONTINUE
IF(MOD(II,IOUT).NE.0) GO TO 72
71 WRITE(6,205) T
DO 73 I=1,NPIPES
K=NPAR(I)+1
WRITE(6,206)
I,(X(I,J),HEAD(I,J),HNEW(I,J),VNEWL(I,J),
$VNEWR(I,J),CAV(I,J),J=1,K)
73 CONTINUE
IF(BYPASS.AND.VNEWR(1,1).GT.0.)
WRITE(6,215)
IF(BYPASS) GO TO 72
IF(RSTART) GO TO 75
QQQ=449.*Q
WRITE(6,212) RPM,QQQ,HPUMP
GO TO 72
75 DO 76 I=1,NSTART
IF(RPMS(I).LE.0.) GO TO 76
QQQ=449.*QP(I)
WRITE(6,212) RPMS(I),QQQ,HPUMP
76 CONTINUE
72 CONTINUE
IF(SAVE) CALL
GRAFCAV(II,NP,IX,INDEX,DELT,VNEWL,VNEWR,HEAD,CAV,
$NPAR,IENTRY)
IF(FAIL) GO TO 400
*
* ** PREPARE FOR
NEXT TIME STEP COMPUTATION **
*
DO 80 I=1,NPIPES
K=NPAR(I)+1
DO 80 J=1,K
VL(I,J)=VNEWL(I,J)
VR(I,J)=VNEWR(I,J)
80 H(I,J)=HNEW(I,J)
IF(T.GT.TMAX) GO TO 400
99 CONTINUE
*
*-------------------------------------------------------------------------------
* END OF TRANSIENT ANALYSIS
*-------------------------------------------------------------------------------
*
* ** CHECK
STABILITY CRITERIA **
*
400
IF(DELTN.GE.(0.9999*DELT)) GO TO 402
DELT=0.999*DELTN
WRITE(*,216) DELT
WRITE(6,216) DELT
GO TO 98
*
* ** WRITE OUT
MAXIMUM AND MINIMUM HEADS AND H-VALUES **
*
402 CONTINUE
WRITE(6,1011) CH12
WRITE(6,207)
DO 401 I=1,NPIPES
WRITE(6,208) I
K=NPAR(I)+1
DO 401 J=1,K
HEADMX=HHIGH(I,J)-(PIPEZ(I)+X(I,J)*L(I)*SINE(I))
HEADMN=
HLOW(I,J)-(PIPEZ(I)+X(I,J)*L(I)*SINE(I))
401 WRITE(6,209)
X(I,J),HEADMX,THIGH(I,J),HEADMN,TLOW(I,J),
$HHIGH(I,J),HLOW(I,J),CAVMAX(I,J),TCAV(I,J)
PMAX=PHMAX/2.31
PMIN=PHMIN/2.31
WRITE(6,210) PHMAX,PMAX,IPMAX,XMAX,TTMAX
WRITE(6,211) PHMIN,PMIN,IPMIN,XMIN,TTMIN
IF(PFILE) WRITE(6,1011) CH12
IF(PFILE) CALL
PROFILE(NPIPES,L,NPAR,PIPEZ,SINE,X,
$HLOW,HHIGH,HSTEAD)
IENTRY=2
IF(HVPRNT) WRITE(6,1011) CH12
IF(HVPRNT) CALL
GRAFCAV(II,NP,IX,INDEX,DELT,VNEWL,VNEWR,HEAD,CAV,
$NPAR,IENTRY)
IENTRY=3
WRITE(6,1011) CH,CH79
IF(PPLOT) WRITE(6,1011) CH12
IF(PPLOT) CALL
GRAFCAV(II,NP,IX,INDEX,DELT,VNEWL,VNEWR,HEAD,CAV,
$NPAR,IENTRY)
IENTRY=4
IF(GRAPH) WRITE(6,1011) CH12
IF(GRAPH) CALL
GRAFCAV(II,NP,IX,INDEX,DELT,VNEWL,VNEWR,HEAD,CAV,
$NPAR,IENTRY)
98 CONTINUE
*
* ******* FORMAT STATEMENTS
*******
*
100 FORMAT(A)
1011 FORMAT(3A)
101
FORMAT(5X,'NUMBER OF PIPES =',I3/5X,'MINIMUM NO. OF PARTS PIPE IS
$DIVIDED INTO =',I3//5X,'RESERVOIR
ELEVATION AT DOWNSTREAM END =',
$F7.1,' FT'/5X,'ELEVATION OF PUMP SUMP
=',F7.1,' FT'//
$5X,'ATMOSPHERIC PRESSURE HEAD =',F5.1,'
FT OF WATER'
$//5X,'ELEVATION OF DOWNSTREAM END OF LAST
PIPE =',F7.1,' FT'//5X,
$'MAXIMUM REAL TIME OF ANALYSIS =',F7.1,'
SEC'//
$5X,'ACCURACY OF ITERATION FOR STEADY
STATE DISCHARGE =',F5.1,
$' GPM'//5X,'PUMP BYPASS DIAMETER =
',F6.2,' IN'/5X,
$'PUMP BYPASS LOSS COEFFICIENT =',F6.2)
1010
FORMAT(/5X,'PRINT OUTPUT SCHEDULE IS IOUT =',I8,I5,4I8/
$5X,'CORRESPONDING INTERVALS (SEC)
=',F4.1,5F8.2)
102
FORMAT(//22X,'PIPE INPUT DATA'/22X,15('-')//5X,'PIPE DIAM-IN LENGT
$H-FT WAVE SPD-FPS PIPEZ-FT F-VALUE VEL-FPS'/
$5X,4('-'),1X,7('-'),1X,9('-'),1X,12('-'),1X,8('-'),1X,7('-'),
$1X,7('-'))
1022
FORMAT(//22X,'PIPE INPUT DATA'/22X,15('-')//5X,'PIPE DIAM-IN LENGT
$H-FT WAVE SPD-FPS PIPEZ-FT C-VALUE
VEL-FPS'/
$5X,4('-'),1X,7('-'),1X,9('-'),1X,12('-'),1X,8('-'),1X,7('-'),
$1X,7('-'))
1020
FORMAT(/5X,'PIPE DELT-SEC PARTS
SINE L/A-SEC INTERPOLATION',1X,
$'A-V
VALVE'/5X,4('-'),1X,8('-'),1X,5('-'),2X,6('-'),1X,7('-'),
$1X,13('-'),1X,9('-'))
103
FORMAT(5X,I3,2X,F6.2,1X,F9.1,6X,F6.0,4X,F6.0,3X,F5.4,3X,
$F5.2)
1033
FORMAT(5X,I3,2X,F6.2,1X,F9.1,6X,F6.0,4X,F6.0,3X,F5.0,3X,F5.2)
1030
FORMAT(5X,I3,2X,F6.3,3X,I3,2X,F7.5,1X,F6.2,5X,F5.3,9X,A)
1040
FORMAT(//5X,'DELTA T =',F9.5,' SEC'//)
104
FORMAT(///15X,'* PUMP DATA *'/)
105 FORMAT(5X,'NUMBER OF PUMPS IN PARALLEL
=',I3/5X,'EACH PUMP HAS',I2
$,' STAGES'/5X,'STEADY STATE PUMP SPEED
=',F7.1,' RPM'/5X,'EACH PUM
$P AND MOTOR UNIT HAS A MOMENT OF INERTIA
OF',F6.1,' LB-FT SQ'//5X
$'PUMP SPEED FOR DATA SHOWN BELOW
IS',F7.1,' RPM')
106
FORMAT(/8X,'Q-GPM H/STAGE-FT BHP/STAGE'/8X,5('-'),2X,10('-'),2X,
$9('-'))
107
FORMAT(6X,F7.1,2X,F8.1,4X,F8.2)
108
FORMAT(/3X,'PUMP NO.',I2,' STARTS AT',F8.2,' SEC AT',F6.0,
$' RPM AND RAMPS TO',F6.0,' RPM IN',F7.2,'
SEC')
1088 FORMAT(//5X,'* STEADY STATE LINE FLOW
RATE =',F9.1,
$' GPM AT A PUMP HEAD =',F7.1,' FT *'/)
201 FORMAT(//10X,A)
2011
FORMAT(10X,A//)
202 FORMAT(21X,14('*')/21X,'* INPUT DATA
*'/21X,14('*')/)
204 FORMAT(/'
PRESSURE HEADS, H-VALUES AND VELOCITIES AS FUNCTIONS
$ OF TIME'/2X,60('-'))
205
FORMAT(//20X,4X,' X HEAD-FT
H-FT VL-FPS VR-FPS CAV-FT3'/
$' TIME =',F9.3,' SEC',4X,'----- -------
----- ------ ------',
$1X,'-------')
206
FORMAT(/14X,'PIPE',I2,4X,F5.3,2F7.0,2X,2F7.2,F8.2/
$(24X,F5.3,2F7.0,2X,2F7.2,F8.2))
207
FORMAT(//29X,27('*')/29X,'* TABLE OF EXTREME VALUES *'/29X,
$27('*')//11X,'X MAX HEAD
TIME MIN HEAD TIME
'
$'MAX H
MIN H CAVMAX TIME'/9X,5('-'),2X,8('-'),2X,
$4('-'),3X,8('-'),2X,4('-'),3X,6('-'),3X,6('-'),3X,6('-'),
$2X,4('-'))
208
FORMAT(1X,'PIPE',I3)
209
FORMAT(9X,F5.3,2X,F7.1,1X,F6.1,3X,F7.1,1X,F6.1,2X,F7.1,2X,
$F7.1,2X,F6.1,1X,F6.1)
210 FORMAT(//'
MAXIMUM HEAD =',F6.1,' FT (',F5.1,' PSI) IN PIPE',
$I3,' AT X =',F5.3,' AT TIME =',F7.2,'
SEC')
211 FORMAT(/'
MINIMUM HEAD =',F6.1,' FT (',F5.1,' PSI) IN PIPE',
$I3,' AT X =',F5.3,' AT TIME =',F7.2,'
SEC')
212
FORMAT(/5X,'PUMP SPEED =',F7.1,' RPM',5X,'PUMP DISCHARGE =',F7.1
$,' GPM EACH',5X,'PUMP HEAD =',F6.1,' FT')
213 FORMAT(//2X,'THE
DIMENSIONED SIZE ',I4,' OF THE ARRAYS X(,), V(,),
$... IS TOO SMALL FOR THE NPARTS YOU
SELECTED.'/2X,'YOU MUST DECREA
$SE NPARTS OR INCREASE THE SIZE OF THE
ARRAYS TO AT LEAST ',I4)
214
FORMAT(////2X,'COLUMN SEPARATION HAS OCCURRED AT',F6.2,' SEC IN PI
$PE',I3,' AT LOCATION',F5.3)
215
FORMAT(/15X,'PUMP HEAD HAS DROPPED TO ZERO AND BYPASS IS OPEN')
216 FORMAT(/'
THE CHECK FOR STABILITY CRITERIA HAS FOUND A VIOLATION.'
$/' THE ANALYSIS MUST BE RERUN WITH A
SMALLER DELTA T.'/
$'
CHANGE DATA FILE TO SET DTNEW =',F7.4,' SEC & RERUN = "T"')
250
FORMAT(/////10X,'YOU HAVE EXCEEDED MAXIMUM DISCHARGE VALUE INPUT.
$EXECUTION IS TERMINATED.')
251
FORMAT(/////10X,'ITERATION TO COMPUTE STEADY STATE DISCHARGE WAS U
$NSUCCESSFUL AFTER 20 TRIES')
260
FORMAT(10X,'YOU HAVE EXCEEDED MAXIMUM BHP VALUE INPUT.')
270
FORMAT(/////10X,'THE SEARCH FOR A VALUE OF Q HAS CYCLED 20 TIMES W
$ITHOUT SUCCESS. EXECUTION IS
TERMINATED.')
271 FORMAT(//'
CONTINUITY CHECK AT THE PUMPS HAS FAILED AT T ='
$,F10.4,' SEC')
790 FORMAT(A)
791 FORMAT(/'
ENTER THE NAME OF YOUR INPUT DATA FILE: '\)
792 FORMAT(/'
ENTER THE NAME OF THE FILE ON WHICH YOU WANT YOUR OUTPUT
$ WRITTEN: '\)
END