Abstract A novel Two-way coupled Euler-Lagrange approach, including bubble-bubble collision, coalescence, with variable bubble radius, and bubbles breakup was applied to simulate air-water bubbly flow in vertical pipes. This approach uses the Continuous Random Walk CRW models for creating the velocity fluctuations according to the given state of turbulent kinetic energy and Dissipation rate at the location of the bubbles which is solved by the k-e turbulence model. This dissertation i) describes the development of the Euler-Lagrange Approach under study, ii) present a study for the two-way coupling effect on both continuous and dispersed phases properties, iii) study the effect of both lift force coefficient and Bubble Induced turbulence on the gas void fraction distributions, iv) study the effect of the bubbles coalescence and breakup on bubbles sized and gas void fraction distributions. And presents the results of the simulations performed under each of these considerations. The two-way coupling process toke the effect of the dispersed phase on the continuous one through inserting source terms in the conservation equations of momentum and Turbulence. And also modify the volume of the computational cell in the Euler solver available for the continuous phase according to the void fraction of the cell. Due to the two-way coupling process some studies needed to be performed like the adjustment of the lift force coefficient and the BIT relation due to the change of the liquid velocity profiles due to the two-way coupling. The bubble-bubble collision was applied in the two-way coupling process. It was found that the effect of considering the collision effect on the void fraction distribution only in the high velocity and high gas holdup cases with small effect. Also the bubble-bubble coalescence was applies as a complementary part of the collision process using the film drainage model of Chesters (1991) that compares the film drainage time with the contact time to calculate the coalescence efficiency. The coalescence model was tested first before applying the breakup, so the bubbles size increased only and this affected on the void fraction distribution negatively. Then the breakup model of Martinez Bazan 1991 was applied to perform the equilibrium in the bubbles breakup and coalescence process. These two final processes of the coalescence and breakup found to consume much time, the reason that did not give us chance for testing many cases with considering both coalescence and breakup. The main investigation point through the development of this work was the BIT term and its effect in the used CRW. This was studied in nearly every phase of the model development to study the effect of the BIT under the different considerations of the bubbles interaction mechanisms. Finally we could conclude a final expression for the BIT that is used successfully with the CRW models under consideration.