Particulate systems are usually modelled by resorting to monovariate population balance equations. In fact, when particles can be described as regular equidimensional objects, such as spheres or cubes, one internal coordinate is enough to describe the evolution of particle size and morphology. In these cases standard numerical methods can be used [1]. However, in many practical cases two or more internal coordinates are needed. For example, when dealing with fractal-like particles, the use of only one internal coordinate is not effective in distinguishing between a compact particle (with fractal dimension close to three) and a ramified one (low fractal dimension). For these cases two internal coordinates must be used, for example particle mass and area, or particle mass and fractal dimension. Even if particles are compact two internal coordinates may be required. For example during precipitation because of aggregation or secondary nucleation, particles are often polycrystalline, namely constituted by more than one crystallite. It is clear that if the total particle mass is taken as the only internal coordinate one would not be able to distinguish between a particle constituted by a few big crystallites and a lot of small ones resulting in two crystals with the same total mass. Also for this case a number density function based on the number of crystallites per particle and on the size of a single crystallite should be used. A convenient way to solve the population balance equation is by using the Quadrature Method of Moments (QMOM) [2] which is based on the simple idea of solving transport equations for the moments of the number density function and it uses a quadrature approximation with N nodes to close the unknown terms. In this work the Direct Quadrature Method of Moments (DQMOM) [2], an improvement of QMOM, is used to solve bivariate population balance equations. This method is again based on a quadrature approximation but it tracks the quadrature approximation directly and can be easily used to solve bivariate population equations. In order to asses its accuracy and stability, DQMOM has been used to describe the evolution of a particulate system undergoing simultaneous nucleation, growth, aggregation, breakage and sintering and its predictions have been compared with those obtained with more sophisticated methods such as Monte Carlo methods. Comparison was carried out in terms of some mixed moments that have particular physical meaning. For example, when the two internal coordinates are particle mass and particle surface area, the mixed moment m10 is the total particle mass, whereas m01 is the total particle area and m00 is the total particle number density. In Fig. 1 DQMOM predictions obtained with a quadrature approximation with three nodes are compared with Monte Carlo simulations for two mixed moments, namely m00 and m02. Comparison is reported for a population of particles undergoing pure coalescence (see Fig. 1a,b) and undergoing aggregation and simultaneous restructuring (see Fig. 1c,d). As it is possible to see the agreement between the very accurate but CPU intensive Monte Carlo method and DQMOM is excellent. Results obtained working under different operating conditions seem to confirm that DQMOM is characterized by high accuracy and by surprisingly low computational costs. In fact, for most of the tested cases, the solution of bivariate population balance equations with DQMOM can be done by tracking from six to nine moments of the number density function, or in other words by resorting to a quadrature approximation with two or three nodes, with a drastic reduction of the CPU time in comparison with Monte Carlo methods.

[1] Ramkrishna D., Population balances, New York: Academic Press (2000). [2] Marchisio D.L., Vigil R.D. and Fox R.O., J. Coll. & Int. Sci., 258, 322-334 (2003). [3] Marchisio D.L. and Fox R.O., J. Aerosol Sci., 36, 43-73 (2005).

Fig. 1. Comparison between Monte Carlo simulations and DQMOM with three nodes, as a function of the dimensionless time. Coalescence: a) m00(t)/m00(t=0); b) m02(t)/m02(t=0). Aggregation and restructuring: c) m00(t)/m00(t=0); d) m02(t)/m02(t=0)