We present a new and computationally efficient methodology using osmotic ensemble Monte Carlo (OEMC) simulation to calculate chemical potential-concentration curves and the solubility of aqueous electrolytes. The method avoids calculations for the solid phase, incorporating readily available data from thermochemical tables that are based on well-defined reference states. It performs simulations of the aqueous solution at a fixed number of water molecules, pressure, temperature, and specified overall electrolyte chemical potential. Insertion/deletion of ions to/from the system is implemented using fractional ions, which are coupled to the system via a coupling parameter lambda that varies between 0 (no interaction between the fractional ions and the other particles in the system) and 1 (full interaction between the fractional ions and the other particles of the system). Transitions between lambda-states are accepted with a probability following from the osmotic ensemble partition function. Biasing weights associated with the lambda-states are used in order to efficiently realize transitions between them; these are determined by means of the Wang-Landau method. We also propose a novel scaling procedure for lambda, which can be used for both nonpolarizable and polarizable models of aqueous electrolyte systems. The approach is readily extended to involve other solvents, multiple electrolytes, and species complexation reactions. The method is illustrated for NaCl, using SPC/E water and several force field models for NaCl from the literature, and the results are compared with experiment at ambient conditions. Good agreement is obtained for the chemical potential-concentration curve and the solubility prediction is reasonable. Future improvements to the predictions will require improved force field models.