Reynolds-averaged Navier-Stokes (RANS) computations of heat transfer involving wall bounded flows at elevated Prandtl numbers typically suffer from a lack of accuracy and/or increased mesh dependency. This can be often attributed to an improper near-wall turbulence modeling and the deficiency of the wall heat transfer models (based on the so called P-functions) that do not properly account for the variation of the turbulent Prandtl number in the wall proximity (y+< 5). As the conductive sub-layer gets significantly thinner than the viscous velocity sub-layer (for Pr >1), treatment of the thermal buffer layer gains importance as well. Various hybrid strategies utilize blending functions dependent on the molecular Prandtl number, which do not necessarily provide a smooth transition from the viscous/conductive sub-layer to the logarithmic region. This work relies on the k-ζ-f turbulence model and the underlying hybrid wall treatment, which is capable of predicting the near-wall momentum and heat transfer with more fidelity, compared to the standard or low-Re variants of the k-z-ε turbulence model. Based on a new DNS database for turbulent flow and heat transfer in a heated pipe (Reτ=360, Pr=1, 10, and 20), a two-layer wall heat transfer model has been formulated. A priori analysis and RANS predictions of the reference heated pipe flow are encouraging, showing improvements of the near-wall heat transfer predictions with respect to accuracy and mesh independence. The potential of the proposed model in real engineering applications is demonstrated in the cooling of electric/hybrid powertrain components, by simulating heat and fluid flow in the e-motor water jacket model.