Computational Study of Zigzag Spacer Design with Elliptical Cross-Section Filaments

The large energy consumption of membrane desalination process has encouraged researchers to explore different spacer designs using Computational Fluid Dynamics (CFD) for maximizing permeate per unit of energy consumed. In previous studies of zigzag spacer designs, the filaments are modeled as circular cross sections in a two-dimensional geometry under the assumption that the flow is oriented normal to the filaments. In this work, we consider the 45° orientation of the flow towards the three-dimensional zigzag spacer unit, which projects the circular cross section of the filament as elliptical in a simplified twodimensional domain. OpenFOAM was used to simulate the mass transfer enhancement in a reverse-osmosis desalination unit employing spiral wound membranes lined with zigzag spacer filaments. Properties that impact the concentration polarization and hence permeate flux were analyzed in the domain with elliptical filaments as well as a domain with circular filaments to draw suitable comparisons. The range of variation in characteristic parameters across the domain between the two different configurations is determined. It was concluded that ignoring the elliptical projection of circular filaments to the flow direction, can introduce significant margin of error in the estimation of mass transfer coefficient.


Introduction
Reverse Osmosis desalination technique has gained pace over the years due to advancement in technology and abundance of sea-water, especially in the Middle East. The process utilizes high pressure fluids being pumped into devices that employ hydrophilic membranes to filter out the water into a permeate channel and retaining the concentrated water in the retentate channel. The excessively high pressures however translate into very high costs. Efforts are being constantly made to maximize the permeate obtained per unit of cost of filtration. Among several forms of modules used to carry out this process, spiral wound membrane (SWM) stands quite popular for industrial use because of its high membrane packing density compared to rectangular channels. SWM feature spacer nets throughout their length which serve two purposes. First, they provide mechanical support to multi-layered retentate and permeate channels and second, create instability in the flow field [1]. The primary ailment of the reverse osmosis procedure, however is the concentration polarization as reported by both Yun et al. [2] and Bahmanyar et al. [3]. Concentration polarization occurs when the concentration at the membrane is higher than that in the bulk fluid [4]. This not only acts as an added layer of resistance for water to flow into the permeate channel but also results in permeate water to flow back into the retentate channel, deteriorating the efficiency of the system. This is counteracted by the disruption of the boundary layer through intended turbulence promoters. The use of turbulence promoters was initiated in 1965 [5] and used for the first time a ladder type configuration of spacers in the channel. He reported significant gains in permeate flux and that the increased energy requirement due to enhanced pressure drop across the channel because of the spacers is outweighed by the increased mass transfer and hence permeate flux obtained as reported by Thomas et al. [5,6] and Da Costa et al. [7]. Yu et al. [8] appreciated the use of baffled channels in increasing mass transfer coefficients and Subramani et al [9] studied the pressure drops and shear rates in spacer filled and open channels and highlighted the usefulness of spacer filled channels. Santos et al. [10] extended the work of Geraldes and furthered identification of most beneficial spacer designs in SWM. They experimentally and numerically conducted experiments with twelve different spacers and reported on the laminar to turbulent transition. Da Costa et al. [11] conducted experiments on several different types and orientation of spacer configuration and reported that incoming flow at an angle to the horizontal of the spacers has greater mass transfer compared to one that comes at zero angle. They utilized circular strands instead of rectangular ones in their studies and concluded that the zigzag configuration shown in Figure 1 is the most efficient configuration of spacer designs. Traditionally, zigzag spacers involve approximating the strands cross-section as circular. This contrasts with the observed elliptical strands in zigzag configuration as depicted by Figure 1. Fimbers et al [12] conducted unsteady flow analyses of the zigzag spacers using circular strands and explored relations between vortex shedding, pressured drops and mass transfer coefficients. Similarly, Schwinge et al [13,14] modeled fluid flow behavior for different filament configurations including zigzag spacers in two dimensions. Wardeh et al [15] investigated the effect of zigzag spacers with circular strands and measured the permeate flux through the membrane.
Three dimensional modelling of spacer strands is very computationally intensive [1] and research is usually limited to simple flow modelling without incorporating mass transfer equations. Iwatsu et al. [16,17] simulated and analyzed the flow inside a cubic driven cavity to conclude that for steady flow at low Reynolds numbers, three dimensional can be reduced to two dimensional without any significant loss in the accuracy of results, however, as the Reynolds Number increases and the flow transitions towards unsteady behaviour, the endwall effects start to dominate which cannot be captured in two dimensional flows. Therefore, the study is only limited to low Reynolds number twodimensional flow.
Solid modelling was used to produce the threedimensional domain which was then reduced to two dimensional through a cross-sectional cut at 45° with the alignment of the strands as shown in Figure 1. It can be observed that the strands need to be modelled as elliptical if a zigzag domain is being simulated. To provide results that are comparable, the radius of the circles in the CC domain is also kept at 1.3mm and this same circle is projected on to a plane at 45° with the horizontal. This results in the minor axis length equal to the original radius, but the major axis is then a function of the radius and the cosine of the projection angle.

Geometry
Santos et al [10] established that CFD analysis can be limited to consecutive longitudinal filaments, hence reducing the computational resources required for the simulations. Therefore, three units of the dimensions shown by Fig. 2 were modelled back to back and the analyses were based on the second unit. This was done to provide sufficient time and length to the flow from the inlet to be completely developed and to prevent the exit boundary condition from impacting the zone of interest.
The case is considered to be incompressible. Density changes and effects of gravity are assumed to be negligible because of the narrow channel [18,19]. The fluid is Newtonian and kinematic viscosity is valued at 1.004 (m^2/s) at a temperature of 20 °C and a diffusion coefficient of 1.99e-09 (m^2/s) [20]. The height of the channel is just a fraction under the sum of the minor diameters of the two strands as shown in Figure 2. The domain was based on the findings of Li et al. [21] who researched on the most efficient dimensioning for zigzag spacer configurations. Furthermore, to determine the flow behaviour in SWM, narrow channels have been used to simulate the flow [14,22]. Ranade and Kumar [23] simulated curvilinear channels lined with spacers to determine the effect of the curvature if any on the fluid flow. Their work reiterated that the flow in straight narrow channels is not significantly different than that in curved channels and was reinforced by Schock and Miquel [24]. Use a two-column format, and set the spacing between the columns at 8 mm. Do not add any page numbers.

Mass Conservation Equation:
(1) Momentum Conservation Equation: (2) Concentration Equation: ( The channel height was taken to be the characteristic length and used for the calculation of Reynolds number (Figure 2.) Table 1 details the test cases run for the purpose of this simulation. Instead of calculating the hydraulic Reynolds number as detailed by Schock and Miquel [24] which is purposeful in the case of three dimensional flow analysis, Channel Reynolds Number is calculated according to the equation [4]: (4) Velocity inlet with a pressure outlet was used. No slip on the wall patches and zero gradient on the outlet were used as velocity boundary conditions.
Outlet was defined as zero atmospheric gauge pressure with Neumann type condition at the inlet and the walls. Similarly, concentration was defined to be 0.006 mol/m^3 at the and Neumann type at the walls and the outlet. Structured meshing through OpenFoam's blockMesh utility was done to entirely capture the curvature of the obstacles. 0.3 Million nodes were used for the simulations.
Open source CFD tool, OpenFOAM [25] was utilized for the simulations. The solver used was SimpleFoam, which is a Semi-Implicit pressure based steady state solver for incompressible flow which assumes constant density, hence, the pressures and shear stresses obtained are density normalized. The default solver is designed to only resolve pressure and velocity fields. Therefore, a new solver was designed which coupled the scalar transport equation with the mass and momentum conservation equations. The following few lines mention briefly the code behind the new solver: { fvScalarMatrix CEqn ( fvm::div(phi, C) -fvm::laplacian(DC, C) ); CEqn.relax(); CEqn.solve(); } Table 2 details the important parameters of the mesh of two configurations. Since the study was comparative in nature, it was imperative that the mesh parameters are controlled and kept similar to allow for valid conclusions to be drawn. Structured meshing was done for this reason. The CC achieved convergence under 6500 iterations and EC converged within 8000 iterations for each case detailed in Table 1.

Results and discussion
The flow over an ellipse in a narrow channel is profoundly different from flow over a circle. The inflection point, where the boundary layer separates from the wall of the channel in EC is noted to be 9.01% upstream of that of CC as observed by the difference in the first x-intercept of the two cases on Figure 3 which measures the shear stress at the bottom wall along Wall B in Figure 2. This difference in the separation points can also be visualized in the streamlines plotted in Figure 6. Even though the flow experiences much stronger adverse pressure gradient as shown by Figure 5, the effect of the higher-pressure gradient is overcome by the greater bulge of the ellipse. The separation bubble in EC is longer compared to that in CC upstream of the spacers. Downstream of the spacer, the flow reattaches earlier in EC and a difference 2.088% is recorded (Figure 3.) This is significant because the separation bubble this time is larger in CC than in EC. Unlike the upstream recirculation zone, this recirculation zone is larger and experiences far greater vorticity. Due to the high velocity of rotation of the fluid particles and the larger radius of motion, the mass transfer phenomena is very significant in this region. Therefore, a larger separation region signifies relatively higher mass transfer and hence an over estimation in flux levels.
Note that the pattern in the shear stresses experience in the two cases is very similar, yet different in values. The velocity and hence the shearing stresses are highest at the constricted section of the channel resulting in peak values of wall shear stresses.
Wall Shear Stresses are given significant importance in the study of reverse-osmosis because they contribute towards increase pressure drop and contrastingly towards reduced concentration polarization. These stresses have the potential to erode away any deposited salt on the membrane thereby increasing flux. Therefore, higher values of wall shear stresses directly contribute towards greater flux in CC. The results also agree with the findings of Fimbres et al. [12] who reported highest wall shear stresses just downstream of the reattachment points due to the converging channel.
Between the CC and EC, vorticity is observed to be higher in CC at almost every location along the membrane as shown by Figure 4. Vorticity is accredited with periodically disrupting the concentration boundary layer which prevents the buildup of concentration along the wall, combating concentration polarization. It can be clearly observed that the angular velocity is higher throughout the wall in CC which directly leads to the observation that the permeate flux would be much higher owing to the higher velocities. This is also reinforced by the wall shear stress plot which shows peaks in regions with higher vorticities.
The stream lines as shown by Figure 6 indicate the formation of strong recirculation zones upstream and downstream of the obstacles. The difference in the separation and reattachment points can be visually observed.  Figure 5 indicates how the density normalized pressure drop across the unit varies between the two cases at different Reynolds numbers. It also shows that the two curves closely follow each other and start to diverge at relatively higher Reynolds numbers. The pressure lines start to diverge after a Reynolds number of 150 is achieved and at 250, an error of 3.15% is observed. Note that this pressure drop is only for flow over one unit cell consisting of a pair of spacers while in actual cases a single desalination module would have millions of such units and the cumulative pressure drop across all the units can take it in the order of mega Pascals, hence, the pressure drop too cannot be accurately predicted through CC at higher Reynolds Numbers. This pressure drop also helps us acknowledge the large energy drainage caused by reverse-osmosis.

Fig. 7: Pressure drop across the unit for both CC and EC
There were also notable differences in the behavior of concentration between the two cases shown by Figure 8 for the line along B in Figure 2. Concentration was relatively lower in the wake region of CC along the wall in comparison to the EC. This results in higher local mass transfer along the wall. The lesser the concentration along the wall as is the case in CC, the greater the permeate flux through the wall and the lesser the chance of any flow back into the retentate channel.
It can be observed from mathematical relationships that trans-membrane pressures are a determinant in permeate flux. All parameters such as vorticity, shear stresses etc in essence alter the local pressure along the membrane wall which in turn determine the flux. According to the mathematical construct given by Kedem and katchalsky [26] and Merten [27], the permeate flux is given by:  Figure 5 clearly indicates that while the pattern for static pressure along the membrane wall for both CC and EC tends to follow a similar pattern, there is a discernible difference between the two with the minimum difference of 5.6% and a maximum difference of 19.3%. Although, a linear relationship is not followed between transmembrane pressures and permeate flux as indicated by Equation 5, it can be concluded from the large difference in recorded local pressures that the permeate flux in CC would be considerably greater.
It can also argued according to the results shown by Figure 3, Figure 4 and Figure 6 that although the separation bubble upstream of the spacer is larger in EC, any gain in mass transfer in EC due to this is counteracted by the largely slowly rotating fluid indicated by the low vorticity and low wall shear stresses in the upstream region. Similarly, the wake region is affected more in CC than in it is in EC. The concentration directly behind the spacer is 5.55% less in CC than in EC. This can be traced back to the higher vorticity in this region in the case of CC.
The concentration boundary layer, measured in millimeters from the bottom wall labeled B in Figure 2 follows a very similar trend in both the cases. Boundary layer height was obtained by evaluating where the concentration changes its value to 0.99 of the bulk concentration value. The results are plotted in Figure 9. A higher boundary layer thickness is recorded in CC with a maximum difference of around 20% which signifies that the region of recirculation is larger in CC and hence a larger amount of water can be channelled closer to the membrane via the vortex systems that are dominant in those regions. However, as the trend suggests in Figure 9, most of the regions in CC and EC offer close proximity to each other and the average boundary layer thickness is observed to be relatively similar.

Conclusion
Both upstream and downstream of the spacer, the CC is observed to enhance the mass transport considerably greater than an EC. This means that some of the previous research done in two dimensional flow will not always yield very accurate results as properties like wall shear stresses, separation bubble length and location of separation and reattachment point can vary up to 9% and local pressures along the membrane wall may even vary by 20% between the two configurations.
The increased vorticity in CC results in greater mixing between the bulk fluid and the fluid present near the membrane walls. The higher vorticity also has a stronger disruption effect on the boundary layer. Since a major use of the spacers is to disrupt the developed boundary layer, a stronger vortex system does this much more efficiently compared to a weaker one with lower rotational speeds. This will impact the overall Sherwood number and hence any conclusions drawn from studies using CC to simulate zigzag spacers.