Summary: Simulation of nano-scale channel flows using a coupled Navier-Stokes/Molecular Dynamics (MD) method is presented. The flow cases serve as examples of the application of a multi-physics computational framework put forward in this work. The framework employs a set of (partially) overlapping sub-domains in which different levels of physical modelling are used to describe the flow. This way, numerical simulations based on the Navier-Stokes equations can be extended to flows in which the continuum and/or Newtonian flow assumptions break down in regions of the domain, by locally increasing the level of detail in the model. Then, the use of multiple levels of physical modelling can reduce the overall computational cost for a given level of fidelity. The present work describes the structure of a parallel computational framework for such simulations, including details of a Navier-Stokes/MD coupling, the convergence behaviour of coupled simulations as well as the parallel implementation. For the cases considered here, micro-scale MD problems are constructed to provide viscous stresses for the Navier-Stokes equations. The first problem is the planar Poiseuille flow, for which the viscous fluxes on each cell face in the finite-volume discretization are evaluated using MD. The second example deals with fully developed three-dimensional channel flow, with molecular level modelling of the shear stresses in a group of cells in the domain corners. An important aspect in using shear stresses evaluated with MD in Navier-Stokes simulations is the scatter in the data due to the sampling of a finite ensemble over a limited interval. In the coupled simulations, this prevents the convergence of the system in terms of the reduction of the norm of the residual vector of the finite-volume discretization of the macro-domain. Solutions to this problem are discussed in the present work, along with an analysis of the effect of number of realizations and sample duration. The averaging of the apparent viscosity for each cell face, i.e. the ratio of the shear stress predicted from MD and the imposed velocity gradient, over a number of macro-scale time steps is shown to be a simple but effective method to reach a good level of convergence of the coupled system. Finally, the parallel efficiency of the developed method is demonstrated.