SUBROUTINE INITAL C... C... MA1 - NH3 ABSORBER C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LS, KYA C... C... COMMON AREA TO LINK SUBROUTINES PARAMETER(NZ=21) COMMON/T/ T, NSTOP, NORUN + /Y/ XA(NZ), YA(NZ) + /F/ XAT(NZ), YAT(NZ) + /S/ XAZ(NZ), YAZ(NZ), YAS(NZ) + /C/ XA1, YA1, XA2, YA2, + ZL, GS, LS, KYA, S C... C... COEFFICIENTS IN POLYNOMIAL FOR EQUILIBRIUM LINE COMMON/E/A(2) C... C... INPUT/OUTPUT (I/O) UNIT NUMBERS COMMON/IO/ NI, NO C... C... GIVEN END COMPOSITIONS (MOLE FRACTIONS) C... C... ENTERING GAS YA1=0.0825D0 C... C... ENTERING LIQUID XA2=0.0D0 C... C... EXITING GAS YA2=0.003D0 C... C... CONVERT PRECEDING MOL FRACTIONS TO MOL RATIOS (BASED ON AIR AND C... WATER AS THE NONTRANSFERRED COMPONENTS) C... C... ENTERING GAS YA1=YA1/(1.0D0-YA1) C... C... ENTERING LIQUID XA2=XA2/(1.0D0-XA2) C... C... EXITING GAS YA2=YA2/(1.0D0-YA2) C... C... FLOW RATES (MOLS/M**2-MIN, SO THE COLUMN CROSS SECTIONAL AREA C... IS REQUIRED) C... C... PI PI=4.0D0*DATAN(1.0D0) C... C... COLUMN DIAMETER (M) D=0.154D0 C... C... COLUMN CROSS SECTIONAL AREA (M**2) S=PI*(D**2)/4.0D0 C... C... GAS VOLUMETRIC FLOW RATE (M**3/MIN) V=11.04D0 C... C... GAS CONSTANT (M**3-ATM/KG MOL-K) (NOTE - ONE KG MOL OF GAS C... OCCUPIES 22.4 M**3 AT 1 ATM AND 273.16 K) R=22.4D0*1.0D0/273.16D0 C... C... TEMPERATURE (K) TK=273.16D0+20.0D0 C... C... PRESSURE (ATM) P=1.0D0 C... C... GAS MOLAR FLOW RATE (KG MOL/MIN-M**2) G=V*P/(R*TK)*(1.0D0/S) C... C... GAS MOLAR FLOW RATE ON SOLUTE FREE BASIS (KG MOLS AIR C... /MIN-M**2) GS=G*(1.0D0/(1.0D0+YA1)) C... C... LIQUID MOLAR FLOW RATE ON SOLUTE FREE BASIS (KG MOLS H20 C... /MIN-M**2) LS=11.18D0/18.0D0*(1.0D0/S) C... C... EXITING LIQUID COMPOSITION XA1=OPLN(YA1,YA2,XA2,GS,LS) C... C... COLUMN LENGTH (M) ZL=3.28D0 C... C... WRITE SUMMARY OF OPERATING CONDITIONS WRITE(NO,1)YA1,XA2,YA2,XA1,LS,GS 1 FORMAT(' Summary of Operating Conditions',//, + ' Entering gas = ',F7.4,/, + ' Entering liquid = ',F7.4,/, + ' Exiting gas = ',F7.4,/, + ' Exiting liquid = ',F7.4,/, + ' Liquid flow rate = ',F7.4,' kg mols H20/min-m**2',/, + ' Gas flow rate = ',F7.4,' kg mols air/min-m**2',/) C... C... COEFFICIENTS IN THE LEAST SQUARES SECOND ORDER POLYNOMIAL C... APPROXIMATION OF THE EQUILIBRIUM DATA (THE COEFFICIENTS ARE C... WRITTEN INTO COMMON/E/ BY THIS CALL TO COEFF) CALL COEFF C... C... MASS TRANSFER COEFFICIENT (KG MOL AIR/M**2-MIN) KYA=73.8D0 C... C... COLUMN HEIGHT USING LM MOL RATIO CONCENTRATION DIFFERENCE WRITE(NO,2) 2 FORMAT(///,' Column height from log mean concentration', + ' differenece') Z=HTLM(GS,KYA,XA1,YA1,XA2,YA2) WRITE(NO,3)Z 3 FORMAT(/,' Column height = ',F6.2,' m',//) C... C... COLUMN HEIGHT USING INTEGRATION OF MOL RATIO CONCENTRATION C... DIFFERENCE WRITE(NO,4) 4 FORMAT(///,' Column height from numerical integration of', + ' concentration difference') Z=HTIN(GS,KYA,YA1,YA2) WRITE(NO,5)Z 5 FORMAT(/,' Column height = ',F6.2,' m',//) C... C... INITIAL CONDITIONS FOR PDE SOLUTION DO 10 I=1,NZ XA(I)=0.0D0 YA(I)=0.0D0 10 CONTINUE RETURN END SUBROUTINE DERV C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LS, KYA C... C... COMMON AREA TO LINK SUBROUTINES PARAMETER(NZ=21) COMMON/T/ T, NSTOP, NORUN + /Y/ XA(NZ), YA(NZ) + /F/ XAT(NZ), YAT(NZ) + /S/ XAZ(NZ), YAZ(NZ), YAS(NZ) + /C/ XA1, YA1, XA2, YA2, + ZL, GS, LS, KYA, S C... C... BOUNDARY CONDITION AT Z = 0 YA(1)=YA1 C... C... BOUNDARY CONDITION AT Z = ZL XA(NZ)=XA2 C... C... XAZ, YAZ CALL DSS020(0.0D0,ZL,NZ,XA,XAZ,-LS) CALL DSS020(0.0D0,ZL,NZ,YA,YAZ, GS) C... C... YA* GAS COMPOSITION (FROM EQUILIBRIUM LINE) DO 1 I=1,NZ YAS(I)=EQUIL(XA(I)) 1 CONTINUE C... C... PDES DO 2 I=1,NZ YAT(I)=-GS*YAZ(I)+KYA*(YAS(I)-YA(I)) XAT(I)= LS*XAZ(I)-KYA*(YAS(I)-YA(I)) 2 CONTINUE YAT(1)=0.0D0 XAT(NZ)=0.0D0 RETURN END SUBROUTINE PRINT(NI,NO) C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LS, KYA C... C... COMMON AREA TO LINK SUBROUTINES PARAMETER(NZ=21) COMMON/T/ T, NSTOP, NORUN + /Y/ XA(NZ), YA(NZ) + /F/ XAT(NZ), YAT(NZ) + /S/ XAZ(NZ), YAZ(NZ), YAS(NZ) + /C/ XA1, YA1, XA2, YA2, + ZL, GS, LS, KYA, S C... C... ARRAYS FOR PLOTTING PARAMETER(NP=21) DIMENSION TP(NP), XAP(NP), YAP(NP) C... C... MONITOR THE OUTPUT WRITE(*,*)T,YA(NZ) C... C... WRITE A HEADING FOR THE SOLUTION IF(T.LT.0.001D0)THEN WRITE(NO,1) 1 FORMAT(/,9X,'t',3X,'Xa(L,t)',3X,'Xa(0,t)', + 3X,'Ya(0,t)',3X,'Ya(L,t)',/) IP=0 END IF C... C... WRITE THE SOLUTION WRITE(NO,2)T,XA(NZ),XA(1),YA(1),YA(NZ) 2 FORMAT(F10.2,4F10.5) C... C... STORE THE SOLUTION FOR PLOTTING IP=IP+1 TP(IP)=T XAP(IP)=XA(1) YAP(IP)=YA(NZ) C... C... PLOT THE SOLUTION IF(IP.LT.NP)RETURN C... C... BOTTOM LIQUID COMPOSITION VS T CALL SPLOTS(1,NP,TP,XAP) WRITE(NO,3) 3 FORMAT(/,' Xa(0,t) vs t') C... C... TOP GAS COMPOSITION VS T CALL SPLOTS(1,NP,TP,YAP) WRITE(NO,4) 4 FORMAT(/,' Ya(L,t) vs t') RETURN END DOUBLE PRECISION FUNCTION OPLN(YA,YA2,XA2,GS,LS) C... C... Function OPLN computes the bulk liquid composition from C... the corresponding gas composition using the equation for C... the operating line based on the top compositions (at end C... 2) C... C... Xa = (Gs/Ls)*(Ya - Ya2) + Xa2 C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LS C... C... EQUATION FOR THE OPERATING LINE BASED ON THE COMPOSITIONS AT C... THE TOP OF THE COLUMN (END 2) XA=(GS/LS)*(YA-YA2)+XA2 C... C... EVALUATE FUNCTION OPLN OPLN=XA RETURN END DOUBLE PRECISION FUNCTION EQUIL(XA) C... C... Function EQUIL computes the gas phase composition in equilibroum C... with the liquid in the column, using the least squares second C... order polynomial evaluated by subroutine COEFF. C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... COEFFICIENTS IN POLYNOMIAL FOR EQUILIBRIUM LINE COMMON/E/A(2) C... C... SECOND ORDER LEAST SQUARES POLYNOMIAL EQUIL=A(1)*XA+A(2)*(XA**2) RETURN END DOUBLE PRECISION FUNCTION HTLM(GS,KYA,XA1,YA1,XA2,YA2) C... C... Function HTLM compuytes the height of a column using the log- C... mean concentration difference, according to the equation C... C... * C... z = (Gs/Kya)*(Ya1 - Ya2)/(Ya - Ya ) C... lm C... where C... * * C... * (Ya1 - Ya1 ) - (Ya2 - Ya2 ) C... (Ya - Ya ) = ---------------------------- C... lm * C... (Ya1 - Ya1 ) C... ln ------------ C... * C... (Ya2 - Ya2 ) C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION KYA C... C... COEFFICIENTS IN POLYNOMIAL FOR EQUILIBRIUM LINE COMMON/E/A(2) C... C... INPUT/OUTPUT (I/O) UNIT NUMBERS COMMON/IO/ NI, NO C... C... EQUILIBRIUM GAS PHASE COMPOSITION AT THE BOTTOM OF THE COLUMN YA1S=EQUIL(XA1) C... C... EQUILIBRIUM GAS PHASE COMPOSITION AT THE TOP OF THE COLUMN YA2S=EQUIL(XA2) C... C... LOG MEAN CONCENTRATION DIFFERENCE DELTA1=YA1-YA1S DELTA2=YA2-YA2S DYLM=(DELTA1-DELTA2)/DLOG(DELTA1/DELTA2) C... C... WRITE LOG MEAN CONCENTRATION DIFFERENCE WRITE(NO,1)DYLM 1 FORMAT(//,' Log-mean concentration difference = ',F9.6,/) C... C... WRITE NUMBER OF TRANSFER UNITS WRITE(NO,2)(YA1-YA2)/DYLM 2 FORMAT(/,' Number of transfer units = ',F6.2,/) C... C... WRITE HEIGHT OF TRANSFER UNIT (FT) WRITE(NO,3)GS/KYA 3 FORMAT(/,' Height of transfer unit = ',F6.3,' m',/) C... C... COLUMN HEIGHT HTLM=(GS/KYA)*(YA1-YA2)/DYLM RETURN END DOUBLE PRECISION FUNCTION HTIN(GS,KYA,YA1,YA2) C... C... Function HTIN computes the height of a column using the design C... integral C... C... Ya2 dYa C... z = (Gs/Kya) INT -------- C... Ya1 Ya - Ya* C... C... C... The integration is done by Simpson's rule, implemented in C... function SIMP C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LS, KYA C... COEFFICIENTS IN POLYNOMIAL FOR EQUILIBRIUM LINE COMMON/E/A(2) C... C... INPUT/OUTPUT (I/O) UNIT NUMBERS COMMON/IO/ NI, NO C... C... EXTERNAL FUNCTION FOR INTEGRAND EXTERNAL YAINT C... C... EVALUATE DESIGN INTEGRAL. NOTE THAT A SIGN REVERSAL IS C... REQUIRED SINCE DYA IN THE INTEGRAL IS NEGATIVE WHILE C... 1/(YA - YA*) IS POSITIVE NP=100 TU=-SIMP(YAINT,YA1,YA2,NP) C... C... WRITE NUMBER OF TRANSFER UNITS WRITE(NO,2)TU 2 FORMAT(//,' Number of transfer units = ',F6.2,/) C... C... WRITE HEIGHT OF TRANSFER UNIT (M) WRITE(NO,3)GS/KYA 3 FORMAT(/,' Height of transfer unit = ',F6.3,' m',/) C... C... COLUMN HEIGHT HTIN=(GS/KYA)*TU RETURN END DOUBLE PRECISION FUNCTION YAINT(YA) C... C... Function YAINT computes the integrand of the integral C... C... C... Ya2 dYa C... INT -------- C... Ya1 Ya - Ya* C... C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LS, KYA COMMON/C/ XA1, YA1, XA2, YA2, + ZL, GS, LS, KYA, S C... C... LIQUID COMPOSITION FROM THE OPERATING LINE XA=OPLN(YA,YA2,XA2,GS,LS) C... C... EQUILIBRIUM GAS COMPOSITION YAS=EQUIL(XA) C... C... INTEGRAND YAINT=1.0D0/(YA-YAS) RETURN END DOUBLE PRECISION FUNCTION SIMP(F,A,B,NP) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... EXTERNAL FUNCTION FOR INTEGRAND EXTERNAL F C... C... CHECK FOR AN EVEN NUMBER OF PANELS IF((NP/2*2).NE.NP)THEN WRITE(*,2)NP 2 FORMAT(//,' NP = ',I3,' (NOT EVEN)') SIMP=0.0D0 RETURN END IF C... C... INTEGRATION INTERVAL DX=(B-A)/DFLOAT(NP) C... C... INITIALIZE SUM SIMP=(F(A)+4.0D0*F(A+DX)+F(B)) C... C... CONTINUE SUM IF(NP.GT.2)THEN DO 1 I=2,NP-1,2 SIMP=SIMP + +2.0D0*F(A+DFLOAT(I )*DX) + +4.0D0*F(A+DFLOAT(I+1)*DX) 1 CONTINUE END IF C... C... INTEGRAL SIMP=(DX/3.0D0)*SIMP RETURN END SUBROUTINE COEFF C... C... Subroutine COEFF tests the least squares second order polynomial C... fit of the NH3-H20-air equilibrium data in Welty, et al, (1984), C... Example 2, pp 680-681. C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... COEFFICIENTS IN POLYNOMIAL FOR EQUILIBRIUM LINE COMMON/E/A(2) C... C... INPUT/OUTPUT (I/O) UNIT NUMBERS COMMON/IO/ NI, NO C... C... ARRAYS FOR FIVE DATA POINTS PARAMETER(N=5) DIMENSION X(N), Y(N), YC(N) C... C... EQUILIBRIUM DATA (FROM WELTY, ET AL) DATA X/0.0164D0, 0.0252D0, 0.0349D0, 0.0455D0, 0.0722D0/ DATA Y/0.021D0 , 0.032D0 , 0.042D0 , 0.053D0 , 0.08D0 / C... C... EVALUATE COEFFICIENTS IN THE LEAST SQUARES, SECOND ORDER C... POLYNOMIAL CALL POLY(X,Y,N) C... C... WRITE THE COEFFICIENTS WRITE(NO,1)A(1),A(2) 1 FORMAT(10X,' a1 = ',F8.4,' a2 = ',F8.4,/) C... C... HEADING FOR THE COMPUTED RESULTS WRITE(NO,2) 2 FORMAT(' i',6X,'XNH3',6X,'YNH3',6X,'YNH3',7X,'o/o',/, + 5X,7X,'tab', 7X,'tab', 6X,'calc',6X,'diff',/) C... C... COMPUTE VALUES OF THE GAS PHASE COMPOSITION FOR COMPARISON WITH C... THE TABULATED VALUES DO 5 I=1,5 C... C... GAS PHASE COMPOSITION FROM THE SECOND ORDER POLYNOMIAL YC(I)=A(1)*X(I)+A(2)*X(I)**2 C... C... COMPUTE AND PRINT THE GAS PHASE COMPOSITIONS DIFF=(YC(I)-Y(I))/Y(I)*100.0D0 WRITE(NO,4)I,X(I),Y(I),YC(I),DIFF 4 FORMAT(I5,3F10.4,F10.1) C... C... NEXT POINT 5 CONTINUE C... C... END OF CALCULATION RETURN END SUBROUTINE POLY(X,Y,N) C... C... Subroutine POLY computes the coefficients in the least squares C... second order polynomial approximating the equilibrium data. C... C... Consider the least squares analysis of data using a polynomial C... model C... C... y = a0 + a1*x + a2*(x**2) + ... (1) C... C... For purposes of discussion, we consider a second order poly- C... nomial. The extension to higher order polynomials is straight- C... forward. C... C... Application of the least squares principle to N data gives C... C... N C... min E(a0,a1,a2) = sum (y(i) - ym(i))**2 C... i=1 C... C... N C... = sum (y(i) - (a0 + a1*x(i) + a2*(x(i)**2)))**2 (2) C... i=1 C... C... We can find the minimum of eq. (2), i.e., the values of a0, a1 C... and a2 which minimize E, by the usual methods of differential C... calculus. C... C... Before proceeding, however, we consider the special case for C... which the model of eq. (1) is required to pass through the C... point y(0) = 0, x(0) = 0. Thus we see a0 = 0 in eq. (1), and C... we therefore have to evaluate only a1 and a2 in eq. (2). C... C... The required partial derivatives (denoted with a "d" because C... of the limitations of the character set) are set to zero C... C... dE/da1 = C... C... N C... 2 sum (y(i) - (a1*x(i) + a2*(x(i)**2)))*(-x(i)) = 0 (3) C... i=1 C... C... C... dE/da2 = C... C... N C... 2 sum (y(i) - (a1*x(i) + a2*(x(i)**2)))*(-x(i)**2) = 0 (4) C... i=1 C... C... Eqs. (3) and (4) can be rearranged to C... C... N N N C... a1 sum x(i)**2 + a2 sum x(i)**3 = sum x(i)*y(i) (5) C... i=1 i=1 i=1 C... C... N N N C... a1 sum x(i)**3 + a2 sum x(i)**4 = sum (x(i)**2)*y(i) (6) C... i=1 i=1 i=1 C... C... Eqs. (5) and (6) are two linear algebraic equations for a1 and C... a2. The following code computes a1 and a2 from a vector of data C... of length N. C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... COEFFICIENTS IN POLYNOMIAL FOR EQUILIBRIUM LINE COMMON/E/A(2) C... C... INPUT/OUTPUT (I/O) UNIT NUMBERS COMMON/IO/ NI, NO C... C... ARRAYS FOR X, Y DIMENSION X(N), Y(N) C... C... COMPUTE THE ELEMENTS OF THE COEFFICIENT MATRIX AND RHS VECTOR OF C... EQS. (5) AND (6) A11=0.0D0 A12=0.0D0 A21=0.0D0 A22=0.0D0 B1=0.0D0 B2=0.0D0 DO 1 I=1,N A11=A11+X(I)**2 A12=A12+X(I)**3 A21=A12 A22=A22+X(I)**4 B1=B1+X(I)*Y(I) B2=B2+(X(I)**2)*Y(I) 1 CONTINUE C... C... COMPUTE A1 AND A2 IN EQ. (1) (BY DETERMINANTS) DET=A11*A22-A12*A21 C... C... TEST IF THE ALGEBRAIC EQUATIONS ARE NUMERICALLY SINGULAR IF(DABS(DET).LT.1.0D-20)THEN C... C... DETERMINANT IS SMALL SO TERMINATE THE CALCULATION. THE C... PRECEDING THRESHOLD OF 1.0D-20 WAS SELECTED ARBITRARILY C... AND IS DETERMINED PRIMARILY BY THE WORDLENGTH OF THE C... COMPUTER WRITE(NO,2)DET 2 FORMAT(' Determinant of coefficient matrix = ',D11.3) STOP END IF C... C... SYSTEM IS NOT NUMERICALLY SINGULAR (ILL-CONDITIONED), SO CONTINUE C... THE SOLUTION OF THE ALGEBRAIC SYSTEM DET1=B1*A22-B2*A12 DET2=B2*A11-B1*A21 C... C... COMPUTE COEFFICIENTS A1 AND A2 FROM THE DETERMINANTS (BY KRAMER'S C... RULE, WHICH IS GENERALLY NOT RECOMMENDED FOR THE SOLUTION OF C... SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS, BUT IS ADEQUATE FOR THIS C... SMALL 2 X 2 PROBLEM) A(1)=DET1/DET A(2)=DET2/DET RETURN END