 Open Access
 Total Downloads : 422
 Authors : Divya Srivastava, Bikash Mohanty, Ravindra Bhargava
 Paper ID : IJERTV2IS120995
 Volume & Issue : Volume 02, Issue 12 (December 2013)
 Published (First Online): 24122013
 ISSN (Online) : 22780181
 Publisher Name : IJERT
 License: This work is licensed under a Creative Commons Attribution 4.0 International License
Modeling and Simulation of MEE System under Fouling Condition
Divya Srivastava1*, Bikash Mohanty2, Ravindra Bhargava2
1Chemical Engineering Department, Meerut Institute of Engineering and Technology, Meerut 250005, India
2Chemical Engineering Department, Indian Institute of Technology Roorkee, Roorkee 247667, India
Abstract
A generalized mathematical model has been developed to study the effect of fouling in Multiple Effect Evaporator (MEE) system used in Indian sugar industry. The developed model has been solved by globally convergent technique for one complete cycle of operation. The model can handle complexities such as exhaust steam inputs in more than one effect, vapor bleeding from desired effects, heat loss from each effect, variations in physical properties, overall heat transfer coefficient (OHTC) through external empirical correlations, and condensate flashing. For the computation of OHTC in the developed model, an empirical correlation, in terms of fouling resistance (Rf) has been developed. The result of the developed model has been validated against industrial data for exit liquor concentration, vapor body temperature and OHTC values. The prediction shows close agreement with the industrial data. Further, in order to improve the steam economy of the MEE system, flash vapors have also been incorporated in the developed model.
Keywords: Mathematical modeling, Globally convergent method, Fouling resistance, Overall heat transfer coefficient, Sugar industry, Multiple effect evaporator
Nomenclature:
Enthalpy of evaporation (kJ/kg K)
A Area of effect (m2)
Fraction of vapor stream for vapor bleed
Boolean constant
BPR Boiling point rise (C)
Bx Concentration of sugar cane juice (Bx)
Cp Specific heat capacity of sugarcane juice (kJ/kg K)
CS,C Condensate flow rate from steam (CS) and vapor (C) (kg/s)
Ff Feed flow rate (kg/s)
h Enthalpy of water at saturated temperature (kJ/kg K)
hl Enthalpy of liquor at saturated temperature (kJ/kg K)
Hs Enthalpy of vapor at superheated temperature (kJ/kg K)
Hv Enthalpy of vapor at saturated temperature (kJ/kg K)
L Mass flow rate of sugarcane juice (kg/s)
Q Rate of heat transfer (kJ/s)
Qloss Heat loss (kJ/s)
Rf Fouling resistance, m2C/kW
SHTS Superheated steam temperature

Temperature (C)
Tf Film temperature
TL Temperature of liquor (C)
Ts Steam temperature
Tv Vapor body temperature (C)

Overall heat transfer coefficient (OHTC) (kW/m2C)

Mass flow rate of vapor due to evaporation (kg/s)
VB Vapor bleed flow rate (kg/s)
VS Mass flow rate of exhaust steam (kg/s)
SVS Scaled mass flow rate of exhaust steam (kg/s)
STs Scaled steam temperature
Subscript
i Effect number
max Maximum
in At inlet condition
17 Effect number

Introduction
Evaporators are vital part of many process industries such as pulp and paper, sugar, caustic, food processing, desalination and other allied industries. The main aim of using Multiple Effect Evaporator (MEE) system is to reduce energy consumption required for evaporation. MEE system itself is a complex system and its complexity further increases if fouling phenomena is incorporated into it. Fouling is the deposition of organic and inorganic materials, present in the liquor, on the heat transfer surface, which further increases thermal resistance and reduces heat transfer coefficient. In all industrial MEE system, the fouling of heat transfer surfaces is a serious problem as it causes loss of production, and also affects the quality of product. The present work is focused on concentration of weak cane sugar juice in MEE system under fouled condition. In general, sugar industries employ feed forward arrangement
[1] and in India weak cane sugar juice is generally concentrated from 13Bx to 60Bx in MEE system.For economic production of sugar crystals, proper design and operation of MEE system is necessary. The cost of manufactured sugar depends mainly on the steam economy of the MEE system. Steam economy can be described as the ratio between total amount of vapor generated and the total amount of steam supplied. It has been found that different researchers [29] have developed models for simulation of MEE systems used in various chemical industries. However, various researchers [1015] have studied simulation of MEE system used in sugar industry but literature on simulation of MEE systems used in sugar industry under fouling condition is very rare. Further, it has also been observed that investigators such as Mwaba 2003[16], Mwaba et al. 2006 [17] and Yu 2003 [18] have tried to study the fouling of heat transfer surfaces in sugar industry. A scrutiny of literature suggests that, very little or almost negligible work is done on simulation of MEE system in sugar industry under fouled condition. Therefore, in the present paper, an attempt has been made to develop a mathematical model which incorporates fouling phenomena in MEE system used in Indian sugar industry.
1.1 Fouling of Heat Transfer Surfaces:
As stated above fouling is the deposition of organic and inorganic materials on the heat transfer surfaces. It is observed that, in cane sugar industry, major
deposit forming minerals are phosphates, sulphates, carbonates, silicates and oxides [1921]. Bruce 2002
[22] suggested that, these deposits accumulates on heat transfer surface or evaporator surface and results in the formation of deposits. Initially a particular type of fouling takes place and then other types dominate. Epstein 1983 [23] discussed various types of fouling such as crystallization fouling, particulate fouling, chemical reaction fouling, corrosion fouling and biological fouling in detail. Further, he also discussed that these fouling types are associated with different events such as initiation, transport, attachment, removal and aging. Sheikholeslami 2007 [24] studied these events associated with fouling in sugar industry in detail and suggested that composite fouling should be added as new fouling type and growth as the new event. Further, it has been suggested by Yu 2003 [18], Yu and Sheikholeslami 2009 [25] that composite fouling is the main mechanism of fouling in cane sugar evaporators. According to Sheikholeslami and Ng 2001 [26] and Yu 2003 [18] composite fouling is a form of fouling which involves more than one type of fouling species or mechanism. These different fouling species interact with each other to produce final deposit. 
Model development:
In the present work, a generalized steady state mathematical model, developed by authors [27, 28] is modified to transient state mathematical model to incorporate the effect of fouling. This is done by adding time, in terms of days as a variable parameter in the developed model. It should be noted that above stated model is capable of incorporating steam input in more than one effect, vapor bleeding from desired effect, heat loss from each effect, variations in physical properties such as boiling point rise (BPR), specific heat capacity and computation of overall heat transfer coefficient (OHTC) under clean condition. Details of this model can be obtained in Srivastava 2011 and Srivastava et al. 2013 [27, 28].
Model developed in the above referred paper is here named as MEE07 and it concentrates the cane sugar juice or weak liquor from 13Bx to 55Bx in a seven effect Rising film feed forward MEE system. The first 2 effects (Effectno. 1 & 2) are Semikestner type and all other effects (Effect no. 37) are of Robert type evaporator. Exhaust steam of superheated nature is introduced in two effects i.e. in effect no. 1 and effect no. 3. Further, it is also observed that vapor bleeding has been done from effect nos. 1, 2, 4, 5 and
6. These bled vapors are used in other parts of industry such as crystallization unit, preheaters etc. For the sake of clarity, schematic diagram of MEE 07 is reproduced in Figure 1.
VB1 VB2 VB4 VB VB6
Vapor from last effect
Steam
Steam
1 2 3 4 5 6 7
Steam
Product
L1 L2 L3
Feed
L4 L5
L6
Vapor stream Liquid stream
Condensate out
Fig. 1 Schematic diagram of MEE07
Industrial data has been collected and empirical correlation for the computation of OHTC (U) under clean condition has been developed for two different types of evaporator [27]. The two developed dimensionless correlations have been reproduced below and given as Eq. 1 and Eq. 2. Further, based on the data collected, empirical correlations are also developed for the computation of heat loss from each effect.
For Semikestner type evaporator:
U = 2.5665 Bx out 1.4553 Tf 0.052 Ff 0.0914
For the present MEE system, one complete cycle of operation consist of 10 days. Fig. 2 shows the trend of OHTC as a function of days for effect no. 1. It can be observed that OHTC values decreases from day 1 to day 10. Other effects also follow the same trend. On the first few days of the cycle, heat transfer surface was clean, so OHTC value was maximum and then falls gradually up to the last day of cycle due to the increase in fouling resistance.
2.6 Fig. 2 OHTC profile from industry data for effect no. 1
Umax
(1)
Bx max
Tmax
Fmax
OHTC, kW/m2C
OHTC, kW/m2C
2.5
For Robert type evaporator:
U
Umax
= 0.1909 Bx out
Bx max
3.153
Tf Tmax
0.298
Ff Fmax
0.0615
2.4
(2)
Where, Umax = 5.71kW/m2C; Bxmax = 56Bx; Tmax = 122C; Fmax = 100kg/s
2.3
0 2 4 6 8 10
Days

Empirical correlation for OHTC of MEE07 under fouling condition
In order to study the effect of fouling, OHTC is computed for one complete cycle of operation. For this purpose, three different sets of industrial data have been collected.
Sugar cane juice contains various organic and inorganic matters which deposits on the heat transfer surface of all effects. As a result, heat transfer coefficient decreases gradually from day 1 to the last day of cycle. The decrease in OHTC value with days reduces the performance of MEE system which was measured in terms of the drop in exit liquor concentration and has been shown in fig. 3. After the
last day of cycle, the whole or partial MEE system needs cleaning of the heat transfer surface to bring it to day1or initial condition.
Further, it has also been observed that, in order to compensate the ill effects of fouling, the industry increases the exhaust steam temperature by 1 or 2C during the intermediate days of cycle of operation. This increase in temperature results in the
improvement of the exit liquor concentration. Table
Fig. 4 Fouling phenomena observed in effect no. 1 of MEE07
Fouling resistance, Rf, m2 0C/kW
Fouling resistance, Rf, m2 0C/kW
0.0120
0.0100
1 shows the trend of exhaust steam temperature (superheated) used by industry.
Days of
operating cycle
1
5
8
Exhaust steam temperature (In effect no. 1)
125C
130C
130C
Exhaust steam temperature (In effect no. 3)
120C
122C
125C
Days of
operating cycle
1
5
8
Exhaust steam temperature (In effect no. 1)
125C
130C
130C
Exhaust steam temperature (In effect no. 3)
120C
122C
125C
Table 1: exhaust steam temperature profile
0.0080
0.0060
0.0040
on 5th day exhaust steam temperature is changed
2 3 4 5 6 7 8
From Fig. 3, it can be seen that, when industry increases exhaust steam temperature on intermediate days i.e. on 5th and 8th day, a rise in exit liquor
concentration is achieved. The industry has adapted this trend so as to increase exit liquor concentration or to bring it at par with the initial exit liquor concentration. However, it has been observed that exit liquor concentration was instead decreasing with
days.
Further, it is also observed that the effect of increasing exhaust steam temperature is observed only up to the first four effects i.e. effect nos. 1 to 4 and last three effects i.e. effect nos. 5 to 7 does not show any changes. To tackle the above difficulty it was thought logical to model Rf with days of operation rather than OHTC for effect no.1 to 4 using a sigmoidal function and it has been given by Eq. 3, whereas, for effect nos. 5 to 7 conventional functional relationship has been used and it has been given by Eq. 4. Moreover, Eq. 3 consists of three constants a, b and c which is functions of steam temperature (Ts), and are computed by Eqs. 3a to 3c respectively.
Fig. 3 Exit liquor concentration obtained from industrial data
=
(3)
1+
56 Where,
Exit concentration, 0Bx
Exit concentration, 0Bx
= 1 + 2
+ 3 2 + 4 3 (3a)
54
= 1 + 2
+ 3 2 + 4 3 (3)
52
= 1 + 2
+ 3 2 + 4 3 (3c)
50
0 2 4 6 8 10 1
Days
An analysis of fouling phenomena associated with evaporators shows that once the steam temperature is increased, vapor body temperature increases and hence fouling resistance increases. Fouling resistance increases in sigmoidal fashion and achieves an
Eq. 4 is the conventional form of equation for computation of Rf in terms of concentration, film temperature, feed flow rate and time which has been expressed in days. Rf values, thus computed from Eqs. 3and 4 for different effects were now added to clean OHTC values obtained from Eqs. 1 and 2. In this way, OHTC values for different days of operation can be obtained.
asymptotic value [17] before the steam temperature
Rf
Bx c
Tf d
Ff e
is increased again. In the present section, efforts have = a daysb out
(4)
been made, to predict the values of fouling resistance (Rf) for different days of the operation of the cycle. Fig. 4 shows the above behavior in effect no.1 of MEE07. All these behavior have added complexity to the simulation.
Rf max
Bx max
Tmax
Fmax

Mathematical model under fouling condition:
In the present work, model developed under clean condition (MEE07) is modified to incorporate the effect of fouling. The incorporation of fouling converts the model from steady state to transient state. The time duration for study is one cycle of
operation and it is 10 days. Based on mass and energy balance over an effect (ith effect) of a MEE system, it is observed that two nonlinear equations
are required to describe an effect [28]. These generalized equations were reproduced here and are given by Eq. 5 and Eq. 6. As in present system, there are 7 effects and each effect is represented by two non linear equations, therefore the present model consists of 14 non linear equations. It should be noted that hese equations are in scaled form and the details of method of scaling is given in Srivastava et al. 2013 [28].
g i = (i SVsi1Lin i1 + 1 i 1 i1
SLi2Lin SLi1Lin i1 + SLi1Lin CPi1
STvi1Ts + BPRi1 SLi Lin Cpi STvi Ts + BPRi
SLi1Lin SLi Lin Hvi Qlossi)/ Lins (5)
g i + 1 = (Ui A1 STi1 Ts STvi Ts + BPRi
i SVsi1Lin i1 1 i 1 i1
SLi2Lin SLi1Lin i1)/ Lin s (6)
The governing equations were now solved for all the days of cycle of operation. While solving the set of equations for 10 days it has been kept in mind that for each day the heat transfer coefficient of all the seven effects of the MEE system changes due to fouling. The variables involved in these governing equations are categorized as specified and unknown variables and were given below:
Unknown variable: Vs1, Vs2, Tv1, Tv2, Tv3, Tv4, Tv5, Tv6, L1, L2, L3, L4, L5, L6
Specified variable: Lin, Cin, TLin, TS1,TS2, TV7, A1, A2, A3, A4, A5, A6, A7, 1, 2, 3, 4, 5, 6, 7, Lout, U1, U2, U3, U4, U5, U6, U7,Ts, SHTS1, SHTS2
To meet the above objectives, the values of OHTC as a function of days of operation and effects of the MEE system are supplied through Eq. 3 and 4. All the input variables for each day of operation were collected from industry and are supplied to computer code through a data file.
A computer program is developed using FORTRAN for simulation of the system. Microsoft Developer Studio is used to compile and run computer code. Main program utilizes thirteen subroutines and two functions. The set of non linear equations were solved using globally convergent method.

Solution Technique:
Globally convergent method has been used for solving set of non linear algebraic equations. It is a modified form of Newtons method. The main drawback with the Newtons method is that, if the initial guess supplied is not close enough to the root, then it may fail to converge. A global method is the one that converges to a solution irrespective of the starting point taken. Newtons method is modified by adding a globally convergent strategy [29], so that with each iteration some progress can be made in achieving solution. The initial guess values for each day of operation were taken from industry data. It should be noted that solver does not pose any convergence problem. Moreover, to check whether convergence was achieved or not, several convergence criteria such as change in function values and variable values were considered along with overall heat balance around all effects.

Results and Discussion:
Before using the model for simulation, it is necessary to test the validity of developed model. For this purpose, results obtained through simulation were validated against three sets of industrial data. Various parameters such as concentration profile, vapor body temperature and OHTC values of all effects obtained from the model were compared with that of industry data. Figure 5 and 6 shows the comparison between simulation results and industrial data for exit liquor concentration and vapor body temperature respectively. Predictions for exit liquor concentration and vapor body temperature were within an error band of Â± 0.8% and Â± 0.3% respectively.
Fig. 5 Comparison of final concentration obtained from industry for all 3 sets with that predicted by the model for one cycle of operation
58
Final concentration, 0Bx
Final concentration, 0Bx
56
54
Industry data (set 1) Model prediction (set1) Industry data (set 2)
Model prediction (set2)
Industry data (set 3) Model prediction (set3)
Industry data (set 1) Model prediction (set1) Industry data (set 2)
Model prediction (set2)
Industry data (set 3) Model prediction (set3)
52
50
48
0 2 4 6 8 10 12
Days
Similarly industrial OHTC and predicted OHTC values were compared and it is observed that the predictions are within a maximum error band of
Â±15%.
Fig. 6(a) Comparison of vapor body temperature obtained from industry with that predicted by the model for effect no. to 4 for one cycle of operation
110
Vapor body temperature, 0C
Vapor body temperature, 0C
105
100
economy has been computed. While computing steam economy, it has been kept in mind, that the final concentration of the liquor is maintained equal to that of without flash vapors. Fig. 7 shows the comparison between steam economy of the MEE system computed by model with and without flash vapors.
95
Effect no. 1 Industry data
Effect no. 1 Model prediction
90 Effect no. 2 Industry data
Effect no. 2 Model prediction
Effect no. 3 Industry data
85 Effect no. 3 Model prediction Effect no. 4 Industry data
Effect no. 4 Model prediction
80
0 2 4 6 8 10 12
Days
Fig. 7 Variation of steam economy of MEE system during the complete cycle with and without flash
2.40
Steam Economy
Steam Economy
2.35
2.30
Fig. 6(b) Comparison of vapor body temperature obtained from industry with that predicted by the model for
effect no. 5 to 7 for one cycle of operation
100
Vapor body temperature, 0C
Vapor body temperature, 0C
80
2.25
2.20
2.15
Steam economy without flash
Steam economy with flash
0 2 4 6 8 10 12
Days
60
Effect no. 5 Industry data
40 Effect no. 5 Model prediction Effect no. 6 Industry data Effect no. 6 Model prediction
Effect no. 7 Industry data
Effect no. 7 Model prediction
20
0 2 4 6 8 10 12
Days
From the above discussion it is clear that developed model is capable of predicting industry data within a acceptable error band. After testing the reliability of the developed model, it is then further used for the betterment of the MEE system.
3.1 Simulation of MEE system with condensate flash
It is a well known fact that fouling deteriorates the heat transfer coefficient of an effect and erodes the evaporation capacity of the MEE system. The consumption of exhaust steam in the system decreases which results in decrease in final concentration of the liquor and the whole MEE system works under poor steam economy. Further, it is also observed that industry does not make use of flash vapors. In this section an attempt has been made to increase steam economy by incorporating primary and secondary flash vapors. Primary and secondary flash tanks are added in MEE07 and the governing equations, based on mass and energy balance, which includes primary and secondary flash vapors were developed and used in the computer code. The governing equations, for all the seven effects with condensate flash vapors were now solved for each day of cycle of operation and steam
From the figure (Fig. 7) it has been clear that one gets an improved steam economy throughout the period of the cycle. The minimum and maximum increase in steam economy is of the order of 8.55%
for 1st day and 8.18% for 10th day respectively.

Conclusion
In the present work, a generalized steady state model, developed by authors for simulation of MEE system used in Indian sugar industry is converted to transient state model to study the effect of fouling. One complete cycle of operation of this industry consists of 10 days. The developed model is capable of handling various complexities of sugar industry such as exhaust steam (saturated/superheated) inputs in more than one effect, vapor bleeding from desired effects, heat loss from each effect, and variations in boiling point rise as well as specific heat capacity, heat transfer coefficient through external empirical correlations, and condensate flashing under clean as well as fouling condition. The developed model predicts exit liquor concentration and vapor body temperature of all effects within an error band f
Â±0.8% and Â±0.3% respectively. Further, an empirical correlation has been developed for the prediction of OHTC of all effects for each day of operation. The model predicts the OHTC values for one complete cycle of operation for all effects within a maximum error band of Â±15%. Based on the accuracy of the results obtained, it can be concluded that developed model can be used for simulation of MEE system. Further, flash vapors are also incorporated in the MEE system to increase the steam economy of the system and simulation results show a minimum and
maximum increase in steam economy of 8.18% and 8.55% respectively.
References

Kern, D.Q. (1950). Process Heat Transfer. McGraw Hill.

Bhargava, R. (2004). Simulation of flat falling film evaporator network. PhD Thesis, Indian Institute of Technology Roorkee, India.

Bhargava, R., Khanam, S., Mohanty, B., & Ray, A.K. (2008a). Selection of optimal feed flow sequence for a multiple effect evaporator system. Computers & Chemical Engineering. 32(10), 22032216.

Bhargava, R., Khanam, S., Mohanty, B., & Ray, A.K. (2008b). Simulation of flat falling film evaporator system for concentration of black liquor. Computers & Chemical Engineering. 32(12), 3213 3223.

Bhargava, R., Khanam, S., Mohanty, B., & Ray, A.K. (2008c). Mathematical model for a multiple effect evaporator system with condensate, feed and product flash and steam splitting. Indian Journal of Chemical Technology. 15(2), 118129.

Bremford, D.J., & MullerSteinhagen, H. (1994). Multiple effect evaporator performance for black liquorI simulation of steady state operation for different evaporator arrangements. APPITA Journal. 47(4), 320326.

ElDessouky, H.T., & Ettouney, H.M. (1999). Multipleeffect evaporation desalination systems: thermal analysis. Desalination. 125(13), 259276.

Khanam, S., & Mohanty, B. (2007a). A Process Integration based Approach for the Analysis of Evaporator System. Chemical Engineering and Technology. 30(12), 16591665.

Khanam, S., & Mohanty, B. (2010). Energy reduction schemes for multiple effect evaporator systems. Applied Energy. 87(4), 11021111.

Jorge, L.M.M., Righetto, A.R., Polli, P.A., Santos, O.A.A., & Maciel, F.R. (2010). Simulation and analysis of a sugarcane juice evaporation system. Journal of Food Engineering. 99(3), 351359.

Kaya D, and Sarac HI. (2007). Mathematical modeling of multipleeffect evaporators and energy economy. Energy, 32 (8)153642.

Lewis, A.E., Khodabocus, F., Dhokun. V., & Khalife. M. (2010). Thermodynamic
simulation and evaluation of sugar refinery evaporators using steady state modeling approach. Applied Thermal Engineering. 30(1415), 21802186.

Pedro, C.C. Pimenta & Sebastiao, F.D.A. (1993). Simulation of multiple effect evaporators. Computers and Chemical Engineering. 17, S153S158.

Radovic, L.R., Tasic, A.Z., Grozanic, D.K., Djordjevic, B.D., & Valent, V.J. (1979). Computer design and analysis of operation of a multiple effect evaporator system in the sugar industry. Industrial and Engineering Chemistry Product Research and Development. 18(2), 318323.

Solberg, D.T. (2004). Online heat transfer measurement and analysis for sugar mill evaporators, M. Sc thesis, Louisiana State University,

Mwaba, M.G. (2003). Analysis of heat exchanger fouling in cane sugar industry, Ph.D Thesis, Eindhoven University.

Mwaba, M.G., Golriz, M.R., & Gu. J. (2006). A semi empirical correlation for crystallization fouling on heat exchange surfaces. Applied Thermal Engineering. 26(4), 440447.

Yu, H. (2003). The mechanisms of composite fouling in Australian sugar mill evaporators by calcium oxalate and amorphous silica, PhD Thesis, The University of New South Wales.

Bajpayee, A. (1992). Corrosion and its prevention in cane sugar industry. PhD Thesis, National Sugar Institute, India.

Honig, P. (1963). Principles of sugar technology. Elsevier Publishing Company.

Walford, S.N., & Walthew, D.C. (1996). Preliminary model for oxalate formation in evaporator scale, in: Proceedings of South African Sugar Technologist Association. 70, 231235.

Bruce, S. (2002). Scale deposit problems in pulp and paper mills. In Proceedings of the African Pulp and Paper Week, International Convention Centre. Durban.

Epstein, N. (1983). Thinking about heat transfer fouling: A 5×5 matrix. Heat Transfer Engineering. 4(1), 4356.

Sheikholeslami, R. (2007). Fouling in membranes and thermal unitsunified approach to its principles, assessment, control and mitigation, Desalinations Publications, First edition.

Yu, H. & Sheikholeslami. R. (2009). Composite fouling on heat exchange surface in Australian sugar mill evaporator. Heat transfer Engineering. 30(13), 1033 1040.

Sheikholeslami, R., & Ng, M. (2001). Calcium sulfate precipitation in the presence of nondominant calcium carbonate: Thermodynamics and kinetics, Industrial and Engineering Chemistry Research. 40(16), 35703578.

Srivastava, D. (2011). Simulation of multiple effect evaporators with and without fouling. PhD thesis, Indian Institute of Technology Roorkee, India

Srivastava, D., Mohanty, B., & Bhargava,
R. (2013). Modeling and Simulation of MEE system used in sugar industry. Chemical Engineering Communication. 200 (8), 10891101.

Press, W.H., Teukolsky, S.A., Vellerling, W.T., & Flannery, B.P. (1992). Numerical Recipies in FORTRANThe Art of Scientific Computing, Second edition, Cambridge University Press,.