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:
http://users.chartertn.net/lpartin
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.";