Dynamic Modelling and Simulation of a Multistage Flash Desalination System

: As the leading thermal desalination method, multistage ﬂash (MSF) desalination plays an important role in obtaining freshwater. Its dynamic modeling and dynamic performance prediction are quite important for the optimal control, real-time optimal operation, maintenance, and fault diagnosis of MSF plants. In this study, a detailed mathematical model of the MSF system, based on the ﬁrst principle and its treatment strategy, was established to obtain transient performance change quickly. Firstly, the whole MSF system was divided into four parts, which are brine heat exchanger, ﬂashing stage room, mixed and split modulate, and physical parameter modulate. Secondly, based on mass, energy, and momentum conservation laws, the dynamic correlation equations were for-mulated and then put together for a simultaneous solution. Next, with the established model, the performance of a brine-recirculation (BR)-MSF plant with 16-stage ﬂash chambers was simulated and compared for validation. Finally, with the validated model and the simultaneous solution method, dynamic simulation and analysis were carried out to respond to the dynamic change of feed seawater temperature, feed seawater concentration, recycle stream mass ﬂow rate, and steam temperature. The dynamic response curves of TBT (top brine temperature), BBT (bottom brine temperature), the temperature of ﬂashing brine at previous stages, and distillate mass ﬂow rate at previous stages were obtained, which speciﬁcally reﬂect the dynamic characteristics of the system. The presented dynamic model and its treatment can provide better analysis for the real-time optimal operation and control of the MSF system to achieve lower operational cost and more stable freshwater quality.


Introduction
With the development of society, the shortage of freshwater resources has gotten worse, and seawater desalination is an important way to solve it [1,2]. So far, among the many seawater desalination methods, multistage flash (MSF) desalination leads the way in the desalination industry because of its high reliability, large capacity of single units, and good water quality [3][4][5][6]. With the MSF desalination process becoming more sophisticated, it is necessary to establish a dynamic model of the MSF process [7][8][9]. Dynamic models can be used to solve problems related to the transient behaviors of the device, such as control strategies, process interaction and troubleshooting, reliability and stability analysis, and startup and shutdown conditions. Therefore, it is of great significance to establish a strict dynamic model and simulation analysis of the system.
Since the 1970s, the dynamic modeling and analysis of MSF have been gradually carried out on the basis of steady-state research. Previously, Glueck and Bradshaw [10] studied the dynamic model of the MSF process and presented a dynamic model, including the energy difference balance between the steam space and the flash chamber distillate, but they did not provide simulation results. Ulrich [11] applied the model to simulate an MSF experimental device containing 6-stage flash chambers and gave the simulation results of the system. However, the simulation results had significant drawbacks in the disturbance Processes 2021, 9,522 3 of 25 (T sea ), feed seawater salt concentration (C F ), recycle stream mass flow rate (W Re ), and steam temperature (T steam ) were simulated and analyzed. For the quick and easy implementation of the dynamic simulation, the dynamic equations are discretized with the Runge-Kutta algorithm and solved with the Interior Point solver in the MATLAB platform. These lay a foundation and guarantee a further understanding of the dynamic change of system state and real-time optimal operation.

Multistage Flash (MSF) Process Description
The MSF desalination system has two structures, which are the one-through multistage flash (OT-MSF) seawater desalination system and the brine-recirculation multistage flash (BR-MSF) seawater desalination system [23][24][25][26]. The brine-recirculation type of MSF process is said to be the industry standard [27], so the dynamic simulation and analysis in this article are based on the BR-MSF system; the flowsheet of the BR-MSF system is shown in Figure 1.
Processes 2021, 9,522 3 of 28 establishes a relatively complete and realistic dynamic mechanism model. Based on the verification and analysis of the model, the dynamic process of feed seawater temperature ( sea T ), feed seawater salt concentration ( F C ), recycle stream mass flow rate ( Re W ), and steam temperature ( steam T ) were simulated and analyzed. For the quick and easy implementation of the dynamic simulation, the dynamic equations are discretized with the Runge-Kutta algorithm and solved with the Interior Point solver in the MATLAB platform. These lay a foundation and guarantee a further understanding of the dynamic change of system state and real-time optimal operation.

Multistage Flash (MSF) Process Description
The MSF desalination system has two structures, which are the one-through multistage flash (OT-MSF) seawater desalination system and the brine-recirculation multistage flash (BR-MSF) seawater desalination system [23][24][25][26]. The brine-recirculation type of MSF process is said to be the industry standard [27], so the dynamic simulation and analysis in this article are based on the BR-MSF system; the flowsheet of the BR-MSF system is shown in Figure 1. The BR-MSF system is composed of a brine heater, a heat rejection section, a heat recovery section, splitters, and mixers. As shown in Figure 1, feed seawater is preheated through the condenser tube of the heat rejection section and divided into two parts; one is returned to the sea and the other is mixed with the recycled brine. The mixed brine enters the end of the heat recovery section and is heated step by step while condensing the flash steam in the flash chamber. Finally, the mixed brine flows out from the first stage of the heat recovery section, enters the brine heater for further heating, and flows into the first stage flash chamber at the highest temperature. As the pressure of the flash chamber decreases step by step, the flash stream enters the flash chamber and evaporates immediately. The generated steam passes through the condensers and drops into the distillate trays after being condensed. The process is repeated until the last stage is reached, the concentrated brine is rejected, the distillate is extracted, and a part of the brine is reused as recycled brine.

Dynamic Mathematical Modeling of the Multistage Flash (MSF) Process
The complete dynamic mathematical model of the BR-MSF desalination system includes four parts: the brine heater module, the flash chamber module, the brine heater module, the splitters and mixers module, and physical parameter equations. Next, the above four parts are modeled based on the following main assumptions: 1. The distillate from whichever stage is salt-free; 2. Noncondensable gases are ignored; The BR-MSF system is composed of a brine heater, a heat rejection section, a heat recovery section, splitters, and mixers. As shown in Figure 1, feed seawater is preheated through the condenser tube of the heat rejection section and divided into two parts; one is returned to the sea and the other is mixed with the recycled brine. The mixed brine enters the end of the heat recovery section and is heated step by step while condensing the flash steam in the flash chamber. Finally, the mixed brine flows out from the first stage of the heat recovery section, enters the brine heater for further heating, and flows into the first stage flash chamber at the highest temperature. As the pressure of the flash chamber decreases step by step, the flash stream enters the flash chamber and evaporates immediately. The generated steam passes through the condensers and drops into the distillate trays after being condensed. The process is repeated until the last stage is reached, the concentrated brine is rejected, the distillate is extracted, and a part of the brine is reused as recycled brine.

Dynamic Mathematical Modeling of the Multistage Flash (MSF) Process
The complete dynamic mathematical model of the BR-MSF desalination system includes four parts: the brine heater module, the flash chamber module, the brine heater module, the splitters and mixers module, and physical parameter equations. Next, the above four parts are modeled based on the following main assumptions: 1.
The distillate from whichever stage is salt-free; 2.
The system is adiabatic; 4.
The influence of the pump is not considered; 5.
The evaporation of freshwater is ignored; 6.
The amount of flashing seawater in the condenser tubes remains constant, and there is no accumulation of salt. The heating steam enters the brine heater at a certain temperature and flow rate, flows through the condenser tubes, and is finally pumped away by the condensing pump. During this period, the cooling brine in the condenser tubes is heated to the top brine temperature (TBT) and enters the first flash chamber for flash. Then, the above process is transformed into a mathematical model, in which the brine heater is generally regarded as the 0-th flash chamber.
Mass balance: Salt mass balance: Assuming that the amount of brine holdup in the condenser tubes remains unchanged and there is no salt accumulation in the condenser tubes, then: where B 0 M refers to the brine holdup in the brine heater, in kg; 0 F C is the salt concentration of flashing seawater leaving the brine heater, in wt %.
Heating steam enthalpy balance: The heating steam enters the brine heater at a certain temperature and flow rate, flows through the condenser tubes, and is finally pumped away by the condensing pump. During this period, the cooling brine in the condenser tubes is heated to the top brine temperature (TBT) and enters the first flash chamber for flash. Then, the above process is transformed into a mathematical model, in which the brine heater is generally regarded as the 0-th flash chamber.
Mass balance: Salt mass balance: where W B0 is the flashing brine mass flow rate leaving the brine heater, in kg·h −1 ; W F1 is the flashing seawater mass flow rate leaving Stage 1, in kg·h −1 ; C B0 is the salt concentration in the flashing seawater leaving the brine heater, in wt %; C F1 is the salt concentration of flashing seawater leaving Stage 1, in wt %.
Assuming that the amount of brine holdup in the condenser tubes remains unchanged and there is no salt accumulation in the condenser tubes, then: where M B0 refers to the brine holdup in the brine heater, in kg; C F0 is the salt concentration of flashing seawater leaving the brine heater, in wt %. Heating steam enthalpy balance: Overall enthalpy balance: where W steam is the steam mass flow rate, in kg·h −1 ; U H is the overall heat transfer coefficient at the brine heater, in kcal·(m 2 ·h· • C) −1 ; A H is the heat transfer area of the brine heater, in m 2 ; T B0 is the temperature of flashing brine leaving the brine heater, in • C; T F1 is the temperature of cooling brine leaving Stage j, in • C; λ steam is the latent heat of the steam, in kcal·kg −1 ; h B0 is the specific enthalpy of the brine heater, in kcal·kg −1 ; h F1 is the specific enthalpy of flashing seawater at Stage j, in kcal·kg −1 . λ steam , h B0 , and h F1 in the above model are functions of concentration and temperature. For specific formulas, please refer to the physical parameter equations.

Model of the Flash Chamber Module
As shown in Figure 3, the j-th stage flash chamber can be divided into four parts: the brine pool, the vapor space, the product tray, and condenser tubes. At the same time, each stage of the flash chamber includes three steps: stream flash, steam condensation, and seawater preheating. Then, based on the four major components and three basic links, dynamic modeling of the flash chamber is carried out. In the process of modeling, thermodynamic loss, the variation of physical parameters with temperature and concentration, the difference of nonideality, and the basic physicochemical phenomenon of an unequal temperature difference between different stages were considered. Considering that there is very little evaporation of freshwater in the distillate tray, the evaporation of freshwater is ignored in this model. The flash chamber module equations in this paper refer to the research of Mazzotti et al. [16].  Flashing brine mass balance: where Flashing brine salt balance: Flashing brine mass balance: where M B j = ρ B j * V represents the flashing brine holdup in the j-th stage flash chamber, in kg; V is the volume of flashing brine holdup, in m 3 ; W Bj is the flashing brine mass flow rate leaving Stage j, in kg·h −1 ; V Bj is the evaporation capacity of brine at Stage j, in kg·h −1 . Flashing brine salt balance: where C Bj is the salt concentration in the flashing brine leaving stage j, in wt %. Flashing brine energy balance: where CP Bj is the heat capacity of brine leaving Stage j, in kcal·(kg· • C) −1 ; T Bj is the temperature of flashing brine leaving Stage j, in • C; h Bj is the specific enthalpy of flashing brine at Stage j, kcal·kg −1 . Distillate mass balance: where M D j = ρ D j * V represents the distillate holdup in the j-th stage flash chamber, in kg; V is the volume of distillate holdup, in m 3 ; W Dj is the distillate mass flow rate leaving Stage j, in kg·h −1 ; F Dj is the steam condensation at Stage j, in kg·h −1 . Distillate energy balance: where h Dj is the specific enthalpy of distillate at Stage j, kcal·kg −1 ; h FDj is the specific enthalpy of condensation at Stage j, in kcal·kg −1 . Temperature relationship between distillate and flashing brine: where BPE represents boiling point elevation, NETD represents the nonequilibrium allowance, and TL represents the temperature loss due to demister and condenser; the specific calculation formula is shown in the physical parameter equations; T Dj is the temperature of distillate leaving Stage j, in • C. Steam mass balance: where M V j = ρ V j * V represents the steam holdup in the j-th stage flash chamber, in kg, and V is the volume of steam holdup, in m 3 . Steam energy balance: where CP vj is the heat capacity of the flashing steam leaving Stage j, in kcal·(kg· • C) −1 ; T vj is the temperature of flashed vapor at Stage j, in • C; h VBj is the specific enthalpy of steam at Stage j, kcal·kg −1 ; h FDj is the specific enthalpy of condensation at Stage j, kcal·kg −1 ; U j is the overall heat transfer coefficient at Stage j, in kcal·(m2 ·h· • C) −1 ; T Fj is the temperature of cooling brine leaving Stage j, in • C. Temperature relationship between distillate and steam: where ∆TL j is the temperature difference loss, in • C. Assuming that the amount of flashing seawater in the condenser tubes remains constant and there is no accumulation of salt, the conservation equations of flashing seawater mass and salt in the condenser tubes are as follows: where M F j represents the flashing seawater holdup in the j-th stage flash chamber, in kg; C Fj is the salt concentration of flashing seawater leaving Stage j, in wt %. Because of the constraint of the condenser tubes, its value remains constant during the flashing process, which can be obtained from the product of the volume of the condenser tubes and the density of the flashing brine. Condenser tubes energy balance: Volume constraint of the flash chamber: where ρ V j is the density of steam at Stage j, kg·m −1 ; ρ Bj is the density of flashing brine at Stage j, kg·m −1 ; ρ Dj is the density of distillate at Stage j, kg·m −1 ; V TUBE is the volume of the condenser tube, m 3 ; N TUBE is the number of condenser tubes.

Model of the Splitters and Mixers Module
Splitters: Equations (20)-(22) express the flow relationship of the three splitters, respectively. W BD is the blowdown mass flow rate, in kg·h −1 . W BN is the flashing brine mass flow rate leaving Stage N, in kg·h −1 . W m is the make-up brine mass flow rate, in kg·h −1 . W F is the flashing seawater mass flow rate to the rejection section, in kg·h −1 . S is the reject recycled mass flow rate, in kg·h −1 . W r is the mass flow rate to the reject seawater splitter, in kg·h −1 . C W is the rejected seawater mass flow rate, in kg·h −1 . Mixers: Equations (23) and (26) express the flow relationship of the two mixers, respectively. Equation (24) expresses the slat balance at one of the mixers. Equations (25) and (27) express the energy balance of the mixers. W R is the cooling brine mass flow rate to the recovery section, in kg·h −1 . C R is the salt concentration in the cooling brine to the recovery section, in wt %. C Re is the recycled brine concentration, in wt %. C m is the salt concentration in make-up water, in wt %. h R is the specific enthalpy of stream to the recovery section, in kcal·kg −1 . h Re is the specific enthalpy of the recycled brine, in kcal·kg −1 . h m is the specific enthalpy of the make-up brine, in kcal·kg −1 . W F is the flashing seawater mass flow rate to the rejection section, kg·h −1 . W S is the seawater mass flow rate, in kg·h −1 . h W F is the specific enthalpy of the brine at the entrance of the rejection section, in kcal·kg −1 . h S is the specific enthalpy of recycled brine at the rejection stage, in kcal·kg −1 . h W S is the specific enthalpy of feed seawater, in kcal·kg −1 .

Physical Parameter Equations
The MSF desalination system exhibits a strong nonlinearity. In order to achieve better simulation performance, the physical and chemical properties of the model variables in the system must be well characterized. The physical parameter equations in this paper refer to the research of Woldai et al. [28] and Helal et al. [29].
Specific heat capacity calculation equation: Equation (28) is integrated between the reference temperature of 32 • F and the boiling temperature in • F to obtain the following equation: Enthalpy calculation equation: Calculation equation of flash temperature difference loss: where where ω j = WF/w j ; w j represents the width of the flash chamber; H j represents the flashing brine level in the j-th stage flash chamber; ∆T B represents the flashing brine temperature difference between the two stages.

Degree of Freedom Analysis
The degree of freedom analysis is an important step before solving the model. It is important to ensure that the model equations have a unique relationship between all inputs and outputs [30]: where N f is the degree of freedom, N v is the number of variables, including input and output, and N e is the number of independent differential and algebraic equations. When N f = 0, it is the best case; the system can correspond to the unique solution. When N f < 0, the number of variables is greater than the number of equations, and additional model equations need to be established. When N f > 0, the number of variables is less than the number of equations, and the system has not determined enough variables.
The flash chamber is divided into four parts for the degree of freedom analysis: there are 4 unknown variables in the flash chamber part, 3 unknown variables in the product tray part, and 3 unknown variables in the condenser tube part. There are 5 unknown variables in the brine heater. In summary, there are 13 unknown variables in each stage of the flash chamber, 5 unknown variables in the brine heater, and 9 unknown variables in the splitters and mixers module. If the BR-MSF device has N stages, the system has (13N + 14) unknown variables.
Regardless of physical parameters such as enthalpy value, specific heat capacity, and heat loss, N v = N e = (13N + 14) in the model. From Equation (39), N f = 0. Therefore, the BR-MSF dynamic model established in this paper is complete and has a unique solution.
It is impractical to directly solve the model with so many dynamic and physical equations. Hence, the dynamic equations were firstly described into algebraic equations by the Runge-Kutta algorithm; then, all the equations, including discrete dynamic equations and physical equations, were put together by simultaneous strategy, and the discretized model, in the form of nonlinear programming, was solved by the Interior Point solver through sparsity treatment in the MATLAB platform in order to carry out the simulation quickly. To fully respond to the dynamic process with different disturbances, the entire simulation time is set as 7200 s, and the interval is set as 20 s. Therefore, the model was discretized into 360 parts in the whole-time domain. Under conventional conditions, one simulation work can be finished within 50 s under the MATLAB platform to solve the discretized model with the Interior Point solver. The platform and model solution strategy lay a foundation and guarantee for real-time optimal operation in future work.

Model Validation
As the research of Rosso et al.
[31] is based on real industrial data, this paper uses the system parameters in the paper by Rosso et al  Tables 1 and 2; the initial operation parameters in the simulation process are shown in Table 3.   Feed Seawater temperature Based on the above parameters, the seawater desalination process is simulated. Figures 4-8 show the process of the TBT, bottom brine temperature (BBT), and flashing brine temperature of the first three stages (T B1 , T B2 , T B3 ), distillate flow rate of the first three stages (W D1 , W D2 , W D3 ), and total distillate flow rate (W DN ) from a given initial value to steady state during the dynamic simulation process. According to the simulation results, the system reaches a steady state at around 3000 s.                 Since the dynamic response of disturbance will drive the multistage flash process from one steady state to another steady state, at that time, the final value of dynamic response should be consistent with the new balanced steady state. Hence, the validation is carried out as follows. Assuming that the level of all stages of flashing brine is 0.457 m, the initial values are calculated based on the results of the steady-state model. Based on the dynamic model established above, the dynamic simulation of the system is carried out under the MATLAB platform.
According to Figures 4-8, the system reaches a steady state around 3000 s. The steady-state value of the system after 3000 s is extracted. Figure 9 shows both the flashing brine temperature after 3000 s during the dynamic process (T Bj ) and the flashing brine temperature within the steady-state simulation results (T 0 Bj ) of Gao Hanhan et al.
[32]. Figure 10 shows both the distillate temperature after 3000 s during the dynamic process (T Dj ) and the distillate temperature during the steady-state process (T 0 Dj ). Figure 11 shows the distillate flow rate after 3000 s during the dynamic process (W Dj ) and the distillate flow rate during the steady-state process at all stages (

Dynamic Simulation and Analysis
In order to study the performance changes of the multistage flash desalination system under different working conditions, based on the dynamic model established above, the influence of sea T , Re W , and steam T on the performance of the system is discussed in this paper.

Dynamic Simulation and Analysis
In order to study the performance changes of the multistage flash desalination system under different working conditions, based on the dynamic model established above, the influence of sea T , Re W , and steam T on the performance of the system is discussed in this paper.

Dynamic Simulation and Analysis
In order to study the performance changes of the multistage flash desalination system under different working conditions, based on the dynamic model established above, the influence of T sea , W Re , and T steam on the performance of the system is discussed in this paper.

Dynamic Response Analysis of Feed Seawater Temperature Change
Seawater temperature is a variable input affected by external disturbances, which will change with the seasons and different times of the day. It is meaningful to analyze the impact of T sea changes on the system. Under the condition that other operating parameters remain unchanged, T sea increases (from 35 to 39 • C) during the simulation process.
As shown in Figure 12, when the system runs normally to 4000 s, a disturbance is caused by the increase of T sea by 300 s; that is, T sea increases from 35 to 39 • C. When the system runs to 4300 s, the disturbance stops and T sea recovers to 35 • C and this lasts until the end of 7200 s.
The dynamic responses of the T sea increase are shown in Figures 13-16. The disturbance time of the T sea increase is 5 min, and the dynamic response time of TBT, BBT, T B1 , T B2 , T B3 , and W D1 , W D2 , W D3 is about 8 min. The results show that the changes of TBT, BBT, and T B1 , T B2 , T B3 are proportional to the change of T sea . BBT is more sensitive to the change in T sea than TBT. However, the change in W D1 , W D2 , W D3 is inversely proportional to the change in T sea , showing a downward trend at first and then an upward trend. The dynamic trends of TBT, BBT, and T B1 , T B2 , T B3 with increasing seawater temperature, are the same as in the paper by Jari et al. [22].

cess.
As shown in Figure 12, when the system runs normally to 4000 s, a disturbance is caused by the increase of sea T by 300 s; that is, sea T increases from 35 to 39 °C. When the system runs to 4300 s, the disturbance stops and sea T recovers to 35 °C and this lasts until the end of 7200 s.     The increase in sea T reduces the overall temperature drop of the device, indicating a decrease in the heat consumption of the flashing brine in the system. According to the law of conservation of energy, the heat lost by the flashing brine will be converted into the heat absorbed by the flashing steam, resulting in less distillate. Therefore, as shown in Figure 16, when the disturbance occurs,

Dynamic Response Analysis of Feed Seawater Salt Concentration Change
F C is one of the important parameters that are susceptible to external influences in the MSF-BR process, and its change will also affect the performance of the entire device.
Under the condition that other operating parameters remain unchanged, F C is increased (from 5.7 to 6.7 wt %) during the simulation process. As shown in Figure 17, when the system runs normally to 4000 s, a disturbance is caused by the increase of F C by 300 s; that is, F C increases from 5.7 to 6.7 wt %. When the system runs to 4300 s, the disturbance stops and F C recovers to 5.7 wt % and this lasts until the end of 7200 s.
Seawater concentration/wt% The dynamic responses of TBT, BBT, and T B1 , T B2 , T B3 are shown in Figures 13-15. The figures show that the rise of the above parameters is smaller than that of the disturbance. The increase in T sea reduces the overall temperature drop of the device, indicating a decrease in the heat consumption of the flashing brine in the system. According to the law of conservation of energy, the heat lost by the flashing brine will be converted into the heat absorbed by the flashing steam, resulting in less distillate. Therefore, as shown in Figure 16, when the disturbance occurs, W D1 , W D2 , W D3 will first decrease and then increase.

Dynamic Response Analysis of Feed Seawater Salt Concentration Change
C F is one of the important parameters that are susceptible to external influences in the MSF-BR process, and its change will also affect the performance of the entire device. Under the condition that other operating parameters remain unchanged, C F is increased (from 5.7 to 6.7 wt %) during the simulation process.
As shown in Figure 17, when the system runs normally to 4000 s, a disturbance is caused by the increase of C F by 300 s; that is, C F increases from 5.7 to 6.7 wt %. When the system runs to 4300 s, the disturbance stops and C F recovers to 5.7 wt % and this lasts until the end of 7200 s.
Under the condition that other operating parameters remain unchanged, is increased (from 5.7 to 6.7 wt %) during the simulation process.
As shown in Figure 17, when the system runs normally to 4000 s, a disturbance is caused by the increase of F C by 300 s; that is, F C increases from 5.7 to 6.7 wt %. When the system runs to 4300 s, the disturbance stops and F C recovers to 5.7 wt % and this lasts until the end of 7200 s. T T T decrease and then increase, with little effect. As shown in Figure 19, the variation trend of BBT is consistent with the trend of the disturbance of F C , with little influence.  T T T decrease and then increase, with little effect. As shown in Figure 19, the variation trend of BBT is consistent with the trend of the disturbance of F C , with little influence.     As shown in Figure 21, due to the first TBT decrease and then increase, the overall flash temperature drop of the system will decrease and then increase accordingly. According to the conservation of energy,

Dynamic Response Analysis of Recycle Stream Mass Flow Rate Change
Re W is one of the important parameters in the MSF-BR process, and it is also one of the adjustable parameters. The study of the influence of Re W on dynamic performance can provide a deeper understanding of the dynamic characteristics of the system and provide a path to dynamic optimization of the system in the future. Under the condition that other operating parameters remain unchanged, Re W is increased (from 6.05 × 10 6 to 6.35 × 10 6 kg/h) during the simulation process. As shown in Figure 22, when the system runs normally to 4000 s, a disturbance is caused by the increase of Re W by 300 s; that is, Re W increases from 6.05 × 10 6 to 6.35 × 10 6 kg/h. When the system runs to 4300 s, the disturbance stops and Re W recovers to 6.05 × C F affects the related physical parameters of the system, thus affecting system performance. It can be seen from Equation (32) that the change in C F will lead to the change in the BPE j value, which will affect the T Bj of the flash chambers. The change in T Bj will affect the flash temperature drop of the whole system and then affect distillate production.
Therefore, as shown in Figures 18-20, the first increase and then decrease of the feed seawater salt concentration made TBT and T B1 , T B2 , T B3 decrease and then increase, with little effect. As shown in Figure 19, the variation trend of BBT is consistent with the trend of the disturbance of C F , with little influence.
As shown in Figure 21, due to the first TBT decrease and then increase, the overall flash temperature drop of the system will decrease and then increase accordingly. According to the conservation of energy, W D1 , W D2 , W D3 also decreases first and then increases.

Dynamic Response Analysis of Recycle Stream Mass Flow Rate Change
W Re is one of the important parameters in the MSF-BR process, and it is also one of the adjustable parameters. The study of the influence of W Re on dynamic performance can provide a deeper understanding of the dynamic characteristics of the system and provide a path to dynamic optimization of the system in the future. Under the condition that other operating parameters remain unchanged, W Re is increased (from 6.05 × 10 6 to 6.35 × 10 6 kg/h) during the simulation process.
As shown in Figure 22, when the system runs normally to 4000 s, a disturbance is caused by the increase of W Re by 300 s; that is, W Re increases from 6.05 × 10 6 to 6.35 × 10 6 kg/h. When the system runs to 4300 s, the disturbance stops and W Re recovers to 6.05 × 10 6 kg/h and this lasts until the end of 7200 s. other operating parameters remain unchanged, Re W is increased (from 6.05 × 10 6 to 6.35 × 10 6 kg/h) during the simulation process. As shown in Figure 22, when the system runs normally to 4000 s, a disturbance is caused by the increase of Re W by 300 s; that is, Re W increases from 6.05 × 10 6 to 6.35 × 10 6 kg/h. When the system runs to 4300 s, the disturbance stops and Re W recovers to 6.05 × 10 6 kg/h and this lasts until the end of 7200 s. The increase in Re W has little effect on BBT. For the later flash chamber of the system, the increase in Re W is equivalent to the increase in Bj W in the flash chamber. Under the effect of the fixed interstage pressure difference of the system, Bj T cannot drop to the previous ideal temperature but may increase instead. As shown in Figure 24, the first increase and then decrease of Re W will cause BBT to decrease first and then increase.  The increase in Re W has little effect on BBT. For the later flash chamber of the system, the increase in Re W is equivalent to the increase in Bj W in the flash chamber. Under the effect of the fixed interstage pressure difference of the system, Bj T cannot drop to the previous ideal temperature but may increase instead. As shown in Figure 24, the first increase and then decrease of Re W will cause BBT to decrease first and then increase.   As shown in Figure 26, because Re W first decreases and then increases, the overall seawater of the device increases first and then decreases. Hence,

Dynamic Response Analysis of Steam Temperature Change
The temperature of steam directly affects the size of TBT and then affects the performance of the whole device. Under the condition that other operating parameters remain unchanged, steam T increases (from 97 to 99 °C) during the simulation process. As shown in Figure 27, when the system runs normally to 4000 s, a disturbance is caused by the increase of steam T ; that is, steam T increases from 97 to 99 °C and this lasts until the end of 7200 s.  As shown in Figure 26, because Re W first decreases and then increases, the overall seawater of the device increases first and then decreases. Hence,

Dynamic Response Analysis of Steam Temperature Change
The temperature of steam directly affects the size of TBT and then affects the performance of the whole device. Under the condition that other operating parameters remain unchanged, steam T increases (from 97 to 99 °C) during the simulation process. As shown in Figure 27, when the system runs normally to 4000 s, a disturbance is caused by the increase of steam T ; that is, steam T increases from 97 to 99 °C and this lasts until the end of 7200 s.  W Re is mixed with the make-up brine mass flow rate (W m ) from the heat rejection section and enters the condenser tubes of the heat recovery section for preheating. After mixing, the cooling brine mass flow rate to the recovery section (W R ) will then enter the brine heater for heating. Obviously, the increase in W Re will reduce the temperature of the cooling brine to the brine heater (T F0 ). As shown in Figures 23-25, the first increase and then decrease in W Re will lead to the decrease and then increase of T B1 , T B2 , T B3 and TBT.
The increase in W Re has little effect on BBT. For the later flash chamber of the system, the increase in W Re is equivalent to the increase in W Bj in the flash chamber. Under the effect of the fixed interstage pressure difference of the system, T Bj cannot drop to the previous ideal temperature but may increase instead. As shown in Figure 24, the first increase and then decrease of W Re will cause BBT to decrease first and then increase.
As shown in Figure 26, because W Re first decreases and then increases, the overall seawater of the device increases first and then decreases. Hence, W D1 , W D2 , W D3 increases first and then decreases.

Dynamic Response Analysis of Steam Temperature Change
The temperature of steam directly affects the size of TBT and then affects the performance of the whole device. Under the condition that other operating parameters remain unchanged, T steam increases (from 97 to 99 • C) during the simulation process.
As shown in Figure 27, when the system runs normally to 4000 s, a disturbance is caused by the increase of T steam ; that is, T steam increases from 97 to 99 • C and this lasts until the end of 7200 s.     Due to the increase of TBT, the flash temperature drop of the system will increase. Hence, as shown in Figure 31, under the condition that other operating parameters remain unchanged,     Due to the increase of TBT, the flash temperature drop of the system will increase. Hence, as shown in Figure 31, under the condition that other operating parameters remain unchanged,   Cooling seawater (W F0 ) enters the brine heater from the outlet of the condenser tubes of the first-stage flash chamber, and W steam releases latent heat to heat the brine in the device. At the outlet of the brine heater, W B0 enters the brine pool of the first-stage flash chamber at TBT to start the flashing. T steam directly affects TBT and T B1 , T B2 , T B3 . As shown in Figures 28-30, with the increase in T steam , both TBT and T B1 , T B2 , T B3 increase, which has a great influence. Compared with TBT, the increase in T steam has little influence on BBT.

Conclusions
Due to the increase of TBT, the flash temperature drop of the system will increase. Hence, as shown in Figure 31, under the condition that other operating parameters remain unchanged, W D1 , W D2 , W D3 will increase.

Conclusions
Multistage flash is one of the most widely used desalination methods, and its dynamic performance is very important for optimal operation, design, maintenance, and fault diagnosis. In this paper, the dynamic modeling and simulation of the multistage flash process were studied for the purpose of easy implementation and real-time operation. Based on more detailed consideration, a whole-process dynamic model, including a brine heat exchanger, a flashing stage room, mixed and split modulates, and physical parameter modulate, was established. For quick simulation, the dynamic equations were discretized and put together with physical equations through a simultaneous strategy and solved using the Interior Point solver, with sparsity treatment under the MATLAB platform. Then, the validation and solving efficiencies were verified by a multistage flash desalination system with 16 stages. In addition, the dynamic response for different key parameters such as feed seawater temperature, feed seawater concentration, recycle stream mass flow rate, and steam temperature were simulated and analyzed. The simulation results can also provide ideas for the combination of the variable changes that make the system operate with more stability. For example, the distillate flow rate decreases with the increase of feedwater concentration but increases with the steam temperature; therefore, it could be reasonable to think that increasing the steam temperature when the feedwater concentration increases is a kind of control strategy to maintain the levels of distillate in the plant. The states of each stage and the key permeance curve of TBT and BBT were also obtained and analyzed under dynamic change. The presented dynamic model and its treatment can provide a better understanding and analysis for real-time optimal operation and control of the MSF system to achieve lower operational cost and more stable freshwater quality.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the privacy.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
In order to better understand the mathematical model of MSF seawater desalination, a symbol description list is added here.

A H
Heat transfer area of the brine heater, m 2 A J Heat transfer area of stage j, m 2 BPE j Boiling point elevation of stage j C B0 Salt concentration in the flashing leaving the brine heater, wt % C Bj Salt concentration in the flashing brine leaving stage j, wt % C F0 Salt concentration of flashing seawater leaving brine heater, wt % C Fj Salt concentration of flashing seawater leaving stage j, wt % C m Salt concentration in make-up water, wt % C R Salt concentration in the cooling brine to the recovery section, wt % C Re Recycle brine concentration, wt % C W Rejected seawater mass flow rate, kg·h −1 CP Bj Heat capacity of brine leaving stage j, kcal·(kg· • C) −1

CP Fj
Heat capacity of flashing seawater leaving stage j, kcal·(kg· • C) −1 CP V j Heat capacity of flashing steam leaving stage j, kcal·(kg· • C) −1 CP RH Heat capacity of cooling brine leaving brine heater, kcal·(kg· • C) −