C**************************************************************** ! This component has been assigned Type Number 204. If that number conflicts with ! another user Type number, you will need to change it and recompile the appropriate ! dll. SUBROUTINE TYPE204(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) !DEC$ATTRIBUTES DLLEXPORT :: TYPE204 C Version from: 12/14/89 C SAMPLE INPUT AND OUTPUT ARE LOCATED AFTER THE NOMENCLATURE C---------------------------------------------------------------- C**** This Type represents a centrifugal pump, which is designed C to run with an electric motor. C The model requires two sets of data: head-flowrate data at C reference speed and efficiency-flowrate data at reference C speed. A linear regression is performed to fit a curve C through those data. The normal equations are solved using C Gaussian elimination with partial pivoting. Newtons method C is implemented to solve a system of nonlinear equations. C**************************************************************** C**************************************************************** C**** Variables: C TL -- Torque [Nm] C N -- Speed [1/s] C Q -- Volumetric flowrate [m^3/s] C H -- total dynamic head [m] C ETA -- pump efficiency C TOFF,QOFF,NOFF,ETAOFF -- T,Q,N,ETA at shutoff condition C TREQ -- required torque to provide flow C NREF -- Reference speed C NINIT -- initial speed for Newton's method C NMAX -- maximum speed of pump C QINIT -- initial flowrate for Newton's method C QNOM -- nominal design flowrate [m^3/hour] C MAXIT -- maximum number of iterations to find the C speed and flowrate in Newton's method C COUNTER -- counting the number of iterations C POFF -- Shaftpower at shutoff condition at reference speed C POFF is input as parameter. If no information is C available, Poff should be set equal some negative number C EPS -- error tolerance in Newton's method C A1-H1 -- coefficients for the H=f(Q,N) and ETA=f(Q,N) relations C S1-S8 -- defined constants in the equations C RHOG -- density, rho, times acceleration const., g C CONST -- specific weight of water(rho*g) devided by two times pi C X -- two dimensional array contains N and Q -- old values C X1=N X2=Q C XNEW -- same as X, but containing new values C F -- two dimensional array contains values of the objectiv function C J -- Jacobian matrix (2*2) contains the derivatives C JINV -- inverse Jacobian matrix C DET -- determinant of the Jacobian matrix C LOOP -- control variable C PHYD -- water power C PMECH -- mechanical (input) power C FAIL -- error message C :0 - everything is okay C :1 - maximum pump speed is exceeded C :2 - newton method could not find proper solution C :3 - newton method did not converge within C specified number of iterations C :4 - chosen pump is to small to match system,i.e. C speed required to provide flow is greater C than maximum pump speed C C**************************************************************** C C SAMPLE INPUT AND OUTPUT C C INPUT SAMPLE VALUES C XIN(1) TL 9.25 C C PARAMETERS C PAR(1) S3 30.48 C PAR(2) S4 421040.0 C PAR(3) NREF 50.0 C PAR(4) NMAX 58.3 C PAR(5) POFF -1.0 C PAR(6) QNOM 15.0 C C OUTPUTS C OUT(1) N 58.19 C OUT(2) TL 9.25 C OUT(3) Q .004918 C OUT(4) H 40.66 C OUT(5) ETA .5806 C OUT(6) PHYD 1964.0 C OUT(7) PMECH 3382.0 C OUT(8) FAIL 0.0 C C****************************************************************** IMPLICIT NONE DOUBLE PRECISION XIN, OUT INTEGER INFO,COUNTER,MAXIT,LOOP,FAIL INTEGER NSTORE,IAV,ISTORE REAL TIME,T,DTDT,ICNTRL,PAR,S REAL TL,N,Q,H,ETA,EPS,RHOG,PHYD,PMECH,QNOM,HOFF REAL NOFF,QOFF,ETAOFF,NREF,POFF,TOFF,TREQ REAL X(2),XNEW(2),J(2,2),JINV(2,2),F(2) REAL S1,S2,S3,S4,S5,S6,S7,S8 REAL A1,B1,C1,E1,F1,G1,H1,PI,CONST REAL DET,DELTA1,DELTA2,NINIT,QINIT,NMAX DIMENSION XIN(1), OUT(8), PAR(6), INFO(15) COMMON /STORE/ NSTORE,IAV,S(5000) COMMON /COEF/ S1,S2,S3,S4,S5,S6,S7,S8,CONST ! Set the version information for TRNSYS IF (INFO(7).EQ.-2) THEN INFO(12) = 15 RETURN 1 ENDIF C---------------------------------------------------------------- C**** First call of component IF (INFO(7).EQ.-1) THEN INFO(6)=8 INFO(9)=0 C**** storage allocation: speed and flowrate from previous C**** iteration are stored in the S-array C**** S(ISTORE) -- NOLD C**** S(ISTORE+1) -- QOLD INFO(10)=2 CALL TYPECK(1,INFO,1,6,0) ISTORE=INFO(10) EPS=1.e-4 MAXIT=100 C**** Set parameters C**** System curve S3=PAR(1) S4=PAR(2) NREF=PAR(3) NMAX=PAR(4) POFF=PAR(5) QNOM=PAR(6) C**** Coefficients for the equations C**** Pump characteristics (obtained from regression) CALL HEADPOLYNOMIAL (A1,B1) C**** Efficiency characteristics CALL EFFPOLYNOMIAL(E1,F1,G1,H1) C**** Const=specific weight of water(rho*g) divided by 2*pi PI=ACOS(-1.) RHOG=9819.8 CONST=RHOG/2./PI C**** Definition of some constant values S1=A1/NREF**2 S2=B1 S5=E1 S6=F1*NREF S7=G1*NREF**2 S8=H1*NREF**3 C**** Initial guess for speed (X1) and flowrate (X2) C**** the initial speed is taken to be the speed where C**** HEADpump=HEADsystem at zero flowrate NINIT=NREF*(S3/A1)**.5 IF(NINIT.GE.NMAX)THEN FAIL=4 ENDIF C**** the initial flowrate is taken to be the flowrate C**** at maximum permitted speed where HEADpump=HEADsystem QINIT=((S3-S1*NMAX**2)/(S2-S4))**.5 X(1)=NINIT X(2)=QINIT C**** The following part evaluates the torque required to C**** provide flow. C**** Evaluation of Q and H at shutoff condition C**** QOFF is choosen to be very small, i.e. QNOM per second C**** divided by 100 QOFF=QNOM/3.6E5 HOFF=A1+B1*QOFF**2 C**** Computation of torque at shutoff condition at reference C**** speed IF(POFF.LT.0)THEN ETAOFF=E1+F1*QOFF+G1*QOFF**2+H1*QOFF**3 TOFF=CONST*QOFF*HOFF/ETAOFF/NREF ELSE TOFF=POFF/(2.*PI*NREF) ENDIF TREQ=TOFF*(NINIT/NREF)**2 FAIL=0 C---------------------------------------------------------------- ELSE C**** First and following calls in timestep FAIL=0 C**** Set inputs TL=XIN(1) X(1)=S(ISTORE) X(2)=S(ISTORE+1) C**** to avoid division by zero, i.e. if torque=0 IF(TL.GT.TREQ)THEN LOOP=0 C**************************************************** 50 CONTINUE COUNTER=1 C**** Iterative process to find the solution of the system of C**** nonlinear equations, using Newton's method 100 CONTINUE C**** set up the objective functions CALL SYSTEM(F,X,TL) C**** set up of the Jacobian matrix CALL JACOBIAN(J,X,TL) C**** To inverse the matrix a direct method via the adjuncte is used C**** Evaluation of the determinant of the Jacobian DET=J(1,1)*J(2,2)-J(2,1)*J(1,2) C**** Evaluation of the inverse of the Jacobian JINV(1,1)=J(2,2)/DET JINV(1,2)=-J(1,2)/DET JINV(2,1)=-J(2,1)/DET JINV(2,2)=J(1,1)/DET C**** Basic algorithm XNEW(1)=X(1)-(JINV(1,1)*F(1)+JINV(1,2)*F(2)) XNEW(2)=X(2)-(JINV(2,1)*F(1)+JINV(2,2)*F(2)) C**** Check on the infinity norm of the solution vector C**** i.e. the maximum value of xnew-xold (=max. error) C**** about convergence DELTA1=ABS(XNEW(1)-X(1)) DELTA2=ABS(XNEW(2)-X(2)) IF (MAX(DELTA1,DELTA2).GT.EPS) THEN IF (COUNTER.LT.MAXIT) THEN COUNTER=COUNTER+1 X(1)=XNEW(1) X(2)=XNEW(2) GOTO 100 ELSE C**** Convergence was not attained after maxit iterations FAIL=3 ENDIF ENDIF C**************************************************** X(1)=XNEW(1) X(2)=XNEW(2) N=XNEW(1) Q=XNEW(2) C**** Check if Q is in the negative range, i.e. the wrong solution C**** is found IF ((Q.LT.0.).OR.(N.LT.0)) THEN C**** to protect an endless loop situation the following C**** control function is introduced IF(LOOP.EQ.0)THEN X(1)=NMAX X(2)=5.*QINIT LOOP=1 GOTO 50 ELSEIF(LOOP.EQ.1)THEN C**** try again with higher limits X(1)=2.*NMAX X(2)=2.*QNOM LOOP=2 GOTO 50 ELSE C**** Newton could not find proper soluion C**** go ahead FAIL=2 ENDIF ENDIF C**** Evaluation of total head and efficiency H=S3+S4*Q**2 ETA=CONST*Q*H/TL/N C**** Check if speed exceeds maximum speed IF(N.GT.NMAX)THEN C**** Flowrate is taken to be the flowrate at maximum C**** permitted speed Q=QINIT FAIL=1 C**** Evaluation of total head H=S3+S4*Q**2 C**** Efficiency is evaluated at NMAX and QINIT ETA=S5+S6*Q/NMAX+S7*Q**2/NMAX**2+S8*Q**3/NMAX**3 C**** Correction of inital guess X(1)=NMAX X(2)=QINIT ENDIF ELSE C**** In this case the affinity laws are valid. The torque C**** or the speed follows a parabolic function N=NREF*(TL/TOFF)**.5 Q=0. H=HOFF*(N/NREF)**2 ETA=0. ENDIF ENDIF C---------------------------------------------------------------- IF((X(1).GT.0).AND.(X(2).GT.0.))THEN S(ISTORE)=X(1) S(ISTORE+1)=X(2) ENDIF PHYD=RHOG*Q*H PMECH=2*PI*TL*N C**** Set outputs OUT(1)=N OUT(2)=TL OUT(3)=Q OUT(4)=H OUT(5)=ETA OUT(6)=PHYD OUT(7)=PMECH OUT(8)=FAIL RETURN 1 END C**************************************************************** SUBROUTINE SYSTEM(F,X,TL) C**** Subroutine provides the objective function of the system C of nonlinear equations C TL,X1,X2 in C F out REAL F(2),X(2) REAL TL,S1,S2,S3,S4,S5,S6,S7,S8,CONST COMMON /COEF/ S1,S2,S3,S4,S5,S6,S7,S8,CONST F(1)=S1*X(1)**2+S2*X(2)**2-S3-S4*X(2)**2 F(2)=CONST/TL*(S1*X(1)**4*X(2)+S2*X(1)**2*X(2)**3) & -S5*X(1)**3-S6*X(1)**2*X(2)-S7*X(1)*X(2)**2-S8* & X(2)**3 END C********************************************************************* SUBROUTINE JACOBIAN(J,X,TL) C**** Subroutine provides all the derivatives for the Jacobian C matrix C X1,X2,TL in C J out REAL J(2,2),X(2) REAL TL,S1,S2,S3,S4,S5,S6,S7,S8,CONST COMMON /COEF/ S1,S2,S3,S4,S5,S6,S7,S8,CONST J(1,1)=2.*S1*X(1) J(1,2)=2.*S2*X(2)-2.*S4*X(2) J(2,1)=CONST/TL*(4.*S1*X(1)**3*X(2)+2.*S2*X(1)*X(2)**3) & -3.*S5*X(1)**2-2.*S6*X(1)*X(2)-S7*X(2)**2 J(2,2)=CONST/TL*(S1*X(1)**4+3.*S2*X(1)**2*X(2)**2)- & S6*X(1)**2-2.*S7*X(1)*X(2)-3.*S8*X(2)**2 END C************************************************************** SUBROUTINE HEADPOLYNOMIAL (A1,B1) C**** Program fits a polynomial through a given number of data C**** points using least squares method C**** the polynmial has the form y=a+bx^2 C**** HEAD -- head of the pump [m] C**** CAP -- Capacity (volumetric flowrate) [m^3/sec] C**** A1,B1 -- Coefficients of thr constructed polynomial C**** MEAN -- Average over the head values C**** RSQUARES -- R^2 measure of curve fit C**** SUM1-SUM6 -- Summation terms occuring in the normal equations C**** Y -- Array containing the fitted values REAL SUM1,SUM2,SUM3,SUM4,SUM5,SUM6,A1,B1 REAL MEAN,RSQUARE REAL HEAD(50),CAP(50),Y(50) INTEGER M C**** Input file contains the data OPEN (12,FILE='HEAD.DAT',STATUS='OLD') REWIND(12) C**** Read in number of data pairs READ(12,*) M C**** Read in the data pairs DO 100 I=1,M READ(12,*) CAP(I),HEAD(I) 100 CONTINUE C**** Evaluation of the sum terms in the normal equations SUM1=0. SUM2=0. SUM3=0. SUM4=0. DO 200 J=1,M SUM1=SUM1+CAP(J)**2 SUM2=SUM2+HEAD(J) SUM3=SUM3+CAP(J)**4 SUM4=SUM4+HEAD(J)*CAP(J)**2 200 CONTINUE C**** Evaluation of the coefficients of the polynomial B1=(M*SUM4-SUM2*SUM1)/(M*SUM3-SUM1**2) A1=(SUM2-B1*SUM1)/M DO 300 L=1,M Y(L)=A1+B1*CAP(L)**2 300 CONTINUE C**** Estimation of the curve fit precission MEAN=SUM2/M SUM5=0. SUM6=0. DO 400 LL=1,M SUM5=SUM5+(Y(LL)-MEAN)**2 SUM6=SUM6+(HEAD(LL)-MEAN)**2 400 CONTINUE C**** R^2 measures the proportion of total variation about the mean C**** explained by the regression RSQUARE=SUM5/SUM6 END C*************************************************************** SUBROUTINE EFFPOLYNOMIAL(E1,F1,G1,H1) C**** Subroutine fits a polynomial through a given number of data C**** points using least squares method C**** the polynmial has the form y=a+bx+cx^2+dx^3 C**** ETA -- pump efficiency C**** CAP -- Capacity (volumetric flowrate) C**** ESUM1-ESUM12 -- Summation terms occuring in the normal equations C**** E1,F1,G1,H1 -- Coefficients of the constructed polynomial C**** MEAN -- Average over the efficiency values C**** RSQUARES -- R^2 measure of curve fit C**** Y -- Array of dimension NN containing the fitted values C**** A -- Array of dimension NN*NN containing the coefficients C**** of the normal equations C**** B -- Array of dimension NN (B-vector) C**** SOL -- Array containing the solution, i.e. the coefficients C**** E1,F1,G1,H1 C**** LINDEX -- index vector (Gaussian elimination) C**** S -- scale vector (Gaussian elimination) INTEGER M,MM,NN PARAMETER (NN=4) REAL ESUM1,ESUM2,ESUM3,ESUM4,ESUM5,ESUM6,ESUM7,ESUM8 REAL ESUM9,SUM10,ESUM11,ESUM12,E1,F1,G1,H1 REAL MEAN,RSQUARE REAL ETA(50),CAP(50),Y(50),A(NN,NN),B(NN),SOL(NN) REAL LINDEX(NN),S(NN) C**** Input file contains the data OPEN (14,FILE='EFF.DAT',STATUS='OLD') REWIND(14) C**** Read in number of data pairs READ(14,*) M MM=M+4 C**** Read in the data pairs DO 100 I=1,M READ(14,*) CAP(I),ETA(I) 100 CONTINUE C**** To force the fitted curve through the origin a couple C**** of data pairs are set equal 0 DO 120 II=M+1,MM CAP(II)=0. ETA(II)=0. 120 CONTINUE C**** Evaluation of the sum terms in the normal equations ESUM1=0. ESUM2=0. ESUM3=0. ESUM4=0. ESUM5=0. ESUM6=0. ESUM7=0. ESUM8=0. ESUM9=0. ESUM10=0. DO 200 J=1,M ESUM1=ESUM1+CAP(J) ESUM2=ESUM2+CAP(J)**2 ESUM3=ESUM3+CAP(J)**3 ESUM4=ESUM4+CAP(J)**4 ESUM5=ESUM5+CAP(J)**5 ESUM6=ESUM6+CAP(J)**6 ESUM7=ESUM7+ETA(J) ESUM8=ESUM8+ETA(J)*CAP(J) ESUM9=ESUM9+ETA(J)*CAP(J)**2 ESUM10=ESUM10+ETA(J)*CAP(J)**3 200 CONTINUE C**** Set up the matrix containing the normal equations A(n,n) C**** In this case it is a (4*4) matrix A(1,1)=REAL(MM) A(2,1)=ESUM1 A(2,2)=ESUM2 A(3,1)=ESUM2 A(3,2)=ESUM3 A(3,3)=ESUM4 A(4,1)=ESUM3 A(4,2)=ESUM4 A(4,3)=ESUM5 A(4,4)=ESUM6 C**** because of the symmetry of the matrix,it can be C**** filled up by DO 250 I=1,4 DO 250 J=1,4 IF(I.NE.J)THEN A(I,J)=A(J,I) ENDIF 250 CONTINUE C**** Set up B vector (dimension=4) B(1)=ESUM7 B(2)=ESUM8 B(3)=ESUM9 B(4)=ESUM10 C**** The linear system of equations (normal equations) is solved C**** using Gaussian elimination with partial pivoting C**** LINDEX & S must occur in the calling sequence because they C**** have a variable dimension. They are not used as output CALL GAUSS(A,B,SOL,LINDEX,S,NN) C**** Output E1=SOL(1) F1=SOL(2) G1=SOL(3) H1=SOL(4) DO 300 L=1,MM Y(L)=E1+F1*CAP(L)+G1*CAP(L)**2+H1*CAP(L)**3 300 CONTINUE C**** Estimation of the curve fit precission MEAN=ESUM7/REAL(MM) ESUM11=0. ESUM12=0. DO 400 LL=1,MM ESUM11=ESUM11+(Y(LL)-MEAN)**2 ESUM12=ESUM12+(ETA(LL)-MEAN)**2 400 CONTINUE C**** R^2 measures the proportion of total variation about the mean C**** explained by the regression RSQUARE=ESUM11/ESUM12 END C******************************************************************* C**** Subroutine Gauss with partial pivoting SUBROUTINE GAUSS(A,B,X,L,S,NN) REAL A(NN,NN),B(NN),L(NN),S(NN),X(NN) C**** Scale factors(s) and pivot(l) DO 30 I=1,NN L(I)=I SMAX=0. DO 20 J=1,NN SMAX=AMAX1(SMAX,ABS(A(I,J))) 20 CONTINUE S(I)=SMAX 30 CONTINUE C**** Look for largest ratio to select pivot row DO 70 K=1,NN-1 RMAX=0. DO 40 I=K,NN R=ABS(A(L(I),K))/S(L(I)) IF (R.GT.RMAX) THEN J=I RMAX=R ENDIF 40 CONTINUE C**** Update pivot vector XK=L(J) L(J)=L(K) L(K)=XK C**** Forward elimination DO 60 I=K+1,NN XMULT=A(L(I),K)/A(L(K),K) DO 50 J=K+1,NN A(L(I),J)=A(L(I),J)-XMULT*A(L(K),J) 50 CONTINUE A(L(I),K)=XMULT 60 CONTINUE 70 CONTINUE C**** Backward substitution DO 100 K=1,NN-1 DO 100 I=K+1,NN B(L(I))=B(L(I))-A(L(I),K)*B(L(K)) 100 CONTINUE X(NN)=B(L(NN))/A(L(NN),NN) DO 150 I=NN-1,1,-1 SUM=B(L(I)) DO 140 J=I+1,NN SUM=SUM-A(L(I),J)*X(J) 140 CONTINUE X(I)=SUM/A(L(I),I) 150 CONTINUE END