Study of the fluid flow pattern in a bubble column reactor for biodiesel production

The applications of bubble column reactor (BCR) are very important as multiphase reactors in process industry. In biodiesel production, one of the reactors that are used is BCR. The advantages of BCR are low operating cost and maintenance due to the compactness and no moving parts. It is important to understand the nature of hydrodynamics and operational parameters to characterize their operation, including gas superficial velocity and bubble rise velocity to make the design and scale-up process. Computational fluid dynamics (CFD) can be used to evaluate the performance of BCR at lower cost compared to the experimental setup for biodiesel production. In this work, a commercial CFD software, FLUENT 14.0 was used for modeling of gas-liquid flow in a BCR for biodiesel production. Multiphase simulations were performed using an Eulerian-Eulerian two-fluid model. In this study, the simulation was conducted by using three different temperature, which are 523 K, 543 K and 563 K. The CFD result predicts the turbulent kinetic energy, gas hold-up and the liquid velocity were fairly good, although the results seem to suggest that further improvement on the interface exchange models and possibly further refine the two-fluid modeling approaches are necessary especially for the liquid velocity and turbulent kinetic energy.


Introduction
Biodiesel is a fuel derived from vegetable oil or animal fat. It is a direct petro-diesel replacement that can significantly reduce gas emissions and provide energy independence [1]. Biodiesel is produced under chemical reactions which are transesterification and esterification. Transesterification is a process where triglycerides (TG) are reacted with alcohol, generally methanol (MeOH) to form fatty acid methyl esters (FAME). There are three steps reactions with the intermediate formation of triglycerides (DG) and monoglycerides (MG) resulting in the production of 3 moles of methyl esters (ME) and 1 mole of glycerol (GL) [2].
There are several types of reactors to produce biodiesel. One of the reactor that is widely used in the biodiesel industry is a bubble column reactor (BCR). They are excellent reactors for biodiesel production which require a large interfacial area for gas-liquid mass transfer and efficient mixing for reacting species. Producing biodiesel by using superheated methanol vapor in the BCR take place under atmospheric pressure and requires high temperature (250°C -300°C), therefore reactor was made of stainless-steel and covered by insulation [3]

Geometry Modeling
In this study, all the CFD simulation was conducted by using ANSYS FLUENT 14.0. For the geometric modeling, the design is constructed by using SOLIDWORKS 2014. The model is then imported to ANSYS Workbench for further analysis. The references of geometry were referred based on previous studies and analysis [2][3]. This interactive process is the first stage in pre-processing before generating the mesh. The dimensions and physical conditions of the BCR are given in Table 1.

Meshing
Meshing is one of the most important parts to determine the efficiency and effectiveness of an analysis. The function of meshing is to perform the Finite Element Analysis. In order to obtain more accurate results, the mesh must be refined. The easiest mesh refinement method is by reducing the element size throughout the modeling domains. After trials and errors, the suitable grid size and mesh type for the BCR simulation were determined which can obtain accurate results and shorten the calculating time. The mesh configurations are shown in Figure 2 and the number of nodes and elements are shown in Table 2.  Finer mesh can capture the flow details efficiently since it has the highest number of nodes and elements. That is the reason why the finer meshes can provide a better result compared to coarse mesh. In addition, the calculation of the problems also can be converged by using a finer mesh.

Boundary Condition
The boundaries of the BCR are either an inlet, outlet, or wall. The inlet is dispersed with a uniform gas velocity and a volume fraction of 1. The velocity at the inlet is 2 m/s. Besides, the outlets have a pressure equivalent to the atmospheric pressure, which is 0.1 MPa. The turbulent kinetic energy and dissipation rate are set to 0 at the inlets and outlets because the eff ects of turbulence at the boundaries are difficult to predict. The walls are specified with a no-slip condition for the gas and liquid phases.
The Euler -Euler model treats both phases as continuous phase, means that the gas phase does not move along the wall, therefore, no-slip boundary condition is applied to the gas phase along the wall. Table 3 summarize the boundary condition for BCR simulation.

Computational Model
The Eulerian multiphase model allows for multiple interacting phases by assuming that each phase behaves as interpenetrating continua [4]. The continuity equation is formulated for each phase without exchange between phases after Reynolds averaging term where uq and αq is the velocity and volume fraction of phase q th respectively.
By comparing to the mass conservation, the momentum conservation is described by the Navier-Stokes equation extended by the phase volume fraction. Such an equation for phase 'q' is as follows: The terms on the right-hand side describe all forces acting on the phase 'q' of a fluid element in the control volume. The effect of virtual mass and lift force in this work has been neglected. From equation (2), τq is the stress tensor of phase q th and calculated as follows: where μq and λq are shear and bulk viscosity of phase q respectively.
The exchange coefficient for gas-liquid flow has been considered as following equation: The momentum exchange due to drag forces is: where CD is the drag coefficient, to compute this coefficient the Schiller-Naumann relation has been used. To apply the effect of magnetic field on the fluid flow, the Lorentz force source term has been added to the momentum equation [5].

Result and Discussion
In this section, some of the characteristics of bubble dynamics were explored. It is now necessary to examine the regime of validity of those analyses. The simulation was conducted under 3 different temperatures; 523 K, 543 K and 563 K. Then the flow regime and velocity profile can be obtained from each temperature. Turbulence modeling is an important factor when modelling two-phase flows.
The standard k-ε model was also tested for turbulent dispersion. (1)

Flow Pattern
The flow pattern inside the reactor can be obtained from each temperature. By using XY Plane as the viewing angle, the production of bubbles and flow pattern inside the reactor can be observed. Figure 3 shows the contours of vapor volume fraction which describes the distribution of methanol bubble in the oil as result of CFD simulations. Red color appears when the volume fraction of vapor is 1 and the blue color shows that vapor volume fraction is 0. It is assumed that the methanol bubble was in vapor form while the oil was in liquid form.
At 2 seconds flow, the bubbles have reached at the top of the reactor and begins to circulate in the reactor. Based on the three conditions from Figure 3, the circulation behavior of the liquid is the same as the model made by Joshi and Sharma which is a Donut model [6]. The bubbles dispersed from the sparger will always move upwards in the form of swarm in a priority direction, and then creating a large liquid vortex structure. As the time increase, the priority direction suddenly changes over into another which could be any direction inside the reactor.

Velocity Vector
To understand the flow characteristic in the reactor, velocity vectors can provide an excellent visualization of the flow and can also show the direction of the flow. In order to simulate the hydrodynamic behavior of the system in turbulent regime, the calculations were started with the laminar model with the liquid at rest [7]. When the liquid velocity in most parts of the reactor was greater than zero, the k and ε fields were initialized and the turbulence model was started. Figure 4 shows the velocity vector at the wall and the center of the reactor at 523 K. From Figure 4, the direction of the vector in the center region is in reverse direction compared to the wall region. This is because most part of the center region in the BCR was filled with the upward liquid flow and it returned to the bottom along the column wall.

Gas Holdup
The effect of gas velocity inlet to gas holdup and surface contact area was analyzed. Gas holdup is the basic parameter indicating the hydrodynamic characteristics of BCR [8]. It affects directly the geometric sizes, and the gas-liquid interfacial, thus affect the mass-transfer rates of BCR. Therefore, it is one of the important parameters for the BCR design. As a result, the contact surface area and gas holdup increase with increasing gas velocity. The effect of temperature on gas holdup is shown in Figure 5.

Velocity Analysis
The velocity analysis based on the graph has been plotted according to the vertical plane and horizontal plane as shown in Figure 6. By plotting the results from both planes, the observation between all conditions can be compared with each other. From the analysis, the best and suitable temperature can be selected. In addition, the graph only focused on the velocity result because the velocity can determine the flow pattern in the BCR either it is highly turbulence or less turbulence regime during the transesterification process. The direction of the velocity flow is upward towards the top of the BCR. Figure 5 shows the plotted velocity against distance at 523 K, 543 K and 563 K.  Figure 6 above, it shows that the velocity of the fluids in mostly increase until it reach 0.14 m reactor level, then the velocity reduces ranging below 0.2 m/s. This might due to the bubble has reached at the entrance of the outlet pipe. After that, the flow velocity begins to increase again until they reach the top of the outlet of the reactor.

Validation
Validation may refer to the comparison between the experimental data and the simulation result based on the study in Ref. [2] and [7]. From the observation, the velocity was resulted in different value of simulation from the experiment. In these problems, the error might cause by a high flow circulation by the fluids in the reactor, thus making the standard k-ε turbulence model cannot deliver the exact predictions for the velocity results. In addition, the geometry drawing of the reactor might be less accurate from the exact dimensions. Therefore, it will affect the result obtained from the simulation. Meanwhile, the errors might come from the experiment itself. During the measurements, some errors might be occurred that affect the experimental result. Figure 7 shows the comparison between experimental results with the CFD simulation in the term of gas holdup.

Conclusion
In this study, the fluid flow inside the BCR has been simulated and analyzed by using CFD software. There are two phases which are triglyceride (liquid) and MeOH (vapor), therefore using Volume of Fluid (VOF) model is suitable in this simulation. Moreover, the standard k-ε model is used to predict the turbulence flow inside the BCR during the mixing process. The results obtained with the 3D version of the k-ε turbulence model are good and acceptable agreement with experimental results. The inclusion of the standard k-ε turbulence model is useful to describe the instantaneous large-scale vertical flow structure correctly.
The effects of inlet gas velocities to the gas holdup were accounted. This parameter is very important to describe the interaction between two-phase fluid flow in the BCR. Many correlations for predicting gas hold up can be found in literatures. Experimental data shows that the total gas hold up is increased with increasing superficial gas velocity. From the simulation result, it is proven that the increasing of inlet gas velocities cause increasing of the gas holdup. The large bubbles have higher rise velocity and decrease the gas hold up.