Application Center - Maplesoft

# Algorithm to find all feasible distillation boundary maps for Ternary Systems

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

Algorithm to Find All Feasible Distillation Boundary Maps (Residue Curve Maps) for Ternary Systems

Reference: Eric J. Peterson and Lee R. Partin, Ind. Eng. Chem. Res. , 36 (1997), 1799-1811

The capabilities of Maple were very beneficial for this specialized research project. It made it easy to quickly program an algorithm to find all of the feasible residue curve maps for three-component systems. A residue curve tracks the liquid composition within a boiling flash. For a description of the research effort, please see my web site:

Step 1: Find all combinations of position numbers [1,2,3,4,5,6,7]

> S := combinat[powerset]([1,2,3,4,5,6,7]);

Step 2: Accept the combinations that include 1, 3 & 5.

> f:=i->if (has(1,S[i]) and has(3,S[i]) and has(5,S[i])) then S[i] end if:
Ss1:=[seq(f(i),i=1..nops(S))];

There are 16 valid combinations.

Steps 3 - 5i: Loop through all permuations of the combinations and apply validity rules to identify valid temperature sequences

> Ss2:=table():
j:=0:
# check all permutations for each of the 16 valid combinations
for i from 1 to 16 do
# Step 3, calculate the permutations
permutations:=combinat[permute](Ss1[i]); # finds the permutations
# check each of the permutations for validity
for k from 1 to nops(permutations) do
valid := FALSE;
order1 := 0;
order2 := 0;
order3 := 0;
order4 := 0;
order5 := 0;
order6 := 0;
order7 := 0;
item:=op(k,permutations);
# Step 4, find the order for each position number
for l from 1 to nops(item) do
if op(l,item) = 1 then order1 := l fi;
if op(l,item) = 2 then order2 := l fi;
if op(l,item) = 3 then order3 := l fi;
if op(l,item) = 4 then order4 := l fi;
if op(l,item) = 5 then order5 := l fi;
if op(l,item) = 6 then order6 := l fi;
if op(l,item) = 7 then order7 := l fi;
od;
# Step 4a, 1 - 3 - 5 must be in rank order by definition
if order1 < order3 and order3 < order5 then
valid := TRUE
fi;
# Step 4b, 2 can not be between 1 and 3 if it is an azeotrope
if order2 > order1 and order2 < order3 then
valid := FALSE
fi;
# 4 can not be between 3 and 5 if it is an azeotrope
if order4 > order3 and order4 < order5 then
valid := FALSE
fi;
# 6 can not be between 1 and 5 if it is an azeotrope
if order6 > order1 and order6 < order5 then
valid := FALSE
fi;
# Step 5, start applying the DBM sketching rules
# Step 5a, determine if position #1 is a node or saddle
out1:=1;
out2:=1;
if order2>0 and order2<order1 then out1:=0 fi;
if order6>0 and order6<order1 then out2:=0 fi;
Lpoint:=`N`;
if out1<>out2 then Lpoint:=`S` fi;
# determine if 3 is a node or saddle
out1:=0;
out2:=1;
if order2>0 and order3<order2 then out1:=1 fi;
if order4>0 and order4<order3 then out2:=0 fi;
Ipoint:=`N`;
if out1<>out2 then Ipoint:=`S` fi;
# determine if 5 is a node or saddle
out1:=0;
out2:=0;
if order4>0 and order5<order4 then out1:=1 fi;
if order6>0 and order5<order6 then out2:=1 fi;
Hpoint:=`N`;
if out1<>out2 then Hpoint:=`S` fi;
# count the pure component nodes and saddles
N1:=0;
S1:=0;
if Lpoint=`N` then N1:=1
else S1:=1 fi;
if Ipoint=`N` then N1:=N1+1
else S1:=S1+1 fi;
if Hpoint=`N` then N1:=N1+1
else S1:=S1+1 fi;
# Step 5b, count the number of binary azeotropes
B:=0;
if order2>0 then B:=1 fi;
if order4>0 then B:=B+1 fi;
if order6>0 then B:=B+1 fi;
# Step 5c, determine the status of the ternary azeotrope
count:=nops(item);
ternary:=0;
if order7>0 then ternary:=1 fi;
if order7=0 then N3:=0
# ternary node if (N1+B)<4
elif (N1+B)<4 then N3:=1
# ternary node if it is the 1st or 2nd lowest T excluding S1's
elif order7=1 then N3:=1
elif order7=2 then N3:=1
elif order7=3 and Lpoint=`S` and order1<3 then N3:=1
elif order7 = 3 and Ipoint=`S` and order3<3 then N3:=1
# ternary node if it is the 1st or 2nd highest T excluding S1's
elif order7=count then N3:=1
elif order7=(count-1) then N3:=1
elif order7 = (count-2) and Hpoint=`S` and
order5 > (count-2) then N3:=1
elif order7 = (count-2) and Ipoint=`S` and
order3 > (count-2) then N3:=1
else N3:=0 fi;
S3:=ternary-N3;
# Step 5d, use equations for N2 and S2
N2:=(2+B-N1-2*N3+2*S3)/2;
S2:=B-N2;
# infeasible if N2 or S2 values are wrong
if not (N2=0 or N2=1 or N2=2 or N2=3) then valid:=FALSE fi;
if not (S2=0 or S2=1 or S2=2 or S2=3) then valid:=FALSE fi;
# Step 5e, check the status of azeotropes
if order2=0 then LIpoint:=`0`
elif order2=1 then LIpoint:=`1`
elif order2=count then LIpoint:=`3`
elif order2<order1 then LIpoint:=`(12)`
else LIpoint:=`(34)` fi;
if order4=0 then IHpoint:=`0`
elif order4=1 then IHpoint:=`1`
elif order4=count then IHpoint:=`3`
elif order4<order3 then IHpoint:=`(12)`
else IHpoint:=`(34)` fi;
if order6=0 then LHpoint:=`0`
elif order6=1 then LHpoint:=`1`
elif order6=count then LHpoint:=`3`
elif order6<order1 then LHpoint:=`(12)`
else LHpoint:=`(34)` fi;
# Step 5f, check for binary azeotropes to fill S2 requirement
Bnodes:=0;
if LIpoint=`1` or LIpoint=`3` then Bnodes:=1 fi;
if IHpoint=`1` or IHpoint=`3` then Bnodes:=Bnodes+1 fi;
if LHpoint=`1` or LHpoint=`3` then Bnodes:=Bnodes+1 fi;
if S2 > (B-Bnodes) then valid:=FALSE fi;
# Step 5g, for ternary saddles, check temperature crosses
cross:=0;
pos1:=0;
pos2:=0;
pos3:=0;
pos4:=0;
pos5:=0;
pos6:=0;
if count = 7 and S3 = 1 then
if order1>order7 then pos1:=1 fi;
if order2>order7 then pos2:=1 fi;
if order3>order7 then pos3:=1 fi;
if order4>order7 then pos4:=1 fi;
if order5>order7 then pos5:=1 fi;
if order6>order7 then pos6:=1 fi;
if pos1<>pos2 then cross:=cross+1 fi;
if pos2<>pos3 then cross:=cross+1 fi;
if pos3<>pos4 then cross:=cross+1 fi;
if pos4<>pos5 then cross:=cross+1 fi;
if pos5<>pos6 then cross:=cross+1 fi;
if cross < 3 then valid:=FALSE fi;
fi;
# Step 5h, further determine the status of binary azeotropes
if (B-Bnodes)>0 and N2=B then
Bnodes:=B;
if LIpoint=`(12)` then LIpoint:=`1`
elif LIpoint=`(34)` then LIpoint:=`3` fi;
if IHpoint=`(12)` then IHpoint:=`1`
elif IHpoint=`(34)` then IHpoint:=`3` fi;
if LHpoint=`(12)` then LHpoint:=`1`
elif LHpoint=`(34)` then LHpoint:=`3` fi;
fi;
# case with remaining saddles,
if (B-Bnodes)>0 and S2=(B-Bnodes) then
if LIpoint=`(12)` then LIpoint:=`2`
elif LIpoint=`(34)` then LIpoint:=`4` fi;
if IHpoint=`(12)` then IHpoint:=`2`
elif IHpoint=`(34)` then IHpoint:=`4` fi;
if LHpoint=`(12)` then LHpoint:=`2`
elif LHpoint=`(34)` then LHpoint:=`4` fi;
fi;
if valid = TRUE then
# Step 5i, determine ternary azeotrope status
ternary:=``;
if N3=1 and order7<4 then ternary:=`-m` fi;
if N3=1 and order7>(count-3) then ternary:=`-M` fi;
if S3=1 then ternary:=`-S` fi;
mapping:=cat(LIpoint,IHpoint,LHpoint,ternary);
# Steps 6, 7 and 8
# Some of the mappings are not definitive and
# further testing is required. You test all of the
# possible maps to see if the temperature ordering
# matches the separatrices.
if mapping=`(34)(12)0` then
mapping:=``;
if order1<order4 then
mapping:=cat(mapping,`320 `) fi;
if order2<order5 then
mapping:=cat(mapping,`410 `) fi;
fi;
if mapping=`(34)(12)3` then
mapping:=``;
if order1<order4 then
mapping:=cat(mapping,`323 `) fi;
if order2<order6 then
mapping:=cat(mapping,`413 `) fi;
fi;
if mapping=`3(12)(34)` then
mapping:=``;
if order6<order2 then
mapping:=cat(mapping,`314 `) fi;
if order1<order4 then
mapping:=cat(mapping,`323 `) fi;
fi;
if mapping=`(34)1(12)` then
mapping:=``;
if order4<order6 then
mapping:=cat(mapping,`312 `) fi;
if order2<order5 then
mapping:=cat(mapping,`411 `) fi;
fi;
if mapping=`(34)(12)1` then
mapping:=``;
if order6<order4 then
mapping:=cat(mapping,`321 `) fi;
if order2<order5 then
mapping:=cat(mapping,`411 `) fi;
fi;
if mapping=`(12)(12)1-S` then
# There are six potential maps for the case. You test
# for all options.
# case 1: 121-S
mapping:=``;
if order7<order1 and order2<order7 and
order7<order4 and order6<order7 then
mapping:=cat(mapping,`121-S `) fi;
# case 2: 121-SH
if order7<order1 and order2<order7 and
order2<order4 and order7<order5 and
order6<order7 then
mapping:=cat(mapping,`121-SH `) fi;
# case 3: 121-SI
if order7<order1 and order2<order7 and
order7<order3 and order6<order4 and
order6<order7 then
mapping:=cat(mapping,`121-SI `) fi;
# case 4: 211-S
if order7<order2 and order4<order7 and
order7<order5 and order6<order7 then
mapping:=cat(mapping,`211-S `) fi;
# case 5: 211-SI
if order6<order2 and order7<order3 and
order4<order7 and order7<order5 and
order6<order7 then
mapping:=cat(mapping,`211-SI `) fi;
# case 6: 211-SL
if order4<order2 and order4<order7 and
order7<order5 and order6<order7 and
order7<order1 then
mapping:=cat(mapping,`211-SL `) fi;
fi;
if mapping=`(12)1(12)-S` then
# case 1: 112-S. You test for temperatures to match
# the separatices of the diagram.
mapping:=``;
if order2<order7 and order7<order3 and
order4<order7 and order7<order6 then
mapping:=cat(mapping,`112-S `) fi;
# case 2: 112-SH
if order2<order6 and order2<order7 and
order7<order3 and order4<order7 and
order7<order5 then
mapping:=cat(mapping,`112-SH `) fi;
# case 3: 112-SL
if order7<order1 and order2<order7 and
order7<order3 and order4<order7 and
order4<order6 then
mapping:=cat(mapping,`112-SL `) fi;
# case 4: 211-S
if order7<order2 and order4<order7 and
order7<order5 and order6<order7 then
mapping:=cat(mapping,`211-S `) fi;
# case 5: 211-SI
if order6<order2 and order7<order3 and
order4<order7 and order7<order5 and
order6<order7 then
mapping:=cat(mapping,`211-SI `) fi;
# case 6: 211-SL
if order4<order2 and order4<order7 and
order7<order5 and order6<order7 and
order7<order1 then
mapping:=cat(mapping,`211-SL `) fi;
fi;
if mapping=`1(12)(12)-S` then
# case 1: 112-S. You test for temperatures to match
# the separatices of the diagram.
mapping:=``;
if order2<order7 and order7<order3 and
order4<order7 and order7<order6 then
mapping:=cat(mapping,`112-S `) fi;
# case 2: 112-SH
if order2<order6 and order2<order7 and
order7<order3 and order4<order7 and
order7<order5 then
mapping:=cat(mapping,`112-SH `) fi;
# case 3: 112-SL
if order7<order1 and order2<order7 and
order7<order3 and order4<order7 and
order4<order6 then
mapping:=cat(mapping,`112-SL `) fi;
# case 4: 121-S
if order7<order1 and order2<order7 and
order7<order4 and order6<order7 then
mapping:=cat(mapping,`121-S `) fi;
# case 5: 121-SH
if order7<order1 and order2<order7 and
order2<order4 and order7<order5 and
order6<order7 then
mapping:=cat(mapping,`121-SH `) fi;
# case 6: 121-SI
if order7<order1 and order2<order7 and
order7<order3 and order6<order4 and
order6<order7 then
mapping:=cat(mapping,`121-SI `) fi;
fi;
if mapping=`(34)(12)(12)-m` and S2=2 then
# case 1: 322-m. You test for temperatures to match
# the separatices of the diagram.
mapping:=``;
if order7<order4 and order7<order6 then
mapping:=cat(mapping,`322-m `) fi;
# case2: 412-m
if order2<order5 and order7<order6 then
mapping:=cat(mapping,`412-m `) fi;
# case3: 421-m
if order2<order5 and order7<order4 then
mapping:=cat(mapping,`421-m `) fi;
fi;
if mapping=`(34)(12)(34)-M` and S2=2 then
# case 1: 324-M. You test for temperatures to match
# the separatices of the diagram.
mapping:=``;
if order1<order4 and order6<order7 then
mapping:=cat(mapping,`324-M `) fi;
# case2: 414-M
if order2<order7 and order6<order7 then
mapping:=cat(mapping,`414-M `) fi;
# case3: 423-M
if order1<order4 and order2<order7 then
mapping:=cat(mapping,`423-M `) fi;
fi;
if mapping=`(34)(34)3-S` then
# case 1: 343-S. You test for temperatures to match
# the separatices of the diagram.
mapping:=``;
if order1<order7 and order7<order2 and
order4<order7 and order7<order6 then
mapping:=cat(mapping,`343-S `) fi;
# case 2: 343-SH
if order1<order7 and order7<order2 and
order4<order2 and order5<order7 and
order7<order6 then
mapping:=cat(mapping,`343-SH `) fi;
# case 3: 343-SI
if order1<order7 and order7<order2 and
order3<order7 and order4<order6 and
order7<order6 then
mapping:=cat(mapping,`343-SI `) fi;
# case 4: 433-S
if order2<order7 and order7<order4 and
order5<order7 and order7<order6 then
mapping:=cat(mapping,`433-S `) fi;
# case 5: 433-SI
if order2<order6 and order3<order7 and
order7<order4 and order5<order7 and
order7<order6 then
mapping:=cat(mapping,`433-SI `) fi;
# case 6: 433-SL
if order2<order4 and order7<order4 and
order5<order7 and order7<order6 and
order1<order7 then
mapping:=cat(mapping,`433-SL `) fi;
fi;
if mapping=`(34)3(34)-S` then
# case 1: 334-S. You test for temperatures to match
# the separatices of the diagram.
mapping:=``;
if order7<order2 and order3<order7 and
order7<order4 and order6<order7 then
mapping:=cat(mapping,`334-S `) fi;
# case 2: 334-SH
if order6<order2 and order7<order2 and
order3<order7 and order7<order4 and
order5<order7 then
mapping:=cat(mapping,`334-SH `) fi;
# case 3: 334-SL
if order1<order7 and order7<order2 and
order3<order7 and order7<order4 and
order6<order4 then
mapping:=cat(mapping,`334-SL `) fi;
# case 4: 433-S
if order2<order7 and order7<order4 and
order5<order7 and order7<order6 then
mapping:=cat(mapping,`433-S `) fi;
# case 5: 433-SI
if order2<order6 and order3<order7 and
order7<order4 and order5<order7 and
order7<order6 then
mapping:=cat(mapping,`433-SI `) fi;
# case 6: 433-SL
if order2<order4 and order7<order4 and
order5<order7 and order7<order6 and
order1<order7 then
mapping:=cat(mapping,`433-SL `) fi;
fi;
if mapping=`3(34)(34)-S` then
# case 1: 334-S. You test for temperatures to match
# the separatices of the diagram.
mapping:=``;
if order7<order2 and order3<order7 and
order7<order4 and order6<order7 then
mapping:=cat(mapping,`334-S `) fi;
# case 2: 334-SH
if order6<order2 and order7<order2 and
order3<order7 and order7<order4 and
order5<order7 then
mapping:=cat(mapping,`334-SH `) fi;
# case 3: 334-SL
if order1<order7 and order7<order2 and
order3<order7 and order7<order4 and
order6<order4 then
mapping:=cat(mapping,`334-SL `) fi;
# case 4: 343-S
if order1<order7 and order7<order2 and
order4<order7 and order7<order6 then
mapping:=cat(mapping,`343-S `) fi;
# case 5: 343-SH
if order1<order7 and order7<order2 and
order4<order2 and order5<order7 and
order7<order6 then
mapping:=cat(mapping,`343-SH `) fi;
# case 6: 343-SI
if order1<order7 and order7<order2 and
order3<order7 and order4<order6 and
order7<order6 then
mapping:=cat(mapping,`343-SI `) fi;
fi;
j:=j+1;
Ss2[j]:=[op(item),mapping];
fi;
od;
od:
# convert the solutions into temperature profile numbers
Ss3:=table():
for i from 1 to nops(op(2,op(1,Ss2))) do
val:=rhs(op(i,op(2,op(1,Ss2))));
val2:=cat(seq(convert(val[m],string),m=1..nops(val)-1));
Ss3[i]:=[op(1,sscanf(val2,`%d`)),val[nops(val)]];
od:
# sort the solution numbers
Fsort:=(x,y)->if op(1,x)<op(1,y) then true else false fi:
solution:=sort(convert(Ss3,list),Fsort):

Results

The solution is stored as a list of lists.

> "There are " || (nops(solution)) || " valid temperature sequences.";

Solution Lists

The first item of each solution list is a valid sequence number. The second item is a character string containing the valid DBM names for that temperature sequence.

> solution;

Valid Temperature Sequences

The valid temperature sequences are extracted from the solution.

> TSequences:=[seq(solution[i][1],i=1..nops(solution))];

Valid DBM Names

The valid DBM names are extracted from the solution.

> DBMs:=sort([op({seq(op(sscanf(solution[i][2],"%s %s %s")),i=1..nops(solution))})]);

> "There are " || (nops(DBMs)) || " valid DBM names.";