Using molecular dynamics simulations, the forces acting on nanoparticles inside polymer brushes are computed. Vertical force profiles are obtained under variation of grafting density, nanoparticle (NP) size, solvent quality, and degree of polydispersity. The force profiles are integrated to obtain the inclusion free energies. If the NP size is fixed, this energy scales with the osmotic pressure, consistent with current theoretical models. These models also predict a scaling of the free energy proportional to the volume of the NP, which is verified in good solvent and at high densities. Otherwise, the power exponent remains lower, and surface tension as a possible cause for the observed deviation is discussed. Polydispersity is shown to reduce the inclusion free energy, while the power law scaling as a function of NP size remains unchanged. Finally, polymer, brushes in Theta-solvent are shown to violate simple scaling predictions within the density regime covered in this work. Using a Flory-Huggins mean-field model, we demonstrate that the universal scaling regime is restricted to very low grafting densities below sigma = 0.01, and the observed deviations are a result of higher order contributions to the virial expansion of the osmotic pressure.