Scope and Limitations of the Mathematical Models Developed for the Forward Feed Multi-Effect Distillation Process—A Review

Abstract

Desalination has become one of the obvious solutions for the global water crisis due to affording high-quality water from seawater and brackish water resources. As a result, there are continuing efforts being made to improve desalination technologies, especially the one producing high-quantity freshwater, i.e., thermal desalination. This improvement must be accomplished via enhancement of process design through optimization which is implicitly dependent on providing a generic process model. Due to the scarcity of a comprehensive review paper for modeling multi-effect distillation (MED) process, this topic is becoming more important. Therefore, this paper intends to capture the evolution of modeling the forward feed MED (most common type) and shed a light on its branches of steady-state and dynamic modeling. The maturity of the models developed for MED will be thoroughly reviewed to clarify the general efforts made highlighting the advantages and disadvantages. Depending on the outputs of this review, the requirements of process development and emerging challengeable matters of modeling will be specified. This, in turn, would afford a possible improvement strategy to gain a reliable and sustainable thermal desalination process.

1. Introduction

Various thermal desalination methods, such as multi-effect distillation (MED) and multi-stage flash (MSF), have been utilized to produce freshwater from brine water. The MED is the oldest large-scale seawater industrial desalination method that is referred to multi-effect boiling (MEB) or multi-effect evaporation (MEE). MED is the utmost effective thermal distillation process thermodynamically, even though it stands in the 2nd place after MSF in the thermal desalination market [1].

Compared to MSF, MED has several advantages, including running the desalination process at low temperatures and normally combined with thermal vapor compression (TVC) unit (external steam provider). Moreover, MED is characterized with high production of freshwater and high heat transfer efficiency. MED also utilizes lower pumping energy compared to MSF [2]. Due to a lower top brine temperature in the MED process, its thermal energy consumption is lower than MSF process. Thus, for a large-scale seawater desalination, MED is a more favorable technology.

TVC is generally used to decrease the external steam subjected to the 1st stage as it uses some of the vapor produced from the last stage. The vapor is compressed in the TVC to produce heat for preheating the upcoming feed. In other words, TVC recovers some of the produced vapor of the last stage as a motive steam, which in turn lowers the heat required for the 1st stage. Therefore, TVC has a positive impact as it increases the system performance ratio ($PR$) besides improving the economic feasibility. In other words, MED normally runs at low top brine temperature in the range 60–70 °C compared to MSF [3]. Thus, the unit capacity of the MED system is much lower than that of the multi-stage flash (MSF) process. Thus, MED has less production capacity compared to MSF. Moreover, the TVC unit steam ejector recovers energy from the low-pressure steam of the last stage of evaporators [4]. Therefore, steam ejectors reduce steam consumption and introduces low-operating pressure conditions. This, in turn, would mitigate the fouling and scaling propensities. Accordingly, multi-effect distillation - thermal vapor compression (MED-TVC) system was deemed most effective in lowering energy consumption, superior efficiency and free from environmental issues compared to MSF [5].

There is a progressive research made to improve the design of MED system as the wider technology of seawater desalination. However, MED design is critically dependent on the establishment of a robust mathematical model. Basically, modeling of an industrial process is a crucial task in correlating the process response pointers with the control variables and a set of fixed parameters. This would help to forecast the variation of process indicators against the expected variation of control variables via simulation. Moreover, accurate models can be successfully used to investigate the proper values of operating conditions to maintain the process at an optimal operation.

Several published studies of the MED system are engrossed on explanation of system features and operational performance. However, studies relating to process modeling progress with thorough analysis are not readily found in the literature. Therefore, the following sections are review on the modeling of MED process based on its two branches of the steady-state and dynamic modeling. The advantages and shortcoming of these models are also depicted in a detailed table and that is followed by some suggestions of methods of improvement and plausible current research direction.

2. Description of Multi-Effect Distillation and Thermal Vapor Compression System

The forward feed MED process can be characterized as an evaporation and condensation process. However, it is originally operated at low pressure that correspond to a decrease in both temperature and pressure from the 1st stage to the last one. The scheme of forward feed MED is basically characterized by supplying the feed seawater to the 1st stage to be evaporated and concentrated to form the brine stream of the next stage; therefore, the high-concentration brine is fed to the last stage. This is a well-established configuration of MED process compared to the parallel one as it works with similar feed flow rate for each stage. Figure 1 demonstrates the three most important parts of the forward feed MED system, including the stages (falling-film evaporators) in series, condenser, and thermo compressor (steam jet ejector). The feed seawater reaches the down condenser to condense part of the vapor leaving the last stage and to preheat the seawater before supplying to TVC with exiting the non-condensable gases (NCG). At the same time, the distillate of the last stage is directed to a demister into two streams, i.e., one stream is directed to the condenser (feed reheater) and then flows to the fresh water line while the remnants portion of vapor is directed to TVC to recompressed and de-superheated to produce a motive stream (superheated steam).

The feed water is sprayed or dispersed in a thin film on the surface of evaporator tubes to boil and evaporate, thereby generating additional vapor. Therefore, it is fair to expect that the disposed brine salinity for the 1st stage is the lowest compared to the last stage. Moreover, the high-concentration stream of the corresponding stage is collected in the brine pool (brine flash box) where it flashes and generates more vapor.

The first stage uses heat extracted from steam from TVC to preheated seawater and form small quantity of water vapor that is deployed to produce heat for the 2nd stage. The generated vapor from the 1st stage discharges its latent heat in the 2nd stage and condensate is formed inside the tubes. Similarly, latent heat is discharged, and a marginal vapor is created in the 2nd stage, which supplies heat to the 3^rd stage. This procedure is repeated for the rest of stages with progressively falling pressure and temperature until the vapor temperature approaches the feed seawater temperature [3]. Hence, the feed is boiled in a series of stages with no demand to provide further heat after the 1st stage. The feed water evaporates outside the tubes and the vapor is conveyed into the tubes of the next stage having lower boiling point temperature and pressure. Here the vapor condenses and vaporizes more feed water [6]. Three pumps are used in this configuration to feed seawater into the last condenser, to reject the brine into sea, and to deliver the product water into product tank.

3. Overview of Existing Models of MED System

3.1. Steady-State Modeling of MED System

Mathematical modeling of any industrial process with advanced computer simulations can deliver a clear insight into the operation of the process and would aid to thoroughly perceive the process conditions to forecast the performance metrics besides optimizing the efficacy of the system.

Exact system modeling is crucial for creating knowledge and awareness for investigating conceivable outcomes for process’s enhancement. Therefore, more than a few MED models have been established. The state-of-the-art of MED steady-state modeling and its evolution is described in brief:

El-Sayed and Silver [8] built the original model for forward feed MED based on several thermodynamic assumptions. The model was able to calculate some process parameters for performance evaluation. Specifically, a formula to examine the ($P_{\mathrm{R}}$) (−) was developed as depicted in Equation (1). Despite its simplicity, the equation gives a speedy and very solid approach to evaluate the performance of the process under existing working conditions; it cannot be utilized for optimization or sensitivity analysis.

$$P_{\mathrm{R}}=\frac{h_{\mathrm{s}}}{\frac{h_{\mathrm{fg}}}{N}+\frac{q_{\mathrm{m},\mathrm{f}}}{q_{\mathrm{m},\mathrm{d}}} c_p \left(\Delta T_{\mathrm{fh}}+\left(\mathsf{\Delta }T\right)\right)+\frac{N-1}{2N} c_p\Delta T_{\mathrm{e}}},$$

(1)

where h_s (kJ/kg) is the enthalpy of vaporization of steam, $N$ the number of stages, and $q_{\mathrm{m},\mathrm{f}}$ and $q_{\mathrm{m},\mathrm{d}}$ (kg/s) are the mass flow rates of feed seawater and distillate, respectively. $c_p$ is the mean specific heat capacity at a fixed pressure (kJ/(kg K)), $\Delta T_{\mathrm{fh}}$ (K) the temperature change between the 1st stage and the feed at exit of the last feed heater, $\sum \left(\mathsf{\Delta }T\right)$ is the sum of boiling point elevation ($\Delta T_{\mathrm{BPE}}$) (K) and temperature variation as a result to pressure loss which is equal to (Σ(ΔT) = $\Delta T_{\mathrm{BPE}}$ + ΔT_Δp), and $\Delta T_{\mathrm{e}}$ (K) the temperature change between two stages. Notwithstanding its uncomplicatedness, Equation (1) is obtained via intensive thermodynamic point of views and is suitable for fast approaching of the $P_{\mathrm{R}}$ (−) and needed transfer areas for FF-MED system under recognized control variables. Nevertheless, it cannot be used to obtain comprehensive information concerning a variety of certain streams or to comprehend system sensitivities to different input parameters. More importantly, the hypothesis of this model is originated from the assumption of some thermodynamic simplifications which include constant fluid properties of specific heat capacity, $\Delta T_{\mathrm{BPE}}$ and latent heat.

An additional equation is given for determining the needed heat transfer surface area as a function of a recognized or presumed overall heat transfer coefficient. Equation (2) was generally used to assess the area demanded for the transfer of heat (Q) (${\mathrm{k}}_{\mathrm{W}}$) between two fluids separated by a thin wall.

$$Q=hA\Delta T_{\mathrm{eff}}.$$

(2)

The total surface area through which the heat is transferred, and an overall surface heat transfer coefficient are ($A$) (m²) and ($h$) (kW/(m² K)), respectively. In addition, $\Delta T_{\mathrm{eff}}$ ((°C) is the effective temperature difference that denotes the logarithmic mean temperature difference ($L_{\mathrm{MTD}}$) (ΔT_LM) (°C) which is normally utilized to design the condenser surface. Basically, this parameter can be used for fluids in both counter flow and parallel flow as given in Equation (3). Moreover, it specified the case of fluids constant temperature in condensation and evaporation.

$$\mathsf{\Delta }{T}_{\mathrm{LM}}=\frac{\left(\Delta T_{\mathrm{a}}-T_{\mathrm{b}}\right)}{\ln \left(\frac{\Delta T_{\mathrm{a}}}{\Delta T_{\mathrm{b}}}\right)},$$

(3)

where $a$ and $b$ identify exchanger ends. $\Delta T_{\mathrm{a}}$ and $\Delta T_{\mathrm{b}}$ (°C) are the temperature difference at $a$ and $b$, respectively.

Hanbury [9] and Minnich [10] introduced several presumptions related to the thermal load to introduce a simplified model for the MEE. In this regard, the following assumptions were considered:

• fixed seawater flow rate for all stages,
• fixed thermal loads in the 2nd and following stages,
• fixed surface heat transfer coefficient for all the stages,
• fixed heat transfer area for all the stages,
• fixed vapor formed in all the stages,
• fixed temperature-drop between the considered stages,
• fixed increase of temperature across the preheaters,
• fixed latent heat of evaporation in all the stages, and
• fixed physical properties of the seawater, distillate, and brine.

The heat supplied to the system to increase the temperature of the feed seawater from $\left(T_{\mathrm{s}}\right)$ to top brine temperature $\left(T_1\right)$ and to boil the vapor in the 1st stage is calculated from

$$Q_{\mathrm{s}}=q_{\mathrm{m},\mathrm{f}}{c}_p\left(T_1-T_{\mathrm{s}}\right)+\frac{q_{\mathrm{m},\mathrm{d}}{\lambda }_1}{N},$$

(4)

where $c_p\left(\mathrm{kJ}/\left(\mathrm{kg}\mathrm{K}\right)\right)$ and ${\lambda }_1$ (kJ/kg) are the specific heat capacity at fixed pressure and latent heat of evaporation, respectively.

The associated $P_{\mathrm{R}}$ (−) is interrelated to the number of stages ($N$) via Equation (5)

$$P_{\mathrm{R}}=N+\frac{q_{\mathrm{m},\mathrm{d}}\lambda }{q_{m,f}{c}_p\left(T_1-T_{\mathrm{p}}\right)},$$

(5)

where $T_p$ (°C) is the temperature of preheater of the 1st stage.

El-Dessouky et al. [11], El-Dessouky and Ettouney [12], and El-Dessouky et al. [13] developed one of the most comprehensive and detailed mathematical model for MED system that involved several sets of heat transfer coefficients, pressure drop and thermodynamic properties. Basically, this model has been built on the previous simple models of El-Dessouky and Assassa [14] and Darwish and El-Hadik [15] (the details are mentioned in

Table 1. Summary of literature review of modeling the forward feed MED process.

No.	Author and Year	Assumptions	The Specific Features and Advantages	Disadvantage
1	El-Sayed and Silver (1980)	Steady-state operation. The simplifying assumptions of mean latent heat, mean specific heat, and mean boiling point rise in an ideal solution are maintained.	Explicit correlations for evaluating the $PR$ and transfer areas were derived. The pressure drop due to friction is modeled in terms of a mean saturation temperature-drop augmenting the stage of boiling point elevation.	Several thermodynamic assumptions were made. The problem of corrosion has an impact on the allowable highest vapor release temperature. It cannot be used to carry out sensitive analysis to several parameters. It cannot be used for optimization or sensitivity analysis.
2	El-Dessouky and Assassa (1985)	Steady-state operation. A fixed heat transfer area in all the stages. Fixed temperature-drop between the stages.	The salinity of the brine in the first stage is lower than the salinity of the feed seawater, permitting for a higher operating temperature.	Very simplified model. The fouling inside and outside the tubes has not been considered.
3	Darwish and El-Hadik (1986)	Steady-state operation. Identical latent heat in all stages, and equivalent specific heat for the brine, distillate, and seawater feed. The temperature-drops between the stages are identical Identical amount of distillate generated from each stage.	Accounted the total the energy consumption. Accounted the calculation of water production cost.	TVC has not been considered. Linear temperature profile is assumed.
4	Hanbury (1995) and Minnich (1995)	Steady-state operation. No feed preheaters. Ignored contribution of distillate flash. The feed flow rate is equal for all stages. Constant temperature difference.	High temperature would produce high heat transfer coefficient.	Very simplified model. The fouling inside and outside the tubes has not been considered.
5	Aly and Marwan in (1997)	Dynamic operation. Full condensation of vapor inside the stage was assumed.	The model was able to capture the transient behavior, including start-up, shutdown, troubleshooting, and load changes.	Ignored the presence of non-condensable gases in the stages, which negatively affect the measurement of heat transfer areas in the evaporators.
6	El-Dessouky et al.(1998)	Steady-state operation. Fixed heat transfer areas in the evaporators and feed preheaters in all stages.	The influence of preheaters and flashing boxes were considered. Considered the influence of vapor leak through the venting system. Consider the distinction in thermodynamic losses from one stage to another. The influence of the physical characteristics of seawater on temperature and salinity was realized. The effect of non-condensable gases on the heat transfer coefficients in the evaporators and the feed preheaters was appraised.	All the feed must be heated to the boiling temperature before boiling commences. New equations for the heat transfer coefficients in evaporators and pre-heaters were developed that considered the relationship between temperature and fouling.
7	El-Dessouky and Ettouney (1999)	Steady-state operation. Ignored heat losses to the surroundings, and salt-free product. No flashed distillate. No feed pre-heating.	The model satisfies the analysis of various configurations of MED process of different number of stages.	Very simplified model that does not require a numerical solver.
8	El-Dessouky et al. (2000) [34]	Steady-state operation. The distillate is salt-free. Fixed and identical heat transfer areas in all stages.	The heat transfer area in each evaporator was modeled by the heat transfer equations as the summation of the area for brine heating and the area for evaporation Explore the influence of the existence of non-condensable gases on the heat transfer coefficients in the evaporators and condenser.	The area considered is the summation of the area of brine heating and the area for evaporation. The influence of temperature on the specific heat transfer area is more obvious at a high number of stages.
9	El-Dessouky and Ettouney (2002)	Steady-state operation. Constant specific heat at an average temperature. Constant thermodynamic losses across all the stages. Steam was employed to evaporate distillate in the 1st stage, not for heating the feed. No vapor flashes in the stages. The temperature change in the stages is identical to the change in condensation and evaporation temperatures. No energy loses to the environment.	High heat-transfer rates. Efficient disengagement of vapors and non-condensable gases. Suitable venting of the non-condensable gases, and simple controlling of fouling and scaling.	Simplified model does not tackle with the actual specification of practical desalination plants. Flashing boxes and preheaters are excluded. Non-condensable gases were not included in mass balance and transport correlations. However, it was used to assess the heat transfer coefficient. Feed enters the 1st stage at its saturation temperature, and, therefore, the steam is not deployed to heat the seawater feed but only to evaporate the distillate.
10	Aly and El Figi (2003)	Steady-state operation. Salt free of the product vapor. Ignoring the heat losses to the environment.	The performance ratio is fully related to the number of stages and marginally related to the top brine temperature. Recognize the influence of the operating conditions on the plant performance and produced fresh water.	Simplified model that ignored the fouling factors on the inside and outside pipes of evaporators and condenser.
11	Dardour et al. (2005)	Dynamic operation. Full condensation of vapor inside the stage was assumed.	A control method was developed for MED process coupled to a high temperature nuclear reactor. The model was already used to optimize the inlet parameter.	Ignored the presence of non-condensable gases in the stages, which negatively impacted the measurement of heat transfer areas in the evaporators.
12	Darwish et al. (2006)	Steady-state operation. Identical vapor generated from each stage. Identical boiling temperature difference between the stages. Identical specific heat for the brine and feed. Identical latent heat and boiling point elevation.	The proposed terms deployed to investigate the efficacy of thermal desalination processes (specific consumed fuel, specific consumed energy, and specific equivalent work) are useful to assess the actual performance of the MED system.	Intensive assumptions of equal mass flow rate of vapor generated in each stage.
13	Al-Sahali and Ettouney (2007)	Steady-state operation. Identical vapor generated from each stage. Identical boiling temperature difference between the stages. Identical specific heat for the feed seawater and brine. Identical latent heat and boiling point elevation.	MED desalination process is modeled with emphasis on process design, energy consumption, and economic facets. Rise of seawater temperature to its saturation one of each stage would retain a high-pitched heat transfer coefficient.	Fixed values of the temperature-drop, heat transfer coefficients, and specific heat. Identical heat transfer area in all the stages.
14	Darwish and Abdulrahim (2008)	Steady-state operation. Equal vapor generated from each stage. Identical boiling temperature change between the stages.	Low values of specific energy consumption and steam temperature are necessitated to activate the MED process. Increasing the number of stages would decrease the temperature and pressure in each stage.	Identical specific heat for the seawater feed exit brine. Identical latent heat and boiling point elevation.
15	Roca et al. (2008)	Dynamic operation. Full condensation of vapor inside the stage.	The model has relaxed the assumption of fixed heat transfer coefficient. The model was able to assess the thermal dynamics of the heater and the freshwater production rate.	Ignored the presence of non-condensable gases in the stages, which negatively impact the measurement of heat transfer areas in the evaporators.
16	Ameri et al. (2009)	Steady-state operation. Identical vapor generated from each stage. Identical boiling temperature difference between the stages. Identical specific heat for the exit brine and feed seawater. Identical latent heat and BPE.	Investigated several design parameters of stages, including number of stages, supplied steam pressure, temperature change of the stages, and feed seawater temperature, on MED process. To segregate the desalination system from environmental variations, entered feed flow rate to the condenser must be lessened under winter season.	MED is an intensive energy consumption that requires to deploy renewable energy systems to ensure the reliability of this technological choice.
17	Kamali et al. (2009)	Steady-state operation. Identical vapor produced from each stage. Identical boiling temperature difference between the stages. Identical specific heat for the seawater feed and exit brine. Identical latent heat and boiling point elevation.	The impact of all input parameters on heat transfer coefficients, pressure and temperature, and performance ratio were appraised based on the operating conditions and design parameters.	Simplified model based on a stringent set of approximation.
18	Khademi et al. (2009)	Steady-state operation. Identical vapor generated from each stage. Identical boiling temperature difference between the stages. Identical specific heat for the seawater feed and exit brine. Identical latent heat and boiling point elevation.	Detailed correlations of the heat exchanger, flash tank, and preheater were developed. The model was deployed to find the optimal values of the input parameters and empower the performance indicators.	The fouling factor has not been implicitly appeared in the model equations.
19	Yılmaz et al. (2012)	Steady-state operation. Identical vapor generated from each stage. Identical boiling temperature difference between the stages. Identical specific heat for the seawater feed and exit brine. Identical latent heat and BPE.	A comprehensive model for forward feed MED operated by renewable energy source was developed. Boiling point elevation and NEA were appraised. Overall heat transfer coefficients of evaporators and condenser are evaluated in terms of stage’s vapor temperature.	The exergy analysis for various compartments of MED process was not engaged in the model.
20	Palenzuela et al. (2012)	Steady-state operation. Identical temperature change across the stages. Identical temperature change across the preheaters. Salt-free distillate from all the stages.	Tow observed equations for the overall heat transfer coefficients in the 1st stage and the preheaters were developed for the first time via experimental operation.	The condensate that exits each stage is judged as saturated liquid. Negligible heat losses. Exergy analysis and thermodynamic losses are not provided in the model.
21	Mistry et al. (2013)	Steady-state operation. Zero salinity of produced water, Vapor exiting a stage passes by a stage preheater to preheat the feed seawater. Vapor and condensate exiting the preheater is subjected to the next stage as a heat source to produce vapor via discharging the remainder of the latent heat of condensation.	The intensive features of the model permit for suitable sensitivities to important variables corresponding to input, operating, and design conditions important for operate in hybrid systems optimization process. Applying of Engineering Equation Solver (EES) and JACOBIAN has simplified and mitigated the requirements assumptions and number of approximations.	Ignored the energy losses to the environment. Assuming of average overall heat transfer coefficient.
22	Druetta et al. (2013)	Steady-state operation. A constant value of the specific heat capacities was considered for each stream. Ignored the pressure drop in the demister and throughout the vapor condensation. Fixed boiling point elevation.	A feasible package to calculate the heat transfer area, energy consumption, and water production cost. Specific equations were presented to estimate the heat transfer areas in evaporators and preheaters and condenser. A complementary optimization study to optimize MED system was associated based on the model developed; therefore, the optimum input parameters and effect of size were provided.	Neglecting thermodynamic losses. No exergy analysis was presented by the model; therefore, the calculation of total energy consumption is not precisely accurate.
23	Mazini et al. (2014)	Dynamic operation. Full condensation of vapor inside the stage.	Full description of input parameters throughout the transient period is illustrated.	Ignored the presence of non-condensable gases in the stages, which negatively affect the measurement of heat transfer areas in the evaporators.
24	Tahir et al. (2015)	Steady-state operation. A constant value of the specific heat capacities was considered for each one of the streams. Ignored the pressure drop in the demister and throughout the vapor condensation.	Overall heat transfer coefficient for each evaporator and preheater are estimated based on two comprehensive equations that evaluated the fouling factor inside the tubes of condenser and preheater and outside the tubes of the evaporators.	The exergy analysis was not been analyzed for various compartment of MED process. BPE was reflected as a fixed value. 1 °C was presented as the temperature losses in piping and during condensation.
25	Azimibavil and Dehkordi (2016)	Dynamic operation. Full condensation of vapor inside the stage.	Improved results regarding the specification of the fluid flow distribution with signifying the cause of high fouling throughout the transient operation and how long it resists. Full description of input parameters throughout the transient period is illustrated. Partial differential equations were entailed to examine the variation of input parameters on both x and y directions.	The sensitivity analysis does not elucidate substantial information of the performance of a practical desalination MED plant. Ignored the presence of non-condensable gases in the stages, which negatively impact the measurement of heat transfer areas in the evaporators. Non-condensable gases have not been considered that negatively impact the estimation of heat transfer areas in the evaporators.
26	Cipollina et al. (2017)	Dynamic operation. Full condensation of vapor inside the stage.	The presence of non-condensable gases in the stages has been considered which has a positive role to measure the heat transfer areas in the evaporators. Complete condensation of vapor inside the effect was not assumed where a partial condensation was replaced.	The exergy analysis was not been analyzed for various compartment of MED process.
27	Filippini et al. (2018)	Steady-state operation. Steam flowrate and temperature are assumed to be known from an upstream process while freshwater production is evaluated. The heat transfer areas for the evaporators and feed preheaters in all stages is equal, according to industrial reality.	Comprehensive model was presented that is easily connected in a hybrid system with Reverse Osmosis process. Comprehensive thermodynamic correlations were applied to evaluate all the relevant thermodynamic characteristics of the MED process in terms of water salinity, temperature, and fouling.	An elevated top brine temperature can cause significant fouling, consequently decreasing the heat exchange coefficient. Freshwater demand and top brine temperature are fixed, while the steam flowrate and properties are evaluated. The exergy analysis was not been analyzed for various compartment of MED process.
28	Khalid at el (2018)	Steady-state operation. The vapor is salt-free since the mist eliminators do not permit brine to be entrained with the created vapor. Fixed specific heat due to a small variation of temperature. Seawater feed enters the 1st stage has 5° of sub cool.	Comprehensive model that relaxed several previous assumptions used to be used in the literature.	Evaporator tubes are totally wetted by sprayed seawater, without exit dry patches. The exergy analysis was not been analyzed for various compartment of MED process.

). They broke down various MED designs comprising the parallel stream, the parallel/cross flow and systems joined with a TVC or mechanical vapor compressor (MVC). They developed a generic model of the MEE to analyze the influence of the operating conditions and design parameters on the freshwater production cost. The model was originally established due to steady-state mass and heat balances combined with the heat transfer rate correlations for any stage. Moreover, the corresponding equations are linked to the $P_{\mathrm{R}}(-)$. However, fixed heat transfer areas for the feed preheaters and evaporators in all stages are assumed. It is noteworthy to mention that the heat transfer formulas utilized in the model used the area examined as the area for evaporation plus the area of brine heating. The results showed that the thermal $P_{\mathrm{R}}$ $(-)$ of the TVC and specific power consumption of the MVC declines at higher top brine temperatures. The $P_{\mathrm{R}}$ $(-)$ is almost independent of top brine temperature and is roughly influenced by the number of stages. In addition, the specific heat transfer area would decline due to raising the temperature of the heating steam. The conversion ratio $(C_{\mathrm{R}}$) ($-$) is shown to be subjected to the brine flow design and it is free from the vapor compression mode. More importantly, they provided specific design relationships to define the influence of top brine temperature and number of stages on the $P_{\mathrm{R}}(-)$, the specific flow rate of cooling water and the specific heat transfer area. The most important correlations of El-Dessouky et al. are as follows:

The total material balance correlations were derived to allocate the brine flow rate exiting stage N, ($q_{\mathrm{m},\mathrm{b}}$) (kg/s) and the feed flow rate ($q_{\mathrm{m},\mathrm{f}}$) (kg/s) as described below

$$q_{\mathrm{m},\mathrm{f}}=q_{\mathrm{m},\mathrm{d}}+q_{\mathrm{m},\mathrm{b}},$$

(6)

$$q_{\mathrm{m},\mathrm{b}}=\frac{q_{\mathrm{m},\mathrm{f}}{X}_{\mathrm{f}}}{X_{\mathrm{b}}}$$

(7)

Substituting Equations (6) into (7), and eliminating $q_{\mathrm{m},\mathrm{f}}$, gives

$$q_{\mathrm{m},\mathrm{b}}=(\frac{X_{\mathrm{f}}}{X_{\mathrm{b}}-X_{\mathrm{f}}})q_{\mathrm{m},\mathrm{d}},$$

(8)

where $f,b,d$ denotes feed, brine, and distillate mass flow rates, respectively. $X$ denotes salinity in ppm. The mass flow rate of the motive superheated steam of the 1st stage $q_{\mathrm{m},\mathrm{s}}$ (kg/s) is estimated based on $q_{\mathrm{m},\mathrm{f}}$ (kg/s), temperature difference between the 1st and 2nd stages, the mass flow rate of vapor and latent heat of the 1st stage ($q_{m,d1}\left(\mathrm{kg}/\mathrm{s}\right), {\lambda }_1$ (kJ/kg)) and corresponding heat loss to the environment $Q$ (W).

$$q_{\mathrm{m},\mathrm{s}}{\lambda }_{\mathrm{s}}=q_{\mathrm{m},\mathrm{f}}{c}_p\left(T_1-T_2\right)+q_{m,d1}{\lambda }_1+Q.$$

(9)

The latent heat $\lambda$ (kJ/kg) depends on the boiling temperature and estimated from

$$\lambda =2589.583+0.9156 T-4.8343\times {10}^{-2} T^2.$$

(10)

Due to presumption of equal thermal loads in all the stages, $q_{\mathrm{m},\mathrm{s}}$ (kg/s) is obtained from

$$q_{\mathrm{m},\mathrm{s}}={\frac{q_{\mathrm{m},\mathrm{d}1}{\lambda }_{\mathrm{v}1}}{{\lambda }_{\mathrm{s}}}}_{ }.$$

(11)

The boiling temperature of the 1st stage is estimated as

$$T_1=T_{v1}+{\left(\Delta T_{BPE1}\right)}_{ }+\Delta T_h.$$

(12)

The $\Delta T_{\mathrm{BPE}}$ (°C) is related to seawater concentration as presented in Equation (13)

$$\Delta T_{BPE}=X\left(B+CX\right)\times {10}^{-3},$$

(13)

where $B$ and $C$ (−) are constants related to brine temperature as follows:

$$B={(6.71+6.43\times {10}^{-2}{T}_b+9.74\times {10}^{-5}{T}_b^2)}^2\times {10}^{-3},$$

(14)

$$C={(22.238+9.59\times {10}^{-3}{T}_b+9.72\times {10}^{-5}{T}_b^2)}^2\times {10}^{-8}.$$

(15)

$T_{\mathrm{v}1}$ (°C) is the saturation temperature of the generated vapor of the 1st stage which is lower than $T_1$ (°C) by the $\Delta T_{\mathrm{BPE}1}$. Specifically, $\Delta T_{\mathrm{BPE}}$ is the growth in the boiling temperature as a result to the dissolved salts in water at a considered pressure. $\Delta T_{\mathrm{h}}$ (°C) is the hydrostatic head.

Equation (16) used to calculate the surface heat transfer coefficient for feed seawater inside the tubes ($h_{\mathrm{i}})$ (W/m² °C)

$$h_i=\frac{(\left(3293.5+T\left(84.24-0.171T\right)-{\rm X}\left(8.471+0.116{\rm X}+0.271T\right)\right)}{\left({\delta }_i/0.17272{)}^{0.2}\right) {(0.656v_v)}^{0.8}\left(\frac{{\delta }_i}{{\delta }_o}\right)},$$

(16)

where $T$ is the temperature (°C). ${\delta }_{\mathrm{i}}$ and ${\delta }_{\mathrm{o}}$ are the inner and outer tube diameters in (m).

The pressure-drop in the demister ($\Delta P_{\mathrm{p}})$ (kPa) is obtained from

$$\Delta P_p=9.583\times {10}^4{({\rho }_{ })}^{1.597}{(v_v)}^{0.710}{(L_p)}^{1.388},$$

(17)

where ${\rho }_{}$ (kg/m³), $v_v$ (m/s), and $L_{\mathrm{p}}$ (m) are mass density of demister, vapor velocity, and length of demister, respectively.

The seawater density ($\rho$) in (kg/m³) is obtained as

$$\rho =1000 \left(A_1{F}_1+A_2{F}_2+A_3{F}_3+A_{4 }{F}_4\right).$$

(18)

The following parameters ($A_1, A_2, A_3, A_4$) are calculated from

$$A_1=4.032219G_1+0.115313G_2+3.26\times {10}^{-4}{G}_3,$$

(19)

$$A_2=-0.108199G_1+1.571\times {10}^{-3}{G}_2-4.23\times {10}^{-4}{G}_3,$$

(20)

$$A_3=-0.012247G_1+1.74\times {10}^{-3}{G}_2-9\times {10}^{-6}{G}_3,$$

(21)

$$A_4=6.92\times {10}^{-4}{G}_1-8.7\times {10}^{-5}{G}_2-5.3\times {10}^{-5}{G}_3,$$

(22)

$$B=\frac{\left(2w_S/1000-150\right)}{150},$$

(23)

$$G_1=0.5, G_2=B, G_3=2B^2-1,$$

(24)

$$A=\frac{\left(\left(2T\right)-200\right)}{160},$$

(25)

$$F_1=0.5,F_2=A,F_3=2A^2-1,F_4=4A^3-3A.$$

(26)

The dynamic viscosity ($\mu$) in ($\frac{\mathrm{kg}}{\mathrm{m}}\mathrm{s}$) is calculated as

$$\mu ={\mu }_{\mathrm{W}}{\mu }_{\mathrm{R}}\times {10}^{-3},$$

(27)

$$\ln \left({\mu }_{\mathrm{W}}\right)=-3.79418+604.129/\left(139.18+T\right),$$

(28)

$${\mu }_{\mathrm{R}}=1+A{w}_{\mathrm{S}}+B{w}_{\mathrm{S}}{}^2,$$

(29)

$$A=1.474\times {10}^{-3}+1.5\times {10}^{-5} T-3.927\times {10}^{-8}{T}^2,$$

(30)

$$B=1.0734\times {10}^{-5}-8.5\times {10}^{-8} T+2.23\times {10}^{-10 }{T}^2.$$

(31)

Water salinity ($w_{\mathrm{S}}$) is the mass fraction of salt measured in (g/kg).

The thermal conductivity of seawater ($\kappa$) in ($\frac{\mathrm{W}}{\mathrm{m}}\mathrm{K}$), is acquired from

$${\log }_{10}\left(\kappa \right)={\log }_{10}(240+2\times {10}^{-4}{w}_{\mathrm{S}})+0.434\left(2.3-\frac{343.5+0.037S_{\mathrm{S}}}{T+273.15}\right){\left(1-\frac{T+273.15}{647.3+0.03w_{\mathrm{S}}}\right)}^{\frac{1}{3}}.$$

(32)

The seawater specific heat capacity at fixed pressure ($c_p$) (kJ/(kg K) is depicted in Equation (33) and is related to water salinity and temperature.

$$c_p=\left(A+BT+C{T}^2+D{T}^3\right)\times {10}^{-3}.$$

(33)

The four constants (A, B, C, and D) are linked to water salinity and calculated from

$$A=4206.8-6.6197w_s+1.2288\times {10}^{-2}{w}_{\mathrm{S}}{ }^2,$$

(34)

$$B=-1.1262+5.4178\times {10}^{-2}{w}_{\mathrm{S}}-2.2719\times {10}^{-4}{w}_{\mathrm{S}}{ }^2,$$

(35)

$$C=1.2026\times {10}^{-2}-5.3566\times {10}^{-4}{w}_{\mathrm{S}}+1.8906\times {10}^{-6}{w}_{\mathrm{S}}{ }^2,$$

(36)

$$D=6.8777\times {10}^{-7}+1.517\times {10}^{-6}{w}_{\mathrm{S}}-4.4268\times {10}^{-9}{w}_{\mathrm{S}}{ }^2.$$

(37)

The model assumes fixed thermal loads in all the stages,

$$Q_1=Q_2=\dots =Q_{\mathrm{N}-1}=Q_{\mathrm{N}}.$$

(38)

Therefore, for the 1st stage and stages from 2 to N,

$$Q_1=q_{\mathrm{m},\mathrm{s}}{\lambda }_{\mathrm{s}},$$

(39)

$$Q_{\mathrm{i}}=q_{\mathrm{m},\mathrm{di}}{\lambda }_{\mathrm{v}i }.$$

(40)

${\lambda }_{\mathrm{v}i}$ (J/kg) is the latent heat of generated vapor at ($T_{\mathrm{i}}-\Delta T_{\mathrm{loss}})$. Equation (40) can be elucidated in the form of Equation (41) to calculate the generated vapor of any stage

$$q_{\mathrm{m},\mathrm{di}}{\mathsf{\lambda }}_{\mathrm{i}}=q_{\mathrm{m},\mathrm{di}-1}{\lambda }_{\mathrm{i}-1}.$$

(41)

$q_{\mathrm{m},\mathrm{di}}$ (kg/s) is the distillate mass flow rate in the corresponding stage.

Moreover, the model presumes equal heat transfer areas for all the stages. Equations (42) and (43) are deployed to obtain the heat transfer area in the 1st and 2nd stages to N, respectively.

$$A_1=\frac{q_{\mathrm{m},\mathrm{d}1}{\lambda }_{\mathrm{v}1,}}{h_1\left(T_{\mathrm{s}}-T_1\right)},$$

(42)

$$A_{\mathrm{i}}=\frac{q_{\mathrm{m},\mathrm{di} }{\mathsf{\lambda }}_{\mathrm{i}}}{h_{\mathrm{i} }\left(T_{\mathrm{i}}-\Delta T_{\mathrm{loss}}\right)}.$$

(43)

The heat transfer area of each condenser $A_{\mathrm{c}}$ (m²) is estimated from

$$A_{\mathrm{c}}=\frac{Q_c}{h_{\mathrm{c} }{\mathsf{\Delta }\mathrm{T}}_{\left(\mathrm{LM}\right)}c}.$$

(44)

$\Delta T_{\mathrm{loss}}$ (°C) denotes the thermodynamic losses in the considered stage which varies between 0.5 to 3 °C, and $h$ (kW/(m² °C)) is the overall heat transfer coefficient. In this regard, the total specific heat transfer area of the whole desalination plant ($s_{\mathrm{A}}$) (m² s/kg) is defined as

$$s_{\mathrm{A}}=\frac{{\sum }_{i=1}^N{A}_{\mathrm{i}}+A_{\mathrm{c}}}{q_{\mathrm{m},\mathrm{d}}}.$$

(45)

The heat transfer surface area of the 1st stage $A_{\mathrm{ev}1}$ (m²) is predicted based on the heat transferred $Q_{\mathrm{ev}1}$ (${\mathrm{k}}_{\mathrm{W}}$), the overall surface heat transfer coefficient $h_{\mathrm{ev}1}$ (kW/(m² °C)) and the difference between steam temperature $T_{\mathrm{s}}$ (°C) and the boiling brine as depicted in Equation (46).

$$A_{\mathrm{ev}1}=\frac{Q_{\mathrm{ev}1}}{h_{\mathrm{ev}1}\left(T_{\mathrm{s}}-T_1\right)}.$$

(46)

Moreover, the thermal load of each stage can be obtained as

$$Q_{\mathrm{i}}=A_{\mathrm{i}}{h}_{\mathrm{i}}\Delta T_{\mathrm{i}}.$$

(47)

Due to the assumption of equal thermal load and heat transfer area in all the stages, therefore

$$\frac{Q_1}{A_1}=\frac{Q_2}{A_2}=\dots =\frac{Q_{N-1}}{A_{\mathrm{N}-1}}=\frac{Q_{\mathrm{N}}}{A_{\mathrm{N}}} ,$$

(48)

$$h_1\Delta T_1=h_2\Delta T_2=\dots =h_{\mathrm{N}-1}\Delta T_{\mathrm{N}-1}=h_{\mathrm{N}}\Delta T_{\mathrm{N}}.$$

(49)

$\Delta T$ (°C) and $h$ (kW/(m² °C)) are the temperature driving force and the overall heat transfer coefficient, respectively. $N$ (−) is the number of stages. Finally, the thermal performance $P_{\mathrm{R}}(-)$ is the most important metrics of the MED desalination system exhibited as

$$P_{\mathrm{R}}=\frac{q_{\mathrm{m},\mathrm{d}}}{q_{\mathrm{m},\mathrm{s}}}.$$

(50)

To summarize, El-Dessouky and Ettouney [16] built up a basic model that requires no numerical solver, and it can anticipate some performance parameters. However, their nonlinear reliance on operative parameters is lost due to many simplifying suppositions, including constant fluid properties, constant thermal load in each stage, no flashed distillate, and no feed pre-heating. It is accepted that the feed seawater goes into the 1st stage at the main stage’s saturation temperature. In other words, steam is utilized just to vaporize distillate, not for heating the feed.

Aly and El-Figi [17] established a steady-state model to examine the performance of the forward feed MED process and came up with the fact that the $P_{\mathrm{R}}(-)$ is meaningfully relevant to number of stages despite of the TBT (top brine temperature). The model is associated with different evaporator areas as well as fixed and equal temperature-drops between stages.

The model developed has been built based on several assumptions:

• salt-free of the generated product water,
• insignificant losses of mass and heat in the vacuum system, and
• no heat losses to the environment.

The model was established to estimate the flow rate, temperature, and heat transfer area of each stage, preheaters, as well as the last condenser.

The mass flow rate of distillate of each stage ($q_{\mathrm{m},\mathrm{d}\left(i\right)}$) (kg/s) and the brine flow rate of the 1st stage were derived based on the material, mass, and energy balance:

$$q_{\mathrm{m},\mathrm{d}\left(i\right)}=\frac{[q_{\mathrm{m},\mathrm{s}\left(\mathrm{i}\right)}{\lambda }_{\left(\mathrm{i}-1\right)}+q_{\mathrm{m},\mathrm{b}\left(\mathrm{i}-1\right)}\left(h_{\left(\mathrm{i}-1\right)}-h_{\mathrm{i}}\right]}{h_{\mathrm{g}\left(\mathrm{i}\right)}-h_{\mathrm{b}\left(\mathrm{i}\right)}},$$

(51)

$$q_{\mathrm{m},\mathrm{b}1}=q_{\mathrm{m},\mathrm{f}}-q_{\mathrm{m},\mathrm{d}1},$$

(52)

where $h$ (kW/(m² K)) is the heat transfer coefficient. The mass flow rate of steam of 1st stage ($q_{\mathrm{m},\mathrm{so}}$) and for each stage ($q_{\mathrm{m},\mathrm{si}}$) (kg/s) are

$$q_{\mathrm{m},\mathrm{so}}=\frac{q_{\mathrm{m},\mathrm{s}1}{\lambda }_s+q_{\mathrm{m},\mathrm{f}}{c}_p\left(T_1-t_1\right)}{{\lambda }_{\mathrm{s}}},$$

(53)

$$q_{\mathrm{m},\mathrm{si}}=q_{\mathrm{m},\mathrm{d}\left(\mathrm{i}-1\right)}-q_{\mathrm{m},\mathrm{d}\left(\mathrm{i}\right)}+q_{\mathrm{m},\mathrm{b}\left(\mathrm{i}\right)}.$$

(54)

The mass flow rate of steam ($q_{\mathrm{m},\mathrm{s}}$) (kg/s) for each preheater used to heat the brine can be calculated from

$$q_{\mathrm{m},\mathrm{s}}=\frac{q_{\mathrm{m},\mathrm{f}}{c}_p\left(T_{\mathrm{i}}-T_{\mathrm{i}+1}\right)}{{\lambda }_{\mathrm{i}}}.$$

(55)

The heat transfer area for each preheater ($A_{\mathrm{P}}$) is calculated based on the $\Delta T_{\mathrm{LM}}$ (°C) as

$$A_{\mathrm{P}}=\frac{q_{\mathrm{m},\mathrm{f}} c_p\left(T_{\mathrm{i}}-T_{\mathrm{i}+1}\right)}{h_{\mathrm{P}\left(\mathrm{i}\right)} \mathsf{\Delta }{T}_{{\mathrm{LM}}_{\left(\mathrm{i}\right)}}}.$$

(56)

The mass flow rate of cooling seawater ($q_{\mathrm{m},\mathrm{cw}})$ (kg/s) in the end condenser is predicted from Equation (57):

$$q_{\mathrm{m},\mathrm{cw}}=\frac{q_{\mathrm{m},\mathrm{d}\left(\mathrm{N}\right)}{\lambda }_{\mathrm{N}} }{c_p\left(T_{\mathrm{co}}-T_{\mathrm{sw}}\right)-q_{\mathrm{m},\mathrm{f} }\mathsf{\Delta }{T}_{{\mathrm{LM}}_{\mathrm{c}}}} =\frac{T_{\mathrm{co}}-T_{\mathrm{sw}}}{\log \left\{\frac{T_{\mathrm{d}\left(\mathrm{N}\right)}-T_{\mathrm{sw}}}{T_{\mathrm{d}\left(\mathrm{N}\right)}-T_{\mathrm{co}}}\right\} },$$

(57)

where $T_{\mathrm{co}}$ and $T_{\mathrm{sw}}$ (°C) are the temperature of cooling water and seawater, respectively.

The heat transfer area of the end condenser ($A_{\mathrm{c}}$) (m²) is estimated as

$$A_{\mathrm{c}}=\frac{q_{\mathrm{m},\mathrm{cw}}{c}_p\left(T_{\mathrm{co}}-T_{\mathrm{sw}}\right)}{h_{\mathrm{c}}{\mathrm{LMTD}}_{\mathrm{c}}}.$$

(58)

Darwish et al. [18] and Darwish and Abdulrahim [5] built a simplified model for MED with assumptions, such as identical mass flow rate of vapor generated in each stage and debated the trade-off between heat transfer area and $P_{\mathrm{R}}(-)$. Moreover, the model had several thermodynamic limitations, such as constant fluid properties and constant heat exchange coefficients. A linear temperature profile is also assumed. With those presumptions, an equation for $P_{\mathrm{R}}(-)$ has been inferred.

$$P_{\mathrm{R}}=\frac{N}{\left(1+N \frac{q_{\mathrm{m},\mathrm{f} }{c}_p\left(T_1-T_{\mathrm{sw}}\right)}{{\dot{q}}_{\mathrm{m},\mathrm{d} }{\lambda }_{\mathrm{s}}}\right)}.$$

(59)

Al-Sahali and Ettouney [3] built a robust model for MED-TVC using a sequential solution approach, relatively to an iterative technique. They considered fixed heat transfer coefficients, specific heat transfer heat and temperature-drop. However, the model does have an innovative approach; it calculates the total temperature-drop ($\Delta T_{\mathrm{T}}$) (°C) from the 1st stage to the last stage via Equation (60)

$$\Delta T_{\mathrm{T}}=T_{\mathrm{s}}-\left(N-1\right) \Delta T_{\mathrm{l}}-T_{\mathrm{b}N},$$

(60)

where $\Delta T_{\mathrm{l}}$ (°C) and $T_{\mathrm{b}N}$ (°C) are temperature losses in each evaporation stage and the disposed brine temperature, respectively. For the 1st stage, the temperature-drop and heat transfer area in the 1st stage $A_{\mathrm{ev}1}$ are considered by Equations (61) and (62), respectively:

$$\Delta T_1=\frac{\Delta T_{\mathrm{T}}}{h_1{\sum }_{i=1}^N \frac{1}{h_{\mathrm{i}}} },$$

(61)

$$A_{\mathrm{ev}1}=\frac{ q_{\mathrm{m},\mathrm{d}1}{\lambda }_1}{h_{\mathrm{ev}1}\left(T_{\mathrm{s}}-T_1\right)}.$$

(62)

The heating mass flow rate of the 1st stage ($M_{\mathrm{s}}$) (kg/s) is

$$M_{\mathrm{s}}=\frac{ q_{\mathrm{m},\mathrm{d}1}{\lambda }_1}{{\lambda }_{\mathrm{s}}}.$$

(63)

The heat transfer area of the end condenser ($A_{\mathrm{c}}$) (m²) and the mass flow rate of cooling seawater ($q_{\mathrm{m},\mathrm{cw}})$ (kg/s) are predicted from Equations (64) and (65):

$$A_{\mathrm{c}}=\frac{ q_{\mathrm{m},\mathrm{dN}}{\mathsf{\lambda }}_{\mathrm{v}N}}{h_{\mathrm{c}}{L}_{{\mathrm{MTD}}_{\mathrm{c}}}},$$

(64)

$$q_{\mathrm{m},\mathrm{dN}}{\mathsf{\lambda }}_{\mathrm{v}N}=q_{\mathrm{m},\mathrm{f} }{q}_{\mathrm{m},\mathrm{cw}}{c}_p\left(T_{\mathrm{f}}-T_{\mathrm{cw}}\right).$$

(65)

Ameri et al. [19] examined the influence of design parameters on MED process conditions by developing a steady-state operational model. The model developed has implemented modifications in the TVC and demisters to improve the precision and model validation. Equations (66) and (67) were used to predict the inside heat transfer coefficient ($h_{\mathrm{i}}$) (W/(m² °C)) and outside heat transfer coefficient ($h_0$) (W/(m² °C)) of the condenser tubes:

$$h_{\mathrm{i}}=\frac{3293.5+T_{\mathrm{mean}}\left(84.24-0.1714T_{\mathrm{mean}}\right)-X\left(8.471+0.116 X+0.2716T_{\mathrm{mean}}\right)}{{(q_{\mathrm{m},\mathrm{di}}/0.17272)}^{0.2} {\left(0.65 v\right)}^{0.8}\left(q_{\mathrm{m},\mathrm{di}}/q_{\mathrm{m},\mathrm{d}0}\right)},$$

(66)

$$h_0=0.729 {[ \frac{g_{{\rho }_{\mathrm{l}}}\left({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{v}}\right) k_{\mathrm{l}}^3\acute{h_{\mathrm{s}}}}{N_{\mathrm{c}} {\mu }_{\mathrm{l}}{q}_{m,d}\left(T_{\mathrm{sal}}-T_{\mathrm{w}}\right)} ]}^{1/4},$$

(67)

$${\acute{h}}_{\mathrm{s}}=h_s+0.68c_{p.l} \left(T_{\mathrm{sat}}-T_{\mathrm{w}}\right),$$

(68)

where $T_{\mathrm{mean}}$ (°C) is the mean temperature inside the tube. $q_{m,d}, q_{\mathrm{m},\mathrm{di}}$, and $q_{\mathrm{m},\mathrm{do}}$ are the tube diameter (m), inner, and outer diameters (m) of the tube, respectively. $h_s$ is the enthalpy of steam (kJ/kg), ${\acute{h}}_{\mathrm{s}}$ (kJ/kg) is the specific enthalpy of evaporation. $N_{\mathrm{c}} (-)$, ${\kappa }_{\mathrm{l}} \left(\frac{\mathrm{W}}{\mathrm{m}\mathrm{K}}\right)$, ${\mu }_{\mathrm{l}} \left(\frac{\mathrm{kg}}{\mathrm{s}}\mathrm{m}\right)$, $T_{\mathrm{sat}}$ (°C), $v$ (m/s), and $g$ (m/s²) are the number of tubes in the condenser, thermal conductivity, dynamic viscosity of the liquid, saturation temperature, velocity and gravitational acceleration, respectively.

Kamali et al. [20] created a simple model of forward feed MED-TVC system to estimate the influence of operating conditions and design parameters on heat transfer coefficients, pressure, and temperature and $P_{\mathrm{R}}(-)$.

The temperature change between all the stages is presumed constant and calculated by

$$\Delta T=\frac{T_{\mathrm{s}}-T_{\mathrm{f}}}{N+1}.$$

(69)

The ratio of total feed flowrate ($q_{V,F1}$) (m³/h) to distillate flowrate ($q_{V,d1}$) (m³/h) of the 1st stage is calculated from

$$\frac{q_{V,F1}}{q_{V,d1}}=\frac{X_{\mathrm{b}1}}{X_{\mathrm{b}1}-X_{\mathrm{F}1}}.$$

(70)

Khademi et al. [21] established a steady-state simulation model for forward feed MED system and used it to simulate and optimize a six-effect evaporator. The model uses detailed equations of three main compartment of the heat exchanger, flash tank, and preheater based on familiar assumptions of MED.

The thermodynamic properties are calculated as follows:

• The $\Delta T_{\mathrm{BPE}}$ is calculated as follows:
  The $BPE$ correlation for the seawater is

  $$\begin{gathered} \Delta T_{\mathrm{BPE}}=[\frac{565.757}{T}-9.81559+1.54739\ln T-\left(\frac{337.178}{T}-6.41981+0.922743\ln T\right)A \\ +\left(32.681/T-0.55368+0.079022\ln T\right)\times A^2] \\ \times \left[A/\left(266919.6/T^2-379.669/T+0.334169\right)\right], \end{gathered}$$

  (71)

  $$A=\frac{19.819X}{1-X}.$$

  (72)
• The saturation temperature $T$ (°C) is

  $$T=\left(42.6776-\frac{3892.7}{\left[\ln \left(\frac{p}{100}\right)-9.48654\right]}\right)-273.15.$$

  (73)
• The saturation pressure p (kPa) is

  $$p=10.17246-0.6167302T+1.832249\times {10}^{-2}{T}^2-1.77376\times {10}^{-4}{T}^3+1.47068\times {10}^{-6}{T}^4.$$

  (74)
• The vapor enthalpy for pure water $H_v$ (kJ/kg) is

  $$H_v=2500.152+1.947T-1.94538\times {10}^{-3}{T}^2.$$

  (75)
• The liquid enthalpy for pure water $H_L$ (kJ/kg) is

  $$H_L=0.58021+4.1519T+3.536659\times {10}^{-4}{T}^2.$$

  (76)
• The latent heat of the water vapor $\lambda$ (kJ/kg) is

  $$\lambda =2589.583+0.9156T-4.8343\times {10}^{-2}{T}^2.$$

  (77)
• The enthalpy equation for the aqueous sodium chloride is

  $$H=A+BT+C{T}^2+D{T}^3+E{T}^4,$$

  (78)

  $$A=\left(0.0005+0.0378X-0.3682X^2-0.6529X^3+2.89X^4\right)\times {10}^3,$$

  (79)

  $$B=\left(4.145-4.974X+4.482X^2+18.31X^3-46.41X^4\right) ,$$

  (80)

  $$C=\left(0.0007-0.0059X+0.0854X^2-0.4961X^3+0.8255X^4\right) ,$$

  (81)

  $$D=\left(-0.0048+0.0639X-0.714X^2+3.273X^3-4.85X^4\right)\times {10}^{-3},$$

  (82)

  $$E=\left(0.02022-0.2432X+2.054X^2-8.211X^3+11.43X^4\right)\times {10}^{-6}.$$

  (83)

Moreover, the temperature-drop per each stage ($\Delta T_{\mathrm{i}}$) is reversibly related to the heat transfer coefficient as follows

$$\Delta T_{\mathrm{i}}=\left(T_{\mathrm{s}}-T_{\mathrm{cond}}\right)\frac{1/h_{\mathrm{i}}}{{\sum }_{i=1}^{N}1/h_{\mathrm{i}}}.$$

(84)

Yılmaz et al. [22] generated a detailed model for a forward feed MED seawater desalination system that operated using renewable energy sources. The usual basic assumptions were made in this model and the main contributions are:

• the change in the thermodynamic losses, including boiling point elevation and non-equilibrium allowance, for each stage is considered.
• the change of overall heat transfer coefficient of evaporators and condenser in terms of stage vapor temperature is considered.

Equation (85) is the correlation for the ${\left(N_{\mathrm{EA}}\right)}_{\mathrm{i}}$ for each stage

$${\left(N_{\mathrm{EA}}\right)}_{\mathrm{i}}=\frac{33 \Delta T_{\mathrm{i}}{}^{0.55}}{T_{\mathrm{v}i}},$$

(85)

where $\Delta T_{\mathrm{i}}$ (°C) is the temperature change of boiling brine between stages ($i-1$) and ($i$), as given in Equation (86). $T_{\mathrm{v}i}$ (°C) is the vapor temperature of stage $i$.

$$\Delta T_{\mathrm{i}}=T_{\mathrm{i}-1}-T_{\mathrm{i}}.$$

(86)

Palenzuela et al. [23] came up with a steady-state model considering a vertically stacked forward feed MED system situated at Platform Solar de Almería (PSA). The system is distinctly designed for distillate distribution to enhance recovered energy, which was comprehensively modeled. The authors for the first time found empirical relationships for the overall heat transfer coefficients in the 1st stage and the preheaters by executing a comprehensive experimental operation in the desalination system. Moreover, other equations were derived to calculate the overall heat transfer coefficient for the preheaters, the cooling seawater outlet temperature, and the vapor temperature of the 1st stage.

The thermal heat transferred from the hot water to seawater of the 1st evaporator of the 1st stage ($Q_h$) is drawn based on the overall heat transfer coefficient of the heater ($h_{\mathrm{h}}$) (kW/(m² °C)), the area of the evaporator ($A_{\mathrm{h}}$) (m²) and the inlet and outlet temperatures of the preheater ($(T_{\mathrm{h}.in}$) and ($T_{\mathrm{h}.out}$) (°C)).

$$Q_h=A_{\mathrm{h}}{h}_{\mathrm{h}} \frac{(T_{\mathrm{h}.in}-T_{\mathrm{b}\left(1\right)})-\left(T_{\mathrm{h}.out}-T_{\mathrm{f}}\right)}{\ln \left(\frac{T_{\mathrm{h}.in}-T_{\mathrm{b}\left(1\right)}}{T_{\mathrm{h}.out}-T_{\mathrm{f}}}\right)}.$$

(87)

In addition, it can be calculated based on the water specific heat ($c_p$) (kJ/(kg °C))

$$Q_{\mathrm{h}}=q_{\mathrm{m},\mathrm{h}}{c}_p \left(T_{\mathrm{h}.in}-T_{\mathrm{h}.out}\right)$$

(88)

The heat transfer for the preheater, given by Equation (89), is based on the heat transfer coefficient ($h_{\mathrm{p}}$) (kW/(m² °C)), the preheater area ($A_{\mathrm{p}}$) (m²), seawater temperature in the exit of the last preheater ($T_{\mathrm{f}}$) (°C) and seawater temperature in the inlet of the first preheater ($T_{\mathrm{p}\left(1\right)}$) (°C). $q_{\mathrm{m},\mathrm{h}}$ is the mass flowrate in the brine heater (kg/s).

$$Q_{\mathrm{p}\left(1\right)}=A_{\mathrm{p}}{h}_{\mathrm{p}}\frac{(T_{\mathrm{v}\left(1\right)}-T_{\mathrm{f}})-(T_{\mathrm{v}\left(1\right)}-T_{\mathrm{p}\left(1\right)})}{\ln \left(\frac{(T_{\mathrm{v}\left(1\right)}-T_{\mathrm{f}})}{(T_{\mathrm{v}\left(1\right)}-T_{\mathrm{p}\left(1\right)})}\right)}$$

(89)

In addition, $Q_{\mathrm{p}\left(1\right)}$ (kW) can be calculated from Equation (90) based on specific heat capacity of seawater at the inlet of first preheater $c_{p\left(1\right)}$ (kJ/(kg $K$)) and outlet of the last one $c_{pf}$ (kJ/(kg K)).

$$Q_{\mathrm{p}\left(1\right)}=q_{\mathrm{m},\mathrm{f}}{c}_{pf}{T}_{\mathrm{f}}-q_{\mathrm{m},\mathrm{f}}{c}_{p\left(1\right)}{T}_{\mathrm{p}\left(1\right)}.$$

(90)

The temperature change across the stages ($\Delta T_{\mathrm{v}}$) (°C) and the preheaters ($\Delta T_{\mathrm{p}}$) (°C) are estimated as

$$\Delta T_{\mathrm{v}}=\frac{T_{\mathrm{v}1}-T_{\mathrm{vN}}}{N-1},$$

(91)

$$\Delta T_{\mathrm{p}}=\frac{T_{\mathrm{f}}-T_{\mathrm{cw},out}}{N-1},$$

(92)

where $T_{\mathrm{cw},out}$ ($K$) is the seawater outlet temperature of the condenser. The heat transfer in the end condenser, ($Q_{\mathrm{c}}$) (kW), is obtained based on the overall heat transfer coefficient of the condenser ($h_{\mathrm{c}}$) (kW/(m² K)) and the area ($A_{\mathrm{c}}$) (m²):

$$Q_{\mathrm{c}}=A_{\mathrm{c}}{h}_{\mathrm{c}} \frac{(T_{\mathrm{v},N}-T_{\mathrm{cw},out})-\left(T_{\mathrm{v},N }-T_{\mathrm{cw}, in}\right)}{\ln \left(\frac{(T_{\mathrm{v},N}-T_{\mathrm{cw},out}}{T_{\mathrm{v},N }-T_{\mathrm{cw}, in}}\right)}.$$

(93)

Moreover, the heat transfer ($Q_{\mathrm{c}}$) (kW) can be calculated from

$$Q_{\mathrm{h}}=(q_{\mathrm{m},\mathrm{f}}+q_{\mathrm{m},\mathrm{CW}}) \left(h_{\mathrm{CW}, out}-h_{\mathrm{CW}, in}\right).$$

(94)

Mistry et al. [24] introduced a detailed and an updated steady-state model for a forward feed MED system. For the first time, a new scheme where the vapor exiting a stage is directed to preheat the coming seawater before being directed to the next stage was presented. More importantly, the $P_{\mathrm{R}}(-)$. and the specific area are related to the number of stages, the heating steam temperature, and the Recovery Ratio ($R_{\mathrm{R}}$) ($-$). The model of Mistry et al. [24] was developed based on fewer assumptions. These include:

• Zero salinity of produced water,
• The physical properties of feed seawater are only related to salinity and temperature,
• Seawater is an incompressible liquid,
• Neglected energy losses to the environment,
• Neglected the non-condensable gasses,
• Distillate vapor is marginally superheated,
• Average overall heat transfer coefficient,
• The temperature only related to heat transfer coefficient in each stage, feed heater, and condenser.

The model provides specific equations for each compartment of the MED system, including the whole system, 1st stage, 2nd stage through n stages, flash box, mixing box, feed heater, condenser, and distillate box. Moreover, the performance parameters of $P_{\mathrm{R}}(-)$, $R_{\mathrm{R}}$ (−), and total specific heat transfer area ($s_{\mathrm{A}}$) (m² s/kg) are calculated based on the following equations:

$$P_{\mathrm{R}}=\frac{q_{\mathrm{m},\mathrm{d}}}{\dot{q}m{,}_{\mathrm{s}}},$$

(95)

$$R_{\mathrm{R}}=\frac{\dot{q}m{,}_{\mathrm{d}}}{q_{\mathrm{m},\mathrm{f}}},$$

(96)

$$s_{\mathrm{A}}=\frac{{\sum }^{\text{}}{A}_{\mathrm{e}}+\sum A_{\mathrm{fh}}+A_{\mathrm{c}}}{q_{\mathrm{m},\mathrm{d}}},$$

(97)

where $A_{\mathrm{e}}$, $A_{\mathrm{fh}},$ and $A_{\mathrm{c}}$ (m²) are the heat transfer area in stage, feed heater, and condenser, respectively.

Druetta et al. [7] developed a non-linear predictive model for forward feed MED process due to mass and energy balances and ignoring thermodynamic losses. In this regard, specific correlations to calculate heat transfer areas in preheaters, evaporation stages and condenser are presented. The model was used to analyze and optimize the system to predict the optimum configuration, operating conditions, and size of each stage. This, in turn, elaborated several flow patterns for the vapor and distillate, which results in a decline of the total specific heat transfer area by 5% in comparison to the conventional configuration. The model developed used several assumptions as follows

• Fixed boiling point elevation,
• neglected the non-condensable gasses,
• neglected the pressure-drop in the demister and throughout the vapor condensation,
• fixed specific heat capacity,
• ignored the influence of fouling and non-condensable gasses on the heat transfer coefficients in the evaporators, preheaters, and the condenser.

More specifically, Druetta et al. [7] assumed that each stage has its unique temperature-drop, distillate production and heat transfer area. To elucidate the system efficiency, the $C_{\mathrm{R}}$ (−), specific cooling water rate ($S_{\mathrm{cw}}$) (mass flowrate of feed/mass flowrate of distillate) ($\frac{\frac{kg}{s} F }{\frac{kg}{s} D}$) and specific heat transfer area ($s_{\mathrm{A}}$) (m² s/(kg of D) are used, besides the $P_{\mathrm{R}}$ (−), as follows:

$$C_{\mathrm{R}}=\frac{q_{\mathrm{m},\mathrm{d}}}{q_{\mathrm{m},\mathrm{f}}},$$

(98)

$$S_{\mathrm{cw}}=\frac{q_{\mathrm{m},\mathrm{Wc}}}{q_{\mathrm{m},\mathrm{d}}}.$$

(99)

Tahir et al. [25] improved a mathematical model for MED system based on the well-known assumptions, including the following:

• 1 °C is considered as the temperature losses in piping and throughout condensation,
• fixed total areas of the stages and preheaters, and
• time-dependent fouling is elaborated to predict the overall heat transfer coefficient.

The overall heat transfer coefficient ($h)$(kW/(m² K)) for each evaporator and preheater are predicted by Equations (100) and (101), respectively. These equations considered fouling factor outside the tubes of the evaporators ($R_{\mathrm{fo}})$ (kW/(m² K)) and inside the tubes of the preheaters and condenser ($R_{\mathrm{f}})$ (kW/(m² K)):

$$\frac{1}{h}=\frac{r_0}{r_{\mathrm{i}} h_{\mathrm{in}}}+R_{\mathrm{fo}}+r_{\mathrm{o}} \frac{\ln \left(\frac{r_{\mathrm{o}}}{r_{\mathrm{i}}}\right)}{{\kappa }_{\mathrm{w}}}+\frac{1}{h_{\mathrm{o}}},$$

(100)

$$\frac{1}{h}=\frac{r_0}{r_{\mathrm{i}} h_{\mathrm{in}}}+R_{\mathrm{f}} \frac{r_{\mathrm{o}}}{r_{\mathrm{i}}}+r_{\mathrm{o}} \frac{\ln \left(\frac{r_{\mathrm{o}}}{r_{\mathrm{i}}}\right)}{{\kappa }_{\mathrm{w}}}+\frac{1}{h_{\mathrm{o}}} ,$$

(101)

$$R_{\mathrm{f}}=0.199 \left(1-\exp \left(\frac{-time}{61.55}\right)\right),$$

(102)

$$h_{\mathrm{o}}=0.0004 \left({\rho }_{\mathrm{sw}}^{2}g \frac{{\kappa }_{\mathrm{sw}}^3}{{\mu }_{\mathrm{sw}}^2}\right)R{e}^{0.2}P{r}^{0.65}{q}^{n0.4},$$

(103)

$$h_{\mathrm{i}}=\left(1+\frac{3.8}{Z^{0.95}}\right)h_{\mathrm{u}},$$

(104)

$$h_{\mathrm{u}}=h_{\mathrm{l}} {\left(1-x\right)}^{0.8},$$

(105)

$$Z={\left(\frac{1}{x}-1\right)}^{0.8}P{r}^{0.4},$$

(106)

$$h_{\mathrm{l}}=0.023 R{e}^{0.8}P{r}^{0.4}\frac{{\kappa }_{\mathrm{LC}}}{2r_{\mathrm{i}}},$$

(107)

where $r_{\mathrm{o}}$ and $r_{\mathrm{i}}$ (m) are the outer and inner radius of the tubes, respectively, ${\kappa }_{\mathrm{w}}$ (W/(m K)) and ${\kappa }_{\mathrm{sw}}$ (W/(m K)) are the thermal conductivity of the tubes and feed seawater, respectively, and $h_{in}$ (kW/(m² K)) and $h_{\mathrm{o}}$ (kW/(m² K)) are the heat transfer coefficients of inner and outer tubes, respectively. $Pr$ (−), $Re(-),x$ (−) are Prandtl number outside the tubes, Reynolds number outside the tubes, and vapor phase mass fraction, respectively. ${\kappa }_{\mathrm{LC}} \left(\frac{\mathrm{W}}{\mathrm{m}}\mathrm{K}\right), {\rho }_{ }\left(\mathrm{kg}/{\mathrm{m}}^3\right), {\mu }_{\mathrm{sw}}\left(\frac{\mathrm{kg}}{\mathrm{m}}\mathrm{s}\right)$ are thermal conductivity in the liquid phase, seawater density, and viscosity, respectively. $g$ is acceleration due to gravity.

A detailed model for MED system was developed by Filippini et al. [26], which included the calculation of energy consumption. Interestingly, thorough thermodynamic equations are deployed to assess all the relevant thermodynamic properties of the process in terms of water salinity, fouling, and temperature. This, in turn, has relaxed the assumptions of fixed thermodynamic properties provided by Darwish et al. [18].

The division of total areas $A_{\mathrm{tot}}$ (m²) of $N$ evaporators, N− feed pre-heaters, and final condenser and product flow rate has generated a specific parameter that does not depend on plant capacity.

$$A{\mathrm{tot}}_s=\frac{A_{\mathrm{tot}}}{q_{\mathrm{m},\mathrm{d}}}.$$

(108)

Another specific parameter was derived by the division of seawater intake ($q_{\mathrm{m},\mathrm{Wc}}$) (kg/s) in the end condenser and the produced fresh water. This parameter also does not depend on the plant capacity as given in Equation (99).

Then, the total energy consumption $e_{\mathrm{s}}$ (kJ/kg) to produce 1 kg of freshwater is given in Equation (109)

$$e_{\mathrm{s}}=\frac{q_{m,s}{\lambda }_{\left(\mathrm{Ts}\right)}}{q_{\mathrm{m},\mathrm{d}}},$$

(109)

Other interesting variables evaluated by the model are:

• temperature direction in the plant: brine temperature, distillate temperature, and feed temperature;
• salinity direction in the plant;
• the thermal load use, i.e., the portion of power provided to the 1st stage is used for evaporation with respect to the portion used for feed heating; and
• the RR, that is, the distillate resulted from 1 kg of seawater.

The relationships between $P_{\mathrm{R}}$ ($-$) and Gain Output Ratio $G_{\mathrm{OR}}$ $(-)$ is stated as

$$P_{\mathrm{R}}=G_{\mathrm{OR}}\frac{2330 \mathrm{kJ}/\mathrm{kg}}{{\lambda }_{\left(\mathrm{Ts}\right)}}.$$

(110)

More specifically, $G_{\mathrm{OR}}$ $(-)$ is well-defined as the mass of distillate to the mass of input steam. In addition, the $P_{\mathrm{R}}$ $(-)$ of TVC $(G_{{\mathrm{OR}}_{\mathrm{TVC}}}$) (−) is related to its $G_{\mathrm{OR}}$ by

$$P_{{\mathrm{R}}_{\mathrm{TVC}}}=G_{{\mathrm{OR}}_{\mathrm{TVC}}}\frac{2330 \mathrm{KJ}/\mathrm{Kg}}{{\lambda }_{\left(\mathrm{Ts}\right)}}.$$

(111)

Khalid at el. [27] presented a mathematical model for a forward feed MED system based on mass and energy balance correlations and making the same assumptions as previous workers to determine the best location the TVC.

The heat transfer surface areas of the 1st stage ($A_1$) (m²) and other proceeding stages $A_{\mathrm{i}}$ (m²) from i = 2 to N are calculated by

$$A_1{h}_1\left(T_{\mathrm{s}}-T_1\right)=q_{,md1}\lambda {\mathrm{v}}_1+q_{\mathrm{m},\mathrm{f}}{c}_{\mathrm{p}}\left(T_1-T_{\mathrm{f}2}\right),$$

(112)

$$A_{\mathrm{i}}{h}_1\left(\Delta T_{\mathrm{i}}-\Delta T_{\mathrm{loss}}\right)=q_{\mathrm{m},\mathrm{di} }\lambda {\mathrm{v}}_i.$$

(113)

The supplied steam flow rate $q_{\mathrm{m},\mathrm{s}}$ (kg/s) of the 1st stage is

$$q_{\mathrm{m},\mathrm{s}}{\mathsf{\lambda }}_{\mathrm{s}}=q_{\mathrm{m},\mathrm{d}1} \lambda {\mathrm{v}}_1+q_{\mathrm{m},\mathrm{f}}{c}_p\left(T_1-T_{\mathrm{f}2}\right).$$

(114)

Equations (116) and (117) are used to signify the energy balance of the down condenser and an estimate of the heat load (${\Phi }_{\mathrm{c}}$) (kW) resulting from the condensation of vapor leaving the last stage:

$$q_{\mathrm{m},\mathrm{dN}}{\lambda }_{\mathrm{v}N}=(q_{\mathrm{m},\mathrm{f}}+q_{\mathrm{m},\mathrm{cw}})c_p\left(T_{\mathrm{f}}-T_{\mathrm{cw}}\right),$$

(115)

$$Q_{\mathrm{c}}=q_{\mathrm{m},\mathrm{dN}}{\lambda }_{\mathrm{v}N}.$$

(116)

3.2. Dynamic Modeling of MED System

The dynamic model is a vital tool that enables one to understand and analyze the process performance in the transient periods, such as start-up and prolonged operation time. It is important to know that the knowledge of transient and corresponding steady-state schemes based on a dynamic model is essential for any industrial process to explore the influence of any expected disturbance of the input parameters on the process behavior via simulation. More importantly, dynamic modeling of any process enables one to forecast how long it takes to regain the steady-state if an unexpected disturbance of an input parameter occurs during steady-state operation. Additionally, an optimal control strategy necessitates the availability of a generic dynamic model. In other words, the improvement of an industrial process, including MED, is dependent on a dynamic simulation model that can forecast the process indicators for a long time of operation. Although several steady-state models were developed for MED system, it can be said that there are so far few attempts in the literature regarding the dynamic modeling. The following exhibits the progress of dynamic model for MED system with a thorough explanation of the most important parts of the models developed.

The first dynamic model of MED system was developed by Aly and Marwan in Reference [28], based on mass, energy, and salt balance correlations for each stage, to study the transient behavior, including start-up, shutdown, troubleshooting (the cause of disorders), and load changes. The overall heat transfer coefficient was presumed to vary linearly with temperature. Moreover, a plausible set of presumptions were outlined, such as:

• non-condensable gases and
• no partial vapor condensation reflected in the tube bundle.

The dynamic modeling equation are implicitly presented within the next research of Mazini et al. [29] (2014).

A computer package was originated by Dardour et al. [30] to simulate the dynamic behavior of forward feed MED system coupled to a high temperature nuclear reactor and to provide a control strategy to optimize the operating conditions. The model has mass and energy balance correlations and complementary heat transfer correlations and physical properties. The dynamic model equations are:

The total liquid and steam mass balance in each stage are represented as

$$\frac{\partial q_{\mathrm{m},1\left(\mathrm{i}\right)}}{\partial t}=q_{\mathrm{m},\mathrm{f}\left(\mathrm{i}\right)}+q_{\mathrm{m},\mathrm{b}\left(\mathrm{i}\right)}-q_{\mathrm{m},\mathrm{s}\left(\mathrm{i}\right)},$$

(117)

$$\frac{\partial q_{\mathrm{m},\mathrm{v}\left(\mathrm{i}\right)}}{\partial t}=q_{\mathrm{m},\mathrm{s}\left(\mathrm{i}\right)}-q_{\mathrm{m},\mathrm{s}\left(\mathrm{i}+1\right)}.$$

(118)

The salt mass balance and energy balance of each evaporator are

$$\frac{\partial x_{1\left(\mathrm{i}\right)} q_{m,1\left(\mathrm{i}\right)}}{\partial t}=X_{\mathrm{f}\left(\mathrm{i}\right)}{q}_{\mathrm{m},\mathrm{f}\left(\mathrm{i}\right)}+X_{\mathrm{b}\left(\mathrm{i}\right)}{q}_{\mathrm{m},\mathrm{b}\left(\mathrm{i}\right)}-X_{1\left(\mathrm{i}\right)}{q}_{\mathrm{m},\mathrm{b}\left(\mathrm{i}+1\right)},$$

(119)

$$\frac{\partial e_{1\left(\mathrm{i}\right)}}{\partial t}=q_{\mathrm{m},\mathrm{f}\left(\mathrm{i}\right)}{h}_{\mathrm{f}\left(\mathrm{i}\right)}+q_{\mathrm{m},\mathrm{b}\left(\mathrm{I}\right)}{h}_{\mathrm{b}\left(\mathrm{i}\right)}+Q_{\left(\mathrm{i}\right)}-q_{\mathrm{m},\mathrm{s}\left(\mathrm{i}\right)}{h}_{\mathrm{s}\left(\mathrm{i}\right)}-q_{\mathrm{m},\mathrm{b}\left(\mathrm{i}+1\right)}{h}_{1\left(\mathrm{i}\right)}.$$

(120)

$X$ represents the salt concentration of the specified stage.

A hybrid forward feed MED plant driven by a low temperature static solar collector field was studied by Roca et al. [31]. To systematically control and optimize the distillate production of this hybrid system, a dynamic model was formed. The model was built based on mass and energy balance correlations considering a number of basic presumptions. These include, no flash vapor generated, no presence of non-condensable gases, fixed physical properties, fixed temperature-drops between stages, and no heat loss. However, the previous assumption of fixed heat transfer coefficient was relaxed in this model. The model is basically developed to forecast the thermal dynamics of the heater and the freshwater production rate. The model equations are as follows:

The mass balance of any stage is

$$\frac{d{q}_{\mathrm{m},\mathrm{b}\left(\mathrm{i}\right)}}{dt}=q_{\mathrm{m},\mathrm{b}\left(\mathrm{i}-1\right)}-q_{\mathrm{m},\mathrm{bi}}-q_{\mathrm{m},\mathrm{v}d},$$

(121)

$$\frac{d\left(M_{\mathrm{b}\left(\mathrm{i}\right)}{X}_{b\left(i\right)}\right)}{dt}=q_{\mathrm{m},\mathrm{b}\left(\mathrm{i}-1\right)}{X}_{\mathrm{b}\left(\mathrm{i}-1\right)}-q_{\mathrm{m},\mathrm{b}\left(\mathrm{i}\right)}{X}_{\mathrm{b}\left(\mathrm{i}\right)}.$$

(122)

The energy balance equations are

$$\frac{d(q_{\mathrm{m},\mathrm{b}\left(\mathrm{i}\right)} h_{\mathrm{b}\left(\mathrm{i}\right)}}{dt}=q_{\mathrm{m},\mathrm{b}\left(\mathrm{i}-1\right)}{h}_{\mathrm{b}\left(\mathrm{i}-1\right)}-q_{m,sw\left(i\right)}{h}_{sw\left(i\right)}-q_{\mathrm{m},\mathrm{v}d}{h}_{\mathrm{v}d}+\left(k_{\mathrm{c}}{A}_{\mathrm{c}}\Delta T\right).$$

(123)

$q_{\mathrm{m},\mathrm{v}d}$ (kg/s) and $q_{\mathrm{m},\mathrm{sw}}$ (kg/s) are the mass flow rates of vapor distillate and feed seawater, respectively.

Mazini et al. [29] created a mathematical dynamic model for MED-TVC system in terms of material, salt, and energy balance of three main compartment, i.e., evaporators, condenser, and TVC. The dynamic equations of MED system were derived for three lumps of each stage, including brine lump, vapor lump, and tube bundle, as described below. The model improved has used the steady-state correlations of El-Dessouky et al. to signify the TVC section. Moreover, the model neglects non-condensable gases and assumes no pressure difference for each stage. The dynamic material and energy balance for the brine lump are

$$\frac{d}{dt}{q}_{\mathrm{m},\mathrm{b}.i}=q_{\mathrm{m},\mathrm{f}}+q_{\mathrm{m},\mathrm{b},i-1}-q_{\mathrm{m},\mathrm{b},i}-q_{\mathrm{m},\mathrm{v},i},$$

(124)

$$\frac{d}{dt}\left(q_{\mathrm{m},\mathrm{b},i} h_{\mathrm{b},i}\right)=q_{\mathrm{m},\mathrm{f}}{h}_f+q_{\mathrm{m},\mathrm{b},i-1}{h}_{\mathrm{b},i-1}-q_{\mathrm{m},\mathrm{b},i}{h}_{\mathrm{b},i}-q_{\mathrm{m},\mathrm{v},i}{h}_{\mathrm{v},i}+Q_{\mathrm{E},i},$$

(125)

where $Q_{\mathrm{E},i}$ (kW) is the heat transfer rate of corresponding evaporator. The salt balance for the ith stage is

$$\frac{d}{dt}\left(q_{\mathrm{m},\mathrm{b},i}{X}_{\mathrm{b},i}\right)=q_{\mathrm{m},\mathrm{f}}{X}_{\mathrm{f}}+q_{\mathrm{m},\mathrm{b},i-1}{X}_{\mathrm{b}.i-1}-q_{\mathrm{m},\mathrm{b},i}{X}_{\mathrm{b},i}.$$

(126)

The dynamic material and energy balances for the vapor lump are

$$\frac{d}{dt}{q}_{\mathrm{m},\mathrm{v},i}=q_{\mathrm{m},\mathrm{e},i}-q_{\mathrm{m},\mathrm{v},i},$$

(127)

$$\frac{d}{dt}\left(M_{\mathrm{v},i}{h}_{\mathrm{v},i}\right)=M_{\mathrm{e},i}{h}_{\mathrm{v},i}-M_{\mathrm{v},i}{h}_{\mathrm{v},i},$$

(128)

where $q_{\mathrm{m},\mathrm{e},i}$(kg/s) is the total rate of vaporization generated in the stage.

To predict the temperature of distillate from each stage, Equation (130) can be used

$$\frac{d}{dt}\left(q_{\mathrm{m},\mathrm{T}}{h}_{\mathrm{T}}\right)=q_{\mathrm{m},\mathrm{v},i-1}\left(h_{\mathrm{v},i-1}-h_{\mathrm{d},i}\right)-Q_{\mathrm{E},i},$$

(129)

where $q_{\mathrm{m},\mathrm{T}}$ (kg), $h_{\mathrm{T}}$ (kJ/kg), and $h_{\mathrm{d},i}$ (kJ/kg) are the vapor mass in the tube bundles, enthalpy of vapor in the tube bundles and enthalpy of distillate exiting the ith stage, respectively.

The mathematical dynamic energy equation of the condenser tubes is

$$\frac{d}{dt}\left(q_{\mathrm{m},\mathrm{con}}{h}_{\mathrm{con}}\right)=Q_{\mathrm{con}}-q_{\mathrm{m},\mathrm{f}}{h}_{\mathrm{f}},$$

(130)

where $q_{\mathrm{m},\mathrm{con}}$ (kg) and $h_{\mathrm{con}}$ (kJ/kg) are the seawater mass in the condenser tubes and enthalpy of feed seawater in the condenser tubes, respectively. $Q_{\mathrm{con}}$ (kW) is the heat transfer rate of condenser.

The dynamic model was then deployed to simulate the dynamic behavior of brine levels, temperature, and salinity after employing a 5% step change in seawater temperature and flow rate. This indicated a response time of the brine level and salinity between 5 to 10 min to rectify the steady-state condition after the applied disturbance. However, a faster response between 1 to 2 min was registered for the temperature.

Azimibavil and Dehkordi [32] developed a dynamic model for forward feed MED system consisting of partial differential equations that covered outside and inside long horizontal tube-bundle based on simplified assumptions. These include:

• the impact of vapor shear stress on the liquid film is not reflected due to insignificant turbulence,
• ignored pressure-drop inside the tubes,
• regular liquid distribution of feed seawater is accomplished on the top row of tubes in the bundle, and
• regular distribution of falling film along the length of each row of tube-bundle.

Azimibavil and Dehkordi [32] studied the dynamic behavior of water and vapor streams in the tube bundle of a corresponding stage. The model uses partial differential equations with the operating variables varying in both x and y directions.

The measured variables are the quality of the vapor, the rate of vapor condensation inside the tubes, the falling film flow rate, water enthalpy in the tubes and the produced vapor rate. The presence of non-condensable gases was also acknowledged because of the insignificant heat transfer rate in the tube bundle with the prediction of the vapor condensation rate.

The mass, momentum, and energy correlations of creeping film on vertical wall in the evaporating phase are summarized below:

$$\frac{\partial {\rho }_{\mathrm{f}}\delta }{\partial t}+\frac{\partial {\rho }_{\mathrm{f}}{v}_{\mathrm{av}}\delta }{\partial y}=-\frac{q_{\mathrm{m},\mathrm{v}}\left(x,N,t\right)/2}{\left(r \theta \mathrm{w}\right)L},$$

(131)

$$\frac{\partial {\rho }_{\mathrm{f}}v}{\partial t}+v\frac{\partial {\rho }_{\mathrm{f} }v}{\partial y}={\mu }_{\mathrm{f}}\frac{{\partial }^{2}v}{\partial x^2}+\Delta \rho g,$$

(132)

$$\frac{\partial {\rho }_{\mathrm{f}}h\delta }{\partial t}+\frac{\partial {\rho }_{\mathrm{f}}h\overline{v}\delta }{\partial y}=\frac{Q\left(x,N,t\right)/2}{\left(r \theta \mathrm{w}\right)L}-\frac{h_{\mathrm{g} }{q}_{\mathrm{m},\mathrm{v}}\left(x,N,t\right)/2}{\left(r \theta \mathrm{w}\right)L},$$

(133)

where $v \left(\mathrm{m}/\mathrm{s}\right), v_{\mathrm{av}} \left(\mathrm{m}/\mathrm{s}\right), \delta \left(\mathrm{m}\right), x\left(\mathrm{m}\right), y\left(\mathrm{m}\right), N,t(-),r\left(\mathrm{m}\right),{\theta }_{\mathrm{w}}(°),L$ (m) are the falling film velocity, average falling film velocity, falling film thickness, streamwise spatial coordinate inside tube, tangential coordinate of falling film flow on tube wall, number of tubes, tube radius, circumferential angle of tube, and tube-bundle length, respectively. ${\mu }_{\mathrm{f}}$ (kg/ms) and $g$ (m/s²) are the dynamic viscosity and gravitational acceleration constant, respectively.

The obtained results of model led to determining the fluid flow distribution and the influence of transient conditions on fouling. This showed that high fouling risk can be present in the transient period with signifying how long it would exist.

A dynamic model for 12 stages forward feed MED-TVC system was presented by Cipollina et al. [33]. The model was applied to thoroughly define the dynamic response of the main variables after applying a specific disturbance. The model constitutes algebraic correlations for each subsection and ordinary differential mass and enthalpy balance equations for single evaporator stage. Moreover, a control equation was deployed to signify the last stage brine pool level. More importantly, the model of Cipollina et al. [33] relaxed the assumption of no presence of non-condensable gases, which negatively impacts the steady-state and dynamic behavior. Basically, the model assumed the following:

• the produced vapor is free of salt,
• adiabatic system of no heat loses,
• ideal gas law is applied,
• Non-condensable gases cannot enter the stages only through inlet seawater,
• Saturated vapor is reached the 1st stage due to superheating the TVC.

Total mass balance equations of the evaporator are

$$\frac{d{q}_{\mathrm{m},\mathrm{vap}}}{d\mathrm{t}}=q_{\mathrm{m},\mathrm{evap}}-q_{m,{\mathrm{vap}}_{\mathrm{out}}}+q_{\mathrm{m},\mathrm{flash}}+q_{m,{\mathrm{N}}_{{\mathrm{CG}}_{\mathrm{shell}}}}+q_{m,{\mathrm{N}}_{{\mathrm{CG}}_{\mathrm{entrained}}}}+q_{m,{\mathrm{N}}_{\mathrm{CG}}\_\mathrm{pool}},$$

(134)

$$\frac{d{m}_{\mathrm{shell}}}{d\mathrm{t}}=q_{\mathrm{m},\mathrm{f}}-q_{\mathrm{m},\mathrm{evap}}-q_{\mathrm{m},\mathrm{br}\_\mathrm{shell}}-q_{m,{\mathrm{N}}_{\mathrm{CG}}\_\mathrm{shell}},$$

(135)

$$\frac{d{m}_{\mathrm{pool}}}{d\mathrm{t}}=q_{\mathrm{m},\mathrm{br}\_\mathrm{shell}}+q_{\mathrm{m},\mathrm{br}\_\mathrm{in}}-q_{\mathrm{m},\mathrm{br}\_\mathrm{out}}-q_{\mathrm{m},\mathrm{flash}}-q_{m,{\mathrm{N}}_{\mathrm{CG}}\_\mathrm{pool}},$$

(136)

where $q_{\mathrm{m},\mathrm{vap}} \left(\mathrm{kg}/\mathrm{s}\right),q_{\mathrm{m},\mathrm{evap}} \left(\mathrm{kg}/\mathrm{s}\right), q_{m,{\mathrm{vap}}_{\mathrm{out}} }\left(\mathrm{kg}/\mathrm{s}\right), q_{\mathrm{m},\mathrm{flash}} \left(\mathrm{kg}/\mathrm{s}\right), q_{m,{\mathrm{N}}_{{\mathrm{CG}}_{\mathrm{shell}}}}\left(\mathrm{kg}/\mathrm{s}\right), q_{m,{\mathrm{N}}_{{\mathrm{CG}}_{\mathrm{entrained}}}}\left(\mathrm{kg}/\mathrm{s}\right), q_{m,{\mathrm{N}}_{{\mathrm{CG}}_{pool}}}\left(\mathrm{kg}/\mathrm{s}\right)$ are mass flow rate of vapor, mass flow rate of vapor generated from tube bundle evaporation surface, vapor leaving from the stage, vapor produced by the flashed brine arrived from the preceding stage, non-condensable gases released from the discrete brine around the tube bundle, non-condensable gases arriving from the previous stage’s vapor, and non-condensable gases released from the brine arriving from the previous stage phase, respectively. $m_{\mathrm{shell}} \left(\mathrm{kg}\right)$, and $m_{\mathrm{pool}} \left(\mathrm{kg}\right)$ are dispersed mass brine around the tube bundle and pool brine, respectively. $q_{\mathrm{m},\mathrm{br}\_\mathrm{shell}} \left(\mathrm{kg}/\mathrm{s}\right), q_{m,{\mathrm{br}}_{\mathrm{in}}}\left(\mathrm{kg}/\mathrm{s}\right)$, and $q_{\mathrm{m},\mathrm{br}\_\mathrm{out}} \left(\mathrm{kg}/\mathrm{s}\right)$ are mass flow rates of dispersed liquid reaching the brine pool from the tube bundle, brine stream arriving from the former stage, and brine stream departing the stage to the following one, respectively.

Mass balance equations for the salt are

$$\frac{d{m}_{\mathrm{shell}}}{dt}=q_{\mathrm{m},\mathrm{f}}{X}_{\mathrm{f}}-q_{\mathrm{m},\mathrm{br}\_\mathrm{shell}}{X}_{\mathrm{shell}},$$

(137)

$$\frac{d{m}_{\mathrm{pool}}}{dt}=q_{\mathrm{m},\mathrm{br}\_\mathrm{shell}}{X}_{\mathrm{shell}}+q_{\mathrm{m},\mathrm{br}\_\mathrm{in}}{X}_{\mathrm{br}\_\mathrm{in}}-q_{\mathrm{m},\mathrm{br}\_\mathrm{out}}{X}_{\mathrm{pool}},$$

(138)

where $m_{\mathrm{shell}}$ ($\mathrm{kg}$) and $m_{\mathrm{pool}} \left(\mathrm{kg}\right)$ are the total mass of salt in the dispersed brine around the tube bundle and pool brine, respectively. $X$ ($\frac{\mathrm{kg}\mathrm{of}\mathrm{non}-\mathrm{condensable}\mathrm{gases}}{\mathrm{kg}\mathrm{of}\mathrm{Liquid}}$) denotes the mass fraction in parts per million of the corresponding stream. The mass balance equations of the non-condensed gases are

$$\frac{dNC{G}_{\mathrm{vap}}}{dt}=q_{m,N_{CG}\_\mathrm{shell}}+q_{m,N_{CG}\_\mathrm{entrained}}+q_{m,N_{CG}\_\mathrm{pool}}-q_{\mathrm{m},\mathrm{vap}\_\mathrm{out}}{Y}_{\mathrm{vap}},$$

(139)

$$\frac{dNC{G}_{\mathrm{shell}}}{dt}=q_{\mathrm{m},\mathrm{f}}{X}_{\mathrm{f}}-q_{m,N_{CG}\_\mathrm{shell}}-q_{\mathrm{m},\mathrm{br}\_\mathrm{shell}}{X}_{\mathrm{shell}},$$

(140)

$$\frac{dNC{G}_{\mathrm{pool}}}{dt}=q_{\mathrm{m},\mathrm{br}\_\mathrm{in}}{X}_{\mathrm{br}\_\mathrm{in}}+q_{\mathrm{m},\mathrm{br}\_\mathrm{shell}} X_{\mathrm{shell}}-q_{\mathrm{m},\mathrm{br}\_\mathrm{out}}{X}_{\mathrm{pool}}-q_{m,N_{CG}\_\mathrm{pool}}.$$

(141)

The mass flow rate of non-condensable gases in the vapor, dispersed brine around the tube bundle, and pool brine phases are denoted as $N_{{CG}_{\mathrm{vap}}}, \left(\mathrm{kg}\right), N_{{CG}_{\mathrm{shell}}}$ ($\mathrm{kg}$), and $N_{{CG}_{\mathrm{pool}}} \left(\mathrm{kg}\right)$, respectively.

Finally, the enthalpy balance correlations in the shell and pool are

$$\frac{d{H}_{\mathrm{shell}}}{dt}=q_{\mathrm{m},\mathrm{f}}{h}_{\mathrm{f}}+Q_{\mathrm{eff}}-q_{\mathrm{m},\mathrm{evap}}{h}_{\mathrm{evap}}-q_{\mathrm{m},\mathrm{br}\_\mathrm{shell}}{h}_{\mathrm{shell}}-q_{m,N_{CG}\_\mathrm{shell}}{h}_{N_{CG}\_\mathrm{shell}},$$

(142)

$$\frac{d{H}_{\mathrm{pool}}}{dt}=q_{\mathrm{m},\mathrm{br}\_\mathrm{shell}}{h}_{\mathrm{shell}}+q_{\mathrm{m},\mathrm{br}\_\mathrm{in}}{h}_{\mathrm{br}\_\mathrm{in}}-q_{m,{\mathrm{br}}_{\mathrm{out}}}{h}_{\mathrm{pool}}-q_{\mathrm{m},\mathrm{flash}}{h}_{\mathrm{flash}}-q_{m,{\mathrm{N}}_{CG}\_\mathrm{pool}}{h}_{{\mathrm{N}}_{CG}\_\mathrm{pool}},$$

(143)

where $H_{\mathrm{shell}} \left(\mathrm{kJ}\right), H_{\mathrm{pool}} \left(\mathrm{kJ}\right)$, and $Q_{\mathrm{eff}} \left(\mathrm{kW}\right)$ are the enthalpy of the dispersed brine around the tube bundle and pool brine and heat transferred through the tube bundle, respectively.

4. Improvement of Modeling of MED System and Its Subsequent Operation

The modeling of MED system is essential to guarantee the implementation of accurate simulation and parametric studies. This would offer a deeper understanding of the process variables leading to successful design resulting in optimum process efficiency at the lowest irreversibility. The following points are some suggestions to improve the MED modeling and process operation.

• The modeling of MED process has been associated with a permanent set of stringent assumptions, such as constant temperature-drop between the stages, fixed distillate production in all the stages, and fixed heat transfer coefficients. It is believed that such assumptions need to be relaxed in the upcoming models to mitigate its passive effect on the model prediction.
• More efforts are needed to understand the thermodynamic limitations of MED process and exergy analysis to narrow the divergence between the model predictions and actual energy level.
• Ignoring the heat losses to the environment throughout the MED process is not a feasible option. Therefore, the empowerment of well-covered and insulated evaporators would resolve the ignored heat losses and possibly allow the utilization of heating inside the evaporators and condenser, thus harvesting a higher freshwater production rate.
• The study of combination of parallel and forward feed configurations, mathematically or experimentally, has not been reported in the literature. It is fair that such integration would mitigate the demerits of each individual configuration and promote distillate production. An accurate model and consequence simulation are needed to realize the feasibility of such design.
• There are many successful theoretical and experimental attempts in the literature of integrated systems of MED and other processes that have resulted in significant energy recovery. However, substantial efforts are required for the current research to lift the thermodynamic efficacy and promote the heat transfer area and reduce the required steam by exploring the most feasible process to be integrated with MED process. In addition, the research needs to be directed towards decreasing corrosion and consequently to reduce the total production cost.

It is fair to expect that the above suggestions can be accordingly deployed in the current MED industry to exceed the operational performance of the conventional desalination plants.

5. Conclusions

This paper has reviewed and discussed the modeling of forward feed MED based on steady-state and dynamic models developed. Specifically, the models were applied to predict the process performance in seawater desalination. Here, the evolution of these models and highlighting the merits and demerits of each model were carried out. A critical appraisal of 29 mathematical models has been carried out, and suggestions for improvement have been proposed. The most significant conclusions made are:

• Two sets of models can be found in the open literature, simple and detailed ones. Simple models can moderately offer acceptable estimation of the main key performance indicators, such as performance ratio, specific heat transfer areas, and flow rate of cooling water. However, detailed models can provide a more precise estimation for a wide range of performance indicators and implicit variables.
• The dynamic modeling of MED system is scarce compared to steady-state modeling.
• It is well agreed that the utilization of several approximations in the modeling of MED process can result in large inconsistencies between the actual data and model predictions.
• The dynamic modeling of MED process is essential in capturing the transient periods, such as start-up and prolonged operation time, which is not possible in steady-state modeling.
• Most of the models developed have discounted the presence of non-condensable gases that marginally influence the estimation of heat transfer areas of the evaporators.
• For MED plants with large number of stages, the assumptions of fixed thermodynamic losses and specially boiling point elevation is not plausible due to increased discrepancy in evaluating the specific heat area.

This comprehensive review would serve as a platform when developing mathematical models in the future since it clarifies the main issues of process modeling to be resolved. This, in turn, would strengthen the prediction of MED process behavior in both steady-state and dynamic operations.