A numerical approach to incorporate the compressibility of bubbles to a two-phase solver with incompressible liquid is presented. The scheme introduces a bubble mass density equation and the ideal gas law (though other thermodynamic relations can also be used) to compute the void fraction for use in the continuous phase equations. Pressure and velocity are strongly coupled on a collocated grid scheme by obtaining face velocities through nonlinear interpolation from nodal values using an algebraic equation. Tight coupling is achieved for pressure and void fraction by careful treatment of the bubble transport and compressibility terms in the resulting pressure Poisson equation. The novel proposed implicit strategy prevents numerical instabilities even at high void fractions and highly compressible bubbles. In addition, the proposed method can capture acoustic waves for void fractions between 0.1% and 90% with errors in speed of sound of less than 2% respect to the speed of sound of isothermal compressible bubbles and water mixtures. 1D, 2D and 3D numerical tests, including the bubbly flow around a surface ship, are performed to demonstrate robustness and range of applications of the proposed scheme. The approach is shown to be stable at high void fractions when other schemes fail. (C) 2018 Elsevier Ltd. All rights reserved.