Binary Batch Distillation

> restart;

We will begin by creating a list of component names

> components:=[Benzene, Toluene];

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;
member(i,components,`pos`);
RETURN(pos);
end:

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);

The number of components is

> nc := nops(components);

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;

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);

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 ( and ), 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;

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

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

where P is the system pressure, 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)];

Vapor pressures can be computed using the Antoine equation:

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

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);

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;

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

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

into which we substitute the K-value model

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

The variables appearing in these equations are

> Vars := indets(DAEqns2,name);

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);

In fact, there is only one degree of freedom since 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;

> DAEqns3 := subs(Pspec,DAEqns2);

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;

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

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

The unknown variables are

> indets(AEqns3,name);

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

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

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

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;

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);

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

The amount of liquid in the still is plotted as a function of in the following command

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

and the temperature as a function of is

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

>