Cocurrent-Flow Model for Gas Separation by Membranes ​​Authors: Housam Binous and Ahmed Bellagi

reference : H . Binous, Gas Permeation Computations with Mathematica®, Chemical Engineering Education, 40 (2), pp . 140 - 144, 2006.
We compute the permeate and residue compositions obtained using a membrane separation module . Membrane area is also computed . Feed is air and is considered as a binary mixture . Cocurrent-flow mixing model is used . We use the shooting technique in order to determine the surface area of the separation module.
Reference:
​
C. J. Geankoplis, Transport Processes and Unit Operations, 3rd edition, Prentice Hall, 1993.
​
W. P. Walawender and S. A. Stern, Separation Science, 7 (5), pp. 553-584, 1972.
The separation factor for this binary mixture (components A and B) that we want to separate by membrane permeation is
*
α
=
P'
A
P'
B
=
10.
In[]:=
α=10;
The pressure of the feed side is
p
h
=190cmHg
while
p
l
=
19 cmHg is the pressure on the permeate side.
In[]:=
p
l
=19;
In[]:=
p
h
=190;
Stage cut is equal to θ=0.2. The material balance is given by the following equation:
In[]:=
θ=0.2;
Mole fraction at end of separation module is
y
i
.
It is related to
x
f
by the following equation:
In[]:=
sol=Solve[y/(1-y)​​(α((xf-
p
l
/
p
h
y)/(1-xf-
p
l
/
p
h
(1-y)))),y]
Out[]=
y
1
18
19+90xf-
361-180xf+8100
2
xf
,y
1
18
19+90xf+
361-180xf+8100
2
xf

In[]:=
yi=y/.sol[[1]];
t= 2.54
-3
10
cm is membrane thickness.
q
O
=(1-θ)
q
f
is reject out flow rate.
q
f
=1
6
10
3
cm
(STP)/s
is feed flow rate.
P'
B
=50
-10
10
3
cm
(STP)cm(s
2
cm
cmHg)
is the permeability of component B. Governing equations for mole fraction x and y in the reject and permeate sides are given by:
L'Hôpital rule gives the derivatives as
A
m
0.
In[]:=
q
f
=110^6;
In[]:=
t=2.54 10^-3;
In[]:=
P'
B
=5010^-10;
In[]:=
xf=0.209;
In[]:=
eq1=D[x[
A
m
],{
A
m
,1}]If[
A
m
0,Evaluate[(
p
l
P'
B
)/(
q
f
t)α(xf
p
h
/
p
l
-yi)(xf-yi)/yi/.t->2.5410^-3/.
p
h
->190/.
p
l
->19/.α->10/.
P'
B
->5010^-10/.
q
f
->810^5],Evaluate[(
p
l
P'
B
)/(
q
f
t)(x[
A
m
]-y[
A
m
])/(y[
A
m
]-xf)((1-x[
A
m
])α(
p
h
/
p
l
x[
A
m
]-y[
A
m
])-x[
A
m
](
p
h
/
p
l
(1-x[
A
m
])-(1-y[
A
m
])))/.t->2.5410^-3/.
p
h
->190/.
p
l
->19/.α->10/.
P'
B
->5010^-10/.
q
f
->810^5]];
In[]:=
eq2=D[y[
A
m
],{
A
m
,1}]If[
A
m
0,Evaluate[(xf-yi)
p
h
/
p
l
(α-yi(α-1))/((
q
f
t)/(
p
l
P'
B
)-((xf-yi)/((
p
l
P'
B
)/(
q
f
t)α(xf
p
h
/
p
l
-yi)(xf-yi)/yi))((α-1)(2yi-
p
h
/
p
l
xf-1)-
p
h
/
p
l
))/.t->2.5410^-3/.
p
h
->190/.
p
l
->19/.α->10/.
P'
B
->5010^-10/.
q
f
->810^5],Evaluate[(
p
l
P'
B
)/(
q
f
t)(x[
A
m
]-y[
A
m
])/(x[
A
m
]-xf)((1-y[
A
m
])α(
p
h
/
p
l
x[
A
m
]-y[
A
m
])-y[
A
m
](
p
h
/
p
l
(1-x[
A
m
])-(1-y[
A
m
])))/.t->2.5410^-3/.
p
h
->190/.
p
l
->19/.α->10/.
P'
B
->5010^-10/.
q
f
->810^5]];
In[]:=
eq3=y[0]yi;
In[]:=
myODEsoln[Ω_]:=NDSolve[{eq1,eq2,eq3,x[0]xf},{x[
A
m
],y[
A
m
]},{
A
m
,0,Ω}]
Use FindRoot to compute surface area of the separation module. Surface area must verify
thematerialbalance.
In[]:=
Off[InterpolatingFunction::dmval]
In[]:=
Off[NDSolve::mxst]
In[]:=
Off[General::spell1]
In[]:=
xend[Ω_?NumericQ]:=Flatten[(x[
A
m
]/.myODEsoln[Ω])/.
A
m
Ω]​​yend[Ω_?NumericQ]:=Flatten[(y[
A
m
]/.myODEsoln[Ω])/.
A
m
Ω]​​
P
sol
=FindRoot[0.2yend[Ω]+0.8xend[Ω]xf,{Ω,3.10^8,5.10^8},MaxIterations1000]
Out[]=
{Ω2.95525×
8
10
}
Surface area is equal to 2.95525
8
10
2
cm
.
In[]:=
P
sol
Out[]=
{Ω2.95525×
8
10
}
Checking the values of the mole fractions x and y at the entrance and exit of the module. Permeate mole fraction is 0.558456 and reject out mole fraction
x
O
= 0.121636.
In[]:=
x[
A
m
]/.myODEsoln[Ω/.
P
sol
]/.
A
m
Ω/.
P
sol
Out[]=
{0.121636}
In[]:=
x[
A
m
]/.myODEsoln[Ω/.
P
sol
]/.
A
m
0
Out[]=
{0.209}
In[]:=
y[
A
m
]/.myODEsoln[Ω/.
P
sol
]/.
A
m
Ω/.
P
sol
Out[]=
{0.558456}
In[]:=
y[
A
m
]/.myODEsoln[Ω/.
P
sol
]/.
A
m
0
Out[]=
{0.654834}
Checking Material Balance: θ
y
P
+(1-θ)
x
O
=
x
f
= 0.209
In[]:=
0.558456 0.2+0.121636 0.8
Out[]=
0.209
Plotting x and y versus surface area of the separation module.
In[]:=
First[y[
A
m
]/.myODEsoln[Ω/.
P
sol
]]/.
A
m
area
Out[]=
InterpolatingFunction
Domain: 0.,2.96×
8
10

Output: scalar
[area]
In[]:=
Plot[%,{area,0.001,Ω/.
P
sol
},PlotStyleRGBColor[1,0,0],Frame->True,GridLines->Automatic,AspectRatio->1,FrameLabel->{Style["Area / cm2",14],Style["permeate composition",14]}]
Out[]=