ABSTRACT (original article): The present paper describes the arc-length continuation method; a mathematical technique that can solve a wide variety of chemical engineering problems. Several case studies involving distillation, reaction engineering, and thermodynamics are treated with the help of this method. The approach is powerful enough to elucidate situations when several turning points occur. Finally, the authors conclude with their experience teaching this technique. Indeed, this methodology was introduced to the graduate students of the nonlinear dynamics course at King Fahd University of Petroleum & Minerals (KFUPM). CITATION (original article): Binous, H., & Shaikh, A. A. (2014). Introduction of the arc-length continuation technique in the chemical engineering graduate program at KFUPM. In Computer Applications in Engineering Education. DOI: https://doi.org/10.1002/cae.21604
Batch rectification is a dynamic separation process in which the composition in the reboiler (still) evolves as the more volatile component is continuously removed in the distillate. This causes the still composition to become progressively leaner in the light component over time. Despite this inherent unsteadiness, many analytical and numerical studies replace the full dynamic description with a sequence of quasi-steady configurations. This approximation forms the basis of the pseudo steady-state assumption used in the present notebook. The key idea is that the internal column dynamics (vapor–liquid equilibrium on trays, mass transfer between stages, and hydraulic adjustments) occur on a much faster timescale than the depletion of the reboiler composition. As a result, at each instant, the column is assumed to be in internal equilibrium corresponding to the current operating conditions, while only the slow variables—such as overall holdup and composition in the still—evolve in time. This separation of timescales justifies treating the column as a sequence of steady-state problems rather than solving a full set of stiff differential-algebraic equations. This is the essence of the pseudo steady-state assumption.
Under this framework, the batch rectification column is described by a set of nonlinear algebraic equations representing stage-to-stage equilibrium, operating lines, and boundary conditions (reboiler and condenser). These equations define a solution manifold: F(u, R) = 0 where u represents stage compositions and R is the reflux ratio. Instead of time, the system is analyzed in parameter space. The reflux ratio is selected as the bifurcation parameter because it strongly controls separation intensity and can induce nonlinear phenomena such as multiplicity of steady states and fold bifurcations.
In summary, the pseudo steady-state assumption reduces a complex dynamic batch distillation process into a tractable sequence of equilibrium states, enabling bifurcation analysis using arc-length continuation. This approach provides deep insight into nonlinear operating regimes, stability limits, and sensitivity to reflux ratio, making it a powerful tool for both theoretical analysis and process design of batch rectification systems.
Two method are successfully compared:
1- the first method uses FindRoot and predictor / corrector iterative approach.
2- the second method is based on Ref. 1 and uses NDSolve in order to solve a system of Differential Algebraic Equations or DAEs.
Under this framework, the batch rectification column is described by a set of nonlinear algebraic equations representing stage-to-stage equilibrium, operating lines, and boundary conditions (reboiler and condenser). These equations define a solution manifold: F(u, R) = 0 where u represents stage compositions and R is the reflux ratio. Instead of time, the system is analyzed in parameter space. The reflux ratio is selected as the bifurcation parameter because it strongly controls separation intensity and can induce nonlinear phenomena such as multiplicity of steady states and fold bifurcations.
In summary, the pseudo steady-state assumption reduces a complex dynamic batch distillation process into a tractable sequence of equilibrium states, enabling bifurcation analysis using arc-length continuation. This approach provides deep insight into nonlinear operating regimes, stability limits, and sensitivity to reflux ratio, making it a powerful tool for both theoretical analysis and process design of batch rectification systems.
Two method are successfully compared:
1- the first method uses FindRoot and predictor / corrector iterative approach.
2- the second method is based on Ref. 1 and uses NDSolve in order to solve a system of Differential Algebraic Equations or DAEs.
Constant relative volatility binary mixture example
alpha=2.2;(*relativevolatility*)xB=0.40;(*stillcomposition*)Nstage=8;(*numberofrectifyingstages*)
Continuation parameters
ds=0.02;nSteps=120;
VLE equation for constant relative volatility mixture
yEq[x_]:=alpha*x/(1+(alpha-1)*x)
Rectifying operating line
yOp[x_,xD_,R_]:=(R/(R+1))*x+xD/(R+1)
Stage equations
stageEquations[vars_,R_]:=Module[{x,y,xD,eqs},x=vars[[1;;Nstage]];y=vars[[Nstage+1;;2Nstage]];xD=vars[[-1]];eqs={};(*equilibriumrelations*)Do[AppendTo[eqs,y[[i]]==yEq[x[[i]]]],{i,1,Nstage}];(*operatinglinerelations*)Do[AppendTo[eqs,y[[i]]==yOp[x[[i+1]],xD,R]],{i,1,Nstage-1}];(*condenser*)AppendTo[eqs,y[[Nstage]]==xD];(*reboilerboundary*)AppendTo[eqs,x[[1]]==xB];eqs]
Variables
vars=Join[Table[x[i],{i,Nstage}],Table[y[i],{i,Nstage}],{xD}];
Initial solution at reflux ration R=1
R0=1.0;guess=Join[Table[0.4+0.05i,{i,Nstage}],Table[0.5+0.05i,{i,Nstage}],{0.85}];sol0=FindRoot[stageEquations[vars,R0],Evaluate@Thread[{vars,guess}],MaxIterations->100];u0=vars/.sol0;
Print["Initial Distillate Composition = ",xD/.sol0];
Initial Distillate Composition = 0.787093
First Method based on FindRoot and predictor / corrector scheme
In[]:=
(*arc-lengthcontinuation*)curve={};AppendTo[curve,{R0,xD/.sol0}];(*tangentinitialization*)dR=1.0;tangent=Normalize[Join[ConstantArray[0.,Length[u0]],{dR}]];uPrev=u0;RPrev=R0;Do[(*PREDICTORSTEP*)pred=Join[uPrev,{RPrev}]+ds*tangent;uGuess=pred[[1;;-2]];RGuess=pred[[-1]];(*ARC-LENGTHSYSTEM*)augEquations[z_]:=Module[{u,RR,physEq,arcEq},u=z[[1;;-2]];RR=z[[-1]];physEq=stageEquations[u,RR]/.Equal->Subtract;arcEq=(z-pred).tangent;Join[physEq,{arcEq}]];zVars=Join[vars,{RR}];zGuess=Join[uGuess,{RGuess}];sol=Quiet@FindRoot[augEquations[zVars]==ConstantArray[0,Length[zVars]],Evaluate@Thread[{zVars,zGuess}],MaxIterations->100];zSol=zVars/.sol;uNew=zSol[[1;;-2]];RNew=zSol[[-1]];(*STORERESULTS*)xDnew=uNew[[-1]];AppendTo[curve,{RNew,xDnew}];Print["Step = ",step," R = ",NumberForm[RNew,4]," xD = ",NumberForm[xDnew,5]];(*UPDATETANGENT*)tangent=Normalize[Join[uNew,{RNew}]-Join[uPrev,{RPrev}]];uPrev=uNew;RPrev=RNew;,{step,1,nSteps}]
Step = 1 R = 1.02 xD = 0.79077
Step = 2 R = 1.038 xD = 0.79399
Step = 3 R = 1.055 xD = 0.7972
Step = 4 R = 1.073 xD = 0.80039
Step = 5 R = 1.09 xD = 0.80356
Step = 6 R = 1.108 xD = 0.80671
Step = 7 R = 1.126 xD = 0.80984
Step = 8 R = 1.143 xD = 0.81294
Step = 9 R = 1.16 xD = 0.81602
Step = 10 R = 1.178 xD = 0.81908
Step = 11 R = 1.195 xD = 0.82211
Step = 12 R = 1.213 xD = 0.82512
Step = 13 R = 1.23 xD = 0.8281
Step = 14 R = 1.247 xD = 0.83106
Step = 15 R = 1.265 xD = 0.83399
Step = 16 R = 1.282 xD = 0.8369
Step = 17 R = 1.299 xD = 0.83978
Step = 18 R = 1.316 xD = 0.84263
Step = 19 R = 1.334 xD = 0.84545
Step = 20 R = 1.351 xD = 0.84824
Step = 21 R = 1.368 xD = 0.85101
Step = 22 R = 1.385 xD = 0.85374
Step = 23 R = 1.402 xD = 0.85645
Step = 24 R = 1.419 xD = 0.85912
Step = 25 R = 1.436 xD = 0.86177
Step = 26 R = 1.453 xD = 0.86438
Step = 27 R = 1.47 xD = 0.86696
Step = 28 R = 1.487 xD = 0.86951
Step = 29 R = 1.504 xD = 0.87203
Step = 30 R = 1.521 xD = 0.87451
Step = 31 R = 1.538 xD = 0.87696
Step = 32 R = 1.555 xD = 0.87938
Step = 33 R = 1.572 xD = 0.88176
Step = 34 R = 1.589 xD = 0.88411
Step = 35 R = 1.606 xD = 0.88642
Step = 36 R = 1.623 xD = 0.8887
Step = 37 R = 1.64 xD = 0.89095
Step = 38 R = 1.657 xD = 0.89316
Step = 39 R = 1.675 xD = 0.89533
Step = 40 R = 1.692 xD = 0.89747
Step = 41 R = 1.709 xD = 0.89957
Step = 42 R = 1.726 xD = 0.90164
Step = 43 R = 1.743 xD = 0.90367
Bifurcation diagram
Second Method, thoroughly described in Ref. 1, is based on NDSolve and solving a system of DAEs
Both methods yield the same bifurcation diagram
References
References
[1] Binous, H., & Shaikh, A. A. (2014). Introduction of the arc-length continuation technique in the chemical engineering graduate program at KFUPM. In Computer Applications in Engineering Education (Vol. 23, Issue 3, pp. 344–351). Wiley. DOI: https://doi.org/10.1002/cae.21604
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Arc-Length Continuation in Chemical Engineering at KFUPM
by Housam Binous and Abdullah Shaikh
Wolfram Community, STAFF PICKS, May 26, 2026
https://community.wolfram.com/groups/-/m/t/3722703
by Housam Binous and Abdullah Shaikh
Wolfram Community, STAFF PICKS, May 26, 2026
https://community.wolfram.com/groups/-/m/t/3722703