An analysis of the two-fluid model in the case of vertical fully developed laminar bubbly flows is conducted. Firstly the phase distribution in the central region of the pipe (where wall effects vanish) is considered. From the model equations an intrinsic length scale L is deduced such that the scaled system reduces to a single equation without parameters. With the aid of this equation some generic properties of the solutions of the model for pipes with diameter greater than about 20L (the usual case, since L is of the order of the bubble radius) are found. We prove that in all physically meaningful solutions an (almost) exact compensation of the applied pressure gradient with the hydrostatic force rho(eff) g occurs (with rho(eff) the effective density and g the gravity). This compensation implies flat void fraction and velocity profiles in the central region not affected by the wall, even when no turbulence effects are accounted for. We then turn to consider the complete problem with a numerical approach, with the effect of the wall dealt via wall forces. The previous mathematical results are confirmed and the near-wall phase distributions and velocity profiles are found. With the numerical code it is also possible to investigate the regime in which the pressure gradient is greater than the weight of the pure liquid, in which case a region of strictly zero void fraction develops surrounding the axis of the pipe (in upward flow of bubbles). Finally, the same code is used to study the effect of reducing the gravity. As g decreases, so does the relative velocity between the phases, making the lift force increasingly dominant. This produces, in upward bubbly flows, narrower and sharper void fraction peaks that also appear closer to the wall. (C) 2003 Elsevier Ltd. All rights reserved.