McCabe-Thiele “stepping” near a pinch point

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

mjm@nd.edu
chemeprof.com

A (potentially interesting) "aside" which shows a relation between McCabe - Thiele calculations near a pinch point and some simple mathematics which has been used to describe intermittency (a popular buzz-word from chaos theory which does describe some real phenomena).

The conclusion is that the number of stages near a pinch point will vary as Sqrt[Rmin/(R-Rmin)], where R is the reflux ratio and Rmin is the minimum reflux ratio.

What happens to stage by stage stepping near a tangency pinch between the operating and equilibrium lines ???

We need some equations for operating and equilibrium lines.  If we note that any line, y= m x +b
could be transformed into Y = X if Y=y+b and X=mx, we can use Y = X to represent any operating line.  

If in general, the equilibrium line is curved, it could be represented as an expansion in powers of x as:

								y = a0 + a1 x + a2 x^2 + a3 x^3 + .......
								
we will truncate it at second order and confine our thinking to small values of x (this can be done with no loss of generality.

(*operating line*)

operate[x_ ]:=x

(*equilibrium curve *)

equilib[x_ ]:= ϵ+x+x^2

Plot[{operate[x],equilib[x]/.ϵ->0,
equilib[x]/.ϵ->.1,equilib[x]/.ϵ->-.1},{x,-.5,.5},
		PlotStyle->{ {  }, {Dashing[{.05,.05}]},
		{Dashing[{.02,.05}]},{Dashing[{.01,.03}]}},
		PlotLabel->{"Tangency when ϵ is 0"}]

Because we are stepping between the operating and equilibrium lines, it is clear that as eps gets closer to 0 (from the plus side), the number of stages will increase.  Can we figure out how much??

ϵ=.1;

calculate steps

numbs = NestList[equilib,-.5,8]

manuipulate the data to get the standard diagram

ans1 = Transpose[{numbs,numbs}]

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

now put it in a form to plot

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

attach one last point to the steps

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

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

ListPlot[steps,PlotJoined->True]

Plot[{equilib[x],operate[x]},{x,-.7,1.7}]

show the diagram

Show[%,%%]

We see that it takes about 8 stages to get past the pinch

ϵ=.05;

calculate steps

numbs = NestList[equilib,-.5,12]

manuipulate the data to get the standard diagram

ans1 = Transpose[{numbs,numbs}]

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

now put it in a form to plot

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

attach one last point to the steps

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

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

ListPlot[steps,PlotJoined->True]

Plot[{equilib[x],operate[x]},{x,-.7,1.7}]

show the diagram

Show[%,%%]

Now what if eps is 0.05

eps=0.05;

NestList[equilib,-.5,12]

Transpose[{%,%}]

Transpose[{ Take[%%,Length[%%]-1],
			Take[%%,{2,Length[%%]}]}]

Flatten[Table[{%%[[i]],%[[i]]},{i,1,Length[%]}],1]

%[[Length[%]]][[2]];

Join[%%,{{%,%}}]

ListPlot[%,PlotJoined->True]

Plot[{equilib[x],operate[x]},{x,-.7,1.7}]

Show[%,%%]

We might guess that N varies as 1/Sqrt[eps] !!!

How could we predict what will be the variation  ??

Note that what we are doing in the stepping process is following the discrete steps:

y = ϵ + x + x^2,
x=y

which is the same as

x(j+1) = x(j) + ϵ + x(j)^2

Note that x(j) is usually small compared to 1, it seems like we are are really dealing with a differential equation:

d x/dz = ϵ + x^2

we can find z, which is proportional to N by performing the integration,

Clear[x,ϵ];

z= Integrate[1/(ϵ+x^2),x]

zz = (z/.x->x1 )- (z/.x->x0)

N[zz/.{ϵ->.1,x0->-.7,x1->1.7}]

N[zz/.{ϵ->.05,x0->-.7,x1->1.7}]

N[zz/.{ϵ->.01,x0->-.7,x1->1.7}]

From which we can see that the differential equation predicts the behavior pretty well .

Over the entire range we get this behavior

Plot[zz/.{x0->-.7,x1->1.7},{ϵ,.5,.01}]

If we do a log - log plot we can see that there is a power relation between N and eps

LogLogPlot[zz/.{x0->-.7,x1->1.7},{ϵ,.2,.01}]

We have guessed that it will be ϵ^-1/2 from the form of the integral, let's see !!

ParametricPlot[{ϵ^-(1/2),zz/.{x0->-.7,x1->1.7}},
				{ϵ,.001,.1}, AspectRatio->.5]

We get a line, so that this is the relation.

What can we do analytically ??

zz

We would like to know how this will vary as eps is varied.  We pick one of the end points as 0 with no loss of generality.

%/.x0->0

This expression looks like:

ArcTan[s]/Sqrt[ϵ]

We know that s will be very large  (i.e.go to infinity) if x1 >> eps so we have:

Limit[%,s->Infinity]

Which is the desired result, consistent with the behavior above:

N stages ~ 1/Sqrt[ϵ]

What does this tell us ??

If ϵ is really this small, you will probably not dare design the device to operate in this region.  
However, if you measure the existence of a tangent pinch, this gives you:

1.	 An idea of how close it really is and also

2.		How things will change if ϵ is varied.  

Note that the value of ϵ is effectively:  (r-rtangent)/rtangent where r is whatever parameter that you are varying that changes the distance between the equilibrium curve and the operating line.

Let's see what all of this has to do with intermittency

We define a function that tells how x1 will vary, using compile speeds things up

mapmcabe = Compile[{x1},eps+x1 + x1^2]

Good idea to initialize the array

Table[xx[i]=0,{i,0,1000}];

some values for the variables

eps=.01;

xx[1]=.01;

xc=.5;

calculate the x's, the second line keeps the result bounded

Do[{	xx[i] = mapmcabe[xx[i-1]],
		If[xx[i-1]>xc,{xx[i]=(1-xx[i])/(1-xc)}]	  
			},
		{i,2,200}]

Then plot the data, the irregular rapid jumps are known as intermittency,

The longest time periods between jumps scale as eps^(-1/2)

ListPlot[Table[xx[i],{i,1,200}],PlotJoined->True,
		PlotRange->{-1,1}]

eps=.002;

Do[{	xx[i] = mapmcabe[xx[i-1]],
		If[xx[i-1]>xc,{xx[i]=(1-xx[i])/(1-xc)}]	  
			},
		{i,2,200}]

ListPlot[Table[xx[i],{i,1,200}],PlotJoined->True,
PlotRange->{-1,1}]

We see that the time period is longer,


if we use more points you can see that the character is unchanged

Do[{	xx[i] = mapmcabe[xx[i-1]],
		If[xx[i-1]>xc,{xx[i]=(1-xx[i])/(1-xc)}]	  
			},
		{i,2,1000}]

ListPlot[Table[xx[i],{i,1,1000}],PlotJoined->True,
PlotRange->{-1,1}]