||With the rapid development of computational power, the status of the numerical simulation has been changing from an engineering tool to a new scientific branch, because the applied scope of the numerical simulation has already been far beyond that of the experiment or theoretical analysis. More and more brand new physical phenomena and theories have been discovered from the numerical simulation. In order to make the numerical simulation as an independent scientific branch, the direct modeling on the discretized space is the basic and natural idea. In this thesis, the direct modelings on the discretized space will be clarified. A new concept, the PDE-based modeling, will be introduced. In the simulation of fluid flow, a direct physical description of the flow evolution on the discretized space is needed. The conservations of mass, momentum and energy in a control volume are the basic physical laws for a fluid system. Then, based on various modeling equations, different schemes basically use different methods to model the flux and source terms in the conservation laws. If a modeling is based on the gas-kinetic equation, we will call it gas-kinetic scheme(GKS). In recent years, high-order numerical methods have been extensively investigated and widely used in computational fluid dynamics(CFD). In this thesis, the state-of-the-art WENO-type initial reconstruction and the gas-kinetic gas evolution model will be used in the construction of a high-order multidimensional GKS for the Navier-Stokes solutions. The spatial and temporal gas evolution is fully coupled in the newly developed high-order methods. In order to distinguish numerical performance due to different flux modeling, the results from the WENO-Godunov and WENO-Steger-Warming methods will be used for comparison. This study demonstrates that, besides high-order initial reconstruction, an accurate gas evolution model or flux function in a high-order scheme is important as well in the capturing of physical solutions. In a real physical flow, the transport, stress deformation, heat conduction, and viscous heating are all coupled in a single gas evolution process. Therefore, it is preferred to develop such a flux model with multi-dimensional wave propagation, and un-splitting treatment of inviscid and viscous terms. In order to capture the physical evolution of a slowly evolving gravitational hydrodynamic system, the numerical scheme needs to be an accurate shock capturing scheme for the description of general time-dependent flow evolution, and to have a well-balanced property to capture the delicate equilibrium balance. In this thesis, based on the gas-kinetic equation a well-balanced gas-kinetic symplecticity-preserving BGK (SP-BGK) scheme is developed. In the construction of such a scheme, the gravitational potential is modeled as a piecewise constant function inside each cell with a potential jump at the cell interface. The design of such a scheme fully uses the energy conservation, Liouville’s theorem, and the symplecticity preserving property of a Hamiltonian flow, which play important roles in the description of particles penetration and reflection from a potential barrier. More importantly, the use of the symplecticity preserving property is crucial in the evaluation of the moments of a post-interaction gas distribution function with a potential jump in terms of the moments of pre-interaction distribution function. The SP-BGK method is a well-balanced shock capturing gas-kinetic scheme for the Navier-Stokes equations under gravitational field.