  The accurate numerical simulation of high Reynolds number incompressible
flows is a challenging topic in computational fluid dynamics. Classical inf-sup
stable methods like the Taylor-Hood element or only $L^2$-conforming
discontinuous Galerkin (DG) methods relax the divergence constraint in the
variational formulation. However, unlike divergence-free methods, this
relaxation leads to a pressure-dependent contribution in the velocity error
which is proportional to the inverse of the viscosity, thus resulting in
methods that lack pressure robustness and have difficulties in preserving
structures at high Reynolds numbers. The present paper addresses the
discretization of the incompressible Navier-Stokes equations with high-order DG
methods in the framework of projection methods. The major focus in this article
is threefold: i) We present a novel postprocessing technique in the projection
step of the splitting scheme that reconstructs the Helmholtz flux in
$H(\text{div})$. In contrast to the previously introduced $H(\text{div})$
postprocessing technique, the resulting velocity field is pointwise
divergence-free in addition to satisfying the discrete continuity equation. ii)
Based on this Helmholtz flux $H(\text{div})$ reconstruction, we obtain a high
order in space, pressure robust splitting scheme as numerical experiments in
this paper demonstrate. iii) With this pressure robust splitting scheme, we
demonstrate that a robust DG method for underresolved turbulent incompressible
flows can be realized.
