This thesis reports research carried out towards numerical simulation of hypersonic flows around space-plane shapes. For high speed flows around conical geometries, a locally conical approximation was introduced, which resulted in locally conical Navier-Stokes equations. In order to achieve accuracy and efficiency for steady state solutions, various methods were investigated. Based on the MacCormack implicit scheme and the Beam-Warming implicit scheme, two implicit procedures were developed to solve the locally conical Navier-Stokes equations (LCNSE). A new implicit boundary treatment was introduced in the MacCormack implicit scheme. The source term in the governing equations was treated explicitly. A simplified Beam-Warming implicit scheme was developed for its application to the LCNSE. Accuracy of the two schemes was investigated. The time step dependence of steady state solution with MacCormack-type schemes was analyzed and a procedure to reduce the error was proposed. To further accelerate the convergence to the steady state, two multigrid methods were applied to the two implicit schemes respectively. An extention of Ni-type multigrid method was developed to accelerate the MacCormack implicit scheme, and the FAS multigrid method was employed to accelerate the simplified Beam-Warming implicit scheme. In parallel, a new approach for fast steady state solution - sparse quasi-Newton method - was proposed to avoid difficulties in linearization associated with implicit schemes for general CFD problems. Formulation was given for three-point and five-point spatial discretization schemes. Preliminary results of a nozzle problem with van Leer's flux splitting and Harten's TVD high shock-resolution schemes illustrated significantly faster convergence to steady state with the sparse quasi-Newton approach than those with corresponding implicit operators of van Leer and Harten. Numerical simulations by solving LCNSE with the two implicit schemes developed in this study were carried out on hypersonic flows around a cone, on the leeside of a delta wing and beneath/over a cone-delta-wing combination. Detailed structures of the complex flow interaction were well predicted including the existence of embedded shock waves and secondary vortices. Comparison with available experimental data was made. Euler solutions were also carried out to compare with the N-S solutions. In the present hypersonic delta wing flow simulation, different phenomena were found than would have been expected from the Miller and Wood classification in the lower speed range. The numerical simulation of hypersonic viscous flows around a cone-delta-wing combination was the first flow field simulation around such a shape representing wing-body interference. It was found that the complexity of the flow field results from the shock-shock, shock-boundary layer and shock-vortex interactions in the flow field. High local heating and its cause were revealed near the corner on both the windward side and the leeward side surfaces of the geometry.