A mesh refinement study and validation of two-fluid model closures were carried out for a bubbling fluidized bed application. The mesh refinement study indicates that a higher degree of mesh refinement is required for atmospheric than for pressurized fluidization. Statistical bubble parameters (bubble frequency, mean bubble rise velocity, mean pierced bubble length and mean bubble volume fraction) were evaluated. The simulated statistical bubble quantities are computed from voidage signals derived from the transient multidimensional solution of two-fluid models. The algorithm for computing these quantities is taken directly from the evaluation program treating the measurement signals. To remedy the long simulation times required to obtain acceptable statistical values, a parallel version of the two-fluid model solver was developed, based on a domain decomposition method for distributed memory computers. A number of problems related to the parallelization are investigated. These are optimal treatment of velocity components on multi-block boundaries, frequency of data exchange at multi-block boundaries, local errors at multi-block boundaries and simulation time requirements.