Application Center - Maplesoft

App Preview:

Binary batch distillation

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


Binary Batch Distillation

> restart;

We will begin by creating a list of component names

> components:=[Benzene, Toluene];

[Maple Math]

In order to shorten the equations that we shall derive we will identify the components using numerical subscripts that represent the position of the component in the list of their names:

> Label := proc(i)
global components;
local pos;

Other methods of labelling the components are possible (first letters, chemical formulae etc); the method used here happens to be simple and widely used. The list of component identities now becomes

> compid := map(Label,components);

[Maple Math]

The number of components is

> nc := nops(components);

[Maple Math]

Now for the equations that model the process. The amount of liquid in the still is given by

> DiffEqns :=Diff(L,x[2])=L/x[2]/(K[2]-1): DiffEqns;

[Maple Math]

The composition of the vapor leaving and the liquid remaining in the still is related by the equilibrium equations

> i:='i': Eqmeqn := i->y[i]=K[i]*x[i]: Eqmeqn(i);

[Maple Math]

The K 's are the so-called K-values or equilibrium ratios, the calculation of which will be discussed below. The K 's are, in general, complicated functions of the composition of both phases ( [Maple Math] and [Maple Math] ), temperature ( T ) and pressure ( P ); the latter variables not appearing explicitly in any of the above expressions. We have a binary system so we need two of these equations. The system of equations is completed by forcing the mole fractions to sum to unity

> i:='i': AEqns := [seq(Eqmeqn(i),i=compid),add(y[i],i=compid)=1,add(x[i],i=compid)=1]: AEqns;

[Maple Math]

The K-values for this system are to be computed using Raoult's law

> RAOULT := i->K[i] = P[i,sat]/P: RAOULT(i);

[Maple Math]

where P is the system pressure, [Maple Math] is the vapor pressure. We form a list of the K-values for all species to be used later.

> Kvalues := [seq(RAOULT(i),i=compid)];

[Maple Math]

Vapor pressures can be computed using the Antoine equation:

> Antoine := (T, A, B, C) -> 10^(A - B / (T + C));

[Maple Math]

Parameters for this equation are as follows:

> params:={AntC[Benzene] = 220.79, AntB[Benzene] = 1211.033, AntA[Benzene] = 6.90565, AntC[Toluene] = 219.482, AntB[Toluene] = 1344.8, AntA[Toluene] = 6.95464}:

which we modify below to use the indexing system employed here.

> Antparams:=subs(seq(components[i]=i,i=compid),params);

[Maple Math]

In order to carry out the numerical compuations we substitute the parameters into the K-value model equations.

> Kval2 := subs(seq(P[i,sat]=Antoine(T,AntA[i],AntB[i],AntC[i]),i=compid),Antparams,Kvalues): Kval2;

[Maple Math]

and the complete set of (differential-algebraic) equations formed by combining the differential and algrbraic equations.

> DAEqns := [DiffEqns,op(AEqns)]: DAEqns;

[Maple Math]

into which we substitute the K-value model

> DAEqns2:=subs(Kval2,DAEqns): DAEqns2;

[Maple Math]

The variables appearing in these equations are

> Vars := indets(DAEqns2,name);

[Maple Math]

Note that this system of DAEs includes both temperature and pressure as implicit variables. The number of degrees of freedom is

> nops(Vars)-nops(DAEqns);

[Maple Math]

In fact, there is only one degree of freedom since [Maple Math] must be specified at the start of the integration. We specify the pressure as 1.2 atmosphere (to be converted to mm Hg) and recreate the DAE system.

> Pspec := P = 1.2*760: Pspec;

[Maple Math]

> DAEqns3 := subs(Pspec,DAEqns2);

[Maple Math]
[Maple Math]

The initial value of all the algebriac variables is determined by solving the only algebraic equations.

> AEqns2:=subs(seq(Kval2[i],i=compid),AEqns): AEqns2;

[Maple Math]

The pressure and initial mole fractions in the liquid are fixed.

> AEqns3:=subs(Pspec,x[2]=0.4,AEqns2): AEqns3;

[Maple Math]

The unknown variables are

> indets(AEqns3,name);

[Maple Math]

We solve these equations using Newton's method.

> read `e:/maple/numerics/nonlin/newton.mpl`:

> start:=Newton(AEqns3,[x[1]=0.3,y[1]=0.5,y[2]=0.5,T=80]);

[Maple Math]

The complete starting point is given by adding the initial values of [Maple Math] and [Maple Math] to this list

> start2:=[x[2]=0.4,L=100,op(start)];

[Maple Math]

Maple out of the box is not capable of solving DAE problems. However, we have implemented a numerical method for solving DAE systems called BESIRK, the code for which is read into Maple now.

> read `e:/maple/numerics/integ/besirk`:

The equations are to be integrated until the liquid contains 80% toluene. The integration range is, therefore:

> xrange := 0.4..0.8;

[Maple Math]

We integrate the DAE system using BESIRK in the next command. Note that BESIRK can take a number of optional commands to limit the step size, or to increase the accuracy of the calculations. None of these options are needed here.

> result:=BESIRK(DAEqns3,start2,xrange);

[Maple Math]

In this example the output is not too extensive (it is an easy problem) and we have included it so we can see the structure. The result is an array with each column represents the sequence of values for ther variable whose name appears at the top of the column. The last values are shown in the bottom row. We make a list that identifies the variables with their final values as follows:

> [seq(result[1,i]=result[linalg[rowdim](result),i],i=1..linalg[coldim](result))];

[Maple Math]

The amount of liquid in the still is plotted as a function of [Maple Math] in the following command

> datatableplot(result,[x[2],L],color=red);

[Maple Plot]

and the temperature as a function of [Maple Math] is

> datatableplot(result,[x[2],T],color=red);

[Maple Plot]