In a paper by Rautenbach and Dahm (see Ref . 1), a quaternary mixture is defined as follows :where the represent the feed mole fractions and are the separation factors (with component 1 being the reference component) .Rautenbach and Dahm compared the exact solution to various pseudo - binary calculations . Here, we successfully reproduce Fig . 10 of Ref . 1 by applying the methodology described in Ref . 2 in order to obtain an exact solution (i . e ., without the pseudo - binary assumption) . Nowadays, pseudo-binary assumption does not find justification since fast computing capabilities and mathematical software can readily tackle the exact problem as show in the present notebook.
()
x
α
1,2,3or4
()
α
i1
i=1,2,3or4
eq1i=/()qAst(-/)/(qBst(-/));eq2i=/()qAst(-/)/(qCst(-/));eq3i=/()qAst(-/)/(qDst(-/));eq4i=+++1;=0.3;=0.5;=0.1;=0.1;qA=1;qB=0.3;qC=0.2;qD=0.1;=1.0;=0.1;soli=FindRoot[{eq1i,eq2i,eq3i,eq4i},{,0.7},{,0.4},{,0.1},{,0.1}];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}];
q
P
y
PAi
q
P
y
PBi
x
fA
p
l
p
h
y
PAi
x
fB
p
l
p
h
y
PBi
q
P
y
PAi
q
P
y
PCi
x
fA
p
l
p
h
y
PAi
x
fC
p
l
p
h
y
PCi
q
P
y
PAi
q
P
y
PDi
x
fA
p
l
p
h
y
PAi
x
fD
p
l
p
h
y
PDi
y
PAi
y
PBi
y
PCi
y
PDi
x
fA
x
fB
x
fC
x
fD
q
f
p
l
p
h
y
PAi
y
PBi
y
PCi
y
PDi
xAp[i_]:=2/stAdd[[i,j]]xA[j];yAp[i_]:=2/stAdd[[i,j]]yA[j]xBp[i_]:=2/stAdd[[i,j]]xB[j];yBp[i_]:=2/stAdd[[i,j]]yB[j];xCp[i_]:=2/stAdd[[i,j]]xC[j];yCp[i_]:=2/stAdd[[i,j]]yC[j];xDp[i_]:=2/stAdd[[i,j]]xD[j];yDp[i_]:=2/stAdd[[i,j]]yD[j];qp[i_]:=2/stAdd[[i,j]]q[j];eq1[i_]:=-yA[i]qp[i](xA[i]-/yA[i])eq2[i_]:=-yB[i]qp[i]qB(xB[i]-/yB[i])eq3[i_]:=-yC[i]qp[i]qC(xC[i]-/yC[i])eq4[i_]:=-yD[i]qp[i]qD(xD[i]-/yD[i])eq5[i_]:=yA[i]qp[i]q[i]xAp[i]+qp[i]xA[i]eq6[i_]:=yB[i]qp[i]q[i]xBp[i]+qp[i]xB[i]eq7[i_]:=yC[i]qp[i]q[i]xCp[i]+qp[i]xC[i]yD[i_]:=1-yA[i]-yB[i]-yC[i]xD[i_]:=1-xA[i]-xB[i]-xC[i]BC={xA[Np+1],xB[Np+1],xC[Np+1],yA[Np+1](/.soli),yB[Np+1](/.soli),yC[Np+1](/.soli),q[Np+1],q[1](1-theta)};unknown=Join[Flatten[Table[{{xA[i],0.5},{yA[i],0.6},{xB[i],0.1},{yB[i],0.3},{xC[i],0.2},{yC[i],0.05},{q[i],0.3}},{i,1,Np+1}],1],{{theta,0.6}}];sys=Flatten@Join[Table[{eq1[i],eq2[i],eq3[i],eq4[i],eq5[i],eq6[i],eq7[i]},{i,2,Np}],{eq5[1],eq6[1],eq7[Np+1],eq1[1],eq2[1],eq3[1],eq4[1]},BC];sol=FindRoot[sys/.st1,unknown,MaxIterations1000,AccuracyGoal20]
Np+1
∑
j=1
Np+1
∑
j=1
Np+1
∑
j=1
Np+1
∑
j=1
Np+1
∑
j=1
Np+1
∑
j=1
Np+1
∑
j=1
Np+1
∑
j=1
Np+1
∑
j=1
p
l
p
h
p
l
p
h
p
l
p
h
p
l
p
h
x
fA
x
fB
x
fC
y
PAi
y
PBi
y
PCi
q
f
q
f
q
f
{xA[1]0.176997,yA[1]0.41559,xB[1]0.558329,yB[1]0.470644,xC[1]0.124668,yC[1]0.0720845,q[1]0.631951,xA[2]0.177755,yA[2]0.416865,xB[2]0.55805,yB[2]0.469729,xC[2]0.124501,yC[2]0.071881,q[2]0.633959,xA[3]0.180011,yA[3]0.420644,xB[3]0.55721,yB[3]0.467012,xC[3]0.124005,yC[3]0.0712793,q[3]0.639951,xA[4]0.183718,yA[4]0.426795,xB[4]0.555804,yB[4]0.462578,xC[4]0.123196,yC[4]0.070305,q[4]0.649835,xA[5]0.188793,yA[5]0.435098,xB[5]0.553829,yB[5]0.456566,xC[5]0.122097,yC[5]0.0689989,q[5]0.663453,xA[6]0.195119,yA[6]0.445266,xB[6]0.551288,yB[6]0.449167,xC[6]0.120741,yC[6]0.0674138,q[6]0.680577,xA[7]0.202544,yA[7]0.45695,xB[7]0.548202,yB[7]0.440614,xC[7]0.119168,yC[7]0.0656113,q[7]0.700905,xA[8]0.210885,yA[8]0.469759,xB[8]0.54461,yB[8]0.431178,xC[8]0.117424,yC[8]0.0636576,q[8]0.724058,xA[9]0.219929,yA[9]0.483283,xB[9]0.540577,yB[9]0.421149,xC[9]0.115558,yC[9]0.0616192,q[9]0.749573,xA[10]0.229441,yA[10]0.497115,xB[10]0.536192,yB[10]0.410828,xC[10]0.113624,yC[10]0.0595599,q[10]0.776909,xA[11]0.23917,yA[11]0.510861,xB[11]0.531566,yB[11]0.400506,xC[11]0.111672,yC[11]0.0575377,q[11]0.805447,xA[12]0.248864,yA[12]0.524169,xB[12]0.526826,yB[12]0.390456,xC[12]0.109753,yC[12]0.0556026,q[12]0.834503,xA[13]0.258271,yA[13]0.536728,xB[13]0.522111,yB[13]0.380922,xC[13]0.107914,yC[13]0.0537962,q[13]0.86334,xA[14]0.267155,yA[14]0.548278,xB[14]0.51756,yB[14]0.372111,xC[14]0.106197,yC[14]0.0521515,q[14]0.891192,xA[15]0.2753,yA[15]0.558611,xB[15]0.513309,yB[15]0.364195,xC[15]0.104639,yC[15]0.0506933,q[15]0.917285,xA[16]0.282515,yA[16]0.567566,xB[16]0.509484,yB[16]0.357311,xC[16]0.10327,yC[16]0.0494394,q[16]0.940869,xA[17]0.288636,yA[17]0.575022,xB[17]0.506198,yB[17]0.351561,xC[17]0.102118,yC[17]0.0484023,q[17]0.961244,xA[18]0.293533,yA[18]0.580895,xB[18]0.503542,yB[18]0.347022,xC[18]0.101202,yC[18]0.0475898,q[18]0.977791,xA[19]0.297102,yA[19]0.585125,xB[19]0.501592,yB[19]0.343746,xC[19]0.100538,yC[19]0.0470069,q[19]0.989997,xA[20]0.299272,yA[20]0.587677,xB[20]0.500401,yB[20]0.341768,xC[20]0.100135,yC[20]0.0466563,q[20]0.997479,xA[21]0.3,yA[21]0.58853,xB[21]0.5,yB[21]0.341107,xC[21]0.1,yC[21]0.0465392,q[21]1.,theta0.368049}
In[]:=
yD[i_]=.
In[]:=
xD[i_]=.
In[]:=
yD[i_][s]:=1-yA[i][s]-yB[i][s]-yC[i][s]xD[i_][s]:=1-xA[i][s]-xB[i][s]-xC[i][s]
st=.;init=Join[Flatten[Table[{xA[i][0](xA[i]/.sol),yA[i][0](xA[i]/.sol),xB[i][0](xB[i]/.sol),yB[i][0](yB[i]/.sol),xC[i][0](xC[i]/.sol),yC[i][0](xC[i]/.sol),q[i][0](q[i]/.sol)},{i,1,Np+1}],1],{theta[0](theta/.sol),st[0]1}];initdev=Join[Flatten[Table[{xA[i]'[0]0.1,yA[i]'[0]0.1,xB[i]'[0]0.1,yB[i]'[0]0.1,xC[i]'[0]0.1,yC[i]'[0]0.1,q[i]'[0]0.1},{i,1,Np+1}],1],{theta'[0]-0.1,st'[0]-1}];aux=Sum[xA[i]'[s]^2+yA[i]'[s]^2+xB[i]'[s]^2+yB[i]'[s]^2+xC[i]'[s]^2+yC[i]'[s]^2+q[i]'[s]^2,{i,1,Np+1}]+theta'[s]^2+st'[s]^21;unknown1=Join[Flatten[Table[{xA[i],yA[i],xB[i],yB[i],xC[i],yC[i],q[i]},{i,1,Np+1}],1],{theta,st}];functions=;
In[]:=
SYS=Join[sys/.functions,{aux},init,initdev];
In[]:=
SOL=Quiet@NDSolve[SYS,unknown1,{s,-10,10},AccuracyGoal5];
In[]:=
solA=First@NSolve[theta[s]yAP+(1-theta[s])xA[1][s]0.3,yAP]
Out[]=
yAP
0.3-1.(1.-1.theta[s])xA[1][s]
theta[s]
In[]:=
solB=First@NSolve[theta[s]yBP+(1-theta[s])xB[1][s]0.5,yBP]
Out[]=
yBP
0.5-1.(1.-1.theta[s])xB[1][s]
theta[s]
In[]:=
solC=First@NSolve[theta[s]yCP+(1-theta[s])xC[1][s]0.1,yCP]
Out[]=
yCP
0.1-1.(1.-1.theta[s])xC[1][s]
theta[s]
In[]:=
pA=ParametricPlot[Evaluate[{theta[s],yAP/.solA}/.SOL[[1]]],{s,-10,10},PlotRange{{0.005,0.995},{0,0.8}},FrameTrue,GridLinesAutomatic,AspectRatio1,PlotStyleThick,FrameTrue,FrameLabel{Style["stage cut",16],Style["permeate mole fraction",16]},PlotPoints1000];
In[]:=
pB=ParametricPlot[Evaluate[{theta[s],yBP/.solB}/.SOL[[1]]],{s,-10,10},PlotRange{{0.005,0.995},{0,1}},FrameTrue,GridLinesAutomatic,AspectRatio1,PlotStyle{Darker@Red,Thick},FrameTrue,PlotPoints1000];
In[]:=
pC=ParametricPlot[Evaluate[{theta[s],yCP/.solC}/.SOL[[1]]],{s,-10,10},PlotRange{{0.005,0.995},{0,1}},FrameTrue,GridLinesAutomatic,AspectRatio1,PlotStyle{Darker@Green,Thick},FrameTrue,PlotPoints1000];
In[]:=
pD=ParametricPlot[Evaluate[{theta[s],1-yAP-yBP-yCP/.solA/.solB/.solC}/.SOL[[1]]],{s,-10,10},PlotRange{{0.005,0.995},{0,1}},FrameTrue,GridLinesAutomatic,AspectRatio1,PlotStyle{Darker@Magenta,Thick},FrameTrue,PlotPoints1000];
In[]:=
Show[pA,pB,pC,pD]
Out[]=
We check our results against Fig . 10 of Ref . 1
Notice the very close match between our colored curves (our calculation) and the data marked with x (data representing the exact solution from Ref. 1)
Notice the very close match between our colored curves (our calculation) and the data marked with x (data representing the exact solution from Ref. 1)
References
References
[1] R . Rautenbach and W . Dahm, The Separation of Multicomponent Mixtures by Gas Permeation, Chem . Eng . Process, vol. 19, issue 4, 1985, pp. 211 - 219.
[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
[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
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Exact simulation of the separation of quaternary mixtures by gas permeation
by Housam Binous and Ahmed Bellagi
Wolfram Community, STAFF PICKS, December 27, 2024
https://community.wolfram.com/groups/-/m/t/3345458
by Housam Binous and Ahmed Bellagi
Wolfram Community, STAFF PICKS, December 27, 2024
https://community.wolfram.com/groups/-/m/t/3345458