This paper treats topology optimization of natural convection problems. A simplified model is suggested to describe the flow of a steady-state incompressible fluid, similar to Darcy’s law for porous media. By neglecting inertia and viscous boundary layers, the flow model is significantly simplified. The fluid flow is coupled to the thermal convection-diffusion equation through the Boussinesq approximation. The coupled non-linear system of equations is discretized with stabilized finite elements and solved in a parallel framework that allows for the optimization of high-resolution three-dimensional problems. A density-based topology optimization approach is used, where the permeability and conductivity of the distributed material is interpolated. Due to the simplified model, the proposed methodology significantly reduces the computational effort required in the optimization. At the same time, it is notably more accurate than even simpler models that rely on Newton’s law of cooling. The methodology discussed herein is applied to the optimization-based design of three-dimensional heat sinks. The final designs are compared with previous work obtained from solving the full set of Navier–Stokes equations, both in terms of design performance and computational cost. The computational time is shown to be decreased to around 5-20% in terms of core-hours, allowing for the possibility of optimizing designs during the workday on a small computational cluster and overnight on a high-end desktop. However, due to the use of a simplified model, the performance of the final designs are evaluated using the full Navier–Stokes equations. This ensures verification of performance, as well as systematic comparison with reference results.