This paper presents topology optimization-based designs of the heater section of a single-phase rectangular Vertical Heater and Vertical Cooler Natural Circulation Loop (NCL) for enhanced flow rate. A 2-D steady-state, laminar incompressible flow is assumed where Navier–Stokes equation is coupled with convection–diffusion equation incorporating Boussinesq approximation. The Grashof numbers considered in this work span a wide range from \(10^3\) to \(10^8\) . The key challenges highlighted include the difficulties associated with solving velocity-based objectives through density-based topology optimization, especially for natural convection problems. Use of standard interpolation methods such as Rational Approximation of Material properties (RAMP) generates porous outputs for better flow rates, and heavy projection schemes are not much successful due to a wide transition zone of intermediate densities leading to suboptimal results. This work proposes a promising solution to find better designs for enhanced flow rates in NCL by interpolation of the buoyancy force term, which is reported for the first time in literature. The new method of buoyancy interpolation has been tested with two different schemes for material interpolation: RAMP and sigmoid. Numerical experimentation showed that with careful calibration, interpolation of buoyancy force with a sigmoid function in sync with projection scheme enhances the chances of the optimization algorithm to find better optimal output. The significance of buoyancy interpolation in density-based topology optimization of natural convection-based cooling systems such as the NCL at high Grashof numbers is clearly evident from this work.