Non-Isothermal Treatment of Oily Waters Using Ceramic Membrane: A Numerical Investigation

: Currently, the oil industry deals with the challenge of produced-water proper disposal, and the membrane-separation technology appears as an important tool on the treatment of these waters. In this sense, this work developed a mathematical model for simulating the oil / water separation by a ceramic membrane. The aim was to investigate the thermal aspects of the separation process via computational ﬂuid dynamic, using the Ansys CFX ® 15 software (15, Ansys, Inc., Canonsburg, PA, USA). Oil concentration, pressure, and velocity distributions, as well as permeation velocity, are presented and analyzed. It was veriﬁed that the mathematical model was capable of accurately representing the studied phenomena and that temperature strongly inﬂuences the ﬂow behavior.


Introduction
In the oil fields, significant amounts of oil-contaminated water are generated which are commonly known as produced water. These waters are provided from the reservoir and water injection, during the secondary recovery process. The oily water is an aqueous solution that is rich in salts, which may contain dispersed or emulsified oil droplets.
Operations at refineries generate a large amount of difficult-to-treat wastewater, rich in organic pollutants. Effluent treatment plants are generally used to treat these waters. The purpose is to change the effluent physical and/or chemical compositions, to enable the water's proper disposal into the environment, according to the actual legislation. During the treatment, dispersed or dissolved hydrocarbons, suspended solids, toxic organic compounds (BTEX-benzene, toluene, ethylbenzene, and xylenes), phenol, and inorganic compounds-such as cyanides, sulfides, metals, phosphates, and nitrogen compounds-are removed [1,2].
Currently, many processes are used to remove contaminants from produced water, such as electrochemical treatment, membrane filtration, biological treatment, flotation, and photocatalysis, among others.
The electrochemical treatment uses an electrical potential to oxidize or reduce substances present in contaminated water; this treatment is considered to be a highly efficient, versatile method with easy automation and control [3]. However, the electrodes are consumed over time and need to be replaced regularly. process via CFD, as the thermal effects have not been well explored in the literature, presenting a new approach to treat oily water and providing data that enable the optimization of the currently used processes. Herein, it is realized a comparison between the two-and three-dimensional approaches of the process, aiming to contribute to the understanding of temperature and geometry influence in the concentration distribution and permeation velocity.

2D Analysis
For the studying of the oil-water separation process, in a two-dimensional approach, we used a tubular membrane with a tube diameter of 0.03 m and a length of 3 m [23], subject to oily water flow ( Figure 1). Due to the membrane angular symmetry, only one longitudinal section was used (YZ plane).
Energies 2020, 13, x FOR PEER REVIEW  3 of 19 literature, presenting a new approach to treat oily water and providing data that enable the optimization of the currently used processes. Herein, it is realized a comparison between the twoand three-dimensional approaches of the process, aiming to contribute to the understanding of temperature and geometry influence in the concentration distribution and permeation velocity.

2D Analysis
For the studying of the oil-water separation process, in a two-dimensional approach, we used a tubular membrane with a tube diameter of 0.03 m and a length of 3 m [23], subject to oily water flow ( Figure 1). Due to the membrane angular symmetry, only one longitudinal section was used (YZ plane). For the 2D numerical simulation, a mesh was generated, using ICEM CFD 15 software (15, Ansys, Inc., Canonsburg, PA, USA), having a total of 205,056 elements and 180,000 nodes, with inflations near the membrane surface, as illustrated in Figure 2.

Three-Dimensional Analysis
In this case, we used the water/oil separation module known as microfiltration, which consists of a tube shell, with diameter De, enclosing a tubular membrane with diameter di. The module has a rectangular input and output ( Figure 3). The main dimensions of the module are shown in Table 1.  For the 2D numerical simulation, a mesh was generated, using ICEM CFD 15 software (15, Ansys, Inc., Canonsburg, PA, USA), having a total of 205,056 elements and 180,000 nodes, with inflations near the membrane surface, as illustrated in Figure 2.
Energies 2020, 13, x FOR PEER REVIEW  3 of 19 literature, presenting a new approach to treat oily water and providing data that enable the optimization of the currently used processes. Herein, it is realized a comparison between the twoand three-dimensional approaches of the process, aiming to contribute to the understanding of temperature and geometry influence in the concentration distribution and permeation velocity.

2D Analysis
For the studying of the oil-water separation process, in a two-dimensional approach, we used a tubular membrane with a tube diameter of 0.03 m and a length of 3 m [23], subject to oily water flow ( Figure 1). Due to the membrane angular symmetry, only one longitudinal section was used (YZ plane). For the 2D numerical simulation, a mesh was generated, using ICEM CFD 15 software (15, Ansys, Inc., Canonsburg, PA, USA), having a total of 205,056 elements and 180,000 nodes, with inflations near the membrane surface, as illustrated in Figure 2.

Three-Dimensional Analysis
In this case, we used the water/oil separation module known as microfiltration, which consists of a tube shell, with diameter De, enclosing a tubular membrane with diameter di. The module has a rectangular input and output ( Figure 3). The main dimensions of the module are shown in Table 1.

Three-Dimensional Analysis
In this case, we used the water/oil separation module known as microfiltration, which consists of a tube shell, with diameter D e , enclosing a tubular membrane with diameter d i . The module has a rectangular input and output ( Figure 3). The main dimensions of the module are shown in Table 1. For the 3D mesh generation, we used the blocking strategy, with the creation of hexahedral elements within the domain, aiming to facilitate the refinement. The resulting final mesh has a total of 217,352 elements and 180,780 nodes, as illustrated in Figure 4.  For the 3D mesh generation, we used the blocking strategy, with the creation of hexahedral elements within the domain, aiming to facilitate the refinement. The resulting final mesh has a total of 217,352 elements and 180,780 nodes, as illustrated in Figure 4. For the 3D mesh generation, we used the blocking strategy, with the creation of hexahedral elements within the domain, aiming to facilitate the refinement. The resulting final mesh has a total of 217,352 elements and 180,780 nodes, as illustrated in Figure 4.

Mathematical Modeling
The water/oil separation process through the tubular membrane was described by using mathematical modeling based on the conservation equations (mass, momentum, and energy) and mass transport. The aim was to evaluate the thermal effect on the fluid-flow behavior.

The Model
For the study of the microfiltration process using a 2D tubular membrane and a shell and tube membrane-separation module, the following considerations were used for the mathematical modeling: (a) Steady-state; (b) Laminar and incompressible flow; (c) Fluids density and viscosity are dependent on the temperature; (d) The local permeation rate was established by the series resistances theory; (e) Concentration layer is considered homogeneous, and the Carman-Kozeny equation is valid; (f) The axial velocity profile is fully developed in the membrane inlet; (g) The obstruction of the porous medium by the oil is neglected; (h) The mass diffusivity of the dispersed phase (oil) was assumed to be constant; (i) The gravitational effect was neglected; (j) The fouling mechanisms were neglected; (k) The temperature range from 35 to 95 • C was used for the analyzes.
From the considerations made, the following equations can be written: (a) Mass conservation equation where ρ is the density, and U is the velocity vector. (b) Linear momentum equation where µ is the solution viscosity, and p is the pressure. (c) Energy conservation equation where H is the enthalpy, and Γ e is the effective thermal diffusivity. (d) Mass transport equation D AB is the mass diffusivity (constant for each Schmidt number) (Equation (5)), and C is the oil concentration.

Boundary Conditions
(a) Inlet condition (2D and 3D analyses) At the equipment inlet is considered fully developed flow, at an initial temperature (T 0 ), with zero radial velocity and axial velocity expressed by Equation (6): Re µ ρ(D e − d i ) (6) where Re is the axial Reynolds number. Oil initial concentration (Equation (7)) and the inlet fluid temperature (Equation (8)) are expressed as follows: (b) Output condition (2D and 3D analyses) At the system outlet was assumed a pressure condition equal to atmospheric (1 atm) and parabolic condition for concentration (Equations (9) and (10)), as follows: (c) Symmetry condition (only for the 2D analysis) The following symmetry conditions were assumed on the tube axis: The following conditions were assumed for the transverse planes:

(d) Porous wall condition (2D and 3D analyses)
A no-slip condition (zero axial velocity) (Equation (17)) was assumed in the membrane surface, as well as an equation responsible for convective heat transport (Equation (18)).
where T ∞ is the fluid temperature (air), T s is the surface membrane temperature, and h c is the convective heat transfer coefficient (W/m 2 K), calculated according to Equation (19).
where k is the thermal conductivity, Nu is the Nusselt number, and L c is the characteristic length. Characteristic length, for flow over horizontal cylinders, is defined as the tube diameter (Equation (20)): To calculate the Nusselt number, for natural convection on a horizontal cylinder, Equation (21) was used [28].
where Pr is Prandtl number, Ra is the Rayleigh number (Equation (22)), and Gr is the Grashof number (Equation (23)).
where µ is the dynamic viscosity, and β is the thermal expansion coefficient. Through the membrane, there is a radial velocity, U y , expressed by the permeation velocity,U w , as follows: where R p is the resistance due to concentration polarization, and R m is the membrane resistance. Mass transport (Equation (25)) [27] was added, to be modeled as a source term. In this case, the following equation was used: where R r is the oil intrinsic retention by the membrane, considered constant and equal to one. This consideration means 100% of oil rejection by the membrane. In Equation (24), the transmembrane pressure, ∆P, can be express as follows: (26) where P p is the average permeate pressure, and P ex is the external pressure to the membrane. The terms R m and R p in Equation (24) are expressed by Equations (27) and (28), respectively.
where e is the membrane thickness, K is the porous media permeability, δ p is the concentration layer thickness, and r p is the resistance (when the concentration layer is considered homogeneous). The parameter r p is calculated by using the Carman-Kozeny equation, as follows: where ε p is the porosity of the concentration polarization, and d p is the mean particle diameter (oil).
where z is the axial coordinate along the membrane, Sc is the Schmidt number, d i is the internal diameter, Re is the axial Reynolds number, and Re w is the permeation Reynolds number. According to the authors [24], this correlation is especially useful, as it allows for the evaluation of the boundary-layer thickness by polarization with a 10% average error, as a function of the operational parameters.

Physical Properties of the Fluid
Thermal conductivity, specific heat, density, and viscosity of the water were calculated by using expressions as a function of the temperature, as reported in the literature.
(a) Thermal conductivity (k) The thermal conductivity of water as a function of temperature is given by Equation (31) [29]: where T is the temperature in Kelvin.

(b) Specific heat (c p )
Specific heat of the water was determined by Equation (32), as follows [30]: where T is the temperature in Kelvin, M w is the molecular mass of water in kg/kmol, and A, B, C, and D are the regression coefficients as reported in Table 2.
Water density as a function of the temperature was determined by using the following equation [25]: where T c is the critical temperature, and the parameters A, B, and N are the regression coefficients, as reported in Table 3.  [30].
where T is the temperature in Kelvin, and the terms A, B, C, and D are the regression coefficients, as shown in Table 4.
The model validation was realized in a previous work [9], by comparing the numerical results obtained in the simulations with the data available in some works reported in the literature [23][24][25][26]31].

Studied Cases
Simulations were performed, by using different fluid inlet temperatures, and the other parameters were kept constant; these values were adopted based on the work developed by Magalhães et al. [9], Damak et al. [23], and Magalhães et al. [32]. Table 5 presents information about the parameters used in the simulations, and Table 6 shows all the cases studied.  Table 6. Conditions used in the simulations.

Cases T 0 ( • C) Membrane Location
Two-dimensional analysis  Figure 5 shows the concentration profiles along the z coordinate (for y = 0.0225 m), for Cases 2 (two-dimensional) and 9 (three-dimensional), both with T 0 = 55 • C. From the analysis of this figure, it is observed that the oil concentration increases along with the axial position, due to the particle accumulation near the membrane surface. Furthermore, it is observed that the geometry of the device considered in the analysis slightly affect the concentration profile, since the difference between the concentrations of the 3D and 2D analyses, along with the radial position, is only 3 × 10 −3 kg/m 3 . When comparing the two-dimensional result with that obtained in a three-dimensional approach, where the permeation occurs on the outer tube surface of the shell and tube membrane separation module, it is verified that the 3D geometry presents higher concentrations, with variations close to the inlet (z = 0 m) and outlet (z = 3 m), presenting inlet and outlet solute concentration equal to 1.000 and 1.0089 kg/m 3 , respectively; this behavior is induced justly by feed inlet and outlet locations of the concentrate in the tridimensional geometry, which is not observed for 2D geometry.  Figure 6 shows the oil-concentration behavior as a function of the longitudinal position over four angles (0°, 90°, 180°, and 270°) on the membrane surface of the module, Cases 5 (T0 = 35 °C) and 8 (T0 = 95 °C). It illustrates that the oil-concentration distribution is affected by the variation temperature, especially at the outlet region (z = 3 m). In Figure 6a, it can be seen that the concentration profile, at 0° and 180° angular positions, have higher oil concentrations, as compared with the results observed in Figure 6b, for the angular positions 90° and 270°, due to the fluid flow behavior around the membrane. Higher concentrations were obtained, as expected, at the impact region between the fluid and membrane, provided by the feed inlet perpendicular position.  angles (0°, 90°, 180°, and 270°) on the membrane surface of the module, Cases 5 (T0 = 35 °C) and 8 (T0 = 95 °C). It illustrates that the oil-concentration distribution is affected by the variation temperature, especially at the outlet region (z = 3 m). In Figure 6a, it can be seen that the concentration profile, at 0° and 180° angular positions, have higher oil concentrations, as compared with the results observed in Figure 6b, for the angular positions 90° and 270°, due to the fluid flow behavior around the membrane. Higher concentrations were obtained, as expected, at the impact region between the fluid and membrane, provided by the feed inlet perpendicular position.  Figure 7 shows the permeation velocity profile as a function of the axial position for Cases 1-4 (two-dimensional cases) and Cases 5-8 (three-dimensional cases). It is verified that, as the temperature is increased, a reduction in the permeation velocity is observed; this behavior is associated with the reduction of the pressure gradient within the equipment, due to the inlet velocity reduction, caused by the variation in the fluid properties. For example, by raising the temperature from 35 to 95 • C, the relationship between viscosity and density (µ/ρ) is reduced by 4.14 × 10 −7 m 2 /s. Furthermore, it was verified that the geometry influences the permeation velocity; the characteristic curves for the three-dimensional cases present lower velocities than the two-dimensional curves. This behavior is associated with the enlargement of the runoff area under the same operating conditions imposed in the 2D analysis. Reinforcing this analysis, it is possible to visualize the temperature influence on permeate mass flow (Table 7). from 35 to 95 °C, the relationship between viscosity and density (μ/ρ) is reduced by 4.14 × 10 m /s. Furthermore, it was verified that the geometry influences the permeation velocity; the characteristic curves for the three-dimensional cases present lower velocities than the two-dimensional curves. This behavior is associated with the enlargement of the runoff area under the same operating conditions imposed in the 2D analysis. Reinforcing this analysis, it is possible to visualize the temperature influence on permeate mass flow (Table 7).        Figures 9-12 show the temperature fields within the separation module for feed temperatures 35, 55, 75, and 95 °C (Cases 5-8). A temperature gradient is observed, at annular space between the shell module and the ceramic membrane, caused by the heat transfer between the heated membrane surface and atmospheric air, by natural convection, thereby cooling the fluid at membrane proximity. Furthermore, higher temperature gradients are verified for z = 2.25 m and 3 m (close to outlet region), due to the gradual heat losses along with the axial position.  Figures 9-12 show the temperature fields within the separation module for feed temperatures 35, 55, 75, and 95 • C (Cases 5-8). A temperature gradient is observed, at annular space between the shell module and the ceramic membrane, caused by the heat transfer between the heated membrane surface and atmospheric air, by natural convection, thereby cooling the fluid at membrane proximity. Furthermore, higher temperature gradients are verified for z = 2.25 m and 3 m (close to outlet region), due to the gradual heat losses along with the axial position.          In Figures 9-12, significant temperature variations can be observed inside the module, due to the pre-established entry conditions, which are equal to ΔT ≈ 2 o C (Case 5), ΔT ≈ 7 o C (Case 6), ΔT ≈ 13 o C (Case 7), and ΔT ≈ 22 o C (Case 8). This behavior is associated with heat losses, by natural convection for the environment, since only the fluid mixture at the inlet has a higher temperature than the equipment wall. This behavior is associated with the oil transport to the proximity of the membrane surface by the shear forces, which is attenuated by increasing the temperature, due to the decrease in dynamic viscosity and density of the fluid, thereby reducing the accumulation of particles along the membrane surface.

Concentration Field
This behavior can be seen more adequately in Figure 17, which represents the oil concentration profile along the z-axis, in y equal 0.017 m, for temperatures ranging from 35 to 95 °C.  In Figures 9-12, significant temperature variations can be observed inside the module, due to the pre-established entry conditions, which are equal to ∆T ≈ 2 • C (Case 5), ∆T ≈ 7 • C (Case 6), ∆T ≈ 13 • C (Case 7), and ∆T ≈ 22 • C (Case 8). This behavior is associated with heat losses, by natural convection for the environment, since only the fluid mixture at the inlet has a higher temperature than the equipment wall. This behavior is associated with the oil transport to the proximity of the membrane surface by the shear forces, which is attenuated by increasing the temperature, due to the decrease in dynamic viscosity and density of the fluid, thereby reducing the accumulation of particles along the membrane surface.

Concentration Field
This behavior can be seen more adequately in Figure 17, which represents the oil concentration profile along the z-axis, in y equal 0.017 m, for temperatures ranging from 35 to 95 • C. In Figures 9-12, significant temperature variations can be observed inside the module, due to the pre-established entry conditions, which are equal to ΔT ≈ 2 o C (Case 5), ΔT ≈ 7 o C (Case 6), ΔT ≈ 13 o C (Case 7), and ΔT ≈ 22 o C (Case 8). This behavior is associated with heat losses, by natural convection for the environment, since only the fluid mixture at the inlet has a higher temperature than the equipment wall. This behavior is associated with the oil transport to the proximity of the membrane surface by the shear forces, which is attenuated by increasing the temperature, due to the decrease in dynamic viscosity and density of the fluid, thereby reducing the accumulation of particles along the membrane surface.

Concentration Field
This behavior can be seen more adequately in Figure 17, which represents the oil concentration profile along the z-axis, in y equal 0.017 m, for temperatures ranging from 35 to 95 °C.

Final Consideration and Limitations of the Model
In this work, the steady-state regime was used, being a relevant approach when it is desired to analyze the behavior of the separation process without the influence of time. As the main objective is to evaluate the influence of temperature on the flow-dynamics behavior, this approach is appropriate for the proposed study. However, we recommend a transient analysis for future works, to determine the time to change or cleaning the membrane.
The model developed was based on the work of Damak et al. [23][24][25] (used to validate this work). Therefore, we decided to work similarly. Then, the flow laminar regime was assumed. In general, this condition has been used in different experimental works. In addition, Equation (30), which was used to describe the concentration polarization, is only valid for R = 300-1000, i.e., under low Reynolds number condition. However, in industrial scale, outside of the membrane, turbulent flow is the dominant flow regime. Thus, studies under turbulent regime are recommended, since the turbulent flow could help to reduce the concentration polarization.
Furthermore, based on Damak et al. [23][24][25], it was decided to work with the series resistance model, considering the membrane resistance and the resistance due to the formation of the concentration polarization, as well as the Carman-Kozeny equation, to describe the pressure drop across the porous medium. Then, the effect of the membrane data and thickness were incorporated into the model, not being necessary to model the porous media explicitly.
Furthermore, when disregarding the pore obstruction and the fouling mechanisms, it is not possible to analyze their influence on the process. These effects are important, but the information is not found easily. However, within the proposed objective, these considerations do not compromise the results, since the model considers the rejection coefficient equal to one.
The temperature range for analysis was chosen based on the physics of the problem (temperatures below 100 °C), since high temperatures promote physical phenomena such as water evaporation and oil degradation. In industrial practice, all devices operate in isothermal conditions. Thus, new researches incorporating thermal effects play an important role.
With respect to the heat-transfer process, some aspects need to be explained. Based on the work of Damak and co-authors [23][24][25], the membrane was modeled entirely by boundary conditions. Then, it was necessary to use an appropriated thermal condition at the surface of the membrane. The fluid flows inside the system under action of one piece of pumping equipment; however, since the lower fluid velocity was found in the device (typically, less than 0.05) and Gr/Re 2 >> 1 [28], the consideration of free convection to analyze heat transfer, it is very appropriate. Furthermore, with respect to the choice of the numerical meshes, some important points must be explained. The number of elements used in 2D analysis was adequate. For this case, all predicted results were validated with the results reported by Damak and co-authors [23][24][25]. This information

Final Consideration and Limitations of the Model
In this work, the steady-state regime was used, being a relevant approach when it is desired to analyze the behavior of the separation process without the influence of time. As the main objective is to evaluate the influence of temperature on the flow-dynamics behavior, this approach is appropriate for the proposed study. However, we recommend a transient analysis for future works, to determine the time to change or cleaning the membrane.
The model developed was based on the work of Damak et al. [23][24][25] (used to validate this work). Therefore, we decided to work similarly. Then, the flow laminar regime was assumed. In general, this condition has been used in different experimental works. In addition, Equation (30), which was used to describe the concentration polarization, is only valid for R = 300-1000, i.e., under low Reynolds number condition. However, in industrial scale, outside of the membrane, turbulent flow is the dominant flow regime. Thus, studies under turbulent regime are recommended, since the turbulent flow could help to reduce the concentration polarization.
Furthermore, based on Damak et al. [23][24][25], it was decided to work with the series resistance model, considering the membrane resistance and the resistance due to the formation of the concentration polarization, as well as the Carman-Kozeny equation, to describe the pressure drop across the porous medium. Then, the effect of the membrane data and thickness were incorporated into the model, not being necessary to model the porous media explicitly.
Furthermore, when disregarding the pore obstruction and the fouling mechanisms, it is not possible to analyze their influence on the process. These effects are important, but the information is not found easily. However, within the proposed objective, these considerations do not compromise the results, since the model considers the rejection coefficient equal to one.
The temperature range for analysis was chosen based on the physics of the problem (temperatures below 100 • C), since high temperatures promote physical phenomena such as water evaporation and oil degradation. In industrial practice, all devices operate in isothermal conditions. Thus, new researches incorporating thermal effects play an important role.
With respect to the heat-transfer process, some aspects need to be explained. Based on the work of Damak and co-authors [23][24][25], the membrane was modeled entirely by boundary conditions. Then, it was necessary to use an appropriated thermal condition at the surface of the membrane. The fluid flows inside the system under action of one piece of pumping equipment; however, since the lower fluid velocity was found in the device (typically, less than 0.05) and Gr/Re 2 >> 1 [28], the consideration of free convection to analyze heat transfer, it is very appropriate. Furthermore, with respect to the choice of the numerical meshes, some important points must be explained. The number of elements used in 2D analysis was adequate. For this case, all predicted Energies 2020, 13,2092 18 of 20 results were validated with the results reported by Damak and co-authors [23][24][25]. This information can be verified by the good agreement between the data, as shown in previous work [9]. It is expected that the number of elements in the 3D mesh would be higher than that used in this work. However, because of the simultaneous solution of the governing equations as applied to multiphase flow in porous and non-porous media and the limitation in hardware (computational equipment), larger computational efforts and time were needed, and numerical convergence problems occurred (overflow), when a larger element number was used. Thus, as an alternative numerical procedure, the number of elements was reduced, sequentially, until we obtained a convergent solution with the highest possible element number of the mesh (217,352 elements).

Conclusions
Based on the numerical results of the water/oil separation process through the tubular membrane, it is possible to conclude the following: (a) The mathematical model was able to predict the fluid behavior during the water/oil separation process with accuracy, providing a good understanding of the flow inside the system; (b) The increase in temperature provided modifications in the fluid properties, reducing the feed velocity and modifying the concentration, pressure, and velocity distributions; (c) The increase in temperature inside the equipment reduces the pressure gradient, due to a decrease in the inlet velocity, directly influencing the fluid viscosity and density ratio, causing a reduction in the permeation velocity.
Author Contributions: All the authors contributed to the development, analysis, writing, and revision of the paper. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by CNPq, CAPES, and FINEP (Brazilian research agencies).