A combined lattice Boltzmann and discrete element approach is proposed for numerical modelling of magnetorheological fluids. In its formulation, the particle dynamics is simulated by the discrete element method, while the fluid field is resolved with the lattice Boltzmann method. The coupling between the fluid and the particles are realized through the hydrodynamic interactions. Procedures for computing magnetic, contact and hydrodynamic forces are discussed in detail. The applicability of the proposed solution procedure is illustrated via a two-stage simulation of a MR fluid problem with four different particle volume fractions. At the first stage, simulations are performed for the particle chain formation upon application of an external magnetic field; and at the second stage, the rheological properties of the MR fluid under different shear loading conditions are investigated with the particle chains established at the first stage as the initial configuration.