A bivariate population balance equation applied to a grinding process is implemented in a model (PBM). The particles are simultaneously characterized by their size and their mechanical strength, expressed here by the minimum energy needed to break them. PBM is solved by the Direct Quadrature Method of Moments (DQMOM). The mixed moments of the distribution are expressed by the quadrature form of the population density defined for one order (N) and incorporating the weights and the abscissas defined for the two properties. The effect of the quadrature order (N = 2,3,4) and the selected set of the 3N moments needed to solve the system on the accuracy of the results is discussed. For a given order of the quadrature, the selected set of the initial mixed moments slightly affects first the weights and abscissas derived from the initial particle distribution. The set of moments also affects the precision of the moments calculated versus time but only those having high orders in relation with the respective range of the solid properties considered. Problems of convergence and significant differences in the predicted mixed moments are also observed when the order of the quadrature is equal to 2. However, the changes of a bivariate distribution versus time applied to a grinding process are well predicted using the DQMOM approach, choosing a number of nodes equal to 3, associated with a smart selection of the moment set, incorporating all the moments of interest. (C) 2016 Elsevier B.V. All rights reserved.