Towards direct numerical simulations of low-Mach number turbulent reacting and two-phase flows using immersed boundaries. Most reacting and two-phase flows of practical interest are turbulent but take place at low Mach numbers or under incompressible conditions. In order to investigate the properties of such complex flows with high accuracy but acceptable computing times, a suitable tool for direct numerical simulations (DNS), called DINOSOARS, has been developed. The present article describes the numerical components and methods implemented in this code, together with a detailed verification and validation phase, and finishes with two examples of full-scale simulations. We hope it might be useful as a “verification and validation guideline” for other researchers working on DNS of reacting flows. Since applications of growing complexity are considered by DNS, a direct boundary immersed boundary method (DB-IBM) has been implemented, allowing a description of arbitrary geometries on a fixed, but possibly refined, Cartesian mesh. A direct force IBM is implemented as well in DINOSOARS in order to resolve large moving spherical particles (much larger than the Kolmogorov scale) on the grid. Particles below the Kolmogorov scale are treated as point particles, taking into account additionally heat and mass transfer with the continuous flow. The efficient parallelization of the code relies on the open-source library 2DECOMP&FFT. The underlying Poisson equation is solved in a fast and accurate manner by FFT, even for non-periodic boundary conditions. The flexibility of DINOSOARS makes it a very promising tool for analyzing a variety of problems and applications involving turbulent reacting and/or two-phase flows.