Nitrogen recovery for various membrane selectivities​
​Authors: Housam Binous and Ahmed Bellagi

Economical nitrogen production process requires high oxygen/nitrogen selectivity. The effect of selectivity on the efficiency of nitrogen recovery is depicted in Figure 8.25 of Ref. 1. Pressure ratio is taken equal to 10.​​
​​We reproduce Fig. 8.25 of Ref. 1 using a cross-flow gas permeation module and the methodology described in detail in Ref. 2.​​References:​[1] Richard W. Baker, Membrane Technology And Applications, Second Edition, NY: John Wiley & Sons, 2004. [2] H . Binous, U . Zahid and A . Bellagi, Modeling for Multi - component Gas Permeation with Orthogonal Collocation and Arc - Length Continuation, Chemistry Africa (2024) . https://doi .org/10.1007/s42250-024-01167-y
In[]:=
eq1i=
y
PAi
/(
y
PBi
)qA(xA[Np+1]-
p
l
/
p
h
y
PAi
)/(qB(xB[Np+1]-
p
l
/
p
h
y
PBi
));​​eq2i=
y
PAi
+
y
PBi
1;​​xAin=0.79;​​qA=1;qB=alpha;​​
q
f
=1.0;​​
p
l
=1/10
p
h
;​​Np=20;​​d[0,0]=(2Np^2+1)/6//N;​​d[Np,Np]=-(2Np^2+1)/6//N;​​yc[i_]:=Cos[iPi/Np]//N;​​d[i_,i_]:=-yc[i]/(2(1-yc[i]^2))//N;​​​​d[j_,k_]:=If[j≠k&&k≠0&&k≠Np&&j≠0&&j≠Np,(-1)^(j+k)/(yc[j]-yc[k])]//N;​​d[Np,0]=-1/2(-1)^Np;​​d[0,Np]=1/2(-1)^Np;​​d[j_,0]:=If[j≠Np&&j≠0,-1/2(-1)^j/(1-yc[j])];​​d[j_,Np]:=If[j≠Np&&j≠0,1/2(-1)^(j+Np)/(1+yc[j])];​​d[0,j_]:=If[j≠Np&&j≠0,2(-1)^j/(1-yc[j])];​​d[Np,j_]:=If[j≠Np&&j≠0,-2(-1)^(j+Np)/(1+yc[j])];​​Add=Table[d[i,j]//N,{i,0,Np},{j,0,Np}];​​xAp[i_]:=2/st
Np+1
∑
j=1
Add[[i,j]]xA[j];​​yAp[i_]:=2/st
Np+1
∑
j=1
Add[[i,j]]yA[j]​​qp[i_]:=2/st
Np+1
∑
j=1
Add[[i,j]]q[j];​​g=1/10.;​​eq1[i_]:=-yA[i]qp[i]qA(xA[i]-
p
l
/
p
h
yA[i]);​​eq2[i_]:=-(1-yA[i])qp[i]qB((1-xA[i])-
p
l
/
p
h
(1-yA[i]));​​eq3[i_]:=yA[i]qp[i]q[i]xAp[i]+qp[i]xA[i];​​yB[i_]:=1-yA[i]​​xB[i_]:=1-xA[i]​​xA[Np+1]=xAin;yA[Np+1]=
y
PAi
;q[Np+1]=
q
f
;q[1]=
q
f
(1-theta);​​unknown=Join[Flatten[Join[Table[{{xA[i],0.8}},{i,1,Np}],Table[{{yA[i],0.2}},{i,1,Np}],Table[{{q[i],0.99
q
f
}},{i,2,Np}]],1],{{theta,0.06},{
y
PAi
,0.97},{
y
PBi
,0.1}}];​​sys=Flatten@Join[Table[{eq1[i],eq2[i],eq3[i]},{i,2,Np}],{eq1i,eq2i,eq1[1],eq2[1],eq3[Np+1]}];
In[]:=
functions={xA[1]xA[1][s],yA[1]yA[1][s],xA[2]xA[2][s],yA[2]yA[2][s],q[2]q[2][s],xA[3]xA[3][s],yA[3]yA[3][s],q[3]q[3][s],xA[4]xA[4][s],yA[4]yA[4][s],q[4]q[4][s],xA[5]xA[5][s],yA[5]yA[5][s],q[5]q[5][s],xA[6]xA[6][s],yA[6]yA[6][s],q[6]q[6][s],xA[7]xA[7][s],yA[7]yA[7][s],q[7]q[7][s],xA[8]xA[8][s],yA[8]yA[8][s],q[8]q[8][s],xA[9]xA[9][s],yA[9]yA[9][s],q[9]q[9][s],xA[10]xA[10][s],yA[10]yA[10][s],q[10]q[10][s],xA[11]xA[11][s],yA[11]yA[11][s],q[11]q[11][s],xA[12]xA[12][s],yA[12]yA[12][s],q[12]q[12][s],xA[13]xA[13][s],yA[13]yA[13][s],q[13]q[13][s],xA[14]xA[14][s],yA[14]yA[14][s],q[14]q[14][s],xA[15]xA[15][s],yA[15]yA[15][s],q[15]q[15][s],xA[16]xA[16][s],yA[16]yA[16][s],q[16]q[16][s],xA[17]xA[17][s],yA[17]yA[17][s],q[17]q[17][s],xA[18]xA[18][s],yA[18]yA[18][s],q[18]q[18][s],xA[19]xA[19][s],yA[19]yA[19][s],q[19]q[19][s],xA[20]xA[20][s],yA[20]yA[20][s],q[20]q[20][s],thetatheta[s],stst[s],
y
PAi

y
PAi
[s],
y
PBi

y
PBi
[s]};

alpha=2

In[]:=
alpha=2;​​sol=FindRoot[sys/.st0.01,unknown,MaxIterations100];​​xA[Np+1]=.;yA[Np+1]=.;q[Np+1]=.;q[1]=.;​​xA[Np+1][s_]:=xAin;yA[Np+1][s_]:=
y
PAi
[s];q[Np+1][s_]:=
q
f
;q[1][s_]:=
q
f
(1-theta[s]);​​init=Join[Flatten[Table[{xA[i][0](xA[i]/.sol),yA[i][0](yA[i]/.sol)},{i,1,Np}],1],Flatten[Table[q[i][0](q[i]/.sol),{i,2,Np}],1],{
y
PAi
[0](
y
PAi
/.sol),
y
PBi
[0](
y
PBi
/.sol),theta[0](theta/.sol),st[0]0.01}];​​initdev=Join[Flatten[Table[{xA[i]'[0]-0.1,yA[i]'[0]-0.1},{i,1,Np}],1],Flatten[Table[{q[i]'[0]-0.1},{i,2,Np}],1],{theta'[0]-0.1,st'[0]-1,
y
PAi
'[0]-0.1,
y
PBi
'[0]-0.1}];​​aux=Sum[xA[i]'[s]^2+yA[i]'[s]^2,{i,1,Np}]+Sum[q[i]'[s]^2,{i,2,Np}]+theta'[s]^2+st'[s]^2+
y
PAi
'[s]^2+
y
PBi
'[s]^21;​​unknown1=Join[Flatten[Table[{xA[i],yA[i]},{i,1,Np}],1],Flatten[Table[q[i],{i,2,Np}],1],{theta,st,
y
PAi
,
y
PBi
}];​​SYS=Join[sys/.functions,{aux},init,initdev];​​SOL=Quiet@NDSolve[SYS,unknown1,{s,-10,10},AccuracyGoal5];​​bounds=SOL[[1]][[2,2]][[1,1]];​​p1=ParametricPlot[Evaluate[{xA[1][s],xA[1][s](1-theta[s])/xAin}/.SOL[[1]]],{s,bounds[[1]],bounds[[2]]-2},PlotRange{{0.9,1},{0,1}},FrameTrue,GridLinesAutomatic,AspectRatio1,PlotStyle{Thick,Dashed,Darker@Red},FrameTrue,FrameLabel{Style["Nitrogen Purity",16],Style["Nitrogen Recovery",16]},PlotLegendsLineLegend[{Darker@Magenta,Darker@Blue,Darker@Yellow,Darker@Green,Darker@Red},{"α=20","α=10","α=6","α=4","α=2"}]];

alpha=4

In[]:=
alpha=4;​​xA[Np+1]=xAin;yA[Np+1]=
y
PAi
;q[Np+1]=
q
f
;q[1]=
q
f
(1-theta);​​sol=FindRoot[sys/.st0.01,unknown,MaxIterations100];​​xA[Np+1]=.;yA[Np+1]=.;q[Np+1]=.;q[1]=.;​​xA[Np+1][s_]:=xAin;yA[Np+1][s_]:=
y
PAi
[s];q[Np+1][s_]:=
q
f
;q[1][s_]:=
q
f
(1-theta[s]);​​init=Join[Flatten[Table[{xA[i][0](xA[i]/.sol),yA[i][0](yA[i]/.sol)},{i,1,Np}],1],Flatten[Table[q[i][0](q[i]/.sol),{i,2,Np}],1],{
y
PAi
[0](
y
PAi
/.sol),
y
PBi
[0](
y
PBi
/.sol),theta[0](theta/.sol),st[0]0.01}];​​initdev=Join[Flatten[Table[{xA[i]'[0]-0.1,yA[i]'[0]-0.1},{i,1,Np}],1],Flatten[Table[{q[i]'[0]-0.1},{i,2,Np}],1],{theta'[0]-0.1,st'[0]-1,
y
PAi
'[0]-0.1,
y
PBi
'[0]-0.1}];​​aux=Sum[xA[i]'[s]^2+yA[i]'[s]^2,{i,1,Np}]+Sum[q[i]'[s]^2,{i,2,Np}]+theta'[s]^2+st'[s]^2+
y
PAi
'[s]^2+
y
PBi
'[s]^21;​​unknown1=Join[Flatten[Table[{xA[i],yA[i]},{i,1,Np}],1],Flatten[Table[q[i],{i,2,Np}],1],{theta,st,
y
PAi
,
y
PBi
}];​​SYS=Join[sys/.functions,{aux},init,initdev];​​SOL=Quiet@NDSolve[SYS,unknown1,{s,-10,10},AccuracyGoal5];​​bounds=SOL[[1]][[2,2]][[1,1]];​​p2=ParametricPlot[Evaluate[{xA[1][s],xA[1][s](1-theta[s])/xAin}/.SOL[[1]]],{s,bounds[[1]]+0.01,bounds[[2]]-0.01},PlotRange{{0.9,1},{0,1}},FrameTrue,GridLinesAutomatic,AspectRatio1,PlotStyle{Thick,Dashed,Darker@Green},FrameTrue,FrameLabel{Style["Nitrogen Purity",16],Style["Nitrogen Recovery",16]}];

alpha=6

In[]:=
alpha=6;​​xA[Np+1]=xAin;yA[Np+1]=
y
PAi
;q[Np+1]=
q
f
;q[1]=
q
f
(1-theta);​​sol=FindRoot[sys/.st0.01,unknown,MaxIterations100];​​xA[Np+1]=.;yA[Np+1]=.;q[Np+1]=.;q[1]=.;​​xA[Np+1][s_]:=xAin;yA[Np+1][s_]:=
y
PAi
[s];q[Np+1][s_]:=
q
f
;q[1][s_]:=
q
f
(1-theta[s]);​​init=Join[Flatten[Table[{xA[i][0](xA[i]/.sol),yA[i][0](yA[i]/.sol)},{i,1,Np}],1],Flatten[Table[q[i][0](q[i]/.sol),{i,2,Np}],1],{
y
PAi
[0](
y
PAi
/.sol),
y
PBi
[0](
y
PBi
/.sol),theta[0](theta/.sol),st[0]0.01}];​​initdev=Join[Flatten[Table[{xA[i]'[0]-0.1,yA[i]'[0]-0.1},{i,1,Np}],1],Flatten[Table[{q[i]'[0]-0.1},{i,2,Np}],1],{theta'[0]-0.1,st'[0]-1,
y
PAi
'[0]-0.1,
y
PBi
'[0]-0.1}];​​aux=Sum[xA[i]'[s]^2+yA[i]'[s]^2,{i,1,Np}]+Sum[q[i]'[s]^2,{i,2,Np}]+theta'[s]^2+st'[s]^2+
y
PAi
'[s]^2+
y
PBi
'[s]^21;​​unknown1=Join[Flatten[Table[{xA[i],yA[i]},{i,1,Np}],1],Flatten[Table[q[i],{i,2,Np}],1],{theta,st,
y
PAi
,
y
PBi
}];​​SYS=Join[sys/.functions,{aux},init,initdev];​​SOL=Quiet@NDSolve[SYS,unknown1,{s,-10,10},AccuracyGoal5];​​bounds=SOL[[1]][[2,2]][[1,1]];​​p3=ParametricPlot[Evaluate[{xA[1][s],xA[1][s](1-theta[s])/xAin}/.SOL[[1]]],{s,bounds[[1]]+0.01,bounds[[2]]-0.01},PlotRange{{0.9,1},{0,1}},FrameTrue,GridLinesAutomatic,AspectRatio1,PlotStyle{Thick,Dashed,Darker@Yellow},FrameTrue,FrameLabel{Style["Nitrogen Purity",16],Style["Nitrogen Recovery",16]}];

alpha=10

alpha=20

Nitrogen recovery versus Nitrogen purity for various separation factors

Comparison between our results and Figure 8.25 of Ref. 1 indicates close agreement, which provides a validation of the methodology described in Ref. 2.