Countercurrent-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 . Countercurrent-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
0
by the following equation:
In[]:=
sol=Solve[y/(1-y)​​(α((xo-
p
l
/
p
h
y)/(1-xo-
p
l
/
p
h
(1-y)))),y]
Out[]=
y
1
18
19+90xo-
361-180xo+8100
2
xo
,y
1
18
19+90xo+
361-180xo+8100
2
xo

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
O
=810^5;
In[]:=
t=2.54 10^-3;
In[]:=
P'
B
=5010^-10;
In[]:=
eq1[a_]:=D[x[
A
m
],{
A
m
,1}]If[
A
m
0,Evaluate[(
p
l
P'
B
)/(
q
O
t)α(xo
p
h
/
p
l
-yi)(xo-yi)/yi/.t->2.5410^-3/.
p
h
->190/.
p
l
->19/.α->10/.
P'
B
->5010^-10/.
q
O
->810^5/.xoa],Evaluate[(
p
l
P'
B
)/(
q
O
t)(x[
A
m
]-y[
A
m
])/(y[
A
m
]-xo)((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
O
->810^5/.xoa]]
In[]:=
eq2[a_]:=D[y[
A
m
],{
A
m
,1}]If[
A
m
0,Evaluate[(xo-yi)
p
h
/
p
l
(α-yi(α-1))/((
q
O
t)/(
p
l
P'
B
)-((xo-yi)/((
p
l
P'
B
)/(
q
O
t)α(xo
p
h
/
p
l
-yi)(xo-yi)/yi))((α-1)(2yi-
p
h
/
p
l
xo-1)-
p
h
/
p
l
))/.t->2.5410^-3/.
p
h
->190/.
p
l
->19/.α->10/.
P'
B
->5010^-10/.
q
O
->810^5/.xoa],Evaluate[(
p
l
P'
B
)/(
q
O
t)(x[
A
m
]-y[
A
m
])/(x[
A
m
]-xo)((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
O
->810^5/.xoa]]
In[]:=
eq3[a_]:=y[0]yi/.xo->a
In[]:=
myODEsoln[Ω_,xo_]:=NDSolve[{eq1[xo],eq2[xo],eq3[xo],x[0]xo},{x[
A
m
],y[
A
m
]},{
A
m
,0,Ω}]
Use FindRoot to compute surface area of the separation module. Surface area and
x
O
must verify
x
f
=0.209andthematerialbalance.
Surface area sign must reversed.
In[]:=
Off[InterpolatingFunction::dmval]
In[]:=
Off[NDSolve::ndsz]
In[]:=
Off[General::spell1]
In[]:=
xend[Ω_?NumericQ,xo_?NumericQ]:=Flatten[(x[
A
m
]/.myODEsoln[Ω,xo])/.
A
m
Ω]​​yend[Ω_?NumericQ,xo_?NumericQ]:=Flatten[(y[
A
m
]/.myODEsoln[Ω,xo])/.
A
m
Ω]​​
P
sol
=FindRoot[{xend[Ω,xo]0.209,0.2yend[Ω,xo]+0.8xo0.209},{Ω,-5.10^8,-1.10^8},{xo,0.05,0.3},MaxIterations1000]
Out[]=
{Ω-2.85971×
8
10
,xo0.117169}
Surface area sign must reversed. Surface area is equal to 2.85971
8
10
2
cm
and reject out mole fraction
x
O
= 0.117169
In[]:=
P
sol
Out[]=
{Ω-2.85971×
8
10
,xo0.117169}
Checking the values of the mole fractions x and y at the entrance and exit of the module. Permeate mole fraction is 0.576323 and reject out mole fraction
x
O
= 0.117169.
In[]:=
x[
A
m
]/.myODEsoln[Ω/.
P
sol
,xo/.
P
sol
]/.
A
m
Ω/.
P
sol
Out[]=
{0.209}
In[]:=
x[
A
m
]/.myODEsoln[Ω/.
P
sol
,xo/.
P
sol
]/.
A
m
0
Out[]=
{0.117169}
In[]:=
y[
A
m
]/.myODEsoln[Ω/.
P
sol
,xo/.
P
sol
]/.
A
m
Ω/.
P
sol
Out[]=
{0.576323}
In[]:=
y[
A
m
]/.myODEsoln[Ω/.
P
sol
,xo/.
P
sol
]/.
A
m
0
Out[]=
{0.461436}
Checking Material Balance: θ
y
P
+(1-θ)
x
O
=
x
f
= 0.209
In[]:=
0.576323 0.2+0.117169 0.8
Out[]=
0.209
Plotting x and y versus surface area of the separation module.
In[]:=
First[y[
A
m
]/.myODEsoln[Ω/.
P
sol
,xo/.
P
sol
]]/.
A
m
area
Out[]=
InterpolatingFunction
Domain: -2.86×
8
10
,0.
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[]=