A two way coupled Immersed boundary method (IBM) and Discrete element method (DEM) is performed in this work to access the behavior of particles in Magnetorheological fluids (MRFs) under the influence of magnetic field. Particle dynamics under linear and non-linear magnetic field is analyzed in different conditions. The results are compared with analytical solution for the case of one-dimensional particle motion. Particle chain formation for mono-dispersion and poly-dispersion (having three different sizes) is simulated when particles are placed in lower (0.0375T) and higher (0.3T) density magnetic fields. Particle placed in a stationary fluid column in non-linear magnetic field, along with gravitational force acting on the particle, is also simulated at the end. A modified version of Foam extend (for IBM), is used to compute the hydrodynamic stress which is transferred to the DEM part. The particle position and velocity are then numerically accessed using in-house developed XDEM code (Extended DEM). The new position and velocity of the particles are transferred back to IBM, and this keeps repeating till the final time step. A module is added to the calculate the effect of magnetic force on the particles in XDEM which is then validated with a benchmark case. For particle under linear magnetic field, the chain formation time for 5 μm size particle diameter and mixture of 5 μm, 4 μm, and 3 μm particle diameters with 10%, 20%, and 30% volume concentration (vc.%) is discussed. It is found that for 0.0375T of field density, the particle chain formation for mono-dispersion with 10 vc.% requires more time while for 30 vc.% it is very less relatively. For poly-dispersion, the particle chain formation depends on the random position of particles and the vc.%. For particle suspended in fluid under non-linear magnetic field, the influence of particle dynamics is sever which can overcome the gravitational force effect.