An inconsistency exists in bubble force models used in the standard Euler-Euler simulations. The bubble force models are typically developed by assuming that the forces act on the bubbles' centers of mass. However, in the standard Euler-Euler model, each bubble force is a function of the local gas volume fraction because the phase-averaging method is used. This inconsistency can lead to gas over-concentration in the center or near the wall of a channel when the bubble diameter is larger than the computational cell size. Besides, a mesh-independent solution may not exist in such cases. In addition, the bubble deformation is not fully considered in the standard Euler-Euler model. In this thesis, a particle-center-averaging method is used to represent the bubble forces as forces that act on the bubbles' centers of mass. A particle-center-averaged Euler-Euler approach for bubbly flow simulations is developed by combining the particle-center-averaged Euler-Euler framework with a Gaussian convolution method. The convolution method is used to convert the phase-averaged and the particle-center-averaged quantities. The remediation of the inconsistency in the standard Euler-Euler model by the particle-center-averaging method is demonstrated using a simplified two-dimensional test case. Bubbly flows in different vertical pipes are used to validate the particle-center-averaged Euler-Euler approach. The bubbly flow simulation results for the particle-center-averaged Euler-Euler model and the standard Euler-Euler model are compared with experimental data. For monodisperse simulations, the particle-center-averaging method alleviates the over-predictions of the gas volume fraction peaks for wall-peaking cases and for finely dispersed flow case. Whereas, no improvement is found in the simulated gas volume fraction for center-peaking cases because the over-prediction caused by the inconsistency has been smoothed by the turbulent dispersion. Moreover, the axial gas and liquid velocities simulated with both Euler-Euler models are similar, which proves that the closure models for bubble forces and turbulence are correctly applied in the particle-center-averaged Euler-Euler model. For fixed polydisperse simulations, the particle-center-averaging method can also alleviate the over-prediction of the gas volume fraction peak in the center or near the wall of a pipe. The axial gas velocities simulated with both Euler-Euler models are about the same. Comparisons are also made for the simulation results of bubbly flows in a cylindrical bubble column and the experimental data. The gas volume fractions and the axial gas velocities simulated with both Euler-Euler models almost coincide with each other, which indicates that the sink and source terms for the continuity equations and the degassing boundary are set correctly in the particle-center-averaged Euler-Euler model. An oblate ellipsoidal bubble shape is considered in the particle-center-averaged Euler-Euler simulations by an anisotropic diffusion. The influence of bubble shape on the simulation results of bubbly pipe flows is investigated. The results show that considering the oblate ellipsoidal bubble shape in simulations can further alleviate the over-predictions of the gas volume fraction peaks for wall peaking cases, but it has little influence on the gas volume fractions of center-peaking cases and the axial gas velocities.