Simulation Study on Direct Contact Membrane Distillation Modules for High-Concentration NaCl Solution

Membrane distillation technology, as a new membrane-based water treatment technology that combines the membrane technology and evaporation process, has the advantages of using low-grade heat, working at atmospheric pressure with simple configuration, etc. In this study, heat and mass transfer were coupled at the membrane surfaces through the user-defined function program. The effects of feed temperature, feed velocity and permeate velocity on temperature polarization were mainly investigated for a high-concentration NaCl solution. The temperature polarization was increased with the increase of feed temperature and the decrease of feed and permeate velocity. The effects of temperature, inlet velocity and solution concentration on the evaporation efficiency of the membrane module for co- and counter-current operations were investigated in detail. The counter-current operation performed better than co-current operation in most cases, except for the condition where the NaCl concentration was relatively low or the module length was long enough. In addition, the optimal membrane thickness for both PVDF and PTFE was studied. The optimal membrane thickness was found in the range of 10 to 20 μm, which corresponded to the highest permeate flux for the selected materials, pore size distribution, and operation conditions. Membrane material with lower thermal conductivity and larger porosity was prone to get higher permeate flux and had larger optimal membrane thickness. Increasing feed velocity or feed temperature could decrease the optimal membrane thickness.


Introduction
High-salt wastewater treatment attracts the attention of researchers and engineers in the field of wastewater treatment for production processes in chemical industries, seawater desalination, water recovery from desulfurization process in power plants, etc. In power plants, a large amount of desulfurization wastewater is produced in the process of flue gas desulfurization (FGD) and, in many countries, FGD wastewater is required to be recovered and reused due to environmental requirements and economic considerations. The process is usually composed of a pretreatment unit, after which wastewater contains mainly salt, and evaporation unit, where most of water is recovered accompanied

Theoretical Background
DCMD is the MD configuration in which both feed and permeate solutions are in direct contact with the hydrophobic porous membrane. The transmembrane temperature gradient serves as the driving force for mass transfer. As shown in Figure 1, heat was transferred from the feed solution to the heat boundary layer, and then the membrane surface, where water was vaporized. The heat loss from the feed solution was composed of the latent heat of vaporization, and the heat conduction of the membrane itself. The latter reduced the temperature difference between the two solutions. After, the vapor was diffused to the permeate side through the membrane. Vapor condensed at the interface. Conductive heat was transferred from the feed side to the permeate side of the membrane and then through the heat boundary layer.
Membranes 2020, 10, x FOR PEER REVIEW 3 of 20 DCMD is the MD configuration in which both feed and permeate solutions are in direct contact with the hydrophobic porous membrane. The transmembrane temperature gradient serves as the driving force for mass transfer. As shown in Figure 1, heat was transferred from the feed solution to the heat boundary layer, and then the membrane surface, where water was vaporized. The heat loss from the feed solution was composed of the latent heat of vaporization, and the heat conduction of the membrane itself. The latter reduced the temperature difference between the two solutions. After, the vapor was diffused to the permeate side through the membrane. Vapor condensed at the interface. Conductive heat was transferred from the feed side to the permeate side of the membrane and then through the heat boundary layer. The temperature polarization was schematically shown in Figure 1. Temperatures, Tfm, Tpm at the membrane-liquid interfaces were different from bulk solution temperature, Tf, Tp, caused by the resistance to the heat flux across the membrane boundary [26]. Furthermore, a concentration boundary layer was also formed near the membrane surface of the liquid layer, so that the membrane surface solute concentration was higher than the bulk solute concentration, i.e., the concentration polarization [18]. Temperature and concentration polarization reduced the mass transfer driving force, leading to a reduction in permeation flux, which needed to be considered in the simulation.
The heat transfer from feed solution to membrane surface of feed side was given by, The heat transfer across the membrane consists of two parts: the latent heat of vaporization and the heat conduction. The latter one was considered as energy loss. The evaporation efficiency (EE) was defined as the ratio of the efficient heat due to vaporization and the total heat transfer through the membrane and was calculated as, where QN and Qc could be calculated as, The temperature polarization was schematically shown in Figure 1. Temperatures, T fm , T pm at the membrane-liquid interfaces were different from bulk solution temperature, T f , T p , caused by the resistance to the heat flux across the membrane boundary [26]. Furthermore, a concentration boundary layer was also formed near the membrane surface of the liquid layer, so that the membrane surface solute concentration was higher than the bulk solute concentration, i.e., the concentration polarization [18]. Temperature and concentration polarization reduced the mass transfer driving force, leading to a reduction in permeation flux, which needed to be considered in the simulation.
The heat transfer from feed solution to membrane surface of feed side was given by, The heat transfer across the membrane consists of two parts: the latent heat of vaporization and the heat conduction. The latter one was considered as energy loss. The evaporation efficiency (EE) was defined as the ratio of the efficient heat due to vaporization and the total heat transfer through the membrane and was calculated as, where Q N and Q c could be calculated as, The heat transfer from the membrane surface of permeate side to permeate solution was given by, According to the conservation of energy, The latent heat of vaporization was obtained by [8], The thermal conductivity was obtained by [18], When T f ranged from 45 to 60 • C, the mean free path in DCMD was from 0.107 to 0.11 µm. The pore size of membrane ranged from 0.2 to 1.0 µm. Kn number ranged from 0.107 to 0.55. The mass transport took places via the combination of Knudsen diffusion and Molecular diffusion [38], and the mass transfer model was given by, The activity coefficient of water under different concentration conditions could be obtained by the following empirical formula [39]: where x s was the molar fraction of solute in the solution. D w−a was evaluated by the empirical formula, given by [40], The water vapor pressure P fm , and P pm were calculated through T fm , T pm by Antoine equation, In the MD process, the temperature polarization was measured with the temperature polarization coefficient (TPC), The concentration polarization was characterized using a concentration polarization coefficient (CPC), which was defined as the ratio of the membrane surface concentration to the feed inlet concentration, Membranes 2020, 10, 179 5 of 18

Geometry and Governing Equations
As shown in Figure 2, the geometry model was a rectangle flow channel which was 100 mm long, 2.5 mm high. The computational domain consisted of permeate and feed channels. Boundary conditions were set as shown in Figure 2. The inlet of permeate and feed channel were set to Velocity-Inlet. The outlet of permeate and feed channel were set to Pressure-Outlet. The membrane materials used in the simulation were shown in Table 1. The membrane module was placed horizontally. The NaCl solution was taken as feed solution, and pure water was used as the cooling substance. T f was set from 45 to 75 • C according to the fact that the desulfurization wastewater temperature was about 50 • C [41]. The temperature difference of feed and permeate solution was from 20 to 50 • C. The flow rate was set from 0.05 to 0.25 m/s. In Section 4.3, in order to investigate the temperature field and EE of co-and counter-current operations obviously, the calculation domain of x direction was extended by three times. conditions were set as shown in Figure 2. The inlet of permeate and feed channel were set to Velocity-Inlet. The outlet of permeate and feed channel were set to Pressure-Outlet. The membrane materials used in the simulation were shown in Table 1. The membrane module was placed horizontally. The NaCl solution was taken as feed solution, and pure water was used as the cooling substance. Tf was set from 45 to 75 °C according to the fact that the desulfurization wastewater temperature was about 50 °C [41]. The temperature difference of feed and permeate solution was from 20 to 50 °C. The flow rate was set from 0.05 to 0.25 m/s. In Section 4.3, in order to investigate the temperature field and EE of co-and counter-current operations obviously, the calculation domain of x direction was extended by three times.  In the current work, quadrilateral structured mesh was generated. In order to describe the fluid properties of the boundary layer near the membrane accurately, finer grids on both sides of the membrane were taken. The first layer of mesh on both sides of the film had a thickness of 10 μm and a growth factor of 1.05. Grids consisted of about 12,000. Grid independent analysis was undertaken to determine the grid sizes by examining the pressure drop of the modules via a hydrodynamic simulation.
The laminar model was chosen to describe the feed and permeate flows in the channels. The second-order upwind scheme was chosen to discrete the equations. The SIMPLE algorithm was used to solve the equations. The residuals of all the variables were set to 10 −6 .
The feed and permeate flow were governed by continuity equation, Navier-Stokes equations and energy equation, where was the viscous stress tensor, u, p, were the velocity vector, pressure, density, respectively, was dynamic viscosity, cp was the fluid heat capacity.
Mass transport was modeled using advection-diffusion equation,  In the current work, quadrilateral structured mesh was generated. In order to describe the fluid properties of the boundary layer near the membrane accurately, finer grids on both sides of the membrane were taken. The first layer of mesh on both sides of the film had a thickness of 10 µm and a growth factor of 1.05. Grids consisted of about 12,000. Grid independent analysis was undertaken to determine the grid sizes by examining the pressure drop of the modules via a hydrodynamic simulation.
The laminar model was chosen to describe the feed and permeate flows in the channels. The second-order upwind scheme was chosen to discrete the equations. The SIMPLE algorithm was used to solve the equations. The residuals of all the variables were set to 10 −6 .
The feed and permeate flow were governed by continuity equation, Navier-Stokes equations and energy equation, where σ was the viscous stress tensor, u, p, ρ were the velocity vector, pressure, density, respectively, µ was dynamic viscosity, c p was the fluid heat capacity.
Membranes 2020, 10, 179 6 of 18 Mass transport was modeled using advection-diffusion equation, where D was the effective mass diffusivity. Effective mass diffusivity D was obtained by [22], where ϑ Na was the equivalent limiting ionic conductance of the sodium and ϑ Cl was the equivalent limiting ionic conductance of the chloride. In order to simulate the transmembrane mass transfer process, it was necessary to load the mass source of water on the first layer grid near the permeate side and the feed side of the membrane. Besides, the transmembrane heat transfer process was achieved by loading the heat transfer source.
where b was the height of the first grid. For the membrane distillation process, the mass source term on the liquid side was not only related to the feed side grid temperature, but also affected by the corresponding permeate side grid temperature, and it was achieved by the user-defined function (UDF). The function of the UDF program was divided into three steps. First, the temperatures of cell center on both sides of the membrane wall was obtained by the UDF program, which were T fm , T pm , respectively. Then, the heat and mass transfer sources were calculated by UDF program through Equations (9)- (14). Finally, the UDF program was used to add the heat flux to membrane wall as boundary condition and mass transfer source as source terms.

Model Validation
The experimental data in the literature [18] were used to verify the model. In Figure 3, permeate fluxes were shown in different feed temperature for pure water and 24.2 wt.% NaCl solution, respectively. The simulation results agreed to the experiment data, which showed the results based on the current simulation model were reliable and the developed simulation methods were suitable for the following analysis.
where ϑNa was the equivalent limiting ionic conductance of the sodium and ϑCl was the equivalent limiting ionic conductance of the chloride. In order to simulate the transmembrane mass transfer process, it was necessary to load the mass source of water on the first layer grid near the permeate side and the feed side of the membrane. Besides, the transmembrane heat transfer process was achieved by loading the heat transfer source.
where b was the height of the first grid. For the membrane distillation process, the mass source term on the liquid side was not only related to the feed side grid temperature, but also affected by the corresponding permeate side grid temperature, and it was achieved by the user-defined function (UDF). The function of the UDF program was divided into three steps. First, the temperatures of cell center on both sides of the membrane wall was obtained by the UDF program, which were Tfm, Tpm, respectively. Then, the heat and mass transfer sources were calculated by UDF program through Equations (9)- (14). Finally, the UDF program was used to add the heat flux to membrane wall as boundary condition and mass transfer source as source terms.

Model Validation
The experimental data in the literature [18] were used to verify the model. In Figure 3, permeate fluxes were shown in different feed temperature for pure water and 24.2 wt.% NaCl solution, respectively. The simulation results agreed to the experiment data, which showed the results based on the current simulation model were reliable and the developed simulation methods were suitable for the following analysis.

Results and Discussion
In the current work, the influence of Tf, vf, vp on the temperature polarization in high-salt solution were simulated and analyzed. The effects of Tf, vf, vp and solution salinity on EE were studied. The optimal membrane thickness was discovered in different operation conditions.

Temperature Polarization
As shown in Figures 4 and 5, there was a decrease in both TPC and Tf over the membrane surface along the x-axis in each Tf. This was because that the temperature boundary layer became thicker

Results and Discussion
In the current work, the influence of T f , v f , v p on the temperature polarization in high-salt solution were simulated and analyzed. The effects of T f , v f , v p and solution salinity on EE were studied. The optimal membrane thickness was discovered in different operation conditions.

Temperature Polarization
As shown in Figures 4 and 5, there was a decrease in both TPC and T f over the membrane surface along the x-axis in each T f . This was because that the temperature boundary layer became thicker along the x direction, and the temperature difference between the membrane surface and the bulk area became larger, which caused a decrease of TPC. With the T f was increased from 45 to 75 • C, the TPC decreased correspondingly from 69.9% (x = 0.01 m) to 62.8% (x = 0.01 m). The variation of temperature difference between two sides of the membrane surface was smaller than that of the temperature difference between the bulk and membrane surface, which led to the decrease of TPC.
Membranes 2020, 10, x FOR PEER REVIEW 7 of 20 temperature difference between two sides of the membrane surface was smaller than that of the temperature difference between the bulk and membrane surface, which led to the decrease of TPC.  Ali's [42] study showed the theoretical and experimental TPC as a function of Reynolds number (Re) in the range of 1000 to 5000. When the Re was below 1000, which was exactly the Re of current study, TPC was around 0.68, showing good coincidence with the current study.
As shown in Figures 6 and 7, the temperature boundary layer became thinner with the increasing v f , v p . When the v f , v p were 0.05 m/s, TPC decreased from 60.36% at inlet to 41.87% at outlet and the temperature difference between inlet and outlet at feed side was 5.86 • C. When v f , v p were increased to 0.25 m/s, TPC decreased from 75.32% at inlet to 54.62% at outlet and the temperature difference between inlet and outlet at feed side was decreased to 4.56 • C. Results suggested that increasing the v f , v p could reduce temperature polarization significantly.

Co-and Counter-Current Operations
In this section, in order to well investigate the temperature field and EE of co-and counter-current operations, the calculation domain of x direction was extended by 3 times. The influence of v f , v p , T f and salinity on EE for co-and counter-current operation was evaluated in detail.   Ali's [42] study showed the theoretical and experimental TPC as a function of Reynolds number (Re) in the range of 1000 to 5000. When the Re was below 1000, which was exactly the Re of current study, TPC was around 0.68, showing good coincidence with the current study.
As shown in Figures 6 and 7, the temperature boundary layer became thinner with the increasing vf, vp. When the vf, vp were 0.05 m/s, TPC decreased from 60.36% at inlet to 41.87% at outlet and the temperature difference between inlet and outlet at feed side was 5.86 °C. When vf, vp were increased to 0.25 m/s, TPC decreased from 75.32% at inlet to 54.62% at outlet and the temperature difference between inlet and outlet at feed side was decreased to 4.56 °C. Results suggested that increasing the vf, vp could reduce temperature polarization significantly.

Temperature Fields
As shown in Figure 8, there was a decrease in T f and an increase in T p along the x-axis in co-and counter-current operations. In the studied cases, the feed side temperature difference between the inlet and outlet for co-and counter-current operations was 2.06 • C and 4.42 • C, respectively, which indicated that counter-current operation performed better than co-current. The temperature distribution of both feed and permeate sides along the x-axis was more homogeneous in counter-current operation, which resulted less thermal stress on the membrane. It was beneficial to prolong the working life of the membrane module. The temperature fields of count-current and co-current were well coincided with experimental data [43].

Co-and Counter-Current Operations
In this section, in order to well investigate the temperature field and EE of co-and countercurrent operations, the calculation domain of x direction was extended by 3 times. The influence of vf, vp, Tf and salinity on EE for co-and counter-current operation was evaluated in detail.

Temperature Fields
As shown in Figure 8, there was a decrease in Tf and an increase in Tp along the x-axis in co-and counter-current operations. In the studied cases, the feed side temperature difference between the inlet and outlet for co-and counter-current operations was 2.06 °C and 4.42 °C, respectively, which indicated that counter-current operation performed better than co-current. The temperature distribution of both feed and permeate sides along the x-axis was more homogeneous in countercurrent operation, which resulted less thermal stress on the membrane. It was beneficial to prolong the working life of the membrane module. The temperature fields of count-current and co-current were well coincided with experimental data [43].

Evaporation Efficiency for Co-and Counter-Current Operations
As shown in Figure 9, the EE was increased with the increase of T f in all cases. It could be concluded that temperature difference was the main influencing factor for EE. Besides, there was a decrease of EE over the membrane surface along the x-axis due to the decrease of temperature. It was essential to control the length of membrane module to avoid the effect. Comparing the co-and counter-current operations in different T f , it was found that the counter-current operation performed better than co-current at upstream at all the studied temperatures. The priority of counter-to co-current operation at 45 • C was more evident than at 75 • C. This phenomenon could be attributed to several factors. In general, the temperature difference between feed and permeate membrane surface in counter-current operation was larger than co-current operation. The higher temperature difference lead to higher EE. According to Lou's study [22], the concentration polarization of counter-current was less than that of co-current operation in the upstream and larger than that of co-current operation in the downstream.
As shown in Figure 10, the decrease in EE with x-axis occurred in all cases due to the reduced temperature difference across the membrane and increased c f . The EE for counter-current operation performed better than that for co-current operation in the same v f , v p . With the v f , v p increasing, the EE of both co-and counter-current increased. The difference of EE between counter-and co-current was decreasing along x-axis. When v f and v p were both 0.05 m/s, the EE for counter-current operation decreased slowly along the x-axis compared to co-current operation. This was due to the fact that the temperature different between the feed and permeate surface along x-axis was almost constant for counter-current operation and gradually decreased for co-current operation.
counter-current operations. In the studied cases, the feed side temperature difference between the inlet and outlet for co-and counter-current operations was 2.06 °C and 4.42 °C, respectively, which indicated that counter-current operation performed better than co-current. The temperature distribution of both feed and permeate sides along the x-axis was more homogeneous in countercurrent operation, which resulted less thermal stress on the membrane. It was beneficial to prolong the working life of the membrane module. The temperature fields of count-current and co-current were well coincided with experimental data [43].

Evaporation Efficiency for Co-and Counter-Current Operations
As shown in Figure 9, the EE was increased with the increase of Tf in all cases. It could be concluded that temperature difference was the main influencing factor for EE. Besides, there was a decrease of EE over the membrane surface along the x-axis due to the decrease of temperature. It was essential to control the length of membrane module to avoid the effect. Comparing the co-and counter-current operations in different Tf, it was found that the counter-current operation performed better than co-current at upstream at all the studied temperatures. The priority of counter-to cocurrent operation at 45 °C was more evident than at 75 °C. This phenomenon could be attributed to several factors. In general, the temperature difference between feed and permeate membrane surface in counter-current operation was larger than co-current operation. The higher temperature difference lead to higher EE. According to Lou's study [22], the concentration polarization of counter-current was less than that of co-current operation in the upstream and larger than that of co-current operation in the downstream. As shown in Figure 10, the decrease in EE with x-axis occurred in all cases due to the reduced temperature difference across the membrane and increased cf. The EE for counter-current operation performed better than that for co-current operation in the same vf, vp. With the vf, vp increasing, the EE of both co-and counter-current increased. The difference of EE between counter-and co-current was decreasing along x-axis. When vf and vp were both 0.05 m/s, the EE for counter-current operation As shown in Figure 11, it was found that the EE was generally low in high-salt solution due to the high concentration resulting in small permeate flux. It was worth noting that the EE for co-current operation started to higher than counter-current when the concentration was 2.9 wt.% at x = 0.15, which suggested that, in low concentrations, the impact of concentration polarization on EE was larger than temperature difference. As shown in Figure 11, it was found that the EE was generally low in high-salt solution due to the high concentration resulting in small permeate flux. It was worth noting that the EE for co-current operation started to higher than counter-current when the concentration was 2.9 wt.% at x = 0.15, which suggested that, in low concentrations, the impact of concentration polarization on EE was larger than temperature difference.

Optimal Membrane Thickness
Theoretically, two competitive effects are associated with membrane thickness. The thin membrane gives low mass transfer resistance and thus provide high transmembrane water flux at fixed operation condition. However, as the membrane thickness decreases, conductive heat loss

Optimal Membrane Thickness
Theoretically, two competitive effects are associated with membrane thickness. The thin membrane gives low mass transfer resistance and thus provide high transmembrane water flux at fixed operation condition. However, as the membrane thickness decreases, conductive heat loss through the membrane increases, leading to a low temperature gradient across the membrane and resulting low flux. Therefore, there must be an optimal thickness balancing the mass transfer resistance and conductive heat losses, and this thickness depends on operating conditions and other membrane characteristics.
Membrane thickness is an important property of membrane module but seldom studied [35,[44][45][46]. Martinez [46] established a resistance-in-series model which was quite suitable to express the change of membrane thickness due to the import of transport resistance. Swaminathan [44] coupled the membrane thickness with the module size together and got the cost-optimal membrane thickness. Wu figured out the optimal membrane thickness both by analytical model and by experiment. The results indicated that the optimal thickness was much thicker than that could be made in experiment. Eykens [35] made a systematic research on membrane thickness. The influence of membrane thickness on DCMD at various salinities was investigated in detail. However, there were no studies involving in comparing of different membranes on membrane thickness and permeate flux.
The current study investigated the influence of membrane pore size, conductivity, tortuosity and porosity through comparing the different membrane materials, i.e., TF450 and HVHP, had the same pore size while the membrane conductivity, tortuosity and porosity were different, GVHP and HVHP had the same membrane conductivity, tortuosity and porosity while the pore size was different. Figure 12 showed the flux as function of membrane thickness for different membrane materials. There existed a threshold of permeate flux. By decreasing the membrane thickness, the permeate flux could be improved until reached this threshold. In Wu's work [45], they found that the flux decreases with the decrease of membrane thickness through experiment, which was well reproduced in our simulations. The optimal membrane thickness and the largest permeate flux varied with different membrane materials. When the membrane thickness was kept as 15 µm, the permeate flux increased with the increasing of pore size for TF1000, TF450 and TF200, that was 20.30, 18.68, and 15.948 kg/(m 2 ·h), respectively. It indicated that the pore size had great impact on the permeate flux. The effect of thermal conductivity and porosity could be analyzed through TF450 and HVHP which had the same pore size, corresponding to permeate flux of 20.30 kg/(m 2 ·h) and 13.8 kg/(m 2 ·h). The thermal conductivity was 0.031 and 0.041 W·m −1 ·K −1 and the porosity was 0.8 and 0.75, respectively. Membrane material which had lower thermal conductivity and larger porosity showed higher permeate flux and thicker optimal membrane thickness. Lower thermal conductivity was helpful to reduce the heat loss, and larger porosity was beneficial to improving permeate flux. However, the optimal membrane thickness at 24.2 wt% was quite different with literature [30] (18 µm in the current work for GVHP and 49 µm in the literature for PP). The reason is that the tortuosity in our study is larger-the larger tortuosity means a thinner optimal membrane thickness. The membrane thermal conductivity is lower, the lower membrane conductivity also corresponds to a thinner optimal membrane thickness. Figure 13 demonstrated the impact of T f on the permeate flux and the membrane thickness for GVHP. Increase of permeate flux was observed when T f increased from 45 to 60 • C. It was clearly seen that the higher T f corresponded to the thinner membrane thickness. When T f was 45 • C, the optimal membrane thickness was 18 µm. However, when T f reached 60 • C, the optimal membrane thickness turned to 10 µm.

Optimal Membrane Thickness for Different T f
reduce the heat loss, and larger porosity was beneficial to improving permeate flux. However, the optimal membrane thickness at 24.2 wt% was quite different with literature [30] (18 μm in the current work for GVHP and 49 μm in the literature for PP). The reason is that the tortuosity in our study is larger-the larger tortuosity means a thinner optimal membrane thickness. The membrane thermal conductivity is lower, the lower membrane conductivity also corresponds to a thinner optimal membrane thickness.  Figure 13 demonstrated the impact of Tf on the permeate flux and the membrane thickness for GVHP. Increase of permeate flux was observed when Tf increased from 45 to 60 °C. It was clearly seen that the higher Tf corresponded to the thinner membrane thickness. When Tf was 45 °C, the optimal membrane thickness was 18 μm. However, when Tf reached 60 °C, the optimal membrane thickness turned to 10 μm. As shown in Figure 14, the permeate flux was improved with the increasing of vf, vp. The optimal membrane thickness differed with different vf, vp. When vf, vp was 0.1 m/s, the optimal membrane thickness was 20 μm, corresponding to permeate flux at 9.51 kg/(m 2 •h), while these data became 15

Optimal Membrane Thickness for Different v f , v p
As shown in Figure 14, the permeate flux was improved with the increasing of v f , v p . The optimal membrane thickness differed with different v f , v p . When v f , v p was 0.1 m/s, the optimal membrane thickness was 20 µm, corresponding to permeate flux at 9.51 kg/(m 2 ·h), while these data became 15 µm corresponding to 13.93 kg/(m 2 ·h) when v f , v p reached 0.2 m/s. The v f , v p affected the optimal membrane thickness to a certain extent, but its impact was small. The developed model could contribute to industrial application of MD, especially for supporting MD module designs in terms of optimization on membrane material and structure, working conditions or operating mode. In the current work, when the pore size varied from 0.2-1 μm, the Knudsen diffusion mechanism should be mainly considered for flux prediction, since the Knudsen diffusion coefficient is 50 times higher than that of molecular diffusion for the average pore size of 0.2 μm, and 10 times higher for the average pore size of 1 μm, when the membrane porosity is fixed. In the industries such as power plants, where a large amount of heat supply is present, the feed temperature should be properly increased to promote flux. In the industries such as for chemical synthesis with a low enthalpy of reaction, the selection of feed and permeate velocities, or pressure loss should be mainly considered for obtaining higher evaporation efficiency, where the model can be used for optimization by applying the conditions of the specific industries. Usually, countercurrent operation performs better than co-current operation in terms of evaporation efficiency.
Counter-current operation mode should be considered as priority. However, when the placement of MD module is geometrically limited and the aspect ratio of the membrane is high, counter-current mode is not necessarily better, the selection of co-current or counter-current mode can be determined by using simulation with the current model.

Conclusions
The simulation model coupling heat and mass transfer at the membrane surface through a userdefined function program was established. The optimal thickness for typical commercialized PVDF and PTFE membranes was studied, and it was found to be decreased with the increase of vf, vp or Tf . The TPC increased with the decrease of Tf or the increase of vf, vp. The counter-current operation The developed model could contribute to industrial application of MD, especially for supporting MD module designs in terms of optimization on membrane material and structure, working conditions or operating mode. In the current work, when the pore size varied from 0.2-1 µm, the Knudsen diffusion mechanism should be mainly considered for flux prediction, since the Knudsen diffusion coefficient is 50 times higher than that of molecular diffusion for the average pore size of 0.2 µm, and 10 times higher for the average pore size of 1 µm, when the membrane porosity is fixed. In the industries such as power plants, where a large amount of heat supply is present, the feed temperature should be properly increased to promote flux. In the industries such as for chemical synthesis with a low enthalpy of reaction, the selection of feed and permeate velocities, or pressure loss should be mainly considered for obtaining higher evaporation efficiency, where the model can be used for optimization by applying the conditions of the specific industries. Usually, counter-current operation performs better than co-current operation in terms of evaporation efficiency. Counter-current operation mode should be considered as priority. However, when the placement of MD module is geometrically limited and the aspect ratio of the membrane is high, counter-current mode is not necessarily better, the selection of co-current or counter-current mode can be determined by using simulation with the current model.

Conclusions
The simulation model coupling heat and mass transfer at the membrane surface through a user-defined function program was established. The optimal thickness for typical commercialized PVDF and PTFE membranes was studied, and it was found to be decreased with the increase of v f , v p or T f . The TPC increased with the decrease of T f or the increase of v f , v p . The counter-current operation mode performed better than co-current operation in most cases except for the condition where the c f was low and the module length was long enough. The difference of EE between co-and counter-current was due to the combined effect of temperature difference and concentration polarization. In addition, membranes with lower thermal conductivity and larger porosity were prone to get higher permeate flux and thinner optimal membrane thickness. The developed model could also provide guidance for MD module designs.

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