Acoustic wave propagation in polydisperse bubbly liquids is relevant to diverse applications, such as ship propellers, underwater explosions, and biomedical applications. The simulation of bubbly liquids can date back to Foldy who presented a general theory. In the linear regime, two frequently used models for bubbly liquids are based on the continuum theory and on the multiple scattering theory. Under the homogenization-based assumption, models based on the volume-averaged equations or on the ensemble-averaged equations are designed to find the solutions of a given two-phase flow. The effective wave number is derived through the linearization of these equations. A second approach to the problem of linear wave propagation utilizes the multiple scattering theory. Bubbles are treated as point-like scatterers, and the total field at any location can be predicted by multiple scattering of scatterers. However, in most of experimental researches, the comparison between the approaches and the experimental results is not satisfactory for frequencies near the peak of phase speed and attenuation. In fact, the discrepancies between measurements and approaches are irregular, and the explanations of these discrepancies need further studying. We indicate that such a discrepancy should be attributed to an implicit assumption in these approaches:the bubbles are spatially uniform distribution and statistically independent of each other. In contrast, the complex bubble structures can be observed in many practical bubbly liquids which have important consequences for the acoustic wave propagation. In this paper, our intent is to model the effect of small bubble cluster on linear-wave propagation in bubbly liquids using the self-consistent method. The quasi-crystal approximation is applied to the self-consistent method, and the effective wave number is derived. According to the experimental results, the small clusters of bubbles often exist in bubbly liquids. Therefore, a three-dimensional random model, the Neyman-Scott point process, is proposed to simulate bubbly liquid with the cluster structure. Using this method, we study the influence of such a phenomenon on acoustic dispersion and attenuation relation. A formula for effective wavenumber in clustered bubbly liquid is derived. Compared with the results from the equation of Commander and ProsperettiJ. Acoust. Soc. Am. 85 732 (1989), our results show that the clustering can suppress peaks in the attenuation and the phase velocity, each of which is a function of frequency. Further, we provide a numerical method. A clustered bubbly liquid is simulated with strict mathematical method and the statistical information is obtained through ratio-unbiased statistical approach. Using such a method, we quantificationally analyze the influence of estimated value on predictions.