This paper treats topology optimization of natural convection problems. A simplified model is suggested to describe the flow of an incompressible fluid in steady state conditions, similar to Darcy's law for fluid flow in porous media. The equations for the fluid flow are 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 a two-material interpolation scheme is applied to both the permeability and conductivity of the distributed material. Due to the simplified model, the proposed methodology allows for a significant reduction of the computational effort required in the optimization. At the same time, it is significantly more accurate than even simpler models that rely on convection boundary conditions based 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 formally compared with results of previous work obtained from solving the full set of Navier-Stokes equations. The results are compared in terms of performance of the optimized designs 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 generating an optimized design during the workday on a small computational cluster and overnight on a high-end desktop.