The numerical investigation of a two-phase bubbly flow inside a long vertical counterflow heat exchanger has been performed in this study. The heat transfer was running from the steam moving inside the tubing to the two-phase bubbly water in the annulus. A two-dimensional steady-state heat transfer and fluid flow model are based on the coupling of a probability density function (PDF) coupled with a k-epsilon turbulence model is determined to be the best computation technique in this large-scale case. The k-epsilon model was used to analyze the carrier fluid (water or steam) and the PDF model to calculate bubble flow. The authors desired to receive practical results of the heat transfer process rather than paying attention to the formation of bubbles, their coalescence, tracking flow regimes, and other turbulence-related features that are useless in industrial applications. The article explored two cases: short (800 m) and long heat exchanger (8000 m). Several thermodynamic parameters were on track: flow velocities, pressure, temperature, and concentration of the bubbles in the radial direction. Einstein's correction equation adjusted the viscosity of the bubbles. Solution converged with non-uniform grid spacing in the longitudinal and radial directions. The results gave satisfactory conclusions about the validity of the obtained data with theoretical and experimental results. This article opens a new way on computational fluid dynamics (CFD) application to a large-scale simulation of long vertical or horizontal pipelines and coaxial heat exchangers with the coupled bubbly flow and heat transfer applicable to geothermal and petroleum industry problems.