Wolfram Cloud Document

ABSTRACT (original article): Many chemical engineering programs teach gas permeation because of its important applications in the oil and gas industry. In fact, it is considered a good replacement for classical separation techniques such as cryogenic technologies, PSA, physical and chemical absorption etc. Gas permeation module modeling can play an important role in both simulation and design of permeation units. In fact, the simulation allows the prediction of the separation performance of gas permeation modules as well as the selection of optimal operating parameters, stage cut and pressure ratio . On the other hand, the design provides useful data such as the module area. Here, we present the theoretical background of this unit operation, and the modeling of a single-stage gas permeation module used to separate multicomponent gas mixtures for the following flow patterns: (1) perfect mixing, (2) one-side mixing, (3) cross flow, (4) co-current flow, (4) counter-current flow. Our numerical simulation methodology is based on Chebyshev orthogonal collocation and arc-length continuation techniques. We apply the proposed methodology to two industrially important systems, which are described in reference (Shindo et al. in Sep Sci Technol 20:445–459, 1985. https:// doi. org/ 10. 1080/ 01496 39850 80606 92): (1) separation of NH3, H2, and N2 gaseous mixture by means of a polyethylene membrane. Such a gas mixture arises in the Haber process. (2) Separation of H2, CH4, CO, N2, and CO2 mixture through a microporous glass membrane. This mixture is important in the partial oxidation of fuels with air and in the steam gasification of biomass. Whenever possible, our results are benchmarked against data obtained from the open literature or predicted using the software Aspen-Plus®. Furthermore, we performed sensitivity analyses where we vary either the stage cut or the pressure ratio in order to investigate the effect of these parameters on both the permeate and retentate compositions as well as on the surface area of the permeation module. Finally, we provide useful tips and computer coding to readers who wish to investigate similar gas permeation problems using our new methodology. CITATION (original article): H. Binous, U. Zahid and A. Bellagi, Modeling for Multi-component Gas Permeation with Orthogonal Collocation and Arc-Length Continuation, Chemistry Africa (2024). DOI: https://doi.org/10.1007/s42250-024-01167-y

This notebook constitutes part B of the Electronic Supplementary Material for Ref. 1.

The coding, based on orthogonal collocation and arc-length continuation, constitutes a model for separation of a five-component mixture (H2, CH4, CO, N2, and CO2) by gas permeation. This mixture was studied previously by Shindo and coworkers (Ref. 2). Here, the case of a counter-current gas permeation module is considered but the model can treat other flow patterns such as (1) perfect mixing, (2) one-side mixing, (3) cross flow and (4) co-current flow.

Dimensionless area and permeate mole fractions versus stage cut are computed and plotted. Finally, we compare our dimensionless area with its counterpart obtained using our Mathematica’s extension of the approach adopted by Wankat (see Ref. 3) and Coker and coworkers (see Ref. 4-6).

Feed composition:

=0.3;

=0.1;

=0.25;

=0.15;

=0.2;

Permeabilities:

qA=1;qB=0.398;qC=0.292;qD=0.288;qE=0.308;

In[]:=

=1.0;

Pressure ratio:

0.1;g=0.1;

Chebyshev orthogonal collocation first we build matrix Add which will allow for the derivatives (calculated at the collocation nodes)

In[]:=

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}];

Derivatives are indicated by qp, xAp, xBp, xCp, xDp, yAp, yBp, yCp, and yDp

In[]:=

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]xBp[i_]:=2/st

Np+1

∑

j=1

Add[[i,j]]xB[j];xCp[i_]:=2/st

Np+1

∑

j=1

Add[[i,j]]xC[j];yBp[i_]:=2/st

Np+1

∑

j=1

Add[[i,j]]yB[j];yCp[i_]:=2/st

Np+1

∑

j=1

Add[[i,j]]yC[j];yDp[i_]:=2/st

Np+1

∑

j=1

Add[[i,j]]yD[j];qp[i_]:=2/st

Np+1

∑

j=1

Add[[i,j]]q[j];xDp[i_]:=2/st

Np+1

∑

j=1

Add[[i,j]]xD[j];xEp[i_]:=2/st

Np+1

∑

j=1

Add[[i,j]]xE[j];

The governing equations y stands for permeate side mol fraction and x for retentate side mol fraction

eq1[i_]:=-xAp[i]q[i]qA(xA[i]-

yA[i])-xA[i](qA(xA[i]-

yA[i])+qB(xB[i]-

yB[i])+qC(xC[i]-

yC[i])+qD(xD[i]-

yD[i])+qE(xE[i]-

yE[i]))eq2[i_]:=-xBp[i]q[i]qB(xB[i]-

yB[i])-xB[i](qA(xA[i]-

yA[i])+qB(xB[i]-

yB[i])+qC(xC[i]-

yC[i])+qD(xD[i]-

yD[i])+qE(xE[i]-

yE[i]))eq3[i_]:=-xCp[i]q[i]qC(xC[i]-

yC[i])-xC[i](qA(xA[i]-

yA[i])+qB(xB[i]-

yB[i])+qC(xC[i]-

yC[i])+qD(xD[i]-

yD[i])+qE(xE[i]-

yE[i]))eq4[i_]:=-xDp[i]q[i]qD(xD[i]-

yD[i])-xD[i](qA(xA[i]-

yA[i])+qB(xB[i]-

yB[i])+qC(xC[i]-

yC[i])+qD(xD[i]-

yD[i])+qE(xE[i]-

yE[i]))eq5[i_]:=qp[i]-(qA(xA[i]-

yA[i])+qB(xB[i]-

yB[i])+qC(xC[i]-

yC[i])+qD(xD[i]-

yD[i])+qE(xE[i]-

yE[i]))eq6[i_]:=(q[i]-(1-theta))yA[i]-(1-theta)xA[1]+q[i]xA[i]eq7[i_]:=(q[i]-(1-theta))yB[i]-(1-theta)xB[1]+q[i]xB[i]eq8[i_]:=(q[i]-(1-theta))yC[i]-(1-theta)xC[1]+q[i]xC[i]eq9[i_]:=(q[i]-(1-theta))yD[i]-(1-theta)xD[1]+q[i]xD[i]

Summation rules:

yE[i_]:=1-yA[i]-yB[i]-yC[i]-yD[i]xE[i_]:=1-xA[i]-xB[i]-xC[i]-xD[i]

Equation in order to compute the permeate-side compositions at node j=1 (facing the outlet retentate stream) for the counter-current permeation module:

eq1i=

PAi

PBi

)qA(xA[1]-

PAi

)/(qB(xB[1]-

PBi

));eq2i=

PAi

PCi

)qA(xA[1]-

PAi

)/(qC(xC[1]-

PCi

));eq3i=

PAi

PDi

)qA(xA[1]-

PAi

)/(qD(xD[1]-

PDi

));eq4i=

PAi

PEi

)qA(xA[1]-

PAi

)/(qE(xE[1]-

PEi

));eq5i=

PAi

PBi

PCi

PDi

PEi

1;

In[]:=

xA[Np+1]=

;xB[Np+1]=

;xC[Np+1]=

;xD[Np+1]=

;

In[]:=

yA[1]=

PAi

;yB[1]=

PBi

;yC[1]=

PCi

;yD[1]=

PDi

;q[Np+1]=

;q[1]=

(1-theta);

Number of unknowns is 185

In[]:=

unknown=Join[Flatten[Join[Table[{{xA[i],0.2},{xB[i],0.1},{xC[i],0.2},{xD[i],0.1}},{i,1,Np}],Table[{{yA[i],0.8},{yB[i],0.3},{yC[i],0.2},{yD[i],0.1}},{i,2,Np+1}],Table[{{q[i],0.6

}},{i,2,Np}]],1],{{theta,0.42},{

PAi

,0.4},{

PBi

,0.1},{

PCi

,0.2},{

PDi

,0.1},{

PEi

,0.15}}];

In[]:=

Length[unknown]

Out[]=

185

Forming the system of equations in total we have 185 equations

In[]:=

sys=Flatten@Join[Table[{eq1[i],eq2[i],eq3[i],eq4[i],eq5[i],eq6[i],eq7[i],eq8[i],eq9[i]},{i,2,Np}],{eq1i,eq2i,eq3i,eq4i,eq5i,eq1[1],eq2[1],eq3[1],eq4[1],eq5[1],eq6[Np+1],eq7[Np+1],eq8[Np+1],eq9[Np+1]}];

In[]:=

Length[sys]

Out[]=

185

The solution we get using FindRoot is in agreement with Ref.1 when g=0.1 and st=1.0

Assume total dimensionless area (st) and fetch for the stage cut (theta) and permeate and retentate composition at outlet

st=1;

In[]:=

sol=FindRoot[sys/.st1,unknown,MaxIterations100]

Out[]=

{xA[1]0.17664,xB[1]0.106743,xC[1]0.300035,xD[1]0.180837,xA[2]0.17741,xB[2]0.106726,xC[2]0.299708,xD[2]0.180633,xA[3]0.179705,xB[3]0.106674,xC[3]0.298736,xD[3]0.180028,xA[4]0.183473,xB[4]0.106579,xC[4]0.297144,xD[4]0.179038,xA[5]0.188629,xB[5]0.106432,xC[5]0.294975,xD[5]0.17769,xA[6]0.195049,xB[6]0.106226,xC[6]0.292289,xD[6]0.176023,xA[7]0.202576,xB[7]0.10595,xC[7]0.289159,xD[7]0.174082,xA[8]0.211017,xB[8]0.105602,xC[8]0.28567,xD[8]0.171924,xA[9]0.220151,xB[9]0.105181,xC[9]0.281921,xD[9]0.169607,xA[10]0.229734,xB[10]0.104696,xC[10]0.278013,xD[10]0.167195,xA[11]0.239511,xB[11]0.104156,xC[11]0.27405,xD[11]0.164754,xA[12]0.249225,xB[12]0.103579,xC[12]0.270137,xD[12]0.162346,xA[13]0.258625,xB[13]0.102985,xC[13]0.26637,xD[13]0.160031,xA[14]0.267478,xB[14]0.102395,xC[14]0.26284,xD[14]0.157864,xA[15]0.275573,xB[15]0.101831,xC[15]0.259626,xD[15]0.155893,xA[16]0.282727,xB[16]0.101314,xC[16]0.256796,xD[16]0.154159,xA[17]0.288785,xB[17]0.100864,xC[17]0.254406,xD[17]0.152696,xA[18]0.293622,xB[18]0.100496,xC[18]0.252503,xD[18]0.151531,xA[19]0.297143,xB[19]0.100224,xC[19]0.25112,xD[19]0.150685,xA[20]0.299283,xB[20]0.100056,xC[20]0.250281,xD[20]0.150172,yA[2]0.372917,yB[2]0.102459,yC[2]0.216732,yD[2]0.128964,yA[3]0.374969,yB[3]0.102243,yC[3]0.21596,yD[3]0.128498,yA[4]0.378321,yB[4]0.101887,yC[4]0.214701,yD[4]0.127737,yA[5]0.382872,yB[5]0.101398,yC[5]0.212994,yD[5]0.126706,yA[6]0.388488,yB[6]0.100786,yC[6]0.210893,yD[6]0.125438,yA[7]0.395001,yB[7]0.100066,yC[7]0.208462,yD[7]0.123971,yA[8]0.402221,yB[8]0.0992532,yC[8]0.205774,yD[8]0.12235,yA[9]0.40994,yB[9]0.0983693,yC[9]0.202909,yD[9]0.120624,yA[10]0.417942,yB[10]0.0974372,yC[10]0.199948,yD[10]0.11884,yA[11]0.426011,yB[11]0.0964814,yC[11]0.19697,yD[11]0.117047,yA[12]0.433941,yB[12]0.0955275,yC[12]0.194051,yD[12]0.115291,yA[13]0.441537,yB[13]0.0946004,yC[13]0.191262,yD[13]0.113614,yA[14]0.448626,yB[14]0.093724,yC[14]0.188666,yD[14]0.112053,yA[15]0.455057,yB[15]0.0929198,yC[15]0.186315,yD[15]0.11064,yA[16]0.460701,yB[16]0.0922071,yC[16]0.184256,yD[16]0.109404,yA[17]0.465454,yB[17]0.0916021,yC[17]0.182524,yD[17]0.108364,yA[18]0.469232,yB[18]0.091118,yC[18]0.18115,yD[18]0.107539,yA[19]0.471973,yB[19]0.090765,yC[19]0.180153,yD[19]0.106941,yA[20]0.473634,yB[20]0.0905504,yC[20]0.17955,yD[20]0.106578,yA[21]0.474191,yB[21]0.0904783,yC[21]0.179348,yD[21]0.106457,q[2]0.587722,q[3]0.594605,q[4]0.605946,q[5]0.621547,q[6]0.641128,q[7]0.664321,q[8]0.690669,q[9]0.719627,q[10]0.750563,q[11]0.782765,q[12]0.815457,q[13]0.847811,q[14]0.878977,q[15]0.908103,q[16]0.93437,q[17]0.957021,q[18]0.975388,q[19]0.98892,q[20]0.997209,theta0.414585,

PAi

0.372227,

PBi

0.102531,

PCi

0.216992,

PDi

0.129121,

PEi

0.179128}

Now we can compute the outlet permeate mole fractions

Arc-length coding

Initial conditions:

Initial conditions for derivatives:

The auxiliary equation:

The dependent variables we aim to find:

The system augmented with the auxiliary equation and the initial conditions all variable should now be a function of s, which is the arc-length

A check of the number of equations and initial conditions will show 3x(185+1)=558

This is just exploitation of our results-Plots and so forth

Applying the method of Wankat (Ref . 3) and Coker (Ref . 4-6 )

References

◼

[1] 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] Y. Shindo, T. Hakuta, H. Yoshitome and H. Inoue H (1985) Calculation methods for multicomponent gas separation by permeation. Sep Sci Technol 20:445–459. https://doi.org/10.1080/01496398508060692

◼

[3] P. C. Wankat, Separation Process Engineering, 2nd edition, Pearson Education, 2006.

◼

[4] D.T. Coker, R. Prabhakar and B.D. Freeman (2003) Tools for teaching: gas separation using polymers. Chem Eng Educ 37:60–67

◼

[5] D.T. Coker, B.D. Freeman and G. K. Fleming (1998) Modeling multicomponent gas separation using hollow-fiber membrane contactors. AIChE J 44:1289–1302. https://doi.org/10.1002/aic.690440607

◼

[6] D.T. Coker, T. Allen, B. D. Freeman and G. K. Fleming GK (1999) Nonisothermal model for gas separation hollow-fiber membranes. AIChE J 45:1451–1468. https://doi.org/10.1002/aic.690450709

CITE THIS NOTEBOOK

Multicomponent gas separation using a counter-current permeation model
by Housam Binous, Umer Zahid & Ahmed Bellagi
Wolfram Community, STAFF PICKS, December 23, 2024
https://community.wolfram.com/groups/-/m/t/3343566