A New Design of Tubular Ceramic Membrane Module for Oily Water Treatment: Multiphase Flow Behavior and Performance Evaluation

Petroleum has been extracted from oil reservoirs using different techniques. This activity is accompanied for a large amount of water and sometimes mixed with gas. This produced water has a high oil concentration and other toxic chemical compounds, thus, it must be treated to be reused or released to environment according to environmental protection regulations. Currently, ceramic membrane technology has been employed in the wastewater treatment, due to its high benefit–cost ratio. In this sense, this work aims to study the oil–water mixture separation process using a new configuration of tubular ceramic membrane module by computational fluid dynamic (ANSYS Fluent software). The proposed model is composed of mass and linear momentum conservation equations coupled to Darcy’s law and SST k-ω turbulence model. Results of the volumetric fraction, pressure, and velocity distribution of the oil and water phases are presented and discussed. The results indicated that the proposed model and new device both have great potential to be used on the water/oil separation process and that the transmembrane pressure remains constant in the axial direction and decreases radially through the membranes, indicating an efficient system that favors the transport of clean water and oil retention.


Introduction
Produced water is the effluent from the oil and gas industry emitted in large quantities when water from underground reservoirs is brought to the surface during activities at sea or on land. It is a large amount of residual water with different characteristics, which carries different contaminants, In the oily water separation process, the forces exerted by pressure, adsorption, and electric charge cause deformation of oil particles with small diameter in such a way that these particles cannot pass through the larger pores in the membrane. This phenomenon occurs because the organic phase (oily phase) in the liquid-liquid dispersion system is separated from the aqueous phase by means of lipophilic or hydrophilic properties of the porous membrane. Knowing that the two liquids are not miscible with each other and that there are differences in the affinity of the membrane for the different liquids, one of these liquids will form a pure liquid layer of a certain thickness on the membrane surface under hydraulic and external forces. Thus, the concentration of the other liquid forms a gradient with the pure liquid layer that stimulates their separation [13].
The concentration polarization phenomenon is established from the formation of an oil layer and other contaminants from the effluent parallel to the membrane surface [14,15]. The increase in the concentration of these retained materials favors a concentration gradient that causes additional resistance to mass transfer across the membrane [16]. In fact, this polarization by concentration of solute on the membrane surface is the primary reason for the decline in flow during the initial period of a membrane separation process. This can occur in conjunction with irreversible incrustations in the membrane and/or reversible formation of an encrustation layer [17].
Increasing the speed of the effluent flow is one of the strategies to enhance the mass transfer rate and decrease polarization by concentration. This effect was clearly observed by Cui et al. [18], when they applied ceramic microfiltration membranes combined with zeolite to treat water contaminated with oil. The authors observed that a higher tangential flow velocity (0.01 m/s) doubled the flow through the membrane without affecting the oil rejection coefficient. In addition, the membrane performance remained excellent, even with high oil concentration in the feed solution (500 mg/L).
Membrane encrustation generally occurs as a result of adsorption of components from the feed solution to the membrane surface [19]. This encrustation reduces the membrane separation performance, but the occurrence of this phenomenon can be mitigated by modifying the membrane surface. For example, Zhou et al. [20] modified commercial microfiltration membranes based on Al 2 O 3 by coating with nanometric-sized ZrO 2 , to reduce clogging of the membrane by oil droplets. The constant flow of the modified membrane was obtained in a very short time, maintaining 88% of the initial flow rate and oil rejection above 97.8%. In addition, the combination of this coating, with the appropriate backwash and transverse flow rate, helps to reduce the obstruction of the membrane by oil droplets.
Another effective method to delay polarization by concentration, reduce fouling and, consequently, increase the permeation flow is the use of pulsed flow and the promotion of turbulent flow [15,21]. Despite the problems with polarization by concentration and fouling, ceramic membranes are more resistant to contamination than organic ones, when it comes to petrochemical effluents [15]. All this potential of ceramic membranes justifies the intensification of scientific works that seek to enhance the treatment of water produced using this device.
Ceramic membranes have high filtration efficiency due to their excellent hydrophilic properties and chemical and hydrothermal stability [13]. Ceramic membranes are considered an economical and promising technology for water treatment and purification due to their ability to achieve almost 100% rejection of dissolved solids. This is due to the fact that the membranes are composed of components with hydrophobic properties (such as Silane agents), which favor only the transport of vapor, while retaining the unwanted solutes in the liquid phase on the feed side [22].
The resistance of the membrane itself is one of the factors that influence the flow performance through the membrane and, necessarily, must be considered, as it will always exist from the beginning to the end of the filtration process. This resistance is due to geometric factors, such as pore size and distribution, thickness and affinity of the membrane surface with the solvent [14,19].
Based on the above, the separation process by ceramic membranes can be understood as being dynamic for different characteristics/constitution of oily effluents and the membrane itself. In addition, the wide variation in several operational parameters results in a different filtration potential and must Membranes 2020, 10, 403 4 of 26 also suit the specific properties of each effluent and type of membrane in operation. Because the great importance, membrane separation processes have been studied on the experimental and theoretical approaches. Under the theoretical point of view, the application of computational fluid dynamics (CFD) techniques has advanced technologically and helped in the understanding and development of new membrane separation processes, but problems such as fouling, concentration polarization, and pore clogging, which reduce the filtration rate and favor the contamination of the permeate, remain challenges to be overcome [23][24][25][26][27][28][29].
In this context, the construction of new design of ceramic membranes, especially with adaptations of the mesh and pore arrangement emerges as one of the first steps for the development of new prototypes of practical and efficient application in the treatment of produced water. In addition to the development of new membrane models, the simulation of the filtration process is essential for establishing the operational characteristics, understanding and maximizing the filtration process. In this premise, it is worth highlighting the behavior of the pressure distribution, velocity and concentration field of the components in the separation module as well as the porosity and permeability in the membrane.
Therefore, in addition to the works reported in the literature, this research aims to study the separation process of produced water (oily water) using a tubular ceramic membrane module and the computational fluid dynamics (CFD) technique. The idea is to evaluation the separation performance and fluid flow behavior inside module and membranes. The following innovative aspect of this research can be highlighted: The study focuses in a new design of the membrane module and in the robustness of the proposed mathematical model that includes coupled phenomena of turbulent multiphase flow inside the module shell and membrane (wall thickness and tubular region), not yet evaluated in the literature.

Description of the Domain under Study
In this study, a shell-tube ceramic membrane filtration module (domain under study) was proposed, which consists of a main cylinder (hull) equipped with a tangential inlet and outlet and four internal cylinders corresponding to the membranes porous ceramics (membranes 1, 2, 3, and 4), as shown in Figure 1.
The separation module acts through nine structural domains, namely: The four volumes of fluids related to the four membranes, the four permeate (internal part of the membrane) and the cylindrical hull, comprised between the internal surfaces, which limit the separator and the outer surfaces of the membranes. The geometric dimensions of this module, used in numerical simulations, are described in Table 1.

Computational Mesh
For a numerical solution of the conservation equations that govern the domain under study, it is necessary to convert the continuous domain into a discrete domain. These control volumes are three-dimensional partitions of the geometric domain, composed of edges and nodal points, in which the discretized governing equations must be solved, to predict the physical phenomenon. The The module entrance is a tangential circular tube through which the contaminated effluent enters the system. The outlet tube, also tangential circular, is intended for the outlet of the concentrate. The filtrate is retained inside the membranes after filtration. The fact that the fluid flow in the membrane is considered not yet included in the literature for this type of separator makes the proposal of this research innovative. Thus, the separation between oil and water occurs as the feed mixture, containing a certain oil concentration, enters tangentially into the separator. Mainly due to the difference in transmembrane pressure, the mixture is forced to cross the membranes and, depending on the operational conditions and membrane characteristics, the oil fraction is retained in the porous structure, and only the aqueous phase easily crosses the membrane, generating the permeate. Finally, the four outlets coming from the membranes receive the permeate flow, while the outlet in the hull receives the concentrate (mixture of water with high oil concentration).

Computational Mesh
For a numerical solution of the conservation equations that govern the domain under study, it is necessary to convert the continuous domain into a discrete domain. These control volumes are three-dimensional partitions of the geometric domain, composed of edges and nodal points, in which the discretized governing equations must be solved, to predict the physical phenomenon. The numerical method causes errors of truncation and idealizations, which decrease, as the mesh build with high number of elements, that is, as the finite limits of solution are reduced in each control volume.
To adequately predict the proposed problem, three hybrid meshes (tetrahedral and hexahedral elements) with different densities of control volumes (elements) were developed, using the ANSYS Meshing ® 15.0 software (Canonsburg, PA, USA), as shown in Table 2. For the hull domain, hexahedral elements close to the outer walls of the membranes and the hull itself were built, and tetrahedral elements were built for the other regions of the study domain. In the membranes, structured meshes were built with only hexahedral elements, which was maintained for the cylindrical domains, referring to the permeate volumes. This procedure was realized using the mesh construction technique named o-grid.  Table 2 shows the parameters of mesh quality, orthogonality (minimum and average), deformation (maximum and average) and number of elements, referring to the three meshes, with 218,704, 508,325, and 853,536 elements, developed throughout the research. It is observed, from the ANSYS Meshing ® software, that all meshes are within the recommended limits for values of deformation (below 0.95) and orthogonality (above 0.1). These meshes were evaluated for the mesh dependence study in relation to the number of elements. Figures 2 and 3 illustrate aspects of mesh 2, with emphasis on the studied domains. Especially in Figure 2, it is possible to observe details in the mesh of the hull, membranes and permeate domains. There is also a cross-section of the entire module, focusing on the membrane mesh, as well as the mesh entrance and exit regions.      The Eulerian-Eulerian multiphase model is capable of modeling multiple phases, treated as separate, yet interactive. The phases can be liquid, gaseous, or solid, in almost any combination. A Eulerian-Eulerian solution treatment is used for each phase separately, even if one of the phases is made up of particles. The model makes no distinction between fluid-fluid and fluid-solid (granular) flow. A single pressure is shared by all phases, and the linear momentum and continuity equations are solved for each phase [30].
The formulation of the model assumes that two or more phases are continuous and interpenetrable. For each phase added to the mathematical model, a variable is introduced: The volumetric fraction (α q ), of the respective phase "q". The volumetric fractions represent the volume occupied by each phase, and the mass and linear momentum conservation laws are satisfied for each phase individually. In each control volume (V), the sum of the volumetric fractions of the phases is equal to 1 (one), as follows. where: Membranes 2020, 10, 403 The effective density of phase q is defined by: where ρ q is the physical density of phase q. Thus, there are three possibilities: α q = 1: Indicates that the volume is completely filled by the q phase. α q = 0: Indicates that the volume does not even have a region filled by the q phase. 0 < α q < 1: Indicates that the volume has an interface between the q phase and one or more phases. Based on the local value of α q , of the n phases existing in the physical process, the properties and the appropriate variables are weighted in each region of the multiphase flow. In this work, it was considered that the mass and linear momentum conservation equations are solved for each of the present phases (continuous and dispersed). For each of the situations, the following considerations were adopted: • Newtonian and incompressible fluid; • Flow in permanent and isothermal regime; • Mass transfer, interfacial momentum and mass sources were neglected; • The non-drag interfacial forces (lift forces, wall lubrication, virtual mass, turbulent dispersion, and solid pressure) were neglected; • There is no transfer of interfacial mass; • The module walls are static and null roughness; • The water stream is considered to be a mixture of water and oil; • The viscosity and density of the mixture and other physical and chemical properties are constant; • The porous medium (ceramic membrane) has an isotropic distribution of pores and permeability; and • There is no reaction or adsorption of the solute on the contact surface in the porous medium.
From the above assumptions, the mass conservation equation for multiphase flow was defined by Equation (4), as follows: where the sub-index "q" represents the phase involved in the two-phase water/oil mixture; α, ρ and → v are the volumetric fraction, density, and velocity vector, respectively.
In turn, the linear momentum conservation equation for multiphase flow was defined by Equation (5), given below: where = τ q is the stress-tensor for the q phase and → R pq is the term for the interface forces. This force depends on friction, pressure, cohesion and other effects and is subject to the following condition: The Fluent software solver uses a simplified interaction term, as follows: Membranes 2020, 10, 403 9 of 26 where K pq = K qp is the interfacial coefficient of momentum transfer; → v p and → v q are the velocities of the p and q phases.
For a two-phase flow, it is assumed that the secondary phase is in the form of drops. This has an impact on how each fluid is assigned to each phase, for example: In flows where there are unequal quantities of two fluids, the predominant fluid must be modeled as the primary fluid, since the sparse fluid is more likely to form droplets or bubbles [31]. The exchange coefficient for bubbling, liquid-liquid or gas-liquid mixtures can be written as follows: is the interfacial area, f is the drag function, which is defined according to the exchange rate model used, and the term τ p is the "particulate relaxation time", defined as follows: where d p is the diameter of the drop.
To determine the drag function f , the Schiller and Naumann model, was used, as follows: where C D , is the drag coefficient, given by: with Re being the relative Reynolds number, defined for the primary phase q and for the secondary phase p, as follows: To describe the turbulent flow within the system, the SST k-ω turbulence model was used. To activate the equations of the standard k-ω and k-ε turbulence models, coupling functions are used in cells near and far from the walls, respectively. The variable k represents the turbulent kinetic energy of flow, while the variable ω represents the energy dissipation rate. The transport of these variables is described by Equations (13) and (14). and where t is time, X is the position vector, u i is the velocity vector, sub-indices i and j represent the components (x, y, and z) of the coordinate axes, so that, if i or j = 1, one has the component in the x-direction; if i or j = 2, we have the component in the y-direction, if i or j = 3, we have the component in the z-direction. In addition, the term G k represents the production of turbulent kinetic energy, G ω represents the generation of ω, and Γ k and Γ ω are the effective diffusivity of k and ω, respectively, due to turbulence.
The variables Y k and Y ω , represent the dissipation of k and ω, the term D ω represents the cross diffusion, and the terms W k and W ω are the source terms of the referred equations.
The effective diffusivities of k and ω (Γ k and Γ ω ) are given by Equations (15) and (16), as follows: and where σ k and σ ω are the turbulent Prandtl numbers given by Equations (17) and (18), below: and In the Equations (17) and (18), F 1 is one of the coupling functions between the models k-ω and k-ε, described in Equation (19), and σ k,1 , σ k,2 , σ ω,1 , and σ ω,2 are constants of the model which are described in Table 3. Table 3. Turbulence model constants k-ω SST (shear-stress transport).

Constants
Numeric Value The F 1 coupling function is given by: where θ 1 is determined by Equation (20).
being d g the distance between the position studied in the domain and the nearest wall. In turn, D + ω , is calculated by Equation (21).

of 26
The term µ t in the Equations (15) and (16) is the turbulent viscosity given as follows: where α 1 = 0.31 [30] and α * is the damping factor relative to the turbulent viscosity, used to correct its value, for low Reynolds numbers, given by Equation (23).
In the Equation (23), Re t is the turbulent Reynolds number, which is calculated by Equation (24), and R k and α ∞ * are model constants defined in Table 1. The parameters Re t e α 0 * , can be obtained by Equations (24) and (25), as follows: and where, being β i,1 and β i,2 described in Table 1.
The term F 2 , from Equation (22), is another coupling function between the models k-ω and k-ε, which is defined by Equation (27). where, The term S of Equation (22) is the magnitude of the fluid deformation rate, which is given by Equation (29). where, The term G k of Equation (13) is defined by Equation (31) and represents the production of turbulent kinetic energy (k): where the term u i u j represents Reynolds stresses, derived from turbulent flow. The term G ω , from Equation (14), represents the production of the turbulent dissipation rate, and is given by: where, being R ω a model constant (Table 1), and α ∞ is defined by Equation (34).
In the Equation (34), α ∞,1 and α ∞,2 are given by: and where λ is a model constant ( Table 3). The term Y k , from Equation (13) represents the dissipation of k, which is defined by Equation (37). where, The term β i * of Equation (38) is given by: where ξ * , R β and β ∞ * are constants in the model ( Table 3).
The term Y ω in Equation (14) represents the dissipation of ω, which can be calculated by Equation (40).

Multiphase Model for Porous Media
In the Fluent software, the porous medium is modeled as a region containing elements of fluid, where the linear momentum equation is modified by the addition of a source term of dissipation. This source term is composed of two parts: One referring to the loss of viscous head (Darcy's law, the first term on the right side of Equation (41)) and one referring to the loss of inertial head (the second term on the right of Equation (41)).
where S i is the source term for the i-th (x, y or z) momentum equation, |v| is the magnitude of the velocity, and M and N are prescribed matrices. This momentum sink contributes to the pressure gradient in the porous cell, creating a pressure drop that is proportional to the fluid velocity in the cell. Finally, homogeneous porous media were defined by Equation (42).
where K i is the permeability, and C 2 is the inertial resistance factor, simplifying the matrices M and N as diagonal matrices with 1 K i and C 2 , respectively, occupying the values on the diagonals of the matrix.
Given the low velocities developed in the volumes relative to the porous medium, the term referring to inertial resistance was neglected in this research.

Boundary Conditions (a) Module Input
The prescribed mass flow condition was established at the input of the separator module. Specifying the mass flow rate allows the total pressure to vary in response to the numerical solution. In this boundary condition, the absolute reference system, the direction of flow (normal to the inlet surface), the turbulence intensity, I = 5%, and the turbulent viscosity ratio R µ = 10, given by Equations (43) and (44) were established: and where u is the rate of velocity fluctuation, and u is the mean velocity of free flow. The values of k and ε are computed as a function of these parameters [30].

(b) Concentrate and Permeate Outlets
The prescribed pressure boundary condition was applied to the permeate and concentrate outlets. There was a zero gauge pressure at the outlets, that is, the ambient pressure. The pressure difference between the input of the separation module and the outputs (atmospheric pressure) drives the flow of produced water through the separation module.

(c) Module Wall and Membrane Surface
The wall conditions were used to connect the fluid and solid regions, surfaces of the module and membranes, in contact, externally, with the volume of the hull and, internally, with the domains related to permeate. Boundary conditions of no-slip and negligible roughness were applied to the external surfaces of the hull, surfaces of the supports, and in the ends of the membranes, hull, and permeate. For the internal and external surfaces of the membranes, the condition of interior wall (open surface) was used, which allowed the flow of produced water through the membranes, where, due to the source term relative to the porous medium, the phases were separated. Table 4 describes the limit regions of the module under study and their respective boundary conditions.

Numerical Solution Method and Simulation Parameters
The discretization method applied in this work was that of finite volumes, the discretization method used by the software, which performs a balance of conservation equations for each control volume, which constitutes the numerical mesh.
For the pressure-velocity coupling, the Coupled algorithm, available in the ANSYS Fluent ® software, was used. The Coupled algorithm solves the continuity equations, based on the results of the calculation of the linear momentum conservation and the pressure, in a coupled way, which gives it, in relation to the segregated solution algorithms, a convergence range, with a smaller number of iterations.
The relaxation factors and the standard convergence criteria of ANSYS Fluent were adopted. The water and oil used in the numerical simulations had densities of 998.2 and 997 kg/m 3 , respectively, and dynamic viscosity of 0.001003 and 0.05 Pa·s, respectively. Based on the work of Cunha [26], the operating variables were standardized in mass flow rate ( . m) 0.5 kg/s, oil concentration at the inlet, (C 0 ), 1 kg/m 3 , mean diameter of oil droplet (d p ), 63 µm, membrane permeability (K), 3 × 10 −13 m 2 , and membrane porosity (ε), 30%.

Plans and Points Studied Inside the Module
To study the characteristics of the filtration process during CFD simulation, in the module's geometry was inserted different Cartesian planes, considering the x (transversal length), y (height), and z (longitudinal length) axis. Then, different planes were selected and evaluated at different points on the x, y, and z planes (Figures 4-6).

Numerical Solution Method and Simulation Parameters
The discretization method applied in this work was that of finite volumes, the discretization method used by the software, which performs a balance of conservation equations for each control volume, which constitutes the numerical mesh.
For the pressure-velocity coupling, the Coupled algorithm, available in the ANSYS Fluent ® software, was used. The Coupled algorithm solves the continuity equations, based on the results of the calculation of the linear momentum conservation and the pressure, in a coupled way, which gives it, in relation to the segregated solution algorithms, a convergence range, with a smaller number of iterations.
The relaxation factors and the standard convergence criteria of ANSYS Fluent were adopted. The water and oil used in the numerical simulations had densities of 998.2 and 997 kg/m 3 , respectively, and dynamic viscosity of 0.001003 and 0.05 Pa·s, respectively. Based on the work of Cunha [26], the operating variables were standardized in mass flow rate (m ) 0.5 kg/s, oil concentration at the inlet, (C ), 1 kg/m 3 , mean diameter of oil droplet (d ), 63 μm, membrane permeability (K), 3 × 10 −13 m 2 , and membrane porosity (ε), 30%.

Plans and Points Studied Inside the Module
To study the characteristics of the filtration process during CFD simulation, in the module's geometry was inserted different Cartesian planes, considering the x (transversal length), y (height), and z (longitudinal length) axis. Then, different planes were selected and evaluated at different points on the x, y, and z planes (Figures 4-6).  The relative volumetric fraction of oil and the pressure were evaluated in the hull regions in the xy planes at z = 0, 30, 75, 120, and 150 mm; in the yz planes, at x = −15, 0, and 15 mm (where x = 0 mm corresponds to the center of the domain) and in the xz planes, at y = −15, 0, and 15 mm (where y = 0 mm corresponds to the center of the domain). The volumetric fraction of oil was also evaluated in the porous medium in the xy planes at z = 30, 75, and 120 mm, and throughout the permeate collecting tube. The mixing fluid velocity was evaluated throughout the porous medium and the permeate collecting tube.
The pressure was also evaluated inside the membranes, considering horizontal lines (x-axis): One line on the right (x ranging from −0.025 m to −0.020 m) and one on the left (x ranging from −0.005 m to −0.010 m), in such a way that the pressure was measured from the external surface (water/oil mixture side) to the inner membrane surface (permeate side); and vertical lines (y-axis): One line in the upper position (y ranging from 0.025 m to 0.020 m) and another in the lower position (y ranging from 0.005 m to 0.010 m). Both lines were evaluated at positions z = 30, 75 and 120 mm from the origin of the membrane (direction of the z axis).       Considering the hull region as a whole, it is observed that the oil volumetric fraction varied from 89.2% to 208.4%. Exactly at the origin of the hull (considering the xy plane when z = 0 mm), it is possible to observe a high volumetric fraction of oil (ranging from 154.4% to 174%) at the location of the module's feed tube (Figure 7b), indicating that the oil is present in the water/oil mixture supplied to the system. Still in the plane of origin of the hull, it is possible to observe a lower oil concentration on the outer sides of the hull. To give you an idea, a volumetric fraction of 1% corresponds to an oil concentration of 9.97 kg/m 3 or 9.97 g/L or 9970 ppm.

Relative Volume Fraction of Oil in the Filtration Module
At 30 mm from the origin of the module (z = 30 mm), an oil volumetric fraction ranging from 98.9 to 209.8% is observed, more evenly distributed in relation to the plane at z = 0 mm and being more concentrated in the region between membranes and close to their surface, as shown in Figure 7c. From 30 mm, along the z axis, the maximum volume fraction of oil gradually reduced to the planes z = 75 mm (198.8%), z = 120 mm (191.1%,) and z = 150 mm (111.2%), indicating that the presence of oil in the hull is reduced, as the fluid goes towards the exit of the module (Figure 7d-f).  Figure 8 shows the distribution of the oil volume fraction at different planes. From the analysis of this figure, it can be seen that the average volumetric fractions of the oil in the parallel planes at different points along the x and y axes varied from 92.4% to 167.9% and from 91.3 to 170.2%, respectively (Figure 8a,b, respectively). Further, it can be seen that the lateral planes in the positions x = −15 mm and 15 mm (Figure 8c,g) and y = −15 mm and 15 mm (Figure 8d,h) present the largest  Figure 8 shows the distribution of the oil volume fraction at different planes. From the analysis of this figure, it can be seen that the average volumetric fractions of the oil in the parallel planes at different points along the x and y axes varied from 92.4% to 167.9% and from 91.3 to 170.2%, respectively (Figure 8a,b, respectively). Further, it can be seen that the lateral planes in the positions x = −15 mm and 15 mm (Figure 8c,g) and y = −15 mm and 15 mm (Figure 8d,h) present the largest maximum oil volume fractions when compared with their respective central planes, x = 0 mm and y = 0 mm (Figure 8e-f).
According to Habert et al. [32], Cunha [26], and Souza [27] this behavior occurs, because the oily fraction is retained mainly on the external surface of the porous medium, as the fluid flow enters the filtration module. Additional resistance to mass transfer occurs due to the establishment of a concentration gradient, leading to a decrease in the permeate flow. Therefore, this is a phenomenon that must be controlled and minimized, as it reduces the permeate flow and can affect the quality of the product [26]. When evaluating the relative oil volumetric fraction in the four porous membranes that make up the filtration module (Figure 9), a slight variation was observed in the maximum volumetric fractions of oil between the membranes as well as in the regions of greater oil presence in each membrane. For example, in membrane 1 (Figure 9a), the relative volumetric fraction of oil varied When evaluating the relative oil volumetric fraction in the four porous membranes that make up the filtration module (Figure 9), a slight variation was observed in the maximum volumetric fractions of oil between the membranes as well as in the regions of greater oil presence in each membrane. For example, in membrane 1 (Figure 9a), the relative volumetric fraction of oil varied from 57.9% to 185.7%, while in membrane 3 (Figure 9c), the oil volumetric fraction varied from 47.3 to 208.7%. In addition, the region with the highest relative volumetric fraction of oil in membrane 1 was close to the exit region of the module, while, in membrane 3, this same region showed the lowest oil concentration. However, the relative oil volumetric fraction field was slightly similar between membranes 2 and 4 ( Figure 9b,d), which, in turn, differed slightly from that observed in membranes 1 and 3. In general, the terminal region (close to the module exit) of membranes 2, 3, and 4 showed lower relative volume volumetric fractions of oil in relation to the region comprising 2/3 of the membrane measured from z = 0 mm.
When analyzing the iso-surfaces of the relative oil volumetric fraction in the permeate (Figure 10), there is an absence of oil in the innermost region of the permeate (permeate outlet) or very low volumetric fraction of oil (α/α 0 = 0.25%), located mainly in the part of the permeate collecting tube, which is close to the feeding region (Figure 10a). On the other hand, in the permeate near the membrane interface, a higher oil concentration is observed (α/α 0 = 75%), located mainly in the permeate outlet ( Figure 10f). In this sense, the fraction of oil retained gradually increased to more distant regions of the permeate collecting tube.
In general, the accumulation of relative oil fractions in the hull, membranes and permeate occurred due to the drag force and the tangential flow direction, which carries the retained particles to the terminal end of the module.
Membranes 2020, 10, x FOR PEER REVIEW 19 of 27 which is close to the feeding region ( Figure 10a). On the other hand, in the permeate near the membrane interface, a higher oil concentration is observed ( / = 75%), located mainly in the permeate outlet (Figure 10f). In this sense, the fraction of oil retained gradually increased to more distant regions of the permeate collecting tube.
In general, the accumulation of relative oil fractions in the hull, membranes and permeate occurred due to the drag force and the tangential flow direction, which carries the retained particles to the terminal end of the module.      Figure 11 shows the pressure distribution in the hull region. From the analysis of this figure, it can be seen that the pressure in the hull was higher in the region of entry and close to this region, reaching a maximum average pressure of 39.08 kPa. On the other hand, at the permeate outlet, the lowest pressure values were observed. This information is ratified when comparing the maximum pressure in the xy plane at z = 0 mm (inlet; Figure 11b) with z = 150 mm (inlet; Figure 11f). Comparing the results obtained in the different planes along the x-axis, at positions z = 0, 30, 75, 120, and 150 mm, a reduction in system pressure is observed that occurs gradually from the entrance to the module outlet.

Pressure in the Filtration Module
Observing the transverse planes at positions z = 0, 30, 75, 120, and 150 mm of the hull, it is noted that the pressure tends to be greater in the lower region of the hull, especially in the region below and between membranes 3 and 4. In the transverse planes in the positions z = 0 mm and z = 30 mm, higher pressures (36.93 and 34.80 kPa, respectively) are observed, located in the lower and lateral regions, while the lowest pressures occur close to the membranes (27.42 and 30.31 kPa, respectively). In the transverse planes in the positions z = 75 mm and z = 120 mm, the minimum and maximum pressure varied from 31.09 to 33.91 kPa and from 31.53 to 33.82 kPa, respectively. In both planes, higher and more distributed pressures can be observed in the lower region of the hull, especially in the region close to the surface of the membranes. In the transverse plane in the position z = 150 mm, the pressure varied from 32.55 to 33.87 kPa, with the highest pressure being located in the lower region, and the lowest pressure, located in the permeate exit region.
The pressure along the membrane was assessed in two horizontal lines (x-axis): One line on the right (x ranging from −0.025 to −0.020 m) and another line on the left (x ranging from −0.005 to −0.010 m), in such a way that this parameter was measured from the outer surface (water/oil mixture side) to the inner membrane surface (permeate side) and evaluated at z = 30, 75, and 120 mm, from the origin of the membrane (direction of the z axis), as shown in Figure 12a. Further, transmembrane pressure was also assessed on two vertical lines (y-axis): A line in the upper position (y ranging from 0.020 to 0.025 m) and another in the lower position (y ranging from 0.005 to 0.010 m), also at z = 30, 75, and 120 mm, from the origin of the membrane (direction of the z axis), as illustrated in Figure 12b.
lowest pressure values were observed. This information is ratified when comparing the maximum pressure in the xy plane at z = 0 mm (inlet; Figure 11b) with z = 150 mm (inlet; Figure 11f). Comparing the results obtained in the different planes along the x-axis, at positions z = 0, 30, 75, 120, and 150 mm, a reduction in system pressure is observed that occurs gradually from the entrance to the module outlet. Observing the transverse planes at positions z = 0, 30, 75, 120, and 150 mm of the hull, it is noted that the pressure tends to be greater in the lower region of the hull, especially in the region below and between membranes 3 and 4. In the transverse planes in the positions z = 0 mm and z = 30 mm, higher pressures (36.93 and 34.80 kPa, respectively) are observed, located in the lower and lateral The pressure in the membrane was practically constant along the z axis (longitudinal length) of the membrane, observed by the similar values in the positions z = 30, 75 and 120 mm of the entrance, and in all the right and left lines (direction x) and top and bottom (y direction). This behavior indicates an efficient system to maintain the same pressure throughout the filtration membrane module. In turn, the pressure differential in the system is a determining factor, to ensure the flow through the membrane, keeping the high filtration process efficiency for a longer period of time [8,11].
origin of the membrane (direction of the z axis), as shown in Figure 12a. Further, transmembrane pressure was also assessed on two vertical lines (y-axis): A line in the upper position (y ranging from 0.020 to 0.025 m) and another in the lower position (y ranging from 0.005 to 0.010 m), also at z = 30, 75, and 120 mm, from the origin of the membrane (direction of the z axis), as illustrated in Figure 12b.   Figure 13 illustrates the velocity distribution of the mixture at different planes within the membrane. Analyzing this figure, it can be seen that the velocity at the surface of the membranes increased in the direction of entry-exit, that is, the farther from the region of entry of the module, the greater the fluid mixture velocity on the membrane surface. In fact, the velocity at the initial position of the membrane surface was zero or very close to zero. This same behavior was observed in the 4 membranes inside the module, however the maximum velocity in membranes 1 and 2 (Figure 13a,b) was higher than the maximum velocity in membranes 3 and 4 (Figure 13c,d).

Mixing Speed in the Filtration Module
The increase in velocity along the membrane surface occurs because the tangential flow runs through the free space inside the module parallel to the membranes surfaces and in the form of a vortex. Thus, at various points along the z-axis, the behavior of the fluids present in the module causes an increase in the fluid mixture velocity that runs parallel to the membrane surfaces. These results are in line with those reported by Cipollina et al. [33], when simulating the flow and temperature fields in a region of a membrane module in spirals, with double layer filament spacer, using CFD tools and considering a 25 × 25 × 3 mm size domain in x, y and z directions, respectively. membrane. Analyzing this figure, it can be seen that the velocity at the surface of the membranes increased in the direction of entry-exit, that is, the farther from the region of entry of the module, the greater the fluid mixture velocity on the membrane surface. In fact, the velocity at the initial position of the membrane surface was zero or very close to zero. This same behavior was observed in the 4 membranes inside the module, however the maximum velocity in membranes 1 and 2 (Figure 13a,b) was higher than the maximum velocity in membranes 3 and 4 (Figure 13c,d). The increase in velocity along the membrane surface occurs because the tangential flow runs through the free space inside the module parallel to the membranes surfaces and in the form of a vortex. Thus, at various points along the z-axis, the behavior of the fluids present in the module causes an increase in the fluid mixture velocity that runs parallel to the membrane surfaces. These results are in line with those reported by Cipollina et al. [33], when simulating the flow and temperature fields in a region of a membrane module in spirals, with double layer filament spacer, using CFD tools and considering a 25 × 25 × 3 mm size domain in x, y and z directions, respectively. The permeate velocity inside the membranes (permeate collecting tube) gradually increased from the entrance region to the exit area of the system (Figure 14). This behavior was common throughout the flow and occurs because, along the longitudinal axis of the membrane (z-axis), the flow is continuously filtered, crossing the membrane and feeding the permeate tube, thereby increasing the mass flow rate of the permeate. Thus, the more mass flow rate of the permeate is added to the collection tube, an increase in the total flow velocity is expected.
Although membrane 1 showed a high maximum velocity on the surface, especially when compared to membranes 3 and 4, the maximum permeate velocity for membrane 1 (Figure 14a) was at least 4 times lower than that observed for membranes 3 ( Figure 14c) and 4 ( Figure 14d). The decrease in the fluid flow through the membranes during the filtration process occurred due to the obstruction of the membrane pores [26,27,32,34] as well as the formation of a concentration polarization layer, which is established on the surface of membranes and which can reduce the permeate volume per unit time [15].
The fluid mixture velocity in the membrane collecting tube (Figure 15) was evaluated in horizontal lines (x-axis, ranging from −0.020 m to −0.010 m) and vertical lines (y-axis, ranging from 0.010 m to 0.020 m) at the positions z = 30, 75, and 120 mm, from the origin of the module (direction of the z axis), as shown in Figure 15a,b. From the analysis of this figure, it appears that the fluid velocity was low at the position z = 30 mm, from the origin of the membrane, increasing in the direction of the z axis, as observed for the positions of z = 75 mm and z = 120 mm. These changes in the fluid mixture velocity inside the tube were very similar between the horizontal and vertical lines. compared to membranes 3 and 4, the maximum permeate velocity for membrane 1 (Figure 14a) was at least 4 times lower than that observed for membranes 3 ( Figure 14c) and 4 (Figure 14d). The decrease in the fluid flow through the membranes during the filtration process occurred due to the obstruction of the membrane pores [26,27,32,34] as well as the formation of a concentration polarization layer, which is established on the surface of membranes and which can reduce the permeate volume per unit time [15]. The fluid mixture velocity in the membrane collecting tube (Figure 15) was evaluated in horizontal lines (x-axis, ranging from −0.020 m to −0.010 m) and vertical lines (y-axis, ranging from 0.010 m to 0.020 m) at the positions z = 30, 75, and 120 mm, from the origin of the module (direction of the z axis), as shown in Figure 15a,b. From the analysis of this figure, it appears that the fluid velocity was low at the position z = 30 mm, from the origin of the membrane, increasing in the direction of the z axis, as observed for the positions of z = 75 mm and z = 120 mm. These changes in the fluid mixture velocity inside the tube were very similar between the horizontal and vertical lines. In the center of the permeate collecting tube (centroid of the x and y axes, −0.015 and 0.015 m, respectively), the permeate velocity was higher and decreased as it approached the inner surface of the membrane (region of membrane-permeate contact). It can be seen that the velocity profile is developing, from the behavior of a laminar flow (parabolic velocity profile) to a turbulent flow (flatter velocity profile).
In the center of the permeate collecting tube (centroid of the x and y axes, −0.015 and 0.015 m, respectively), the permeate velocity was higher and decreased as it approached the inner surface of the membrane (region of membrane-permeate contact). It can be seen that the velocity profile is developing, from the behavior of a laminar flow (parabolic velocity profile) to a turbulent flow (flatter velocity profile).

Conclusions
Currently, ceramic membrane has been considered one of the most promising technologies in the treatment of contaminated water. This is due to its excellent mechanical, chemical, and biological characteristics. In this research, simulations related to the process of separating a mixture of water and oil using a ceramic membrane module and the CFD technique were carried out under different operating conditions. From the results obtained, it can be concluded that: (a) The proposed mathematical modeling proved be effective in predicting the fluid dynamic behavior of the oil and water phases inside the separation module in the separation process; (b) the new configuration of the ceramic membrane module showed major performance compared to one ceramic membrane alone operating in the same conditions;

Conclusions
Currently, ceramic membrane has been considered one of the most promising technologies in the treatment of contaminated water. This is due to its excellent mechanical, chemical, and biological characteristics. In this research, simulations related to the process of separating a mixture of water and oil using a ceramic membrane module and the CFD technique were carried out under different operating conditions. From the results obtained, it can be concluded that: (a) The proposed mathematical modeling proved be effective in predicting the fluid dynamic behavior of the oil and water phases inside the separation module in the separation process; (b) the new configuration of the ceramic membrane module showed major performance compared to one ceramic membrane alone operating in the same conditions; (c) the volumetric fraction of oil increases from the entrance of the module to the exit of the concentrate and decreases inside the membrane, from the external surface to the exit of the permeate; (d) the pressure inside the module decreases from the feed inlet to the outlets of the concentrate and permeate, being higher in the vicinity of the lower inner surface of the hull; and (e) the transmembrane pressure was almost constant along the longitudinal length of the membranes, indicating an efficient system in maintaining the same pressure throughout the filtration module, favoring the mass transfer and oil retention.
Under the practical point of view, the contribution of this research is clearly seen: The proposed membrane module can be used in the analysis of different materials (ceramic or polymeric membranes, oil and water or solid particle and water, for examples); at different operating conditions of feed flow rate (low or high), fluid mixture composition (high or low solute concentration), membrane physical parameters (porosity and permeability), membrane geometrical parameters (length, diameter, and wall thickness), and at different liquid-liquid and liquid-solid separation processes (microfiltration, ultrafiltration, nanofiltration, and reverse osmosis). This way, it is possible to analyze different physical situations and to assist industrial engineers in making appropriate and safety decisions on this topic, according to the level of filter required and type of applications.