Enhancing nano-scale computational fluid dynamics with molecular pre-simulations: unsteady problems and design optimisation

We demonstrate that a computational ﬂuid dynamics (CFD) model enhanced with molecular-level infor- mation can accurately predict unsteady nano-scale ﬂows in non-trivial geometries, while being efﬁcient enough to be used for design optimisation. We ﬁrst consider a converging–diverging nano-scale channel driven by a time-varying body force. The time-dependent mass ﬂow rate predicted by our enhanced CFD agrees well with a full molecular dynamics (MD) simulation of the same conﬁguration, and is achieved at a fraction of the computational cost. Conventional CFD predictions of the same case are wholly inade-quate. We then demonstrate the application of enhanced CFD as a design optimisation tool on a bifurcat- ing two-dimensional channel, with the target of maximising mass ﬂow rate for a ﬁxed total volume and applied pressure. At macro scales the optimised geometry agrees well with Murray’s Law for optimal branching of vascular networks; however, at nanoscales, the optimum result deviates from Murray’s Law, and a corrected equation is presented. (cid:2) 2015 The Authors. Published by Elsevier Ltd. ThisisanopenaccessarticleundertheCCBYlicense(http://


Introduction
Many emerging applications of nanofluidic technology take advantage of different physical effects that dominate at small scales; examples can be found in air and water purification [1,2], and in micro chemical reactors [3,4]. The design of these technologies would be greatly facilitated by being able to perform numerical simulations that predict mass flow rates and heat transfer. Computational fluid dynamics (CFD) is regularly used to model and create optimal every-day engineering designs efficiently. However, the assumptions used to derive the continuum fluid equations become invalid in highly-confined systems, making the equations inaccurate. On the other hand, molecular dynamics (MD) can be used to perform highly detailed simulations of nano-scale systems; it has been successfully used to study the behaviour of protein folding [5], crystal formation [6] and chemical reactions [7]. The drawback is that MD is extremely computationally intensive, especially when used to model systems comprising hundreds of thousands of molecules that would be required for engineering applications.
The continuum fluid assumptions become inaccurate for gas flows as the smallest characteristic scale of the geometry (e.g. channel height) approaches the mean distance between molecular collisions (i.e. the mean free path) [8]. When modelling dense liquids (as we do in this paper) there is not a well-defined condition for when the fluid assumptions become inaccurate. However, it appears that they fail when water is confined in channels of width $1-2 nm (see [9] and references therein), and MD simulations have been used to show that Lennard-Jones fluids confined in geometries of $2-3 nm still show continuum behaviour [10][11][12]. At the nano-scale the fluid molecules form layers parallel to an interface, which causes the strain rate to vary rapidly within several molecular diameters [13]. These large variations mean the stress no longer has local linear behaviour [14,15].
Despite the complexities of fluid behaviour at the nano-scale, it has recently been shown that useful predictions from CFD can be obtained for some simple geometries, if appropriate fluid state models, viscosity relationships, and slip models are extracted from an MD pre-simulation [16]. In [17], a similar approach is used to obtain CFD predictions of flow through a nanotube; with results agreeing well with full MD simulations.
However, what remains to be tested is the robustness of nanoscale CFD when applied to more complex engineering calculations. In this paper, we test if our enhanced CFD model is robust enough to predict flow behaviour in a non-trivial geometry (a convergingdiverging channel), while using various forms of applied forcing to generate unsteadiness within it. As a demonstration of its efficiency, we go on to apply the enhanced CFD to the design optimisation of a small fluidic network.
The paper is structured as follows. In the next section we summarise the MD pre-simulation procedure and the CFD model used. We then use these models to perform unsteady simulations of a converging-diverging channel, where the width of the channel is close to the continuum limit. Our model is then applied to simulate an industrially-relevant problem of flow through a bifurcating channel. We show that to optimally design the channels the slip velocity at walls must be taken into account. The paper then ends with a summary.

The MD pre-simulations
As in [16], we employ preliminary molecular dynamics (MD) simulations to obtain fluid properties and boundary conditions that enable the effective use of a Navier-Stokes fluid solver for nano-scale applications. This approach can be classed as a ''sequential molecular-continuum hybrid method'' (see [18] for a review of hybrid methods), where 'sequential' refers to the fact that the MD is performed in advance of, and so independent of, the continuumfluid solver. Fig. 1 (far left) shows a schematic of the MD pre-simulation domain; it is symmetrical about its centrelines and uses periodic boundary conditions in the streamwise direction (in the x-direction) and into the page (in the z-direction). The domain has bulk, shear and interface zones (as labelled) for measuring state, constitutive and boundary properties, respectively. Pressure and density are measured in the bulk zone. In addition to this, in the bulk zone an artificial streamwise body force ðF x Þ is applied (Fig. 1, centre left), which creates a velocity profile in the domain similar to that illustrated (Fig. 1, centre right). We assume that the equation of state in the bulk zone is unaffected by the magnitude of strain rate generated. In the shear zone the fluid is, therefore, subject to a constant shear stress, s xy , directly resulting from the bulk-zone forcing. A linear flow velocity profile is developed in the shear zone, and this is least-squares fitted to obtain a strain rate and shear viscosity coefficient, l. Any significant density oscillations associated with molecular layering are confined to the interface zone. In this zone we calculate what we term the 'CFD surface displacement', d, which is the distance that a CFD wall/surface needs to be displaced from the centres of surface atoms in order to accurately represent the boundary of the fluid (as opposed to the boundary of the solid); see d in Fig. 1. We take this displacement to be the distance from the centre of the surface wall atoms to where the fluid density becomes at least 10% of the bulk, i.e. q P aq bulk , where a ¼ 0:1.
Note, the surface displacement is quite insensitive to the percentage of the bulk density chosen as the threshold, since the density increases from zero to well above the bulk density over a very short distance. For example, had we chosen the threshold to be at 20% of the bulk density, the surface displacement would have only been 1-2% larger, for a typical case.
The linear velocity profile obtained in the shear zone is extrapolated into the interface zone to find the apparent slip length, n, as defined from the CFD surface (see Fig. 1, centre right).
The molecular dynamics pre-simulations, and the full-scale MD simulations used for benchmarking, are performed using the mdFoam solver [19][20][21][22] that is implemented within the OpenFOAM libraries [23]. For the test cases considered in this paper we have adopted a simple Lennard-Jones (LJ) fluid model (at 292.8 K), where the solid LJ wall atoms are fixed/frozen [24]. However, the methodology is general to any given molecular model. For full details of the molecular pre-simulation domain and the molecular dynamics parameters used, the reader is referred to [16].  Fig. 1. Schematic of molecular dynamics pre-simulation for extracting fluid dynamic properties that are essential inputs to an enhanced CFD solver for nano-scale flows. Fig. 2(a) shows MD pre-simulation measurements of pressure, obtained from the standard Irving-Kirkwood expression [25], varying with the mass density. The MD pre-simulation results are leastsquares-fitted to a 2nd order polynomial. This then serves as an equation of state within the enhanced CFD solver to connect the mass continuity equation to the momentum equation. In this case the polynomial is p ¼ 0:001559q 2 À 3:387q þ 2020:6.

Line of symmetry
The strain-rate is extracted from the MD shear zone by a leastsquares linear fit to the relaxed and time-averaged velocity profile. The applied shear stress is measured using the Irving-Kirkwood equation, from which we obtain a dynamic shear viscosity coefficient for the LJ fluid at a given bulk density. The viscosity coefficients measured from our MD pre-simulations of Lennard-Jones argon are shown in Fig. 2(b). A least-squares polynomial fit of 2nd order in density is also plotted: l ¼ 7:96 Â 10 À10 q 2 À 1:774 Â 10 À6 q þ 0:001106. This is then used in our enhanced CFD simulations to close the momentum equation. Note, due to the breakdown of the continuum assumption and the existence of nonlocal stress, this state-dependent viscosity becomes only approximate when applied to a nano-confined fluid.
The surface displacement d defines the location of the CFD boundaries relative to the atomic (actual) walls. If d varies substantially with density (or any other fluid property), the geometry of the enhanced CFD domain becomes dependent on the CFD solution itself. However, for the fluid/solid combination considered in this paper, over the density ranges considered, d is effectively constant, as seen in Fig. 3.
In certain cases the value of d will itself be dependent on the geometry, particularly for high curvatures, such as around sharp corners and obstructions. It is beyond the scope of the current work to attempt to accommodate these influences, while noting that, later, we obtain good agreement with full MD simulations without doing so. To tackle geometry-dependent flow properties (including surface displacement) would dramatically increase the parameter space that the pre-simulations would be required to supply information for; in fact, for such problems a 'concurrent' hybrid approach is likely to be more efficient.
As the spatial-scale of the geometry increases, the relative significance of the surface displacement reduces. We can develop a simple gauge of its impact by considering the percentage that it modifies the mass flow rate in a simple channel in two limiting cases: assuming no-slip at the walls (i.e. n ¼ 0); and for very high slip (i.e. n ) h, where h is the channel width). In the no-slip case, for Poiseuille flow, the mass flow rate is proportional to the cube of the channel width; the percentage difference of using the surface displacement is then For the cases in Section 3, where the channel width varies, is $28-44%. For high-slip cases, where the velocity profile becomes pluglike, the mass flow rate becomes proportional to the square of the channel width, giving a percentage difference: Considering again the cases in Section 3, would be $20-32%; i.e.
the impact of the surface displacement is likely to be very significant regardless of the degree of velocity slip. Based on the estimates of Eqs. (1) and (2), the impact of a surface displacement d $ 0:2 nm will only be less than 1% (i.e. negligible) for channels greater than 75-100 nm. Liquid slip velocity at surfaces is calculated using the Navier slip condition: where n is the slip-length and _ c is the shear-rate at the bounding surface. The least-squares-fitted linear velocity profile is used to calculate the slip-length (as defined from the CFD surface). Based on the strain-rate/slip-length relationship proposed in [24], and assuming a linear dependence on density, a least-squares fit is performed to the following equation: where q is the density, _ c c is the critical shear rate (see [24]), and   3) and (4) is directly introduced as a Robin boundary condition in the enhanced CFD solver.

The enhanced CFD model
We use the laminar, compressible flow solver sonicLiquidFoam, which we have modified to (a) accommodate a nonlinear equation of state, (b) allow a density-dependent viscosity, and (c) incorporate slip boundary conditions of the form given in Eq. (4). A compressible solver is used despite the very low Mach numbers, since significant compressibility can occur in micro and nano geometries due to very high viscous pressure losses [27,28].

Unsteady simulations
We now simulate the unsteady flow behaviour of a Lennard-Jones fluid along a converging-diverging channel; a case chosen to demonstrate the robustness of the enhanced CFD model when applied to non-trivial flow problems.
Owing to the lack of detailed and reliable experimental flow measurements at the nano-scale, in this section we compare our enhanced CFD predictions with full-scale MD simulation results. This comparison is intended to test whether enhanced CFD can produce flow field solutions of comparable accuracy to full MD in complex nano-scale geometries, without the need for ad hoc corrections, and at only a fraction of the cost of full MD.

Cases
We consider a two-dimensional geometry: a convergingdiverging channel with a smoothly varying height in the streamwise direction. A gravity-type force is applied to the fluid to generate an unsteady/transient flow. As test cases, we choose flows that exhibit non-continuum behaviour (e.g. slip at surfaces), and do not contain a significant bulk flow region, i.e., the width of the channel is at the 2-3 nm continuum-fluid limit for a Lennard-Jones fluid.
The converging-diverging channel is shown in Fig. 5 and has a length l ¼ 68 nm in the streamwise direction x, a depth of 5.44 nm, and heights of 3.4 nm and 2.04 nm at the inlet and throat sections, respectively. The channel is periodic in both the streamwise and spanwise direction. The height between top and bottom walls hðxÞ varies in the streamwise direction according to a sinusoidal function, where 4a ¼ 1:36 nm is the change of height from inlet to throat, and h inlet is the height of the channel at the inlet. The full MD domain is divided into 200 bins in the x-direction of bin-width dx ¼ 0:34 nm, and the instantaneous mass flow rate and density are measured in each bin. In the enhanced CFD domain, we define a plane across the channel at equivalent positions, and sum the mass fluxes from each cell the plane crosses, at each time-step, to get the instantaneous data. Dependency studies on the mesh resolution and on the time step showed that 50,000 cells and a time step of 21.6 fs were more than sufficient to obtain converged CFD solutions.
All the flows start from rest, then a time-varying gravity force F g ðtÞ is applied. We consider four different forces applied to the fluid: where t is the simulation time; 3. Long oscillations: an unsteady oscillating gravity force of the same amplitude, but with a larger period T ¼ 10:8 ns, i.e.

Results for the four cases
To test the reliability of our CFD predictions that have MD presimulation input, we compare results with full-domain molecular dynamics calculations (referred to as 'full MD'). To test whether our enhanced CFD model is an improvement over conventional CFD, we also compare results with predictions from compressible CFD with no-slip at the wall and without incorporating a CFD surface offset (referred to as 'no-slip CFD'). We also compare with incompressible CFD with slip incorporated but no surface displacement (referred to below as 'incomp. slip CFD').
In Fig. 7 we plot the mass flow rate variation with time in a single bin near the inlet of the channel for each case. We see that in all cases the enhanced CFD model is able to accurately predict how the mass flow rate changes in time. Fig. 7(a), in which a constant force is applied throughout the channel, shows that the CFD reaches steady state at the same time as the MD simulation, and that a similar final mass flow rate is reached. There are, however, substantial differences between the enhanced CFD, the no-slip CFD, and the incomp. slip CFD results. The oscillations that are observed at the early times in Fig. 7(a) in the enhanced CFD results and also the MD data are due to an acoustic response of the nano channel to impulse forcing. A first estimate of the natural acoustic period is obtained by T ¼ l=c ¼ 0:07 ns (where c is the speed of sound). This corresponds reasonably closely with the observed kinks in the mass flow rate. Fig. 7(b) shows the results when an oscillating force with period 0.22 ns is applied. The mass flow rate in the enhanced CFD oscillates with the right frequency, the correct amplitude, and is also in phase with the full MD results. The no-slip CFD, on the other hand, appears to have the correct frequency but the amplitude is incorrect and it is oscillating out of phase, while the incomp. slip CFD is in phase but overpredicts the amplitude.
In Fig. 7(c) we have an oscillating force with period 10.8 ns. The mass flux in the enhanced CFD oscillates with the right frequency, correct amplitude, and is in phase, whereas the no-slip CFD appears to have the correct frequency, and is oscillating in phase, but the amplitude is incorrect. In Fig. 7(d) the period of the oscillating force increases from 0.22 ns to 10.8 ns; even in this more elaborate case, the enhanced CFD prediction is accurate. Table 1 provides an indication of the computational cost for the full-domain MD simulations. The longest simulations presented in this paper ran in parallel (on 24 CPUs) for 48 days. The enhanced CFD itself has negligible computational cost by comparison, although the MD pre-simulations require the computational resources indicated in the last row of Table 1. However, these pre-simulations need only be performed once for a particular fluid/solid combination, and then can be used for any number of flow geometries thereafter.

Design optimisation
We now demonstrate how the enhanced CFD model can be used in design optimisation problems at the nanoscale. The example we choose is the optimal design of a bifurcating nano-channel network (see Fig. 8); such a design exploration would not be feasible using full MD simulations. The problem is to find the optimal widths of the channels in a bifurcating channel (i.e. those that give greatest mass flow rate), for a constant pressure difference Dp between the inlet and the outlets, and a constant volume V. At the macro scale the solution to this problem is given by Murray's Law [29,30], which was first derived using the Hagen-Poiseuille Law to minimise the power required to sustain the flow of blood through vessels. It has also since been shown to describe the water transport though biological vessels in plants [31], and at the micro scale it has been used to optimally design MEMS devices with rectangular or trapezoidal cross sections [32]. For a 2D two-level network, like the geometry in Fig. 8, Murray's Law is where h 0 is the width of the inlet parent channel, and h 1 to h N are the widths of the outlet daughter channels. For a symmetric bifurcating channel with N ¼ 2 and h 1 ¼ h 2 , the optimum ratio of channel widths is then given by (d)   The optimisation we perform, with the given constraints of constant volume and fixed pressure difference, is a linear search on channel width (equal increments) to find the maximum mass flow rate. We choose a volume of 1100 nm 3 , a pressure difference of 10 MPa, channel lengths l 0 ¼ l 1 ¼ 75 nm, a junction length l j ¼ 20 nm and junction width, h j ¼ 4 nm. The volume V can be calculated as: If this geometry was optimised using MD, each simulation would take approximately 30 days, whereas each enhanced CFD simulation takes approximately 500 s to perform. Fig. 9 shows the results from this optimisation with our enhanced CFD model used on a micro-scale channel and on a nano-scale channel. We see that for a micro scale channel, the optimum width occurs when h 2 0 =2h 2 1 ¼ 1: this is the expected result according to Murray's Law. At the nano-scale, however, we observe a significant deviation from the standard Murray's Law, which is now discussed.
A deviation from the standard Murray's Law has been noted for rarefied gases [33] but has not so far been demonstrated for a liquid. To uncover the origin of this deviation we derive Murray's Law using Poiseuille's equations with Navier slip at the walls, i.e. uðhÞ ¼ uðÀhÞ ¼ n du dy where n is the slip length, and the velocity is at a maximum at y ¼ 0 i.e. du dy j y¼0 ¼ 0. The mass flow rate is then where _ m is the mass flow rate, q is the density, l is the dynamic viscosity and Dp is the pressure difference between the inlet of the parent channel and the outlet of the daughter channel. Murray's Law is found by minimising the power P required to maintain flow, which for flow through a channel is where b is a constant of proportionality. By eliminating Dp with Eq. (11) in this equation and differentiating, we find that when the power is minimised the mass flow rate is where k ¼ 2=3 ffiffiffiffiffiffiffiffiffiffiffi ffi qb=l p . For a symmetric bifurcating channel, the mass flow rate through the parent channel must equal the total mass flow rate through the daughter channels, i.e. _ m 0 ¼ 2 _ m 1 , therefore, the optimal ratio of channel widths becomes h 2 0 It is clear that when h 0 ; h 1 ) n this becomes Eq. (9), as expected. It can also be noted that when the flow becomes plug-like, i.e. h 0 ; h 1 ( n, this ratio becomes h 2 0 =2h 2 1 ¼ 2 1=3 . We can now use Eqs. (10) and (14) to calculate the expected value of h 2 0 =2h 2 1 . When comparing this to the optimum found by the enhanced CFD we get excellent agreement, as highlighted in Fig. 9. This shows that the slip at the walls is the important factor in the deviation from the expected optimum. A CFD model that includes an accurate model of the wall-fluid interaction is, therefore, potentially very important in the design of nano-scale devices.

Summary
We have shown that a CFD model enhanced with data from MD pre-simulations is capable of making accurate predictions of unsteady liquid flow along a converging-diverging channel that has a width close to the expected continuum-fluid limit. This enhanced CFD approach is far more accurate than conventional CFD calculations, and significantly more computationally efficient than full MD simulations.
We have also demonstrated the enhanced CFD approach applied to a design optimisation problem: that of a bifurcating nanofluidic network. The widths of channels in the network should be optimised to maximise the mass flow rate through the network, for a fixed pressure drop and network volume. We have shown that slip at the nano-scale can have a very significant effect on the optimum channel dimensions, and we have derived an analytical equation which corrects the well-known Murray's Law. This is one of many possible cases where nano-scale flow effects modify the optimal design of nanofluidic systems when compared with their macroscopic counterparts.

Murray's Law
Slip adjusted Murray's Law Fig. 9. The normalised mass flow rate for different ratios of mother/daughter channels widths for micro-scale channels (squares) where the dash-dotted line is the expected optimum, and for nano-scale channels (circles) where the dashed line is the expected optimum.