A CFD Comparative Study of Bubbling Fluidized Bed Behavior with Thermal Effects Using the Open-Source Platforms MFiX and OpenFOAM

This work studies the performance of two open-source CFD codes, OpenFOAM and MFiX, to address bubbling fluidized bed system at different temperature and heat transfer conditions. Both codes are used to predict two parameters that are relevant for the design of fluidized units: the minimum fluidization velocity as a function of the temperature of the bed and wall-to-bed heat transfer coefficient from a lateral wall and from internal tubes. Although the CFD solvers are structuraly similar, there are some key differences (available models, meshing techniques, and balance formulations) that are often translated into differences in the fields prediction. The computational results are compared between both codes and against the experimental data. The minimum fluidization velocity can be correctly predicted with both codes at different temperatures while, in general, for the heat transfer and the fluidization patterns, MFiX shows slightly more accurate results compared to OpenFOAM but with low versatility for meshing curved geometries which might translate into higher computational costs for the same level of accuracy.


Introduction
For the design of fluidized bed systems, the minimum fluidization velocity is arguably the most important variable [1,2] and can be generally defined as the minimum superficial velocity at which the pressure drop through the bed is equal to the bed weight per unit cross-section. A large amount of experimental work has been carried out on this parameter, and many correlations have been proposed for its prediction in the literature [3]. Regarding fluidization with heat transfer, the thermal uniformity is one of the main features of bubbling fluidized beds. This condition is caused by the presence of gas bubbles that induce a high amount of solids recirculation. The same mechanism produces high heat transfer coefficients towards submerged objects, establishing thermal gradients in a narrow region close to the surface of the object. In this sense, the internals are incorporated into fluidized beds for different purposes. In some cases, they are incorporated to add or extract heat from the bed using vertical or horizontal tubes (FBHE, Fluidized bed heat exchanger). In other cases, they are incorporated to prevent the growth of bubbles and, in this way, influence their average size, determining lower speeds during their ascent and eventual passage through tube bundles located above [4]. In any case it is necessary to evaluate the bed-to-surface heat transfer coefficient beforehand to carry out the design of the equipment [5]. The estimation of the minimum fluidization velocity can be carried out by employing correlations. Pattipati and Wen [6] showed that the minimum fluidization velocity (U m f ) is a function of temperature and can be correlated with the properties of the fluidizing gas that depend on the same variable. During their experiments, they observed that U m f decreases when the temperature increases for diameters of sand particles smaller than 2 mm, while the opposite occurs for particles of greater diameter. Likewise, the authors also concluded that the correlation of Wen and Yu [7], developed at room temperature, was valid for the predictions of U m f at elevated temperatures. Regarding the heat transfer coefficient, its estimation through correlations is not so simple. The wall-to-bed heat transfer coefficient h is the result of a combined mechanism of convection and radiation for both gas (interstitial and within bubbles) and for particles.
For the evaluation of both of these parameters, Computational Fluid Dynamics (CFD) techniques come as a nonexpensive tool complementing, and sometimes even replacing, the experimental approaches. The determination of U m f using CFD in systems at high temperatures has been studied by various authors. Gosavi et al. [8] studied systems with temperatures between 30-600 • C for lithium titanate (Li 2 TiO 3 ) spherical particles, belonging to group B of the Geldart classification [9], with air as the fluidizing agent. The simulations were developed in two dimensions using the Eulerian Two-Fluids Model (TFM). The model predicted the minimum fluidization velocity with 95% of accuracy when compared to experimental observations. Additionally, the authors concluded that the model is capable of predicting the decrease in U m f with increasing temperature. Shao et al. [2] used a 3D model with an Eulerian-Lagrangian approach to predict the minimum fluidization velocity at high pressure and temperature, with ranges between 0.1 and 4 MPa for pressure and 25 and 800 • C for temperature. The model was validated with experimental values reported in the bibliography. The authors concluded that the CFD model is suitable for the prediction of U m f and that it is also an inexpensive and fast option, compared to the determination of U m f experimentally. On the other hand, the study of wall-to-bed heat transfer has been studied by different authors using CFD [10]. Moreover, a phenomenological heterogeneous model to predict the heat transfer rates between bubbling fluidized beds and immersed surfaces was reported by Mazza et al. [11,12]. One aspect to consider in modeling fluidized beds with heat transfer to or from surfaces using TFM is that the thermal conductivities of the fluid phase and the solid phase (κ g and κ s ) should be interpreted as effective transport coefficients [13]. The direct use of the molecular thermal conductivities of the solid and the gas results in an overestimation of the energy transferred [14]. The effective conductivity model used in most of the cases reported in the literature is the Zehner and Schlünder model [15], commonly regarded as the standard approach [16]. Another relevant issue comes from the high degree of refinement in the heat exchange zone necessary for the correct resolution of the temperature field and, thus, the heat transfer. The first authors to carry out simulations with these characteristics were Gidaspow and Syamlal [17] and later Kuipers et al. [18].
This work seeks to determine both of these relevant parameters using two widespread open-source codes for CFD simulation: MFiX [19] and OpenFOAM [20]. Both of these codes are available and free so any trained user can download the software, install it on a personal computer or work station and use them for the study and design of fluidized bed systems. Therefore, evaluating the performance of both codes for addressing fluidized beds with heat transfer becomes specially important. It is the purpose of the present work to determine the accuracy of these codes and draw some conclusions and recommendation when they are used to predict the minimum fluidization velocity and heat transfer coefficient for different arrangements, validating the results with experimental data available in the literature.

Computational Model
This section describes the continuum equations that are part of the Two-Fluid Model (TFM) implemented in the open-source codes OpenFOAM [20] and MFiX [19,21]. For the sake of simplicity, the equations and models presented below are assumed to be formulated similarly in both codes and comments are made upon the differences.

Continuity Equations
The mass conservation equations for both phases can be written as: In practice, only one of the phase volume fraction is solved, and the volume fraction of the remaining phase is computed by considering: In addition, the sum of both equations gives rise to the continuity equation of the mixture, which is written as: This is only true when both phases are considered to be incompressible. Since the coupling between velocity and pressure is performed in a segregated manner, Equation (4) is used alongside the momentum equations to formulate an equation for the pressure field, following the general structure of the SIMPLE algorithm for multiphase flows [22][23][24][25].

Momentum Balance
The momentum balance for both phases may be written as: This general formulation, particularly the momentum balance for the solids phase, is based on the work of Ishii [26]. Here, the stress tensors may be written as: The interphase momentum transfer is given by the drag forces, and the drag coefficient is computed based on the Gidaspow model [17]: where:

Granular Rheology
The current model is based on treating both phases as an interpenetrating continua. Therefore, under this approach, the rheology of the granular phase needs to be properly modeled. For the low concentration of particles, the kinetic theory of granular flow [27] brings closure to the equations by introducing the granular temperature field (θ) which is used to compute the granular phase viscosity and obeys an energy balance equation: The parameters involved are defined as [27][28][29]: For high concentrations, the grains are in contact each other and rubbing and friction take place. For these conditions, the frictional theory based on soils mechanics [30,31] serves as a modeling approach for the solids pressure and solids viscosity: Here, the frictional pressure is computed following the approach used in MFiX [32], while the solids viscosity is computed following the work of [31]:

Internal Energy Balance
Both phases obey an internal energy balance which predicts that the rate of change of internal energy is equal to the changes due to convection, diffusion and heat transfer between phases. This might be written as: Here, the thermal conductivities are not a property of each phase material but an effective conductivity based on the current phase concentration and can be computed based on the model Bauer and Schlünder [33]: where and Then, On the other hand, the heat transfer between phases is calculated based on Gunn's correlation [34]:

Numerical Method
The aforementioned models are solved using MFiX v21.3.2 [19,21] and OpenFOAM v20.12 [20]. All of them are already available in the standard distribution of MFiX, while in OpenFOAM, the heat transfer model between phases Equation (28) and the effective conductivity model of Bauer and Schlünder (Equations (22)-(27)) were implemented for this work.
Both computational codes use the same approach for addressing the mathematical model. They are based on the Finite Volume Method (FVM) where both phases are treated as incompressible [35,36], and a SIMPLE-based algorithm [37] is used for the segregated coupling of pressure and the velocities of each phase. OpenFOAM allows one to perform iterations to enforce the mass balance within following the approach of PISO [38]. Moreover, the momentum equations are coupled based on the Partial Elimination Algorithm (PEA) [39,40].
It is worth mentioning that both codes use different meshing techniques. While Open-FOAM have a dedicated mesher tool and can import grids generated by other softwares, MFiX relies only on its own mesher tool, which is based on generating structured grids and the cut-cell technique for addressing curved surfaces. This difference becomes very relevant for addressing industrial-scale problems with curved surfaces. In general, for these situations, a uniformly highly refined grid might become unaffordable; therefore, MFiX would rely on a coarser grid which, in the presence of wall heat transfer effects, might not be enough. This issue is addressed with practical examples in the following section.

Results and Discussion
The following tests are selected based on the availability of experimental data but also with the intention of having simple geometries to validate the numerical approach.
The physical properties and parameters used for each test are summarized in Table 1, and the numerical parameters and mathematical models involved in these cases are included in Table 2. The first test case is based on the experimental setup of Subramani et al. [41]. In this work, minimum fluidization velocities were determined with the bed at different temperatures and filled with Geldart B particles. The experiments were carried out on a cylindrical bed made of silica glass with an internal diameter of 2.8 cm and a length of 25 cm, and the temperatures ranged from 273 to 973 K. The air is preheated before entering the bed at the corresponding temperature.
For the computational simulations, a mesh convergence analysis was performed for each software, resulting in a o-grid type of mesh consisting of 44,000 cells for OpenFOAM and a grid of 35,000 structured cells for MFiX, based on the cut-cell technique. These refinements were selected following an a priori analysis of mesh convergence and were proven to produce a good balance between the computational costs involved and the accuracy of the numerical solution for these conditions. All the physical and numerical parameters involved for the simulations are described on Tables 1 and 2. Both codes required around one hour of overall computational time in a single CPU to obtain a statistically steady solution of the pressure field (each point on Figure 1). Figure 1 shows the fluidization curves obtained with MFiX. Each point corresponds to the pressure drop obtained for a simulation with a fixed superficial velocity. Here, it can be observed that the qualitative trend of having smaller values of U m f as the temperature of the bed is increased. The U m f prediction with both codes is shown in Figure 2 along with the experimental results of Subramani et al. [41]. The values shown correspond to a graphical intersection between a linear fitting of the pressure drop values of the packed bed region and the fluidized region of the fluidization curves. Both codes show a slight underestimation of the U m f but, in general, in good agreement with the experimental results with a maximum error of 10%. Different topics contribute to generating the differences observed on the predicted values of the minimum fluidization velocity from OpenFOAM and MFiX. Without excluding some others, it must be mentioned that the momentum balance formulations are not strictly identical in both software. In addition, even if the coupling between phases is based on the Partial Elimination Algorithm [39] in both codes, there are still some differences in the formulations. Namely, as explained in Section 2.5, the algorithm in OpenFOAM is designed based on (a multiphase version of) the PIMPLE method, which is a combination of SIMPLE [37] and PISO [38], unlike MFiX, which uses the SIMPLE method directly. Using two or more PISO inner iterations per SIMPLE iteration enforces the mass balance per time step, increasing the convergence of the segregated coupling between pressure and the phase velocities, which can also be achieved by modifying the pressure under-relaxation, as is performed in MFiX.

Test 2: Heat Transfer from a Vertical Wall
In this test, the heat transfer coefficient is estimated based on simulations for a problem based on the experimental setup of Yusuf et al. [42]. The problem consists of a pseudo-2D fluidized bed with a jet inlet of high velocity (U = 16.6 m/s) in the bottom part in contact with the lateral wall, as shown in Figure 3. All the walls are adiabatic except the lateral right wall, which is at 333 K and the air inlet is at 293 K. The rest of the inlet at the bottom of the bed is set at minimum fluidization velocity (U = 0.18 m/s). In the experiment, this condition is usually achieved by using an air distributor consisting of a perforated plate in the whole base of the bed except for the jet inlet part. In the simulations, this is modeled by imposing a fixed velocity which is calculated by dividing its value by local phase fraction. The dimensions of the bed is 0.2 m of width, 0.7 m of height and 0.025 m of thickness, and the solids phase consists of glass spherical particles of 0.491 mm of diameter. The rest of the parameters for the simulation are summarized in Table 2.
The grid used for both codes consists of uniform refinement in the vertical direction and a linear grading of refinement in the horizontal direction with smaller cells closer to the hot wall, as shown schematically in Figure 3. Table 3 shows different grid refinements and how the solution is affected by it.  Table 3 shows the grid refinement analysis for MFiX, resulting in the adoption of mesh 3 for this test. The difference of meshes between codes translated into different overall computational times (although a uniform time-step of 1×10 −5 s was considered for both cases). OpenFOAM required around 10 h of computational time to simulate 2 s, while MFiX needed around 4 h. The eruption of the first bubbles with both codes are shown in Figures 4 and 5. Here, it can be observed that the hydrodynamics predicted by both codes is clearly different. MFiX predicts a more compact bed with bubbles only produced above the jet, while OpenFOAM predicts small bubbles above the region that is at minimal fluidization conditions, which agrees with the expected behavior for Geldart B particles. In addition, compared to MFiX, OpenFOAM predicts a bigger first bubble above the jet, more splashing of particles once the first bubble erupts and a layer of solids in contact to the wall while the first bubble is moving upwards. This behavior can be observed in Figure 6, which shows the time-averaged solids fraction field for both codes.
The results are compared to the experimental observations and numerical predictions of Yusuf et al. [42]. The numerical results shown here correspond to the same modeling of the thermal conductivity of the phases (as described in Section 2). An argument to explain the differences between the experimental and numerical predictions can be related to the low sampling frequency during the experiment, which might filter the peaks observed numerically. Another reason might be related to the use of a conductivity model that is meant for the bulk of the bed. In any case, differences in the heat transfer predicted by both codes are to be expected given the different flow patterns shown in Figures 4 and 5.   Table 4 shows a time-averaged value of the heat transfer coefficient (between t = 1 and 2 s) where, in spite of the thermal conductivity model adopted, MFiX predicts a heat transfer coefficient that is close to the experimental measures. The simulation of Yusuf et al. shows much higher time-averaged values of the heat transfer, while OpenFOAM results fall in between. These differences might be correlated to the hydrodynamic behavior observed with both codes. It is expected that having a layer of particles in contact to the wall, as predicted by OpenFOAM, increases the effective phase conductivity and, therefore, increases the local heat transfer. On the other hand, although MFiX instantaneous heat transfer predictions do not agree completely with the experiment, the local time-averaged heat transfer coefficient is very similar.

Test 3: Heat Transfer from Submerged Tubes
This test is based on the work of Kim et al. [5]. Experiments were carried out in a 3D fluidized bed (0.48 × 0.6 × 0.34 m). A tube bundle in a triangular arrangement (pitch length 0.08 m), with each tube of 0.34 m length and 25.4 mm outside diameter, is located within the particulate bed (as shown in Figure 8a). A central tube wall is set a constant temperature of 333 K, where a thermal probe is located to evaluate the heat transfer between the tube and the bed. Sand particles are considered for the experiment and the simulations, all the numerical and physical parameters are summarized in Tables 1 and 2. A grid sensitivity analysis is performed a priori for both codes based on a meshconverged fields evaluation (see Table 5). Moreover, the meshing technique of each code is different, meaning it is not possible to evaluate the performance of both codes using the same FVM grid. Nevertheless, results are compared using the coarser refinements for each code upon which the heat transfer coefficient between the hot tube and the bed do not change significantly for a higher level of refinement. For MFiX, a uniform structured grid of 3,133,440 hexahedral cells where the boundary cells are truncated so that they conform to the boundary surface (cut-cell technique), as shown in Figure 8b, is used. For OpenFOAM, the refinement at which the heat transfer coefficient converged to a fixed value consists of three levels of refinement around the tubes with cells of 2 mm in contact with the tubes (as shown in Figure 8c and a maximum cell size of 1 cm far from the tubes bank region. The mesh is generated with snappyHexMesh, and the amount of cells is 330,152. Here, it is important to mention that it is not possible to make a further refinement close to nonplanar surface boundaries with the MFiX mesher. This implies that a uniform refined grid in the whole domain is necessary to accurately predict the field gradients of velocity and temperature near the tubes, which increase the computational costs relative to OpenFOAM. OpenFOAM required around 1 day of overall computational time running in parallel in four CPUs, while MFiX required around 5 days.   Figure 9 shows an instantaneous solids fraction distribution predicted by OpenFOAM and MFiX. Here, OpenFOAM shows a more expanded bed with only a few defined bubbles. MFiX, unlike OpenFOAM, shows clearly defined bubbles with regions of particles at maximum packing. In addition, smaller bubbles appear above the distributor, and larger bubbles move upwards around the tubes bundle. The instantaneous local heat transfer coefficient is computed according to Equation (29). Then, a time-and surface-averaged over the surface of the tube is computed. Figure 10 shows this result for the experiments and with both CFD codes. Here it can be seen that both code seem to moderately overpredict the heat transfer. This might be due to the need of a near-wall effective conductivity model. Moreover, while OpenFOAM seems to follow the general trend of heat transfer as a function of the fluidization velocity, MFiX shows almost no dependence of the heat transfer on velocity. Although a highly refined mesh was used for this problem, it is likely that this problem requires an even higher refinement near the hot tube for MFiX. This issue becomes relevant considering that the MFiX mesher does not allow for a selective refinement near curved surfaces and a uniform highly refined mesh would be necessary to capture the thermal gradients close to the active heat transfer surfaces.

Conclusions
This work analyzes the performance of the open-source CFD codes MFiX and Open-FOAM for predicting heat transfer and minimum fluidization velocities in bubbling fluidized beds. Both codes use the Two-Fluids Model coupled with the KTGF and Frictional theory for the rheological closure and include energy balances for each phase. Expressions for particle-to-fluid heat transfer coefficient and for stagnant thermal conductivity were implemented in OpenFOAM to simulate the thermal behavior.
Values of the minimum fluidization velocity and its dependence on the temperature are appropriately predicted by both codes. Regarding the wall-to-bed heat transfer coefficient estimation, both codes using the same models predict slightly different fluid-dynamic patterns which eventually have an impact on the heat transfer. For the case of the heat transfer from a lateral wall in a pseudo-2D system, MFiX predicts bubbles that erupt above the air jet with little amount of solids within and almost no bubbles in the rest of the bed, while OpenFOAM predicts a much more chaotic fluidization with small bubbles in the width of the bed above the distributor. Compared to MFiX, OpenFOAM predicts a bigger main bubble above the jet with a layer of particles that is in contact to the hot wall most of the time. This different behavior affects the heat transfer prediction since it modifies the instantaneous volumetric distribution of phases and the effective conductivities. The MFiX results are in close agreement with the experimental data while OpenFOAM requires more refinement near the wall to achieve mesh-converged fields and the heat transfer coefficient is overestimated. Nonetheless, these results are closer to the experiment values than that of the simulations made by the authors, using the same physical models. Regarding the heat transfer from a tube in an immersed tube bundle, both codes seem to overpredict the time-averaged heat transfer coefficient for different fluidization velocities. This is in agreement with the results of Test 2 and suggests the need of a near-wall conductivity model. Nonetheless, OpenFOAM predicts the same trend of the experimental observations of heat transfer for a superficial gas velocity value around 1.25 times of minimum fluidization velocity. In this regard, MFiX is not able to reproduce the same trend. Moreover, the meshing technique available in MFiX does not allow for a selective refinement close to the curved surfaces (i.e., tube wall) where large thermal gradients arise which translates into a need of a very large amount of cells in total to be able to capture the high gradients around the hot tubes. This is not the case of OpenFOAM, which can reproduce similar results with far fewer cells than MFiX, therefore saving a lot of computational cost.
In general, both codes are able to predict global hydrodynamic patterns in fluidized beds and how they are influenced by thermal effects. Regarding the CFD predictions in problems involving wall heat transfer, MFiX results, compared to OpenFOAM, present a high level of accuracy with the experimental data for simple geometries involving planar boundaries. However, for the simulation of large-scale systems with nonplanar walls with heat transfer (such as tube bundles immersed in a fludized medium where high thermal gradients are expected), MFiX becomes hindered by its own meshing tool by not allowing a relative refinement in the domain. It should be borne in mind that MFiX was originally conceived for these types of applications, involving multiphase flow in fluidized conditions, while OpenFOAM is a general CFD multipurpose platform with a much broader scope for fluid dynamics. In the context of the present applications, it is expected that future developments in MFiX will be directed to the meshing tools, while in OpenFOAM, the efforts should be focused on the accuracy of the hydrodynamics in bubbling fluidization.   Prandtl number (Pr = µ g c p,g /κ g,0 ) s Subindex for solid phase g Subindex for gas phase