A computational method for the simulation of hypersonic nozzles is discussed. The method uses an excluded volume equation of state to represent high-pressure effects in the reservoir. Vibrational energy relaxation and chemical reactions may be simulated with finite-rate models. The Spalart-Allmaras turbulence model with a compressibility correction is used to model the turbulent boundary layer on the nozzle wall. Typical high-pressure nozzle flows require extreme grid stretching near the nozzle surface to resolve the boundary layer. The numerics are designed to work well with cell-aspect ratios as high as 105. A low-dissipation form of Steger-Warming flux vector splitting is modified to compute the fluxes for the non-ideal equation of state. The equations are integrated in time using the implicit Data-Parallel Line Relaxation method that is very effective on highly stretched grids. Non-dimensional time steps of 106, corresponding to 5-10 μsec may be used, resulting in machine zero convergence in approximately 2000 iterations. Comparisons with previous simulations and experimental data are generally good, though there are still concerns with comparisons to some AEDC Tunnel 9 data.