Chapter 7 - Tubular Reactor with Axial Dispersion

Exa 7.1 Page 122

from __future__ import division
from numpy import zeros
e1=1 ; e2=1                                                      #assumed values
u=1 ; D=10**-4 ; k=1 ; C_a_in=1 ; delta_x=10/50                         #given data
cf_ca1_n1=-2*D/delta_x**2-3*u/delta_x-k-2*u**2/D                 #co-efficient of C-A1 at node 1
cf_ca2_n1=2*D/delta_x**2+u/delta_x
cf_da1_n1=-(2*u**2/D+2*u/delta_x)*C_a_in                         #right hand side co-efficient
cf_ca1_n2=D/delta_x**2+u/delta_x
cf_ca2_n2=-2*D/delta_x**2-u/delta_x-k
cf_ca3_n2=D/delta_x**2
cf_da1_n2=0
cf_ca2_n3=cf_ca1_n2
cf_ca3_n3=cf_ca2_n2
cf_ca4_n3=cf_ca3_n2
cf_da1_n3=0
cf_ca50_n51=2*D/delta_x**2+u/delta_x                           #co-efficient of C-A50 at node 51
cf_ca51_n51=-2*D/delta_x**2-u/delta_x-k
cf_da51_n51=0
a=[0]
for i in range(2,51):
    a.append(cf_ca1_n2)
    
a[49]=cf_ca2_n1
c=[cf_ca2_n1]
for i in range(2,51):
    c.append(cf_ca3_n2)

d=[cf_da1_n1]
for i in range(2,52):
    d.append(cf_da1_n2)
x=[]
for i in range(1,52):
    x.append(0)

b=[cf_ca1_n1]
x1=[]
for i in range(2,52):
    b.append(cf_ca2_n2)
while e1>1e-6 and e2>1e-6:
    for i in range(1,52):
        x1.append(x[i-1])

    i=1 ;  n=51 ; Beta=[b[i-1]]
    Gamma=[d[i-1]/Beta[i-1]]
    i1=i+1,
    for j in range(2,52):
        Beta.append(b[j-1-1]-a[j-1-1]*c[j-2-1]/Beta[j-2-1])
        Gamma.append((d[j-1]-a[j-1-1]*Gamma[j-2-1])/Beta[j-1-1])
    x=zeros(n)
    x[-1]=Gamma[n-1]
    n1=n-i
    for k in range(1,n1+1):
        j=n-k ; x[j-1]=Gamma[j-1]-c[j-1]*x[j]/Beta[j-1]
    
    e1=abs(x[(41)]-x1[(41)]); e2=abs(x[(17)]-x1[(17)])

print "First Order Reaction-50 parts - "
for i in range(1,52):
    print x[i-1]

First Order Reaction-50 parts - 
0.999699990083
3.12062838745e-11
0.000249850053626
8.66524048103e-08
0.000208139054698
7.22363767651e-08
0.000173451600052
6.01978238317e-08
0.0001445450073
5.01655498325e-08
0.000120455845492
4.18052053996e-08
0.000100381265215
3.48381549557e-08
8.36522160046e-05
2.90321989598e-08
6.97111480663e-05
2.41938351073e-08
5.80934301186e-05
2.01618092385e-08
4.84118640498e-05
1.68017410206e-08
4.03437802862e-05
1.40016452881e-08
3.36202837822e-05
1.16681997736e-08
2.80172922214e-05
9.72363484122e-09
2.33480677471e-05
8.10314156079e-09
1.94569933173e-05
6.75271174067e-09
1.62143862631e-05
5.62733792942e-09
1.35121761931e-05
4.6895133967e-09
1.12603031968e-05
3.90798209982e-09
9.38371630679e-06
3.25669697484e-09
7.81987218174e-06
2.71395183372e-09
6.51665064666e-06
2.26165793522e-09
5.43061761928e-06
1.88474111898e-09
4.52557752837e-06
1.57063941036e-09
3.77136697934e-06
1.30953827722e-09
3.14441964484e-06

Exa 7.2 Page 131

from __future__ import division
from numpy import zeros
e1=1 ; e2=1
u=1 ; D=10**-4 ; k=1 ; C_a_in=1 ; delta_x=10/20
cff_ca2_n1=2*D/delta_x**2+u/delta_x                       #co-efficient of C-A2 at node 1
cff_da1_n1=-(2*u**2/D+2*u/delta_x)*C_a_in                 #right hand side co-efficient
cff_ca1_n2=D/delta_x**2+u/delta_x
cf_ca1_n2=D/delta_x**2+u/delta_x
cff_ca3_n2=D/delta_x**2                                    #co-efficient of C-A3 at node 2
cf_ca3_n2=D/delta_x**2
cff_da1_n2=0
cff_ca2_n3=cf_ca1_n2
cff_ca4_n3=cf_ca3_n2
cff_da1_n3=0
cff_ca20_n21=2*D/delta_x**2+u/delta_x                        #co-efficient of C-A20 at node 21
cff_da21_n21=0
a=[0]

for i in range(2,21):
    a.append(cff_ca1_n2)

a.append(cff_ca2_n1) ; c=[cff_ca2_n1]
c=[0]
for i in range(2,21):
    c.append(cff_ca3_n2)

d=[cff_da1_n1]
for i in range(2,22):
    d.append(cff_da1_n1)
       
x=[] ; x1=[]
for i in range(1,22):
    x.append(0)    

while e1>1e-6 and e2>1e-6:
    for i in range(1,22):
        x1.append(x[i-1])
    cff_ca1_n1=-2*D/delta_x**2-3*u/delta_x-x1[0]-2*u**2/D    #             //main diagonal elements dependence on conc.
    b=[cff_ca1_n1]
    for i in range(2,22):
        b.append(-2*D/delta_x**2-u/delta_x-x[i-1])

    #solving by TDMA method
    i=1 ;  n=21 ; Beta=[b[i-1]]
    Gamma=[d[i-1]/Beta[i-1]]
    i1=i+1
    for j in range(i1,n+1):
        Beta.append((b[j-1]-a[j-1]*c[j-2])/Beta[j-2])
        Gamma.append((d[j-1]-a[j-1]*Gamma[j-2])/Beta[j-1])
    x=zeros(n)
    x[-1]=Gamma[n-1]
    n1=n-i
    for k in range(1,n1+1):
        j=n-k
        x[j-1]=Gamma[j-1]-c[j-1]*x[j]/Beta[j-1]
    
    e1=abs(x[0]-x1[0]) ;  e2=abs(x[20]-x1[20])

print "Second Order Reaction-20 parts"    
for i in range(1,22):
    print x[i-1]
    

Second Order Reaction-20 parts
0.999899990007
-99989995.0013
-7.99620241884
-499875443627.0
-3.99760071124
-2.50037661937e+19
-3.99784153685
-1.25099285374e+31
-3.99749207478
-6.26163384772e+46
-3.99765851293
-3.13583861982e+66
-3.99757331924
-1.57118962175e+90
-3.9976154302
-7.87634360632e+117
-3.99759424268
-3.95035153098e+149
-3.99760512762
-1.98227988436e+185
-3.98527206101e-36