Assignment – 5

Gas absorption again!

This time, again it was gas absorption but with some changes. This time we had to perform the gas absorption calculations with taking carier gas into account.

This time it was a different story, I had my Chemical Engineering Laboratory yesterday and we had a surprise viva based on our experiments and I did well. After the lab, I was tired but had a meeting with the student editor of our institute’s journal (The Bombay Technologist). We had a cup of Tea and that refreshed me up again. I reached home late as usual, had food and it was already midnight. I started the assignment and

As usual I did this at 3 AM!

https://github.com/chiragmule1729/3-AM-Assignments.git

But this time I tried something different!

Reference –

Problem Statement

Seperation of C3 and C4 from lower alkane hydrocarbons

Formulas Used

ln(x)=A+BT+ClnT+DTln(x)=A+BT+ClnT+DT H=1xH=1x

  • T in K
  • H in atm
  • x= mole fraction of solute in water at soltue partial pressure of 1 atm

In [204]:

import scipy
import scipy.integrate as scint
import scipy.optimize as scopt
import matplotlib.pyplot as plt
import pandas as pd
pi=scipy.pi
e=scipy.e

Data

In [205]:

T=25+273.16 #in K
P=101325 #in Pa
#1=methane,#2=Water,#3=propane,#4=butane
#Vapour pressure formula
#ln(P)=c1+c2/T+c3lnT+c4*T^c5
c1=[73.649,83.107,76.945]
c2=[-7258.2,-6486.2,-6729.8]
c3=[-7.3037,-9.2194,-8.179]
c4=[4.1653e-6,6.9844e-6,5.3017e-6]
c5=[2,2,2]
R=8.314472
A=-181.587
B=8632.13
C=24.7981
D=0
H=e**(-1*(A+B/T+C*scipy.log(T)+D*T))*101325
G=30 #mol/s
L=40 #mol/s
Ng1in=100 #mol/s
Ng2in=0
Ng3in=12
Ng4in=15
Nl1in=0
Nl2in=40
Nl3in=0
Nl4in=0
Di=0.5 #column diamter in m
Qg=(Ng1in*18+Ng2in*18+Ng3in*44+Ng4in*58)/(1.184*1e3)
ug=Qg/(0.25*pi*Di**2)
ka=0.0269*ug**0.82
Z=10

In [206]:

def hc(b):
    return e**(c1[b]+c2[b]/T+c3[b]*scipy.log(T)+c4[b]*T**c5[b])

In [207]:

def model(SV,z,obj):
    e=obj.e
    Di=obj.Di
    L=obj.L
    G=obj.G
    P=obj.P
    T=obj.T
    H=obj.H
    Z=obj.Z
    A=obj.A
    B=obj.B
    C=obj.C
    D=obj.D
    [Ng1,Ng2,Ng3,Ng4,Nl1,Nl2,Nl3,Nl4]=SV
    dNg1bydz=-ka/(R*T)*(P*Ng1/(Ng1+Ng2+Ng3+Ng4)-H*Nl1/(Nl1+Nl2+Nl3+Nl4))
    dNg2bydz=-ka/(R*T)*(P*Ng2/(Ng1+Ng2+Ng3+Ng4)-hc(0))
    dNg3bydz=-ka/(R*T)*(P*Ng3/(Ng1+Ng2+Ng3+Ng4)-hc(1))
    dNg4bydz=-ka/(R*T)*(P*Ng4/(Ng1+Ng2+Ng3+Ng4)-hc(2))
    dNl1bydz=ka/(R*T)*(P*Ng1/(Ng1+Ng2+Ng3+Ng4)-H*Nl1/(Nl1+Nl2+Nl3+Nl4))
    dNl2bydz=-ka/(R*T)*(P*Ng2/(Ng1+Ng2+Ng3+Ng4)-hc(0))
    dNl3bydz=-ka/(R*T)*(P*Ng3/(Ng1+Ng2+Ng3+Ng4)-hc(1))
    dNl4bydz=-ka/(R*T)*(P*Ng4/(Ng1+Ng2+Ng3+Ng4)-hc(2))
    return [dNg1bydz,dNg2bydz,dNg3bydz,dNg4bydz,dNl1bydz,dNl2bydz,dNl3bydz,dNl4bydz]

In [208]:

def get_boundary_residuals(SV0,obj):
    z=scipy.array([0,obj.Z])
    trial_solution = scint.odeint(model, SV0, z, args = (obj,))
    SVi = trial_solution[0]
    SVb = trial_solution[-1]
    error = []
    k = 0
    for ic in obj.initial_conditions:
        error.append(ic - SVi[k])
        k += 1
    for bc in obj.boundary_conditions:
        error.append(bc - SVb[k])
        k += 1
    return error

Applying Boundary Value problem

In [209]:

class absorption:
    def __init__(self):
        self.e=2.718281828
        self.L=5
        self.G=5
        self.Di=0.5
        self.Z=10
        self.T=25+273.16
        self.P=101325
        self.A=-181.587
        self.B=8632.13 
        self.C=24.7981
        self.D=0
        self.H=e**(-1*(A+B/T+C*scipy.log(T)+D*T))*101325
        self.Ng1in=30
        self.Ng2in=0
        self.Ng3in=12
        self.Ng4in=15
        self.Nl1in=0
        self.Nl2in=40
        self.Nl3in=0
        self.Nl4in=0
    def solve(self):
        self.x=scipy.linspace(0,Z,10)
        self.initial_conditions=[self.Ng1in,self.Ng2in,self.Ng3in,self.Ng4in]
        self.boundary_conditions=[self.Nl1in,self.Nl2in,self.Nl3in,self.Nl4in]
        SV0=self.initial_conditions+self.boundary_conditions
        shooting_method_solution = scopt.least_squares(get_boundary_residuals, SV0,args = (self,))
        SV0=shooting_method_solution.x
        self.solution = scint.odeint(model, SV0,self.x,args = (self,))
        self.solution=pd.DataFrame({"h(m)":self.x,
                "Ng1":self.solution[:,0],
                "Ng2":self.solution[:,1],
                "Ng3":self.solution[:,2],
                "Ng4":self.solution[:,3],
                "Nl1":self.solution[:,4],
                "Nl2":self.solution[:,5],
                "Nl3":self.solution[:,6],
                "Nl4":self.solution[:,7],})
    def plot(self):
        fig = plt.figure()
        ay = fig.add_subplot(121)
        ay.plot(self.solution["h(m)"],self.solution["Ng1"],'r',label="Methane")
        ay.plot(self.solution["h(m)"],self.solution["Ng2"],'g',label="Water")
        ay.plot(self.solution["h(m)"],self.solution["Ng3"],'b',label="Propane")
        ay.plot(self.solution["h(m)"],self.solution["Ng4"],'y',label="Butane")
        ay.legend()
        ay.title.set_text("Flow rates in the gas phase")
        ay.xaxis.label.set_text("Height")
        ay.yaxis.label.set_text("moles per second")
        ax = fig.add_subplot(122)
        ax.plot(self.solution["h(m)"],self.solution["Nl1"],'r',label="Methane")
        ax.plot(self.solution["h(m)"],self.solution["Nl2"],'g',label="Water")
        ax.plot(self.solution["h(m)"],self.solution["Nl3"],'b',label="Propane")
        ax.plot(self.solution["h(m)"],self.solution["Nl4"],'y',label="Butane")
        ax.legend()
        ax.title.set_text("Flow rates in the liquid phase")
        ax.xaxis.label.set_text("Height")

In [210]:

ccgas=absorption()
ccgas.Z=10

In [211]:

ccgas.solve()

Flowsheet

In [212]:

ccgas.solution

Out[212]:

h(m)Ng1Ng2Ng3Ng4Nl1Nl2Nl3Nl4
00.00000030.0000008.069046e-1012.00000015.000000-0.13288338.7910295.579589e+001.154830e+01
11.11111129.8667702.988101e-0111.12325912.7736070.00034739.0898394.702849e+009.321910e+00
22.22222229.8667705.424077e-0110.30477410.8289140.00034739.3334373.884363e+007.377217e+00
33.33333329.8667707.366282e-019.5496219.1501910.00034639.5276573.129210e+005.698494e+00
44.44444429.8667718.876953e-018.8612867.7182510.00034639.6787242.440876e+004.266554e+00
55.55555629.8667711.001966e+008.2414716.5111690.00034639.7929951.821061e+003.059472e+00
66.66666729.8667711.085668e+007.6900425.5052200.00034639.8766971.269632e+002.053523e+00
77.77777829.8667711.144655e+007.2051404.6759310.00034539.9356847.847296e-011.224234e+00
88.88888929.8667721.184220e+006.7834253.9991020.00034539.9752493.630142e-015.474054e-01
910.00000029.8667721.208971e+006.4204113.4516970.00034540.000000-6.502894e-09-9.126016e-09

In [213]:

ccgas.plot()

Assignment – 4

Multicomponent Absorption

Multicomponent gas absorption

What exactly are we doing?

We are solving for the initial value problem There is a gas mixture consists of “Acetaldehyde” & “Acetone”.The mixture is scrubbed using pure water. (Heat balances are ignored). Pure water comes from the top & gases come from bottom. No reaction is taking place.

Formulae used

https://github.com/chiragmule1729/3-AM-Assignments.git

How are we solving?

We are solving for the initial value boundary residuals. Residual goes to zero.In [1]:

import scipy
import scipy.integrate as scint
import scipy.optimize as scopt
import matplotlib.pyplot as plt
import pandas as pd
from scipy.misc import derivative

In [2]:

Ls=200 #molar water flow rate
Gi=150 #molar gas flow rate input
D=0.5
y1i=0.05 #mole fractions
y2i=0.05 #mole fractions
Gs=(1-y1i-y2i)*Gi
M=Gs/Gi*28.81+y1i*30+y2i*28
Qg=Gi*M*1e-3/1.225
S=scipy.pi*0.25*D**2
ug=Qg/S
kLa=0.0269*ug**0.82
kGa=0.0005
H1=40.54
H2=20.13
H=3

In [3]:

#finding x1i
def findx1i(x1,y1,x1i):
    x1i=(y1+kLa/kGa*x1)/(kLa/kGa+H1)
    return x1i

In [4]:

#finding x2i
def findx2i(x2,y2,x2i):
    x2i=(y2+kLa/kGa*x2)/(kLa/kGa+H2)
    return x2i

In [5]:

def model(SV,z,obj):
    [x1,x2,y1,y2]=SV
    D=obj.D
    S=obj.S
    Ls=obj.Ls
    H1=obj.H1
    H2=obj.H2
    H=obj.H
    x1i=0.005
    x2i=0.005
    dx1bydz=kLa*S/Ls*(1-x1)**2*scipy.log((1-x1)/(1-findx1i(x1,y1,x1i)))
    dx2bydz=kLa*S/Ls*(1-x2)**2*scipy.log((1-x2)/(1-findx2i(x2,y2,x2i)))
    dy1bydz=-kGa*S/Gs*(1-y1)**2*scipy.log((1-y1)/(1-H1*x1i))
    dy2bydz=-kGa*S/Gs*(1-y2)**2*scipy.log((1-y2)/(1-H2*x2i))
    return [dx1bydz,dx2bydz,dy1bydz,dy2bydz]

In [6]:

def get_boundary_residuals(SV0, obj):
    z = scipy.array([0, obj.H])
    trial_solution = scint.odeint(model, SV0, z, args = (obj,))
    SVi = trial_solution[0]
    SVb = trial_solution[-1]
    error = []
    k = 0
    for ic in obj.initial_conditions:
        error.append(ic - SVi[k])
        k += 1
    for bc in obj.boundary_conditions:
        error.append(bc - SVb[k])
        k += 1
    return error

In [7]:

class absorption:
    def __init__(self):
        self.D=0.5
        self.S=scipy.pi*0.25*D**2
        self.Ls=250
        self.H1=299
        self.H2=354
        self.H=3
        self.y1in=0.04
        self.y2in=0.04
        self.x1in=0
        self.x2in=0
    def solve(self):
        self.x=scipy.linspace(0,self.H,10)
        self.initial_conditions=[self.y1in,self.y2in]
        self.boundary_conditions=[self.x1in,self.x2in]
        SV0=self.initial_conditions + self.boundary_conditions
        shooting_method_solution = scopt.least_squares(get_boundary_residuals,
                                                      SV0,
                                                      args = (self,))
        SV0=shooting_method_solution.x
        solution = scint.odeint(model, SV0, self.x, args = (self,))
        self.solution=pd.DataFrame({
            "x":self.x,
            "x1":solution[:,0],
            "x2":solution[:,1],
            "y1":solution[:,2],
            "y2":solution[:,3],
            
        }
        
        )
    def plotx(self):
        fig = plt.figure()
        ax = fig.add_subplot(111)
        ax.plot(self.solution["x"], 
                self.solution["x1"],
                'c', label = "Acetaldehyde")
        ax.plot(self.solution["x"],
               self.solution["x2"],
               'g', label="Acetone")
        ax.legend()
        ax.xaxis.label.set_text("Height in m")
        ax.yaxis.label.set_text("mole fraction in gas")
        ax.title.set_text("Concentration Profiles")
    def ploty(self):
        fig = plt.figure()
        ax = fig.add_subplot(111)
        ax.plot(self.solution["x"], 
                self.solution["y1"],
                'c', label = "Acetaldehyde")
        ax.plot(self.solution["x"],
               self.solution["y2"],
               'g', label="Acetone")
        ax.legend()
        ax.xaxis.label.set_text("Height in m")
        ax.yaxis.label.set_text("mole fraction in liquid")
        ax.title.set_text("Concentration Profiles")

In [8]:

cchex1 = absorption()
cchex1.H = 5

In [12]:

cchex1.solve()

In [10]:

cchex1.solution

Out[10]:

xx1x2y1y2
00.0000000.0400000.0400002.556887e-069.503465e-07
10.5555560.0400000.0400002.272789e-068.447526e-07
21.1111110.0399990.0400001.988692e-067.391587e-07
31.6666670.0399990.0400001.704593e-066.335648e-07
42.2222220.0399990.0399991.420495e-065.279708e-07
52.7777780.0399980.0399991.136397e-064.223768e-07
63.3333330.0399980.0399998.522978e-073.167826e-07
73.8888890.0399980.0399995.681987e-072.111885e-07
84.4444440.0399970.0399992.840993e-071.055942e-07
95.0000000.0399970.039999-4.116221e-13-9.680519e-14

In [11]:

cchex1.plotx(),cchex1.ploty()

Out[11]:

(None, None)

In [ ]:

Assignment – 3

Copper Disc Dissolution

Copper Disc Dissolution

In the Chemical Engineerig Laboratory, last week I performed the experiment of “Copper Disc Dissolution”. We also had to do the calculations at home for this experiment. In the same week, we got to learn ‘Pandas’ and as it is a good tool for computation, I used it to perform the calculations of my experiment.

This week was really really very hectic for me, we were having Funtech which is our intra college festival and being a part of Technological Association, I was a part of organising team. During last 7 days, I also had 2 company interviews for in plant training and there were thousand more things which I don’t even remember now. For the last 6 days, I was literally running all the time, used to go home very late (by 12.00-12.30 in the night) and I used to leave again in the early morning by 6.45 am. But, whatever I was doing, I was enjoying each and every part of it and there was a thrill in that. Last night was the final night and it was one of the best nights of the college, the fest ended well with a bang! We all went to have dinner together and then I took a train to go home. I took a quick nap in the travel 😋 because I knew after going home I had this assignment still pending of which today was the submission and I will have to stay up and do it othwrwise I will lose marks.

And as usual, I did this at 3 AM !

https://github.com/chiragmule1729/3-AM-Assignments.git

In [16]:

import scipy
import math
import pandas as pd
import matplotlib.pyplot as plt

Observation Table

In [48]:

example = pd.DataFrame(
{
    "Speed":[0,354.3,425.5,777.8,999.9,1183.4],
    "TitrationI":[31.2,31.2,37.9,42.1,48.4,50.0],
    "TitrationF":[31.2,37.9,42.1,48.4,50.0,52.1],
    "Disc diameter":[4.532,4.530,4.525,4.510,4.500,4.492],
    "Thickness":[0.63,0.608,0.60,0.57,0.551,0.541],   
},
    index = ["initial",1,2,3,4,5]
)
example

Out[48]:

SpeedTitrationITitrationFDisc diameterThickness
initial0.031.231.24.5320.630
1354.331.237.94.5300.608
2425.537.942.14.5250.600
3777.842.148.44.5100.570
4999.948.450.04.5000.551
51183.450.052.14.4920.541

Calculation

In [55]:

def Bi(row,a,b,c):
    ti=row["TitrationI"]
    bi=(a*(b-ti))/c
    return bi

In [56]:

example["Bi"]=example.apply(Bi,axis=1, args=(0.1,53.5,10))

In [57]:

example

Out[57]:

SpeedTitrationITitrationFDisc diameterThicknessBiBf
initial0.031.231.24.5320.6300.2230.223
1354.331.237.94.5300.6080.2230.156
2425.537.942.14.5250.6000.1560.114
3777.842.148.44.5100.5700.1140.051
4999.948.450.04.5000.5510.0510.035
51183.450.052.14.4920.5410.0350.014

In [58]:

def Bf(row,a,b,c):
    tf=row["TitrationF"]
    bf=(a*(b-tf))/c
    return bf

In [59]:

example["Bf"]=example.apply(Bf,axis=1, args=(0.1,53.5,10))

In [60]:

example

Out[60]:

SpeedTitrationITitrationFDisc diameterThicknessBiBf
initial0.031.231.24.5320.6300.2230.223
1354.331.237.94.5300.6080.2230.156
2425.537.942.14.5250.6000.1560.114
3777.842.148.44.5100.5700.1140.051
4999.948.450.04.5000.5510.0510.035
51183.450.052.14.4920.5410.0350.014

In [61]:

def ksle(row):
    bi=row["Bi"]
    bf=row["Bf"]
    ksle=248.139*math.log(bi/bf)
    return ksle

In [62]:

example["ksle"]=example.apply(ksle,axis=1)

In [63]:

example

Out[63]:

SpeedTitrationITitrationFDisc diameterThicknessBiBfksle
initial0.031.231.24.5320.6300.2230.2230.000000
1354.331.237.94.5300.6080.2230.15688.663976
2425.537.942.14.5250.6000.1560.11477.830673
3777.842.148.44.5100.5700.1140.051199.596266
4999.948.450.04.5000.5510.0510.03593.418768
51183.450.052.14.4920.5410.0350.014227.367466

In [68]:

def Re(row,a,b):
    d=row["Disc diameter"]
    n=row["Speed"]
    Re=(d*d*n*a)/(100*60*b)
    return Re

In [69]:

example["Re"]=example.apply(Re,axis=1, args=(1000,0.00119))

In [70]:

example

Out[70]:

SpeedTitrationITitrationFDisc diameterThicknessBiBfksleRe
initial0.031.231.24.5320.6300.2230.2230.0000000.000000e+00
1354.331.237.94.5300.6080.2230.15688.6639761.018285e+06
2425.537.942.14.5250.6000.1560.11477.8306731.220221e+06
3777.842.148.44.5100.5700.1140.051199.5962662.215760e+06
4999.948.450.04.5000.5510.0510.03593.4187682.835851e+06
51183.450.052.14.4920.5410.0350.014227.3674663.344359e+06

In [79]:

def Sh(row):
    Re=row["Re"]
    Sh=(0.03387*0.62*(Re**0.5))/(0.565)
    return Sh

In [80]:

example["Sh"]=example.apply(Sh,axis=1)

In [81]:

example

Out[81]:

SpeedTitrationITitrationFDisc diameterThicknessBiBfksleReSh
initial0.031.231.24.5320.6300.2230.2230.0000000.000000e+000.000000
1354.331.237.94.5300.6080.2230.15688.6639761.018285e+0637.505340
2425.537.942.14.5250.6000.1560.11477.8306731.220221e+0641.056101
3777.842.148.44.5100.5700.1140.051199.5962662.215760e+0655.324799
4999.948.450.04.5000.5510.0510.03593.4187682.835851e+0662.589305
51183.450.052.14.4920.5410.0350.014227.3674663.344359e+0667.969623

In [88]:

def kslp(row):
    d=row["Disc diameter"]
    Sh=row["Sh"]
    kslp=(0.0907)*Sh*100/d
    return kslp

In [89]:

example["kslp"]=example.apply(kslp,axis=1)

Calculation Table

In [90]:

example

Out[90]:

SpeedTitrationITitrationFDisc diameterThicknessBiBfksleReShkslp
initial0.031.231.24.5320.6300.2230.2230.0000000.000000e+000.0000000.000000
1354.331.237.94.5300.6080.2230.15688.6639761.018285e+0637.50534075.093474
2425.537.942.14.5250.6000.1560.11477.8306731.220221e+0641.05610182.293665
3777.842.148.44.5100.5700.1140.051199.5962662.215760e+0655.324799111.262955
4999.948.450.04.5000.5510.0510.03593.4187682.835851e+0662.589305126.152222
51183.450.052.14.4920.5410.0350.014227.3674663.344359e+0667.969623137.240535

In [108]:

n=[0,354.3,425.5,777.8,999.9,1183.4]

In [109]:

ksle=[88.663976,77.830673,199.596266,93.418768,227.367466]

In [110]:

kslp=[0,75.093474,82.293665,111.262955,126.152222,137.240535]

Graph

In [112]:

plt.plot(n,kslp,'g-',label='kslp vs rpm')
plt.ylabel('kslp')
plt.xlabel('rpm')
plt.legend(loc='best')
plt.show()

In [116]:

example.to_excel("Copper Disc Dissolution.xlsx", sheet_name="17CHE123")

Results and Conclusion

From the above graph, we can observe –

  • Predicted Mass transfer coefficient and experimental coefficient fairly match
  • Rate of mass transfer increases with increase in RPM
  • The reaction is phase-4 regime

Assignment -2

CONNECTING TANKS SYSTEM ANALYSIS

In the Chemical Process Control course, I recently studied the concept of “Connecting Tanks System”. In the class, we were told to solve this analytically.It was a simple system but for complicated systems, it will get difficult to solve the ODEs analytically. At the same time, I was having a course on python and my first assignment was solving any ODE using python. I learnt the basic syntax and solved a general ODE just when it striked in my mind that I can solve the Connecting Tank system very easily with python. So, I tried this at 3 AM !

What’s the question

There are two tanks (Tank 1 and Tank 2). water comes from source and enters Tank 1, goes to Tank 2 and then goes out from Tank 2. R1, R2 are resistances of the pipes. Let output flow be proportional to height of the water for both the tanks. We want to know how the system will respond on the step change in the input flow of water.In [1]:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

In [2]:

# Defining the function
def control(z,t,Fi):
    R1=1 #Resistance of outlet of tank 1
    R2=2 #Resistance of outlet of tank 2
    T1=5 #Tau for tank 1, A=5
    T2=10 #Tau for tank 2, A=5
    h1 = z[0] #deiation Level in tank 1 wrt time
    h2 = z[1] #deviation level in tank 2 wrt time
    dh1bydt = (-h1 + Fi*R1)/T1
    dh2bydt = (-h2 + h1*(R2/R1))/T2
    dzdt = [dh1bydt,dh2bydt]
    return dzdt

In [3]:

# initial condition at t=0
z0 = [0,0] #Because h1 & h2 are deviation variables

In [4]:

# time points
t = np.linspace(0,200,400)

In [5]:

# step change in input
Fi = np.zeros(400)
# change = 5.0
Fi[51:] = 5.0

In [6]:

# store solution
x = np.empty_like(t)
y = np.empty_like(t)
# record initial conditions
x[0] = z0[0]
y[0] = z0[1]

In [7]:

# solve ODE
for i in range(1,400):
    # span for next time step
    tspan = [t[i-1],t[i]]
    # solve for next step
    z = odeint(control,z0,tspan,args=(Fi[i],))
    # store solution for plotting
    x[i] = z[1][0]
    y[i] = z[1][1]
    # next initial condition
    z0 = z[1]

In [8]:

# plot results
plt.plot(t,Fi,'g:',label='input')
plt.plot(t,x,'b-',label='Tank 1 output')
plt.plot(t,y,'r--',label='Tank 2 output')
plt.ylabel('variables')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()

To refer the live jupyter file visit the link https://github.com/chiragmule1729/3-AM-Assignments.githttps://github.com/chiragmule1729/3-AM-Assignments.git

Design a site like this with WordPress.com
Get started