In Direct Numerical Simulations (DNS), resolving thermal boundary layers around inclusions such as bubbles, droplets, or particles is particularly challenging with affordable grid resolutions, especially at moderate to high Prandtl numbers. In Interface Capturing and Tracking methods, quantities resolved on the Eulerian grid near the interface are often smeared, making it difficult to achieve converged predictions of heat transfer, particularly for the temperature and concentration gradients along the interface normal. To address this, several studies introduced subresolution techniques using analytical profiles to represent boundary layer variations below the grid scale. These methods were later refined to ensure strict energy conservation between grid and subgrid regions, but are restricted to situations with large-scale separation between the modeled boundary layer and the bubble’s shape or velocity profile. A recent alternative, suitable for moderate Prandtl numbers (tested range 1 to 5) and in regimes where analytical solutions fail, resolves a one-dimensional profile in the normal direction. This 1D approach reconstructs tangential temperature variations from the resolved grid to account for all convective and diffusive contributions but may introduce artificial heat-transfer dispersion, motivating the fully three-dimensional methodology developed here. The proposed approach lies between early Boundary-Fitted Methods, and hybrid formulations that overlay a local spherical grid around each fixed particle onto a global Cartesian mesh. To handle deformable bubbles in the local problem, the interface and surrounding domain are discretised in parametric coordinates to enable efficient scalar resolution around bubbles. A steady-state convection–diffusion equation is solved with a fixed interface geometry and prescribed velocity field. This subresolution method can be coupled with a Cartesian grid via time-dependent boundary conditions, aiming to reduce DNS requirements to fewer than 16 cells per bubble diameter. Validation on simplified configurations up to Peclet numbers of 100, from spherical to ellipsoidal shapes, under constant scalar boundary conditions, shows promising agreement with reference solutions, supporting integration into a DNS solver.