In this thesis, I investigate the viscoelastic properties of magnetorheological fluids (MRFs) through a simulation method that combines three established methods: molecular dynamics for the simulation of MRF particles, the Lattice-Boltzmann method for the simulation of the carrier fluid, and the Immersed Boundary method describing the coupling between the MRF particles and the carrier fluid. I first describe the method, rank it among the other methods for MRF simulation, and present the simulation step in detail. Finally, I use the method to describe the viscoelasticity of MRFs in shear mode of operation. I focus on the stress-strain curve, creep, and the dynamics of MRF under shear flow. Some of the numerical results are also compared with experimental results from the literature.