The route to reliable quantum nanoelectronic devices hinges on precise control of the electrostatic environment. For this reason, accurate methods for electrostatic simulations are essential in the design process. The most widespread methods for this purpose are the Thomas-Fermi (TF) approximation, which provides quick approximate results, and the Schrödinger-Poisson (SP) method, which better takes into account quantum mechanical effects. The mentioned methods suffer from relevant shortcomings: the TF method fails to take into account quantum confinement effects that are crucial in heterostructures, while the SP method suffers severe scalability problems. This paper outlines the application of an orbital-free approach inspired by density functional theory. By introducing gradient terms in the kinetic energy functional, our proposed method incorporates corrections to the electronic density due to quantum confinement while it preserves the scalability of a theory that can be expressed as a functional minimization problem. This method offers a new approach to addressing large-scale electrostatic simulations of quantum nanoelectronic devices.