２成分の蒸留計算に比べて、３成分以上の蒸留計算はあまり取り扱われませんが、実際の蒸留塔においては多成分の蒸留がほとんどを占めています。
このページでは、FORTRANによるプログラムを利用して計算した多成分の蒸留例を紹介します。
目次
3成分の蒸留計算
このページでは、以下のプログラム（FORTRAN）を用いて蒸留計算を行いました。
PROGRAM MAIN
IMPLICIT NONE
INTEGER I,J,K,L,II,JJ,KK,LL
INTEGER IW,ITE,INFO
INTEGER NCOMP,NPLATE,NPLATE_FEED,N_JACOB
CHARACTER(16) input_data
CHARACTER(18) CHARA(10)
REAL(8) P0
REAL(8) FEED_STREAM,TOP,REF,BTM
REAL(8) delta
REAL(8) ErrorTolerance
REAL(8) allError
INTEGER,ALLOCATABLE,DIMENSION(:) :: IPIV
REAL(8),ALLOCATABLE,DIMENSION(:) :: FEED,LIQ,VAP
REAL(8),ALLOCATABLE,DIMENSION(:) :: D_LIQ,D_VAP
REAL(8),ALLOCATABLE,DIMENSION(:) :: LIQ_PRIME,VAP_PRIME
REAL(8),ALLOCATABLE,DIMENSION(:) :: TEMP,P,D_TEMP
REAL(8),ALLOCATABLE,DIMENSION(:) :: xLIQ_FEED
REAL(8),ALLOCATABLE,DIMENSION(:,:) :: xLIQ,D_xLIQ
REAL(8),ALLOCATABLE,DIMENSION(:,:) :: ALPHA,GAM,D_ALPHA
REAL(8),ALLOCATABLE,DIMENSION(:,:) :: ANTOINE,WILLSON,CP
REAL(8),ALLOCATABLE,DIMENSION(:) :: TB
REAL(8),ALLOCATABLE,DIMENSION(:,:) :: JACOB
REAL(8),ALLOCATABLE,DIMENSION(:) :: FX
REAL(8),ALLOCATABLE,DIMENSION(:) :: E,D_E
REAL(8),ALLOCATABLE,DIMENSION(:) :: H_LIQ,H_VAP,H_FEED
REAL(8),ALLOCATABLE,DIMENSION(:) :: D_H_LIQ,D_H_VAP
REAL(8),ALLOCATABLE,DIMENSION(:,:) :: M_LIQ,M_VAP,M_FEED
REAL(8),ALLOCATABLE,DIMENSION(:,:) :: D_M_LIQ,D_M_VAP
REAL(8),ALLOCATABLE,DIMENSION(:) :: M_TOP
REAL(8),ALLOCATABLE,DIMENSION(:,:) :: M,Q,D_M,D_Q
REAL(8),ALLOCATABLE,DIMENSION(:) :: working
REAL(8) COMP_FEED(10)
INTEGER MATERIAL(10)
IW = 6
P0 = 1013.25d0
delta = 0.010d0
ErrorTolerance = 0.10d0**(5.0d0)
CALL GETARG(1,input_data)
CALL READ_INPUT(input_data,NCOMP,NPLATE,NPLATE_FEED,
1 FEED_STREAM,TOP,REF,BTM,MATERIAL,COMP_FEED,CHARA)
ALLOCATE (IPIV((2*NCOMP+1)*NPLATE))
ALLOCATE (FEED(NPLATE),LIQ(NPLATE),VAP(NPLATE))
ALLOCATE (D_LIQ(NPLATE),D_VAP(NPLATE))
ALLOCATE (ALPHA(NCOMP,NPLATE),D_ALPHA(NCOMP,NPLATE))
ALLOCATE (GAM(NCOMP,NPLATE))
ALLOCATE (TEMP(NPLATE),D_TEMP(NPLATE),P(NPLATE))
ALLOCATE (JACOB((2*NCOMP+1)*NPLATE,(2*NCOMP+1)*NPLATE))
ALLOCATE (FX((2*NCOMP+1)*NPLATE))
ALLOCATE (xLIQ(NCOMP,NPLATE),D_xLIQ(NCOMP,NPLATE))
ALLOCATE (xLIQ_FEED(NCOMP))
ALLOCATE (ANTOINE(3,NCOMP),WILLSON(NCOMP,NCOMP))
ALLOCATE (CP(3,NCOMP),TB(NCOMP))
ALLOCATE (H_LIQ(NPLATE),H_VAP(NPLATE),H_FEED(NPLATE))
ALLOCATE (D_H_LIQ(NPLATE),D_H_VAP(NPLATE))
ALLOCATE (M_LIQ(NCOMP,NPLATE),M_VAP(NCOMP,NPLATE))
ALLOCATE (D_M_LIQ(NCOMP,NPLATE),D_M_VAP(NCOMP,NPLATE))
ALLOCATE (M_FEED(NCOMP,NPLATE),M_TOP(NCOMP))
ALLOCATE (E(NPLATE),D_E(NPLATE))
ALLOCATE (M(NCOMP,NPLATE),Q(NCOMP,NPLATE))
ALLOCATE (D_M(NCOMP,NPLATE),D_Q(NCOMP,NPLATE))
ALLOCATE (working(NPLATE))
CLOSE(1)
WRITE(IW,*) ""
WRITE(IW,'(a10,10a10)') "",("----------",I=1,7)
WRITE(IW,'(15x,a35)') "CALCULATION PROPERTY"
WRITE(IW,'(a10,10a10)') "",("----------",I=1,7)
WRITE(IW,*) ""
WRITE(IW,'(a20)') " "
WRITE(IW,'(14x,a16)') "COLUMN PROPERTY"
WRITE(IW,'(23x,a20,i10)') " NPLATE = ",NPLATE
WRITE(IW,'(23x,a20,i10)') " FEED_PLATE = ",NPLATE_FEED
WRITE(IW,'(a20)') " "
WRITE(IW,'(13x,a25)') "STREAM PROPERTY(mol/hr)"
WRITE(IW,'(23x,a20,f10.2)') "FEED = ",FEED_STREAM
WRITE(IW,'(23x,a20,f10.2)') "TOP = ",TOP
WRITE(IW,'(23x,a20,f10.2)') "REFLUX = ",REF*TOP
WRITE(IW,'(a20)') " "
WRITE(IW,'(12x,a16,18x,a12)') "FEED MATERIAL","(mol ratio)"
DO I = 1, NCOMP
WRITE(IW,'(29x,a10,a4,f10.5)')
1 CHARA(MATERIAL(I))," = ",COMP_FEED(I)
ENDDO
CALL SET_PROPERTY(NCOMP,MATERIAL,ANTOINE,WILLSON,CP,TB)
CALL Calc_PXY(ANTOINE,WILLSON,P0,NComp)
CALL SET_xLIQ_START(NCOMP,NPLATE,COMP_FEED,xLIQ_FEED,xLIQ,
1 TB,NPLATE_FEED,FEED_STREAM,TOP,BTM)
CALL SET_INITIAL(NCOMP,NPLATE,NPLATE_FEED,FEED_STREAM,
1 COMP_FEED,TB,xLIQ_FEED,xLIQ,VAP,LIQ,FEED,TOP,BTM,REF,
1 M_VAP,M_LIQ,M_FEED,M_TOP,TEMP)
CALL CALC_ALPHA(NCOMP,NPLATE,ANTOINE,WILLSON,P0,TEMP,xLIQ,ALPHA)
N_JACOB = (2*NCOMP+1)*NPLATE
WRITE(IW,*) ""
WRITE(IW,'(a10,10a10)') "",("----------",I=1,7)
WRITE(IW,'(15x,a35)') "START CALCULATION"
WRITE(IW,'(a10,10a10)') "",("----------",I=1,7)
WRITE(IW,*) ""
WRITE(IW,'(20x,a10,x,a16)') "ITERATION","ERROR"
WRITE(IW,*) ""
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCC
CCC LOOP
CCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
DO ITE = 1,2000
CALL CALC_ALPHA(NCOMP,NPLATE,ANTOINE,WILLSON,P0
1 ,TEMP,xLIQ,ALPHA)
CALL CALC_ERROR(NPLATE,NCOMP,TEMP,REF,TOP,BTM,
1 M_VAP,M_LIQ,VAP,LIQ,M_FEED,M_TOP,ALPHA,CP,E,M,Q)
DO II = 1,(2*NCOMP+1)*NPLATE
DO JJ = 1,(2*NCOMP+1)*NPLATE
JACOB(II,JJ) = 0.0d0
ENDDO
ENDDO
DO II = 1,(2*NCOMP+1)*NPLATE
FX(II) = 0.0d0
ENDDO
DO JJ = 1,(2*NCOMP+1)*NPLATE
CALL CALC_DEL_VEC(JJ,NCOMP,NPLATE,delta
1 ,TEMP,M_VAP,M_LIQ,D_TEMP,D_M_VAP,D_M_LIQ)
CALL CALC_LIQ_VAP_xLIQ(NCOMP,NPLATE
1 ,D_VAP,D_LIQ,D_M_VAP,D_M_LIQ,D_xLIQ)
CALL CALC_ALPHA(NCOMP,NPLATE,ANTOINE,WILLSON,P0
1 ,D_TEMP,D_xLIQ,D_ALPHA)
CALL CALC_ERROR(NPLATE,NCOMP,TEMP,REF,TOP,BTM
1 ,D_M_VAP,D_M_LIQ,D_VAP,D_LIQ,M_FEED,M_TOP,D_ALPHA,CP
1 ,D_E,D_M,D_Q)
CALL CALC_JACOB_COL(JJ,NCOMP,NPLATE
1 ,E,M,Q,D_E,D_M,D_Q,delta,JACOB,FX)
ENDDO
CALL DGESV(N_JACOB,1,JACOB,N_JACOB,IPIV,FX,N_JACOB,INFO)
CALL CALC_NEW_VEC(NCOMP,NPLATE
1 ,TEMP,M_VAP,M_LIQ,M_TOP,VAP,LIQ,FX,xLIQ,TOP)
DO J = 1,NPLATE
working(J) = abs(E(J))
DO I = 1,NCOMP
working(J) = working(J) + abs(M(I,J)) + abs(Q(I,J))
ENDDO
ENDDO
allError = 0.0d0
DO J = 1,NPLATE
allError = allError + working(J)
ENDDO
IF(allError .LE. REAL(NPLATE) * ErrorTolerance) THEN
GOTO 100
ELSE IF(ITE<10) THEN
GOTO 200
ELSE IF(ITE .LE. 2000 .AND. MOD(ITE,10) .EQ.0) THEN
GOTO 200
ELSE
GOTO 300
ENDIF
200 WRITE(IW,'(15x,i10,10x,e15.5)') ITE,allError/real(NPLATE)
300 CONTINUE
CCCCCCCCCCCCCCCCCCCC
CCCC
CCCC ITE LOOP END
ENDDO
CCCC
CCCCCCCCCCCCCCCCCCCC
100 WRITE(IW,'(15x,i10,10x,e15.5,a23)') ITE,allError/real(NPLATE)
1 ,"< 1.0 E-05 (Converged)"
CALL WRITE_RESULT(IW,NCOMP,NPLATE,NPLATE_FEED,FEED_STREAM,
1 MATERIAL,CHARA,COMP_FEED,VAP,LIQ,FEED,TOP,TEMP,
1 xLIQ,xLIQ_FEED,ALPHA)
END PROGRAM
SUBROUTINE CALC_ALPHA
1 (NCOMP,NPLATE,ANTOINE,WILLSON,P0,TEMP,xLIQ,ALPHA)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
REAL(8) xLIQ(NCOMP,NPLATE)
REAL(8) ALPHA(NCOMP,NPLATE)
REAL(8) ANTOINE(3,NCOMP)
REAL(8) WILLSON(NCOMP,NCOMP)
REAL(8) P(NCOMP)
REAL(8) GAM(NCOMP,NPLATE)
REAL(8) TEMP(NPLATE)
DO JJ = 1, NPLATE
DO I= 1, NCOMP
workingFirst = 0.0d0
DO J = 1, NCOMP
workingFirst = workingFirst + xLIQ(J,JJ) * WILLSON(I,J)
ENDDO
workingThird = 0.0d0
DO K = 1, NCOMP
workingBunbo = 0.0d0
DO J = 1, NCOMP
workingBunbo = workingBunbo
1 + xLIQ(J,JJ)*WILLSON(K,J)
ENDDO
workingThird = workingThird
1 + xLIQ(K,JJ)*WILLSON(K,I)/workingBunbo
ENDDO
GAM(I,JJ) =
1 EXP(-LOG(workingFirst) + 1.0d0 -workingThird)
ENDDO
ENDDO
DO J = 1, NPLATE
DO I = 1, NCOMP
P(I) = ANTOINE(1,I) - ANTOINE(2,I)/(TEMP(J)+ANTOINE(3,I))
P(I) = (10.0d0**(P(I))) * P0 / 760.0d0
ENDDO
DO I = 1 , NCOMP
ALPHA(I,J) = GAM(I,J)*P(I)/P0
ENDDO
c DO I = 1 , NCOMP
c ALPHA(I,J) = P(I)/P0
c ENDDO
ENDDO
ENDSUBROUTINE
SUBROUTINE SET_PROPERTY(NCOMP,MATERIAL,ANTOINE,WILLSON,CP,TB)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
REAL(8) DB_ANTOINE(3,10),DB_WILLSON(10,10)
REAL(8) DB_CP(3,10),DB_TB(10)
REAL(8) ANTOINE(3,NCOMP),WILLSON(NCOMP,NCOMP)
REAL(8) CP(3,NCOMP),TB(NCOMP)
INTEGER MATERIAL(10)
DB_ANTOINE(1,1) = 7.87863d0
DB_ANTOINE(2,1) = 1473.110d0
DB_ANTOINE(3,1) = 230.0d0
DB_ANTOINE(1,2) = 8.04494d0
DB_ANTOINE(2,2) = 1554.300d0
DB_ANTOINE(3,2) = 222.64990
DB_ANTOINE(1,3) = 7.96680d0
DB_ANTOINE(2,3) = 1668.209d0
DB_ANTOINE(3,3) = 228.000d0
DB_ANTOINE(1,4) = 7.99733d0
DB_ANTOINE(2,4) = 1569.699d0
DB_ANTOINE(3,4) = 209.50000d0
DB_ANTOINE(1,5) = 6.66040d0
DB_ANTOINE(2,5) = 813.055d0
DB_ANTOINE(3,5) = 132.92900d0
DB_ANTOINE(1,6) = 9.13620d0
DB_ANTOINE(2,6) = 2443.000d0
DB_ANTOINE(3,6) = 273.010d0
DB_ANTOINE(1,7) = 8.13590d0
DB_ANTOINE(2,7) = 1582.399d0
DB_ANTOINE(3,7) = 218.89900
DB_ANTOINE(1,8) = 7.02440d0
DB_ANTOINE(2,8) = 1161.000d0
DB_ANTOINE(3,8) = 224.000d0
DB_ANTOINE(1,9) = 6.97420d0
DB_ANTOINE(2,9) = 1209.600d0
DB_ANTOINE(3,9) = 216.00000d0
DB_ANTOINE(1,10) = 7.09800d0
DB_ANTOINE(2,10) = 1238.600d0
DB_ANTOINE(3,10) = 217.00000d0
DB_WILLSON(1,1) = 1.00000d0
DB_WILLSON(1,2) = 0.99990d0
DB_WILLSON(1,3) = 0.38208d0
DB_WILLSON(1,4) = 1.00000d0
DB_WILLSON(1,5) = 2.33400d0
DB_WILLSON(1,6) = 0.24260d0
DB_WILLSON(1,7) = 0.17570d0
DB_WILLSON(1,8) = 0.14900d0
DB_WILLSON(1,9) = 0.38470d0
DB_WILLSON(1,10) = 0.73290d0
DB_WILLSON(2,1) = 1.00000d0
DB_WILLSON(2,2) = 1.00000d0
DB_WILLSON(2,3) = 0.05840d0
DB_WILLSON(2,4) = 1.00000d0
DB_WILLSON(2,5) = 1.00000d0
DB_WILLSON(2,6) = 1.30100d0
DB_WILLSON(2,7) = 1.00000d0
DB_WILLSON(2,8) = 0.78930d0
DB_WILLSON(2,9) = 0.33240d0
DB_WILLSON(2,10) = 0.73268d0
DB_WILLSON(3,1) = 1.09178d0
DB_WILLSON(3,2) = 0.99900d0
DB_WILLSON(3,3) = 1.00000d0
DB_WILLSON(3,4) = 0.71000d0
DB_WILLSON(3,5) = 0.78600d0
DB_WILLSON(3,6) = 0.71900d0
DB_WILLSON(3,7) = 0.95190d0
DB_WILLSON(3,8) = 0.37930d0
DB_WILLSON(3,9) = 0.43380d0
DB_WILLSON(3,10) = 0.28579d0
DB_WILLSON(4,1) = 1.00000d0
DB_WILLSON(4,2) = 1.00000d0
DB_WILLSON(4,3) = 0.03030d0
DB_WILLSON(4,4) = 1.00000d0
DB_WILLSON(4,5) = 1.89200d0
DB_WILLSON(4,6) = 1.00000d0
DB_WILLSON(4,7) = 1.00000d0
DB_WILLSON(4,8) = 0.64770d0
DB_WILLSON(4,9) = 0.66280d0
DB_WILLSON(4,10) = 0.61890d0
DB_WILLSON(5,1) = 0.28021d0
DB_WILLSON(5,2) = 1.00000d0
DB_WILLSON(5,3) = 0.07761d0
DB_WILLSON(5,4) = 0.36230d0
DB_WILLSON(5,5) = 1.00000d0
DB_WILLSON(5,6) = 1.00000d0
DB_WILLSON(5,7) = 1.00000d0
DB_WILLSON(5,8) = 0.53420d0
DB_WILLSON(5,9) = 0.77660d0
DB_WILLSON(5,10) = 0.71728d0
DB_WILLSON(6,1) = 1.88550d0
DB_WILLSON(6,2) = 0.73660d0
DB_WILLSON(6,3) = 0.02129d0
DB_WILLSON(6,4) = 1.00000d0
DB_WILLSON(6,5) = 1.00000d0
DB_WILLSON(6,6) = 1.00000d0
DB_WILLSON(6,7) = 1.00000d0
DB_WILLSON(6,8) = 0.94060d0
DB_WILLSON(6,9) = 0.74523d0
DB_WILLSON(6,10) = 0.62100d0
DB_WILLSON(7,1) = 2.99900d0
DB_WILLSON(7,2) = 1.00000d0
DB_WILLSON(7,3) = 0.04919d0
DB_WILLSON(7,4) = 1.00000d0
DB_WILLSON(7,5) = 1.00000d0
DB_WILLSON(7,6) = 1.00000d0
DB_WILLSON(7,7) = 1.00000d0
DB_WILLSON(7,8) = 1.00000d0
DB_WILLSON(7,9) = 0.74523d0
DB_WILLSON(7,10) = 0.62100d0
DB_WILLSON(8,1) = 2.07000d0
DB_WILLSON(8,2) = 0.70593d0
DB_WILLSON(8,3) = 0.21130d0
DB_WILLSON(8,4) = 0.74020d0
DB_WILLSON(8,5) = 0.88495d0
DB_WILLSON(8,6) = 0.53050d0
DB_WILLSON(8,7) = 1.00000d0
DB_WILLSON(8,8) = 1.00000d0
DB_WILLSON(8,9) = 1.56800d0
DB_WILLSON(8,10) = 1.01000d0
DB_WILLSON(9,1) = 0.92521d0
DB_WILLSON(9,2) = 1.14630d0
DB_WILLSON(9,3) = 0.02512d0
DB_WILLSON(9,4) = 0.73726d0
DB_WILLSON(9,5) = 0.77660d0
DB_WILLSON(9,6) = 0.92100d0
DB_WILLSON(9,7) = 0.92100d0
DB_WILLSON(9,8) = 0.54900d0
DB_WILLSON(9,9) = 1.00000d0
DB_WILLSON(9,10) = 0.76624d0
DB_WILLSON(10,1) = 0.17922d0
DB_WILLSON(10,2) = 0.52250d0
DB_WILLSON(10,3) = 0.05525d0
DB_WILLSON(10,4) = 0.65770d0
DB_WILLSON(10,5) = 0.74420d0
DB_WILLSON(10,6) = 0.67630d0
DB_WILLSON(10,7) = 0.67630d0
DB_WILLSON(10,8) = 0.92500d0
DB_WILLSON(10,9) = 0.76624d0
DB_WILLSON(10,10)= 1.00000d0
DB_CP(1,1) = 37.50d0
DB_CP(2,1) = 0.04154d0
DB_CP(3,1) = 0.0d0
DB_CP(1,2) = 38.6d0
DB_CP(2,2) = 0.05972d0
DB_CP(3,2) = 0.0d0
DB_CP(1,3) = 44.0d0
DB_CP(2,3) = 0.03343d0
DB_CP(3,3) = 0.0d0
DB_CP(1,4) = 11.08d0
DB_CP(2,4) = 0.07853d0
DB_CP(3,4) = 0.0d0
DB_CP(1,5) = 11.08d00
DB_CP(2,5) = 0.08173d0
DB_CP(3,5) = 0.0d0
DB_CP(1,6) = 45.88d0
DB_CP(2,6) = 0.09743d0
DB_CP(3,6) = 0.0d0
DB_CP(1,7) = 45.88d0
DB_CP(2,7) = 0.10533d0
DB_CP(3,7) = 0.0d0
DB_CP(1,8) = 29.04d0
DB_CP(2,8) = 0.07375d0
DB_CP(3,8) = 0.0d0
DB_CP(1,9) = 31.99d0
DB_CP(2,9) = 0.09246d0
DB_CP(3,9) = 0.0d0
DB_CP(1,10) = 34.25d0
DB_CP(2,10) = 0.10453d0
DB_CP(3,10) = 0.0d0
DB_TB(1) = 64.7d0
DB_TB(2) = 78.4d0
DB_TB(3) = 100.0d0
DB_TB(4) = 97.0d0
DB_TB(5) = 82.3d0
DB_TB(6) = 111.7d0
DB_TB(7) = 82.3d0
DB_TB(8) = 56.0d0
DB_TB(9) = 79.6d0
DB_TB(10) = 77.1d0
DO I = 1,NCOMP
ANTOINE(1,I) = DB_ANTOINE(1,MATERIAL(I))
ANTOINE(2,I) = DB_ANTOINE(2,MATERIAL(I))
ANTOINE(3,I) = DB_ANTOINE(3,MATERIAL(I))
ENDDO
DO I = 1,NCOMP
DO J = 1,NCOMP
WILLSON(I,J) = DB_WILLSON(MATERIAL(I),MATERIAL(J))
ENDDO
ENDDO
DO I = 1,NCOMP
CP(1,I) = DB_CP(1,MATERIAL(I))
CP(2,I) = DB_CP(2,MATERIAL(I))
CP(3,I) = DB_CP(3,MATERIAL(I))
ENDDO
DO I = 1,NCOMP
TB(I) = DB_TB(MATERIAL(I))
TB(I) = DB_TB(MATERIAL(I))
TB(I) = DB_TB(MATERIAL(I))
ENDDO
ENDSUBROUTINE
SUBROUTINE READ_INPUT(input_data,NCOMP,NPLATE,NPLATE_FEED,
1 FEED_STREAM,TOP,REF,BTM,MATERIAL,COMP_FEED,CHARA)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
REAL(8) work
REAL(8) COMP_FEED(10)
INTEGER MATERIAL(10)
CHARACTER(16) input_data
CHARACTER(16) c_work
CHARACTER(18) CHARA(10)
OPEN(1,FILE=input_data)
READ(1,*) c_work, NCOMP
IF(NCOMP.EQ.2) THEN
READ(1,*) c_work,
1 MATERIAL(1),MATERIAL(2)
READ(1,*) c_work, COMP_FEED(1),COMP_FEED(2)
ELSEIF(NCOMP.EQ.3) THEN
READ(1,*) c_work,
1 MATERIAL(1),MATERIAL(2),MATERIAL(3)
READ(1,*) c_work, COMP_FEED(1),COMP_FEED(2),COMP_FEED(3)
ELSEIF(NCOMP.EQ.4) THEN
READ(1,*) c_work,
1 MATERIAL(1),MATERIAL(2),MATERIAL(3)
1 ,MATERIAL(4)
READ(1,*) c_work, COMP_FEED(1),COMP_FEED(2),COMP_FEED(3)
1 ,COMP_FEED(4)
ELSEIF(NCOMP.EQ.5) THEN
READ(1,*) c_work,
1 MATERIAL(1),MATERIAL(2),MATERIAL(3)
1 ,MATERIAL(4),MATERIAL(5)
READ(1,*) c_work, COMP_FEED(1),COMP_FEED(2),COMP_FEED(3)
1 ,COMP_FEED(4),COMP_FEED(5)
ELSE
WRITE(IW,*) "ERROR 2<=NCOMP<=5"
STOP
ENDIF
work = 0.0d0
DO I=1,NCOMP
work = work + COMP_FEED(I)
ENDDO
DO I = 1,NCOMP
COMP_FEED(I) = COMP_FEED(I) / work
ENDDO
DO I=1,NCOMP
IF(MATERIAL(I).GT.10) THEN
WRITE(IW,*) "ERROR MATERIAL_NUMBER <= 10"
STOP
ENDIF
ENDDO
READ(1,*) c_work,NPLATE
READ(1,*) c_work,NPLATE_FEED
READ(1,*) c_work,FEED_STREAM
READ(1,*) c_work,TOP
READ(1,*) c_work,REF
BTM = FEED_STREAM - TOP
CHARA(1) = "MeOH"
CHARA(2) = "EtOH"
CHARA(3) = "WATER"
CHARA(4) = "1-PrOH"
CHARA(5) = "2-PrOH"
CHARA(6) = "1-BuOH"
CHARA(7) = "t-BuOH"
CHARA(8) = "Acetone"
CHARA(9) = "MEK"
CHARA(10)= "AcOH"
ENDSUBROUTINE
SUBROUTINE SET_INITIAL(NCOMP,NPLATE,NPLATE_FEED,FEED_STREAM,
1 COMP_FEED,TB,xLIQ_FEED,xLIQ,VAP,LIQ,FEED,TOP,BTM,REF,
1 M_VAP,M_LIQ,M_FEED,M_TOP,TEMP)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
REAL(8) FEED(NPLATE),LIQ(NPLATE),VAP(NPLATE)
REAL(8) xLIQ_FEED(NCOMP),xLIQ(NCOMP,NPLATE)
REAL(8) M_VAP(NCOMP,NPLATE),M_LIQ(NCOMP,NPLATE)
REAL(8) M_FEED(NCOMP,NPLATE),M_TOP(NCOMP)
REAL(8) TEMP(NPLATE)
REAL(8) COMP_FEED(10)
REAL(8) TB(NCOMP)
DO I = 1, NPLATE
FEED(I) = 0.0d0
ENDDO
FEED(NPLATE_FEED)=FEED_STREAM
VAP(1) = 0.0d0
DO I = 2, NPLATE
VAP(I) = (REF + 1.0d0) * TOP
ENDDO
DO I = 1, NPLATE-1
IF(I.EQ.1) THEN
LIQ(1) = REF * TOP
ELSE
LIQ(I) = LIQ(I-1) + FEED(I)
ENDIF
ENDDO
LIQ(NPLATE) = FEED_STREAM-TOP
T_MIN = TB(1)
T_MAX = TB(1)
DO I = 2, NCOMP
IF(TB(I).LT.T_MIN) THEN
T_MIN = TB(I)
ENDIF
IF(TB(I).GT.T_MAX) THEN
T_MAX = TB(I)
ENDIF
ENDDO
DO I = 1, NPLATE
TEMP(I) = T_MIN + (T_MAX - T_MIN) / NPLATE * real(I)
ENDDO
DO J = 1,NPLATE
DO I = 1,NCOMP
M_VAP (I,J) = xLIQ(I,J)*VAP(J)
M_LIQ (I,J) = xLIQ(I,J)*LIQ(J)
M_FEED(I,J) = xLIQ_FEED(I)*FEED(J)
ENDDO
ENDDO
DO I = 1,NCOMP
M_TOP(I) = TOP*xLIQ(I,1)
ENDDO
DO J = 1,NPLATE
LIQ(J) = 0.0d0
VAP(J) = 0.0d0
ENDDO
DO J = 1,NPLATE
DO I = 1,NCOMP
LIQ(J) = LIQ(J) + M_LIQ(I,J)
VAP(J) = VAP(J) + M_VAP(I,J)
ENDDO
ENDDO
DO I = 1,NCOMP
DO J = 1,NPLATE
xLIQ(I,J) = M_LIQ(I,J)/LIQ(J)
ENDDO
ENDDO
ENDSUBROUTINE
SUBROUTINE CALC_ERROR(NPLATE,NCOMP,TEMP,REF,TOP,BTM,
1 M_VAP,M_LIQ,VAP,LIQ,M_FEED,M_TOP,ALPHA,CP,E,M,Q)
IMPLICIT NONE
INTEGER I,J
INTEGER NPLATE,NCOMP
DOUBLE PRECISION H_VAP(NPLATE),H_LIQ(NPLATE),H_FEED(NPLATE)
DOUBLE PRECISION M_VAP(NCOMP,NPLATE),M_LIQ(NCOMP,NPLATE)
DOUBLE PRECISION M_TOP(NCOMP),M_FEED(NCOMP,NPLATE)
DOUBLE PRECISION ALPHA(NCOMP,NPLATE)
DOUBLE PRECISION LIQ(NPLATE),VAP(NPLATE)
DOUBLE PRECISION TEMP(NPLATE)
DOUBLE PRECISION CP(3,NCOMP)
DOUBLE PRECISION E(NPLATE),M(NCOMP,NPLATE),Q(NCOMP,NPLATE)
DOUBLE PRECISION REF,TOP,BTM
DO J = 1, NPLATE
H_VAP(J) = 0.0d0
H_LIQ(J) = 0.0d0
H_FEED(J) = 0.0d0
DO I = 1, NCOMP
H_VAP(J) = H_VAP(J) + M_VAP(I,J)*(CP(1,I)+TEMP(J)*CP(2,I))
H_LIQ(J) = H_LIQ(J) + M_LIQ(I,J) * CP(3,I)
H_FEED(J) = H_FEED(J)+ M_FEED(I,J)* CP(3,I)
ENDDO
ENDDO
DO I=1,NCOMP
M(I,1) = M_VAP(I,1) + M_LIQ(I,1) + M_TOP(I)
1 - (M_VAP(I,2))
DO J = 2,NPLATE-1
M(I,J) = M_VAP(I,J) + M_LIQ(I,J)
1 - (M_VAP(I,J+1) + M_LIQ(I,J-1) + M_FEED(I,J))
ENDDO
M(I,NPLATE) = M_VAP(I,NPLATE) + M_LIQ(I,NPLATE)
1 - M_LIQ(I,NPLATE-1)
ENDDO
E(1) = LIQ(1) - REF * ( TOP)
E(NPLATE) = LIQ(NPLATE) - BTM
DO I = 2, NPLATE-1
E(I) = H_VAP(I) + H_LIQ(I)
1 - (H_VAP(I+1) + H_LIQ(I-1) + H_FEED(I))
ENDDO
E(NPLATE) = LIQ(NPLATE) - BTM
Q(1,1) = LIQ(1)
DO I = 1,NCOMP
Q(1,1) = Q(1,1) - ALPHA(I,1)*M_LIQ(I,1)
ENDDO
DO I = 2,NCOMP
Q(I,1) = M_LIQ(I,1)*VAP(2) - LIQ(1)*M_VAP(I,2)
ENDDO
DO J = 2,NPLATE
DO I = 1,NCOMP
Q(I,J) = ALPHA(I,J)*VAP(J)/LIQ(J)*M_LIQ(I,J) - M_VAP(I,J)
ENDDO
ENDDO
ENDSUBROUTINE
SUBROUTINE CALC_LIQ_VAP_xLIQ(NCOMP,NPLATE
1 ,VAP,LIQ,M_VAP,M_LIQ,xLIQ)
IMPLICIT NONE
INTEGER I,J
INTEGER NPLATE,NCOMP
DOUBLE PRECISION VAP(NPLATE),LIQ(NPLATE)
DOUBLE PRECISION M_VAP(NCOMP,NPLATE),M_LIQ(NCOMP,NPLATE)
DOUBLE PRECISION xLIQ(NCOMP,NPLATE)
DO J = 1,NPLATE
LIQ(J) = 0.0d0
VAP(J) = 0.0d0
ENDDO
DO J = 1,NPLATE
DO I = 1,NCOMP
LIQ(J) = LIQ(J) + M_LIQ(I,J)
VAP(J) = VAP(J) + M_VAP(I,J)
ENDDO
ENDDO
DO I = 1,NCOMP
DO J = 1,NPLATE
xLIQ(I,J) = M_LIQ(I,J)/LIQ(J)
ENDDO
ENDDO
ENDSUBROUTINE
SUBROUTINE CALC_JACOB_COL(JJ,NCOMP,NPLATE,E,M,Q,D_E,D_M,D_Q,
1 delta,JACOB,FX)
IMPLICIT NONE
INTEGER I,J,II,JJ,KK,LL,JJJ
INTEGER NPLATE,NCOMP
DOUBLE PRECISION E(NPLATE)
DOUBLE PRECISION M(NCOMP,NPLATE),Q(NCOMP,NPLATE)
DOUBLE PRECISION D_E(NPLATE)
DOUBLE PRECISION D_M(NCOMP,NPLATE),D_Q(NCOMP,NPLATE)
DOUBLE PRECISION xLIQ(NCOMP,NPLATE)
DOUBLE PRECISION JACOB((2*NCOMP+1)*NPLATE,(2*NCOMP+1)*NPLATE)
DOUBLE PRECISION FX((2*NCOMP+1)*NPLATE)
DOUBLE PRECISION delta
c
IF(JJ.LE.NPLATE) THEN
JJJ = JJ - NPLATE
JJJ = (JJJ+NCOMP-1)/NCOMP
ELSE IF(JJ.LE.((NCOMP+1)*NPLATE)) THEN
II = mod((JJ-NPLATE),NCOMP)
IF(II.EQ.0) THEN
II = NCOMP
ENDIF
JJJ = JJ - NPLATE
JJJ = (JJJ+NCOMP-1)/NCOMP
ELSE
II = mod((JJ-NPLATE),NCOMP)
IF(II.EQ.0) THEN
II = NCOMP
ENDIF
JJJ = JJ - NPLATE - NPLATE * NCOMP
JJJ = (JJJ+NCOMP-1)/NCOMP
c
ENDIF
IF(JJ.LE.NPLATE) THEN
DO KK = 1,NPLATE
JACOB(KK,JJ)
1 = (D_E(KK)-E(KK)) /delta
DO LL = 1,NCOMP
JACOB(NPLATE+NCOMP*(KK-1)+LL,JJ)
1 = (D_M(LL,KK)-M(LL,KK))/delta
ENDDO
DO LL = 1,NCOMP
JACOB((NCOMP+1)*NPLATE+NCOMP*(KK-1)+LL,JJ)
1 = (D_Q(LL,KK)-Q(LL,KK))/delta
ENDDO
ENDDO
FX(JJ ) = -E(JJ)
ELSE IF(JJ.LE.((NCOMP+1)*NPLATE)) THEN
DO KK = 1,NPLATE
JACOB(KK, JJ)
1 = (D_E(KK)-E(KK)) /delta
DO LL = 1,NCOMP
JACOB(NPLATE+NCOMP*(KK-1)+LL,JJ)
1 = (D_M(LL,KK)-M(LL,KK))/delta
ENDDO
DO LL = 1,NCOMP
JACOB((NCOMP+1)*NPLATE+NCOMP*(KK-1)+LL ,JJ)
1 = (D_Q(LL,KK)-Q(LL,KK))/delta
ENDDO
ENDDO
FX(JJ) = -M(II,JJJ)
c
c
c
ELSE
DO KK = 1,NPLATE
JACOB(KK, JJ)
1 = (D_E(KK)-E(KK)) /delta
DO LL = 1,NCOMP
JACOB(NPLATE+NCOMP*(KK-1)+LL
1 ,JJ)
1 = (D_M(LL,KK)-M(LL,KK))/delta
ENDDO
DO LL = 1,NCOMP
JACOB((NCOMP+1)*NPLATE+NCOMP*(KK-1)+LL
1 ,JJ)
1 = (D_Q(LL,KK)-Q(LL,KK))/delta
ENDDO
ENDDO
FX(JJ) = -Q(II,JJJ)
ENDIF
ENDSUBROUTINE
SUBROUTINE CALC_NEW_VEC(NCOMP,NPLATE
1 ,TEMP,M_VAP,M_LIQ,M_TOP,VAP,LIQ,FX,xLIQ,TOP)
IMPLICIT NONE
INTEGER I,J,II,JJ,KK,LL,JJJ
INTEGER NPLATE,NCOMP
DOUBLE PRECISION xLIQ(NCOMP,NPLATE)
DOUBLE PRECISION M_VAP(NCOMP,NPLATE),M_LIQ(NCOMP,NPLATE)
DOUBLE PRECISION M_TOP(NCOMP)
DOUBLE PRECISION TEMP(NPLATE)
DOUBLE PRECISION VAP(NPLATE),LIQ(NPLATE)
DOUBLE PRECISION FX((2*NCOMP+1)*NPLATE)
DOUBLE PRECISION TOP
DO J = 1,NPLATE
TEMP(J) = TEMP(J)
1 + FX(J)/10.0d0
ENDDO
DO I = 1,NCOMP
DO J = 1,NPLATE
M_VAP(I,J) = M_VAP(I,J)
1 + FX(NPLATE +NCOMP*(J-1)+I)/10.0d0
ENDDO
ENDDO
DO I = 1,NCOMP
DO J = 1,NPLATE
M_LIQ(I,J) = M_LIQ(I,J)
1 + FX((NCOMP+1)*NPLATE+NCOMP*(J-1)+I)/10.0d0
ENDDO
ENDDO
DO I = 1, NCOMP
M_TOP(I) = TOP * xLIQ(I,1)
ENDDO
DO I = 1,NCOMP
DO J = 1,NPLATE
IF(M_VAP(I,J).LE.0.0d0) THEN
M_VAP(I,J) = 0.0d0
ENDIF
IF(M_LIQ(I,J).LE.0.0d0) THEN
M_LIQ(I,J) = 0.0d0
ENDIF
M_VAP(I,1) = 0.0d0
ENDDO
ENDDO
DO J = 1,NPLATE
LIQ(J) = 0.0d0
VAP(J) = 0.0d0
ENDDO
DO J = 1,NPLATE
DO I = 1,NCOMP
LIQ(J) = LIQ(J) + M_LIQ(I,J)
VAP(J) = VAP(J) + M_VAP(I,J)
ENDDO
ENDDO
DO I = 1,NCOMP
DO J = 1,NPLATE
xLIQ(I,J) = M_LIQ(I,J)/LIQ(J)
ENDDO
ENDDO
ENDSUBROUTINE
SUBROUTINE CALC_DEL_VEC(JJ,NCOMP,NPLATE,delta
1 ,TEMP,M_VAP,M_LIQ,D_TEMP,D_M_VAP,D_M_LIQ)
IMPLICIT NONE
INTEGER I,J,II,JJ,KK,LL,JJJ
INTEGER NPLATE,NCOMP
DOUBLE PRECISION M_VAP(NCOMP,NPLATE),M_LIQ(NCOMP,NPLATE)
DOUBLE PRECISION TEMP(NPLATE)
DOUBLE PRECISION D_M_VAP(NCOMP,NPLATE),D_M_LIQ(NCOMP,NPLATE)
DOUBLE PRECISION D_TEMP(NPLATE)
DOUBLE PRECISION delta
DO I = 1,NCOMP
DO J = 1,NPLATE
D_M_VAP(I,J) = M_VAP(I,J)
D_M_LIQ(I,J) = M_LIQ(I,J)
ENDDO
ENDDO
DO J = 1,NPLATE
D_TEMP(J) = TEMP(J)
ENDDO
IF(JJ.LE.NPLATE) THEN
D_TEMP(JJ) = D_TEMP(JJ) + delta
ELSE IF(JJ.LE.((NCOMP+1)*NPLATE)) THEN
II = mod((JJ-NPLATE),NCOMP)
IF(II.EQ.0) THEN
II = NCOMP
ENDIF
JJJ = JJ - NPLATE
JJJ = (JJJ+NCOMP-1)/NCOMP
D_M_VAP(II,JJJ) = D_M_VAP(II,JJJ) + delta
ELSE
II = mod((JJ-NPLATE),NCOMP)
IF(II.EQ.0) THEN
II = NCOMP
ENDIF
JJJ = JJ - NPLATE - NPLATE * NCOMP
JJJ = (JJJ+NCOMP-1)/NCOMP
D_M_LIQ(II,JJJ) = D_M_LIQ(II,JJJ) + delta
ENDIF
ENDSUBROUTINE
SUBROUTINE SET_xLIQ_START(NCOMP,NPLATE,COMP_FEED,xLIQ_FEED,xLIQ,
1 TB,NPLATE_FEED,FEED_STREAM,TOP,BTM)
IMPLICIT NONE
INTEGER I,J,K,LL,JJJ
INTEGER NPLATE,NCOMP,NPLATE_FEED
INTEGER N_TB(NCOMP)
DOUBLE PRECISION xLIQ(NCOMP,NPLATE),xLIQ_FEED(NCOMP)
DOUBLE PRECISION COMP_FEED(NCOMP)
DOUBLE PRECISION TB(NCOMP),TB_work(NCOMP)
DOUBLE PRECISION FEED_STREAM,TOP,BTM
DOUBLE PRECISION T_SAVE
DO I = 1, NCOMP
xLIQ_FEED(I) = COMP_FEED(I)
ENDDO
DO I = 1, NPLATE
DO J = 1, NCOMP
xLIQ(J,I) = xLIQ_FEED(J)
ENDDO
ENDDO
DO I = 1 , NCOMP
N_TB(I) = I
TB_work(I) = TB(I)
ENDDO
DO I = 1, NCOMP-1
DO J = I+1, NCOMP
IF(TB_work(I).GT.TB_work(J)) THEN
T_SAVE = TB_work(I)
TB_work(I) = TB_work(J)
TB_work(J) = T_SAVE
K = N_TB(I)
N_TB(I) = N_TB(J)
N_TB(J) = K
ENDIF
ENDDO
ENDDO
DO I = 1, NCOMP/2
J = N_TB(I)
xLIQ(J,1) = FEED_STREAM/TOP*xLIQ_FEED(J)
IF(xLIQ(J,1).GE.1.0d0) THEN
xLIQ(J,1) = 1.0d0
ENDIF
xLIQ(J,NPLATE) = 0.00d0
ENDDO
DO I = NCOMP/2+1,NCOMP
J = N_TB(I)
xLIQ(J,1) = 0.0d0
xLIQ(J,NPLATE) = FEED_STREAM/BTM*xLIQ_FEED(J)
IF(xLIQ(J,NPLATE).GE.1.0d0) THEN
xLIQ(J,NPLATE) = 1.0d0
ENDIF
ENDDO
DO J = 2,NPLATE_FEED
DO I = 1,NCOMP
xLIQ(I,J) = (xLIQ(I,NPLATE_FEED)-xLIQ(I,1))
1 / real(NPLATE_FEED-1)*real(J) + xLIQ(I,1)
ENDDO
ENDDO
DO J = NPLATE_FEED,NPLATE-1
DO I = 1,NCOMP
xLIQ(I,J) = (xLIQ(I,NPLATE)-xLIQ(I,NPLATE_FEED))
1 / real(NPLATE-NPLATE_FEED)*real(J-NPLATE_FEED) + xLIQ_FEED(I)
ENDDO
ENDDO
ENDSUBROUTINE
SUBROUTINE WRITE_RESULT(IW,NCOMP,NPLATE,NPLATE_FEED,FEED_STREAM,
1 MATERIAL,CHARA,COMP_FEED,VAP,LIQ,FEED,TOP,TEMP,
1 xLIQ,xLIQ_FEED,ALPHA)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
REAL(8) FEED(NPLATE),LIQ(NPLATE),VAP(NPLATE)
REAL(8) xLIQ_FEED(NCOMP),xLIQ(NCOMP,NPLATE)
REAL(8) ALPHA(NCOMP,NPLATE)
REAL(8) TEMP(NPLATE)
REAL(8) COMP_FEED(10)
INTEGER MATERIAL(10)
CHARACTER(18) CHARA(10)
WRITE(IW,*) ""
WRITE(IW,'(a10,10a10)') "",("----------",I=1,7)
WRITE(IW,'(15x,a35)') "END CALCULATION"
WRITE(IW,'(a10,10a10)') "",("----------",I=1,7)
WRITE(IW,*) ""
WRITE(IW,*) ""
WRITE(IW,'(10x,a17)') "STREAM (mol/hr)"
WRITE(IW,*) ""
WRITE(IW,'(4a20)') "FEED ","TOP OUT", "BTM OUT","REFLUX"
WRITE(IW,'(4(10x,f10.5))') FEED_STREAM,TOP,LIQ(NPLATE),LIQ(1)
WRITE(IW,*) ""
WRITE(IW,*) ""
WRITE(IW,'(10x,a40)') "CONCENTRATION of COMPONENT (mol ratio)"
WRITE(IW,*) ""
WRITE(IW,'(25x,3a15)') "FEED","TOP","BTM"
DO I = 1,NCOMP
WRITE(IW,'(13x,a15,3f15.5)')
1 CHARA(MATERIAL(I)),COMP_FEED(I),xLIQ(I,1),xLIQ(I,NPLATE)
ENDDO
WRITE(IW,*) ""
WRITE(IW,*) ""
WRITE(IW,'(11x,a20)') "CONDITION IN COLUMN"
WRITE(IW,*) ""
IF(NCOMP.EQ.2) THEN
WRITE(IW,'(10x,a10,a20,15x,3a10)')
1 "PLATE","xLIQ","TEMP","VAP","LIQ"
ELSEIF(NCOMP.EQ.3) THEN
WRITE(IW,'(10x,a10,a20,25x,3a10)')
1 "PLATE","xLIQ","TEMP","VAP","LIQ"
ELSEIF(NCOMP.EQ.4) THEN
WRITE(IW,'(10x,a10,a20,35x,3a10)')
1 "PLATE","xLIQ","TEMP","VAP","LIQ"
ELSEIF(NCOMP.EQ.5) THEN
WRITE(IW,'(10x,a10,a20,45x,3a10)')
1 "PLATE","xLIQ","TEMP","VAP","LIQ"
ENDIF
WRITE(IW,'(30x,5(a9,x))')
1 (CHARA(MATERIAL(J)),J=1,NCOMP)
IF(NCOMP.EQ.2) THEN
WRITE(IW,'(76x,f10.5)') TOP
ELSEIF(NCOMP.EQ.3) THEN
WRITE(IW,'(86x,f10.5)') TOP
ELSEIF(NCOMP.EQ.4) THEN
WRITE(IW,'(96x,f10.5)') TOP
ELSEIF(NCOMP.EQ.5) THEN
WRITE(IW,'(106x,f10.5)') TOP
ENDIF
DO I = 1, NPLATE
IF(NCOMP.EQ.2) THEN
IF(I.EQ.1) THEN
WRITE(IW,'(10x,i10,a,5x,2f10.5,10x,f10.5,a10,f10.5)')
1 I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),"--------",LIQ(I)
ELSEIF(I.EQ.NPLATE_FEED) THEN
WRITE(IW,'(10x,a7,i3,a,5x,2f10.5,10x,3f10.5)')
1 "FEED>>>",I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),VAP(I),LIQ(I)
ELSE
WRITE(IW,'(10x,i10,a,5x,2f10.5,10x,3f10.5)')
1 I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),VAP(I),LIQ(I)
ENDIF
ELSEIF(NCOMP.EQ.3) THEN
IF(I.EQ.1) THEN
WRITE(IW,'(10x,i10,a,5x,3f10.5,10x,f10.5,a10,f10.5)')
1 I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),"--------",LIQ(I)
ELSEIF(I.EQ.NPLATE_FEED) THEN
WRITE(IW,'(10x,a7,i3,a,5x,3f10.5,10x,3f10.5)')
1 "FEED>>>",I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),VAP(I),LIQ(I)
ELSE
WRITE(IW,'(10x,i10,a,5x,3f10.5,10x,3f10.5)')
1 I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),VAP(I),LIQ(I)
ENDIF
ELSEIF(NCOMP.EQ.4) THEN
IF(I.EQ.1) THEN
WRITE(IW,'(10x,i10,a,5x,4f10.5,10x,f10.5,a10,f10.5)')
1 I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),"--------",LIQ(I)
ELSEIF(I.EQ.NPLATE_FEED) THEN
WRITE(IW,'(10x,a7,i3,a,5x,4f10.5,10x,3f10.5)')
1 "FEED>>>",I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),VAP(I),LIQ(I)
ELSE
WRITE(IW,'(10x,i10,a,5x,4f10.5,10x,3f10.5)')
1 I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),VAP(I),LIQ(I)
ENDIF
ELSEIF(NCOMP.EQ.5) THEN
IF(I.EQ.1) THEN
WRITE(IW,'(10x,i10,a,5x,5f10.5,10x,f10.5,a10,f10.5)')
1 I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),"--------",LIQ(I)
ELSEIF(I.EQ.NPLATE_FEED) THEN
WRITE(IW,'(10x,a7,i3,a,5x,5f10.5,10x,3f10.5)')
1 "FEED>>>",I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),VAP(I),LIQ(I)
ELSE
WRITE(IW,'(10x,i10,a,5x,5f10.5,10x,3f10.5)')
1 I,"|",(xLIQ(J,I),J=1,NCOMP),TEMP(I),VAP(I),LIQ(I)
ENDIF
ENDIF
ENDDO
ENDSUBROUTINE
SUBROUTINE Calc_PXY(ANTOINE,WILLSON,PI0,NComp)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
REAL(8) ANTOINE(3,NComp),WILLSON(NComp,NComp)
REAL(8) K_BP(NComp)
REAL(8) xLiqWork(NComp)
IW = 6
WRITE(IW,*) ""
WRITE(IW,'(a10,10a10)') "",("----------",I=1,10)
WRITE(IW,'(15x,a5,a45,a35)') "|","Calc BoilingPoint","|"
WRITE(IW,'(a10,10a10)') "",("----------",I=1,10)
WRITE(IW,*) ""
WRITE(IW,'(8a10,a30)') "I","xLiq","xGas","K_Value",
1 "J","xLiq","xGas","K_Valuc","BP"
DO I = 1, NComp
DO J = I+1, NComp
WRITE(IW,*) ""
DO II = 1, NComp
xLiqWork(II) = 0.00d0
ENDDO
xLiqWork(I) = 0.00d0
xLiqWork(J) = 0.00d0
DO II = 0, 1000
xLiqWork(I) = 0.001d0 * II
xLiqWork(J) = 1.0d0 - xLiqWork(I)
CALL CalcBoilingPoint
1 (xLiqWork,ANTOINE,WILLSON,PI0,NComp,TB,K_BP)
IF(MOD(II,10).EQ.0) THEN
WRITE(IW,'(2(i10,3f10.3),10x,f20.5)')
1 I,xLiqWork(I),K_BP(I)*xLiqWork(I),K_BP(I)
1 ,J,xLiqWork(J),K_BP(J)*xLiqWork(J),K_BP(J)
1 ,TB
ENDIF
ENDDO
ENDDO
ENDDO
ENDSUBROUTINE
SUBROUTINE CalcBoilingPoint(x,ANTO,WILL,PI0,NComp_,TB,K_)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
REAL(8) x(NComp_)
REAL(8) K_(NComp_)
REAL(8) ANTO(3,NComp_)
REAL(8) WILL(NComp_,NComp_)
REAL(8) P(NComp_)
REAL(8) T_STEP(2),ERROR(2),T_c(2)
REAL(8) GAM(NComp_)
REAL(8) T_DIIS(3,3)
DO I= 1, NComp_
workingFirst = 0.0d0
DO J = 1, NComp_
workingFirst = workingFirst + x(J) * WILL(I,J)
ENDDO
workingThird = 0.0d0
DO K = 1, NComp_
workingBunbo = 0.0d0
DO J = 1, NComp_
workingBunbo = workingBunbo + x(J)*WILL(K,J)
ENDDO
workingThird = workingThird + x(K)*WILL(K,I)/workingBunbo
ENDDO
GAM(I) = -log(workingFirst) + 1.0d0 - workingThird
GAM(I) = exp(GAM(I))
ENDDO
TEMP = 200.0d0
DO ITE_TLOOP = 1,15
IF(ITE_TLOOP < 3 ) THEN
T_new = TEMP - 10.0d0
ELSE IF(ITE_TLOOP >= 4) THEN
TEMP = T_new
ENDIF
TEMP = T_new
T_STEP(2) = T_STEP(1)
T_STEP(1) = TEMP
DO I = 1,NCOMP_
P(I) = ANTO(1,I) - ANTO(2,I)/( TEMP + ANTO(3,I))
P(I) = (10.0d0**(P(I))) * 1013.25d0 / 760.0d0
ENDDO
CALL CalcTotalPressure(x,P,GAM,PI,NCOMP_)
ERROR(2) = ERROR(1)
ERROR(1) = PI-PI0
IF(abs(ERROR(1)) < 1.0d0/(10.0d0**9)) THEN
GOTO 100
ENDIF
IF(ITE_TLOOP>1) THEN
DO I = 1,2
T_DIIS(3,I) = -1.0d0
T_DIIS(I,3) = -1.0d0
ENDDO
T_DIIS(3,3) = 0.0d0
DO I = 1,2
DO J = 1,2
T_DIIS(I,J) = ERROR(I)*ERROR(J)
ENDDO
ENDDO
det = T_DIIS(1,1)*T_DIIS(2,2)*T_DIIS(3,3)
1 + T_DIIS(1,2)*T_DIIS(2,3)*T_DIIS(3,1)
1 + T_DIIS(1,3)*T_DIIS(2,1)*T_DIIS(3,2)
1 - T_DIIS(1,3)*T_DIIS(2,2)*T_DIIS(3,1)
1 - T_DIIS(1,2)*T_DIIS(2,1)*T_DIIS(3,3)
1 - T_DIIS(1,1)*T_DIIS(2,3)*T_DIIS(3,2)
T_c(1) = -(T_DIIS(1,2)*T_DIIS(2,3) - T_DIIS(1,3)*T_DIIS(2,2))/det
T_c(2) = -(T_DIIS(1,3)*T_DIIS(2,1) - T_DIIS(1,1)*T_DIIS(2,3))/det
T_new = T_c(1)*T_STEP(1) + T_c(2)*T_STEP(2)
ENDIF
ENDDO
100 CONTINUE
TB = TEMP
DO I = 1,NCOMP_
P(I) = ANTO(1,I) - ANTO(2,I)/( TB + ANTO(3,I))
P(I) = (10.0d0**(P(I))) * 1013.25d0 / 760.0d0
ENDDO
DO I = 1, NComp_
c K_(I) = P(I) * GAM(I) / PI0
K_(I) = GAM(I)
ENDDO
END SUBROUTINE
SUBROUTINE CalcTotalPressure
1 (XComp,PRESSURE,GAMMA_,TOTALPRESSURE,N)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
REAL(8) XComp(N)
REAL(8) PRESSURE(N)
REAL(8) GAMMA_(N)
TOTALPRESSURE = 0.0d0
DO I = 1, N
TOTALPRESSURE = TOTALPRESSURE + XComp(I)*GAMMA_(I)*PRESSURE(I)
ENDDO
END SUBROUTINE
プログラムの利用方法
STEP
プログラムをダウンロード後、拡張子を「.f」に変更
サイトの仕様上、拡張子「f」ではアップロードできなかったので、「txt」ファイルでダウンロードできるようにしています。プログラムの中身を確認後、問題なければ拡張子を「f」に変更してください。プログラム利用後に問題が起きても当サイトでは責任は負いません。
STEP
「COLUMN.f」をコンパイルして実行ファイルを作成する
コンパイルは次のコードで実行できます。（fortranのコンパイラとlapack環境が必要です。）
gfortran COLUMN.f -llapack -lblasSTEP
実行ファイルにinputファイルを読み込ませて計算を実行させる
./a.out input実行ファイルに計算条件を記載したinputを読み込ませることで蒸留計算が実行されます。inputファイルの例は次です。
NComp , 4
Comp , 1,2,3,8
xLiq_FEED , 20,30,10,40
NSTAIRS , 21
NSTAIR_FEED , 13
FEED_IN , 10.0
TOP , 1.5
REF_RATIO , 2.0
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1) = "MEOH" 64.7[C]
2) = "ETOH" 78.4[C]
3) = "WATER" 100.0[C]
4) = "1-Propanol" 97.0[C]
5) = "2-Propanol" 82.3[C]
6) = "1-Butanol" 111.7[C]
7) = "t-ButylAlcol" 82.3[C]
8) = "Acetone" 56.0[C]
9) = "MetylEtylKetone" 79.6[C]
10)= "EtyleAcetete" 77.1[C]
あわせて読みたい
【COCO/ChemSepで】自作プログラムの正しさを確認しました
別ページにある通り、私は自作した蒸留計算プログラム「COLUMN.f」を持っていますが、正しい計算結果を出力しているかが気になっていました。 そこで、無料の化学工学シ…
まずは3成分（メタノール/エタノール/水）の蒸留を考えてみましょう。
微小成分の水を塔底から排出し、塔頂からメタノールを取り出すプロセスとします。
| 成分 | メタノール/エタノール/水 |
| FEED組成(mol比) | 0.45/0.50/0.05 |
| 蒸留段数 | 21 |
| FEED段 | 13 |
| FEED流量(mol/hr) | 10.0 |
| 蒸留塔上部の流出量(mol/hr) | 1.5 |
| 還流比 | 6.0 |
プログラミングを利用して蒸留計算を行った結果は下の図のようになります。
蒸留の計算結果は縦軸に蒸留段数（上が塔頂、下が塔底）を取り、組成と塔内温度を横軸にとるのが一般的です。
慣れてくるとこの軸のほうが見やすく感じるようになったよ。
メタノール/エタノール/水の混合物の場合、水とエタノールのモル濃度が塔底に向かって単純増加、メタノールのモル濃度が塔頂に向かって単純増加していることが分かります。
また、温度は塔頂が約68℃、塔底が約74℃とほぼメタノールとエタノールの沸点に一致しています。
蒸留計算がうまく収束しないと、計算結果の温度分布が単純増加せずに計算が終了することがありますので、温度プロファイルは必ず確認が必要です。
4成分の蒸留計算
次に、MeOH、1-Propanol、2-Butanol、水に対する 蒸留計算を行ってみましょう。
この場合、組成が単純増加、減少しなくなりますよ。
| 成分数 | 4 |
| 成分 | MEOH/2-PrOH/1-BuOH/H2O |
| FEED組成(mol比) | 0.2,0.5,0.25,0.05 |
| 蒸留段数 | 21 |
| FEED段 | 13 |
| FEED流量(mol/hr) | 10.0 |
| 蒸留塔上部の流出量(mol/hr) | 1.5 |
| 還流比 | 2.0 |
計算結果をグラフ化するとこちらになります。
蒸留計算を行うとFEED段の上部と下部では濃度勾配が大きく異なることが分かります。
これはFEED段下部では液流量が増えるためです。
また、2-Propanolは出口モル比は上部で0.4以下、下部で0.5程度ですが、蒸留塔の中上段では0.7程度まで濃縮しています。
蒸留塔内部では一部の成分が濃縮して、副反応による新たな不純物の生成や反応暴走が起こり得るから、注意が必要だね。
蒸留塔における第三成分の濃縮による事故は1991年6月のライオン株式会社の事故などがあります。
共沸系の蒸留計算
最後にMeOH-WATERとEtOH-WATERの蒸留を行い、蒸留塔内の共沸を見てみましょう。
| 成分数 | 2 |
| 成分 | MEOH(またはETOH)/H2O |
| FEED組成(mol比) | 0.75,0.25 |
| 蒸留段数 | 21 |
| FEED段 | 17 |
| FEED流量(mol/hr) | 10.0 |
| 蒸留塔上部の流出量(mol/hr) | 1.0 |
| 還流比 | 5.0 |
MeOH-WATER系では、素早くMeOH 100%液が得られている一方で、EtOH-WATER系はEtOH 88mol%程度からなかなか濃縮が進んでいません。
これは共沸現象と呼ばれており、実プラントでの蒸留においても挙動の予測が難しく、生産における課題にしばしばなります。