Due to the practical importance of boiling heat transfer, various attempts have been made to date to formulate algorithms for computer simulations of heat transfer between vapor films and/or bubbles departing from heated walls and the surrounding liquid. However, due to problem complexity, the majority of testing, verification and applications of computational models have been limited to relatively simple problems, typically 1D or 2D. The purpose of this paper is to present the results of a study on the development of a novel approach to the complete three-dimensional (3D) modeling of simultaneous evaporation and condensation occurring at subcooled boiling flow conditions in heated channels. The focus of the proposed method has been on a multiple-option model of interfacial heat transfer, which allows for an accurate evaluation of the interfacial heat flux between vapor bubbles and the surrounding subcooled or superheated liquid. The modeling of a variety of heating/cooling situations has been made possible by using appropriate newly-developed expressions for a local, position- and bubble-shape-dependent, Nusselt Number at the bubble/liquid interface. To assure the required accuracy, the new heat transfer model was extensively tested and verified in a stand-alone fashion before it has been numerically coupled with a Level-Set (LS) type solver and implemented in the NPHASE/CMFD computer code. The complete model has been tested again and experimentally validated. The stand-alone tests included model verification against the theoretical solution of the Stefan problem, and the near-analytical solution to the model of time-dependent condensation rate around a 2D bubble. As a starting point to a three-dimensional analysis, a DNS-type model has been developed to evaluate the local Nusselt Number for subcooled/superheated liquid flow around bubbles of different (fixed) shapes. The results have been validated against the available existing data. The formulae for the local Nusselt Number which have been deduced as a result of the follow-up analysis have then been implemented in the overall model of evaporation/condensation heat transfer around flowing bubbles of varying shapes. The results of NPHASE/CMFD-based 3D simulations have been validated against experimental data for a single bubble, and a good agreement has been obtained. Subsequently, simulations for multiple bubbles at subcooled boiling flow conditions have been conducted to illustrate the wide range of applicability of the proposed approach, and to demonstrate the ability of the current computational framework to model micro-scale mechanisms behind complicated nonequilibrium boiling processes. It is envisioned that future uses of the first-principle methodology discussed in this paper as virtual experiments should facilitate the development and verification of various multiphase closure laws for RANS-scale models for a variety of boiling heat transfer situations. (C) 2018 Elsevier Ltd. All rights reserved.