We consider discretization-based homogenization of micro-heterogeneous stationary heat conduction problems within the wider framework of Variationally Multiscale methods. The solution procedure for the fine-scale fluctuations is based on the local mesh-size in the macroscopic finite element problem. This allows for bridging the extremes -- from single scale analysis to computational homogenization -- seamlessly within one finite element simulation. One crucial component in two-scale analysis, often denoted "Finite Element squared" (FE${}^2$) is the prolongation of the macroscopic fields onto the subscale problem. In computational homogenization on Statistical Volume Elements (SVEs), it is well known that temperature (Dirichlet) and flux (Neumann) boundary conditions result in upper and lower bounds, respectively, on the effective heat conductivity. However, in the standard (primal) variational format, it is only the Dirichlet boundary conditions that satisfy the conformity requirement in the standard (primal) variational form. Hence, other choices of boundary conditions (e.g. flux or periodic boundary conditions) are not directly applicable if one wants to consider connected adjacent SVE's without assuming (infinite) separation of scales. In order to allow for flux boundary conditions, while satisfying the conformity requirement, we propose a novel mixed dual variational format. In this format, the finite element multiscale problem, sometimes denoted "Finite Element squared" (FE2), (i) incorporates the flux boundary conditions when computational homogenization is adopted, and (ii) is a conforming approximation when the subscale fluctuation is fully resolved within each macroscopic element. We are thus able to compute the FE approximation using either the (standard) primal format and temperature boundary conditions on SVE's or the dual format and flux boundary conditions on SVE's while allowing for the transition from homogenization to single-scale analysis. We compare the numerical performance of the two different variational settings (and consequent boundary conditions) in the FE2 context with different choices of discretization on the macro scale and choices of "SVE-quadrature". In particular, we compare the convergence properties for FE-refinement on both "scales" and discuss the possibilities to utilize the two approximations a strategy to control the approximation error.