A numerical method is tailored for the solution of a reduced set of model equations developed for the description of reactive flows in chemical reactors. The numerical method synthesis intends to utilize experience both in chemical engineering solving dispersion models containing complex reaction kinetics, phase equilibria and non-ideal thermodynamics, and in fluid dynamics applying classical numerical algorithms constructed for pure flow calculations. The two types of model equations involved in reactive flow simulations reflect very different physics. A modular method is therefore suggested enabling a split between the flow- and chemistry model parts. In this way more optimal solution methods can be adopted for each operator, in contrast to more traditional computational fluid dynamics (CFD) methods where all equations (and operators) are normally solved by the numerical solvers originally intended for pure flow calculations. In addition, emphasis has been put on modularity enabling the same model framework to be used both for the more traditional 1D and 2D dispersion models with simple ID flow calculations, as well as for the more elaborate 2D and 3D reactive flow calculations. In the flow part of the algorithm a fractional-step algorithm is constructed for the simulation of dynamic reactive flows in chemical reactors. The numerical scheme is based on the compressible transport equations in the low-Mach-number limit. The method can handle real gas (and liquid) mixtures with variable density as well as constant density fluids. The velocity field is advanced using a projection scheme, which consists of a partial convection-diffusion update followed by a pressure correction step with an intermediate an-elastic filter. A variable density corrector step is implemented for variable density systems in order to couple the evolution of the density and the velocity fields. In the chemistry part of the model all scalar fields are updated using Strang-type operator-split integration steps that combine several explicit convection and semi-implicit diffusion transport operators with a suitable solver tailored for non-linear sink/source terms. Diffusion terms are discretized with a second order central difference scheme in space. Convection terms are discretized. with a second order TVD scheme in space. The temperature advection is discretized using a second order upwind scheme. The chemical reaction part of the model is discretized by an implicit Euler approximation, and the resulting set of algebraic equations is solved using a Broyden subroutine. The other source terms are solved using an explicit Euler approximation. The performance and behavior of the operator-split scheme are assessed based on simulations of two industrial chemical processes (i.e. the synthesis gas and methanol production processes) performed in multi-tube fixed bed reactors. These processes are important parts of the Statoil methanol plant at Tjeldbergodden in mid Norway. Both ID and 2D pseudo-homogeneous dispersion models (i.e. heat and mass balances) with prescribed velocity, mixture density and total pressure profiles are used describing the reactor performance. The predicted profiles for both the methanol and synthesis gas processes were in accordance with results reported in the literature. An oscillatory radial void fraction distribution was then implemented in the 2D model. It was found that the non-uniform void fraction distribution had no significant effect on the temperature and mole fraction profiles. A ID CFD simulation was pet-formed to evaluate the effect of the variations in velocity, pressure and mixture d The changes in these variables significantly effect the composition and conversion in the reactor. A 2D CFD model was then developed for further studies of the multi-tube fixed bed reactor for the synthesis gas process, analyzing the influence of non-uniform void fraction distributions on the flow and chemical conversion. The non-uniform void fraction distribution induces a significant reduction in axial pressure drop and a much higher fluid velocity close to the wall. However, the influence on the temperature and mole fraction profiles was hardly noticeable. Mass- and enthalpy budgets were implemented for the heat and mass quantities in the program. These budgets show that the enthalpy and both the mixture and component mass balances are fulfilled in the simulations. The numerical CFD algorithm developed in this paper has been found to be both computationally stable and accurate. Computational efficiency analysis shows that the Poisson solver requires about 75% of the total computational time. The computational time spend on the chemistry part of the model is about 15% of the total time. About 10% of the cost used on the chemistry part is spend on the chemistry solver.