A hybrid molecular–continuum method for unsteady compressible multiscale flows

Abstract

We present an internal-flow multiscale method (‘unsteady-IMM’) for compressible, time-varying/unsteady flow problems in nano-confined high-aspect-ratio geometries. The IMM is a hybrid molecular–continuum method that provides accurate flow predictions at macroscopic scales because local microscopic corrections to the continuum-fluid formulation are generated by spatially and temporally distributed molecular simulations. Exploiting separation in both time and length scales enables orders of magnitude computational savings, far greater than seen in other hybrid methods. We apply the unsteady-IMM to a converging–diverging channel flow problem with various time- and length-scale separations. Comparisons are made with a full molecular simulation wherever possible; the level of accuracy of the hybrid solution is excellent in most cases. We demonstrate that the sensitivity of the accuracy of a solution to the macro–micro time-stepping, as well as the computational speed-up over a full molecular simulation, is dependent on the degree of scale separation that exists in a problem. For the largest channel lengths considered in this paper, a speed-up of six orders of magnitude has been obtained, compared with a notional full molecular simulation.

1. Introduction

From rapid point-of-care health diagnostics to alleviating the global problem of diminishing fresh water supplies, nanofluidics is expected to be an important component of many future technological applications. For example, nanotubes a few nanometres in diameter but millimetres in length could be used in reverse osmosis desalination and other highly selective filtration membranes.

Nanofluidic applications generally have a few factors in common.

1. (i) The flow occurs in nano-confined geometries (e.g. pores, tubes, channels etc.) that are typically of large aspect ratio. The characteristic length in the streamwise direction can be several orders larger than a characteristic length perpendicular to the flow, such as the pore width or diameter. For example, aligned nanotube membranes for desalination (Holt et al. 2006) or carbon-dioxide storage (Mantzalis, Asproulis & Drikakis 2011); nano/micro heat exchangers; lab-on-a-chip; and lubrication/tribology applications (Heo et al. 2005).

2. (ii) The flow is dominantly non-continuum. For example, molecular layering exists throughout the majority of the flow (Somers & Davis 1992), stress and viscosity can no longer be treated as local quantities (Todd, Hansen & Daivis 2008), and velocity slip can be present at surfaces of various hydrophobic and hydrophilic materials (Lee, Leese & Mattia 2012).

3. (iii) Conventional continuum flow models, such as the standard Navier–Stokes–Fourier equations, on their own are no longer applicable. Even with significant modifications, at very high confinement the models are still questionable (Holland et al. 2014).

4. (iv) The flow is highly multiscale. Continuum variables (such as velocity) vary, in the direction of the confinement, on a scale similar to molecular scales, so no scale separation exists between them. However, in the direction of flow there generally does exist separation between the scale over which the continuum variables vary and the microscopic molecular scale.

Molecular dynamics (MD) is currently the best method available to simulate these types of flows accurately. The fluid and solid parts of the problem are modelled using molecular constituents (atoms and molecules), whose behaviour in space and time is described by Newton’s equations of motion and intermolecular potentials derived from quantum mechanics or experiment. The problem with MD, however, is that it is very computationally intensive. Current processing power prohibits a full molecular calculation of systems any larger than a few nanometres for a few nanoseconds of problem time. A hybrid molecular–continuum approach is therefore an attractive alternative: it aims to combine the best attributes of both MD (i.e. molecular detail) and a continuum-fluid formulation (i.e. computational cheapness).

Various hybrid methods that couple continuum with molecular solvers seamlessly within the same simulation have been proposed. These can be broadly categorised into either domain decomposition (DD) (O’Connell & Thompson 1995; Hadjiconstantinou & Patera 1997; Nie et al. 2004) or heterogeneous multiscale methods (HMM) (Ren & E 2005; Yasuda & Yamamoto 2008; Kessler, Oran & Kaplan 2010; Asproulis, Kalweit & Drikakis 2012; Borg, Lockerby & Reese 2013a ). DD methods are applicable to configurations where the domain can be split distinctly into molecular and continuum-fluid regions, with overlapping regions placed at the molecular–continuum interfaces where two-way coupling based on fluxes or state properties occurs. Typically a Schwarz alternating method is used to converge the solutions in the disparate regions iteratively (Hadjiconstantinou & Patera 1997). A review of DD methods for dense fluids can be found in Mohamed & Mohamad (2010). As an alternative, the HMM is more effectively applied to rheological flows, or problems in which a constitutive relation and boundary information do not exist. The HMM places a continuum-fluid solver grid over the whole domain, and microscopic simulations dispersed at the nodes of the computational grid provide the missing information.

For flow configurations of high aspect ratio, both DD and HMM approaches are too computationally expensive, and in fact HMM is much less accurate than a full molecular simulation. This is discussed in more detail in Borg, Lockerby & Reese (2013c ), Patronis et al. (2013) and Patronis & Lockerby (2014). The internal-flow multiscale method (IMM) of Borg et al. (2013c ) was developed instead for these kind of geometries, and for isothermal, incompressible, steady flows. The method consists of iteratively solving a one-dimensional (1D) steady continuity equation; MD simulations are applied at regularly spaced nodes along the streamwise direction of the domain, with heights chosen to match the local geometry of the full channel (see figure 1). Measurements of mass flow rates are extracted from these MD subdomains, which are then used to evaluate macro continuity errors. Based on these errors the equivalent pressure gradients (i.e. body forces) applied to each MD simulation are adjusted. A converged solution is obtained after only a few (typically two or three) iterations. As the fluid and wall conditions are defined from a microscopic (intermolecular) perspective in the micro elements, non-continuum models for slip, density layering, and non-local viscosity and stress are not required in the macro description. The flow rate is the sole parameter measured in the micro solutions, and provides an adequate and sufficient micro-to-macro coupling. However, several problems with this simple approach have become evident, in particular, the fluid compressibility is not incorporated and the method is only applicable to steady flows.

Figure 1. Schematic of the IMM treatment of an arbitrary channel with high aspect ratio in the streamwise direction $s$ : a macro model is used to solve the full domain, and micro refinement is applied at nodes of the macro grid using parallel-wall MD simulations with periodic boundary conditions.

Substantial fluid compressibility was evident in the flow simulations in Borg et al. (2013c ), even for a liquid at extremely low Mach numbers (Gad-El-Hak 2005). In these highly confined flows, viscous-related pressure losses can significantly compress the fluid. Patronis et al. (2013) include density as a coupling variable for rarefied gas flows, and Patronis & Lockerby (2014) have been able to deal with temperature variations.

In this present paper we describe recent advances towards an extended IMM framework capable of simulating, for the first time, time-dependent and compressible flow problems accurately and efficiently. In these flow problems there can be both spatial and temporal scale separation; our method, termed ‘unsteady-IMM’, can exploit both and thereby achieve a multiplicative effect on the computational savings. The temporal part of the unsteady-IMM is posed in the general time-stepping framework described by Lockerby et al. (2013); this involves advancing the macro and micro models in time separately, with information exchange organised dynamically, depending on the level of time-scale separation.

This paper is structured as follows. In § 2 we introduce the multiscale flow problem in high-aspect-ratio geometries. Section 3 gives a detailed description of the new unsteady-IMM hybrid scheme. In § 4 we demonstrate hybrid simulations of a converging–diverging channel conveying a Lennard-Jones liquid. We conclude in § 5.

2. Multiscale flow in confined high-aspect-ratio geometries

2.1. Length-scale separation

We present a simulation technique for problems that are scale-separated in the flow direction (i.e. between the scales over which continuum variables vary and the microscopic molecular scales), but not in the direction perpendicular to the flow. For this type of high-aspect-ratio flow geometry we can evaluate the degree of streamwise length-scale separation by a dimensionless number:

(2.1) $$\begin{eqnarray}S_{L}(s)=\left|\frac{\text{d}h(s)}{\text{d}s}\right|^{-1},\end{eqnarray}$$

applied locally in the streamwise direction $s$ , where $h(s)$ is the variation of the channel height along $s$ . Equation (2.1) can also be extended to include streamwise variations in hydrodynamic variables and to channel curvature (see Borg et al. 2013c ; Patronis et al. 2013; Patronis & Lockerby 2014).

Provided $S_{L}\gg 1$ it can safely be assumed that, in small streamwise sections of the channel, the walls are approximately parallel to the centre streamline. This is a local parallel-flow assumption, and it enables the representation of small sections of the channel geometry by micro subdomains with exactly parallel walls and periodicity applied in the streamwise direction (Borg et al. 2013c ). Consequently, this is a convenient set-up for MD ensembles. See figure 1 for a schematic of the IMM treatment of an example channel.

The number of micro elements, which we denote as ${\it\Pi}$ , required to solve a given macro flow problem accurately is highly dependent on $S_{L}$ : more MD elements are required when $S_{L}\rightarrow 1$ , and fewer are required when $S_{L}\rightarrow \infty$ . In circumstances where there is no length-scale separation in some parts of the channel flow configuration, i.e. where $S_{L}\lesssim$ 1 at an arbitrary $s$ , the parallel-flow assumption no longer holds and a suitable alternative is to model the local geometry by an exact MD representation that is then coupled to the other IMM micro elements. This has been demonstrated in Borg, Lockerby & Reese (2013b ) and Stephenson et al. (2014) for junction components in network channel flows.

The choice of ${\it\Pi}$ affects the accuracy (which can be determined through successive refinement) and computational cost. The micro gearing dimensionless parameter is usually a good indicator of how much computationally cheaper a hybrid simulation is than a full-MD simulation:

(2.2) $$\begin{eqnarray}g_{L}=\frac{{\rm\Delta}s}{{\it\delta}s},\end{eqnarray}$$

where ${\rm\Delta}s\approx L/{\it\Pi}$ is the separation between micro elements (i.e. a macro cell size), $L$ is the length of the channel and ${\it\delta}s$ is the streamwise size of a micro element (see figure 1). Computational savings over a full-MD simulation are obtained when $g_{L}>1$ . In MD, ${\it\delta}s$ has to be larger than around $3\times$ the intermolecular potential cut-off radius in order to prevent finite size effects from affecting the results; so ${\it\delta}s=4.08~\text{nm}$ is chosen for all of the cases in this paper. The separation ${\rm\Delta}s$ of micro elements is assumed uniform along $s$ for numerical simplicity.

2.2. Time-scale separation

Up until now the IMM could not be applied to unsteady flow problems. Unsteadiness presents a truly multiscale problem, as there are potentially varying degrees of scale separation in both time and space that can be exploited. The degree of time-scale separation between the macro (continuum) and micro (molecular) models can be evaluated at a time $t$ as

(2.3) $$\begin{eqnarray}S_{T}(t)=\frac{T_{\mathit{macro}}(t)}{T_{\mathit{micro}}},\end{eqnarray}$$

where $T_{\mathit{macro}}(t)$ is the time-scale characterising the temporal variation in the macro variables (e.g. the period of an oscillatory flow, or the time for a macro variable to relax to a steady state) and $T_{\mathit{micro}}$ is the characteristic time taken for the variables in the micro model to relax to any step change in the macro variables (e.g. the start-up time from rest).

From a numerical perspective, time is discretised using time steps ${\rm\Delta}t$ and ${\it\delta}t$ for the macro and micro solvers, respectively. The maximum time step that is allowed for the MD micro solver is based entirely on stability, and is taken to be ${\it\delta}t=5~\text{fs}$ for all cases in this paper. The macroscopic time step is expressed as

(2.4) $$\begin{eqnarray}{\rm\Delta}t=\frac{T_{\mathit{macro}}(t)}{n_{\mathit{macro}}},\end{eqnarray}$$

where $n_{\mathit{macro}}$ is a large number that: (a) captures the macroscopic temporal variation of the hydrodynamic properties, (b) satisfies the Courant–Friedrichs–Lewy (CFL) condition given by $v_{s}{\rm\Delta}t/{\rm\Delta}s<1$ in one dimension, where $v_{s}$ is the streamwise velocity, and (c) assures the hybrid simulation is stable. What we notice from (2.4) is that ${\rm\Delta}t$ may change depending on $T_{\mathit{macro}}(t)$ , and hence the scale separation $S_{T}(t)$ . To enable this adaptability we need a generalised time-stepping procedure for the macro and micro time scales throughout the hybrid simulation.

The Continuous micro solution–Asynchronous, Intermittent coupling (CAI) time-stepping method of Lockerby et al. (2013) exchanges information between micro and macro models at well-defined ‘intervals’; the circles in figure 2 show the points at which the exchange occurs. A ‘micro interval’ consists of $N$ MD time steps (i.e.  $N{\it\delta}t$ ), while a ‘macro interval’ consists of one macro time step ${\rm\Delta}t$ . In situations with large time-scale separation, i.e.  $S_{T}(t)\gg 1$ , the macro interval can be set much larger than the micro interval (see the right-hand side of figure 2).

Figure 2. The CAI time-stepping scheme (Lockerby et al. 2013). Horizontal arrows indicate the incremental time-stepping of the macro (top) and micro (bottom) solver, while arrows indicate the direction of information exchange.

The physical approximation and gain in efficiency stem from asynchronously coupling the macro and micro models. The macro model advances further in time, each time step, than the micro model, but they are coupled as if they are of equal times. This enables the macro model to advance more quickly than it would otherwise, but the overall hybrid method is still acceptably accurate whenever the variations in macro variables are slow compared with the relaxation time scale of the micro model (i.e. if the macro and micro models are time-scale separated).

Of course, whenever time-scale separation does not exist this temporal gearing cannot be applied and the intervals must be made equal, as shown in the left-hand side of figure 2. In cases where the degree of scale separation in the flow varies (such as a sudden start-up flow evolving to steady state), the macro time step ${\rm\Delta}t$ as well as $N$ can vary with the time-local $S_{T}(t)$ , as explained below.

An indication of the computational speed-up provided by the CAI time-stepping method over a fully coupled time-stepping procedure (in which ${\rm\Delta}t=N{\it\delta}t$ ) can therefore be obtained using a micro gearing parameter for time:

(2.5) $$\begin{eqnarray}g_{T}=\frac{{\rm\Delta}t}{N{\it\delta}t},\end{eqnarray}$$

which is the ratio of macro to micro coupling intervals. Only when $g_{T}>1$ (i.e.  $S_{T}>1$ ) can there be computational savings over a fully coupled scheme.

As stated above, for flows of interest the scale separation may not necessarily be constant in time, hence the generality of $S_{T}(t)$ in (2.3). The CAI method enables adaptable gearing throughout the simulation, with a relationship between $g_{T}$ and $S_{T}$ suggested in Lockerby et al. (2013) as

(2.6) $$\begin{eqnarray}g_{T}=K_{g}(S_{T}-1)+1,\end{eqnarray}$$

where $K_{g}$ is a proportional gain that acts as a control on how aggressively scale separation is exploited at the expense of accuracy. Assuming ${\rm\Delta}t$ is known and dictated by the macroscopic variations, (2.6) can be substituted into (2.5) to obtain an expression for an adaptive choice of $N$ :

(2.7) $$\begin{eqnarray}N=\frac{{\rm\Delta}t}{{\it\delta}t[K_{g}(S_{T}-1)+1]}.\end{eqnarray}$$

3. The unsteady-IMM

Our hybrid method discretises the domain with ${\it\Pi}$ uniformly distributed micro elements, and advances the macro and the micro models separately on their own time line. Coupling information is exchanged between the macro and micro models at adaptive time intervals ${\rm\Delta}t$ and $N{\it\delta}t$ , respectively.

3.1. Macro model

The macro model consists of the unsteady, isothermal, compressible mass and momentum conservation equations:

(3.1) $$\begin{eqnarray}\frac{\partial {\it\rho}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\it\rho}\boldsymbol{V}=0,\end{eqnarray}$$

and

(3.2) $$\begin{eqnarray}{\it\rho}\frac{\partial \boldsymbol{V}}{\partial t}+{\it\rho}\boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{V}=-\boldsymbol{{\rm\nabla}}p+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}+\boldsymbol{F}_{s},\end{eqnarray}$$

where ${\it\rho}$ is the mass density, $\boldsymbol{V}$ is the convective velocity, $p$ is the hydrostatic pressure, $\unicode[STIX]{x1D64F}$ the shear stress tensor and $\boldsymbol{F}_{s}$ is a source of flow in terms of force per unit volume. In these equations $t$ refers to the macroscopic time scale. For low Reynolds number flows of high aspect ratio, (3.1) and (3.2) become

(3.3) $$\begin{eqnarray}\frac{\partial {\it\rho}}{\partial t}+\frac{\partial {\it\rho}v_{s}}{\partial s}=0,\end{eqnarray}$$

and

(3.4) $$\begin{eqnarray}{\it\rho}\frac{\partial v_{s}}{\partial t}=-\frac{\partial p}{\partial s}+\frac{\partial \unicode[STIX]{x1D64F}_{sy}}{\partial y}+\frac{\partial \unicode[STIX]{x1D64F}_{sz}}{\partial z}+F_{s},\end{eqnarray}$$

respectively, where $v_{s}$ and $F_{s}$ are the velocity and forcing in the streamwise directions. It is worth noting here that (3.4) contains no constitutive approximation (e.g. Navier–Stokes), transport model (viscosity) or slip boundary conditions; it is maintained in its general form.

Integrating the continuity equation (3.3) over the channel cross-sectional area $A$ we obtain

(3.5) $$\begin{eqnarray}\frac{\partial {\it\rho}}{\partial t}+\frac{1}{A}\frac{\partial {\dot{m}}}{\partial s}=0,\end{eqnarray}$$

where ${\dot{m}}$ is the mass flow rate and ${\it\rho}$ is the mass density; both variables are a function of time $t$ and streamwise location $s$ .

Equation (3.5) is the main macro model, applied in the streamwise direction, with ${\it\rho}(s,t)$ being the unknown macro state variable. The $\text{micro}\rightarrow \text{macro}$ coupling in this equation is imposed by extracting the mass flow rate ${\dot{m}}(s)$ from the micro state as explained in § 3.2.

For the periodic problems of this paper, the spatial derivative in (3.5) is evaluated from values at subdomain locations using a trigonometric interpolation through all subdomains, as described in Patronis et al. (2013).

3.2. Micro model

The micro model consists of standard MD (Allen & Tildesley 1987) that solves for the molecular movements and interactions in each subdomain independently (see the example of a subdomain in figure 1):

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}_{i}=\frac{\text{d}\boldsymbol{r}_{i}}{\text{d}t}, & \displaystyle\end{eqnarray}$$

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle m_{i}\frac{\text{d}\boldsymbol{v}_{i}}{\text{d}t}=\mathop{\sum }_{j}-\boldsymbol{{\rm\nabla}}\mathscr{V}+\boldsymbol{f}^{\mathit{ext}}, & \displaystyle\end{eqnarray}$$

where $m_{i}$ , $\boldsymbol{r}_{i}$ and $\boldsymbol{v}_{i}$ are the mass, position and velocity of the $i$ th molecule. On the right-hand side of (3.7), the first term is the net intermolecular force due to neighbouring molecules and the second term represents an external forcing applied to all liquid molecules. Without full macroscopic coupling, the latter is just the source term:

(3.8) $$\begin{eqnarray}\boldsymbol{f}^{\mathit{ext}}=\frac{F_{s}}{{\it\rho}_{n}}\hat{\boldsymbol{n}}_{s},\end{eqnarray}$$

where ${\it\rho}_{n}$ is the number density and $\hat{\boldsymbol{n}}_{s}$ is a unit vector parallel to the walls. Details of the interaction potential and parameters are given in appendix A, since these are independent of the hybrid method proposed here.

Each subdomain is periodic in the streamwise direction, so the mass in each subdomain is conserved. With an applied external forcing, the flow rate in one MD subdomain is therefore time-accurate and can be measured:

(3.9) $$\begin{eqnarray}{\dot{m}}(t)=\frac{1}{{\it\delta}sN}\mathop{\sum }_{{\it\tau}}^{N}\mathop{\sum }_{i}m_{i}\boldsymbol{v}_{i}\boldsymbol{\cdot }\hat{\boldsymbol{n}}_{s},\end{eqnarray}$$

where the first summation is over the time period $t\rightarrow t+N{\it\delta}t$ and the second summation is over fluid molecules in the MD subdomain.

The micro subdomain with forcing given by (3.8) does not, however, contain the complete $\text{macro}\rightarrow \text{micro}$ momentum coupling. The problem with a MD set-up that contains parallel walls and periodicity applied in the streamwise direction is that it cannot support a streamwise pressure gradient. Therefore the micro momentum balance of one MD subdomain, i.e. 

(3.10) $$\begin{eqnarray}{\it\rho}\frac{\partial v_{s}}{\partial t}=\frac{\partial \unicode[STIX]{x1D64F}_{sy}}{\partial y}+\frac{\partial \unicode[STIX]{x1D64F}_{sz}}{\partial z}+F_{s},\end{eqnarray}$$

is inconsistent with the macro balance of (3.4).

For coupled simulations to be accurate, the conservation of momentum must hold at both macro and micro scales. To produce in the micro model the hydrodynamically equivalent effect of the missing local pressure gradient, we apply an artificial body force ${\it\Phi}$ , which we call the pressure gradient correction, as an additional external forcing to (3.8) in the corresponding subdomain, i.e. 

(3.11) $$\begin{eqnarray}\boldsymbol{f}^{\mathit{ext}}=\frac{F_{s}}{{\it\rho}_{n}}\hat{\boldsymbol{n}}_{s}+\frac{{\it\Phi}}{{\it\rho}_{n}}\hat{\boldsymbol{n}}_{s},\end{eqnarray}$$

where

(3.12) $$\begin{eqnarray}{\it\Phi}=-\frac{\text{d}p}{\text{d}s}.\end{eqnarray}$$

Equations (3.11) and (3.12) therefore now ensure an appropriate momentum conservation; so the local macroscopic pressure gradient comprises the $\text{macro}\rightarrow \text{micro}$ coupling.

3.3. Closure requirements

The streamwise density variation ${\it\rho}(s)$ , which is obtained by solving the continuity equation (3.5), is imposed on the MD subdomains by a particle insertion/deletion protocol (in this case, the USHER algorithm of Delgado-Buscalioni & Coveney (2003)). The fluid density in each micro element is therefore updated on-the-fly during the simulation. However, as there is a change in streamwise density, the streamwise pressure $p(s)$ must be updated first before solving the coupling equation (3.12). This closure is resolved by defining an equation of state relating density to pressure, i.e. of the form $p({\it\rho})$ . While this equation could be obtained from reliable literature sources, any slight mismatch between the macro and micro fluid descriptions, for example when high confinement exists, can produce inaccuracies in the hybrid predictions when compared with the full-MD solution. Instead, we measure the macro pressure–density relationship from each MD subdomain in a presimulation, as outlined in Borg, Lockerby & Reese (2014), and fit a least-squares polynomial of the form

(3.13) $$\begin{eqnarray}p({\it\rho})=\mathop{\sum }_{j=0}^{M}a_{j}{\it\rho}^{(M-j)},\end{eqnarray}$$

where $a_{j}$ are the fit coefficients (see appendix B), and $M=4$ is the degree of the polynomial.

3.4. Algorithm overview

The hybrid solver has been written in OpenFOAM, a C++ toolbox that contains classes for creating computational fluid dynamics solvers and, recently, also MD solvers (Macpherson & Reese 2008; Ritos et al. 2013). The hybrid algorithm runs each MD subdomain in serial on an assigned central processing unit (CPU) core. As such, the micro subdomains run in parallel, and parallel communication of data (i.e. the ${\it\Pi}$ measured mass flow rates) is only required at each coupling interval, making this a very efficient parallelisation scheme. This approach also guarantees that the computational cost for a hybrid simulation is essentially independent of the number of subdomains ( ${\it\Pi}$ ), provided there are ${\it\Pi}$ CPU cores available. In the simulations in this paper, for example, we choose ${\it\Pi}=5$ , which is smaller even than the number of core processors found on many desktop computers.

Figure 3. Converging–diverging channel case: (a) the full-MD simulation, (b) the length-scale separation $S_{L}(s)$ taken from (2.1) and (c) schematic of the unsteady-IMM approach.

The procedure of the algorithm is briefly outlined here, but is also detailed in appendix C in proper algorithmic format. The forcings imposed on micro solvers are given by (3.11); initially ${\it\Phi}=0$ for each subdomain. After $N$ MD time steps, the mean mass flow rate is measured in each subdomain using (3.9) and a low-order polynomial is least-squares fitted to the data to obtain a spatially smoothed ${\dot{m}}(s)$ before trigonometric interpolation; this prefiltering step removes potentially destabilising MD noise. Time-stepping is applied as described in (2.3)–(2.7), to obtain new macro and micro time intervals. Next, the macro equations (3.5), (3.13) and (3.12) are solved numerically in this order to give the new values of macro densities ${\it\rho}(s)$ , pressures $p(s)$ and pressure gradient corrections ${\it\Phi}(s)$ , respectively, at the next macro time interval. The new pressure-gradient correction ${\it\Phi}(s)$ provides the $\text{macro}\rightarrow \text{micro}$ condition, so the forcings are updated in (3.11) and the hybrid algorithm proceeds back to the MD subdomains to obtain new mass flow rates for the next time interval. The number of molecules in each micro subdomain is modified using USHER at the beginning of each MD run to match the new local macroscopic density ${\it\rho}(s)$ .

4. Results

We apply our unsteady-IMM hybrid method to the flow of liquid argon along a periodic converging–diverging channel that has a smoothly varying height in the streamwise direction. A time-dependent gravity-type force is applied to the fluid to generate an unsteady/transient flow, although pressure-driven or boundary-driven flows can also be treated (Patronis et al. 2013). The solution dependence on ${\it\Pi}$ has been reported by Patronis et al. (2013) and Borg et al. (2013c ), and is not repeated here: we choose a value of ${\it\Pi}$ that is comfortably sufficient for accurate solutions in the geometry we consider. Our hybrid results are benchmarked against a full-MD simulation to assess computational expense.

Figure 4. Gravity force varying with time for various time-scale-separated channel flow cases: (a) case A, startup flow, (b) case B, oscillating force with small scale separation $S_{T}=2$ , (c) case C, oscillating force with large scale separation $S_{T}=100$ and (d) case D, oscillating force with mixed scale separation $S_{T}=2\rightarrow 100$ .

4.1. Validation

The full-MD converging–diverging channel shown in figure 3(a) has a length $L=68~\text{nm}$ in the streamwise direction $s$ , a depth of ${\rm\Delta}z=5.44~\text{nm}$ , and heights of 3.4 and 2.04 nm at the inlet and throat sections, respectively. The channel is periodic in both the streamwise and spanwise directions. The height between the top and bottom walls $h(s)$ varies in the streamwise direction according to a sinusoidal function, in order to make the variation in the scale separation gradual. For this channel geometry:

(4.1) $$\begin{eqnarray}h(s)=2a\left[\cos \left(\frac{2{\rm\pi}s}{L}\right)-1\right]+h_{\mathit{inlet}},\end{eqnarray}$$

where $4a=1.36~\text{nm}$ is the difference in height from inlet to throat and $h_{\mathit{inlet}}$ is the inlet height. The variation of length-scale separation in the streamwise direction $S_{L}(s)$ is obtained from (2.1) and shown in figure 3(b). The minimum $S_{L}$ for this configuration is ${\sim}$ 16.

All flows in both the full-MD and the hybrid simulations start from rest, and have the same spatial density and velocity distributions ${\it\rho}(s)$ , $v_{s}(s)$ . Then a time-varying gravity force $F_{g}(t)$ is applied in the $s$ direction. Four types of time-scale-separated flows are simulated (see figure 4).

1. (a) No time-scale separation. A start-up flow problem where an instantaneous gravity force of $F_{g}=0.487~\text{pN}$ is applied. This case is required to obtain an estimate of $T_{\mathit{micro}}\approx 108~\text{ps}$ .

2. (b) Small time-scale separation, $S_{T}=2$ . An oscillating gravity force is applied with amplitude $F_{g}=0.487~\text{pN}$ and period $T_{\mathit{macro}}=0.22~\text{ns}$ .

3. (c) Large time-scale separation, $S_{T}=100$ . An oscillating gravity force of the same amplitude is applied, but with a larger period $T_{\mathit{macro}}=10.8~\text{ns}$ .

4. (d) Mixed time-scale separation, $S_{T}=2\rightarrow 100$ . An oscillating gravity force of the same amplitude is applied, but with gradually increasing period: $0.22\rightarrow 10.8~\text{ns}$ .

The hybrid solver set-up of the converging–diverging channel is shown in figure 3(c), and consists of ${\it\Pi}=5$ micro elements distributed uniformly at streamwise locations:

(4.2) $$\begin{eqnarray}s_{i}=L(i-1)/{\it\Pi},\end{eqnarray}$$

where $i=1,2,\dots ,{\it\Pi}$ are the indices of the micro elements. The channel height of each micro element is a mapping from the macro domain geometry, which can be obtained by substituting the subdomain locations $s_{i}$ (4.2), in the channel geometry (4.1). In the streamwise direction all micro elements have a width of ${\it\delta}s=4.08~\text{nm}$ , which gives a micro subdomain separation (macro spacing) of ${\rm\Delta}s=13.6~\text{nm}$ , and a consequent computational speed-up estimate due to spatial gearing of $g_{L}=3.33$ . Table 1 gives an overview of the three time-scale-separated cases we consider, each with three different temporal gearings $g_{T}$ we test.

In order to improve the MD statistics, the unsteady-IMM uses 10 independent realisations for each micro element. Mass flow rate measurements taken in each micro interval ( $N{\it\delta}t$ ) are therefore further averaged across 10 ensembles. In the full-MD simulations we use either block-averaging or the wavelet proper orthogonal decomposition (WPOD) method (Zimón et al. 2014) to reduce noise.

In figure 5 we present our results for the streamwise variation of mass flow rate ${\dot{m}}(s,t)$ for various time instances during the start-up problem (Case A), and in figure 6 we present the mass flow rate, density and average velocity in each subdomain location as they vary in time. We observe very good qualitative agreement between the full-MD and the hybrid results. The mass flow rate measurements in figure 6 pick up an acoustic response at the beginning of the start-up problem, which is caused by the impulse forcing. Reassuringly, these features are resolved in both the full-MD and the unsteady-IMM solutions.

Figure 5. Full-MD and hybrid predictions of mass flow rates ${\dot{m}}(s,t)$ along the streamwise direction $(s)$ for various instantaneous times $(t)$ in the start-up problem (Case A).

Figure 6. Transient results for the start-up flow (Case A) at each micro subdomain location: (a) mass flow rate, (b) mass density and (c) average velocity are recorded from each MD subdomain location and compared with the result extracted from a corresponding bin of width 2.9 nm in the full-MD simulation.

Table 1. The micro gearing $g_{T}$ and the scale-separation sensitivity parameter $K_{g}$ for the oscillatory-forcing cases considered.

In figure 7(a–c) we show mass flow rates varying with time for the three oscillating cases B-D, respectively; we plot here only the results for the inlet MD subdomain of the converging–diverging channel (i.e.  $i=0$ , $s=0$ ) for comparison.

Figure 7. Predicted mass flow rate ( ${\dot{m}}$ ) varying with time ( $t$ ) at the inlet of the converging–diverging channel ( $s=0$ ) for all oscillatory-forcing cases: (a) Case B, (b) Case C and (c) Case D.

The hybrid results for the high-frequency oscillating flow (Case B) match well with the full-MD solution when $g_{T}=1$ (i.e. spatial gearing in the hybrid, but no temporal gearing), as shown in figure 7(a). In this case, as well as in Case A, the time-scale separation is low ( $S_{T}\leqslant 2$ ), meaning that the microscopic time scale is of similar magnitude to the oscillating period and so computational savings due to hybridisation in time cannot be obtained without loss of accuracy. To demonstrate this we include the simulation results for increasing values of $g_{T}$ in figure 7(a), that show that for less-time-scale-separated cases the time accuracy of the solution is highly sensitive to $g_{T}$ greater than one. Figure 8 shows the actual percentage error between the peaks of the full-MD and those of the unsteady-IMM simulation. For example, a change of gearing from $g_{T}=1$ to $g_{T}=2$ produces a ${\sim}$ 30 % error in the peaks of the mass flow rate, which grows to 60 % error when $g_{T}=4$ .

Figure 8. Percentage error in the peak mass flow rate values for the oscillating flow Cases B and C when our hybrid results are compared with the full-MD solutions. Sensitivity of this error to the gearing is clearly shown by the dotted lines (curve fits to the data) for the two time-scale-separated cases.

In figure 9 we present subdomain measurements of mass flow rate, density and velocity as a function of time in the unsteady-IMM ( $g_{T}=1$ ) solution of Case B. The hybrid method agrees very well with the full-MD simulation across all subdomain locations, and captures the temporal oscillations of mass density.

Figure 9. Transient results for the high-frequency oscillating unsteady flow (Case B) at each micro subdomain location: (a) mass flow rate, (b) mass density and (c) average velocity are recorded from each MD subdomain location and compared with the result extracted from a corresponding bin of width 2.9 nm in the full-MD simulation.

For the low-frequency oscillating flow (Case C), the time-scale separation is much larger ( $S_{T}=100$ ), and so the sensitivity to $g_{T}$ is much less, as shown qualitatively in figure 7(c) and quantitatively in figure 8. In the latter figure this sensitivity is contrasted with the less-scale-separated case ( $S_{T}=2$ ). Clearly, a reasonable accuracy ( ${<}$ 10 % error) would be obtained in Case C for $g_{T}<20$ . To make the analysis complete we also present the full range of subdomain results for Case C in figure 10. We discuss the observed secondary oscillations in § 4.2 below.

Figure 10. Transient results for the low-frequency oscillating unsteady flow (Case C) at each micro subdomain location: (a) mass flow rate, (b) mass density and (c) average velocity are recorded from each MD subdomain location and compared with the result extracted from a corresponding bin of width 2.9 nm in the full-MD simulation. The density profiles have been smoothed using block-averaging to minimise the short time-scale artefacts in the plots.

The final case, which consists of a mixed frequency oscillating flow (Case D), demonstrates the usefulness of the hybrid scheme: now the time-scale separation is a function of time. As such, the gearing parameter $g_{T}(t)$ adapts on-the-fly during the solution procedure, depending on the local time-scale separation $S_{T}(t)$ . The proportional gain parameter $K_{g}$ for these cases is indicated in table 1, and remains constant in each simulation. The resulting time variation of the gearing $g_{T}(t)$ and macro time-step ${\rm\Delta}t(t)$ are shown in figure 11(a,b), respectively, for the three chosen values of $K_{g}$ . Figure 7(c) clearly shows that these adaptive simulations give excellent results across the range of time-scale separation, while figure 12 shows the full range of results for $K_{g}=0.1$ . The high frequencies of the flow are resolved very accurately by the unsteady-IMM as it chooses a small gearing, and the lower frequencies of the flow because a higher gearing is chosen.

Figure 11. (a) Time micro gearing $(g_{T})$ and (b) macro time step ( ${\rm\Delta}t$ ) varying with time ( $t$ ), for the variable frequency oscillating flow problem in Case D. The gearing is an indication of the computational savings over a fully coupled case. Note that the discontinuity in plot (a) occurs because $N$ in (2.7) is an integer, and the nearest integer is calculated.

Figure 12. Transient results for the mixed frequency oscillating unsteady flow (Case D) at each micro subdomain location: (a) mass flow rate, (b) mass density and (c) average velocity are recorded from each MD subdomain location and compared with the result extracted from a corresponding bin of width 2.9 nm in the full-MD simulation. The density profiles have been smoothed using block-averaging to minimise the short time-scale artefacts in the plots.

4.2. Short time-scale artefacts

In general, our multiscale simulation approach balances micro-scale resolution with macro-scale prediction. We capture relevant small-scale phenomena in order to predict large-scale phenomena. However, the results in figure 7(b) highlight an important cautionary example of how temporal gearing can create erroneous artefacts on short time scales (i.e. short-period oscillations), which might otherwise give the impression of being physical.

The source of this high-frequency oscillation is, in fact, acoustic. To estimate the natural acoustic frequency of the flow geometry, consider a micro and a macro model for an inviscid flow in a straight channel of length $L$ :

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\rho}}{\partial t}+\frac{\partial {\it\rho}u}{\partial s}=0,\quad (\text{macro}) & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\rho}u}{\partial t}+\frac{\partial p}{\partial {\it\rho}}\,\frac{\partial {\it\rho}}{\partial s}=0,\quad (\text{micro}). & \displaystyle\end{eqnarray}$$

Coupling the two models produces a wave equation:

(4.5) $$\begin{eqnarray}\frac{\partial ^{2}{\it\rho}}{\partial t^{2}}=c^{2}\frac{\partial ^{2}{\it\rho}}{\partial s^{2}},\end{eqnarray}$$

where the wave speed $c=\sqrt{\partial p/\partial {\it\rho}}$ , i.e. the speed of sound. The fundamental acoustic frequency is then determined by the length of the channel: $f=c/L$ .

For the situation where temporal gearing is applied, the micro and macro equations are modified:

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle g_{T}\frac{\partial {\it\rho}}{\partial t}+\frac{\partial {\it\rho}u}{\partial s}=0,\quad (\text{macro}) & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\rho}u}{\partial t}+\frac{\partial p}{\partial {\it\rho}}\,\frac{\partial {\it\rho}}{\partial s}=0,\quad (\text{micro}) & \displaystyle\end{eqnarray}$$

which combine to produce a fundamental acoustic frequency with a dependence on the temporal gearing:

(4.8) $$\begin{eqnarray}f=\frac{c}{L\sqrt{g_{T}}}.\end{eqnarray}$$

In figure 13 this expression (derived for a straight channel) is compared with the frequencies of the high-frequency oscillations observed in our converging–diverging channel (extracted from the results of Case A and Case B using a Fourier transform); there is a clear correlation.

Figure 13. Frequency of secondary oscillations varying with micro gearing. The secondary frequencies are taken from fast Fourier transforms of the mass flow rate measurements in the first micro subdomain, for Cases A and C. Comparisons are made with the frequency prediction from (4.8). Vertical lines indicate the range of density fluctuations observed in the hybrid simulations, which affect the speed of sound $c$ in (4.8).

However, this does not explain why such acoustic frequencies are not observed in the full-MD simulation (which is without temporal gearing). The reason is that temporal gearing reduces the ratio of the acoustic period to the viscous development time. This means that ‘echoes’ in the flowfield, generated by flow disturbances and reinforced by the periodicity of the domain, are less prone to being damped out by viscosity. In other words, the Wormersley number (the ratio of the oscillatory inertia force to the shear forces) is increased when we employ temporal multiscaling.

The existence of such oscillations, which are not a numerical artefact but are an artefact of the temporal gearing, reinforce the need for multiscale simulations to go hand-in-hand with gearing sensitivity studies (see § 4.4) in order to identify whether a short time-scale phenomenon is physical or not. But there are other practical issues that also arise. The high-frequency oscillations mean the particle insertion protocol (here, the USHER algorithm) has to work very hard to track high-frequency changes in density. Insertion of molecules in MD has a prohibitive computational overhead when the fluid is dense, as discussed in Delgado-Buscalioni & Coveney (2003) and Borg et al. (2014).

4.3. Computational savings

The primary motivation for adopting any molecular-continuum scheme is to obtain accurate results at a lower processing cost than a full-MD simulation. First, we present in this section the measured computational speed-up of our hybrid simulations. This is obtained by dividing the total time spent for the full-MD simulation ( $t_{\mathit{comp}}^{\mathit{F}}$ ) by the time taken for the corresponding hybrid simulation ( $t_{\mathit{comp}}^{\mathit{H}}$ ), i.e. ( $t_{\mathit{comp}}^{\mathit{F}}/t_{\mathit{comp}}^{\mathit{H}}$ ). The speed-ups for all of the cases above are shown in figure 14. In order to make a fair comparison, the computational times for the full-MD and the hybrid simulations have been normalised by the number of allocated processors, as well as by the processing time required to reach the same signal-to-noise ratio.

Figure 14. Computational speed-up varying with micro gearing (in time) for all four Cases A–D relative to the full-MD simulation. Symbols denote exact speed-up measurements taken from processor clock times, while dashed lines are an estimate of the speed-up given by $g_{L}\times g_{T}$ . Note, Case D has a varying gearing, so for this plot we select the largest gearing.

Naturally, the computational speed-up is modest for cases where accuracy is more sensitive to the micro gearing $g_{T}$ , such as in the less-time-scale-separated Cases A and B. However, a sizeable speed-up is achieved for the more scale-separated Cases C and D. Note that the computational savings in figure 14 should be compared with the level of accuracy in figure 8 in order to make a practical compromise. (However, in practice, when full-MD is too costly for a benchmark case, the accuracy of a hybrid result can be assessed by running spatial and temporal gearing independence studies.) The ability to apply gearing is an important feature of the hybrid method proposed here, and it is one that cannot be exploited in a full-MD approach.

A useful task before running any simulations is to obtain an estimate of how much speed-up a hybrid simulation will potentially give. This will help decide which simulation approach to use: a full-MD or a hybrid approach. The spatial and temporal gearings, $g_{L}$ and $g_{T}$ respectively, can be used as an a priori indication of speed-up, with the total overall speed-up given by their product: $g_{L}\times g_{T}$ . To demonstrate how well this fits with the actual speed-ups we found, we plot this measure as dashed lines in figure 14. The agreement between the estimation and the actual speed-up is reasonable for most cases. That there is some computational overhead in the unsteady-IMM simulations of the more scale-separated cases is probably due to the cost incurred by the USHER algorithm in resolving the acoustic oscillations in density, as explained in § 4.2.

4.4. Highly scale-separated flows

In the previous tests, the geometry of the converging–diverging channel, as well as the time scales of the macroscopic variables, were chosen such that a full-MD simulation would be able to validate our hybrid results. We now demonstrate the full potential of the unsteady-IMM method by showing its applicability to cases where much larger length- and time-scale separations exist, and show that speed-ups of orders of magnitude can be obtained. In such cases it is not practical to run a full-MD simulation, but we perform a gearing-dependency study to provide confidence in the predictions (as is done similarly for mesh selection in computational fluid dynamics, for example).

The converging–diverging channel with the geometry described by (4.1) is made successively longer in the streamwise direction for three cases: $L=1~{\rm\mu}\text{m}$ , $L=10~{\rm\mu}\text{m}$ and $L=100~{\rm\mu}\text{m}$ . These are typical length scales for flows, e.g. through carbon nanotube membranes (e.g.  $L=34{-}126~{\rm\mu}\text{m}$ in Majumder et al. (2005)). The test problem we carry out in this section is a start-up flow, as in figure 4(a), where the gravity forcing in the three cases is set to $F_{G}=0.24~\text{pN}$ , $F_{G}=0.015~\text{pN}$ and $F_{G}=3.89~\text{fN}$ , respectively.

The start-up problem for relatively long channels exhibits a variety of time scales. At the beginning of the simulation the dynamics are fast, and the mass flow rate approaches the steady-state value very quickly, within ${\sim}$ 200 ps. In this time period the time-scale separation number is small, $S_{T}=1$ , and so in all cases we use a tight coupling between macro and micro time scales, i.e.  $g_{T}=1$ . The remainder of the simulation has much slower dynamics, in particular, the pressure and density have much longer development times that are dependent on the length of the channel. This long development time in micro/nano channels has been observed in, e.g., Patronis et al. (2013). During this time the scale separation is much larger, and the macroscopic time intervals can be made much greater than the micro time intervals without much loss in accuracy. This produces orders of magnitude of computational savings. For the first channel length case, $L=1~{\rm\mu}\text{m}$ , we find that the pressures in the channel reach a steady-state relatively quickly, and so we maintain a temporal gearing of $g_{T}=1$ throughout the simulation. For the two longer channels, the development times are significantly longer and we allow the gearing to vary adaptively, depending on the case, as listed in table 2.

Table 2. Computational speed-up estimates given by $g=g_{L}\times g_{T}$ for the two long channels.

The pressure results evolving in time in all three cases are shown in figure 15. The hybrid simulations for the shortest channel case ( $L=1~{\rm\mu}\text{m}$ ) have no temporal gearing ( $g_{T}=1$ ), so the profiles are presented here for comparison with the other two cases. The profiles in the medium-length channel ( $L=10~{\rm\mu}\text{m}$ ) contain features of the shortest channel results across all of the subdomains, but naturally these occur over much longer time scales. What is also noteworthy is that the gearing chosen here does not produce any noticeable effect on the accuracy. For the pressure profiles in the longest channel ( $L=100~{\rm\mu}\text{m}$ ), some deviations are observed across the three gearings chosen, but this may be because of the noisy (negligible) flow rates obtained in these simulations.

Figure 15. Pressure $p$ (MPa) varying with time $t$ (ns) in each micro subdomain for the three different lengths of the converging–diverging channel.

Table 2 also gives estimates of the total speed-up $g$ that our hybrid simulations would obtain relative to full-MD simulations: the two longer cases are estimated to have speed-ups of $O(10^{4})$ and $O(10^{6})$ , respectively. This is many orders of magnitude larger than speed-ups reported for other hybrid methods that use MD as their micro solver, i.e.  $O(10^{1})$ (Hadjiconstantinou & Patera 1997). Our approach therefore represents a positive step towards practical simulation of nano-channel flows in realistically large geometries, over longer time scales.

5. Conclusions

We have presented an unsteady internal multiscale flow method (‘unsteady-IMM’) for simulating high-aspect-ratio flows with nano-scale confinement. Both length- and time-scale separation are treated. What is needed to demonstrate that a hybridisation scheme is accurate is that the hybrid model provides the same results as a full-resolution simulation over a range of conditions; what is needed to demonstrate that it is useful is that it obtains these results at a fraction of the computational cost of the full-resolution simulation. We have demonstrated both of these in the current paper. We have validated the unsteady-IMM method against full-MD simulations of a converging–diverging channel for various time-scale-separated flows. In all cases we established the sensitivity of the accuracy to the temporal gearing, which we showed is dependent on the time-scale separation $S_{T}$ . We also demonstrated the potential of the method for highly scale-separated flows for which full-MD simulations are not possible.

There is much further work possible. Here, $S_{T}$ has been estimated at the outset of the problem. Future work could involve a step that determines $S_{T}$ on-the-fly. Noise inherent in the MD subdomain simulations is also another problem. This could be resolved by smoothing techniques such as proper orthogonal decomposition (POD) or wavelet POD (WPOD) methods (Zimón et al. 2014). The unsteady-IMM method is general, and so may be used with other micro models, such as the direct simulation Monte Carlo technique for modelling micro-channel rarefied gas flows.

Evidently our method is limited by geometrical factors: unsteady-IMM may become less efficient as the channel width becomes larger, i.e. when the flow is no longer confined and non-continuum effects no longer dominate. In such cases, however, scale separation may exist in the direction normal to the wall and so an alternative hybrid approach would be the HMM-FWC method of Borg et al. (2013a ). We also expect some limitations to the method as the width of the channel becomes much smaller, for example, for strands of molecules flowing in an aquaporin. In this case, however, a full-MD simulation may be tractable.