A highly efficient numerical scheme has been incorporated to solve a traditional two-fluid model for dense gas-solid flow in large scale fluidized beds. Difficulties associated with the numerical solution and boundary condition enforcement especially in the azimuthally direction and at the axis of the cylindrical domain are discussed in detail. Higher order discretization schemes with deferred correction approach have been implemented for the convection terms. The p-e algorithm (Van der Hoef et al., 2006) has been implemented to deal with densely packed regions in the fluidized bed. A modified SIMPLE algorithm is used to solve the pressure and volume fraction corrections, and a projection method is used to obtain solutions of the momentum e...