' PROGRAM SHRINKIM.BAS
' PROGRAM WRITTEN BY INGEMAR BJERLE MARCH 2002
' Two parallel reactions in fixed bed reacting according to shrinking core.
' This problem cannot be solved with the assumption of a concentration
' profile moving at a constant velocity through the bed that often is done.
' Instead the partial differential equation must be solved for both
' rate equations. Below the program showing the solution is given.
' In case a constant moving concentration profile is assumed the ratio
' C/C0 is constant across the bed area and that is also valid for q/q0.
' When the partial differential equation is applied for each of the
' rate equations the concentration C within a differential volume element
' is allowed to vary across the bed area. After leaving an element
' an average concentration is calculated. The implicit method is used.
' For those interested a 40 pages pdf-file is available. In this also 
' solutions with Matlab and Scientist are given.
' e-mail: ingemar.bjerle@telia.com


    CLS
    DIM QA1(101): DIM QA0(101):DIM CA1(101):DIM CA0(101)
    DIM QB1(101): DIM QB0(101):DIM CB1(101):DIM CB0(101)
    open "c:\IMP-80-20" for output as #1 'WRITE TO FILE

    NN=10
    K2=50    'to high value makes the solution unstable
    DX=1/NN  'step along the reactor

    INPUT "Max fixbed conc. of A                qa0(st=0.01)=";QA0
    INPUT "Max fixbed conc. of B              qb0(0.01-0.04)=";QB0
    INPUT "Max gas phase conc. of A and B   CA0=CB0(st=0.01)=";CA0
    INPUT "    Rate constant A                    Kva(st=80)=";KVA
    INPUT "    Rate constant B                    kvb(st=20)=";KVB


    PRINT""
     XX=QB0/QA0:CB0=CA0
     EPS=.3:  B=1:  TET=.15:RB=2000
     FOR I=1 TO NN: CA0(I)=0:CB0(I)=0:QA0(I)=0:QB0(I)=0: NEXT I
     QA1(0)=1:CA1(0)=1:QB1(0)=1:CB1(0)=1
     K3A=RB*QA0/CA0/EPS/B
     K3B=RB*QB0/CB0/EPS/B
     DT=K2*DX*TET*EPS               'time step
     K1A=DT*B*(1-EPS)*KVA*CA0/RB/QA0
     K1B=DT*B*(1-EPS)*KVB*CB0/RB/QB0

   FOR I=1 TO NN: CA1(I)=0:CB1(I)=0: NEXT I' GUESS OF THE ROOT VALUES AT START

   [REACTA]      REM  *** Reaction A ***
      CA5=CA1(NN-10)
      FOR K=1 TO NN  'SOLUTION OF FIRST PART DIFF EQ BY THE IMPLICIT METHOD
       QA1(K)=QA0(K)+K1A*CA0(K)*(1-QA0(K))  'ALGORITM FOR RATE EQ
       CA1(K)=(CA1(K-1)*K2+CA0(K)+K3A*QA0(K)-QA1(K)*K3A)/(1+K2)
      NEXT K
      IF ABS((CA5-CA1(NN-10))/CA1(NN-10))>.001 GOTO [REACTA]'  CONVERGENCE CHECK

   [REACTB]     REM *** Reaction B ***
      CB5=CB1(NN-10)
      FOR K=1 TO NN
       QB1(K)=QB0(K)+K1B*CB0(K)*(1-QB0(K))
       CB1(K)=(CB1(K-1)*K2+CB0(K)+K3B*QB0(K)-QB1(K)*K3B)/(1+K2)
      NEXT K
      IF ABS((CB5-CB1(NN-10))/CB1(NN-10))>.001 GOTO [REACTB]

      REM *** Uppgrading level ***
      FOR K=1 TO NN
       CA0(K)=(CA1(K)+XX*CB1(K))/(1+XX)
       QA0(K)=QA1(K):QB0(K)=QB1(K)
       CB0(K)=CA0(K)
      NEXT K

     REM **/ Print ***
     N=N+1
     N1=N1+1
     IF N1>=9.99/DT THEN [PRINT]
     GOTO [REACTA]

  [PRINT]
     PRINT USING("####.#",N*DT);"  ";
     PRINT USING("####.###",CA0(1));
     PRINT USING("####.###",CA0(NN));
     PRINT USING("####.###", QB1(NN));
     PRINT USING("####.###",QA1(NN))
     PRINT #1,USING("####.###",N*DT);"    "; CA0(NN)
     N1=0
     IF CA0(NN)>.99 THEN [CLOSE]
     GOTO [REACTA] 

  [CLOSE]
     CLOSE #1

     PRINT"  time     CA0(1) CA0(N)  QB1(N)  QA1(N)"
     PRINT "kva= "; KVA;"  kvb= ";KVB;"  qb0/qa0= ";XX
     PRINT" PRESS ENTER TO END"
     INPUT R$
     PRINT "   **** END ****"
     END