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: the 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; 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; 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.
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 ( region) and a continuum-coupled-to-atomistic region ( 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 the mass flow region; other researchers [14, 19] brought forward the buffer 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 using Couette 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 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 the 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 Section 2 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 Figures 1 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 plane, while a three-dimensional simulation is performed in the atomistic regions with axis as the extension direction.
2.2. Atomistic Region and Continuum Region
In this section, we introduce numerical methods used in the atomistic region and the continuum region. Section 2.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 bywhere is the distance between atom and , 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 , , and , where is the atom mass and a well-defined liquid phase of Argon with , , and the dynamic viscosity , where is the Boltzmann constant and is the characteristic time of the Lennard-Jones potential. We choose a cutoff distance to save the computation time. The wall atoms are modeled by two (111) planes with the direction of the lattice along orientation. For the two-dimensional continuum simulation, there is no flow in direction. Therefore, in the atomistic region, and directions are under periodic boundary conditions while 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 direction which is perpendicular to the bulk flow direction and dose not influence the bulk flow velocity. The equation of motion for atom is given bywhere the summation donates all pair-interaction forces on atom and is the damping ratio. denotes a Gaussian distributed force with zero mean and variance of . We choose the damping ratio of Langevin thermostat 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 .
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 bywhere is the bulk flow velocity, 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 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, the region, the region, the buffer region, the nonperiodic boundary force region (npbf region), and the mass flow region. The latter two are often unified to the control region. The schematic of the overlap region is shown in Figure 3.
2.3.1. Region
The region is located on the outer boundary of the continuum region centered on HSI in Figure 3. In this region, the 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 and direction to match the mesh cells in the continuum region. For the th bin, the boundary condition for the continuum velocity u is given by the following average equation:where is the velocity of the th atom in the bin , is the number of atoms in the bin , 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 sampling points:
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:where is the calculated velocity on the HSI using CAM sampling technique, n is the normal vector to the boundary, is an element of the boundary, and is the whole boundary of the continuum region.
2.3.2. Region
The 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 region, the region can be divided into bins in and direction to match the mesh cells in continuum region.
The CCLD requires that the mean atomistic velocities in bin should be equal to the average continuum velocities in the bin ; that iswhere is the number of atoms in the bin . In this paper, the volume of the 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 ; that iswhere is , , or 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 in and direction while in 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].
Each atom traveling through 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 and use a feedback force to apply on atoms in the missing interactions region.
Firstly, we calculate the average external force experienced by bin ; that iswhere is the interaction force between atoms and , is the number of atoms in bin , and is the average external force experienced by atoms in bin . 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 . This reversed force is given as follows, which is constructed with the simple exponential smoothing, including a smoothing parameter :
In each MD time step, we calculate in npbf region bins, use (10) to construct the feedback force , 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 the mass flux across the MD-continuum interface, we bring forward the mass flow region and use it to deal with the atom insertion and deletion. The mass flow region can be split into several bins to meet with continuum cells. In one continuum time step , the number of atoms to be inserted into or deleted from the mass flow bin is given bywhere is the bin area normal to direction and n is the face normal vector, pointed into the atomistic region. is the continuum velocity, is the liquid density, and is the atom mass. If is negative, certain atoms should be deleted from the mass flow region and if is positive, certain atoms should be inserted into the mass 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 . 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 and the velocity normal to .
2.3.5. Buffer Region
The buffer region is located between the region and the region or between the control region and the region. We use it to alleviate the fluctuation of the artificial operations on atoms. In the buffer region, there is no atom insertion or deletion and no constraint on atom velocities. But we have to narrow the width of the buffer 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 the buffer 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 , integration time steps of Navier-Stokes equation , and sampling average time . As we have mentioned before, is chosen as in our paper. In order to solve NS equation accurately, firstly, must be far less than characteristic time on the mesh, that is, diffusion time , where is the dynamic viscosity and . Secondly, it must be larger than the decay time of velocity autocorrelation function so as to reduce thermal noise from HSI [6]. Finally, it should meet the CFL condition [34], . The CFD time step we choose in this paper is , where is 100. The time advancing mechanism in the HAC simulation is shown in Figure 6, which is sequential coupling.
In each coupling simulation time step, CFD advances a certain time, for example, one , and transfers continuum velocities to the atomistic region through CCLD in the operation. Then, MD advances the same time interval, that is, , samples and averages atom velocities, and then passes them to the continuum region in the operation, thus finishing one cycle of HAC simulation.
Indeed, due to the small time step , 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 Section 3.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 and Layer 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 Section 2.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 the layer number (LN) and the layer width (LW) in Table 1.
For the completeness of comparison, we add an alternative configuration which includes four layers, that is, the control region (the npbf region and the mass flow region), the buffer region, the region, and the region. Those configurations that contain the depress 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, the mass 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. The npbf region is located within a cutoff distance near the MD-continuum interface. The widths of these two remain unchanged.
For the 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 the buffer 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 the buffer region or to remove it. For the 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 Section 4.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: and
In Section 3.1, we have talked about spatial average in the region, while in this section we consider the temporal average parameters of the 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 in the region and the time interval for averaging is 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): the number of CFD time steps corresponding to one HAC data exchange operation , and (ii): the number of sampling MD points in one data exchange operation(iii): the interval of sampling points in one data exchange operation
In 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 time steps in one ; that is, , , and .
In Section 4.3, we will discuss and test the proper sampling number and exchange times 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 the npbf region near the MD-continuum interface within a cutoff radius distance. Then, we refine the npbf region into npbf bins with width of on the basis of CFD mesh width in direction; that is, , where is the finer parameter. Issa and Poesio [13] only provide a nonperiodic boundary force model with an empirical smoothing parameter . 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 bywhere is the distance between the two walls, 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 is as used by Yen et al. [8]. The sliding velocity is at . 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 , , and 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 , , and , respectively, and the combination of functional regions has 4 layers including the control region, the region, the buffer region, and the region. The continuum region is divided into cells for numerical calculation; that is, . The three-dimensional MD region is also divided into bins matching the cells in the continuum region, which is bins and in direction.
Initially, the mean fluid is zero in the whole simulation domain. At , the upper wall in the continuum region begins to move at , 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.
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 Table 2. The total simulation time of the HAC method is only of the full MD simulation which exhibits considerable efficiency.
4.1.2. Slip Sudden Start Couette Flow
The second case is the slip sudden start Couette flow, with wall-fluid interaction parameters, , , and , as used by Thompson and Troian [29]. The velocity of sliding wall is . In this case, the height of the channel is 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 and , 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.
4.2. Effects of Functional Region Configurations on Couette Flows
Based on the discussion of the mutable variables mentioned in Section 3.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 and the sliding velocity is 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 Section 3.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.
Firstly, we consider the tests with the no-slip boundary condition, with the same domain decomposition method in Section 4.1.1, and with the extended channel height of . The width of functional regions is in direction. Therefore, the total width of the overlap region in these five combinations ranges from to . The cumulative mean error between the results of the HAC simulation and the analytical solution is defined aswhere is the number of points in each resulting line, with the value of 50 when . 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 the average 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 minimum average 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.
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 is and the width of the overlap region is varying from to , the result of “5-layer” combination obtains the minimum average 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 . 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 Section 4.1.2. The slip length is defined as [37]
In our model Couette flow, the simplified definition of (14) is given by [29]where is the shearing rate and is the channel height. We take the slip length at the stationary wall and at the sliding wall 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:and the unified relative error between them is given by
The full MD results of the slip length at stationary wall and sliding wall are and respectively, and the ideal slip length calculated by (15) is .
The unified relative error of and 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.
In 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 at under the no-slip and slip boundary conditions. The width of each functional region is the same: therefore, the width of functional regions is in “3-layer” combination, in “4-layer” combination, and in “5-layer” combination.
Using the same data processing method as the previous test, we depict the cumulative mean errors among five combinations and average deviation as shown in Figure 10.
From Figure 10, the result of “4-layer-1” combination achieves the minimum average 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.
Based 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 minimum average 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 a buffer region located between the region and the region, which can effectively relax the data exchanging between the continuum region and the atomistic region, while the buffer region between the control region and the 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 Section 4.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 region, the buffer region, and the region separately. The settings of different widths are listed in Table 6. We choose another two situations from Section 4.2.1 for comparison.
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 and average deviation are depicted in Figure 11.
From Figure 11, the results show that under the no-slip boundary condition “wider ” 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 “wider ” achieves the best accuracy among these five situations with error of 0.30942.
In 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 “wider ,” when enlarging the region, the CCLD in the region has more atoms to apply on. Therefore, the continuum region of “wider ” transfers data more efficiently and “wider ” predicts the slip length with the minimum error.
From the above two tests, “wider ” that enlarges the 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 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 Section 3.2 under the no-slip and slip boundary conditions with channel height as shown in Figures 1 and 2.
4.3.1. Effects of Different Sampling Number on Couette Flows
We firstly set and and vary the sampling number under the no-slip and slip boundary conditions. The parameter list is shown in Table 8.
The first test is under the no-slip boundary and we depict the cumulative mean errors during the five time intervals and average deviation as shown in Figure 12. The result of “N-32” with sample number 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.
Just as the conclusion of Section 4.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 times influence the total number of the and the operations in one HAC simulation. We discuss the proper times of data exchanging in the current section.
Firstly, we change the value of , sample, and average over all MD time steps in one under the no-slip boundary condition to check the effects of on Couette flows. The parameters are listed in Table 10.
In 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 and average deviation are plotted in Figure 13.
We can figure out that, when the changes, the accuracy of the first five situations is lower than the last one with . In Figure 13, the accuracy of “ex-2” and “ex-4” is very much closer to “compare” with only 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 Table 11 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.
In the above two tests, the total number of sampling points is different. As we have discussed in Section 4.3, sampling number with 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.
We perform the tests under the no-slip and slip boundary conditions and plot the normalized errors under the no-slip boundary condition in Figure 14 and the relative errors under the slip boundary condition in Table 13.
From Figure 14 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 and 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 other 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.
Figure 15 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.
(a)
(b)
Figure 17 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 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 use Couette 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.
For and , 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 a buffer region located between the region and the region which can effectively relax the exchanging data between the continuum region and the atomistic region, while the buffer region located between the control region and the 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 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.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
The authors would like to thank the National Key Research and Development Program of China (no. 2016YFB0201301), Science Challenge Project (no. JCKY2016212A502), and the Open Fund from State Key Laboratory of High Performance Computing (nos. 201303-01, 201503-01, and 201503-02) for funding.