The CFD modeling based on the two-fluid model (TFM) was used on studying the novel indirectly heated bubbling fluidized bed steam reformer (IHBFB-SR). This is a collaborative project between Petrogas and TU Delft for testing the advance’s reactor configuration, which indirectly supplies the heat via radiant burner on the center toward the surrounding bed, thereby improving heat transfer efficiency and reduced the losses. The present work aimed to observe the hydrodynamic and heat transfer of the reactor by first employing air as the fluidizing gas and corundum (Geldart B size) as the bed material. The minimum fluidization condition, bubbles development, and voidage profile are the main objectives for the hydrodynamic simulation. In the heat transfer study, the effect of radiative heat transfer, bubbles and voidage profiles, and different radiative models (P1 and DO) on the heat transfer mechanism were examined. The 2D and 3D models were built, and three drag models: Gidaspow, adjusted Syamlal, and EMMS/Bubbling was employed. The simulation results were then compared to the experimental data obtained, such as minimum fluidization flowrate, pressure drop, bed expansion, and temperature profile on a specific flow rate. Initially, a grid independency test was conducted using five different grid sizes. It is concluded that the appropriate grid size for simulating the IHBFB-SR with a bed material particle size of 0.496 mm should be at least 7.5 mm or 15 times the corundum particle size. The present research used 2.5 mm or five times dp. The minimum fluidization obtained was in the range of 14-16 kg/h based on both 2D and 3D simulation. Nonetheless, if primarily refers to the 3D model results, the minimum fluidization condition should lie around 14 kg/h. Three drag models have also been compared. It was found that the adjusted Syamlal gives the closest result compared to the experimental data. Nevertheless, the drag modification, based only on minimum fluidization conditions like modified Syamlal, tends to overpredict the drag coefficient on the entire range of solid volume fraction. There are also no significant differences in the bubbles or voidage profiles among those three drag models. Adjusted Syamlal has a slightly larger bubble while Gidaspow and EMMS bubbling has a bit smaller one. The expected better result of using the EMMS bubbling drag model does not appear to have a considerable impact on the case of Geldart B or larger particles. There was also an underprediction of pressure drop and bed expansion on the simulation. Two major factors were the absence of proper particle shape representation through a sphericity factor and the lack of precise simulation of particle size distribution. Only a perfect rounded sphere of corundum with a sphericity factor of one and one uniform size of the particle was assumed. In the case of heat transfer simulation, bubbles and voidage effect firstly studied. It was found that the increase of the bubbles frequency and size, as represented on the voidage profile, would improve the heat transfer process indicated by the increase of the heat flux. The bubbles’ occurrence on the bed plays a critical role in the mass transfer’s improvement from and to the vicinity of the radiant burner wall, thus increasing heat transfer. In contrast with the voidage, the increase of superficial gas velocity does not directly influence the heat flux. It was also proved that the radiative heat transfer improved the overall heat flux in this bubbling fluidized bed reactor by about 16.11 %. Though it appears small, this contribution fits the range of other proven works, with a similar operating environment and particle size used. There are two different radiative models performed in the simulation: first-order spherical harmonics method (P1) and discrete ordinate method (DO). Both models presented a similar trend, with P1, has a slightly higher magnitude. However, the P1 model shows a peculiar result of having a strange lower temperature lower on the bottom part of the bed. On the contrary, that is not the case for the DO, which is well known for its accuracy but a higher computational cost demand. It is then concluded that for the present setup, the DO model could perform better. Lastly, the overall heat transfer process was investigated. Since no specific experimental data available for validating the heat transfer properties, only the final steady temperature values of five different thermocouples were used. Comparing these two, it was found that the present model did not give satisfactory results due to overprediction of air temperature above the bed and underprediction of bed temperature at the same time. Some improvements are required, as will be presented in the recommendations section.