A three-dimensional modified discrete element method (mDEM) in which the effective particle surface for the hydrodynamic drag force and the disturbance of neighboring particles on the flow field are taken into account is proposed to simulate the deformation and breakup process of large aggregates in flows. First, the dynamic shape factor of rectangular aggregates is simulated and compared with the corresponding experiment, showing that the method can predict the behavior of aggregates in fluids quantitatively. Secondly the method is applied to simulate the breakup of large particle-cluster and cluster-cluster aggregates in shear and elongational flows. It is found that a power law relation holds between the average number of particles in broken fragments and the intensity of flow field. In the case of simple shear flow, the value of is approximated by the following universal function of N-DA the ratio of the representative hydrodynamic drag force and the adhesive force, independently of the number, the size and the size distribution of constituent particles, and the minimum separation distance between particle surfaces. This correlation agrees well with the experimental results reported elsewhere. = 27.9 x N-DA(-0.872) It is also predicted that elongational flow is more effective to break up aggregates than the simple shear flow under usual flow conditions.