This work investigates hydrogen adsorption onto various graphene flakes such as coronene and coronene-like as suitable models of graphene within the framework of the DFT-B3LYP method. The non-local van der Waals (vdW) density functional (B3LYP-D3) method is used for both structural geometry optimization and total energy estimations. Calculations were carried out for a hydrogen molecule above a coronene surface with both conventional and vdW corrected DFT to investigate how these approaches perform in the case of hydrogen adsorption on a graphene surface. Our first-principles results within the B3LYP-D3/def2-TZVPP model show that hydrogen physisorbs on a coronene surface with an adsorption energy of -5.013 (kJ mol(-1)) which is in good agreement with the experimental value. The influence of the basis set and graphene flake size were also evaluated and the results indicate that these slightly affect the adsorption properties. We found also that it is crucial to use non-local dispersion interactions to get accurate results for hydrogen adsorption on a graphene surface. Furthermore, the co-adsorption of H2 molecules onto the graphene surface was investigated. The results obtained at the B3LYP-D3/def2-TZVP level show that H2 molecules can be physisorbed on both sides of the graphene layer with adsorption properties similar to those for a single surface. Finally, we showed that H2 molecules might be bound to the graphene surface via a bilayer adsorption scheme with weak adsorption energy. Charge population and electron density analysis confirm the weak binding nature of the system under consideration.