The mesoscale model for gas–liquid flow based on the energy minimization multi-scale method, which is denoted by the GL-EMMS model, demonstrates its capability of predicting the regime transition in bubble columns. Since the dominant mechanisms and stability conditions are important for this mesoscale model, it is desirable to reveal more details through the computational fluid dynamics (CFD) work. This work proposes a gradual drag correlation based on the GL-EMMS model, and then incorporates it into the Euler–Euler (E–E) and Euler–Lagrange (E–L) frameworks respectively. Cases with different gas distributors are simulated via the opensource platform of OpenFOAM. A simplified meshing treatment of the inlet boundary is adopted and has proven effective to reducing the computing load. The E–E simulation predicts the total gas holdup and the radial distribution of liquid velocity well, but obvious deviation is observed in the radial distribution of gas holdup, especially at lower superficial gas velocities for the single-orifice aeration. Nonetheless, this shortage can be overcome saliently by the E–L method. These verified simulations enable further possible iterative investigation between the GL-EMMS model and the CFD work, and also inspire some meaningful simulating strategy for bubble columns.