A reliable and fast numerical scheme is crucial for the 1D simulation of blood flow. In this paper, a 1D blood flow model is incorporated with a Kelvin-Voigt viscoelastic constitutive relation of wall. This lead to a nonlinear hyperbolic-parabolic system, which is then solved with four numerical schemes, namely: MacCormack, Taylor-Galerkin, second order finite volume and local discontinuous Galerkin. The numerical schemes are tested on an uniform vessel, a simple bifurcation and a network with 55 arteries. In all of the cases, the numerical solutions are checked favorably against analytic, semi-analytic solutions or clinical observations. Among the numerical schemes, comparisons are made in four aspects: the accuracy, the ability to capture shock-like phenomena, the computation speed and the complexity of the implementation. The suitable conditions for the application of each scheme are discussed.