The velocity-vorticity formulation for the solution of unsteady three-dimensional incompressible flows by Kim, Moin and Moser (JFM 1987) is extended for low Mach number combustion. A key advantage of the method is the elimination of the pressure from the momentum equations and, as a result, the errors and complications of pressure boundary conditions inherent in pressure-splitting algorithms common in primitive variable-based CFD solvers. Another merit of the proposed formulation is its efficiency for horizontally homogeneous flows where Fourier expansions can be used for discretization of horizontal derivatives. The resulting elliptic system comprises two evolution equations for two dependent variables and two Poisson type equations. Spatial discretization in the vertical is performed using sixth-order accurate compact finite difference and time stepping is carried out using third order implicit-explicit Runge-Kutta scheme. Open boundary conditions are used at the top boundary and free-slip, no-flux conditions are employed at the bottom boundary. Test cases involving the simulation of compressible Taylor Green vortex flows using a parallel MPI based FORTRAN code are briefly discussed and future research directions are summarized.