McCabe-Thiele analysis for a pinch at top of column as a function of α.

This notebook has been written in  Mathematica by 

Mark J. McCready
Professor of Chemical and Biomolecular Engineering
University of Notre Dame
Notre Dame, IN 46556
USA

mjm@nd.edu
http://chemeprof.com

It is copyrighted to the extent allowed by which ever laws pertain to the World Wide Web and the Internet.

I would hope that as a professional courtesy, this notice remain visible to other users. 
There is no charge for copying and dissemination   

Version: 6/29/2021

This notebook addresses how the expected number of trays (stages) would change for a tough separation as the relative volatility is varied in a range of values slightly above 1.

Using  α as the relative volatility,  the number of trays varies as 1/(α-1), or 1/Log[α], which are about the same for α close to 1.   This is assuming constant α and is for the limit infinite reflux ratio.  This result, with the same α scaling, is similar to the Fenske equation (King, Separation Processes, 1971, p 456) which is derived for constant α and total reflux but not near a pinch.

Equilibrium relation

The definition of relative volatility is

eq1 = y1/x1/( y2/x2) == α

eq2 = eq1/.{y2-> 1- y1, x2-> 1-x1}

For this analysis we want to linearize this expression around x1=1.

why1 = y1/.Solve[eq2,y1][[1]]

whyexpand = Series[why1,{x1,1,1}]

Let’s check the expansion:

Plot[{(((x1 α)/(1-x1+x1 α)))/.α->2, (1+((x1-1)/α))/.α->2},{x1,0,1}]

The range of linear behavior is (obviously) wider if α is closer to 1.

Plot[{(((x1 α)/(1-x1+x1 α)))/.α->1.1, (1+((x1-1)/α))/.α->1.1},{x1,0,1}]

eqgraph = Plot[{(((x1 α)/(1-x1+x1 α)))/.α->1.3, (1+((x1-1)/α))/.α->1.3},{x1,.6,1}]

So the equilibrium line at the top is y = 1/α + 1-1/α

operating line

op1 = y -y0 == L/G ( x- x0)

whyop = y/. Solve[op1,y][[1]]

Since y0=x0

whyop1 = whyop/.x0->y0

opgraph = Plot[whyop1/.{y0->.99,L->.95,G->1},{x,.6,.99},PlotStyle->Red]

Show[eqgraph,opgraph]

So there is some accessible range for the linear approximation.  It looks like the L/G will need to be pretty close to unity if α is this small.

Infinite reflux ratio

For an infinite reflux ratio, each step “starts”on the 45 degree line, (y=x).  This makes calculating successive steps just a matter of using a simple recursion relation.  Each step is just a “move” at constant y from the y=x line to the left until the x value of the intersection with the equilibrium line is found.

To do this we solve the equilibrium relation for x

exeq = x/.Solve[ y == 1+((x-1)/α),x][[1]]

To do the stepping, we would start at y = y0 (say =.99) and calculate x.  For α = 1.3 we get

exeq/.{y->.99, α-> 1.3}

Mathematica can apply a function recursively using the “NestList” command.  Let’s do 5 steps.

exp1 =NestList[1-α+α # &,y0,5]

exp2=exp1/.{y0-> .99, α->1.3}

To make a plot lets try 10 steps.  We want to mimic the McCabe - Thiele graphical construction.

numbs =( NestList[1-α+α # &,y0,10])/.{y0-> .99, α->1.3}

Make a list of pairs of these numbers for the 45 line intersections:

ans1 = Transpose[{numbs,numbs}]

Now get the eq line intersections which are pairs of a the previous “45” value as y and the new “x” from the 1-α + α y function.  You might have to trust me!

ans2 = Transpose[{ Take[numbs,{2,Length[numbs]}],Take[numbs,Length[numbs]-1]
			}]

Now get the eq and “45” pairs alternating in 1 list for plotting:

ans3=Flatten[Table[{ans1[[i]],ans2[[i]]},{i,1,Length[ans2]}],1]

One more point at the bottom

ans4=ans3[[Length[ans3]]][[2]]

Put this one into the list

steps=Join[ans3,{{ans4,ans4}}]

Plot this list which is the steps:

stepplot = ListPlot[steps,Joined->True]

Plot the equilibrium line

eqline = Plot[(1+((x-1)/α))/.α->1.3,{x,.86,1},PlotStyle->Green]

plot y = x

whyex = Plot[x,{x,.86,1},PlotStyle->Red]

fullresult =Show[stepplot,whyex,eqline]

We see that the 10 steps get a mixture at ~.86 up to .99.

How will this change with α?

To make a plot lets try 10 steps.  We want to mimic the McCabe - Thiele graphical construction.

numbs =( NestList[1-α+α # &,y0,10])/.{y0-> .99, α->1.2}

Make a list of pairs of these numbers for the 45 line intersections:

ans1 = Transpose[{numbs,numbs}]

Now get the eq line intersections which are pairs of a the previous “45” value as y and the new “x” from the 1-α + α y function.  You might have to trust me!

ans2 = Transpose[{ Take[numbs,{2,Length[numbs]}],Take[numbs,Length[numbs]-1]
			}]

Now get the eq and “45” pairs alternating in 1 list for plotting:

ans3=Flatten[Table[{ans1[[i]],ans2[[i]]},{i,1,Length[ans2]}],1]

One more point at the bottom

ans4=ans3[[Length[ans3]]][[2]]

Put this one into the list

steps=Join[ans3,{{ans4,ans4}}]

Plot this list which is the steps:

stepplot = ListPlot[steps,Joined->True]

Plot the equilibrium line

eqline = Plot[(1+((x-1)/α))/.α->1.2,{x,.86,1},PlotStyle->Green]

plot y = x

whyex = Plot[x,{x,.86,1},PlotStyle->Red]

fullresult =Show[stepplot,whyex,eqline]

We see that the 10 steps get a mixture at ~.93up to .99 if α is reduced.

Can we pick a range and α and solve for the needed stages?

ExpandAll[NestList[1-α+α # &,y0,5]]

We see that there is a starting composition that would work for any chosen number of stages.  For, say 5 stages, if we started with x = .975, we could get to .99.

1-(α^5)+y0 (α^5)/.{y0->.99,α->1.2}

Since the result of the recurring function is simple, we can pick a desired starting mixture and calculate the number stages to get to, say, y=0.99.

We just the value of “n” after we pick α, x0 and a starting x.    Let’s start at .86 to match the first example.

FindRoot[((1-α^n + x0 α^n) /.{x0->.99, α->1.3})== .86, {n,20}]

The result is 5 stages as we expected!

Variation of n with α.

alphaplot=Plot[n/.FindRoot[((1-α^n + x0 α^n) /.{x0->.99})== .86, {n,20}][[1]],{α,1.01,1.3},PlotRange->{0,300}, AxesLabel->{"α","number of trays"}]

We can guess the functional form:

guessplot = Plot[2.7/(α-1),{α,1.01,1.3},PlotStyle->Red,PlotRange->All]

Show[alphaplot,guessplot]

Sort of amazing!    Perhaps there is some explicit formula?

Solve[(1-α^n + x0 α^n)== β,n]

Mathematica knows it!   We see that the result is a version of the “Fenske” equation.  (King, Separation Processes, 1971, p 456)

guessplot2 = Plot[(((Log[((-1+β)/(-1+x0))])/(Log[α])))/.{β->.86,x0->.99},{α,1.01,1.3},PlotStyle->Orange,PlotRange->All]

Sure enough it works!

Show[alphaplot,guessplot,guessplot2]

How could I guess so well?

Series[((Log[((-1+β)/(-1+x0))])/(Log[α])),{α,1,2}]

What is the numerical value of the numerator!

Log[((-1+β)/(-1+x0))]/.{β->.86,x0->.99}

Close to what I guessed from the intersection at α=1.01.

Have to say this is good fun!

Now what if want to allow for a finite reflux ratio?

For this case there is more work to do and the result is much less elegant.

The calculations for when the operating line is not the 45 degree seem to be best done in a loop.  We pick “g” for the intersections with the equilibrium line (although we could use the complete relation if we wished.) and “f” for the intersections with the operating line.
Again we start at x0=y0=.99

g[0] =.99

f[0] = .99

Here is the loop for α=1.3.  We do 10 steps.  First for infinite reflux ratio.

Making operating line 45.

Do[g[i] = 1-1.3+1.3 g[i-1];
f[i] =    1/1  g[i]+(1-1) .99,{i,1,10}]

The equilibrium intersections are:

ans11=Table[{g[i],g[i-1]},{i,1,10}]

The operating line intersections are:

ans22=Table[{g[i],f[i]},{i,1,10}]

As before

ans3=Flatten[Table[{ans11[[i]],ans22[[i]]},{i,1,Length[ans22]}],1]

Need the first point

Join[{{.99,.99}},ans3]

test1 = ListPlot[%,Joined->True]

eqtest = Plot[(1+((x-1)/α))/.α->1.3,{x,.86,1},PlotStyle->Green]

optest = Plot[(((L x+G y0-L y0)/G))/.{L->1,G->1,y0->.99},{x,.86,.99},PlotStyle->Red]

testans =Show[test1,eqtest,optest]

from above

%445

You can trust me that it matches above”

Show[testans,%445]

Making operating line an arbitrary value

L=.9, G=1

Clear[f,g]

g[0] = .99

f[0] =.99

Do[g[i] = 1-1.3+1.3 f[i-1];
f[i] =    .9/1  g[i]+(1-.9) .99,{i,1,10}]

ans111=Table[{g[i],f[i-1]},{i,1,10}]

ans222=Table[{g[i],f[i]},{i,1,10}]

ans31=Flatten[Table[{ans111[[i]],ans222[[i]]},{i,1,Length[ans222]}],1]

ans41 = Join[{{.99,.99}},ans31]

steps2 = ListPlot[ans41,Joined->True]

eq2plot = Plot[(1+((x-1)/α))/.α->1.3,{x,.86,1},PlotStyle->Green]

optest = Plot[(((L x+G y0-L y0)/G))/.{L->.9 ,G->1,y0->.99},{x,.86,.99},PlotStyle->Red]

Show[steps2,eq2plot,optest]

So we could only go from 0.92 to 0.99 at this L/G.    

We can make it a little worse.

Making operating line an arbitrary value

L=.85, G=1

Clear[f,g]

g[0] = .99

f[0] =.99

Do[g[i] = 1-1.3+1.3 f[i-1];
f[i] =    .85/1  g[i]+(1-.85) .99,{i,1,10}]

ans111=Table[{g[i],f[i-1]},{i,1,10}]

ans222=Table[{g[i],f[i]},{i,1,10}]

ans31=Flatten[Table[{ans111[[i]],ans222[[i]]},{i,1,Length[ans222]}],1]

ans41 = Join[{{.99,.99}},ans31]

steps2 = ListPlot[ans41,Joined->True]

eq2plot = Plot[(1+((x-1)/α))/.α->1.3,{x,.86,1},PlotStyle->Green]

optest = Plot[(((L x+G y0-L y0)/G))/.{L->.85 ,G->1,y0->.99},{x,.86,.99},PlotStyle->Red]

Show[steps2,eq2plot,optest]

So we could only go from 0.94 to 0.99 at this L/G.    

How does the result vary with α?

Clear[f,g]

Do[g[i] = 1-α+α f[i-1];
f[i] =    m  g[i]+(1-m) y0,{i,1,10}]

We can again see the infinite reflux result for α-=1.3, ~.86 and as L/G decreases, the separation degree decreases rapidly.

Plot[f[10]/.{α->1.3,f[0]->.99, y0->.99},{m,.7,1},PlotLegends->{"α->1.3"},AxesLabel->{"L/G","composition"},PlotRange->{.85,1}]

Plot[f[10]/.{α->1.2,f[0]->.99, y0->.99},{m,.7,1},PlotLegends->{"α->1.2"},AxesLabel->{"L/G","composition"},PlotStyle->Red]

Plot[f[10]/.{α->1.1,f[0]->.99, y0->.99},{m,.7,1},PlotLegends->{"α->1.1"},AxesLabel->{"L/G","composition"},PlotStyle->Green]

We see the substantial effect of reflux ratio and that the effect of α is most pronounced near infinite reflux ratio!

Show[%%%,%%,%]

Unfortunately, we can’t make any progress with a solution for “n”.  We see that for the 5th tray, all powers of n appear,

Expand[f[5]]

With the previous result recovered:

%/.m->1

So let’s declare victory with the result that “n” varies as 1/(α-1) or more specifically as ((Log[((-1+β)/(-1+x0))])/(Log[α]))!