Discontinuous Galerkin methods have been shown to be reliable for solving convection-dominated problems, including the incompressible Navier-Stokes equations, while also being well-suited for modern high-performance computing hardware. Matrix-free methods, which compute the action of the discontinuous Galerkin operator by evaluating the underlying finite-element integrals, are attractive for large-scale problems due to their reduced memory use and memory traffic. At low to moderate polynomial degrees, evaluation schemes with local dense matrices of tabulated shape functions prove to be efficient in terms of throughput achieved by the operator. To improve node-level performance of the operator on simplex elements, explicit data parallelism is leveraged by utilizing Single Instruction Multiple Data capabilities. Further, replacing the matrix-vector product in the evaluation step with a dense matrix-matrix product renders the operator compute-bound for cubic elements, achieving approximately 60% of the theoretical peak performance. Since unstructured tetrahedral grids generated by mesh generators often involve element orderings with poor data locality, a hierarchical reordering strategy is developed to improve data locality of the vector access. To demonstrate the effectiveness of the approach, it is compared to matrix-based implementations, including cases involving curvilinear elements and a modified integration scheme for the Poisson operator. The matrix-free implementation outperforms the matrix-based approach for quadratic and cubic elements by a factor of two to three. To solve the Navier-Stokes equations, a high-order dual-splitting scheme is employed. The pressure Poisson equation, discretized with discontinuous Galerkin methods as the natural ansatz space for incompressible flows, is preconditioned by a hybrid multigrid scheme. This preconditioner combines auxiliary continuous finite-element spaces, p-refinement, and h-refinement, while employing algebraic multigrid on the coarse level. Within the preconditioner, the matrix-free implementation transitions to a matrix-based approach where deemed more efficient. Finally, the fluid solver is coupled with a scalar transport problem, demonstrating its applicability to biomedical problems.