Academic investigations digging into the methane flow mechanisms at the nanoscale, closely related to development of shale gas reservoirs, had attracted tremendous interest in the past decade. At the same time, a good understanding of the complex essence remains challenging, while the broad theoretical scope, as well as application value, possesses great attraction. In this work, with the help of molecular dynamics methods nested in LAMMPS software, a fundamental framework is established to mimic the nanoconfined fluid flow through realistic organic shale matrix. Denoting evident discrepancy with existed contributions, shale matrix in this work is composed of specific number of kerogen molecules, rather than simple carbon-based nanotube. Recently, promotion efforts have been implemented in the academic community with the use of kerogen molecules, however, gas flow simulations are still lacking, and the pore shape in the current papers is always hypothesized as slit pores. The pore-geometry assumption seriously conflicts with the general observation phenomenon according to the advanced laboratory experiments, such as SEM image, AFM technology, that the organic pores tend to have circular pore geometry. In order to fill the knowledge gap, the circular nanopore with desirable pore size surrounded by kerogen molecules is constructed at first. The organic nanopore with various thermal maturity can be obtained by altering the kerogen molecular type, expecting to achieve more physically and theoretically similar to the realistic shale matrix. After that, methane flow simulation is performed by utilization of non-equilibrium molecular dynamics, the methane density as well as velocity distribution under different displacement pressures are depicted. Furthermore, detailed discussion with respect to the simulation results is provided. Results show that (a) displacement pressure acts as a dominant role affecting methane flow velocity and, however, fails to affect methane density distribution, a behavior mainly controlled by molecular–wall interactions; (b) the velocity distribution feature appears to be in line with the parabolic law under high atmosphere pressure, which can be attributed to small Knudsen number; (c) the simulation time will be prolonged with larger displacement pressure imposed on nanoconfined methane. Accordingly, this work can provide profound basis for accurate evaluation of nanoconfined gas flow behavior through shale matrix.