# Coupling Strategies Investigation of Hybrid Atomistic-Continuum Method Based on State Variable Coupling

Authors: Qian Wang; Xiao-Guang Ren; Xin-Hai Xu; Chao Li; Hong-Yu Ji; Xue-Jun Yang
Journal: Advances in Materials Science and Engineering
 (2017)

---

## Abstract

Different configurations of coupling strategies influence greatly the accuracy and convergence of the simulation results in the hybrid atomistic-continuum method. This study aims to quantitatively investigate this effect and offer the guidance on how to choose the proper configuration of coupling strategies in the hybrid atomistic-continuum method. We first propose a hybrid molecular dynamics- (MD-) continuum solver in LAMMPS and OpenFOAM that exchanges state variables between the atomistic region and the continuum region and evaluate different configurations of coupling strategies using the sudden start Couette flow, aiming to find the preferable configuration that delivers better accuracy and efficiency. The major findings are as follows:(1) the C→A region plays the most important role in the overlap region and the “4-layer-1” combination achieves the best precision with a fixed width of the overlap region; (2) the data exchanging operation only needs a few sampling points closer to the occasions of interactions and decreasing the coupling exchange operations can reduce the computational load with acceptable errors; (3) the nonperiodic boundary force model with a smoothing parameter of 0.1 and a finer parameter of 20 can not only achieve the minimum disturbance near the MD-continuum interface but also keep the simulation precision.


---

## Body

## 1. Introduction

In recent years, with the rapid development of nanotechnology, microscale/nanoscale devices such as microelectromechanical system (MEMS) devices and lab-on-a-chip devices have been widely used. The fluid flows in these devices involve a broad range of scales from atomistic scale to macroscopic scale [1]. Generally speaking, fluid simulation based on the continuum assumption uses Navier-Stokes (NS) equations to investigate fluid dynamics under macroscopic scale. As the characteristic scale decreases, the fluid flows at the microscales/nanoscales exhibit quite different properties with the flows at macroscale, such as the invalidity of the continuum assumption [2] and increased viscosity in nanochannels [3]. Molecular dynamics (MD), one of widely used microfluid simulation methods, resolves fluid features at microscales/nanoscales. However, its computation-intensive feature brings heavy loads to both simulation time and memory usage, limiting the simulation scale nanometer in length scale and nanosecond in time scale. In order to simulate physical problems with a large length scale and to capture microscopic physical phenomena, many multiscale simulation methods have been proposed. For solid simulation, the “bridging domain” and “bridging scale” method [4, 5] uses Lagrangian multipliers and solution projection method to seamlessly couple two solvers in different scales with few unusual effects. For dense liquid simulation, the hybrid atomistic-continuum method (HAC) has been raised [6–8]. HAC applies molecular dynamics in regions where the atomistic description is needed, for example, boundary regions and corner vortex regions, while using the continuum method in the remaining regions to obtain both computation efficiency and simulation accuracy. In this paper, we focus on the simulation of the dense fluid.There are two types of coupling approaches available: the flux-based method [9, 10] and the state variable-based method [6]. The former uses mass flux, momentum flux, and energy flux to exchange data between the continuum method and molecular dynamics, thus satisfying conservation laws naturally. The latter couples these two methods using mass and momentum state variables to perform simulation. For incompressible problems, the computational cost involved in calculating the molecular fluxes across a surface is much higher than the calculation of state variables only [11]. Therefore, we focus on the state variable-based HAC method in this paper.The first hybrid method combining molecular dynamics with the continuum method for dense fluid was proposed by O’Connell and Thompson [6]. For the one-dimensional Couette flow, the domain was split into an atomistic region and a continuum region, using an overlap region to alleviate dramatic density oscillation and couple the results of these two regions. The overlap region contains a nonperiodic boundary force region (npbf region), an atomistic-coupled-to-continuum region (A→C region) and a continuum-coupled-to-atomistic region (C→A region). However, this method has a limitation that it does not cope with the mass transfer across the MD-continuum interface.Later, HAC models emerged differing in forms of coupling strategies, boundary conditions extraction, and nonperiodic boundary force models [12–16]. In order to cope with the mass flux transfer, some researchers [17, 18] introduced themass flow region; other researchers [14, 19] brought forward thebuffer region to further relax the fluctuations between the MD results and the continuum results; still other researchers [20–22] proposed different expressions for the nonperiodic boundary force. More detailed review is given by Mohamed and Mohamad [23].Nevertheless, over the past 20 years, several issues still remained to be investigated in the hybrid atomistic-continuum method based on state variables. Strategy configurations of the spatial coupling, temporal coupling, and associated parameters have significant influences on simulation efficiency and accuracy. The existing research simply configured these strategies from their own point of view, while detailed analysis of these coupling strategies has never been performed. Firstly, for the HAC methods based on domain decomposition, spatial configurations of the overlap region in the existing approaches differ from each other and mostly set functional regions with the same width. Only Yen et al. [8] explored the appropriate size of the pure MD region and the overlap region. Secondly, on the occasions of data exchanging, existing approaches generally use the average of all sampling points of subsequent MD steps, to alleviate thermal noise from finite space and time sampling in the MD region. However, the increased quantity of samples brings better elimination of thermal noise and results in time lagging on the average results transferred to the continuum region. It is important to explore the effects of different quantities of sampling points and occasions for data exchanging on the convergence of the coupled simulation. Finally, when dealing with nonperiodic boundaries in the MD region, Issa and Poesio [13] proposed the FoReV algorithm and empirically configured a smoothing parameter. But when coupled to the continuum method, the effects of choosing different parameters on the local liquid structure near the MD-continuum interface have to be investigated.Based on the above analysis, we could conclude that there exist several coupling strategies worthy of further investigation. In this paper, we design a domain decomposition type of hybrid MD-continuum solver using open source software LAMMPS [24] and OpenFOAM [25] and investigate the coupling strategy issues of the HAC simulation usingCouette channel flow as the model flow. The main contributions of this paper are summarized as follows:(i)
We analyze the effects of choosing different spatial strategies for configuring the overlap region on the model flow. The finding is that, when given the equal width of the functional regions, “5-layer” combination obtains the best numerical precision as the width of the overlap region varies. By contrast, “4-layer-1” combination is proved to be the best settings while given the fixed width of the overlap region. Apart from that, we also find out that enlarging the C→A region could result in a better simulation accuracy in the model flow.(ii)
We investigate the more efficient temporal strategies for data exchanging. The efficiency is obtained by the study on the quantity of samples and the time points for data exchanging. The practical conclusion is that data exchanging in theA→C operation only needs a few sampling points which are close to the occasions of interactions to guarantee modeling efficiency and to reduce the times of sampling. With acceptable errors, we also find out that timely data exchanging performs better than the other settings, but decreasing the coupling exchange operations can further reduce the computational load.(iii)
We conduct analysis of parameters for the nonperiodic boundary force model on the model flow. We add a finer parameter to the force model in order to apply the FoReV algorithm in the HAC model effectively. Results indicate that under the domain decomposition along the flow direction, the proper combination of the smoothing parameter and the finer parameter can not only achieve the minimum disturbance on the local structure near the MD-continuum interface but also keep the simulation precision.The remainder of the paper is organized as follows: the hybrid atomistic-continuum simulation methodology is presented in Section2 and the discussion of mutable parameters for coupling strategies is described in Section 3. In Section 4, we apply these coupling parameters to the benchmark problems and compare the convergence and accuracy of the results of numerical tests. Section 5 concludes the paper.


## 2. Model Configuration and Methodology

In this section, we propose the MD-continuum solver based on the HAC physical model and coupling strategies using LAMMPS and OpenFOAM. OpenFOAM servers as the main framework which is highly modular and elegant extendibility [26] and LAMMPS is built as a library to be called from. Section 2.1 introduces the decomposition of the simulation domain; Section 2.2 includes the numerical methods of the atomistic region and the continuum region; Section 2.3 gives the configuration of the overlap region and associated coupling operations and Section 2.4 introduces the temporal coupling.


### 2.1. Domain Decomposition

In this paper, we investigate coupling strategy issues contained in the HAC simulation using Couette flows as the model flows which are under incompressible constant temperature condition. We take two kinds of boundary conditions into account, that is, the no-slip and the slip boundary condition. The computational domain and the coordinate system in the current HAC simulation are shown in Figures1 and 2. Under the no-slip boundary condition, the simulation domain is split into an atomistic region, a continuum region, and an overlap region (left), while there are two atomistic regions and two overlap regions under the slip boundary condition (right). The outer boundary of the continuum region which resides in the atomistic region is called hybrid solution interface (HSI). The atomistic regions include liquid fluid atom regions and wall atom regions and are located near the wall regions in order to provide accurate boundary conditions to the continuum part. In the current study, a two-dimensional simulation is performed in the continuum region, that is, only in xy plane, while a three-dimensional simulation is performed in the atomistic regions with z axis as the extension direction.Figure 1
Simulation domain decomposition and coordinates of the no-slip boundary condition;H is the height of the channel.Figure 2
Simulation domain decomposition and coordinates of the slip boundary condition;H is the height of the channel.


### 2.2. Atomistic Region and Continuum Region

In this section, we introduce numerical methods used in the atomistic region and the continuum region. Section2.2.1 shows the molecular dynamics simulation method and physical parameters and Section 2.2.2 gives the continuum solution using the finite volume method (FVM).


#### 2.2.1. Atomistic Region

The truncated and shifted Lennard-Jones (LJ) potential is used to model the interactions between fluid atoms as well as wall atoms in the atomistic region. The potential is given by(1)ϕr=4ϵσr12-σr6-σrc12+σrc6,r≤rc,0,r>rc,where rij is the distance between atom i and j, rc is the cutoff radius, and ϵ and σ are the characteristic molecular energy and the molecular length scale, respectively. In this study, we choose liquid Argon as the model liquid flow with LJ parameters σ=0.34nm, ϵ=1.67×10-21J, and m=6.63×10-26kg, where m is the atom mass and a well-defined liquid phase of Argon with TkBϵ-1=1.1, ρσ3=0.81, and the dynamic viscosity μ=2.14ϵτσ-3, where kB is the Boltzmann constant and τ=(mσ2/ϵ)1/2 is the characteristic time of the Lennard-Jones potential. We choose a cutoff distance rc=2.2σ to save the computation time. The wall atoms are modeled by two (111) planes with the x direction of the lattice along [112¯] orientation. For the two-dimensional continuum simulation, there is no flow in z direction. Therefore, in the atomistic region, x and z directions are under periodic boundary conditions while y direction is not.The constant temperature in the simulation system is maintained by Langevin thermostat [28, 29] that couples the system to a thermal reservoir through the addition of Gaussian noise and frictional terms. The thermostat only applies to z direction which is perpendicular to the bulk flow direction and dose not influence the bulk flow velocity. The equation of motion for i atom is given by(2)my¨i=-∑j≠i∂ϕ∂yj-mΓy˙i+ηi,where the summation donates all pair-interaction forces on atom i and Γ is the damping ratio. ηi denotes a Gaussian distributed force with zero mean and variance of 2mkBTΓ. We choose the damping ratio of Langevin thermostat Γ=1.0τ-1 here. Initially, the liquid atoms are arranged at the lattice and each velocity component of atoms conforms to the Gaussian distribution. The Newtonian equations for atoms are integrated with the velocity Verlet algorithm with time step Δt=0.005τ.


#### 2.2.2. Continuum Region

In the continuum region, we model the system using the incompressible Navier-Stokes equation. The NS equation and the continuity equation are given by(3)∇·u=0,∂u∂t+∇·uu=-1ρ∇P+ν∇2u,where u is the bulk flow velocity, P is the pressure, and ν is the kinematic viscosity. The equations above are solved numerically using the finite volume method. Density in continuum region is the same as the atomistic region and periodic boundary condition is applied in x direction. We solve this two-dimensional NS equation using the PISO algorithm with the icoFoam solver in OpenFOAM [25].


### 2.3. Overlap Region

We categorize the existing configurations of the overlap region into five functional regions, that is, theA→C region, the C→A region, thebuffer region, the nonperiodic boundary force region (npbf region), and themass flow region. The latter two are often unified to thecontrol region. The schematic of the overlap region is shown in Figure 3.Figure 3
Detailed schematic diagram of the overlap region.


#### 2.3.1.A→C Region

TheA→C region is located on the outer boundary of the continuum region centered on HSI in Figure 3. In this region, the A→C operation transfers the spatial and temporal average of atom velocities to the continuum region as the new velocity boundary. This region can be further divided into bins in x and y direction to match the mesh cells in the continuum region. For the ith bin, the boundary condition for the continuum velocityu is given by the following average equation:(4)ui=1Ni∑kNivk,where vk is the velocity of the kth atom in the bin i, Ni is the number of atoms in the bin i, and the bracket represents the temporal average. The most widely used sampling average methods are SAM and CAM [30]. We use the CAM method in our simulation with the following equation, and the temporal average is performed over S sampling points:(5)uiCAM=∑j=1S∑i=1Nkvk∑j=1SNk.However, sampling over finite spatial and temporal scales incurs statistical fluctuations. The HSI is the portion of the continuum region that receives boundary data from the atomistic region. Due to inherent statistical fluctuations in this data, the boundary conditions in the complete continuum region boundary may not exactly mass conserve. The most common remedy [31] is to apply a correction factor as follows:(6)vHSI·ncorrected=vHSI·n-∫ϕvϕ·ndS∫HSIdS,where vHSI is the calculated velocity on the HSI using CAM sampling technique,n is the normal vector to the boundary, dS is an element of the boundary, and ϕ is the whole boundary of the continuum region.


#### 2.3.2.C→A Region

TheC→A region transfers the velocities of continuum mesh cells to the velocities of atoms located in them through the crude constraint Lagrangian dynamics (CCLD) [6] and provides boundary conditions for the atomistic region. Similar to the A→C region, the C→A region can be divided into bins in x and y direction to match the mesh cells in continuum region.The CCLD requires that the mean atomistic velocities in bini should be equal to the average continuum velocities in the bin i; that is(7)u¯i=1Ni∑pNivp,where Ni is the number of atoms in the bin i. In this paper, the volume of the C→A region bins is equal to that of cells in the continuum region. Through CCLD, we can get the new velocities of atoms located in the bin i; that is(8)x˙jα=vjα+ξαujα-1Ni∑kNivkα,where α is x, y, or z direction and ξ is the constraint strength. The constraint strength is 0.01 as used by O’Connell and Thompson [6] while 1 by Nie et al. [7, 8]. In this paper, we set ξ=1.0 in x and y direction while ξ=0 in z direction.


#### 2.3.3. Nonperiodic Boundary Force Region

In order to prevent atoms from drifting away from the atomistic region, keep the quantity of atoms unchanged, and remedy the nonperiodic boundary in the HAC simulation, the existing research had proposed many external force models [7, 20, 21, 27]. Figure 4 shows the atoms missing interaction force near the MD-continuum interface. In this paper, we model an external force for these atoms to alleviate density fluctuation due to the nonperiodic boundary condition and reflecting wall based on the FoReV algorithm proposed by Issa and Poesio [13].Figure 4
Missing interaction region due to the nonperiodic boundary at the MD-continuum interfaceymax.Each atom traveling throughymax will be reflected back into the atomistic region with the same displacement across the interface and reversed velocity. This guarantees a constant number of atoms in the atomistic region. We divide the region near the MD-continuum interface into several bins. In Figure 5, we summarize the total external force experienced by atoms residing in bin k and use a feedback force to apply on atoms in the missing interactions region.Figure 5
Detailed schematic of thenpbf region residing in thecontrol region.Firstly, we calculate the average external force experienced by bink; that is(9)Fbin,k=∑i∈k,j∉kFijNk,where Fij is the interaction force between atoms i and j, Nk is the number of atoms in bin k, and Fbin,k is the average external force experienced by atoms in bin k. Through a feedback mechanism, we apply a normal reversed force on each atom that resides in the cutoff distance near the MD-continuum interface, which makes Fybin,k=0. This reversed force -Fb,kn+1 is given as follows, which is constructed with the simple exponential smoothing, including a smoothing parameter α:(10)Fb,kn+1=αFbin,kn+1+1-αFb,kn.In each MD time step, we calculateFbin,k innpbf region bins, use (10) to construct the feedback force Fb,k, and apply it to associated atoms to remedy nonperiodic boundary effect.


#### 2.3.4. Mass Flow Region

For the incompressible condition, the total number of atoms in the atomistic region should be kept unchanged. In order to simulate themass flux across the MD-continuum interface, we bring forward themass flow region and use it to deal with the atom insertion and deletion. Themass flow region can be split into several bins to meet with continuum cells. In one continuum time step ΔtCFD, the number of atoms to be inserted into or deleted from the mass flow bin i is given by(11)n=AiρuiαΔtCFDm,where Ai=(Δx×Δz)·n is the bin area normal to y direction andn is the face normal vector, pointed into the atomistic region. uiα is the continuum velocity, ρ is the liquid density, and m is the atom mass. If n is negative, certain atoms should be deleted from themass flow region and if n is positive, certain atoms should be inserted into themass flow region. Since an atom is inseparable, the nearest integer is taken and the remaining fraction is included at the next insert/delete operation. In Section 2.3.3, the reflecting wall boundary condition prevents atoms from drifting away freely, so the amount of atoms in the atomistic region can only be changed through the insert and delete operation.The atom insertion is via the USHER algorithm [32]. USHER algorithm tries to find a point in a certain bin with the potential that is equal to the average potential of the bin. The initial insertion point is randomly chosen and updated by Newton-Raphson method. The atom deletion is to remove those atoms near the MD-continuum interface ymax. The main thought of the algorithm of deletion is to choose atoms that are most probable to leave the atomistic region using the distance to ymax and the velocity normal to ymax.


#### 2.3.5. Buffer Region

Thebuffer region is located between the A→C region and the C→A region or between thecontrol region and the C→A region. We use it to alleviate the fluctuation of the artificial operations on atoms. In thebuffer region, there is no atom insertion or deletion and no constraint on atom velocities. But we have to narrow the width of thebuffer region to decrease the computational load but keep a certain width to relax the results of the atomistic region and the continuum region before being coupled together. We will discuss how to set thebuffer region in later sections.


### 2.4. Temporal Coupling

There are three time variables to be considered in the HAC simulation [33]: integration time steps of Newtonian equation ΔtMD, integration time steps of Navier-Stokes equation ΔtCFD, and sampling average time Δtave. As we have mentioned before, ΔtMD is chosen as 0.005τ in our paper. In order to solve NS equation accurately, firstly, ΔtCFD must be far less than characteristic time on the mesh, that is, diffusion time ρΔxΔy/μ, where μ is the dynamic viscosity and μ=2.14ϵτσ-3. Secondly, it must be larger than the decay time of velocity autocorrelation function tvv=0.14ρ-1T-1/2 so as to reduce thermal noise from HSI [6]. Finally, it should meet the CFL condition [34], uflowΔtCFD<Δx/2. The CFD time step we choose in this paper is ΔtCFD=M×ΔtMD=0.5τ, where M is 100. The time advancing mechanism in the HAC simulation is shown in Figure 6, which is sequential coupling.Figure 6
Time coupling and advancing mechanism in the current HAC method. In one coupling cycle, it includes four steps, that is, the CFD advancing, theC→A operation, the MD advancing, and the A→C operation.In each coupling simulation time step, CFD advances a certain time, for example, oneΔtCFD, and transfers continuum velocities to the atomistic region through CCLD in the C→A operation. Then, MD advances the same time interval, that is, M×ΔtMD, samples and averages atom velocities, and then passes them to the continuum region in the A→C operation, thus finishing one cycle of HAC simulation.Indeed, due to the small time stepΔtMD, the HAC simulation time of a benchmark case will take a long time. However, the efficiency of the HAC method is defined by comparing the full MD simulation for the same scale of the benchmark case. The HAC method only applies atomistic simulation in part of the simulation domain. Obviously the HAC method is much more efficient than the full MD method with simulation of the full domain.


## 2.1. Domain Decomposition

In this paper, we investigate coupling strategy issues contained in the HAC simulation using Couette flows as the model flows which are under incompressible constant temperature condition. We take two kinds of boundary conditions into account, that is, the no-slip and the slip boundary condition. The computational domain and the coordinate system in the current HAC simulation are shown in Figures1 and 2. Under the no-slip boundary condition, the simulation domain is split into an atomistic region, a continuum region, and an overlap region (left), while there are two atomistic regions and two overlap regions under the slip boundary condition (right). The outer boundary of the continuum region which resides in the atomistic region is called hybrid solution interface (HSI). The atomistic regions include liquid fluid atom regions and wall atom regions and are located near the wall regions in order to provide accurate boundary conditions to the continuum part. In the current study, a two-dimensional simulation is performed in the continuum region, that is, only in xy plane, while a three-dimensional simulation is performed in the atomistic regions with z axis as the extension direction.Figure 1
Simulation domain decomposition and coordinates of the no-slip boundary condition;H is the height of the channel.Figure 2
Simulation domain decomposition and coordinates of the slip boundary condition;H is the height of the channel.


## 2.2. Atomistic Region and Continuum Region

In this section, we introduce numerical methods used in the atomistic region and the continuum region. Section2.2.1 shows the molecular dynamics simulation method and physical parameters and Section 2.2.2 gives the continuum solution using the finite volume method (FVM).


### 2.2.1. Atomistic Region

The truncated and shifted Lennard-Jones (LJ) potential is used to model the interactions between fluid atoms as well as wall atoms in the atomistic region. The potential is given by(1)ϕr=4ϵσr12-σr6-σrc12+σrc6,r≤rc,0,r>rc,where rij is the distance between atom i and j, rc is the cutoff radius, and ϵ and σ are the characteristic molecular energy and the molecular length scale, respectively. In this study, we choose liquid Argon as the model liquid flow with LJ parameters σ=0.34nm, ϵ=1.67×10-21J, and m=6.63×10-26kg, where m is the atom mass and a well-defined liquid phase of Argon with TkBϵ-1=1.1, ρσ3=0.81, and the dynamic viscosity μ=2.14ϵτσ-3, where kB is the Boltzmann constant and τ=(mσ2/ϵ)1/2 is the characteristic time of the Lennard-Jones potential. We choose a cutoff distance rc=2.2σ to save the computation time. The wall atoms are modeled by two (111) planes with the x direction of the lattice along [112¯] orientation. For the two-dimensional continuum simulation, there is no flow in z direction. Therefore, in the atomistic region, x and z directions are under periodic boundary conditions while y direction is not.The constant temperature in the simulation system is maintained by Langevin thermostat [28, 29] that couples the system to a thermal reservoir through the addition of Gaussian noise and frictional terms. The thermostat only applies to z direction which is perpendicular to the bulk flow direction and dose not influence the bulk flow velocity. The equation of motion for i atom is given by(2)my¨i=-∑j≠i∂ϕ∂yj-mΓy˙i+ηi,where the summation donates all pair-interaction forces on atom i and Γ is the damping ratio. ηi denotes a Gaussian distributed force with zero mean and variance of 2mkBTΓ. We choose the damping ratio of Langevin thermostat Γ=1.0τ-1 here. Initially, the liquid atoms are arranged at the lattice and each velocity component of atoms conforms to the Gaussian distribution. The Newtonian equations for atoms are integrated with the velocity Verlet algorithm with time step Δt=0.005τ.


### 2.2.2. Continuum Region

In the continuum region, we model the system using the incompressible Navier-Stokes equation. The NS equation and the continuity equation are given by(3)∇·u=0,∂u∂t+∇·uu=-1ρ∇P+ν∇2u,where u is the bulk flow velocity, P is the pressure, and ν is the kinematic viscosity. The equations above are solved numerically using the finite volume method. Density in continuum region is the same as the atomistic region and periodic boundary condition is applied in x direction. We solve this two-dimensional NS equation using the PISO algorithm with the icoFoam solver in OpenFOAM [25].


## 2.2.1. Atomistic Region

The truncated and shifted Lennard-Jones (LJ) potential is used to model the interactions between fluid atoms as well as wall atoms in the atomistic region. The potential is given by(1)ϕr=4ϵσr12-σr6-σrc12+σrc6,r≤rc,0,r>rc,where rij is the distance between atom i and j, rc is the cutoff radius, and ϵ and σ are the characteristic molecular energy and the molecular length scale, respectively. In this study, we choose liquid Argon as the model liquid flow with LJ parameters σ=0.34nm, ϵ=1.67×10-21J, and m=6.63×10-26kg, where m is the atom mass and a well-defined liquid phase of Argon with TkBϵ-1=1.1, ρσ3=0.81, and the dynamic viscosity μ=2.14ϵτσ-3, where kB is the Boltzmann constant and τ=(mσ2/ϵ)1/2 is the characteristic time of the Lennard-Jones potential. We choose a cutoff distance rc=2.2σ to save the computation time. The wall atoms are modeled by two (111) planes with the x direction of the lattice along [112¯] orientation. For the two-dimensional continuum simulation, there is no flow in z direction. Therefore, in the atomistic region, x and z directions are under periodic boundary conditions while y direction is not.The constant temperature in the simulation system is maintained by Langevin thermostat [28, 29] that couples the system to a thermal reservoir through the addition of Gaussian noise and frictional terms. The thermostat only applies to z direction which is perpendicular to the bulk flow direction and dose not influence the bulk flow velocity. The equation of motion for i atom is given by(2)my¨i=-∑j≠i∂ϕ∂yj-mΓy˙i+ηi,where the summation donates all pair-interaction forces on atom i and Γ is the damping ratio. ηi denotes a Gaussian distributed force with zero mean and variance of 2mkBTΓ. We choose the damping ratio of Langevin thermostat Γ=1.0τ-1 here. Initially, the liquid atoms are arranged at the lattice and each velocity component of atoms conforms to the Gaussian distribution. The Newtonian equations for atoms are integrated with the velocity Verlet algorithm with time step Δt=0.005τ.


## 2.2.2. Continuum Region

In the continuum region, we model the system using the incompressible Navier-Stokes equation. The NS equation and the continuity equation are given by(3)∇·u=0,∂u∂t+∇·uu=-1ρ∇P+ν∇2u,where u is the bulk flow velocity, P is the pressure, and ν is the kinematic viscosity. The equations above are solved numerically using the finite volume method. Density in continuum region is the same as the atomistic region and periodic boundary condition is applied in x direction. We solve this two-dimensional NS equation using the PISO algorithm with the icoFoam solver in OpenFOAM [25].


## 2.3. Overlap Region

We categorize the existing configurations of the overlap region into five functional regions, that is, theA→C region, the C→A region, thebuffer region, the nonperiodic boundary force region (npbf region), and themass flow region. The latter two are often unified to thecontrol region. The schematic of the overlap region is shown in Figure 3.Figure 3
Detailed schematic diagram of the overlap region.


### 2.3.1.A→C Region

TheA→C region is located on the outer boundary of the continuum region centered on HSI in Figure 3. In this region, the A→C operation transfers the spatial and temporal average of atom velocities to the continuum region as the new velocity boundary. This region can be further divided into bins in x and y direction to match the mesh cells in the continuum region. For the ith bin, the boundary condition for the continuum velocityu is given by the following average equation:(4)ui=1Ni∑kNivk,where vk is the velocity of the kth atom in the bin i, Ni is the number of atoms in the bin i, and the bracket represents the temporal average. The most widely used sampling average methods are SAM and CAM [30]. We use the CAM method in our simulation with the following equation, and the temporal average is performed over S sampling points:(5)uiCAM=∑j=1S∑i=1Nkvk∑j=1SNk.However, sampling over finite spatial and temporal scales incurs statistical fluctuations. The HSI is the portion of the continuum region that receives boundary data from the atomistic region. Due to inherent statistical fluctuations in this data, the boundary conditions in the complete continuum region boundary may not exactly mass conserve. The most common remedy [31] is to apply a correction factor as follows:(6)vHSI·ncorrected=vHSI·n-∫ϕvϕ·ndS∫HSIdS,where vHSI is the calculated velocity on the HSI using CAM sampling technique,n is the normal vector to the boundary, dS is an element of the boundary, and ϕ is the whole boundary of the continuum region.


### 2.3.2.C→A Region

TheC→A region transfers the velocities of continuum mesh cells to the velocities of atoms located in them through the crude constraint Lagrangian dynamics (CCLD) [6] and provides boundary conditions for the atomistic region. Similar to the A→C region, the C→A region can be divided into bins in x and y direction to match the mesh cells in continuum region.The CCLD requires that the mean atomistic velocities in bini should be equal to the average continuum velocities in the bin i; that is(7)u¯i=1Ni∑pNivp,where Ni is the number of atoms in the bin i. In this paper, the volume of the C→A region bins is equal to that of cells in the continuum region. Through CCLD, we can get the new velocities of atoms located in the bin i; that is(8)x˙jα=vjα+ξαujα-1Ni∑kNivkα,where α is x, y, or z direction and ξ is the constraint strength. The constraint strength is 0.01 as used by O’Connell and Thompson [6] while 1 by Nie et al. [7, 8]. In this paper, we set ξ=1.0 in x and y direction while ξ=0 in z direction.


### 2.3.3. Nonperiodic Boundary Force Region

In order to prevent atoms from drifting away from the atomistic region, keep the quantity of atoms unchanged, and remedy the nonperiodic boundary in the HAC simulation, the existing research had proposed many external force models [7, 20, 21, 27]. Figure 4 shows the atoms missing interaction force near the MD-continuum interface. In this paper, we model an external force for these atoms to alleviate density fluctuation due to the nonperiodic boundary condition and reflecting wall based on the FoReV algorithm proposed by Issa and Poesio [13].Figure 4
Missing interaction region due to the nonperiodic boundary at the MD-continuum interfaceymax.Each atom traveling throughymax will be reflected back into the atomistic region with the same displacement across the interface and reversed velocity. This guarantees a constant number of atoms in the atomistic region. We divide the region near the MD-continuum interface into several bins. In Figure 5, we summarize the total external force experienced by atoms residing in bin k and use a feedback force to apply on atoms in the missing interactions region.Figure 5
Detailed schematic of thenpbf region residing in thecontrol region.Firstly, we calculate the average external force experienced by bink; that is(9)Fbin,k=∑i∈k,j∉kFijNk,where Fij is the interaction force between atoms i and j, Nk is the number of atoms in bin k, and Fbin,k is the average external force experienced by atoms in bin k. Through a feedback mechanism, we apply a normal reversed force on each atom that resides in the cutoff distance near the MD-continuum interface, which makes Fybin,k=0. This reversed force -Fb,kn+1 is given as follows, which is constructed with the simple exponential smoothing, including a smoothing parameter α:(10)Fb,kn+1=αFbin,kn+1+1-αFb,kn.In each MD time step, we calculateFbin,k innpbf region bins, use (10) to construct the feedback force Fb,k, and apply it to associated atoms to remedy nonperiodic boundary effect.


### 2.3.4. Mass Flow Region

For the incompressible condition, the total number of atoms in the atomistic region should be kept unchanged. In order to simulate themass flux across the MD-continuum interface, we bring forward themass flow region and use it to deal with the atom insertion and deletion. Themass flow region can be split into several bins to meet with continuum cells. In one continuum time step ΔtCFD, the number of atoms to be inserted into or deleted from the mass flow bin i is given by(11)n=AiρuiαΔtCFDm,where Ai=(Δx×Δz)·n is the bin area normal to y direction andn is the face normal vector, pointed into the atomistic region. uiα is the continuum velocity, ρ is the liquid density, and m is the atom mass. If n is negative, certain atoms should be deleted from themass flow region and if n is positive, certain atoms should be inserted into themass flow region. Since an atom is inseparable, the nearest integer is taken and the remaining fraction is included at the next insert/delete operation. In Section 2.3.3, the reflecting wall boundary condition prevents atoms from drifting away freely, so the amount of atoms in the atomistic region can only be changed through the insert and delete operation.The atom insertion is via the USHER algorithm [32]. USHER algorithm tries to find a point in a certain bin with the potential that is equal to the average potential of the bin. The initial insertion point is randomly chosen and updated by Newton-Raphson method. The atom deletion is to remove those atoms near the MD-continuum interface ymax. The main thought of the algorithm of deletion is to choose atoms that are most probable to leave the atomistic region using the distance to ymax and the velocity normal to ymax.


### 2.3.5. Buffer Region

Thebuffer region is located between the A→C region and the C→A region or between thecontrol region and the C→A region. We use it to alleviate the fluctuation of the artificial operations on atoms. In thebuffer region, there is no atom insertion or deletion and no constraint on atom velocities. But we have to narrow the width of thebuffer region to decrease the computational load but keep a certain width to relax the results of the atomistic region and the continuum region before being coupled together. We will discuss how to set thebuffer region in later sections.


## 2.3.1.A→C Region

TheA→C region is located on the outer boundary of the continuum region centered on HSI in Figure 3. In this region, the A→C operation transfers the spatial and temporal average of atom velocities to the continuum region as the new velocity boundary. This region can be further divided into bins in x and y direction to match the mesh cells in the continuum region. For the ith bin, the boundary condition for the continuum velocityu is given by the following average equation:(4)ui=1Ni∑kNivk,where vk is the velocity of the kth atom in the bin i, Ni is the number of atoms in the bin i, and the bracket represents the temporal average. The most widely used sampling average methods are SAM and CAM [30]. We use the CAM method in our simulation with the following equation, and the temporal average is performed over S sampling points:(5)uiCAM=∑j=1S∑i=1Nkvk∑j=1SNk.However, sampling over finite spatial and temporal scales incurs statistical fluctuations. The HSI is the portion of the continuum region that receives boundary data from the atomistic region. Due to inherent statistical fluctuations in this data, the boundary conditions in the complete continuum region boundary may not exactly mass conserve. The most common remedy [31] is to apply a correction factor as follows:(6)vHSI·ncorrected=vHSI·n-∫ϕvϕ·ndS∫HSIdS,where vHSI is the calculated velocity on the HSI using CAM sampling technique,n is the normal vector to the boundary, dS is an element of the boundary, and ϕ is the whole boundary of the continuum region.


## 2.3.2.C→A Region

TheC→A region transfers the velocities of continuum mesh cells to the velocities of atoms located in them through the crude constraint Lagrangian dynamics (CCLD) [6] and provides boundary conditions for the atomistic region. Similar to the A→C region, the C→A region can be divided into bins in x and y direction to match the mesh cells in continuum region.The CCLD requires that the mean atomistic velocities in bini should be equal to the average continuum velocities in the bin i; that is(7)u¯i=1Ni∑pNivp,where Ni is the number of atoms in the bin i. In this paper, the volume of the C→A region bins is equal to that of cells in the continuum region. Through CCLD, we can get the new velocities of atoms located in the bin i; that is(8)x˙jα=vjα+ξαujα-1Ni∑kNivkα,where α is x, y, or z direction and ξ is the constraint strength. The constraint strength is 0.01 as used by O’Connell and Thompson [6] while 1 by Nie et al. [7, 8]. In this paper, we set ξ=1.0 in x and y direction while ξ=0 in z direction.


## 2.3.3. Nonperiodic Boundary Force Region

In order to prevent atoms from drifting away from the atomistic region, keep the quantity of atoms unchanged, and remedy the nonperiodic boundary in the HAC simulation, the existing research had proposed many external force models [7, 20, 21, 27]. Figure 4 shows the atoms missing interaction force near the MD-continuum interface. In this paper, we model an external force for these atoms to alleviate density fluctuation due to the nonperiodic boundary condition and reflecting wall based on the FoReV algorithm proposed by Issa and Poesio [13].Figure 4
Missing interaction region due to the nonperiodic boundary at the MD-continuum interfaceymax.Each atom traveling throughymax will be reflected back into the atomistic region with the same displacement across the interface and reversed velocity. This guarantees a constant number of atoms in the atomistic region. We divide the region near the MD-continuum interface into several bins. In Figure 5, we summarize the total external force experienced by atoms residing in bin k and use a feedback force to apply on atoms in the missing interactions region.Figure 5
Detailed schematic of thenpbf region residing in thecontrol region.Firstly, we calculate the average external force experienced by bink; that is(9)Fbin,k=∑i∈k,j∉kFijNk,where Fij is the interaction force between atoms i and j, Nk is the number of atoms in bin k, and Fbin,k is the average external force experienced by atoms in bin k. Through a feedback mechanism, we apply a normal reversed force on each atom that resides in the cutoff distance near the MD-continuum interface, which makes Fybin,k=0. This reversed force -Fb,kn+1 is given as follows, which is constructed with the simple exponential smoothing, including a smoothing parameter α:(10)Fb,kn+1=αFbin,kn+1+1-αFb,kn.In each MD time step, we calculateFbin,k innpbf region bins, use (10) to construct the feedback force Fb,k, and apply it to associated atoms to remedy nonperiodic boundary effect.


## 2.3.4. Mass Flow Region

For the incompressible condition, the total number of atoms in the atomistic region should be kept unchanged. In order to simulate themass flux across the MD-continuum interface, we bring forward themass flow region and use it to deal with the atom insertion and deletion. Themass flow region can be split into several bins to meet with continuum cells. In one continuum time step ΔtCFD, the number of atoms to be inserted into or deleted from the mass flow bin i is given by(11)n=AiρuiαΔtCFDm,where Ai=(Δx×Δz)·n is the bin area normal to y direction andn is the face normal vector, pointed into the atomistic region. uiα is the continuum velocity, ρ is the liquid density, and m is the atom mass. If n is negative, certain atoms should be deleted from themass flow region and if n is positive, certain atoms should be inserted into themass flow region. Since an atom is inseparable, the nearest integer is taken and the remaining fraction is included at the next insert/delete operation. In Section 2.3.3, the reflecting wall boundary condition prevents atoms from drifting away freely, so the amount of atoms in the atomistic region can only be changed through the insert and delete operation.The atom insertion is via the USHER algorithm [32]. USHER algorithm tries to find a point in a certain bin with the potential that is equal to the average potential of the bin. The initial insertion point is randomly chosen and updated by Newton-Raphson method. The atom deletion is to remove those atoms near the MD-continuum interface ymax. The main thought of the algorithm of deletion is to choose atoms that are most probable to leave the atomistic region using the distance to ymax and the velocity normal to ymax.


## 2.3.5. Buffer Region

Thebuffer region is located between the A→C region and the C→A region or between thecontrol region and the C→A region. We use it to alleviate the fluctuation of the artificial operations on atoms. In thebuffer region, there is no atom insertion or deletion and no constraint on atom velocities. But we have to narrow the width of thebuffer region to decrease the computational load but keep a certain width to relax the results of the atomistic region and the continuum region before being coupled together. We will discuss how to set thebuffer region in later sections.


## 2.4. Temporal Coupling

There are three time variables to be considered in the HAC simulation [33]: integration time steps of Newtonian equation ΔtMD, integration time steps of Navier-Stokes equation ΔtCFD, and sampling average time Δtave. As we have mentioned before, ΔtMD is chosen as 0.005τ in our paper. In order to solve NS equation accurately, firstly, ΔtCFD must be far less than characteristic time on the mesh, that is, diffusion time ρΔxΔy/μ, where μ is the dynamic viscosity and μ=2.14ϵτσ-3. Secondly, it must be larger than the decay time of velocity autocorrelation function tvv=0.14ρ-1T-1/2 so as to reduce thermal noise from HSI [6]. Finally, it should meet the CFL condition [34], uflowΔtCFD<Δx/2. The CFD time step we choose in this paper is ΔtCFD=M×ΔtMD=0.5τ, where M is 100. The time advancing mechanism in the HAC simulation is shown in Figure 6, which is sequential coupling.Figure 6
Time coupling and advancing mechanism in the current HAC method. In one coupling cycle, it includes four steps, that is, the CFD advancing, theC→A operation, the MD advancing, and the A→C operation.In each coupling simulation time step, CFD advances a certain time, for example, oneΔtCFD, and transfers continuum velocities to the atomistic region through CCLD in the C→A operation. Then, MD advances the same time interval, that is, M×ΔtMD, samples and averages atom velocities, and then passes them to the continuum region in the A→C operation, thus finishing one cycle of HAC simulation.Indeed, due to the small time stepΔtMD, the HAC simulation time of a benchmark case will take a long time. However, the efficiency of the HAC method is defined by comparing the full MD simulation for the same scale of the benchmark case. The HAC method only applies atomistic simulation in part of the simulation domain. Obviously the HAC method is much more efficient than the full MD method with simulation of the full domain.


## 3. Mutable Parameters in Coupling Strategies

Concentrating on the influences of configuring different coupling strategies on accuracy and efficiency of the HAC simulation, in this section, the coupling strategies are embodied into the parameters listed as follows: configurations of the functional regions in Section3.1, variables of data exchanging in Section 3.2, and coefficients of the nonperiodic boundary force model in Section 3.3.


### 3.1. Parameters of Functional Regions:Layer Number andLayer Width

In the HAC simulation, the physical results of the atomistic region and the continuum region should keep consistent in the overlap region. The legal configuration of the overlap region can exchange data between continuum solver and atomistic solver correctly and alleviate the unusual effects due to artificial operations. The width of each component of the overlap region will lead to different simulation accuracy. Therefore, the configuration of the functional regions must be carefully designed and tested.Generally speaking, the overlap region should be located at a certain distance away from the solid wall. In Section2.3, we mention five functional regions to carry out data coupling. In the existing research, there are different combinations of functional regions and different widths of them. We summarize thelayer number (LN) and thelayer width (LW) in Table 1.Table 1
Configurations of functional regions in the overlap region.
Research LN Configuration details (region) LW O’Connell and Thompson 1995 [6] 3 npbf, A→C, C→A Equal Nie et al. 2004 [7] 3 control (npbf, mass flow), A→C, C→A Equal Zhou et al. 2014 [27] 3 control (npbf, mass flow, C→A), buffer, A→C Equal Cosden and Lukes 2013 [18] 4 mass flow, C→A, buffer, A→C Equal Sun et al. 2010 [19] 4 depress (extending to the continuum region, npbf), C→A, buffer, A→C Equal Yen et al. 2007 [8] 5 depress (extending to the continuum region, npbf), buffer, C→A, buffer, A→C Equal Ko et al. 2010 [14] 5 control (npbf, mass flow), buffer, C→A, buffer, A→C EqualFor the completeness of comparison, we add an alternative configuration which includes four layers, that is, thecontrol region (thenpbf region and themass flow region), thebuffer region, the C→A region, and the A→C region. Those configurations that contain thedepress region extending into the continuum region in Table 1 are not considered in the current comparison. We can also see that the functional regions have equal width in the previous research.In this paper, themass flow region is located at the top of the atomistic region, which is near the MD-continuum interface, and its width is equal to one CFD mesh cell width. Thenpbf region is located within a cutoff distance near the MD-continuum interface. The widths of these two remain unchanged.For theC→A region, the larger its width, the stronger the control exerted on the atomistic region by the continuum region and the more obvious effect CCLD has on the data exchanging in this region. For thebuffer region, its function is to alleviate the fluctuations among other kernel functional regions. Increasing its width means increasing the computation load, so it is better to minimize the width of thebuffer region or to remove it. For the A→C region, due to the existence of statistical thermal fluctuations in the atomistic region, the scope of it apparently influences the sampling average results and thus influences its capability to provide accurate boundary conditions to the continuum region.On the premise of accurate HAC simulation, different combinations of functional regions and different widths of them may have influences on simulation results. In Section4.2, we will test and discuss such influences and give suggestions about how to configure the overlap region properly.


### 3.2. Parameters of the Data Exchanging Operation:N and ex_steps

In Section3.1, we have talked about spatial average in the A→C region, while in this section we consider the temporal average parameters of the A→C operation. Generally speaking, in order to correctly extract macroscopic velocity, averages are taken over microscopic variables within a control volume and during time, which is called the binning method [35]. Specifically, in our HAC model, the control volume is the bin i in the A→C region and the time interval for averaging is Δtave for the continuum boundary extraction. Timely and correctly exchanging sampling atomistic data to the continuum solver is a kernel operation in the HAC method. We must pay attention that the temporal average over all sample points will introduce time lagging boundary results and would not exchange coupling information timely. Due to the overload of sampling operations, we should also question whether decreasing the times of coupling data exchanging operations could further reduce the sampling times and the computational load within acceptable errors.There are several parameters involved in the temporal average operation:(i)
ex_steps: the number of CFD time steps corresponding to one HAC data exchange operation Δtave, and Δtave=ex_steps×ΔtCFD(ii)
N: the number of sampling MD points in one data exchange operation(iii)
interval: the interval of sampling points in one data exchange operationIn the previous research [7, 8], data was often exchanged between the continuum region and the atomistic region in each CFD time step, and temporal average is performed over all M time steps in one ΔtCFD; that is, ex_steps=1, N=M, and interval=1.In Section4.3, we will discuss and test the proper sampling number N and exchange times ex_steps in one data exchange operation. Furthermore, we will provide guidance on the macroscopic velocity extraction.


### 3.3. Parameters of the Nonperiodic Boundary Force Model:α and β

Due to the lack of interaction with missing atoms, the atoms near the nonperiodic boundary will show unphysical effects. A nonperiodic boundary is applied to render density inhomogeneous and force nonuniform near the MD-continuum interface. In the existing HAC models, the nonperiodic boundary is often remedied with an external force model on those atoms near the boundary. This remedy operation will keep the pressure in the atomistic region correctly and alleviate the density fluctuation.We construct thenpbf region near the MD-continuum interface within a cutoff radius distance. Then, we refine thenpbf region intonpbf bins with width of binynpbf on the basis of CFD mesh width Lycell in y direction; that is, Lycell=β×binynpbf, where β is the finer parameter. Issa and Poesio [13] only provide a nonperiodic boundary force model with an empirical smoothing parameter α=1×10-4. A proper combination of the smoothing parameter α and the finer parameter β to minimize disturbance on the local liquid structure near the MD-continuum interface is worth investigating. In Section 4.4, we will provide numerical verifications on different combinations of α and β and discuss the appropriate values of them.


## 3.1. Parameters of Functional Regions:Layer Number andLayer Width

In the HAC simulation, the physical results of the atomistic region and the continuum region should keep consistent in the overlap region. The legal configuration of the overlap region can exchange data between continuum solver and atomistic solver correctly and alleviate the unusual effects due to artificial operations. The width of each component of the overlap region will lead to different simulation accuracy. Therefore, the configuration of the functional regions must be carefully designed and tested.Generally speaking, the overlap region should be located at a certain distance away from the solid wall. In Section2.3, we mention five functional regions to carry out data coupling. In the existing research, there are different combinations of functional regions and different widths of them. We summarize thelayer number (LN) and thelayer width (LW) in Table 1.Table 1
Configurations of functional regions in the overlap region.
Research LN Configuration details (region) LW O’Connell and Thompson 1995 [6] 3 npbf, A→C, C→A Equal Nie et al. 2004 [7] 3 control (npbf, mass flow), A→C, C→A Equal Zhou et al. 2014 [27] 3 control (npbf, mass flow, C→A), buffer, A→C Equal Cosden and Lukes 2013 [18] 4 mass flow, C→A, buffer, A→C Equal Sun et al. 2010 [19] 4 depress (extending to the continuum region, npbf), C→A, buffer, A→C Equal Yen et al. 2007 [8] 5 depress (extending to the continuum region, npbf), buffer, C→A, buffer, A→C Equal Ko et al. 2010 [14] 5 control (npbf, mass flow), buffer, C→A, buffer, A→C EqualFor the completeness of comparison, we add an alternative configuration which includes four layers, that is, thecontrol region (thenpbf region and themass flow region), thebuffer region, the C→A region, and the A→C region. Those configurations that contain thedepress region extending into the continuum region in Table 1 are not considered in the current comparison. We can also see that the functional regions have equal width in the previous research.In this paper, themass flow region is located at the top of the atomistic region, which is near the MD-continuum interface, and its width is equal to one CFD mesh cell width. Thenpbf region is located within a cutoff distance near the MD-continuum interface. The widths of these two remain unchanged.For theC→A region, the larger its width, the stronger the control exerted on the atomistic region by the continuum region and the more obvious effect CCLD has on the data exchanging in this region. For thebuffer region, its function is to alleviate the fluctuations among other kernel functional regions. Increasing its width means increasing the computation load, so it is better to minimize the width of thebuffer region or to remove it. For the A→C region, due to the existence of statistical thermal fluctuations in the atomistic region, the scope of it apparently influences the sampling average results and thus influences its capability to provide accurate boundary conditions to the continuum region.On the premise of accurate HAC simulation, different combinations of functional regions and different widths of them may have influences on simulation results. In Section4.2, we will test and discuss such influences and give suggestions about how to configure the overlap region properly.


## 3.2. Parameters of the Data Exchanging Operation:N and ex_steps

In Section3.1, we have talked about spatial average in the A→C region, while in this section we consider the temporal average parameters of the A→C operation. Generally speaking, in order to correctly extract macroscopic velocity, averages are taken over microscopic variables within a control volume and during time, which is called the binning method [35]. Specifically, in our HAC model, the control volume is the bin i in the A→C region and the time interval for averaging is Δtave for the continuum boundary extraction. Timely and correctly exchanging sampling atomistic data to the continuum solver is a kernel operation in the HAC method. We must pay attention that the temporal average over all sample points will introduce time lagging boundary results and would not exchange coupling information timely. Due to the overload of sampling operations, we should also question whether decreasing the times of coupling data exchanging operations could further reduce the sampling times and the computational load within acceptable errors.There are several parameters involved in the temporal average operation:(i)
ex_steps: the number of CFD time steps corresponding to one HAC data exchange operation Δtave, and Δtave=ex_steps×ΔtCFD(ii)
N: the number of sampling MD points in one data exchange operation(iii)
interval: the interval of sampling points in one data exchange operationIn the previous research [7, 8], data was often exchanged between the continuum region and the atomistic region in each CFD time step, and temporal average is performed over all M time steps in one ΔtCFD; that is, ex_steps=1, N=M, and interval=1.In Section4.3, we will discuss and test the proper sampling number N and exchange times ex_steps in one data exchange operation. Furthermore, we will provide guidance on the macroscopic velocity extraction.


## 3.3. Parameters of the Nonperiodic Boundary Force Model:α and β

Due to the lack of interaction with missing atoms, the atoms near the nonperiodic boundary will show unphysical effects. A nonperiodic boundary is applied to render density inhomogeneous and force nonuniform near the MD-continuum interface. In the existing HAC models, the nonperiodic boundary is often remedied with an external force model on those atoms near the boundary. This remedy operation will keep the pressure in the atomistic region correctly and alleviate the density fluctuation.We construct thenpbf region near the MD-continuum interface within a cutoff radius distance. Then, we refine thenpbf region intonpbf bins with width of binynpbf on the basis of CFD mesh width Lycell in y direction; that is, Lycell=β×binynpbf, where β is the finer parameter. Issa and Poesio [13] only provide a nonperiodic boundary force model with an empirical smoothing parameter α=1×10-4. A proper combination of the smoothing parameter α and the finer parameter β to minimize disturbance on the local liquid structure near the MD-continuum interface is worth investigating. In Section 4.4, we will provide numerical verifications on different combinations of α and β and discuss the appropriate values of them.


## 4. Results and Discussion

In this section, we use the model Couette flow with either no-slip or slip boundary condition to verify the effects of choosing different parameters in coupling strategies on the HAC simulation. We run the simulations on a high performance cluster where each computing node contains 12 Intel Xeon E5-2620 2.10 GHz CPU cores and a total main memory of 16 GB [36]. Section 4.1 validates the correctness of our HAC solver; Section 4.2 tests the effects of different spatial configurations of the overlap region on the a Couette flow; Section 4.3 investigates the effects of different temporal strategies of data exchanging on the model flow and Section 4.4 discusses the proper combination of parameters in the nonperiodic boundary force model.


### 4.1. Verification and Efficiency

In this section, two benchmark cases are carried out to test the validity and practical performance of our HAC solver by comparing with either the analytical solution or the full MD results.


#### 4.1.1. No-Slip Sudden Start Couette Flow

We firstly consider the typical no-slip sudden start Couette flow proposed by O’Connell and Thompson [6] as shown in Figure 1. The analytical solution for a sudden start Couette flow is given by(12)uy,t=uwallyLy+2uwallπ∑n=1∞cos⁡nπn×sinnπyLyexp⁡-νn2π2tLy2,where Ly is the distance between the two walls, uwall is the sliding velocity, and ν is the kinematic viscosity. We compare the resulting velocity profiles from our HAC model with the analytical solution.The height of the channel isH=44σ as used by Yen et al. [8]. The sliding velocity is Uw=1.0σ/τ at y/H=Ly. The simulation domain is spilt into the atomistic region near the lower wall and the continuum region near the upper wall. In this case, the no-slip boundary is set for the sliding wall in the continuum region as well as the stationary wall in the atomistic region. The wall-fluid interaction parameters σwf/σ=1.0, ϵwf/ϵ=0.6, and ρwf/ρ=1.0 are employed as used by Thompson and Troian [29]. The height of the pure continuum region, the pure MD region, and the overlap region are 24σ, 12σ, and 8σ, respectively, and the combination of functional regions has 4 layers including thecontrol region, the C→A region, thebuffer region, and the A→C region. The continuum region is divided into nx×ny=5×16 cells for numerical calculation; that is, Δx×Δy=3σ×2σ. The three-dimensional MD region is also divided into bins matching the cells in the continuum region, which is nx×ny×nz=5×11×1 bins and 10σ in z direction.Initially, the mean fluid is zero in the whole simulation domain. Att=0, the upper wall in the continuum region begins to move at Uw=1.0σ/τ, while the lower stationary wall in the atomistic region keeps still. The results are then averaged over the five time intervals as indicated in Figure 7. The transient velocity profiles of our HAC model match well with the analytical solution, especially in the overlap region, and the steady state profile is linear as expected. Therefore, this case can be used as to demonstrate the correctness of the proposed HAC model.Figure 7
Velocity evolution profiles averaged over five time intervals compared with the analytical solution.We also want to compare the efficiency of our HAC method to the full MD simulation using this benchmark case in serial mode. We list the detailed simulation time consuming in Table2. The total simulation time of the HAC method is only 39.25% of the full MD simulation which exhibits considerable efficiency.Table 2
Detailed time consuming for the HAC method and the full MD simulation.
Method MD domain (σ3) Number of particles CFD domain (σ3) Time (s) HAC 15×20×10 2592 15×32×10 3649.32 Full MD 15 × 44 × 10 5616 — 9295.87


#### 4.1.2. Slip Sudden Start Couette Flow

The second case is the slip sudden start Couette flow, with wall-fluid interaction parameters,σwf/σ=0.75, ϵwf/ϵ=0.6, and ρwf/ρ=4.0, as used by Thompson and Troian [29]. The velocity of sliding wall is Uw=1.0σ/τ. In this case, the height of the channel is H=44σ and there are two MD regions near the sliding wall and the stationary wall, respectively, while the continuum region is in the middle part of the channel as shown in Figure 2. The dimensions of the two MD regions and the continuum region are 22σ and 20σ, respectively. The partitions of cells and bins are the same as the previous case. We run independent realization of the same system ten times with the same configuration for both the HAC simulation and the full MD simulation for better thermal noise decrease as with Nie et al. [7]. The comparison of the evolutionary velocity profiles predicted by our HAC simulation and the full MD is presented in Figure 8. As we can see, the results of the two solutions agree quite well with each other with a small discrepancy in the evolution, and the deviation diminishes at the final steady state.Figure 8
Velocity evolution profiles averaged over four time intervals compared with the full MD results.


### 4.2. Effects of Functional Region Configurations on Couette Flows

Based on the discussion of the mutable variables mentioned in Section3.1, we test and discuss the effects of different combinations and different widths of the functional regions on the convergence and accuracy of Couette flow simulation here. The channel height of Couette flow is H=100σ and the sliding velocity is Uw=1.0σ/τ for both the no-slip and slip boundary conditions as shown in Figure 1 and Figure 2, while other test conditions are the same as those in Section 4.1.1 and Section 4.1.2.


#### 4.2.1. Effects of Different Functional Region Combinations on Couette Flows

In Section3.1, we have summarized the existing configurations of the overlap region. Five of the combinations are presented in Table 3 for the comparison of the influences on the simulation of Couette flow.Table 3
Details of different combinations of functional regions to be compared with.
Type of combinations Detailed combinations of functional regions 3-layer-1 control (npbf, mass flow), A→C, C→A 3-layer-2 control (npbf, mass flow, C→A), buffer, A→C 4-layer-1 control (npbf, mass flow), C→A, buffer, A→C 4-layer-2 control (npbf, mass flow), buffer, C→A, A→C 5-layer control (npbf, mass flow), buffer, C→A, buffer,A→CFirstly, we consider the tests with the no-slip boundary condition, with the same domain decomposition method in Section4.1.1, and with the extended channel height of 100σ. The width of functional regions is 2σ in y direction. Therefore, the total width of the overlap region in these five combinations ranges from 6σ to 10σ. The cumulative mean error between the results of the HAC simulation and the analytical solution is defined as(13)err=∑i=1NUiHAC-Uianalysis/UianalysisN,where N is the number of points in each resulting line, with the value of 50 when H=100σ. For it is a time-stepping simulation, we take five time intervals into account. In order to clearly distinguish between these five combinations, we depict the results of five intervals separately and average the total difference over these five intervals to depict theaverage deviation among these configurations using the first combination as the base, that is, “3-layer-1” as shown in Figure 9. In the following sections, we will focus on the minimumaverage deviation, that is, the last figure in Figure 9, and use the same kind of data conversion to clearly explain the distinction among different configurations.Figure 9
Cumulative mean errors during five time intervals (0–50)τ, (50–100)τ, (100–150)τ, (1000–1500)τ, and (7000–8000)τ and final average deviation.In the above figures, we find that, at the beginning of the HAC simulation, due to the thermal noise in the atomistic region, the results in the first three intervals deviate from the analytical solution greatly. However, this deviation vanishes when the simulation process reaches the final steady state.When the width of functional regions is2σ and the width of the overlap region is varying from 6σ to 10σ, the result of “5-layer” combination obtains the minimumaverage deviation compared with those of the other four. It is because of that “5-layer” is configured with the most atoms and the widest overlap region with the width of 10σ. In such combination, the fluctuations between the results from the two different numerical methods can be sufficiently alleviated, similar to the conclusion given by Yen et al. [8].Secondly, we consider the slip boundary condition of the above five combinations with parameters same as those of Section4.1.2. The slip length is defined as [37](14)Ls=Uw∂u/∂nw.In our model Couette flow, the simplified definition of (14) is given by [29](15)Ls=Uw/γ˙-H2,where γ˙ is the shearing rate and H is the channel height. We take the slip length at the stationary wall Lstationary and at the sliding wall Lsliding into consideration when the simulations reach the steady state and define the relative errors between the HAC simulation results and the full MD results as follows:(16)errstationary=LengthstationaryHAC-LengthstationaryfullMDLengthstationaryfullMD,errsliding=LengthslidingHAC-LengthslidingfullMDLengthslidingfullMDand the unified relative error between them is given by(17)errHAC=errstationary+errsliding2.The full MD results of the slip length at stationary wall and sliding wall areLstationary=2.6409σ and Lsliding=2.3619σ respectively, and the ideal slip length calculated by (15) is Ls=2.5014σ.The unified relative error ofLstationary and Lsliding is shown in Table 4. The HAC simulation results of these five combinations are all different from the full MD results, while the result of “5-layer” combination obtains the minimum error with 0.3735, which corresponds to the conclusion of that under the no-slip boundary condition.Table 4
Slip length prediction of five combinations and unified relative errors with different widths of the overlap region.
Situation 3-layer-1 3-layer-2 4-layer-1 4-layer-2 5-layer Lstationary 4.3012σ 4.8011σ 4.0739σ 4.9841σ 4.0947σ Lsliding 2.6735σ 1.7497σ 3.2287σ 3.5236σ 2.826σ errHAC 0.3804 0.5386 0.4548 0.68955 0.3735In the above two tests, the widths of the overlap regions are not the same with each other. Next, we carry out another two tests with the width of the overlap region fixed at10σ under the no-slip and slip boundary conditions. The width of each functional region is the same: therefore, the width of functional regions is 3.33σ in “3-layer” combination, 2.5σ in “4-layer” combination, and 2σ in “5-layer” combination.Using the same data processing method as the previous test, we depict the cumulative mean errors among five combinations andaverage deviation as shown in Figure 10.Figure 10
Cumulative mean errors during five time intervals (0–50)τ, (50–100)τ, (100–150)τ, (1000–1500)τ, and (7000–8000)τ and final average deviation.From Figure10, the result of “4-layer-1” combination achieves the minimumaverage derivation compared with “5-layer” one and also has a better time-stepping performance than “5-layer” one. We also test the slip boundary condition, and the relative errors of the slip length are listed in Table 5. The result of “4-layer-1” combination gains the slightest deviation from the full MD results with error of 0.3599.Table 5
Slip length prediction of five combinations and unified relative errors with the fixed width of the overlap region.
Situation 3-layer-1 3-layer-2 4-layer-1 4-layer-2 5-layer Lstationary 4.8303σ 4.7316σ 4.1067σ 4.4892σ 4.0947σ Lsliding 2.9541σ 3.8346σ 2.7511σ 3.2528σ 2.826σ errHAC 0.53985 0.7076 0.3599 0.53855 0.3735Based on the above four tests, we draw the following conclusion: under the condition of the equal width of the functional regions, when the width of the overlap region varies, “5-layer” combination obtains the best accuracy with minimumaverage deviation, while when the width of the overlap region is fixed, “4-layer-1” combination is the best configuration of these five combinations, which has a reasonable width of the functional regions and alleviates the fluctuations sufficiently. Furthermore, we draw the conclusion that it is reasonable to set abuffer region located between the C→A region and the A→C region, which can effectively relax the data exchanging between the continuum region and the atomistic region, while thebuffer region between thecontrol region and the C→A region could not be set so as to relieve the computational load.


#### 4.2.2. Effects of Different Widths of Functional Regions on Couette Flows

In Section4.2.1, the widths of functional regions are the same, while in this section we discuss the effects of different widths of functional regions on the Couette flow simulation. Following the discussion results of the previous section, we use “4-layer-1” type of combination to deploy the overlap region with the fixed width but to widen the C→A region, thebuffer region, and the A→C region separately. The settings of different widths are listed in Table 6. We choose another two situations from Section 4.2.1 for comparison.Table 6
Configuration situations of different widths of the functional regions.
Situation C → A region buffer region A → C region wider C → A 4σ 2σ 2σ wider buffer 2σ 4σ 2σ wider A → C 2σ 2σ 4σ compare-1 2σ 2σ 2σ compare-2 2.5σ 2.5σ 2.5σWe simulate the tests under the no-slip and slip boundary conditions and compare the simulation accuracy and convergence of these five situations. Under the no-slip boundary condition, the cumulative mean errors andaverage deviation are depicted in Figure 11.Figure 11
Cumulative mean errors during five time intervals (0–50)τ, (50–100)τ, (100–150)τ, (1000–1500)τ, and (7000–8000)τ and final average deviation.From Figure11, the results show that under the no-slip boundary condition “widerC→A” has a better accuracy than “compare-2.” For the slip boundary condition, the relative errors are listed in Table 7. It also shows that the “widerC→A” achieves the best accuracy among these five situations with error of 0.30942.Table 7
Slip length prediction of five situations and unified relative errors with different widths of the functional regions.
Situation wider C → A wider buffer wider A → C compare-1 compare-2 Lstationary 4.0375σ 5.2954σ 4.9102σ 4.0739σ 4.1067σ Lsliding 2.1705σ 3.2045σ 1.7241σ 3.2287σ 2.7511σ errHAC 0.30942 0.680925 0.5647 0.4548 0.3599In our model flow under the slip boundary condition, the atomistic regions provide more accuracy boundary conditions for the continuum region, while the continuum region just serves as the transmission container to transfer data between these two atomistic regions. In “widerC→A,” when enlarging the C→A region, the CCLD in the C→A region has more atoms to apply on. Therefore, the continuum region of “widerC→A” transfers data more efficiently and “widerC→A” predicts the slip length with the minimum error.From the above two tests, “widerC→A” that enlarges the C→A region performs best in matching with the analytical solution and predicting the slip length. Conclusions could be drawn that under the fixed width of the overlap region widening the C→A region improves simulation accuracy.


### 4.3. Effects of Sampling Parameters on Couette Flows

In this section, we test different groups of sampling and averaging parameters mentioned in Section3.2 under the no-slip and slip boundary conditions with channel height H=100σ as shown in Figures 1 and 2.


#### 4.3.1. Effects of Different Sampling Number on Couette Flows

We firstly setex_steps=1 and interval=1 and vary the sampling number N under the no-slip and slip boundary conditions. The parameter list is shown in Table 8.Table 8
Parameter list of different sampling numbers.
Situation N-1 N-2 N-4 N-8 N-16 N-32 N-64 compare N 1 2 4 8 16 32 64 100The first test is under the no-slip boundary and we depict the cumulative mean errors during the five time intervals andaverage deviation as shown in Figure 12. The result of “N-32” with sample number N=32 near the occasion of data exchanging obtains the minimum cumulative mean error. The second test is performed under the slip boundary condition and the relative errors are listed in Table 9.Table 9
Slip length prediction of eight situations and unified relative errors with different sampling numbers.
Situation N-1 N-2 N-4 N-8 N-16 N-32 N-64 compare Lstationary 4.1917σ 4.7359σ 4.4798σ 4.1362σ 4.1975σ 3.8371σ 4.4089σ 4.1067σ Lsliding 3.0493σ 2.459σ 3.1436σ 3.1479σ 3.1785σ 2.751σ 2.7515σ 2.7511σ errHAC 0.43915 0.4172 0.51365 0.4495 0.46755 0.30885 0.4172 0.3599Figure 12
Cumulative mean errors during five time intervals (0–50)τ, (50–100)τ, (100–150)τ, (1000–1500)τ, and (7000–8000)τ and final average deviation.Just as the conclusion of Section4.2.2 shows, under the current model flow, sampling points closer to the occasion of data exchanging provide more important messages for the simulation evolution. Therefore, “N-32” provides better data exchanging boundary conditions while ensuring the quantity of sampling points compared with other situations.From these two tests, we find that “N-32”, which provides 32 sampling points, achieves the best performance with the least error of 0.30885, offers the most effective information of the atomistic region, and also reduces the sampling times.


#### 4.3.2. Effects of Different Data Exchanging Times on Couette Flows

The data exchanging timesex_steps influence the total number of the C→A and the A→C operations in one HAC simulation. We discuss the proper times of data exchanging in the current section.Firstly, we change the value ofex_steps, sample, and average over all MD time steps in one Δtave under the no-slip boundary condition to check the effects of ex_steps on Couette flows. The parameters are listed in Table 10.Table 10
Parameter list of different data exchanging times and sampling numbers.
Situation ex-2 ex-4 ex-8 ex-16 ex-32 compare ex_steps 2 4 8 16 32 1 N 200 400 800 1600 3200 100In this section, the time intervals are multiples of the data exchanging time, and we consider four time intervals in this test. The cumulative mean errors andaverage deviation are plotted in Figure 13.Figure 13
Cumulative mean errors during four time intervals (64–128)τ, (256–512)τ, (960–1472)τ, and (6976–8000)τ and final average deviation.We can figure out that, when theex_steps changes, the accuracy of the first five situations is lower than the last one with ex_steps=1. In Figure 13, the accuracy of “ex-2” and “ex-4” is very much closer to “compare” with only 1% difference, while the total times of data exchanging operations are a half or a quarter of the last one.Secondly, we test the slip boundary condition, and the results are listed in Table11 with different sampling numbers; “compare” reaches the highest accuracy, but “ex-2” and “ex-4” have only deviated from “compare” with a percentage of 5, that is, 0.36691 and 0.3993, respectively.Table 11
Slip length prediction of six situations and unified relative errors with differentex_steps and N.
Situation ex-2 ex-4 ex-8 ex-16 ex-32 compare Lstationary 4.092σ 4.3663σ 4.5941σ 5.2711σ 4.5607σ 4.1067σ Lsliding 2.7973σ 2.708σ 2.3235σ 2.0588σ 2.5892σ 2.7511σ errHAC 0.36691 0.39994 0.37795 0.5621 0.41155 0.3599In the above two tests, the total number of sampling points is different. As we have discussed in Section4.3, sampling number with N=32 performs better than all MD time steps averaging. Next, we take another two cases into account with the same sampling points. The parameters are shown in Table 12.Table 12
Parameter list of different exchanging times but of equal sampling number.
Situation ex-2 ex-4 ex-8 ex-16 ex-32 compare ex_steps 2 4 8 16 32 1 N 32 32 32 32 32 32We perform the tests under the no-slip and slip boundary conditions and plot the normalized errors under the no-slip boundary condition in Figure14 and the relative errors under the slip boundary condition in Table 13.Table 13
Slip length prediction of six situations and unified relative errors with differentex_steps and equal N.
Situation ex-2 ex-4 ex-8 ex-16 ex-32 compare Lstationary 3.8645σ 3.8896σ 4.8903σ 4.1258σ 4.3519σ 3.8371σ Lsliding 2.8012σ 2.7598σ 2.2585σ 2.6894σ 3.4859σ 2.751σ errHAC 0.32645 0.32065 0.4478 0.3505 0.5619 0.30885Figure 14
Cumulative mean errors during four time intervals (64–128)τ, (256–512)τ, (960–1472)τ, and (6976–8000)τ and final average deviation.From Figure14 and Table 13, we can see that the first two situations perform better than “ex-8,” “ex-16,” and “ex-32” but worse than the last one under the no-slip boundary condition as well as the slip boundary condition. But the first two situations are worse than the last situation only by 1% and 2% in these two tests with error of 0.32645 and 0.32065, respectively.From the results of all the tests, we can conclude that timely data exchanging between the atomistic region and the continuum region performs better than otherex_steps settings. But with acceptable errors, we can choose a half or a quarter of the times of data exchanging to reduce the computational load. Therefore, “ex-2” and “ex-4” can be better choices.


### 4.4. Effects of Parameters of the Nonperiodic Boundary Force Model on Couette Flows

The effects of the smoothing parameterα and the finer parameter β on Couette flow simulation are discussed in this section. We use two tests under the no-slip and slip boundary conditions to measure these effects as shown in Figures 1 and 2. The combination of parameters are listed in Table 14. These 16 situations are tested to pick out the best combination of parameters with the least disturbance on the local liquid structure near the MD-continuum interface.Table 14
Combinations of parameters in the nonperiodic boundary force model with 16 situations.
Situation STNI-1 STNI-2 STNI-3 STNI-4 α 1×10-4 0.001 0.01 0.1 β 20 20 20 20  Situation STNI-5 STNI-6 STNI-7 STNI-8  α 1 × 1 0 - 4 0.001 0.01 0.1 β 10 10 10 10  Situation STNI-9 STNI-10 STNI-11 STNI-12  α 1 × 1 0 - 4 0.001 0.01 0.1 β 5 5 5 5  Situation STNI-13 STNI-14 STNI-15 STNI-16  α 1 × 1 0 - 4 0.001 0.01 0.1 β 1 1 1 1Figure15 shows the normalized cumulative mean errors under the no-slip boundary condition. Figures 16 and 17 show the relative errors of the slip length and density near the MD-continuum interface under the slip boundary condition.Figure 15
Differences among the sixteen situations using the results of “STNI-1” as the normalized base (a) and summation errors over the five time intervals (b).
(a) (b)Figure 16
Relative errors of slip length under 16 situations.Figure 17
Relative errors of density under 16 situations; the base density isρσ3=0.81.Figure17 shows that, with the decrease of the finer parameter β, the density near the MD-continuum interface drifts away from the base increasingly and “STNI-4” achieves the minimum disturbance on the local liquid structure. For Figure 15, “STNI-4” performs best in matching with the analytical solution and in Figure 16 “STNI-4” predicts the slip length better than the other situations except “STNI-6” and “STNI-10”, with 4% difference.Based on the above analysis, we can conclude that, in our current model flow, the smoothing parameterα and the finer parameter β can be chosen as 0.1 and 20, respectively, thus providing the best accuracy and prediction capability.


## 4.1. Verification and Efficiency

In this section, two benchmark cases are carried out to test the validity and practical performance of our HAC solver by comparing with either the analytical solution or the full MD results.


### 4.1.1. No-Slip Sudden Start Couette Flow

We firstly consider the typical no-slip sudden start Couette flow proposed by O’Connell and Thompson [6] as shown in Figure 1. The analytical solution for a sudden start Couette flow is given by(12)uy,t=uwallyLy+2uwallπ∑n=1∞cos⁡nπn×sinnπyLyexp⁡-νn2π2tLy2,where Ly is the distance between the two walls, uwall is the sliding velocity, and ν is the kinematic viscosity. We compare the resulting velocity profiles from our HAC model with the analytical solution.The height of the channel isH=44σ as used by Yen et al. [8]. The sliding velocity is Uw=1.0σ/τ at y/H=Ly. The simulation domain is spilt into the atomistic region near the lower wall and the continuum region near the upper wall. In this case, the no-slip boundary is set for the sliding wall in the continuum region as well as the stationary wall in the atomistic region. The wall-fluid interaction parameters σwf/σ=1.0, ϵwf/ϵ=0.6, and ρwf/ρ=1.0 are employed as used by Thompson and Troian [29]. The height of the pure continuum region, the pure MD region, and the overlap region are 24σ, 12σ, and 8σ, respectively, and the combination of functional regions has 4 layers including thecontrol region, the C→A region, thebuffer region, and the A→C region. The continuum region is divided into nx×ny=5×16 cells for numerical calculation; that is, Δx×Δy=3σ×2σ. The three-dimensional MD region is also divided into bins matching the cells in the continuum region, which is nx×ny×nz=5×11×1 bins and 10σ in z direction.Initially, the mean fluid is zero in the whole simulation domain. Att=0, the upper wall in the continuum region begins to move at Uw=1.0σ/τ, while the lower stationary wall in the atomistic region keeps still. The results are then averaged over the five time intervals as indicated in Figure 7. The transient velocity profiles of our HAC model match well with the analytical solution, especially in the overlap region, and the steady state profile is linear as expected. Therefore, this case can be used as to demonstrate the correctness of the proposed HAC model.Figure 7
Velocity evolution profiles averaged over five time intervals compared with the analytical solution.We also want to compare the efficiency of our HAC method to the full MD simulation using this benchmark case in serial mode. We list the detailed simulation time consuming in Table2. The total simulation time of the HAC method is only 39.25% of the full MD simulation which exhibits considerable efficiency.Table 2
Detailed time consuming for the HAC method and the full MD simulation.
Method MD domain (σ3) Number of particles CFD domain (σ3) Time (s) HAC 15×20×10 2592 15×32×10 3649.32 Full MD 15 × 44 × 10 5616 — 9295.87


### 4.1.2. Slip Sudden Start Couette Flow

The second case is the slip sudden start Couette flow, with wall-fluid interaction parameters,σwf/σ=0.75, ϵwf/ϵ=0.6, and ρwf/ρ=4.0, as used by Thompson and Troian [29]. The velocity of sliding wall is Uw=1.0σ/τ. In this case, the height of the channel is H=44σ and there are two MD regions near the sliding wall and the stationary wall, respectively, while the continuum region is in the middle part of the channel as shown in Figure 2. The dimensions of the two MD regions and the continuum region are 22σ and 20σ, respectively. The partitions of cells and bins are the same as the previous case. We run independent realization of the same system ten times with the same configuration for both the HAC simulation and the full MD simulation for better thermal noise decrease as with Nie et al. [7]. The comparison of the evolutionary velocity profiles predicted by our HAC simulation and the full MD is presented in Figure 8. As we can see, the results of the two solutions agree quite well with each other with a small discrepancy in the evolution, and the deviation diminishes at the final steady state.Figure 8
Velocity evolution profiles averaged over four time intervals compared with the full MD results.


## 4.1.1. No-Slip Sudden Start Couette Flow

We firstly consider the typical no-slip sudden start Couette flow proposed by O’Connell and Thompson [6] as shown in Figure 1. The analytical solution for a sudden start Couette flow is given by(12)uy,t=uwallyLy+2uwallπ∑n=1∞cos⁡nπn×sinnπyLyexp⁡-νn2π2tLy2,where Ly is the distance between the two walls, uwall is the sliding velocity, and ν is the kinematic viscosity. We compare the resulting velocity profiles from our HAC model with the analytical solution.The height of the channel isH=44σ as used by Yen et al. [8]. The sliding velocity is Uw=1.0σ/τ at y/H=Ly. The simulation domain is spilt into the atomistic region near the lower wall and the continuum region near the upper wall. In this case, the no-slip boundary is set for the sliding wall in the continuum region as well as the stationary wall in the atomistic region. The wall-fluid interaction parameters σwf/σ=1.0, ϵwf/ϵ=0.6, and ρwf/ρ=1.0 are employed as used by Thompson and Troian [29]. The height of the pure continuum region, the pure MD region, and the overlap region are 24σ, 12σ, and 8σ, respectively, and the combination of functional regions has 4 layers including thecontrol region, the C→A region, thebuffer region, and the A→C region. The continuum region is divided into nx×ny=5×16 cells for numerical calculation; that is, Δx×Δy=3σ×2σ. The three-dimensional MD region is also divided into bins matching the cells in the continuum region, which is nx×ny×nz=5×11×1 bins and 10σ in z direction.Initially, the mean fluid is zero in the whole simulation domain. Att=0, the upper wall in the continuum region begins to move at Uw=1.0σ/τ, while the lower stationary wall in the atomistic region keeps still. The results are then averaged over the five time intervals as indicated in Figure 7. The transient velocity profiles of our HAC model match well with the analytical solution, especially in the overlap region, and the steady state profile is linear as expected. Therefore, this case can be used as to demonstrate the correctness of the proposed HAC model.Figure 7
Velocity evolution profiles averaged over five time intervals compared with the analytical solution.We also want to compare the efficiency of our HAC method to the full MD simulation using this benchmark case in serial mode. We list the detailed simulation time consuming in Table2. The total simulation time of the HAC method is only 39.25% of the full MD simulation which exhibits considerable efficiency.Table 2
Detailed time consuming for the HAC method and the full MD simulation.
Method MD domain (σ3) Number of particles CFD domain (σ3) Time (s) HAC 15×20×10 2592 15×32×10 3649.32 Full MD 15 × 44 × 10 5616 — 9295.87


## 4.1.2. Slip Sudden Start Couette Flow

The second case is the slip sudden start Couette flow, with wall-fluid interaction parameters,σwf/σ=0.75, ϵwf/ϵ=0.6, and ρwf/ρ=4.0, as used by Thompson and Troian [29]. The velocity of sliding wall is Uw=1.0σ/τ. In this case, the height of the channel is H=44σ and there are two MD regions near the sliding wall and the stationary wall, respectively, while the continuum region is in the middle part of the channel as shown in Figure 2. The dimensions of the two MD regions and the continuum region are 22σ and 20σ, respectively. The partitions of cells and bins are the same as the previous case. We run independent realization of the same system ten times with the same configuration for both the HAC simulation and the full MD simulation for better thermal noise decrease as with Nie et al. [7]. The comparison of the evolutionary velocity profiles predicted by our HAC simulation and the full MD is presented in Figure 8. As we can see, the results of the two solutions agree quite well with each other with a small discrepancy in the evolution, and the deviation diminishes at the final steady state.Figure 8
Velocity evolution profiles averaged over four time intervals compared with the full MD results.


## 4.2. Effects of Functional Region Configurations on Couette Flows

Based on the discussion of the mutable variables mentioned in Section3.1, we test and discuss the effects of different combinations and different widths of the functional regions on the convergence and accuracy of Couette flow simulation here. The channel height of Couette flow is H=100σ and the sliding velocity is Uw=1.0σ/τ for both the no-slip and slip boundary conditions as shown in Figure 1 and Figure 2, while other test conditions are the same as those in Section 4.1.1 and Section 4.1.2.


### 4.2.1. Effects of Different Functional Region Combinations on Couette Flows

In Section3.1, we have summarized the existing configurations of the overlap region. Five of the combinations are presented in Table 3 for the comparison of the influences on the simulation of Couette flow.Table 3
Details of different combinations of functional regions to be compared with.
Type of combinations Detailed combinations of functional regions 3-layer-1 control (npbf, mass flow), A→C, C→A 3-layer-2 control (npbf, mass flow, C→A), buffer, A→C 4-layer-1 control (npbf, mass flow), C→A, buffer, A→C 4-layer-2 control (npbf, mass flow), buffer, C→A, A→C 5-layer control (npbf, mass flow), buffer, C→A, buffer,A→CFirstly, we consider the tests with the no-slip boundary condition, with the same domain decomposition method in Section4.1.1, and with the extended channel height of 100σ. The width of functional regions is 2σ in y direction. Therefore, the total width of the overlap region in these five combinations ranges from 6σ to 10σ. The cumulative mean error between the results of the HAC simulation and the analytical solution is defined as(13)err=∑i=1NUiHAC-Uianalysis/UianalysisN,where N is the number of points in each resulting line, with the value of 50 when H=100σ. For it is a time-stepping simulation, we take five time intervals into account. In order to clearly distinguish between these five combinations, we depict the results of five intervals separately and average the total difference over these five intervals to depict theaverage deviation among these configurations using the first combination as the base, that is, “3-layer-1” as shown in Figure 9. In the following sections, we will focus on the minimumaverage deviation, that is, the last figure in Figure 9, and use the same kind of data conversion to clearly explain the distinction among different configurations.Figure 9
Cumulative mean errors during five time intervals (0–50)τ, (50–100)τ, (100–150)τ, (1000–1500)τ, and (7000–8000)τ and final average deviation.In the above figures, we find that, at the beginning of the HAC simulation, due to the thermal noise in the atomistic region, the results in the first three intervals deviate from the analytical solution greatly. However, this deviation vanishes when the simulation process reaches the final steady state.When the width of functional regions is2σ and the width of the overlap region is varying from 6σ to 10σ, the result of “5-layer” combination obtains the minimumaverage deviation compared with those of the other four. It is because of that “5-layer” is configured with the most atoms and the widest overlap region with the width of 10σ. In such combination, the fluctuations between the results from the two different numerical methods can be sufficiently alleviated, similar to the conclusion given by Yen et al. [8].Secondly, we consider the slip boundary condition of the above five combinations with parameters same as those of Section4.1.2. The slip length is defined as [37](14)Ls=Uw∂u/∂nw.In our model Couette flow, the simplified definition of (14) is given by [29](15)Ls=Uw/γ˙-H2,where γ˙ is the shearing rate and H is the channel height. We take the slip length at the stationary wall Lstationary and at the sliding wall Lsliding into consideration when the simulations reach the steady state and define the relative errors between the HAC simulation results and the full MD results as follows:(16)errstationary=LengthstationaryHAC-LengthstationaryfullMDLengthstationaryfullMD,errsliding=LengthslidingHAC-LengthslidingfullMDLengthslidingfullMDand the unified relative error between them is given by(17)errHAC=errstationary+errsliding2.The full MD results of the slip length at stationary wall and sliding wall areLstationary=2.6409σ and Lsliding=2.3619σ respectively, and the ideal slip length calculated by (15) is Ls=2.5014σ.The unified relative error ofLstationary and Lsliding is shown in Table 4. The HAC simulation results of these five combinations are all different from the full MD results, while the result of “5-layer” combination obtains the minimum error with 0.3735, which corresponds to the conclusion of that under the no-slip boundary condition.Table 4
Slip length prediction of five combinations and unified relative errors with different widths of the overlap region.
Situation 3-layer-1 3-layer-2 4-layer-1 4-layer-2 5-layer Lstationary 4.3012σ 4.8011σ 4.0739σ 4.9841σ 4.0947σ Lsliding 2.6735σ 1.7497σ 3.2287σ 3.5236σ 2.826σ errHAC 0.3804 0.5386 0.4548 0.68955 0.3735In the above two tests, the widths of the overlap regions are not the same with each other. Next, we carry out another two tests with the width of the overlap region fixed at10σ under the no-slip and slip boundary conditions. The width of each functional region is the same: therefore, the width of functional regions is 3.33σ in “3-layer” combination, 2.5σ in “4-layer” combination, and 2σ in “5-layer” combination.Using the same data processing method as the previous test, we depict the cumulative mean errors among five combinations andaverage deviation as shown in Figure 10.Figure 10
Cumulative mean errors during five time intervals (0–50)τ, (50–100)τ, (100–150)τ, (1000–1500)τ, and (7000–8000)τ and final average deviation.From Figure10, the result of “4-layer-1” combination achieves the minimumaverage derivation compared with “5-layer” one and also has a better time-stepping performance than “5-layer” one. We also test the slip boundary condition, and the relative errors of the slip length are listed in Table 5. The result of “4-layer-1” combination gains the slightest deviation from the full MD results with error of 0.3599.Table 5
Slip length prediction of five combinations and unified relative errors with the fixed width of the overlap region.
Situation 3-layer-1 3-layer-2 4-layer-1 4-layer-2 5-layer Lstationary 4.8303σ 4.7316σ 4.1067σ 4.4892σ 4.0947σ Lsliding 2.9541σ 3.8346σ 2.7511σ 3.2528σ 2.826σ errHAC 0.53985 0.7076 0.3599 0.53855 0.3735Based on the above four tests, we draw the following conclusion: under the condition of the equal width of the functional regions, when the width of the overlap region varies, “5-layer” combination obtains the best accuracy with minimumaverage deviation, while when the width of the overlap region is fixed, “4-layer-1” combination is the best configuration of these five combinations, which has a reasonable width of the functional regions and alleviates the fluctuations sufficiently. Furthermore, we draw the conclusion that it is reasonable to set abuffer region located between the C→A region and the A→C region, which can effectively relax the data exchanging between the continuum region and the atomistic region, while thebuffer region between thecontrol region and the C→A region could not be set so as to relieve the computational load.


### 4.2.2. Effects of Different Widths of Functional Regions on Couette Flows

In Section4.2.1, the widths of functional regions are the same, while in this section we discuss the effects of different widths of functional regions on the Couette flow simulation. Following the discussion results of the previous section, we use “4-layer-1” type of combination to deploy the overlap region with the fixed width but to widen the C→A region, thebuffer region, and the A→C region separately. The settings of different widths are listed in Table 6. We choose another two situations from Section 4.2.1 for comparison.Table 6
Configuration situations of different widths of the functional regions.
Situation C → A region buffer region A → C region wider C → A 4σ 2σ 2σ wider buffer 2σ 4σ 2σ wider A → C 2σ 2σ 4σ compare-1 2σ 2σ 2σ compare-2 2.5σ 2.5σ 2.5σWe simulate the tests under the no-slip and slip boundary conditions and compare the simulation accuracy and convergence of these five situations. Under the no-slip boundary condition, the cumulative mean errors andaverage deviation are depicted in Figure 11.Figure 11
Cumulative mean errors during five time intervals (0–50)τ, (50–100)τ, (100–150)τ, (1000–1500)τ, and (7000–8000)τ and final average deviation.From Figure11, the results show that under the no-slip boundary condition “widerC→A” has a better accuracy than “compare-2.” For the slip boundary condition, the relative errors are listed in Table 7. It also shows that the “widerC→A” achieves the best accuracy among these five situations with error of 0.30942.Table 7
Slip length prediction of five situations and unified relative errors with different widths of the functional regions.
Situation wider C → A wider buffer wider A → C compare-1 compare-2 Lstationary 4.0375σ 5.2954σ 4.9102σ 4.0739σ 4.1067σ Lsliding 2.1705σ 3.2045σ 1.7241σ 3.2287σ 2.7511σ errHAC 0.30942 0.680925 0.5647 0.4548 0.3599In our model flow under the slip boundary condition, the atomistic regions provide more accuracy boundary conditions for the continuum region, while the continuum region just serves as the transmission container to transfer data between these two atomistic regions. In “widerC→A,” when enlarging the C→A region, the CCLD in the C→A region has more atoms to apply on. Therefore, the continuum region of “widerC→A” transfers data more efficiently and “widerC→A” predicts the slip length with the minimum error.From the above two tests, “widerC→A” that enlarges the C→A region performs best in matching with the analytical solution and predicting the slip length. Conclusions could be drawn that under the fixed width of the overlap region widening the C→A region improves simulation accuracy.


## 4.2.1. Effects of Different Functional Region Combinations on Couette Flows

In Section3.1, we have summarized the existing configurations of the overlap region. Five of the combinations are presented in Table 3 for the comparison of the influences on the simulation of Couette flow.Table 3
Details of different combinations of functional regions to be compared with.
Type of combinations Detailed combinations of functional regions 3-layer-1 control (npbf, mass flow), A→C, C→A 3-layer-2 control (npbf, mass flow, C→A), buffer, A→C 4-layer-1 control (npbf, mass flow), C→A, buffer, A→C 4-layer-2 control (npbf, mass flow), buffer, C→A, A→C 5-layer control (npbf, mass flow), buffer, C→A, buffer,A→CFirstly, we consider the tests with the no-slip boundary condition, with the same domain decomposition method in Section4.1.1, and with the extended channel height of 100σ. The width of functional regions is 2σ in y direction. Therefore, the total width of the overlap region in these five combinations ranges from 6σ to 10σ. The cumulative mean error between the results of the HAC simulation and the analytical solution is defined as(13)err=∑i=1NUiHAC-Uianalysis/UianalysisN,where N is the number of points in each resulting line, with the value of 50 when H=100σ. For it is a time-stepping simulation, we take five time intervals into account. In order to clearly distinguish between these five combinations, we depict the results of five intervals separately and average the total difference over these five intervals to depict theaverage deviation among these configurations using the first combination as the base, that is, “3-layer-1” as shown in Figure 9. In the following sections, we will focus on the minimumaverage deviation, that is, the last figure in Figure 9, and use the same kind of data conversion to clearly explain the distinction among different configurations.Figure 9
Cumulative mean errors during five time intervals (0–50)τ, (50–100)τ, (100–150)τ, (1000–1500)τ, and (7000–8000)τ and final average deviation.In the above figures, we find that, at the beginning of the HAC simulation, due to the thermal noise in the atomistic region, the results in the first three intervals deviate from the analytical solution greatly. However, this deviation vanishes when the simulation process reaches the final steady state.When the width of functional regions is2σ and the width of the overlap region is varying from 6σ to 10σ, the result of “5-layer” combination obtains the minimumaverage deviation compared with those of the other four. It is because of that “5-layer” is configured with the most atoms and the widest overlap region with the width of 10σ. In such combination, the fluctuations between the results from the two different numerical methods can be sufficiently alleviated, similar to the conclusion given by Yen et al. [8].Secondly, we consider the slip boundary condition of the above five combinations with parameters same as those of Section4.1.2. The slip length is defined as [37](14)Ls=Uw∂u/∂nw.In our model Couette flow, the simplified definition of (14) is given by [29](15)Ls=Uw/γ˙-H2,where γ˙ is the shearing rate and H is the channel height. We take the slip length at the stationary wall Lstationary and at the sliding wall Lsliding into consideration when the simulations reach the steady state and define the relative errors between the HAC simulation results and the full MD results as follows:(16)errstationary=LengthstationaryHAC-LengthstationaryfullMDLengthstationaryfullMD,errsliding=LengthslidingHAC-LengthslidingfullMDLengthslidingfullMDand the unified relative error between them is given by(17)errHAC=errstationary+errsliding2.The full MD results of the slip length at stationary wall and sliding wall areLstationary=2.6409σ and Lsliding=2.3619σ respectively, and the ideal slip length calculated by (15) is Ls=2.5014σ.The unified relative error ofLstationary and Lsliding is shown in Table 4. The HAC simulation results of these five combinations are all different from the full MD results, while the result of “5-layer” combination obtains the minimum error with 0.3735, which corresponds to the conclusion of that under the no-slip boundary condition.Table 4
Slip length prediction of five combinations and unified relative errors with different widths of the overlap region.
Situation 3-layer-1 3-layer-2 4-layer-1 4-layer-2 5-layer Lstationary 4.3012σ 4.8011σ 4.0739σ 4.9841σ 4.0947σ Lsliding 2.6735σ 1.7497σ 3.2287σ 3.5236σ 2.826σ errHAC 0.3804 0.5386 0.4548 0.68955 0.3735In the above two tests, the widths of the overlap regions are not the same with each other. Next, we carry out another two tests with the width of the overlap region fixed at10σ under the no-slip and slip boundary conditions. The width of each functional region is the same: therefore, the width of functional regions is 3.33σ in “3-layer” combination, 2.5σ in “4-layer” combination, and 2σ in “5-layer” combination.Using the same data processing method as the previous test, we depict the cumulative mean errors among five combinations andaverage deviation as shown in Figure 10.Figure 10
Cumulative mean errors during five time intervals (0–50)τ, (50–100)τ, (100–150)τ, (1000–1500)τ, and (7000–8000)τ and final average deviation.From Figure10, the result of “4-layer-1” combination achieves the minimumaverage derivation compared with “5-layer” one and also has a better time-stepping performance than “5-layer” one. We also test the slip boundary condition, and the relative errors of the slip length are listed in Table 5. The result of “4-layer-1” combination gains the slightest deviation from the full MD results with error of 0.3599.Table 5
Slip length prediction of five combinations and unified relative errors with the fixed width of the overlap region.
Situation 3-layer-1 3-layer-2 4-layer-1 4-layer-2 5-layer Lstationary 4.8303σ 4.7316σ 4.1067σ 4.4892σ 4.0947σ Lsliding 2.9541σ 3.8346σ 2.7511σ 3.2528σ 2.826σ errHAC 0.53985 0.7076 0.3599 0.53855 0.3735Based on the above four tests, we draw the following conclusion: under the condition of the equal width of the functional regions, when the width of the overlap region varies, “5-layer” combination obtains the best accuracy with minimumaverage deviation, while when the width of the overlap region is fixed, “4-layer-1” combination is the best configuration of these five combinations, which has a reasonable width of the functional regions and alleviates the fluctuations sufficiently. Furthermore, we draw the conclusion that it is reasonable to set abuffer region located between the C→A region and the A→C region, which can effectively relax the data exchanging between the continuum region and the atomistic region, while thebuffer region between thecontrol region and the C→A region could not be set so as to relieve the computational load.


## 4.2.2. Effects of Different Widths of Functional Regions on Couette Flows

In Section4.2.1, the widths of functional regions are the same, while in this section we discuss the effects of different widths of functional regions on the Couette flow simulation. Following the discussion results of the previous section, we use “4-layer-1” type of combination to deploy the overlap region with the fixed width but to widen the C→A region, thebuffer region, and the A→C region separately. The settings of different widths are listed in Table 6. We choose another two situations from Section 4.2.1 for comparison.Table 6
Configuration situations of different widths of the functional regions.
Situation C → A region buffer region A → C region wider C → A 4σ 2σ 2σ wider buffer 2σ 4σ 2σ wider A → C 2σ 2σ 4σ compare-1 2σ 2σ 2σ compare-2 2.5σ 2.5σ 2.5σWe simulate the tests under the no-slip and slip boundary conditions and compare the simulation accuracy and convergence of these five situations. Under the no-slip boundary condition, the cumulative mean errors andaverage deviation are depicted in Figure 11.Figure 11
Cumulative mean errors during five time intervals (0–50)τ, (50–100)τ, (100–150)τ, (1000–1500)τ, and (7000–8000)τ and final average deviation.From Figure11, the results show that under the no-slip boundary condition “widerC→A” has a better accuracy than “compare-2.” For the slip boundary condition, the relative errors are listed in Table 7. It also shows that the “widerC→A” achieves the best accuracy among these five situations with error of 0.30942.Table 7
Slip length prediction of five situations and unified relative errors with different widths of the functional regions.
Situation wider C → A wider buffer wider A → C compare-1 compare-2 Lstationary 4.0375σ 5.2954σ 4.9102σ 4.0739σ 4.1067σ Lsliding 2.1705σ 3.2045σ 1.7241σ 3.2287σ 2.7511σ errHAC 0.30942 0.680925 0.5647 0.4548 0.3599In our model flow under the slip boundary condition, the atomistic regions provide more accuracy boundary conditions for the continuum region, while the continuum region just serves as the transmission container to transfer data between these two atomistic regions. In “widerC→A,” when enlarging the C→A region, the CCLD in the C→A region has more atoms to apply on. Therefore, the continuum region of “widerC→A” transfers data more efficiently and “widerC→A” predicts the slip length with the minimum error.From the above two tests, “widerC→A” that enlarges the C→A region performs best in matching with the analytical solution and predicting the slip length. Conclusions could be drawn that under the fixed width of the overlap region widening the C→A region improves simulation accuracy.


## 4.3. Effects of Sampling Parameters on Couette Flows

In this section, we test different groups of sampling and averaging parameters mentioned in Section3.2 under the no-slip and slip boundary conditions with channel height H=100σ as shown in Figures 1 and 2.


### 4.3.1. Effects of Different Sampling Number on Couette Flows

We firstly setex_steps=1 and interval=1 and vary the sampling number N under the no-slip and slip boundary conditions. The parameter list is shown in Table 8.Table 8
Parameter list of different sampling numbers.
Situation N-1 N-2 N-4 N-8 N-16 N-32 N-64 compare N 1 2 4 8 16 32 64 100The first test is under the no-slip boundary and we depict the cumulative mean errors during the five time intervals andaverage deviation as shown in Figure 12. The result of “N-32” with sample number N=32 near the occasion of data exchanging obtains the minimum cumulative mean error. The second test is performed under the slip boundary condition and the relative errors are listed in Table 9.Table 9
Slip length prediction of eight situations and unified relative errors with different sampling numbers.
Situation N-1 N-2 N-4 N-8 N-16 N-32 N-64 compare Lstationary 4.1917σ 4.7359σ 4.4798σ 4.1362σ 4.1975σ 3.8371σ 4.4089σ 4.1067σ Lsliding 3.0493σ 2.459σ 3.1436σ 3.1479σ 3.1785σ 2.751σ 2.7515σ 2.7511σ errHAC 0.43915 0.4172 0.51365 0.4495 0.46755 0.30885 0.4172 0.3599Figure 12
Cumulative mean errors during five time intervals (0–50)τ, (50–100)τ, (100–150)τ, (1000–1500)τ, and (7000–8000)τ and final average deviation.Just as the conclusion of Section4.2.2 shows, under the current model flow, sampling points closer to the occasion of data exchanging provide more important messages for the simulation evolution. Therefore, “N-32” provides better data exchanging boundary conditions while ensuring the quantity of sampling points compared with other situations.From these two tests, we find that “N-32”, which provides 32 sampling points, achieves the best performance with the least error of 0.30885, offers the most effective information of the atomistic region, and also reduces the sampling times.


### 4.3.2. Effects of Different Data Exchanging Times on Couette Flows

The data exchanging timesex_steps influence the total number of the C→A and the A→C operations in one HAC simulation. We discuss the proper times of data exchanging in the current section.Firstly, we change the value ofex_steps, sample, and average over all MD time steps in one Δtave under the no-slip boundary condition to check the effects of ex_steps on Couette flows. The parameters are listed in Table 10.Table 10
Parameter list of different data exchanging times and sampling numbers.
Situation ex-2 ex-4 ex-8 ex-16 ex-32 compare ex_steps 2 4 8 16 32 1 N 200 400 800 1600 3200 100In this section, the time intervals are multiples of the data exchanging time, and we consider four time intervals in this test. The cumulative mean errors andaverage deviation are plotted in Figure 13.Figure 13
Cumulative mean errors during four time intervals (64–128)τ, (256–512)τ, (960–1472)τ, and (6976–8000)τ and final average deviation.We can figure out that, when theex_steps changes, the accuracy of the first five situations is lower than the last one with ex_steps=1. In Figure 13, the accuracy of “ex-2” and “ex-4” is very much closer to “compare” with only 1% difference, while the total times of data exchanging operations are a half or a quarter of the last one.Secondly, we test the slip boundary condition, and the results are listed in Table11 with different sampling numbers; “compare” reaches the highest accuracy, but “ex-2” and “ex-4” have only deviated from “compare” with a percentage of 5, that is, 0.36691 and 0.3993, respectively.Table 11
Slip length prediction of six situations and unified relative errors with differentex_steps and N.
Situation ex-2 ex-4 ex-8 ex-16 ex-32 compare Lstationary 4.092σ 4.3663σ 4.5941σ 5.2711σ 4.5607σ 4.1067σ Lsliding 2.7973σ 2.708σ 2.3235σ 2.0588σ 2.5892σ 2.7511σ errHAC 0.36691 0.39994 0.37795 0.5621 0.41155 0.3599In the above two tests, the total number of sampling points is different. As we have discussed in Section4.3, sampling number with N=32 performs better than all MD time steps averaging. Next, we take another two cases into account with the same sampling points. The parameters are shown in Table 12.Table 12
Parameter list of different exchanging times but of equal sampling number.
Situation ex-2 ex-4 ex-8 ex-16 ex-32 compare ex_steps 2 4 8 16 32 1 N 32 32 32 32 32 32We perform the tests under the no-slip and slip boundary conditions and plot the normalized errors under the no-slip boundary condition in Figure14 and the relative errors under the slip boundary condition in Table 13.Table 13
Slip length prediction of six situations and unified relative errors with differentex_steps and equal N.
Situation ex-2 ex-4 ex-8 ex-16 ex-32 compare Lstationary 3.8645σ 3.8896σ 4.8903σ 4.1258σ 4.3519σ 3.8371σ Lsliding 2.8012σ 2.7598σ 2.2585σ 2.6894σ 3.4859σ 2.751σ errHAC 0.32645 0.32065 0.4478 0.3505 0.5619 0.30885Figure 14
Cumulative mean errors during four time intervals (64–128)τ, (256–512)τ, (960–1472)τ, and (6976–8000)τ and final average deviation.From Figure14 and Table 13, we can see that the first two situations perform better than “ex-8,” “ex-16,” and “ex-32” but worse than the last one under the no-slip boundary condition as well as the slip boundary condition. But the first two situations are worse than the last situation only by 1% and 2% in these two tests with error of 0.32645 and 0.32065, respectively.From the results of all the tests, we can conclude that timely data exchanging between the atomistic region and the continuum region performs better than otherex_steps settings. But with acceptable errors, we can choose a half or a quarter of the times of data exchanging to reduce the computational load. Therefore, “ex-2” and “ex-4” can be better choices.


## 4.3.1. Effects of Different Sampling Number on Couette Flows

We firstly setex_steps=1 and interval=1 and vary the sampling number N under the no-slip and slip boundary conditions. The parameter list is shown in Table 8.Table 8
Parameter list of different sampling numbers.
Situation N-1 N-2 N-4 N-8 N-16 N-32 N-64 compare N 1 2 4 8 16 32 64 100The first test is under the no-slip boundary and we depict the cumulative mean errors during the five time intervals andaverage deviation as shown in Figure 12. The result of “N-32” with sample number N=32 near the occasion of data exchanging obtains the minimum cumulative mean error. The second test is performed under the slip boundary condition and the relative errors are listed in Table 9.Table 9
Slip length prediction of eight situations and unified relative errors with different sampling numbers.
Situation N-1 N-2 N-4 N-8 N-16 N-32 N-64 compare Lstationary 4.1917σ 4.7359σ 4.4798σ 4.1362σ 4.1975σ 3.8371σ 4.4089σ 4.1067σ Lsliding 3.0493σ 2.459σ 3.1436σ 3.1479σ 3.1785σ 2.751σ 2.7515σ 2.7511σ errHAC 0.43915 0.4172 0.51365 0.4495 0.46755 0.30885 0.4172 0.3599Figure 12
Cumulative mean errors during five time intervals (0–50)τ, (50–100)τ, (100–150)τ, (1000–1500)τ, and (7000–8000)τ and final average deviation.Just as the conclusion of Section4.2.2 shows, under the current model flow, sampling points closer to the occasion of data exchanging provide more important messages for the simulation evolution. Therefore, “N-32” provides better data exchanging boundary conditions while ensuring the quantity of sampling points compared with other situations.From these two tests, we find that “N-32”, which provides 32 sampling points, achieves the best performance with the least error of 0.30885, offers the most effective information of the atomistic region, and also reduces the sampling times.


## 4.3.2. Effects of Different Data Exchanging Times on Couette Flows

The data exchanging timesex_steps influence the total number of the C→A and the A→C operations in one HAC simulation. We discuss the proper times of data exchanging in the current section.Firstly, we change the value ofex_steps, sample, and average over all MD time steps in one Δtave under the no-slip boundary condition to check the effects of ex_steps on Couette flows. The parameters are listed in Table 10.Table 10
Parameter list of different data exchanging times and sampling numbers.
Situation ex-2 ex-4 ex-8 ex-16 ex-32 compare ex_steps 2 4 8 16 32 1 N 200 400 800 1600 3200 100In this section, the time intervals are multiples of the data exchanging time, and we consider four time intervals in this test. The cumulative mean errors andaverage deviation are plotted in Figure 13.Figure 13
Cumulative mean errors during four time intervals (64–128)τ, (256–512)τ, (960–1472)τ, and (6976–8000)τ and final average deviation.We can figure out that, when theex_steps changes, the accuracy of the first five situations is lower than the last one with ex_steps=1. In Figure 13, the accuracy of “ex-2” and “ex-4” is very much closer to “compare” with only 1% difference, while the total times of data exchanging operations are a half or a quarter of the last one.Secondly, we test the slip boundary condition, and the results are listed in Table11 with different sampling numbers; “compare” reaches the highest accuracy, but “ex-2” and “ex-4” have only deviated from “compare” with a percentage of 5, that is, 0.36691 and 0.3993, respectively.Table 11
Slip length prediction of six situations and unified relative errors with differentex_steps and N.
Situation ex-2 ex-4 ex-8 ex-16 ex-32 compare Lstationary 4.092σ 4.3663σ 4.5941σ 5.2711σ 4.5607σ 4.1067σ Lsliding 2.7973σ 2.708σ 2.3235σ 2.0588σ 2.5892σ 2.7511σ errHAC 0.36691 0.39994 0.37795 0.5621 0.41155 0.3599In the above two tests, the total number of sampling points is different. As we have discussed in Section4.3, sampling number with N=32 performs better than all MD time steps averaging. Next, we take another two cases into account with the same sampling points. The parameters are shown in Table 12.Table 12
Parameter list of different exchanging times but of equal sampling number.
Situation ex-2 ex-4 ex-8 ex-16 ex-32 compare ex_steps 2 4 8 16 32 1 N 32 32 32 32 32 32We perform the tests under the no-slip and slip boundary conditions and plot the normalized errors under the no-slip boundary condition in Figure14 and the relative errors under the slip boundary condition in Table 13.Table 13
Slip length prediction of six situations and unified relative errors with differentex_steps and equal N.
Situation ex-2 ex-4 ex-8 ex-16 ex-32 compare Lstationary 3.8645σ 3.8896σ 4.8903σ 4.1258σ 4.3519σ 3.8371σ Lsliding 2.8012σ 2.7598σ 2.2585σ 2.6894σ 3.4859σ 2.751σ errHAC 0.32645 0.32065 0.4478 0.3505 0.5619 0.30885Figure 14
Cumulative mean errors during four time intervals (64–128)τ, (256–512)τ, (960–1472)τ, and (6976–8000)τ and final average deviation.From Figure14 and Table 13, we can see that the first two situations perform better than “ex-8,” “ex-16,” and “ex-32” but worse than the last one under the no-slip boundary condition as well as the slip boundary condition. But the first two situations are worse than the last situation only by 1% and 2% in these two tests with error of 0.32645 and 0.32065, respectively.From the results of all the tests, we can conclude that timely data exchanging between the atomistic region and the continuum region performs better than otherex_steps settings. But with acceptable errors, we can choose a half or a quarter of the times of data exchanging to reduce the computational load. Therefore, “ex-2” and “ex-4” can be better choices.


## 4.4. Effects of Parameters of the Nonperiodic Boundary Force Model on Couette Flows

The effects of the smoothing parameterα and the finer parameter β on Couette flow simulation are discussed in this section. We use two tests under the no-slip and slip boundary conditions to measure these effects as shown in Figures 1 and 2. The combination of parameters are listed in Table 14. These 16 situations are tested to pick out the best combination of parameters with the least disturbance on the local liquid structure near the MD-continuum interface.Table 14
Combinations of parameters in the nonperiodic boundary force model with 16 situations.
Situation STNI-1 STNI-2 STNI-3 STNI-4 α 1×10-4 0.001 0.01 0.1 β 20 20 20 20  Situation STNI-5 STNI-6 STNI-7 STNI-8  α 1 × 1 0 - 4 0.001 0.01 0.1 β 10 10 10 10  Situation STNI-9 STNI-10 STNI-11 STNI-12  α 1 × 1 0 - 4 0.001 0.01 0.1 β 5 5 5 5  Situation STNI-13 STNI-14 STNI-15 STNI-16  α 1 × 1 0 - 4 0.001 0.01 0.1 β 1 1 1 1Figure15 shows the normalized cumulative mean errors under the no-slip boundary condition. Figures 16 and 17 show the relative errors of the slip length and density near the MD-continuum interface under the slip boundary condition.Figure 15
Differences among the sixteen situations using the results of “STNI-1” as the normalized base (a) and summation errors over the five time intervals (b).
(a) (b)Figure 16
Relative errors of slip length under 16 situations.Figure 17
Relative errors of density under 16 situations; the base density isρσ3=0.81.Figure17 shows that, with the decrease of the finer parameter β, the density near the MD-continuum interface drifts away from the base increasingly and “STNI-4” achieves the minimum disturbance on the local liquid structure. For Figure 15, “STNI-4” performs best in matching with the analytical solution and in Figure 16 “STNI-4” predicts the slip length better than the other situations except “STNI-6” and “STNI-10”, with 4% difference.Based on the above analysis, we can conclude that, in our current model flow, the smoothing parameterα and the finer parameter β can be chosen as 0.1 and 20, respectively, thus providing the best accuracy and prediction capability.


## 5. Conclusions

In this paper, we design a domain decomposition type of hybrid MD-continuum solver, using open source software LAMMPS and OpenFOAM. We useCouette channel flow as our model flow to investigate the coupling strategy issues. Focusing on the fixed channel height and the sliding velocity, under the no-slip and slip boundary conditions, we make a deep analysis of different combinations and different widths of the functional regions, different parameters of data exchanging, and various combinations of parameters in the nonperiodic boundary force model.Forlayernumber and layerwidth, we find that under the condition of the equal width of functional regions, when the width of the overlap region is varied, “5-layer” combination obtains the best accuracy, while when the width of the overlap region is fixed, “4-layer-1” combination is the best setting of these five combinations, which has a reasonable width of functional regions and alleviates the fluctuations sufficiently. Furthermore, we can give the conclusion that it is reasonable to set abuffer region located between the C→A region and the A→C region which can effectively relax the exchanging data between the continuum region and the atomistic region, while thebuffer region located between thecontrol region and the C→A region can be removed to save the computational load. We can also figure out that under the fixed width of the overlap region, widening the C→A region gets a better simulation result in our model flow.As to the sampling parameters of temporal average, the present results disclose that the data exchanging operation only needs a few sampling points which are close to the occasions of interactions to guarantee modeling efficiency and to reduce the sampling times, and also timely data exchanging between the atomistic region and the continuum region performs best beyond other settings, but with acceptable errors; a half or a quarter of the data exchanging times can be chosen to reduce the computational load. The discussion of parameters of the nonperiodic boundary force model on the model flow shows that, under the domain decomposition along the flow direction, the smoothing parameter value of 0.1 and the finer parameter value of 20 can achieve the minimum disturbance on the local structure near the MD-continuum interface and keep the simulation accuracy.In this paper, we mainly focus on the HAC method based on geometrical coupling. There are other coupling techniques, such as embedded coupling [38]. The microscopic method can be used to calibrate the parameters in the continuum mechanics model [39, 40]. In our future work, we aim to extend the simulation power of our framework to support more kinds of multiscale coupling.



---

*Source: 1014636-2017-02-14.xml*