We validate and test two algorithms for the time integration of the Boussinesq form of the Navier—Stokes equations within the Large Eddy Simulation (LES) methodology for turbulent flows. The algorithms are implemented in the OpenFOAM framework. From one side, we have implemented an energy-conserving incremental-pressure Runge–Kutta (RK4) projection method for the solution of the Navier–Stokes equations together with a dynamic Lagrangian mixed model for momentum and scalar subgrid-scale (SGS) fluxes; from the other side we revisit the PISO algorithm present in OpenFOAM (pisoFoam) in conjunction with the dynamic eddy-viscosity model for SGS momentum fluxes and a Reynolds Analogy for the scalar SGS fluxes, and used for the study of turbulent channel flows and buoyancy-driven flows. In both cases the validity of the anisotropic filter function, suited for non-homogeneous hexahedral meshes, has been studied and proven to be useful for industrial LES. Preliminary tests on energy-conservation properties of the algorithms studied (without the inclusion of the subgrid-scale models) show the superiority of RK4 over pisoFoam, which exhibits dissipative features. We carried out additional tests for wall-bounded channel flow and for Rayleigh–Bènard convection in the turbulent regime, by running LES using both algorithms. Results show the RK4 algorithm together with the dynamic Lagrangian mixed model gives better results in the cases analyzed for both first- and second-order statistics. On the other hand, the dissipative features of pisoFoam detected in the previous tests reflect in a less accurate evaluation of the statistics of the turbulent field, although the presence of the subgrid-scale model improves the quality of the results compared to a correspondent coarse direct numerical simulation. In case of Rayleigh–Bénard convection, the results of pisoFoam improve with increasing values of Rayleigh number, and this may be attributed to the Reynolds Analogy used for the subgrid-scale temperature fluxes. Finally, we point out that the present analysis holds for hexahedral meshes. More research is need for extension of the methods proposed to general unstructured grids.