Wolfram Cloud Document

ABSTRACT (original article): The present paper describes the arc-length continuation method; a mathematical technique that can solve a wide variety of chemical engineering problems. Several case studies involving distillation, reaction engineering, and thermodynamics are treated with the help of this method. The approach is powerful enough to elucidate situations when several turning points occur. Finally, the authors conclude with their experience teaching this technique. Indeed, this methodology was introduced to the graduate students of the nonlinear dynamics course at King Fahd University of Petroleum & Minerals (KFUPM). CITATION (original article): Binous, H., & Shaikh, A. A. (2014). Introduction of the arc-length continuation technique in the chemical engineering graduate program at KFUPM. In Computer Applications in Engineering Education (Vol. 23, Issue 3, pp. 344–351). Wiley. DOI: https://doi.org/10.1002/cae.21604

Consider an ideal equimolar binary mixture of benzene and toluene at

760mmHg

. This mixture is fed to a 10-stage column with a total condenser and partial reboiler. The feed stage is stage 7. Column stages are counted from the top. The feed flow rate is taken equal to

10kmol/hr

. We set the values of the reflux ratio

R=3

, and the reboil ratio

s=2

We plot the top (in blue) and bottom (in red) product purities (benzene mole fraction in the distillate stream and toluene mole fraction in the bottom stream) versus the thermal quality of two-phase feed mixtures (

0<q<1

) as well as superheated feed mixtures (

q<0

). It is possible to extend the code to include sub-cooled (

q>1

) feed mixtures. The code uses a rigorous approach, which takes into account both energy and mass balances. The computation is based on the arc length continuation technique. Finally, our results are compared to those given by HYSYS (http://www.aspentech.com/hysys) when

0<q<1

. Both data sets present good agreement.

In[]:=

Rr=3.0;(*refluxratio*)sr=2.0;(*reboilratio*)

P=760;(*totalpressureinmmHg*)

vapor pressure using Antoine equations:

Aa1=6.90565;Ba1=1211.033;Ca1=220.79;Aa2=6.95464;Ba2=1344.8;Ca2=219.482;Psat1=10^(Aa1-Ba1/(T+Ca1));Psat2=10^(Aa2-Ba2/(T+Ca2));

equilibrium constants assuming Raoult’s law:

K1=Psat1/P/.TTF;K2=Psat2/P/.TTF;

equimolar feed mixture:

z1=0.5;z2=0.5;Tbpoint=T/.FindRoot[z1Psat1+z2Psat2P,{T,70}];Tdpoint=T/.FindRoot[z1P/Psat1+z2P/Psat21,{T,80}];

=138.2(80.1-80.1)+30820-96.3(80.1-80.1);

=167.5(110.6-80.1)+33330-138.2(110.6-80.1);

mixture liquid and vapor enthalpies:

H[i_]:=y[i](

+96.3(T[i]-80.1))+(1-y[i])(

+138.2(T[i]-80.1));h[i_]:=x[i](138.2(T[i]-80.1))+(1-x[i])(167.5(T[i]-80.1));h0=y[1](138.2(T[0]-80.1))+(1-y[1])(167.5(T[0]-80.1));hF=If[Tbpoint≤TF≤Tdpoint,(1-phi)(x1z(138.2(TF-80.1))+(1-x1z)(167.5(TF-80.1)))+(phi)(y1z(

+96.3(TF-80.1))+(1-y1z)(

+138.2(TF-80.1))),If[TF≤Tbpoint,(z1(138.2(TF-80.1))+(z2)(167.5(TF-80.1))),z1(

+96.3(TF-80.1))+(z1)(

+138.2(TF-80.1))]];Hd=z1(

+96.3(Tdpoint-80.1))+(z1)(

+138.2(Tdpoint-80.1));hb=z1(138.2(Tbpoint-80.1))+(z2)(167.5(Tbpoint-80.1));eq1a=y1z(K1)z1/(1+(K1-1)phi1);eq2a=x1zz1/(1+(K1-1)phi1);eq3a=y2z(K2)z2/(1+(K2-1)phi1);eq4a=x2zz2/(1+(K2-1)phi1);PSat1[i_]:=10^(Aa1-Ba1/(T[i]+Ca1));PSat2[i_]:=10^(Aa2-Ba2/(T[i]+Ca2));F=10;(*feedflowrate*)zF=0.5;(*feedcomposition*)

MESH equations:

equations for the total condenser:

eq3=Di+L0V1;eq4=Dih0+L0h0V1H[1]+Qc;eq5=V1+L1V2+L0;eq6=V1y[1]+L1x[1]V2y[2]+L0y[1];eq7=y[1]Px[1]PSat1[1];eq8=V1H[1]+L1h[1]V2H[2]+L0h0;eq9=V2+L2V3+L1;eq10=V2y[2]+L2x[2]V3y[3]+L1x[1];eq11=y[2]Px[2]PSat1[2];eq12=V2H[2]+L2h[2]V3H[3]+L1h[1];eq13=V3+L3V4+L2;eq14=V3y[3]+L3x[3]V4y[4]+L2x[2];eq15=y[3]Px[3]PSat1[3];eq16=V3H[3]+L3h[3]V4H[4]+L2h[2];eq17=V4+L4V5+L3;eq18=V4y[4]+L4x[4]V5y[5]+L3x[3];eq19=y[4]Px[4]PSat1[4];eq20=V4H[4]+L4h[4]V5H[5]+L3h[3];eq21=V5+L5V6+L4;eq22=V5y[5]+L5x[5]V6y[6]+L4x[4];eq23=y[5]Px[5]PSat1[5];eq24=V5H[5]+L5h[5]V6H[6]+L4h[4];eq25=V6+L6V7+L5;eq26=V6y[6]+L6x[6]V7y[7]+L5x[5];eq27=y[6]Px[6]PSat1[6];eq28=V6H[6]+L6h[6]V7H[7]+L5h[5];

equations for the feed stage:

eq29=V7+L7V8+L6+F;eq30=V7y[7]+L7x[7]V8y[8]+L6x[6]+FzF;eq31=y[7]Px[7]PSat1[7];eq32=V7H[7]+L7h[7]V8H[8]+L6h[6]+FhF;eq33=V8+L8V9+L7;eq34=V8y[8]+L8x[8]V9y[9]+L7x[7];eq35=y[8]Px[8]PSat1[8];eq36=V8H[8]+L8h[8]V9H[9]+L7h[7];eq37=V9+L9V10+L8;eq38=V9y[9]+L9x[9]V10y[10]+L8x[8];eq39=y[9]Px[9]PSat1[9];eq40=V9H[9]+L9h[9]V10H[10]+L8h[8];eq41=V10+L10V11+L9;eq42=V10y[10]+L10x[10]V11y[11]+L9x[9];eq43=y[10]Px[10]PSat1[10];eq44=V10H[10]+L10h[10]V11H[11]+L9h[9];

equations for the partial reboiler:

eq45=V11+BL10;eq46=V11y[11]+Bx[11]L10x[10];eq47=y[11]Px[11]PSat1[11];eq48=V11H[11]+Bh[11]L10h[10]+Qr;eq49=L0/DiRr;(*refluxratio*)eq50=V11/Bsr;(*reboilratio*)

bubble temperature equations:

eq51=PSat1[1]x[1]+PSat2[1](1-x[1])P;eq52=PSat1[2]x[2]+PSat2[2](1-x[2])P;eq53=PSat1[3]x[3]+PSat2[3](1-x[3])P;eq54=PSat1[4]x[4]+PSat2[4](1-x[4])P;eq55=PSat1[5]x[5]+PSat2[5](1-x[5])P;eq56=PSat1[6]x[6]+PSat2[6](1-x[6])P;eq57=PSat1[7]x[7]+PSat2[7](1-x[7])P;eq58=PSat1[8]x[8]+PSat2[8](1-x[8])P;eq59=PSat1[9]x[9]+PSat2[9](1-x[9])P;eq60=PSat1[10]x[10]+PSat2[10](1-x[10])P;eq61=PSat1[11]x[11]+PSat2[11](1-x[11])P;eq62=PSat1[0]y[1]+PSat2[0](1-y[1])P;

Fetch an initial condition for vapor fraction equal to 1.01 or equivalently a feed thermal quality equal to 0.01

In[]:=

phi1=1.01;

Auxiliary equation

Initial conditions

Governing equations

Unknowns

Using NDSolve to find the solution of the DAE (differential algebraic system of equations)

Result obtained using arc-length continuation

Data obtained using HYSYS

Reference

[1] Binous, H. and Shaikh, A.A. (2015), Introduction of the arc-length continuation technique in the chemical engineering graduate program at KFUPM. Computer Applications in Engineering Education, 23: 344-351. https://doi.org/10.1002/cae.21604

CITE THIS NOTEBOOK

Distillation-column product purity as function of feed thermal quality
by Housam Binous and Ahmed Bellagi
Wolfram Community, STAFF PICKS, March 24, 2024
https://community.wolfram.com/groups/-/m/t/3146544