This paper proposes wall function models to simulate the heat transfer around a cylinder in cross flow with an isothermal and rough surface. The selected case has similitudes with aircraft wing icing: the ice roughness shape, height and distribution. Moreover, the flow is somewhat similar to that found on iced airfoil; and the surface is isothermal like when icing. The Reynolds-Averaged Navier-Stokes, turbulence, energy and mass conservation 7-equation system is solved by two Computational Fluid Dynamics (CFD) codes. To represent accurately the effects of roughness on the heat transfer, the present authors had to modify both codes and to propose new thermal wall functions for them. In addition, it was implemented a momentum wall function that is not so common in CFD codes but it is a standard in aircraft icing simulation. Basically, the present paper is focused on presenting the simulation results improvement obtained by the implementation of thermal wall functions that also considered the effect thermal resistance of viscous sub-layer on convective heat transfer. The numerical results were compared to experimental data of heat transfer around cylinders with an isothermally heated surface and equivalent sand grain roughness height of Ks/D = 900 ⋅ 10-⁵. For the selected case, the flow regime was trans-critical at Reynolds numbers Re = ue ⋅ D/v = 2.2 ⋅ 10⁵. Due to significant blockage effects present in the experimental data, the tunnel walls were also meshed and simulated. In sum, the implementation into CFD codes was considered adequate because the results were close to experimental data around the whole cylinder surface, not only along the leading edge or the separated region. The results accuracy was improved when compared with CFD factory original models. However, the results indicate that further works on wall functions and their validation are required before using CFD in wing aircraft icing.