The simulation is running on MatLab, and I have checked my code for all possible errors. The system used is cyclohexane-n-heptane-toluene. The column is running at total reflux.
Dr. Carlos I checked for the same, couldn't find it. Actually to give you a deeper insight, the vapour rates are continously decreasing and becoming negative, which in turn is affecting the composition. Then, it is decreasing too and the composition in the top tray also became negative, which is causing the Bubble-Point sub-routine to dysfunction and hence the appearance of NaN.
And Dr. Samuel, other people have successfully modelled the system using the odesolving method that I am employing!
Can you prove that your solutions should remain non-negative?
Without looking at the equations, I can't tell for sure, but it looks like once one of the component becomes negative, you create NaNs out of a square root or a log. Your problem is then to keep all your components non negative. Apart from the perhaps brutal NonNegative option in Matlab, you can play with the relative and the absolute tolerance (RelTol and AbsTol). RelTol controls the error in your large solution components and AbsTol controls the error of your small components. If you want to keep the solutions positive, then decrease AbsTol. For finite time spans, you should manage to keep your solutions non-negative.
If this doesn't work, go back to equation system. The equations might be correct, but there could be something wrong with the parameter values, or the initia conditions.
Differential equations involved in Batch Distillation are reputed to be stiff due to different dynamic in the boiler (low inertia) and plates (highly dynamic). Integrating these equations together causes some numerical problems. So the best way to avoid them is to use Gear's method a Backword Differential Formula implemented in LSODI package available in Fortran. May be there is an equivalent in matlab.