In this paper the use of computational fluid dynamics (CFD) for the calculation of flow fields in bubble columns is explained. The local bubble size distribution is considered with the aid of a simplified balance equation for the average bubble volume in bubbly flow. Models are developed for the rate of bubble break-up and coalescence based on physical principals. The flow fields in cylindrical bubble columns without internals are calculated using the Euler-Euler method. The small and large bubble fraction are considered as pseudocontinuous phases in addition to the liquid phase. The calculated flow fields are characterised by several large scale vortices. The local volume fractions of gas and liquid are very inhomogeneous and highly time dependent. The calculated volume fractions, velocities and bubble size distributions agree well with experimental results for bubble columns up to 0.3 m in diameter.