A rigorous model is presented for simulating the gas-particle two phase flow and the involved chemical reaction in bubbling fluidized beds (BFB). The present model employs the recently developed (c(2)) over bar-epsilon(c) closures for the turbulent mass transfer equation, thus enables the prediction of species concentration distributions without relying on empirical turbulent Schmidt number. The Two-Fluid model (TMF) is used to represent the fluid mechanics and the coupled chemical reaction. In mathematical expressions, the k(g) - epsilon(g) - k(p) - epsilon(p) - Theta equations for the description of phases turbulence is adopted. To test the validity of the present model, simulations are conducted for reactive multiphase flow in a BFB. The model predictions show good agreement with the experimental data in literatures. Part II of these serial works will further extend the application of this model to simulating a gas combustion process with strong heat effect in BFB. (C) 2018 Elsevier B.V. All rights reserved.