The boiling vapor-liquid flow at the shell side of shell-and tube heat exchangers was simulated by Euler-Euler transient 3D method in this paper. The mass and heat transfers between the two-phase fluid and heated wall for the subcooled boiling phenomenon were described by the Rensselaer Polytechnic Institute model (RPI model), while the steam condensation within the subcooled liquid was described by the Lee model. Firstly, different turbulence and interfacial force models were evaluated by comparing with the experimental data of Bartolomej (1982). It was found that the turbulence models have minor influence on the temperature and vapor volume fraction distributions. As the bubble size in the subcooled boiling process is small (usually <1 mm), the velocity slip between the vapor bubbles and the liquid is not so important. The simulation results using different drag force models are similar, and the Tomiyama model offers relatively better predictions. The non-drag forces could not significantly improve the accuracy in our simulations. Then the gas-liquid' boiling flow at the shell side of shell-and-tube heat exchangers was then simulated. It was found that the water temperature increases almost linearly near the inlet zone, and the increase speed was slowed down when the bulk temperature approached to the saturated point as the boiling process happened more frequently and consumed much heat. The heat exchangers with the triangle and square configurations have similar temperature and vapor distributions. Further analyses for those two kinds of tube configurations are needed.