CFD Investigation of Spacer-Filled Channels for Membrane Distillation

The membrane distillation (MD) process for water desalination is affected by temperature polarization, which reduces the driving force and the efficiency of the process. To counteract this phenomenon, spacer-filled channels are used, which enhance mixing and heat transfer but also cause higher pressure drops. Therefore, in the design of MD modules, the choice of the spacer is crucial for process efficiency. In the present work, different overlapped spacers are investigated by computational fluid dynamics (CFD) and results are compared with experiments carried out with thermochromic liquid crystals (TLC). Results are reported for different flow attack angles and for Reynolds numbers (Re) ranging from ~200 to ~800. A good qualitative agreement between simulations and experiments can be observed for the areal distribution of the normalized heat transfer coefficient. Trends of the average heat transfer coefficient are reported as functions of Re for the geometries investigated, thus providing the basis for CFD-based correlations to be used in higher-scale process models.


Introduction
In regions where the existing freshwater sources do not meet the water demand, production of potable water from saline and brackish waters is needed [1]. Among the desalination processes, membrane distillation (MD) has recently attracted significant interest due to its features, which make MD a process easily scalable to small-medium size [2][3][4]. MD is a thermally driven membrane-based process where the transport of water vapor is realized through a micro-porous hydrophobic membrane, which leads to a theoretical 100% rejection of salts. The separation process is driven by the vapor pressure difference at the membrane/liquid interfaces due to the temperature gradient existing between the two sides of the membrane. MD is characterized by its compactness and low operational pressure requirement as well as low working temperature. These features make it possible to power MD by low grade heat (waste heat) or renewable energy and boost its applicability in remote areas and small scale production [5][6][7]. Another important distinct feature of MD is its ability to desalinate highly saline brines and wastewater. MD is the most appropriate option to increase the recovery ratio of multi-stage flash (MSF) and reverse osmosis (RO) plants due to it its ability to treat feeds with high salinity and impurities. One main advantage of MD technology is that evaporation and condensation surfaces are tightly packed, thus leading to a compact system with low capital cost per unit product. The classification of MD systems is related to the adopted condensation methods. In general, MD systems are classified into four basic simple configurations designed relative to the permeate vapor condensation method used. These include direct contact membrane distillation (DCMD), air gap [18,36,37], was not exactly equal to 2d (twice the filament diameter), i.e., 4 mm, but slightly less (3.8 mm) because a certain amount of interpenetration existed. Thus, the pitch to channel height ratio was P/H = 2.63.   Both in the experiments and in the computations, the bottom wall was assumed adiabatic and heat transfer occurring only through the top one (active wall), in which distributions of temperature, wall heat flux, and heat transfer coefficient were assessed. This last quantity was locally defined as: in which Tbulk is the local fluid bulk temperature, while q" and Twall are the local wall heat flux and wall temperature, respectively. For the reasons discussed in [23], an average heat transfer coefficient is better defined not as the surface average h of the local h, but rather as in which the symbol   denotes averaging operation on the active wall. The Reynolds number (Re) is defined as:  Both in the experiments and in the computations, the bottom wall was assumed adiabatic and heat transfer occurring only through the top one (active wall), in which distributions of temperature, wall heat flux, and heat transfer coefficient were assessed. This last quantity was locally defined as: (1) in which T bulk is the local fluid bulk temperature, while q and T wall are the local wall heat flux and wall temperature, respectively. For the reasons discussed in [23], an average heat transfer coefficient is better defined not as the surface average h of the local h, but rather as (2) in which the symbol denotes averaging operation on the active wall. The Reynolds number (Re) is defined as: where D h is the hydraulic diameter, conventionally assumed to be equal to twice the channel height (7.6 mm) as in a void plane channel of infinite width. U is the approach, or superficial, velocity, defined as the velocity that the fluid would have if the channel were void of any spacer; it is equal to the ratio between the flow rate and the passage area U = Q/A p (in the experiments A p = 9.12 × 10 −4 m 2 ). ν is the kinematic viscosity, which was assumed equal to 6.78 × 10 −7 m 2 /s in the cases investigated (water at 43 • C). For each configuration, the values of the flow rate investigated and the corresponding values of the Reynolds number are those given in Table 2. Flow rates were selected based on the experiments' limitations (the range of the flowmeters, measurement uncertainty at low flow rates, and pressure build up at the Plexiglas channel), and to enable comparison of results with several studies conducted in the literature by other researchers within the same range of values.

Experimental Set-Up and Procedures
The experimental results were obtained in a previous study [35] using the setup shown in Figure 2. It is composed of two narrow (~4 mm height) hot and cold water channels separated by a 1 mm thick polycarbonate sheet. These channels are enclosed by an upper and a lower 2 cm thick transparent Plexiglass sheet. Test spacers were custom designed and fabricated using 3D printing with overall dimensions of 24 cm × 24 cm × 0.38 cm. The spacers were designed to fit and fill the hot channel while touching a thermochromic liquid crystals (TLC) sheet, which was made to adhere to the lower surface of the polycarbonate sheet. TLC sheets adopted had a color play from 35 to 40 • C. The hot and cold water temperatures were measured using RTD sensors located at the inlet and outlet of the two channels, while the flow rates were measured using flow sensors and rotameters. Julabo temperature-controlled heating and cooling circulators were used to control the hot and cold water temperatures at the channel inlets. Further details on experimental procedure and data processing can be found in [35]. The color distribution of the TLC at the upper surface of the hot channel was monitored by a high resolution camera, and images (1920 × 1080 pixels) were captured for different feed temperatures and flow rates (from 1 to 4 L/min), while the flow rate in the cold channel was kept at 8 L/min. where Dh is the hydraulic diameter, conventionally assumed to be equal to twice the channel height (7.6 mm) as in a void plane channel of infinite width. U is the approach, or superficial, velocity, defined as the velocity that the fluid would have if the channel were void of any spacer; it is equal to the ratio between the flow rate and the passage area U = Q/Ap (in the experiments Ap = 9.12 × 10 −4 m 2 ). ν is the kinematic viscosity, which was assumed equal to 6.78 × 10 −7 m 2 /s in the cases investigated (water at ∼43 °C). For each configuration, the values of the flow rate investigated and the corresponding values of the Reynolds number are those given in Table 2. Flow rates were selected based on the experiments' limitations (the range of the flowmeters, measurement uncertainty at low flow rates, and pressure build up at the Plexiglas channel), and to enable comparison of results with several studies conducted in the literature by other researchers within the same range of values.

Experimental Set-Up and Procedures
The experimental results were obtained in a previous study [35] using the setup shown in Figure 2. It is composed of two narrow (∼4 mm height) hot and cold water channels separated by a 1 mm thick polycarbonate sheet. These channels are enclosed by an upper and a lower 2 cm thick transparent Plexiglass sheet. Test spacers were custom designed and fabricated using 3D printing with overall dimensions of 24 cm × 24 cm × 0.38 cm. The spacers were designed to fit and fill the hot channel while touching a thermochromic liquid crystals (TLC) sheet, which was made to adhere to the lower surface of the polycarbonate sheet. TLC sheets adopted had a color play from 35 to 40 °C. The hot and cold water temperatures were measured using RTD sensors located at the inlet and outlet of the two channels, while the flow rates were measured using flow sensors and rotameters. Julabo temperaturecontrolled heating and cooling circulators were used to control the hot and cold water temperatures at the channel inlets. Further details on experimental procedure and data processing can be found in [35]. The color distribution of the TLC at the upper surface of the hot channel was monitored by a high resolution camera, and images (1920 × 1080 pixels) were captured for different feed temperatures and flow rates (from 1 to 4 L/min), while the flow rate in the cold channel was kept at 8 L/min.  These images were processed by first converting them to the HSV (hue, saturation, value) format. It is known that the hue component of this format correlates with temperature but requires prior calibration. At the beginning of each spacer experiment, a TLC calibration test was preliminarily performed and a best fit polynomial curve was obtained for the temperature as a function of hue. During the calibration test both the upper and lower channels were fed by the same hot water source. In the real experiments, the cold water was set to 30 • C and the hot water was set to 43 • C, allowing a 13 • C difference between the two channels for heat transfer. Hot and cold water inlet temperatures were selected after several pre-tests to ensure that the temperature distributions on the TCL surface were clearly visible and fell within the color play range for all the flow rates.
It should be observed that in order to characterize the temperature polarization phenomenon, it is sufficient to reproduce the thermal field and convective heat transfer occurring in the spacer-filled channel only. This can be achieved by imposing at the active wall a realistic thermal boundary condition, representative of those occurring in MD modules, independent of whether there are or not an actual membrane and a vapor flux.

CFD Investigations
The problem was described by the continuity, momentum (Navier-Stokes), and energy equations for water at T = 43 • C and p = 1 bar, which were solved by the ANSYS-CFX ® code (ANSYS, USA). In the simulations, the unit periodic cell approach was used [22,34,38], which simulates flow and heat exchange phenomena in periodic lattices under fully developed conditions (i.e., at sufficient distance from inlets and outlets). A source term, accounting for the large-scale temperature gradient, was implemented in the energy equation, and a body force per unit volume, accounting for the large-scale pressure gradient, in the momentum equations. In each run, this latter term was dynamically adjusted to obtain the required Re value (corresponding to one of those achieved in the experiments).
As reported in the literature, for geometries similar to those examined here [19,26], the fluid flow becomes unsteady for Re > 350. For this reason, steady-state laminar simulations were carried out for Re = 205 and 305, while, for 410 ≤ Re ≤ 820, the shear stress transport (SST) turbulence model was used, which is a blend between the k-ω model near the walls and the k-ε model in the outer region. As demonstrated in previous works [29], for the present transitional flows, ω-based models, which fully resolve the near-wall layer, are preferable to ε-based models, which make use of wall functions; among them, the SST model gives the most accurate predictions in terms of both distributions and average values of the heat transfer coefficient.
The computational domains (different according to the intrinsic angle) are shown in Figure 3, while details of the computational finite volume grid for case (c) are shown in Figure 4. Only four corner blocks were meshed with tetrahedra; the remaining volume was meshed with mapped hexahedra or with sweep mode in the regions adjacent to the cylindrical surfaces. In regard to the boundary conditions, periodicity was imposed at the walls perpendicular to the filaments; at the other walls a no slip condition was set. The fluid-filament interfaces and the bottom wall were assumed adiabatic, while at the top wall a mixed condition (Robin) was imposed by prescribing an external temperature of 30 • C and an external heat transfer coefficient equal to λ/δ, λ being the fluid's thermal conductivity and δ the channel thickness.
hexahedra or with sweep mode in the regions adjacent to the cylindrical surfaces. In regard to the boundary conditions, periodicity was imposed at the walls perpendicular to the filaments; at the other walls a no slip condition was set. The fluid-filament interfaces and the bottom wall were assumed adiabatic, while at the top wall a mixed condition (Robin) was imposed by prescribing an external temperature of 30 °C and an external heat transfer coefficient equal to λ/δ, λ being the fluid's thermal conductivity and δ the channel thickness.

CFD Prediction of Temperature and Heat Transfer Coefficient Distributions and Comparison with Experiments
For the sake of brevity, results are illustrated here in detail only for the two geometrical configurations in Figure 5. Configuration 1 is characterized by an intrinsic angle θ between filaments of 60° and a flow attack angle α of 30°, while configuration 2 is characterized by θ = 90° and α = 60°.
The solid profiles indicate the spacer filaments, which touch the upper, thermally active, wall. The results of the whole computational campaign (aiming at the determination of the overall dependence of h on Re for each geometrical configuration, to be used in process modelling applications) are reported in the final section. For each configuration, three Reynolds numbers were considered: the lowest (205), an intermediate one (410), and the highest (820) of those reported in Table 2. Distributions of h, normalized by the mean value h, are shown. Experimental distributions were obtained by postprocessing the TLC images, as described in Section 2.2, and averaging the results over 3 × 3 unit cells to reduce dispersion.

CFD Prediction of Temperature and Heat Transfer Coefficient Distributions and Comparison with Experiments
For the sake of brevity, results are illustrated here in detail only for the two geometrical configurations in Figure 5. Configuration 1 is characterized by an intrinsic angle θ between filaments of 60 • and a flow attack angle α of 30 • , while configuration 2 is characterized by θ = 90 • and α = 60 • . The solid profiles indicate the spacer filaments, which touch the upper, thermally active, wall. The results of the whole computational campaign (aiming at the determination of the overall dependence of h on Re for each geometrical configuration, to be used in process modelling applications) are reported in the final section.

CFD Prediction of Temperature and Heat Transfer Coefficient Distributions and Comparison with Experiments
For the sake of brevity, results are illustrated here in detail only for the two geometrical configurations in Figure 5. Configuration 1 is characterized by an intrinsic angle θ between filaments of 60° and a flow attack angle α of 30°, while configuration 2 is characterized by θ = 90° and α = 60°.
The solid profiles indicate the spacer filaments, which touch the upper, thermally active, wall. The results of the whole computational campaign (aiming at the determination of the overall dependence of h on Re for each geometrical configuration, to be used in process modelling applications) are reported in the final section. For each configuration, three Reynolds numbers were considered: the lowest (205), an intermediate one (410), and the highest (820) of those reported in Table 2. Distributions of h, normalized by the mean value h , are shown. Experimental distributions were obtained by post-processing the TLC images, as described in Section 2.2, and averaging the results over 3 × 3 unit cells to reduce dispersion.

Configuration 1
Results for configuration 1 are reported in Figure 6. For clarity purposes, the insets show the direction of the flow and the arrangement of the filaments. Results for configuration 1 are reported in Figure 6. For clarity purposes, the insets show the direction of the flow and the arrangement of the filaments. The comparison shows a fair overall agreement between simulations and experiments in the heat transfer coefficient distribution. At the lowest Re, the h distribution was roughly symmetric between the upstream (left) and downstream (right) halves of the cell, showing that inertial effects were negligible. As Re increased, a stripe of low h values grew shortly after the upstream upper filament, associated with a region of separated flow (as confirmed by an analysis of the CFDpredicted flow fields). At the same time, a band of high h values appeared shortly before the downstream upper filament, in correspondence with the reattachment of the separated shear layer. negligible. As Re increased, a stripe of low h values grew shortly after the upstream upper filament, associated with a region of separated flow (as confirmed by an analysis of the CFD-predicted flow fields). At the same time, a band of high h values appeared shortly before the downstream upper filament, in correspondence with the reattachment of the separated shear layer. Some discrepancy existed between experiments and CFD results, which predicted stripes of high h, narrower, longer, and inclined with respect to the filaments. Note also that experimental results showed residual fluctuations, not eliminated by the spatial averaging performed (which only involved nine unit cells); CFD simulations, by their nature themselves, yielded time-averaged values of h (and of any other variable) and thus appeared smooth and regular.

Configuration 2
Results for configuration 2 are reported in Figure 7. As in the previous Figure 6, the inset shows the direction of the flow and the arrangement of the filaments. Some discrepancy existed between experiments and CFD results, which predicted stripes of high h, narrower, longer, and inclined with respect to the filaments. Note also that experimental results showed residual fluctuations, not eliminated by the spatial averaging performed (which only involved nine unit cells); CFD simulations, by their nature themselves, yielded time-averaged values of h (and of any other variable) and thus appeared smooth and regular.

Configuration 2
Results for configuration 2 are reported in Figure 7. As in the previous Figure 6, the inset shows the direction of the flow and the arrangement of the filaments. In this case only a partial agreement between simulation and experiments could be registered. The numerical simulations predicted, both under laminar flow conditions (Figure 7b) and using the turbulence model (Figure 7d,f), a band where the most intense heat transfer occurs, located shortly In this case only a partial agreement between simulation and experiments could be registered. The numerical simulations predicted, both under laminar flow conditions (Figure 7b) and using the turbulence model (Figure 7d,f), a band where the most intense heat transfer occurs, located shortly before the downstream (upper) spacer filament, and a secondary band located shortly upstream of it (i.e., closer to the center of the unit cell). These features were not present in the experimental results, which showed a broad maximum of h only slightly upstream of the upper filament and centered on the opposite filament (the one that did not touch the active upper wall).
General levels and the overall distribution of h/h avg , however, were satisfactorily reproduced. As in the previous configuration, experimental results exhibited a considerable amount of residual fluctuations, not completely damped by the 9-cell spatial averaging performed, while CFD results were smooth and regular. before the downstream (upper) spacer filament, and a secondary band located shortly upstream of it (i.e., closer to the center of the unit cell). These features were not present in the experimental results, which showed a broad maximum of h only slightly upstream of the upper filament and centered on the opposite filament (the one that did not touch the active upper wall). General levels and the overall distribution of h/havg, however, were satisfactorily reproduced. As in the previous configuration, experimental results exhibited a considerable amount of residual fluctuations, not completely damped by the 9-cell spatial averaging performed, while CFD results were smooth and regular.   A fair agreement between the experiments and the CFD simulations could be noticed. CFD results for configuration 1 (Figure 8) overpredicted havg by 13% when Re was lower than 600; at higher Re, a slight underprediction was registered (14%). In Figure 9, results for configuration 2 showed that CFD overestimates havg by 14%, with the maximum difference between experiments and simulations (25% overprediction) observed at Re ≈ 500-600. before the downstream (upper) spacer filament, and a secondary band located shortly upstream of it (i.e., closer to the center of the unit cell). These features were not present in the experimental results, which showed a broad maximum of h only slightly upstream of the upper filament and centered on the opposite filament (the one that did not touch the active upper wall). General levels and the overall distribution of h/havg, however, were satisfactorily reproduced. As in the previous configuration, experimental results exhibited a considerable amount of residual fluctuations, not completely damped by the 9-cell spatial averaging performed, while CFD results were smooth and regular.   A fair agreement between the experiments and the CFD simulations could be noticed. CFD results for configuration 1 (Figure 8) overpredicted havg by 13% when Re was lower than 600; at higher Re, a slight underprediction was registered (14%). In Figure 9, results for configuration 2 showed that CFD overestimates havg by 14%, with the maximum difference between experiments and simulations (25% overprediction) observed at Re ≈ 500-600. A fair agreement between the experiments and the CFD simulations could be noticed. CFD results for configuration 1 (Figure 8) overpredicted h avg by~13% when Re was lower than~600; at higher Re, a slight underprediction was registered (~14%). In Figure 9, results for configuration 2 showed that CFD overestimates h avg by~14%, with the maximum difference between experiments and simulations (~25% overprediction) observed at Re ≈ 500-600.

Summary Results for All Geometries Investigated
In Figures 10 and 11, the CFD results obtained for all the configurations reported in Table 1 are shown. In Figure 10, h avg is plotted as a function of Reynolds number for the three configurations in which the flow direction bisects the intrinsic angle θ, i.e., θ = 30 • , α = 15 • ; θ = 60 • , α = 30 • ; and θ = 90 • , α = 45 • . Curves of h avg against Re, for the configuration with θ = 90 • , are reported in Figure 11 considering flow attack angles values lower (graph a) or higher (graph b) than 45 • .
From Figures 10 and 11, it is possible to identify the configuration with θ = 90 • , α = 60 • as that providing the highest heat transfer coefficient.

Summary Results for All Geometries Investigated
In Figures 10 and 11, the CFD results obtained for all the configurations reported in Table 1 are shown. In Figure 10, havg is plotted as a function of Reynolds number for the three configurations in which the flow direction bisects the intrinsic angle θ, i.e., θ = 30°, α = 15°; θ = 60°, α = 30°; and θ = 90°, α = 45°. Curves of havg against Re, for the configuration with θ = 90°, are reported in Figure 11 considering flow attack angles values lower (graph a) or higher (graph b) than 45°.
From Figures 10 and 11, it is possible to identify the configuration with θ = 90°, α = 60° as that providing the highest heat transfer coefficient.

Summary Results for All Geometries Investigated
In Figures 10 and 11, the CFD results obtained for all the configurations reported in Table 1 are shown. In Figure 10, havg is plotted as a function of Reynolds number for the three configurations in which the flow direction bisects the intrinsic angle θ, i.e., θ = 30°, α = 15°; θ = 60°, α = 30°; and θ = 90°, α = 45°. Curves of havg against Re, for the configuration with θ = 90°, are reported in Figure 11 considering flow attack angles values lower (graph a) or higher (graph b) than 45°.
From Figures 10 and 11, it is possible to identify the configuration with θ = 90°, α = 60° as that providing the highest heat transfer coefficient.

Conclusions
Flow and heat transfer were predicted by computational fluid dynamics for channels provided with spacers consisting of two overlapped layers with a fixed pitch to height ratio (~2.63) and different values of the intrinsic angle between the filaments and of the flow attack angle. The Reynolds number ranged between 200 and 820.
In a previous work, experiments were carried out for the same configurations and Reynolds numbers by using thermochromic liquid crystals (TLC), which provide the distribution of the heat transfer coefficient h on the active wall.
For two selected cases, distributions of h obtained by CFD were compared with the experimental distributions, obtained by post-processing the real images of the TLC. From the comparison, a fair agreement concerning the shape of the h distributions was registered; the match was better at low Reynolds number, when no turbulence model was used, and worse at high Re, when simulations relied on the SST turbulence model.
For the same test cases, the average heat transfer coefficient was also fairly well predicted by CFD simulations, with discrepancies of the order of 10-20%, which can be regarded as only minor in view of the complex geometry and flow structure.
Finally, the results of the whole computational campaign were reported in terms of average heat transfer coefficient as a function of the Reynolds number and can be used in process modelling applications involving spacers with the same geometrical features as those investigated here.
In MD systems, geometrical configurations yielding the highest possible heat transfer coefficient h should be adopted in order to minimize the membrane surface (the largest component of plant cost). With reference to the geometries investigated in this work, for Reynolds numbers above 400-500, this optimal configuration seems to be characterized by an intrinsic angle θ of 90 • and a flow attack angle α of 60 • .