CFD simulation and PIV validation of the gas/liquid behavior in an UASB reactor

As the world population increases, the need to develop more efficient wastewater treatment systems requires the use of new technologies. Software aided project and optimization of bioreactors and bioprocesses have become a matter of interest in recent years, especially due to the advance in the state-of-the-art of computational resources. This work aimed to perform gas/liquid numerical simulations using the Fluent 16.2 software and to validate this model through Particle Image Velocimetry (PIV) and shadow imaging techniques. Eulerian-Eulerian, laminar, tridimensional and transient simulations were carried out. The results for the mass imbalance for the gas and liquid phases, gas volumetric fraction, gas velocity, bubble size, liquid magnitude and upflow velocity and the velocity profiles for the liquid phase were successfully validated against experimental data. Concerning the dispersed phase, it was found a difference of 4.37% for the gas volumetric fraction between experiments and simulations. Simulated results showed a difference for the bubble mean velocity of 1.73% when compared with shadow imaging results. No coalescence was observed along the experiments, and the flow regime was characterized as dispersed bubble flow. Regarding the liquid phase, it was found a difference of 3.2% for the mean velocity, between simulated and PIV results. Simulated and experimental velocity profiles showed a better agreement at the center of the reactor. Some differences were observed in those profiles, due to geometry simplifications assumed in order to get a better mesh. Considering the good agreement between simulation and experiments, the model was considered validated.


INTRODUCTION
Anaerobic digestion based processes have become a popular kind of bioprocess for treating wastewaters, especially due to their ability to combine the efficient treatment of effluents and the generation of renewable energy. In the late 1970s, a reactor known as upflow anaerobic sludge blanket (UASB) was developed and became widely popular, mainly due to its relatively simple design allied to high treatment efficiencies (SOUZA, 1986;LETTINGA;POL, 1991).
Considering the current need for renewable energy sources as well as the increase in the wastewater generation, the development of new technologies, which allow higher treatment efficiencies and rates, has drawn the attention of researchers from different fields. Most of the works on that field use experimental studies to assess the influence of changes in the process control and operation, e.g. the use of recycle streams, series reactors, changes in the hydraulic retention time (HRT) (EL-SHEIKH et al., 2011;INTANOO et al., 2014;RIZVI et al., 2015;FERRAZ JUNIOR et al., 2016). Changes on the biomass source and the supplementation of the process aiming higher biogas production are also under study by some authors (JANKE et al., 2016;BARROS;DUDA;OLIVEIRA, 2016).
Even though UASB reactors are considered simple in means of construction and operation, predicting the hydrodynamic behavior in biogas producing bioreactors is challenging, considering that it consists of a complex multiphase flow, formed by a liquid, a solid and a gas phase. Nonetheless, the understanding of the hydrodynamic behavior in UASB reactors plays an important role in the improvement of these reactors technology. The biogas yield is dependent on the mixing quality within the reactor, once it induces to a higher contact time between biomass and wastewater.
On the other hand, in cases where the upflow liquid velocity is low, the auto-mixture within the system is promoted mainly due to the gas velocity. In other words, the gas/liquid momentum transfer in UASB reactors has an influence on the overall anaerobic process of the reactor. In general, the interfacial momentum transfer between a gas and a liquid can be due to these three forces:

Computational fluid dynamics model
Considering that, the main anaerobic reactions take place at the lower and middle parts of the reactor, the volume control for the simulations was assumed as showed in Figure 1. Due to the need of a high mesh quality, simplifications concerning the influent distribution system were assumed; also, the conical region of the reactor was not simulated. The commercial software ICEM CFD was used for the domain discretization through structured blocks meshing techniques, and hexahedral elements, in the three-dimensional space.
Taking into consideration the work of Lima et al. (2011), which showed that the flow in a similar UASB reactor could be considered symmetric, as well as the computational costs implicated in simulating a full reactor, it was simulated only half of the reactor (longitudinal section). Aiming to obtain a mesh independent result for the simulations, a mesh independence study was performed.
UASB reactor gas/liquid CFD simulation and PIV validation Thus, three meshes with different numbers of elements were generated: mesh 1 (1,035,064 elements), mesh 2 (528,000 elements) and mesh 3 (248,720 elements). The Grid Convergence Index (GCI) method, proposed by Roache (1994), was used to assess the mesh independence through the analysis of two performance parameters. Therefore, the volume upflow velocity for the gas and the gas volumetric fraction within the reactor were chosen as those parameters. Considering the GCI results as well as the computational costs involved (unpublished results), the mesh 2 (528,000 elements) was chosen for this study.
The commercial software package Fluent 16.2 was used in order to perform the numerical simulations. A two-phase, Eulerian-Eulerian For the E-L model, higher is the dispersed phase volumetric fraction, higher is the computational effort needed, while in the E-E approach, this effort remains constant. Even though the gas is present in a small volumetric fraction, in future works, it is intended to add the solid fraction to the simulation of the UASB reactor, thus increasing significantly the dispersed phase volumetric fraction.
A Reynolds number of 7.45 was calculated using the designed mean upflow velocity for the liquid inside the UASB reactor (2.525 × 10 -5 m.s -1 ), the water properties at 293.15 K and the dimensions of the tubular section of the reactor (0.3 m of diameter). Thus, a laminar model was used to perform the simulations. The liquid phase was considered to be water at 293.15 K, and the standards properties from the software database were used. Considering the biogas as a mixture of 65% methane (CH 4 ) and 35% carbon dioxide (CO 2 ), the gas phase properties were calculated using the harmonic mean of their properties. Table 1 shows the physical properties used in the simulations for both: dispersed and continuous phase.
Concerning the simulations, it was not considered the mass and energy transfer among the phases. The movement and mass conservation for each phase were governed by the fluid-fluid moment conservation Equations 1 and 2, which are solved by Fluent (ANSYS, 2015): Where: q = continuous phase; p = dispersed phase; τ = phase stress-strain tensor; According to many authors, such as Lima et al. (2011), Cheng andZhu (2005), Baten, Ellenberger and Krishna (2003), and Simcik et al. (2011), the drag force is the main responsible for the interfacial momentum transfer between the phases. Therefore, in this study, only the drag force was considered.
Depending on the kind of the flow and on the characteristics of the dispersed and continuous phases, a correlation for the drag force coefficient (CD) should be chosen. As the gas bubbles were spherical and the distribution was sparse, the model proposed by Schiller and Naumann (1933) was selected, and the CD was calculated using Equation 3: Where the Reynolds number for the bubble is calculated using A liquid inlet velocity of 7.11 × 10 -3 m.s -1 and a gas inlet velocity of 1.12 × 10 -3 m.s -1 were assumed, resulting in a liquid flow rate of 6.43 × 10 -3 m 3 .h -1 and a gas flow rate of 4.15 × 10 -2 Nm 3 .h -1 . At the gas outlet, a degassing boundary condition was used, thus allowing only the gas to flow through it. The absolute pressure (Pabs) is the sum of the pressure at the outlet (Pstat) and the operation pressure (Pref), as showed in Equation 5. Therefore, the pressure at the liquid outlet boundary was assumed to be 0 Pa, and an operation pressure of 101,325 Pa was used, thus the absolute pressure would correspond to the atmospheric one at the liquid outlet boundary.
A symmetry boundary condition was assumed at the center of the reactor, the software assumes zero flux of all quantities across that type of boundary. A non-slipping boundary condition was assumed at the walls, it considers that the velocities at these in the X, Y and Z directions are equal to zero.
The CFD model is solved using the phase coupled SIMPLE algorithm in order to solve the velocity pressure coupled differential equations. Aiming to discretize the governing equations into algebraic ones in the final solution, the second order upwind scheme was used for the momentum and the QUICK for the volume fraction.
The adaptive time stepping method was used along the simulations, together with a maximum number of iterations per time step of 25. As convergence criteria, residuals of 1e -4 were adopted, and the stability of some control variables was assessed. Transient simulations were performed.

Reactor
For the experimental validation of the CFD model, the PIV technique was applied using a small-scale UASB reactor. It was previously designed for treating sugarcane vinasse with a chemical oxygen demand (COD) between 15,000 and 65,000 mg O 2 .L -1 (CORTEZ; FREIRE; ROSILLO-CALLE, 1996) and thus generating biogas. A reactor 2.12 m high and with 0.3 m of diameter was used. A usable volume of 140 L was assumed. It was built in poly(methyl)methacrylate (PMMA), thus allowing an internal view of the flow within the reactor. The laser beam generated by the PIV can suffer some refraction effects when it reaches the curved wall of the UASB reactor. In order to reduce these effects, the reactor was placed in a hexahedral vessel, which was filled with tap water.
The liquid entered the reactor through a distributor with eight pipes, placed at the bottom of the reactor, which was designed in order to generate a better mixing between gas and liquid, as proposed by Bastiani et al. (2016). Four concentric rings were placed at the bottom of the reactor, where the gas was inserted. According to Bastiani et al. (2016), this was the alternative that showed the best results on mimicking the effect of the gas being generated inside the reactor. When using PIV techniques, it is important to avoid turbid fluids, in order to achieve better results for the images of the flow fields. Therefore, the use of wastewater and biomass to generate the biogas on the reactor would have a strong interference on the results. The reactor was filled with tap water, and the gas used for the experiments was compressed air at room temperature. The water flow was controlled using a hydraulic pump manufactured by Netzsch, together with a frequency inverter (Weg -CFW 10). A liquid flow of 6.42 L.h -1 was assumed.
For the gas flow control, a rotameter -Applitech ® 1900 -, in combination with a differential pressure controller, was used to maintain a flow of 0.041 Nm 3 ·h -1 .

Particle Image Velocimetry and shadowgraphy set up
In order to gather relevant and accurate data for the validation of the CFD model, a PIV system was used, which is a non-intrusive method that allows the flow visualization and the particle movement analysis.
This method was used to assess the flow velocity and profile for the liquid phase. It consists in the illumination of a section of the flow in two near instants of time, using a laser. The position of tracing particles is then registered using two cameras, and the data is processed in specific software for the analysis of the flow.
In this work, a stereoscopic 3D PIV acquired from LaVision was used. It consists of two cameras ImagerProSX 5M, with Nikon lenses.
It was necessary the use of a 540 nm cut-off filter on the lenses, to avoid possible damages caused by the laser light scattering initiated by bubbles. The laser generator used was the Quantel Double-Pulsed Nd:YAG EverGreen EVG 00200. Fluorescent 20-50 μm tracers were seed into the flow in order to allow the flow analysis. The images collected were processed in a computer using a processor Intel ® Xeon E5620, a video card Zotac/NVIDIA GeForce GT 610, operational system Windows 7 and a programmable timing unit (PTU) VZ Trigger PTU V9.0. The DaVis ® UASB reactor gas/liquid CFD simulation and PIV validation 8.2.2 software was used for the images' treatment and processing, in order to obtain flow fields. The PIV measurements were preceded by a calibration using a calibration plate. This plate was placed inside the reactor and then the cameras were focused on the center of it, aiming to mark on the software three pre-determined points, in order to proceed with the software aided calibration. Following this calibration, the reactor was seed with fluorescent particles. The images were collected using the configuration showed in Figure 2.
The cameras and the laser were placed in vertical rails, allowing capturing the images throughout all the reactor height. The cameras were placed at 60 and 120º within the reactor central longitudinal plane for the tridimensional results optimization, and the laser was placed in a position that it could lighten this plane. Considering that the reactor height was much bigger than the lenses capturing field, the reactor was divided into seven sections for the images capture. In order to measure the particle displacement, two images are captured for each camera, one for each laser pulse. The delay time between the laser pulses changed along the sections, considering that higher is the flow velocity, smaller should be the delay time in order to allow a particle displacement of around 10 pixels between the two pulses.
Aiming to obtain more precise results for the averaged vector field, it is important to obtain as many images for each section as possible, allowing the producing of a high number of instantaneous vector fields to calculate a better time averaged result. On the other hand, a high number of images demands a longer processing and post processing time. Therefore, it was decided to capture 200 images of the flow for each section of the reactor. For each image, four frames were obtained, one frame each camera for each laser pulse, thus allowing the post processing of the images in order to obtain the instantaneous vector fields, that would be used to further calculate the time averaged vector fields.
Aiming to gather data about the gas phase, the PIV system was changed to a shadowgraphy configuration, which allows capturing the gas bubble size and velocity. The same hardware and software were used with a different configuration, but no tracer particles were fed into the reactor. Figure 3 shows a scheme of the shadow imaging system used in this work.
The shadowgraphy technique consists of a pulsed backlight illumination (a dispersed laser light obtained using a laser diffuser), which is aligned with a camera. The laser illuminates the reactor section in two different moments, and the camera captures the "shadow" of the particles of interest (the gas bubbles in this work) in each laser pulse. Then the software is used to calculate the displacement of the observed bubbles between laser pulses, obtaining information about bubbles velocity and position. Over this study, the shadowgraphy measurement of 11 vertical positions was performed, and the results were then tabulated.

RESULTS AND DISCUSSION
Considering the designed values for the gas flow rate and the liquid flow rate, the mass flow rate at the inlet for the gas and the liquid phases were calculated, using the density of each phase. Table 2 sums up the results for the mass flow rates and the mass imbalances for the analytical and computational cases. It shows a difference for the gas mass flow   D' Bastiani, C. et al. rate of 0.39%, and for the liquid mass flow rate of 0.05% between simulated and analytical results at the inlet. It represents a not significant difference, thus showing that the model is suitable for the simulation.
The calculated values for the mass imbalance results also presented a good one, considering the flow conditions. The simulation showed a gas volume fraction within the reactor of 0.0007627044. In order to validate this result, the volume of gas inside the reactor was experimentally determined. The amount of gas inside the reactor can be considered, as the difference in the volume occupied by the fluids with and without the gas. Aiming to measure that, the pump of the liquid was turned on, and the height of the reactor occupied by the liquid was measured. Then the gas was injected into the reactor, and the height was measured again. The difference between those heights was used to calculate the volume occupied by the gas and thus the gas volumetric fraction. An amount of 9.47192 × 10 -5 Nm 3 of gas in the reactor was obtained or a volumetric fraction of 0.0007307548, a difference of 4.37% between experimental and simulated results. Considering that, it is possible to say that the numerical and experimental results for the gas volumetric fraction show a good agreement.
Considering the reactor dimensions and the gas and liquid flow rates, the superficial velocities were calculated using Equation 6: Where: Uis = phase i superficial velocity (m·s -1 ); Qi = phase i flow rate; A = transverse area of the reactor.
The gas superficial velocity calculated was 1.6 × 10 -4 m.s -1 , and the liquid superficial velocity was 2.5 × 10 -5 m.s -1 . Mishima and Ishii (1983) developed new flow-regime transition criteria for upward two-phase flow in vertical tubes, and presented a flow-regime map based on the gas and liquid superficial velocities.
Although a different tube diameter was used in the present study, it was observed a bubbly flow for the calculated superficial velocities, as shown in Figure 4, in accordance with the Mishima and Ishii (1983) flow-regime map. In bubble columns, it was observed by Chen, Reese and Fan (1994), and Díaz, Montes and Galán (2008) mainly two flow regimes: the dispersed bubble and the coalesced bubble. Figure 4 shows the bubble distribution in a 70 × 60 mm region of the reactor, where is possible to visualize the dispersed bubble flow regime. The image was collected using shadowgraphy.
Concerning the shadow imagining experiments, the results for the bubble size for different vertical positions were obtained and are showed in Figure 5, together with their standard deviation. Results show that the diameter is not increased proportionally with the vertical position.
Therefore, it is possible to conclude that there is no bubble coalescence within the reactor. According to Othman, Hamzan andTerry (2011), andDíaz, Montes andGalán (2008). At low Ugs there is no bubble coalescence, due to the dispersed bubble regime.
Both, low gas and liquid flow rates within UASB reactors are a consequence of effluent COD, reactor design and organic loading applied. Higher is the COD to be removed in the treatment and the organic loading applied, higher is the HRT of the effluent, thus the liquid flow rate is reduced. On the other hand, the gas flow rate increases with higher COD and organic loading, due to the increment on the gas yield. Thus, the flow regime in UASB reactors is also a function of the aforementioned characteristics (COD, reactor design and organic loading).  In previous works developed by one of the authors (unpublished results), three CFD simulations were performed using different bubble sizes: 1 mm, 2 mm and 3 mm. The same simulation settings previously presented in this work, as well as the same reactor, were used in these simulations. The results showed a correlation between bubble size and gas velocity, which is presented in Equation 7: Regarding the liquid phase, the velocity profiles obtained from the PIV experiments for different vertical positions along the reactor were compared with the CFD results. Figure 6 shows the PIV and the CFD results for the mean velocity for a plane along the central axis.    the walls, causing the area for the liquid to flow to be reduced and, like this, the liquid velocity to be higher than the calculated superficial velocity. While at the center of the reactor, there is a smaller gas volumetric fraction, so the liquid flows downwards through that space ( Figure 8A). At the top of the reactor, the gas is led by the deflectors towards its outlet, where it leaves the volume control ( Figure 8B), as it is expected for UASB reactors.
Throughout this study, numerical and experimental results for different parameters were obtained, such as: mass imbalance for both phases, gas volumetric fraction, gas velocity, bubble size, liquid magnitude and upflow velocity and the velocity profiles for the liquid phase. Considering the complexity of the flow, the differences between numerical and experimental results presented over this study are adequate to validate the numerical model, thus allowing further analyses and optimizations.

CONCLUSION
According to the results presented in this study, one can conclude that the use of PIV techniques allowed the successful validation of the numerical model proposed for the simulation of the gas/liquid interactions in an UASB reactor. It was possible to notice that there was no bubble coalescence within the reactor. Furthermore, it was observed the existence of liquid internal recirculation, driven by the gas upflow within the reactor. Differences of less than 5% between numerical models and validation results were observed, except for the upflow liquid velocity, which showed a difference of 9.18%.
The good agreement between simulation and validation results will allow the use of the described numerical model for further gas/ liquid simulations, and the inclusion of the solid phase in forthcoming simulations of UASB reactors.