The Euler-Euler model is widely used in bubbly flow simulations up to industrial dimensions. The standard Euler-Euler model is based on the phase-averaging method. After averaging, the bubble forces in the field equations are functions of the local void fraction. In simulations, when the bubble diameter is larger than the computational cell spacing, the forces can transport the gas belonging to the same bubble in different directions. By contrast, a closure model for the bubble force is typically developed based on the assumption that the force is a resultant force that acts on the bubble’s center of mass. This inconsistency can lead to a nonphysical gas concentration in the center of a channel or near the channel wall if the bubble diameter is larger than the cell spacing. The purpose of the present contribution is to develop an Euler-Euler model where the bubble force consistency is recovered for two-phase flow simulations where the diameter of the disperse phase can be larger than the cell spacing. Such an Euler-Euler model is developed by combining an existing particle-center-averaged Euler-Euler framework with a Gaussian convolution method. To validate this Euler-Euler approach, a comparison is made with experimental data for the bubbly flows in two different vertical pipes. The results show that the proposed Euler-Euler model recovers the bubble force consistency and alleviates the overprediction of the void fraction peak near the wall, while its simulation results in the axial gas and liquid velocity and the liquid turbulence kinetic energy are similar to the results of the standard Euler-Euler model.