Fluid simulations with atomistic resolution: a hybrid multiscale method with field-wise coupling

We present a new hybrid method for simulating dense ﬂuid systems that exhibit multiscale behaviour, in particular, systems in which a Navier-Stokes model may not be valid in parts of the computational domain. We apply molecular dynamics as a local microscopic reﬁnement for correcting the Navier-Stokes constitutive approximation in the bulk of the domain, as well as providing a direct measurement of velocity slip at bounding surfaces. Our hybrid approach di ↵ ers from existing techniques, such as the heterogeneous multiscale method (HMM), in some fundamental respects. In our method, the individual molecular solvers, which provide information to the macro model, are not coupled with the continuum grid at nodes (i.e. point-wise coupling), instead coupling occurs over distributed heterogeneous ﬁelds (here referred to as ﬁeld-wise coupling). This a ↵ ords two major advantages. Whereas point-wise coupled HMM is limited to regions of ﬂow that are highly scale-separated in all spatial directions (i.e. where the state of non-equilibrium in the ﬂuid can be adequately described by a single strain tensor and temperature gradient vector), our ﬁeld-wise coupled HMM has no such limitations and so can be applied to ﬂows with arbitrarily-varying degrees of scale separation (e.g. ﬂow from a large reservoir into a nano channel). The second major advantage is that the position of molecular elements does not need to be collocated with nodes of the continuum grid, allowing the resolution of the microscopic correction to be adjusted independently of the resolution of the continuum model. This means the computational cost and accuracy of the molecular correction can be independently controlled and optimised. The macroscopic constraints on the individual molecular solvers are artiﬁcial body-force distributions, used in conjunction with standard periodicity. We test our hybrid method on the Poiseuille ﬂow problem for both Newtonian (Lennard-Jones) and non-Newtonian (FENE) ﬂuids. The multiscale results are validated with full-scale molecular dynamics simulations of the same case. Very close agreement is obtained for all cases, with as few as two micro elements required to accurately capture both the Newtonian and non-Newtonian ﬂowﬁelds. The multiscale method converges very quickly (within 3–4 iterations) and is an order of magnitude more computationally e  cient than the full-scale simulation.


Introduction
In an important class of fluid dynamics problems the overall physics is multiscale, to the extent that microscopic processes strongly dictate the macroscopic behaviour. For example, granular flows in avalanches, shockwaves in high-speed rarefied flows, rain filtration through soil, weather changes etc. The present paper is specifically related to micro-and nano-scale fluid dynamics, where atomistic processes occurring over pico-or nano-scales determine the bulk e↵ects occurring over microand milli-scales, both temporally and spatially. The development of future technologies depending on nano-and micro fluidics requires methods that resolve multiscale phenomena accurately and e ciently.
Molecular Dynamics (MD) is now a recognised computational tool for accurately modelling fluid behaviour at atomistic scales as an ensemble of discrete molecules that obey fundamental laws of physics. In essence, a flow problem can be entirely described by such a microscopic model. While this is suitable for simulating small systems over short time-scales, such as water flows through carbon-nanotubes, in nano and micro engineering applications, it becomes computationally intractable to resolve all the important space and time scales involved.
In such flow problems, the continuum conservation laws can be used as an an accurate, though incomplete, model for macroscopic flow behaviour. Molecular dynamics can then be employed to replace, or close equations, in parts of the computational domain where a classical macroscopic constitutive or boundary model does not exist (e.g. liquids next to surfaces, rheological fluids, chemical reactions, etc). This is the so-called hybrid methodology, in which "the best of both worlds" in terms of computational cost and accuracy can be attained by unifying disparate macroscopic and microscopic solvers in a coupled simulation.
The most common approach to hybrid methods for dense fluids is domain-decomposition (DD) [1][2][3][4][5][6]. As the name implies, the computational domain is divided uniquely into macro and micro parts, providing only an overlapping region for matching hydrodynamic flux or state properties that act as Dirichlet or Neumann-type boundary conditions at the hybrid interface. Figs. 1(a) and (b) illustrate the DD method applied to a simple Poiseuille flow problem. Despite providing a reasonable rationale for modelling near-wall flows, DD schemes have some disadvantages. First, DD can only be used for flows that have a bulk flow region for which a constitutive model is known and is accurate. Furthermore, while it may be reasonably straightforward to segment a domain explicitly into micro/macro subdomains if it is a simple 1D system, it may be more challenging if the hybrid interface has more complex features [7]. Another issue is that of applying acceptable 'open-system' or 'macro-micro' boundary conditions at the non-periodic MD interface, which is a topic that is not well developed in the literature. Finally, DD methods do not generally work well for systems that exhibit varying degrees of scale separation, and may be of comparable computational cost to a full MD simulation of the case [8,9].
A di↵erent hybrid methodology that exploits the existence of scale separation between molecular and continuum processes was proposed by E and co-workers [10,11]. This was dubbed the Heterogeneous Multiscale Method (HMM), and is illustrated in Fig. 1(c). The central idea of HMM is that continuum governing equations for mass, momentum and energy conservation are able to predict the overall macroscopic behaviour, with only the missing/unknown fluid-constitutive and boundary information (such as shear-stress and slip) being provided by small independent MD elements distributed at the nodes of the computational mesh. These micro elements can be either (a) bulk nodes or (b) nodes on the boundary wall. Prior to supplying this information, every MD element is constrained using continuum fluid properties obtained at the nodal point, such as the strain tensor or temperature gradient. This point-wise coupling [12] between the computational grid and the individual molecular elements is valid for scale-separated conditions, where continuum property variations can be assumed to be small over the scale of the micro element (which has a small but finite size). In such conditions, the HMM approach naturally decouples processes occurring at di↵erent scales and can thus achieve major computational savings over a full MD simulation. The HMM has shown its potential in dense fluid simulations by various authors [12][13][14]. While these hybrid simulations all follow the original methodology proposed in [13] (as described above), they di↵er from each other in the way bulk-type micro elements are constrained from the local strain-rate (which is a second rank tensor) supplied by the continuum solution. For example, Ren and E [13] adopt a concept similar to the Lees and Edwards shear boundary condition [15] to impose a two dimensional strain-rate. In principle each MD box is set to deform dynamically every MD time-step, based on the local macroscopic velocity field. A disadvantage of this deformation method is that the MD box may become too skewed, and the authors solve this by expensive reinitialisation steps of the MD box from the final configuration state. Yasuda and Yamamoto [14] instead apply a rotational transformation on the strain-rate tensor, to change it from the continuum co-ordinate system to a co-ordinate system more convenient for MD simulations. In particular, the transformation ensures that the diagonal terms of the new strain-rate tensor all become zero, which when applied to 2D isotropic problems gives only one non-zero strain-rate term in the o↵-diagonal. This means that the one-dimensional Lees and Edwards boundary method can be conveniently applied without problems. The shear-stress tensor that is measured from the micro elements is then supplied back to the macro model after being transformed to the original coordinate system. Asproulis et al. [12] strain the MD box using the Parrinello-Rahman deformation technique [16] in conjunction with the SLLOD algorithm [17]. This can be applied more generally to three dimensional shearing, because the local strain rate is applied directly to the equations of motion of all molecules.
However, the second type of micro elements in HMM simulations provide the local wall boundarycondition information, so the MD boxes are prohibited from adopting the same box-deformation methods as bulk micro elements. Instead, bounding boxes of wall micro elements are kept fixed, and velocity-type boundary conditions only are applied in a 'controlling region' next to the non-periodic boundary, in conjunction with periodic boundary conditions in the flow direction. In essence, there are many similarities of HMM wall micro elements with DD schemes, such as that proposed by Hadjiconstantinou and Patera [2], i.e. a 'Maxwell-Demon' is applied only in the control region and it resamples new molecular velocities from a Maxwell-Boltzmann distribution at the target continuum temperature and local velocity supplied from the overlaying continuum solution.
The biggest drawback of the HMM is that it is only suitable for completely scale-separated flow systems, due to the point-wise macro-micro coupling. For example, in simulations of water, a typical minimum MD element size is x micro =5 nm (to avoid wrap-around e↵ects). For the HMM method to be worthwhile and not be more expensive than full MD, the continuum grid spacing ( x macro ) must then be at least greater than 5 nm, which represents a very severe restriction on the minimum scale of hydrodynamic variable variation that can be resolved with HMM in nano/micro-flows. Such a restriction precludes the use of HMM for technological applications that exhibit a mixed degree of scale separation. For example, where regions of the domain have characteristic scales much larger than the minimum MD scale (and where large savings can be made) and other regions have characteristic scales comparable to or smaller than the MD subdomains (e.g. a carbon nano tube connecting two reservoirs for filtration applications [18]).
Another significant drawback in HMM of the micro elements and macro nodes being collocated (i.e. there being a micro element at every continuum node) is that the micro and macro resolution are intrinsically linked. If the number of micro Molecular Dynamics simulations are to be reduced, the macro resolution must be reduced, and vice versa. This may seem like a reasonable trade-o↵, given that generally in computational methods greater resolution comes with a computational penalty. However, in many cases this macro-micro link is unnecessarily restrictive and limiting. This is because, generally, micro MD elements provide the missing information to the macro model relating to stresses and heat-flux. But such quantities typically vary more slowly than the macro hydrodynamic variables of velocity, temperature and density, and thus will not necessarily require the same resolution. Consider the Navier-Stokes-Fourier constitutive relations: temperature and velocity variations are of a higher spatial-derivative order than the related heat flux and stress, respectively (e.g. in classical Poiseuille flow, the velocity profile is parabolic while the stress variation is linear). Although the Navier-Stokes-Fourier constitutive relations are not expected to hold fully at the nano scale, it is likely that the same tendency will exist. This means that, in general, lower resolution is required for the component of the simulation that is providing information on correcting stress and heat flux (the micro elements) than is required of the component of the simulation resolving hydrodynamic variables (the macro grid).
The approach we present in this paper is a major modification to standard HMM to circumvent the limitations described above (i.e. to improve e ciency and increase the range of applicability). This modification is in the form of a Field-Wise Coupling (FWC) of the macro grid and micro elements (in shorthand we refer to this method as HMM-FWC). The micro elements in the method now represent fields correlating directly with identically sized continuum-flow subregions, and they can be positioned independently of the continuum grid; see Fig. 1(d). This enables application to flows with mixed degrees of scale separation (where hydrodynamic variables can vary substantially over the dimensions of a micro element) and allows the position of the micro elements to be optimised separately to the macro grid.
The remainder of this paper is structured as follows. In Section 2 we describe the new multiscale method in detail, providing a step-by-step account of the algorithm. In Section 3 we test the method on a simple and canonical flow problem for Newtonian and non-Newtonian fluids. In Section 4 we draw conclusions and discuss opportunities for further work.

A Heterogeneous Multiscale Method with Field-Wise Coupling (HMM-FWC)
As in conventional HMM, we use the same fundamental notion that a continuum-fluid conservation model is applicable in the entire flow domain. However, instead of using MD to replace constitutive and boundary field data in the continuum model, we solve the macro description using a modified Navier-Stokes formulation; locally constrained micro simulations enter the hybrid formulation to provide constitutive and boundary-slip corrections only where needed.
For the current work we assume steady, incompressible and isothermal conditions. The governing three-dimensional macro-conservation equations for mass and momentum equations are: where U = (u, v, w) is the fluid velocity, ⇢ is the mass density, p is the pressure, F is a body force per unit volume, µ is an assumed constant viscosity coe cient, and is a tensor field containing a correction to the Navier-Stokes constitutive model: where T is the deviatoric stress tensor and e = 1 2 (rU + rU T ). The basic strategy of our multiscale method is to calculate an interpolated stress-correction field over the entire domain from evaluation of Eqn. (3) at the distributed MD micro elements. From this, U and p are then found from a solution of Eqns. (1) and (2).
The primary characteristic of our method distinguishing it from HMM as presented by Ren and E [13] is that coupling does not occur at every continuum node, and is applied over fields rather than points (compare Figs. 1(c) and (d)). As discussed in the previous section, the assumption underpinning point-wise coupling in HMM is not valid when the hydrodynamic variation occurs at the same scale as the micro elements. In such cases, a direct physical correlation between macro and micro space is needed, requiring a flowfield of arbitrary complexity (e.g. a strain-rate field) to be simulated within each micro element. The imposition of this field on a micro element is not possible using techniques developed in the HMM literature, and so we propose a new constraint procedure that is consistent for both bulk and wall micro elements (explained in §2.3 and §2.4 below).
For the dual purpose of providing a detailed worked example, and also to prove the concept for a canonical flow, application of the method to an isothermal 1D Poiseuille flow geometry is described below. However, our underlying method is general, and may be extended to higher dimensions and other flow applications. We address the issues associated with advancing the method to two and three dimensions in the conclusions section.

Overview of the algorithm
Our multiscale method is an iterative scheme that combines solutions from macro and micro solvers until a converged global solution is reached. The flow diagram of each sub-process in one multiscale iteration is illustrated in Fig. 2, and briefly described as follows: 1. Macro solution (see §2.2) A solution to the continuum Eqns. (1) and (2) is obtained for U and p over the entire domain, with a stress-correction field ( ) taken from the previous iteration. Initially, =0.
2. Macro-to-micro projections (see §2. 3) The macro solution is projected onto all micro subdomains. These projected flowfields serve as target conditions for the individual MD simulations to recreate. Each micro subdomain is divided into a 'core' region and a 'constrained' region. In the core region, the velocity field projection is a one-to-one mapping from the macro solution.
In the constrained regions, however, the macro solution is modified so as to be compatible with the periodic boundary conditions of the MD elements; as such, it is only data obtained in the core regions that is physically accurate. a predetermined period 2 .

5.
Micro-to-macro compression (see §2.6) Field measurements of U and T are obtained from each micro element, and a least-squares polynomial fit is performed to extract continuous data.
6. Stress-correction (see §2.7) The stress correction field is calculated from each micro element using Eqn. (3) and interpolated over the entire macro domain, to be used in Step 1 of the next iteration. Other stress-correction fields are calculated for the constrained regions of the micro subdomains, for use in Step 3 of the next iteration.
In the following sections we explain these six steps in detail, for the worked example of Poiseuille flow.

Macro solution (Step 1)
For the Poiseuille-flow configuration of Fig. 1(a), the continuity Eqn. (1) is automatically satisfied and the stress-corrected x-momentum Eqn. (2) simplifies to: where x is a coordinate in the streamwise direction, y is a coordinate across the channel, F x is a body force driving the flow (which is constant for Poiseuille flow), and (= xy ) is a shear component of the stress-correction tensor: where ⌧ xy is the true shear stress. Eqn. (4) with (5) has no constitutive approximation and therefore holds for both molecular and continuum modelling treatments of the Poiseuille problem. In fact, the coupled equations are applicable even when the viscosity is not known exactly; an estimate for µ will su ce, since the stress-correction term accommodates for any inaccuracies in the Navier-Stokes constitutive relation. A solution to Eqn. (4) can be obtained for the Poiseuille flow case by integrating twice in the y-direction, over half the channel height H, giving a semi-analytical expression for the velocity profile u(y): There are two unknowns in this expression: the stress correction field (y), and the slip velocity u B . Initially both these are zero, but in subsequent iterations these properties are evaluated from the MD micro elements. We explain how these properties are measured and computed in §2.6 below.
For more complex problems, where an analytical expression is not possible, the macro solution could be obtained using standard computational fluid dynamics (CFD), such as the finite volume method. Note, in the molecular simulations the body force per unit volume F x is converted to an actual force using: f x = F x /⇢ n , where ⇢ n is the number density.

Macro-micro projection (Step 2)
The aim of the field-wise coupling approach is for the individual micro elements (see Fig. 3) to reproduce the same flowfield that occurs over the corresponding region in the continuum domain. To this end the estimated macroscopic velocity field u(y), which is obtained from the stress-corrected momentum Eqn. (4), is projected onto the domain of each micro element. This projection is indicated by the red dashed line in Fig. 4(b). This projection is then modified in the 'constrained' regions, in order to be consistent with the periodic boundary conditions 3 of the MD solver (see solid blue line in Fig. 4b). In these constrained regions, where the velocity profile has been artificially modified, body forces are applied (calculated in the next step, §2.4) in order to reproduce the modified velocity profile. No artificial body-forcing is applied in the core regions.
The flow velocity projection (denoted by the superscript ⇤ ) in the core region of the ith micro element (u ⇤ i ) is obtained directly from the macro solution, i.e.   where y 0 i is a new local co-ordinate system used to identify the core region, as indicated in Fig. 4. The height of the core region is denoted by h core . In the constrained region the velocity projection is modified with the intention of bridging the velocities at the two sides of the core region with a continuous velocity field (as well as continuous first and second derivative velocity fields) that pass through the top and bottom periodic boundaries.
To preserve continuity over the two interfaces of velocity and its first and second derivative, a sixth order polynomial is used: where y 00 i is a di↵erent local co-ordinate system used to identify the constrained region (see Fig.  4), h cs is the height of the constrained region, and a k,i are the polynomial coe cients which are uniquely determined through the following six boundary conditions: For the micro elements touching the walls of the Poiseuille flow problem, there are issues of non-periodicity. All that is then required to be able to use periodic boundary conditions is for the wall micro element to be a symmetrical flow domain with the line of symmetry running through the constrained region or wall, as shown in Fig. 5(a). The velocity projection for the wall micro element, both for the core and constrained regions are identical to those in the bulk and can be similarly calculated using Eqns. (7) and (8), as illustrated in Fig. 5(b).

Micro constraints (Step 3)
The modified flowfield of the constrained regions can be generated using an artificial body force f ⇤ x , as shown in Figs. 4(c) and 5(c) for bulk and wall micro elements respectively. To calculate this force distribution we substitute the target velocity field (8) into the flux-corrected Navier-Stokes Eqn. (4): where cs,i is the stress correction field in the constrained region of the ith subdomain and is calculated in the previous iteration (Step 6, §2.7). Note, in the initial iteration cs =0. In the core regions of the micro subdomains, molecular dynamics is performed without any artificial forcing imposed. As such, the forcing in the core region takes that of the Poiseuille driving force only, i.e.

Micro solution (Step 4)
All our micro elements are simulated using the open source non-equilibrium Molecular Dynamics code, mdFoam [7,19,20], which is implemented within OpenFOAM's C++ libraries [21]. Detached micro elements located in the bulk of the flow domain (bulk elements, see Fig. 3) each consist of a fully-periodic domain that is filled solely with fluid molecules at the target density and temperature defined by the problem. The micro element next to the channel bounding surface (wall elements, see Fig. 3) also include wall-material molecules, which are modelled as rigid atoms with parameters chosen from [22] to generate a reasonable slip behaviour.
For the period of the simulation, all fluid molecules in each micro element evolve according to Newton's equations of motion: where k = (1, . . . , N) is a particular molecule in the system of N molecules, m k is the molecule mass, r k = (x k , y k , z k ) is its position and v k = (u k , v k , w k ) is its velocity. The total intermolecular force, f 0 k in Eqn. (13) is computed from pair forces with neighbouring molecules by: where U (r kj ) is the interaction potential between fluid-fluid or solid-fluid interactions. In this paper we simulate for simplicity monatomic liquid argon using the 12-6 Lennard-Jones (LJ) potential (although water or other molecules can also be simulated straightforwardly): where ✏ = 1.6568 ⇥ 10 21 J and = 0.34 nm are the potential's characteristic energy and length scales, and r kj = |r k r j | is the separation of two arbitrary molecules (j, k) within a cut-o↵ radius, r cut = 1.36 nm. The equations of motion are numerically discretised in space and time using a timestep t which is 5.4 fs for the LJ simulations in §3.1 and 2.2 fs for the FENE chain simulations in §3. 2.
The term f ext k in Eqn. (13) is a locally-applied external force derived from the constraint step (Step 3, §2.4). In this fluid problem, f ext (y) = f ⇤ x,i (y)n x is a force distribution applied in the constrained region according to the macro-micro constraint (10), and f ext = f xnx is the macro force applied in the core region according to Eqn. (11), wheren x is the unit vector in the streamwise xdirection. In order to maintain isothermal conditions, the heat introduced by the induced shearing is removed using a velocity unbiased Berendsen thermostat [23] in bins of 0.68 nm thickness in the y-direction.
The size of the constrained region is picked for convenience to be the same as the core region. This produces strain rates in both regions that are similar magnitude, and thus places similar requirements on the thermostat. Almost certainly this region could be shrunk for greater computational savings with no penalty in accuracy. However, if the region is chosen to be excessively small (e.g. ⇠1 molecular diameter width), an extremely large rate of strain is applied locally, which in turn will generate an amount of heat extremely di cult to remove homogeneously using the thermostat.

Micro compression (Step 5)
After all micro elements reach a relaxed steady-state, measurements are obtained using a cumulative averaging technique to reduce noise. Each micro element is divided into spatially-oriented bins in the y-direction in order to resolve the velocity and shear-stress profiles. Velocity in each bin is measured using the Cumulative Averaging Method (CAM) [24], while the stress tensor field is measured using the Irving-Kirkwood relationship [25]. A least-squares polynomial fit to the data is performed, which helps further to reduce noise. The fit produces a continuous function that avoids stability issues arising from supplying highly fluctuating data to the macro solver. A least squares fit is applied to an N th order polynomial for the velocity profile in the core region, and an M th order polynomial for the velocity profile in the constrained region: and where b k,i and c k,i are the coe cients of the polynomials used in the core micro region and constrained region respectively. An estimate of the new slip velocity u B for input to the macro solution (6) is taken directly from the compressed wall micro-element solution (16), at y 0 i = 0. Shear-stress least-squares fits in the core and constrained regions use Eqns. (16) and (17) as well, but with one polynomial degree less, i.e.
Values for N and M are taken to be 3 and 6 respectively for all our hybrid simulations.

Stress correction (Step 6)
The stress correction fields in both core i,core (y 0 ) and constrained i,cs (y 00 ) regions can be calculated using Eqn. (5) for each micro domain: where angle bracket terms on the right hand side are obtained from Eqns. (16)- (19). The final step is to interpolate the flux correction field from the core region of all micro elements to cover the entire macro domain, i.e. over what we refer to as the 'macro gaps', as shown in Fig.  6. To do this we use a Qth-order polynomial: where g k,i are the coe cients of the polynomial, h i,gap is the distance between the core region of two adjacent micro elements, and y 000 i is a gap co-ordinate system as shown in Fig. 6. The order of the polynomial a↵ects the overall noise characteristics of the algorithm, and the optimal choice depends on the case being simulated, as will be discussed later. The new stress correction field is now a continuous function (y) that is used in the macro solution of the next iteration (i.e in Eqn. 6).

Iteration and convergence
Updated continuum solutions u(y) are therefore obtained at each iteration of the algorithm, until convergence u l 1 ! u l , where l is the iteration index. Here we use a simple convergence expression based on solutions of u of the following form: where N P are the number of macro nodes and ⇣ tol is a pre-specified tolerance parameter.
A A ' Figure 6: Illustration of macro gaps (dark regions) over which the stress correction is interpolated.

Validation test for a Newtonian fluid
The HMM-FWC algorithm is now tested on the Poiseuille flow problem of a Lennard-Jones fluid, such that the choice of state point and maximum strain rate gives a Newtonian response. The capability of the method is shown by intentionally choosing the viscosity that is input to the macro solution to be ⇠40% smaller than the actual viscosity at the target density and temperature of the problem. We show how our method corrects the mismatch of viscosity (i.e. that the micro solution corrects the macro solution) indirectly through the stress-correction field.
The Poiseuille channel height, L=162 nm, is chosen small enough in this test case so that results from our hybrid simulation can be compared and validated against a full MD simulation (which is very computationally expensive). Computational savings over the full MD will be far greater for larger and more realistic nano and micro flow systems.
The description of the fluid (state, potentials, thermostat) and wall (packing, lattice-structure, potential) need to be the same for both full MD and hybrid simulations for a fair comparison of accuracy to be made. The potential properties for the liquid-liquid and wall-liquid interactions are taken from [22], with the intention to generate slip at the solid-liquid interface. The values for these are, l l = 3.4 ⇥ 10 10 m, ✏ l l = 1.65678 ⇥ 10 21 J and w l = 2.55 ⇥ 10 10 m, ✏ w l = 0.33 ⇥ 10 21 J in Eqn. (15). The solid wall is chosen to be a simple cubic lattice structure, with number density 5.5 ⇥ ⇢ n , where ⇢ n = 1437 m 3 is the fluid number density. The mass of liquid atoms is m l = 6.6904 ⇥ 10 26 kg. In the full MD description, there are 103,593 liquid molecules and 9,464 wall molecules. The macro force f x = 4.872 fN and Berendsen thermostat at T = 292.8K are applied to the equations of motion of all molecules until a steady-state solution is reached (after ⇠ 300 ns). The steady-state full velocity profile is then averaged over 2 ns. The hybrid setup of the channel is shown in Fig. 7 and consists of the corrected Navier-Stokes description, Eqn. (6), with an under-estimated viscosity µ est = 0.137 ⇥ 10 3 kg/ms, that is solved numerically over 4801 grid nodes on one side of a symmetric nano-channel. The macro spacing is x macro = 17 ⇥ 10 12 m. Three separate hybrid simulations with varying numbers of micro elements (⇧ = 1, 2, 3) are run for the same case. For all these MD simulations we choose the dimensions h core = h cs = 6.8 nm and x = z = 5.44 nm, where the latter is exactly four times the Lennard-Jones cut-o↵. The wall description in the micro elements are identical to that in the full MD simulation, including the number of molecules and the potentials. There are 12,823 and 8,645 liquid molecules in the wall and bulk micro elements respectively. The gaps between micro elements for the three cases are displayed in Fig. 7. For this case, a value of Q=2 is su cient for the interpolant polynomial in Eqn. 21 because the flow is Newtonian in the whole computational domain, so we expect a linear stress-correction. A higher polynomial order can also be used, but this increases noise-related oscillations in the overall result. A tolerance value of ⇣ tol = 0.01 was chosen for this test case to determine convergence, which was obtained in all three hybrid simulations after only 3 iterations. The converged results for velocity profiles from the hybrid simulations are presented in Fig. 8. Velocity profiles from all three hybrid cases match very well with the steady-state velocity profile measured from the full MD simulation. Also plotted in the figures are the velocity profiles predicted by the no-slip Navier-Stokes solution (with the initial estimated viscosity µ est ). We confirm from these results that the 40% lower viscosity that is chosen as an initial guess produces a 40% error in velocity in the middle of the channel, but does not a↵ect the final hybrid solution.
We present the stress correction field (y) (as described in §2.7) from these multiscale simulations in Fig. 9. In this example, a non-zero is generated entirely due to the estimated viscosity being inaccurate. Comparison with a full MD simulation for is not possible because stress correction applies only to the multiscale method. However, in order to provide an indirect calculation we input the stress h⌧ xy i F and velocity hui F measurements obtained from the full-scale simulation inside Eqn. (20) to obtain . In all three hybrid cases there is good agreement with the full MD estimate of , and the fields are approximately linear. The linear profile is expected, given that the simple LJ fluid we chose in this example remains Newtonian. Fluctuations in the calculation of in the hybrid simulations stem from the noise inherent in the micro MD simulations.
In Fig. 10 we present an overview of the most important accuracy and convergence results of our hybrid simulations. No major instabilities in the hybrid algorithm were observed when simulations were run for longer than the converged iteration l=3. Fig. 10(a) shows the convergence parameter ⇣, which is calculated from the velocity profiles at two successive iterations using Eqn. (22); after the third iteration, all three cases lie below the pre-specified tolerance value ⇣ tol = 0.01. In Fig. 10(b) we present a measure of accuracy for the overall velocity profile for all three hybrid simulations when compared to the full-scale simulation, using a root-mean-square (R.M.S.) estimate at each iteration step: where u F , u H are the full-MD and hybrid velocity solutions respectively, and N p are the number of macro nodes. The predicted R.M.S. errors are very close to zero for all cases. Similarly in Fig. 10(c) we present the percentage errors and show that these are less than 1% for ⇧=2,3 and less than 2% for ⇧=1. There are only marginal di↵erences between the three hybrid cases, and so we can be confident that the solution of the Poiseuille flow of a Newtonian fluid could be obtained using just one wall micro element (⇧=1), with additional elements not improving the accuracy of the prediction substantially. Note that in conventional HMM (using point-wise coupling) the minimum number of subdomains is equal to the number of grid points in the macro domain (here, it would be over 4000). The slip velocity u B in Eqn. (6) and the velocity in the midpoint of the channel are presented in Fig. 10(d) and (e) respectively. Again, there is excellent agreement between the full-scale MD simulation and the hybrid results for all three cases.
To verify that we get an accurate prediction of viscosity from the hybrid simulation, we assume that the stress correction has the same form as the Navier-Stokes relationship, i.e., whereμ is the missing contribution of viscosity. Eqn. (24) can be numerically integrated to calculate a corrected viscosity: which is then plotted against the viscosity measured in the full MD simulation in Fig. 10(f). Again, there is good agreement in the viscosity prediction µ corr between our hybrid approach and the full MD solution.

Validation test for a non-Newtonian fluid
There are many applications in the polymer processing industry where accurate predictions of flows of various rheologically-complex fluids, such as polymer solutions and polymer melts, are important [12,26]. One of the main challenges associated with modelling these fluids is that even if the non-linear relationship between shear-stress and strain-rate is known, we are still left with a system of partial di↵erential equations which are often too complex to solve numerically [26]. An alternative to continuum fluid modelling is to adopt a model of the fluid's molecular constituents, and solve the flow using non-equilibrium MD simulations [26,27].
In this section we demonstrate the capability of our multiscale method for a polymer fluid. The method is properly tested by choosing operational strains in a regime where the relationship with shear-stress is highly non-linear, thus imposing a non-linear variation of viscosity throughout the entire computational domain.
The test fluid in this case consists of a pseudoplastic polymer that is modelled artificially in MD using dumbbell chains of N B connected 'beads' (that replace the notion of atoms) [26]. Every bead in the chain is connected to another via a finitely extensible nonlinear elastic (FENE) potential, which is expressed as: where k = 0.43 kg s 2 is the FENE spring constant and r 0 = 5.1Å is the finite extensibility of a pair of connected beads. These properties are chosen from Refs [12,28] in order to minimise bond cross-overs. Here we choose a fluid with N B = 10 beads per chain, at a bead density of ⇢ n = 1447 m 3 and temperature T = 292.8 K. A Lennard-Jones potential (15) is also applied in a similar way as before between all neighbouring beads within the cut-o↵ distance.
The same Poiseuille flow case as in §3.1 is used here to test our multiscale method for non-Newtonian flows. Properties of the setup remain the same as before, except when explicitly stated. The dimensions x = z = 6.8 nm are made larger for both the full MD and hybrid cases in order to allow enough distance for one chain not to interact with itself across the periodic boundaries if fully extended (⇠4.6 nm). As a result, the number of molecules in the full MD was increased to 166,274 (with 2,592 wall molecules), 22,666 molecules in the wall micro element, and 13,526 molecules in the bulk micro element. The macro force was also increased to f x = 24.36 fN in order to operate the fluid in a non-linear stress-strain regime, as discussed above. To control the level of slip at the wall we have reduced the wall density to the same as the fluid, and modified the wall-liquid LJ parameters to w l = 3.4 ⇥ 10 10 m, ✏ w l = 0.995 ⇥ 10 21 J, as suggested in [22].
The full MD case is set up by first initialising the channel height with the target number of beads corresponding to the target density, followed by a short simulation until all chains are grown to the correct bead-length. The gravity-type forcing, i.e. f x = 24.36 fN is then applied until a steady-state is reached (after ⇠130 ns), followed by a short simulation period for measurement (lasting ⇠1 ns).
For the multiscale approach, we run three hybrid simulations, as before, for ⇧ = 1, 2, 3 micro elements. The estimated viscosity input to the macro solution, µ est = 0.683 ⇥ 10 3 kg/ms, is taken at very small strain rates (see Fig. 12(b)), and is initially uniform throughout the domain. Our multiscale method corrects for this approximation through the stress correction field , from which a more accurate variation of e↵ective viscosity can then be calculated. The variation of is expected to be non-linear so the polynomials in the macro gaps are of degree Q=4, as defined in Eqn. (21). All hybrid simulations converge in 3-4 iterations, for a tolerance of ⇣ tol = 0.02. The channel velocity profiles are presented in Fig. 11(a), while the overall error plots are displayed in Fig. 11(b). Excellent agreement is found between the full-scale MD simulation and the hybrid simulations for ⇧=2 and ⇧=3 micro elements (less than 1%). However, less good agreement is found for the hybrid case with just one wall micro element (⇧=1), which oscillates between 2-7% of the full scale result; inaccuracies are caused by the lack of micro resolution in the bulk, whereas near-wall predictions (such as slip and stress) match just as well as the higher ⇧ cases. This can also be observed in the overall stress-correction field in Fig. 11(c), which, unlike before, now varies non-linearly. Note that there is no actual stress correction in the full MD simulation, but as before, we evaluate it using the equivalent of Eqn. (20), using the shear stress h⌧ xy i F and velocity hui F profiles measured from the full MD simulation.
To evaluate an 'e↵ective viscosity' field µ e↵ for both hybrid and full scale simulations, we use the Oswald-de Waele relationship: The full MD viscosity field is evaluated directly from this expression using the stress h⌧ xy i F and velocity hui F measurements. For the hybrid simulation, Eqn. (27) must first be substituted in the corrected Navier-Stokes Eqn. (5) to give: Calculations of e↵ective viscosity are plotted in Fig. 11(d). Again, there is good agreement between the full scale MD and hybrid calculations adopting ⇧=2,3 micro elements, but poor agreement for the ⇧=1 case in the bulk of the channel. Good agreement is also obtained with the standard power-law relationship: where a is the flow behaviour index (a < 1 for pseudoplastic fluids) and K is the flow consistency index. The values for a and K are obtained using both bulk and near-wall MD (pure) simulations of the same dumbell fluid and thermodynamic state, under various strain-rates. E↵ective viscosity vs strain-rate and shear-stress vs strain-rate are plotted in Figs. 12 (a) and (b), respectively. In these simulations we obtained fit values for a = 0.275 and K = 5.3 ⇥ 10 7 Pa s a .

Computational savings
In this section we present the computational savings that our HMM-FWC method has a↵orded over the full-scale MD simulations, by computing the total processing time P = total number of time-steps ⇥ average processing time per time-step, for both approaches. For the full MD simulation this is simply: where t c,F is the average clock time per MD time-step, t tran,F is the transient-time required for the flow to develop from rest to steady-state, t av,F is the averaging-time required to reduce noise from molecular measurements and t is the MD time-step. For the hybrid method, it is: where subscripts W and B refer to the wall and bulk micro elements respectively, and I is the number of iterations at which the hybrid algorithm converges (I=3 in §3.1 and I=4 in §3.2). The cumulative averaging times are the same for all MD cases, so t av,F = t av,W = t av,B = t av (t av ⇡ 9 ns in §3.1, while t av ⇡ 2 ns in §3. of the MD domains. For example, in the Newtonian case §3.1, the transient time in the hybrid simulation is t tran,W = t tran,B = 2 ns, while the full MD su↵ers from a very long start-up time t tran,F ⇡ 300 ns. 4 Not having to simulate this long transient period is a major computational advantage of our multiscale approach for studying steady-state flows. In the non-Newtonian case t tran,W = t tran,B = 1 ns, while t tran,F ⇡ 130 ns. The ratio of Eqns. (30) to (31) gives the total computational speed-up factor, c = P F /P H . For the cases in §3.1, the average clock times measured for the wall and bulk micro elements were: t c,W =0.0326 s and t c,B =0.015 s respectively, while the full MD simulation measured a clock time of t c,F =0.251 s. So this means 74⇥, 50⇥ and 38⇥ speed-ups for the hybrid cases with ⇧ = 1, 2 and 3 micro elements respectively. For the cases in §3.2, the average clock times were: t c,W =0.059 s, t c,B =0.033 s and t c,F =0.595 s. This means 100⇥, 66⇥ and 49⇥ speed-ups for the non-Newtonian cases with ⇧ = 1, 2 and 3 micro elements respectively (I = 4). Clearly our hybrid method is very computationally e cient.

Conclusions
We have proposed a multiscale method that di↵ers from the Heterogeneous Multiscale Method (HMM) [13] by employing a field-wise macro/micro coupling (as opposed to a point-wise coupling at every macroscopic node). This enables a greater range of scales over which our method can be applied, which is essential for the modelling of practical nano and micro flow systems. The ability to decouple the positions of the micro elements from the macro grid nodes a↵ords major computational savings (in both examples presented here, excellent results were obtained for as few as two micro elements). The field-wise coupling method is implemented using a novel macro/micro constraint procedure that involves calculating stress-corrections to the standard Navier-Stokes equations, and projections of flowfields of arbitrary complexity on bulk-and wall-type micro elements.
We have simulated both Newtonian and non-Newtonian Poiseuille flow, and the multiscale method gives flow predictions in close agreement with the full-scale MD simulation of the same case. The method converged very quickly, within 3-4 iterations, which is an attractive feature for extending it to study quasi-steady transient problems, as described in [29]. While the Newtonian flow case showed that having just one wall micro element is su cient to produce accuracies within 1-2% of the full-scale result, we found that in the non-Newtonian case at least one bulk micro element and one wall micro element were required to capture the significant variation of viscosity produced by a FENE fluid under high strain-rates. The multiscale method is approximately 74⇥ and 66⇥ more computationally e cient than the full-scale simulation for the Newtonian and non-Newtonian cases respectively.
The HMM-FWC method presented in this paper can be extended to two-and three-dimensional flow problems. However, the method becomes more complicated and there are some particular issues that need to be addressed carefully. In 3D flows, for example, MD micro elements are constructed with a constraining region that surrounds the core region completely (instead of the top-bottom approach in 1D flows), and where the micro constraint forces are a 3D varying field. The accurate application of this force field is paramount, and thus a fine 3D mesh is required, which is computationally costly. Similar issues arise for the macro-to-micro projection step and the micro compression step -high order interpolants are required to provide velocity projections and fits through measurements of shear-stress and velocity taken from the MD core region. In the projection step there is also the issue of preserving continuity at the core-constrained region (which is not necessarily straightforward). Furthermore, in the MD simulations there is the exacerbated problem of noise, since smaller bins are required to capture 3D gradients, which may lead to instabilities in the hybrid algorithm if the measured fields are not su ciently averaged. For the one-dimensional examples considered in the current paper, mass is automatically conserved, leaving only the momentum conservation Eqn. (2) needing to be solved macroscopically. For more general flow problems, a macroscopic continuity equation must also be solved alongside an equation of state relating density to pressure. The induced pressure variations would then be coupled to the momentum conservation equation using a distributed force field. The development of the HMM-FWC method to higher-dimension flow problems is the subject of future work.