For simulation of gas?solid bubbling fluidized bed, the bubble-based energy minimization multi-scale (EMMS) model was extended to the sub-grid level through a two-step scheme. The modified EMMS drag is not only related to the voidage, but also associated with velocities within a computing grid. For a specified operating condition, the emulsion voidage was obtained under the optimization of a macro-level stability condition. Then the structure parameters and corresponding drag coefficient were solved by a set of conservation equations with the emulsion voidage correlation at the sub-grid level. Finally, this new drag was coupled with the two-fluid model (TFM) through a user-defined function (UDF) to simulate bubbling fluidized bed with Geldart A particles. It was shown that this new drag could well describe heterogeneous flows of bubbling fluidized beds and reasonably predict the axial profile of average solid concentration with a dilute-top and dense-bottom distribution. The values of axial average solid concentration are close to the experimental data with a relative deviation of 6.4%. Obviously, the radial solid concentration has higher value near the wall than at the center, showing good agreement with experiment at different heights (0.6, 0.8 and 1.1 m) except at the bottom (0.4 m) with a relative deviation of less than 14%. By comparison, the Gidaspow drag model underestimates radial solid concentration with a relative deviation of 91% at four bed heights.