In this thesis, turbulent Bose-Einstein condensates are modeled and simulated based on solutions to the semi-classical Gross-Pitaevskii equation, using a mean-field description. The initial configuration consists of a Gaussian free gas distribution, which is then propagated in imaginary time, under the influence of a harmonic trapping potential, to reach a ground state. During this process, turbulence is induced in the system by phase imprinting. The condensate is left to freely expand, from a variety of different turbulent states, in an expanding coordinate system. The time steps are iterated using a 4th-order Runge-Kutta scheme and the spatial derivatives are evaluated using finite differences. Motivated by the recent experiments of Bagnato et al., a variety of different turbulent state variations are considered, in order to reproduce the experimentally observed self-similar free expansion of turbulent ellipsoidal condensates; including variations in vortex orientation, separation and density. The experimental condensates are modeled on anisotropic two-dimensional condensates with equally-oriented fcc vortex lattices. The results show that the rotation of the bulk condensate, the rate at which the vortices tend to an Abrikosov lattice configuration and the self-similarity of the free expansion, all increase with vortex density. Extrapolations of these trends agree with the experimental results.