Reaction Equilibrium for Multiple Gas Phase Reactions

We begin by creating a list of component identities and to calculate the number of components.

> components:=[A,B,C,D,X,Y,Z]; nc:=nops(components);

Next, we create a list of reactions:

> Reactions:=[R1,R2,R3];

Each reaction may be expressed as a Maple equation:

> Reaction[R1] :=A+B=C+D;

> Reaction[R2] := B+C=X+Y;

> Reaction[R3] :=A+X=Z;

The stoichiometric cofficients can be deduced from the set of reactions as follows:

> printlevel:=2:
for r in Reactions do
for i in components do
nu[i,r]:=coeff(lhs(Reaction[r])-rhs(Reaction[r]),i);
od: od;
printlevel:=1:

We need mole balances for each component:

> r:='r':
for i in components do
CMB[i] := C[i,0] = C[i] + add(nu[i,r]*e[r],r=Reactions);
print(CMB[i]);
od:

where we use e to denote the extent of reaction. The initial concentration is denoted with a second subscript of zero.

The reaction equilibrium constants may be expressed as

> i:='i': for r in Reactions do
Kdef[r]:=K[r] = mul ((C[i])^(-nu[i,r]),i=components);
print(Kdef[r]);
od:

We create a set of equations to be solved

> Eqns := {seq(CMB[i],i=components),seq(Kdef[r],r=Reactions)};

and a set of unknown variables:

> Vars := indets(Eqns,name);

The number of degrees of freedom follows directly

> DegFree := nops(Vars)-nops(Eqns);

Case 1

> Specs[1] := C[A,0]=1.5;

> Specs[2] := C[B,0]=1.5;

> Specs[3] := seq(C[i,0]=0,i=[C,D,X,Y,Z]);

At the conditions of interest the reaction equilibrium coefficient has the following values

> Kvals := {K[R1] = 1, K[R2]=2.63, K[R3]=5};

We now augment the set of equations with the set of specifications

> AllEqns:=Eqns union {seq(Specs[k],k=1..3),op(Kvals)};

and ask Maple to solve the problem

> result:=solve(AllEqns,Vars);

Maple found several solutions to this problem without much effort. The number of solutions can be counted using the nops command (note that in this case we must enclose result in braces).

> nops({result});

It is easy to see, however, that not all of the solutions are physically meaningful; some of them contain negative concentrations, while others are complex. Only one of the solutions is physically meaningful (all quantities being non-negative). We can pick this solution from the set of 8 by hand or we can be more sophisticated and ask Maple to find it for us.

> result2:=remove(hastype,{result},negative);

Note that we did not have to rearrange the equations into a form more suitable for numerical computation. We did not have to provide initial guesses of the unknown variables, and we did not have to work out the independent equations since we could do that with Maple as well. Of course, it is possible to compute a purely numerical solution using the Newton's method package in the share library and demonstrated below.

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

Newton requires two arguments, a list of equations to be solved and a list of equations that provide the initial values. The list of eqns is constructed from the set used in the last example

> Eqns:=[op(AllEqns)];

The startingpoint is constructed in a rather simple-minded way (that is effective here). We make a list of the variables in the equations:

> Vars := [op(indets(Eqns,name))];

and give all of them the initial value of 0.1.

> Start := [seq(x=0.1,x=Vars)];

> result:=Newton(Eqns,Start);

>

This is the result obtained earlier. As an exercise we suggest that readers try different starting points. It may be advisable in some cases to arrange the equations in a form that does not involve possible divisions by zero (see the equilibrium equations).