We present new algorithms for simulating pressure gradient flows in molecular dynamics (MD) simulations. Nano-channel inlet and outlet non-periodic boundary conditions are implemented using hydrodynamic state reservoirs and flux boundary models at arbitrary boundaries of the domain geometry. We demonstrate the new method in a complex nano-mixer configuration and for droplet deformation in extensional flow channels. The technique which we propose is applicable to any complex nano-channel configuration, and may serve as a useful tool in engineering design of nano-scale applications

Borg, Matthew Karl

Nicholls, William

Reese, Jason

English

Name not available

MD boundary conditions for pressure gradient flows : nano-mixing and nano-droplet deformation in extensional flows

University of Strathclyde Institutional Repository

Strathprints Institutional RepositoryBorg, Matthew Karl and Nicholls, William and Reese, Jason (2010) MD boundary conditionsfor pressure gradient flows : nano-mixing and nano-droplet deformation in extensional flows.[Proceedings Paper]Strathprints is designed to allow users to access the research output of the University of Strathclyde.Copyright c© and Moral Rights for the papers on this site are retained by the individual authorsand/or other copyright owners. You may not engage in further distribution of the material for anyprofitmaking activities or any commercial gain. You may freely distribute both the url (http://strathprints.strath.ac.uk/) and the content of this paper for research or study, educational, ornot-for-profit purposes without prior permission or charge.Any correspondence concerning this service should be sent to Strathprints administrator:mailto:strathprints@strath.ac.ukhttp://strathprints.strath.ac.uk/MD boundary conditions for pressure gradient flows: nano-mixing, andnano-droplet deformation in extensional flowsMatthew K. Borg*, William D. Nicholls and Jason M. ReeseDepartment of Mechanical Engineering,University of Strathclyde, Glasgow G1 1XJ, UK* matthew.borg@strath.ac.ukABSTRACTWe present new algorithms for simulating pressure-gradient flows in molecular dynamics (MD) simulations.Nano-channel inlet and outlet non-periodic boundaryconditions are implemented using hydrodynamic state-reservoirs and flux boundary models at arbitrary bound-aries of the domain geometry. We demonstrate the newmethod in a complex nano-mixer configuration and fordroplet deformation in extensional flow channels. Thetechnique which we propose is applicable to any complexnano-channel configuration, and may serve as a usefultool in engineering design of nano-scale applications.Keywords: Nanofluidics, Molecular Dynamics, Hy-brid Continuum-Molecular, Boundary Conditions1 INTRODUCTIONPeriodic boundary conditions (PBCs) [1], [2] are byfar the most frequently used boundary condition in molec-ular dynamics (MD) simulations, since they are easyto implement numerically and they conserve mass, mo-mentum and energy at the boundaries. Despite theseadvantages, PBCs have geometry and periodicity con-straints — only simple (typically cubic), symmetric sys-tems close to equilibrium may be simulated. In more re-alistic engineering problems, however, it is necessary tosimulate non-equilibrium MD within arbitrary geome-tries [3], [4] entailing non-periodic boundary conditions(NPBCs) [5], [6]. An important fluid mechanics problemfor the development of nanoelectromechanical systems(NEMS) using MD as the simulation tool, is pressuregradient flows in nano-channels that encompass complexthree-dimensional geometrical features, such as bends,constrictions, and multi-inlet and outlet regions of dif-ferent sizes. The application of the gravitational method[7] concomitant with periodic boundary conditions isambiguous for simulating these cases, and generally in-applicable for applying pressure gradients, mainly owingto the streamwise pressure and density homogeneitieswhich the method imposes. Another fluid mechanicsproblem in which development of NPBCs are importantis in the implementation of multi-scale hybrid methods,for example, at the coupling interface between the con-tinuum and MD sub-domains.2 METHODWe consider a domain of arbitrary geometry, definedby an unstructured polyhedral mesh, as typically usedin finite-volume Computational Fluid Dynamics (CFD).The MD fluid occupies the domain geometry and con-sists of “molecules” that interact through a pair-wisepotential U(rij), where rij = |rij | and rij = ri − rj isthe separation vector between a pair of molecules (i, j).The standard shifted Lennard-Jones (LJ) 12-6 poten-tial and the Weeks Chandler Anderson (WCA) potential[2] are used in this paper. The positions ri and veloci-ties vi of molecules evolve according to standard Newto-nian dynamics, using a molecular time-integration step∆tm = 10 fs.boundary regionnffrcutrcutFlow of mass, momentum and energy at a boundary face, fPMD Mesharbitrary non-periodic boundary, bExternal regioncontrolling zonebuffer zonecontrol cell, PFigure 1: Schematic showing the implementation of ournon-periodic boundary conditions at the inlet of a nano-channel using MD. Boundary conditions consist of im-posing flow of molecular properties in a buffer regionnext to an arbitrary boundary b, in addition to feedback-loop control of macroscopic (continuum) properties inthe adjacent controlling zone.An arbitrary non-periodic boundary of the MD do-main may be defined as an “inlet” or “outlet”. Atan inlet, we construct small controlling regions, definedby a group of cells on the mesh (see Figure 1). Con-trol of hydrodynamic and thermodynamic state prop-erties is established by directly applying perturbationsto molecules residing within an arbitrary control cell.This process is facilitated by the ‘cell-occupancy’ datastructure [3], in which each cell on the MD mesh storesan inexpensive link to occupant molecules. Fields oftarget continuum properties (typically, velocity, densityand temperature) are prescribed to the state controller.Control is then performed using a simple feedback loopcontrol system, in which the sampled property (aver-aged over the preceding time interval ∆t ∼ 1 ps) inan arbitrary control cell is first compared with its tar-get property, and then the ‘error’ is imposed within thesame cell over the subsequent ∆t.To manipulate mass density in an arbitrary controlcell P , the number of molecules is controlled uniformlyover t→ tn = t+∆t:∆NP (t→ tn) = NINT([ρtP (tn)− 〈ρP (t)〉] VPmi), (1)where VP is the volume of cell P , mi is the mass ofthe molecule being inserted, ρtP (tn) is the target massdensity and 〈ρP (t)〉 is the time- and cell-averaged den-sity. The NINT(x) function is required to apply thenearest integer, since only whole molecules may be in-troduced (∆NP > 0) or removed (∆NP < 0). An ar-bitrary molecule is inserted within the cell using theUSHER algorithm [8] so that its initial potential energyis close to the average potential energy of the cell 〈UP 〉,while an existing molecule is deleted based on its closestpotential energy to 〈UP 〉.1Continuum velocity is controlled by applying an ex-ternal force to all molecules i within the control cell atuniform control steps within t→ tnfexti (tk) =∆uP (t → tn)ncont∆tmmi, (2)where t ≤ tk < tn, {k = 1, 2, . . . , ncont}, ncont are thenumber of control steps and∆uP (t→ tn) = λ(utP (tn)− 〈uP (t)〉). (3)In Equation (3), λ ∼ 0.8 is a modifiable parameter thatdefines the rate at which the velocities of molecules oc-cupying cell P are accelerated towards the target value.Control of temperature involves an adaptation of thevelocity-rescaling thermostat [2]. The method scales thevelocities of molecules within the control cell P using ascaling factorχP (tk) =(T tP (tn)〈TP (t)〉)1/(2ncont). (4)1The potential energy of one molecule is given by Ui =12∑Nmolj=1( 6=i)U(rij), where Nmol is the number of molecules withininteraction range.To complete the notion of non-periodic boundaryconditions, inlets include also the flow of mass, momen-tum and energy at the termination edge of the bound-ary. For mass and momentum flow, molecules are in-serted/deleted (using a similar procedure to the densitycontroller) at an arbitrary face f of the boundary tomatch the target mass flux, m˙f = (ρtut)f · Af , whereAf = Af nˆf is the face-normal surface area vector and(ρtut)f is the face interpolated target mass flux den-sity. Energy flux is implemented by applying externalforces on molecules residing within interaction distance(rcut = 0.85 nm) of the boundary so as to include theeffect of the missing molecules beyond the boundary [9],[5], [6].An outlet of a nano-channel is modelled to controlthe flux exiting the system that contributes to a pressuregradient across the domain. A vacuum boundary condi-tion [10], that is, a boundary model that removes instan-taneously the molecules colliding with the boundary, hasbeen shown to produce highly non-linear pressure gradi-ents, as well as uncontrollable (typically, very low) pres-sure values at the outlet. In our method, the boundaryis instead modelled by a diffuse thermal wall, in whichcolliding molecules are reflected back into the domainwith a velocity resampled from a Maxwell-Boltzmanndistribution [11]. The mass flux at the same boundaryis then controlled at uniform control steps using the re-lationship:m˙f (t) = −VZnfaces∆t(ρtZ(t)− 〈ρZ(t)〉), (5)where Z denotes a large spatial zone, placed adjacent tothe boundary, typically covering the entire domain size,and nfaces are the number of faces on the outlet bound-ary; ρtZ(t) is a phenomenological parameter which de-termines the pressure/density gradient occurring withinthe target zone. Additionally to this, an energy fluxboundary model is applied at the outlet.3 RESULTSIn order to test these models, we investigate flowsof liquid argon molecules, which have a characteristiclength scale of σ = 0.34 nm, a characteristic energy ofǫ = 120kb = 1.65678× 10−21 J, where kb is the Boltz-mann constant, and a mass m = 6.69×10−26 kg. A cut-off radius of rcut = 2.5σ = 0.85 nm is used. We presentour results in reduced units: time, t∗ = t√ǫ/mσ2, num-ber density ρ∗ = ρσ3, temperature T ∗ = T (kb/ǫ) andvelocity u∗ = u√m/ǫ.3.1 Simple 2D Poiseuille flowWe verify the new inlet and outlet boundary con-ditions by applying them to a simple two-dimensionalisothermal MD flow channel. The nano-channel of di-mensions (Lx = 175σ, Ly = 30σ, Lz = 20σ), consists ofa top and bottom MD face-centre cubic wall of thick-ness 5σ, and periodic boundary conditions in the cross-channel direction. The inlet region is set to control prop-erties uniformly at values of ρ∗ = 0.6, T ∗ = 2.4 andu∗ = (0.5, 0, 0), while at the outlet we choose a targetdensity of ρ∗Z = 0.5. Results show a linear pressure anddensity gradient along the channel, as seen in Figure 2. 0.4 0.5 0.6 0.7 0  50  100  150DENSITY  (ρ*)STREAMWISE DIRECTION  (x*)inletlinear fit 1 1.5 2 2.5 3 3.5PRESSURE  (p*)linear fitFigure 2: Centre-line profiles of pressure p∗ and densityρ∗ measured from the Poiseuille flow MD simulation.3.2 3D complex mixing channelWe now consider a micro-scale mixing channel ofthree-inlet one-outlet design, taken from Hertzog et al.[12], and reduce its scale to nanometers so that reason-able mixing time-scales may be simulated using molec-ular dynamics. Our modified nano-mixer geometry isconstructed using the following route, which we gener-ally adopt for other complex geometry cases. Initiallythe part is drawn in Pro/ENGINEERr, a commercialcomputer-aided design (CAD) drawing application, andexported in STEP format (.stp) to GAMBITr (a mesh-ing utility normally used for FLUENTr CFD). Thegeometry is meshed using hexahedral cells, and sub-sequently exported in mesh format (.msh) into Open-FOAM 2, where it is filled with molecules of two species,Fluid I and Fluid II, as shown in Figure 3.The two fluids are essentially isotopes of argon; theLennard-Jones potential is used for all I-I, II-I and II-IIfluid interactions. Both fluids have identical properties2OpenFOAM is an open-source toolbox of C++ libraries forcomputational physics in which our MD code is implemented. Itmay be downloaded freely from [13].but different identification (id) numbers, so that mix-ing can be observed and measured. The non-periodicboundary conditions at the inlet A (Fig 3), supply moleculesof Fluid I at a constant rate, u∗A = (0, 0, 0.25), while in-lets B1, B2 supply molecules of Fluid II at the same rate:u∗B1 = (0.07, 0, 0), u∗B2 = (−0.07, 0, 0). The NPBCs atthe inlets are of the same description as in the previouscase for Poiseuille flow, with target density and temper-ature set to ρ∗ = 0.6 and T ∗ = 2.4 respectively. Atthe outlet C, the mass flux boundary model removesmolecules of Fluid I and II equally, based on the totalzonal density ρ∗Z = 0.55. A Berendsen thermostat is ap-plied in the whole domain to keep the local temperatureclose to T ∗ = 2.4. The values for inlet velocities andthe zonal density were selected using a trial and errorprocedure until satisfactory mixing was observed at theoutlet of the channel (see Figure 4).Inlet AInlet B1Outlet C10nmxyzFluid IFluid IIInlet B2Figure 3: The nano-mixer MD case shown in its mixedstate, i.e., after time t∗ = 2000. Inlet NPBCs are ap-plied at A,B1,B2 boundaries, depicted by the shadedregions, while an outlet NPBC is applied at boundaryC. Note: only 30,000 random molecules are shown outof the ∼200,000 molecules that occupy the simulationdomain. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0  10  20  30  40  50  60DENSITY,  ρ∗OUTLET CHANNEL WIDTH, x*Fluid IFluid IIFluid I+IISAMPLING135σFigure 4: The partial-density cross-channel distributionfor both fluids, sampled at the outlet of the mixer.3.3 Nano-droplet deformation inextensional flow channelsWe demonstrate the application of our non-periodicboundary conditions to extensional-type flow geometriesin order to investigate emulsion droplet deformation andbreakup. Such studies are important in nanofluidic andmicrofluidic technologies such as heat-exchange devices,lab-on-a-chip components and ink-jet printing.The domain geometry consists of a cross-shaped nano-channel similar to the mixing channel above. Inlet NPBCsare applied at opposing boundaries (A1,A2), includingcontrol of density at ρ∗ = 0.6, temperature at T ∗ = 2.4,and velocity constrained to point inwards to the domain,i.e. u∗A1 = (0,−0.25, 0) and u∗A2 = (0, 0.25, 0). At theoutlets (B1,B2), two independent mass flux boundarymodels are applied at a value of density ρ∗Z = 0.57,where Z is their respective outlet channel zone. Theseboundary conditions are imposed in order to create anextensional shear flow at the throat of the channel whichacts to elongate, deform and breakup the three-dimensionalviscous droplets that will pass through the channel.(a) Throat Case A (b) Throat Case B (c) Throat Case CFigure 5: The three trial throat geometries used for theextensional droplet flow cases. Solvent molecules areexcluded from the snapshots for clarity.Droplets of 6nm diameter are inserted close to theinlet boundary A1 at a rate of 2.314× 109 s−1, so thatmultiple droplets do not interact during their motionthrough the channel. The interaction potentials are pur-posely chosen to prevent the molecules in the droplet (d)from diffusing into the solvent (s). A simple method toimplement this is to use the repulsive WCA potential forthe droplet-solvent pair-interactions with the followingmodified characteristics σ∗sd = 1.4 and ǫ∗sd = 0.8. Thedroplet-droplet and solvent-solvent interactions are cho-sen to be the common shifted LJ potential with standardcharacteristics σ∗dd = σ∗ss = 1.0 and ǫ∗dd = ǫ∗ss = 1.0.Deformation and breakup characteristics are gener-ally dependent on a number of factors, such as dropletsize, boundary conditions and geometry of the throat.Here we investigate droplet deformation caused by thegeometry of the throats and examine the three cases inFigure 5.The results show that breakup into two smaller dropletsoccurs more repeatedly in the throat case C, in whichthe throat geometry at the exit is narrower than the twoother cases. In this case, however, all droplets that exitthe throat tend to stick to the boundary of the channeldue to the attraction characteristics of the wall and lowspeeds which are present close to the wall. In the othercases, droplets sticking to the wall are rare, due to thelarger exit throat sizes.4 CONCLUSIONSWe have presented new algorithms for applying inletand outlet non-periodic boundary conditions (NPBCs)in molecular dynamics simulations, in order to set uppressure gradient flows in arbitrary domain geometries.In other work we show how the NPBCs may also beapplied to couple a continuum solver with MD, in a hy-brid simulation. The boundary conditions are testedin a simple Poiseuille flow channel, in a complex 3Dmixing-channel and in a cross-shaped extensional flownano-channel (to investigate nano-droplet deformationat the throat). The boundary method we propose is par-allelised and generally applicable to any complex nano-channel configuration, and so may serve as a useful toolin engineering design and simulation for nano-scale ap-plications.5 ACKNOWLEDGMENTSThis work is funded in the UK by the James WeirFoundation (MKB), the Institution of Mechanical En-gineers (WDN) and the Royal Society of Edinburgh /Scottish Government (JMR).REFERENCES[1] D. C. Rapaport. The Art of Molecular DynamicsSimulation. Cambridge University Press, 2nd edi-tion, 2004.[2] M. P. Allen and D. J. Tildesley. Computer Simula-tion of Liquids. Oxford University Press, 1987.[3] G. B. Macpherson and J. M. Reese. Molecular dy-namics in arbitrary geometries: Parallel evaluationof pair forces. Molecular Simulation, 34(1):97–115,Jan 2008.[4] G. B. Macpherson, M. K. Borg, and J. M. Reese.Generation of initial molecular dynamics configura-tions in arbitrary geometries and in parallel. Molec-ular Simulation, 33(15):1199–1212, December 2007.[5] E. M. Kotsalis, J. H. Walther, and P. Koumout-sakos. Control of density fluctuations in atomistic-continuum simulations of dense liquids. PhysicalReview E, 76:016709, 2007.[6] M. K. Borg, G. B. Macpherson, and J. M. Reese.Controllers for imposing continuum-to-molecularboundary conditions in arbitrary fluid flow geome-tries. Accepted by Molecular Simulation, 2010.[7] J. Koplik, J. R. Banavar, and J. F. Willemsen.Molecular dynamics of Poiseuille flow and movingcontact lines. Physical Review Letters, 60:1282–1285, March 1988.[8] R. Delgado-Buscalioni and P. V. Coveney.USHER: An algorithm for particle insertionin dense fluids. Journal of Chemical Physics,119(2):978–987, 2003.[9] T. Werder, J. H. Walther, and P. Koumoutsakos.Hybrid atomistic-continuum method for the simu-lation of dense fluid flows. Journal of Computa-tional Physics, 205(1):373–390, May 2005.[10] M. Sun and C. Ebner. Molecular-dynamics simu-lation of compressible fluid flow in two-dimensionalchannels. Physical Review A, 46(8):4813–4818, Oct1992.[11] G. Ciccotti and A. Tenenbaum. Canonical ensembleand nonequilibrium states by molecular dynamics.Journal of Statistical Physics, 23(6):767–772, Dec1980.[12] D. E. Hertzog, B. I., B. Mohammadi, O. Bakajin,and J. G. Santiago. Optimization of a microfluidicmixer for studying protein folding kinetics. Jour-nal of Analytical Chemistry, 78(13):4299–4306, Jul2006.[13] OpenFOAM: The open source CFD toolbox. Avail-able online: www.openfoam.org.

An algorithm for particle insertion in dense ﬂuids.

Canonical ensemble and nonequilibrium states by molecular dynamics.

Computer Simulation of Liquids.

Control of density ﬂuctuations in atomisticcontinuum simulations of dense liquids. Physical Review E,

Controllers for imposing continuum-to-molecular boundary conditions in arbitrary ﬂuid ﬂow geometries. Accepted by Molecular Simulation,

Generation of initial molecular dynamics configurations in arbitrary geometries and in parallel.

Generation of initial molecular dynamics conﬁgurations in arbitrary geometries and in parallel. Molecular Simulation,

Hybrid atomistic-continuum method for the simulation of dense ﬂuid ﬂows.

Molecular dynamics in arbitrary geometries: Parallel evaluation of pair forces.

Molecular dynamics of Poiseuille ﬂow and moving contact lines. Physical Review Letters,

Molecular-dynamics simulation of compressible ﬂuid ﬂow in two-dimensional channels.

OpenFOAM: The open source CFD toolbox.

Optimization of a microﬂuidic mixer for studying protein folding kinetics.

The Art of Molecular Dynamics Simulation.

http://strathprints.strath.ac.uk/32105/1/Reese_JM_Pure_MD_boundary_conditions_for_pressure_gradient_flows_Nano_mixing_and_nano_droplet_deformation_in_extensional_flows_2011.pdf

MD boundary conditions for pressure gradient flows : nano-mixing and nano-droplet deformation in extensional flows

Abstract

Similar works

Full text

Available Versions

Name not available

University of Strathclyde Institutional Repository