Molecular dynamics simulations in hybrid particle-continuum schemes: Pitfalls and caveats

TL;DR

This paper examines the challenges of molecular dynamics in hybrid multiscale methods, offering insights on simulation strategies, thermostat choices, finite-size effects, and proposing a new reduced-order coupling technique to improve efficiency.

Contribution

It provides a detailed analysis of molecular simulation pitfalls in HMM and introduces a novel reduced-order method for coupling microscopic and macroscopic models.

Findings

Simulating many small systems is preferable to few large ones.

Simple isokinetic thermostats are generally sufficient.

The proposed reduced-order technique efficiently approximates non-linear stress-strain relations.

Abstract

Heterogeneous multiscale methods (HMM) combine molecular accuracy of particle-based simulations with the computational efficiency of continuum descriptions to model flow in soft matter liquids. In these schemes, molecular simulations typically pose a computational bottleneck, which we investigate in detail in this study. We find that it is preferable to simulate many small systems as opposed to a few large systems, and that a choice of a simple isokinetic thermostat is typically sufficient while thermostats such as Lowe-Andersen allow for simulations at elevated viscosity. We discuss suitable choices for time steps and finite-size effects which arise in the limit of very small simulation boxes. We also argue that if colloidal systems are considered as opposed to atomistic systems, the gap between microscopic and macroscopic simulations regarding time and length scales is significantly smaller. We also propose a novel reduced-order technique for the coupling to the macroscopic solver, which allows us to approximate a non-linear stress-strain relation efficiently and thus further reduce computational effort of microscopic simulations.

Equations

$${\sigma }_{\alpha \beta }=-\frac{1}{V}\left(\sum _i^N\left(m_i{v}_{i,\alpha }{v}_{i,\beta }\right)+\sum _i^N\sum _{j>i}^N\left(r_{ij,\alpha }{F}_{ij,\beta }\right)\right),$$

$$p=-\frac{\text{Tr}\sigma }{3}.$$

$$\eta =\left|\frac{{\sigma }_{13}}{\dot{\gamma }}\right|,$$

$$V_{WCA}(r)=4\epsilon \left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6+\frac{1}{4}\right],r<r_C=\sqrt[6]{2}\sigma .$$

$${\mathbf{v}}_{ij}^{\prime }={\xi }_{ij}\sqrt{2k_{B}T}{\hat{r}}_{ij}$$

$${\mathbf{v}}_i^{\prime }={\mathbf{v}}_i+{\mathbf{\Delta }}_{ij}$$

$${\mathbf{v}}_j^{\prime }={\mathbf{v}}_j-{\mathbf{\Delta }}_{ij},$$

$$2{\mathbf{\Delta }}_{ij}={\hat{r}}_{ij}\left({\mathbf{v}}_{ij}^{\prime }-{\mathbf{v}}_{ij}\right)\cdot {\hat{r}}_{ij}$$

$$\frac{3}{2}N{k}_{B}T=\frac{1}{2}\sum _i^N{m}_i{v}_i^2$$

$$\lambda =\sqrt{\frac{T_0}{T}}$$

$$\dot{\mathbf{q}}=\mathbf{p}+▽\mathbf{u}\cdot \mathbf{q}$$

$$\dot{\mathbf{p}}=\frac{\mathbf{f}}{m}-▽\mathbf{u}\cdot \mathbf{p}$$

$$\bigtriangledown\mathbf{u}=\left(\begin{array}[]{ccc}0&0&0\\ 0&0&0\\ \dot{\gamma}&0&0\end{array}\right)$$

$${\mathbf{q}}^{n+\frac{1}{2}}$$

$${\mathbf{p}}^{n+\frac{1}{2}}$$

$${\mathbf{q}}^{n+1}$$

$${\mathbf{F}}^{n+1}$$

$${\mathbf{p}}^{n+1}$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)=\frac{1}{Re}\nabla \cdot \mathbf{\sigma }+\mathbf{g},$$

$$\mathbf{u}={\mathbf{u}}_D,$$

$$\mathbf{u}(t=0)={\mathbf{u}}^{(0)}$$

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)=-\nabla p+\frac{1}{Re}(\mu \Delta \mathbf{u})+\mathbf{g}\text{in}\Omega .$$

$$\displaystyle\tilde{\mathbf{u}}=\Delta t\big{(}-\mathbf{u}^{(n)}\cdot\nabla\mathbf{u}^{(n)}+\frac{1}{Re}\nabla\cdot(\boldsymbol{\sigma}^{(n)}/\rho)+\mathbf{g}^{(n)}/\rho\big{)},\quad$$

$$\displaystyle\Delta\bar{p}^{(n+1)}=\nabla\cdot\big{(}\tilde{\mathbf{u}}/\Delta t\big{)},$$

$$\displaystyle\frac{\partial\bar{p}}{\partial\mathbf{n}}^{(n+1)}=\mathbf{n}\cdot\big{(}\tilde{\mathbf{u}}/\Delta t\big{)},$$

$$\widetilde{\stackrel{~}{\mathbf{u}}}=\tilde{\mathbf{u}}-\Delta t\nabla {\bar{p}}^{(n+1)},$$

$${\mathbf{u}}^{(n+1)}={\mathbf{u}}^{(n)}+\widetilde{\stackrel{~}{\mathbf{u}}},$$

$$\tilde{\mathbf{u}}_{AB}=\Delta t\sum_{j=0}^{J-1}\beta_{j}\big{(}-\mathbf{u}^{(n-j)}\cdot\nabla\mathbf{u}^{(n-j)}+\frac{1}{Re}\nabla\cdot(\boldsymbol{\sigma}^{(n-j)}/\rho)+\mathbf{g}^{(n-j)}/\rho\big{)},$$

$$\mathbf{S}=\left(\begin{array}{cc}\frac{\partial u}{\partial x}& S_{12}\\ S_{12}& \frac{\partial v}{\partial y}\end{array}\right),$$

$$\mathbf{S}=\left(\begin{array}{ccc}\frac{\partial u}{\partial x}& 0& S_{13}\\ 0& 0& 0\\ S_{13}& 0& \frac{\partial w}{\partial z}\end{array}\right).$$

$$\mathbf{\Theta }=\left(\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right),$$

$${\mathbf{S}}^{\prime }=\mathbf{\Theta }\mathbf{S}{\mathbf{\Theta }}^T=\left(\begin{array}{cc}0& \dot{\gamma }/2\\ \dot{\gamma }/2& 0\end{array}\right).$$

$${\mathbf{S}}^{\prime }=\mathbf{P}\mathbf{\Lambda }{\mathbf{P}}^T=\mathbf{P}\left(\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right){\mathbf{P}}^T,$$

$${\mathbf{S}}^{\prime }=\frac{1}{2}\left(\begin{array}{cc}-1& 1\\ 1& 1\end{array}\right)\left(\begin{array}{cc}-\dot{\gamma }/2& 0\\ 0& \dot{\gamma }/2\end{array}\right)\left(\begin{array}{cc}-1& 1\\ 1& 1\end{array}\right)\left[=\left(\begin{array}{cc}0& \dot{\gamma }/2\\ \dot{\gamma }/2& 0\end{array}\right)\right]$$

Full text

††journal: Computer Physics Communications

1 Introduction

Modeling and computational simulation of soft matter liquids remains a challenging problem because these fluids may exhibit complex non-Newtonian effects, such as shear-thinning/thickening, viscoelasticity or flow-induced phase transition. Such complex behavior is attributed to microstructure changes in fluids when a system is subject to an external mechanical shear force [28, 29]. Therefore, computational modeling of soft-matter fluids has to necessarily take into account microscopic effects in order to obtain reliable numerical solutions.

Clearly, the most accurate description of soft-matter fluids can be obtained by the molecular dynamics (MD). However, such microscale description is computationally inefficient, if large scale regions in space and time need to be simulated. To overcome this restriction and to obtain practically tractable simulation techniques hybrid molecular-continuum methods have been proposed in the literature aiming in combining the best attributes of both parts: the molecular accuracy with the computational efficiency of continuum models.

Bridging the large range of dynamically coupled scales is a fundamental challenge that is a driving force in the development of new mathematical algorithms. In general, hybrid models can be divided in two groups: based on the Eulerian-Lagrangian decomposition or on domain decomposition. In the first type the Lagrangian-type particles are embedded in the Eulerian fluid description, see, e.g., [14, 34, 42]. The second type of the methods is based on the domain decomposition into a small accurate atomistic region embedded into a coarser macrosopic model, see, e.g., [15]. In the literature we can find several hybrid models combining particle dynamics with the macroscopic continuum model, see, e.g., the hybrid heterogeneous multiscale methods described in [8, 6, 9, 11, 33, 34, 42], the triple-decker atomistic-mesoscopic-continuum method [15], the seamless multiscale methods [7, 10], the equation-free multiscale methods [22, 23] or the internal-flow multiscale method [2, 3]. In [24] a overview of multiscale flow simulations using particles is presented. The essential question that arises in building a coupled multiscale method is how micro- and macroscopic models are linked together, i.e., how projection/lift (or compression/reconstruction) operators are defined and implemented.

On the artificial boundary of the particle domain embedded into the macroscopic domain following typical strategies of constraint dynamics can be found in the literature: the Maxwell buffer [20], the relaxation dynamics [31], the least constraint dynamics [30], and the flux imposition [17]. Truncation of the microscopic domain is realized by imposing suitable boundary conditions. Non-periodic boundary conditions involve particle insertions and deletions, special wall reflections and body force terms, see [15] for more details. The deformation of the boxes mimics the time evolution of the control volume element in continuum and requires an adaption of standard Lees-Edwards periodic boundary conditions [40].

In order to extract mean flow field information from particle-based simulations averaging needs to be performed after a specified number of time steps. For example, the required rheological information for the stress tensor is calculated using the Irving-Kirkwood expression [21] and passed to the macroscopic continuum model.

In fact all of these techniques can be considered as hybrid particle-continuum methods under the statistical influence of microscale effects since coefficients in coarse-grained equations are estimated from data that are obtained from microscale simulations. As demonstrated in [3] the sensitivity of the accuracy of a solution, as well as the computational speed-up over a full molecular simulation, is dependent on the degree of scale separation that exists in a problem. For the case when processes occurring on a small scale are only loosely coupled with the behavior on a much larger scale and the so-called scale separation in the flow direction occurs, the hybrid multiscale schemes can be successfully applied, see [2, 3, 8, 9, 14, 33, 34, 41, 42, 43] and the references therein.

A complementary approach, the fluctuating hydrodynamics goes beyond the mean flow field of the hybrid simulation. In this case, the statistical influence of microscale effects is explicitly taken into account in the macroscopic flow equations leading to the stochastic partial differential models, such as the Landau-Lifshitz-Navier-Stokes system [4, 5, 19]. We refer reader to a recent study [1], where the errors in the fluctuations, due to both the truncation of the domain and the constraint dynamics performed on the artificial boundary are analysed for hybrid shear flow simulations.

In our recent paper [14] we have developed a novel hybrid multiscale method that is based on the combination of the discontinuous Galerkin (dG) method and molecular dynamics (MD) in order to simulate complex fluids, such as colloids in a Newtonian solvent. It has been shown that the method can be applied successfully to complex fluids when scale separation occurs and we can assume that the statistical influence of the microscale can be controlled on the macroscale. Our dG-MD hybrid method combines the following advantages (i) for macroscopic flow equations the dG method is applied which allows more flexible discretization including per-cell momentum conservation, (ii) the reduced order techniques are included in order to control the number of needed but computationally expensive MD simulations. The main goal of the present paper is to focus on the molecular dynamics part, which typically poses the bottleneck in these hybrid molecular-continuum approaches. We will discuss strategies, which minimize the computational effort in the particle-based simulations and discuss optimum choices for thermostats, time steps and relate time and length scales from simulations to experiments. Moreover we investigate the coupling of the microscopic simulation data to the macroscopic flow solver and propose a novel reduced-order strategy based on the combination of the proper orthogonal decomposition, the regularized least-square approximation and a suitable greedy algorithm to approximate the unknown nonlinear stress-strain function efficiently. This is the first time that the reduced-order technique is used in the context of hybrid simulation methods. As a test case we investigate Couette and Poiseuille flow in two and three dimensions.

2 Microscale (particle-based) simulations

In non-equilibrium Molecular Dynamics (MD) [18], we simulate colloidal particles with a coarse-grained model. Since we simulate a one dimensional flow, standard Lees-Edwards periodic boundary conditions [25] are applied. In our case, we shear in $z$-direction, and apply the velocity-Verlet algorithm to solve the equations of motion [38].

The stress-tensor is calculated via the Irving-Kirkwood formula [21], using the peculiar velocities of particles $\mathbf{v}_{i}=\tilde{\mathbf{v}}_{i}-\mathbf{v}_{i,S}$, where $=\tilde{\mathbf{v}}_{i}$ is the total and $\mathbf{v}_{i,S}$ the streaming velocity of particle $i$, respectively.

[TABLE]

where $\mathbf{r}_{ij}$ and $F_{ij}$ are the distance and the force between particle $i$ and $j$. The pressure corresponds to

[TABLE]

The dynamic viscosity is

[TABLE]

where $\dot{\gamma}$ is the shear rate. Note that $\sigma_{13}$ will be used later on to couple microscopic simulations to the macroscopic solver.

Colloids are treated as hard spheres. The interaction of two particles is simulated with a Weeks-Chandler-Andersen potential (WCA), which corresponds to the repulsive part of the Lennard-Jones potential:

[TABLE]

Two particles in the interaction radius $r_{C}$ reject each other. Outside this radius $r_{C}$, the potential is [math]. In formula (4) $\varepsilon$ corresponds to the well depth of the Lennard-Jones potential. In our simulations, we set $\varepsilon=1k_{B}T$. $\sigma$ is the typical diameter of a colloid and is used as the length scale of the MD simulation.

It is worth noting that if particles represent colloids, $\sigma\approx 1~{}\mu m$ and typical time scales $t_{MD}$ are in the order of $10^{-4}~{}s$, which is close to experimental relaxation times. If particles represent atoms instead ($\sigma\approx 10^{-10}~{}m$), the gap between time scales of the simulation and experimental (macroscopic) time scale spans many orders of magnitude ($t_{MD}\approx 10^{-12}s$). See Appendix A.

Shearing a system leads to friction, which heats up the sheared system. Since we want to simulate in the NVT-ensemble, we need to cool the system down to an initial temperature $T_{0}$. In the following we discuss two thermostats. The Lowe-Andersen thermostat (LAT) [27] conserves momentum and is Galilean invariant. In this thermostat, after solving the equations of motion via the velocity-Verlet algorithm, particles undergo a ”bath” collision with a probability of $\Gamma\Delta t$ (with $0<\Gamma\Delta t\leq 1$).

Two particles within an interaction range $r_{ij}<r_{C}$ have the possibility to get a ”kick” along the center of mass, where the ”kick” is taken from a Maxwell-Boltzmann distribution.

[TABLE]

where $\xi_{ij}$ is a random number with unit variance and $\hat{r}_{ij}$ the normalized particle distance. The new velocities for particles $i$ and $j$ after a bath collision are:

[TABLE]

and

[TABLE]

where

[TABLE]

which leads to momentum conservation. In the following we set the bath collision probability to $\Gamma\cdot\Delta t=0.001$.

The second thermostat we would like to discuss is the simple isokinetic (ISO) thermostat [42], a popular choice amongst hybrid simulation schemes. After each velocity-Verlet integration step, we compute the temperature of a system of $N$ particles:

[TABLE]

and rescale the velocities with the factor

[TABLE]

where $T_{0}$ is the temperature at which we want to simulate. (Here, $T_{0}=1\frac{\varepsilon}{k_{B}}$.)

Additionally, we implement the SLLOD algorithm [42], which applies the shear profile to the equations of motion. Like this, we reach the steady state even faster. The equations of motion change to

[TABLE]

[TABLE]

where $\bigtriangledown\mathbf{u}$ is the matrix representation of the applied stress:

[TABLE]

since we shear the $xy$-plane in the $z$-direction.

Here, we extend the equations of motion for the half step, e.g., one equation for the generalized coordinates $\mathbf{q}^{n+\frac{1}{2}}$ due to the new terms of the SLLOD method in (11) and (12). These equations will then be needed for the full step. For a given $\Delta\mathbf{u}$, the half steps look like:

[TABLE]

where $n$ is the current time step, $n+\frac{1}{2}$ is the time of the half step and $n+1$ is the full step. The index $x$ in $\mathbf{q}$ and $\mathbf{p}$ correspond to the $x$ coordinate of the vector.

Analogously to the half step, we have to modify the equations of motion for $\mathbf{q}$, which includes the half step $\mathbf{q}^{n+\frac{1}{2}}$, consistent with the time level of $\mathbf{p}^{n+\frac{1}{2}}$:

[TABLE]

This method takes advantage of the Crank-Nicolson method. As pointed out in [26] it allows larger time steps and yields one order smaller errors.

Figure 2 shows averaged shear profiles obtained with (a) the Lowe-Andersen thermostat and (b) the isokinetic thermostat with SLLOD in agreement with the shear rate imposed by Lees-Ewards boundary and SLLOD conditions. Note that SLLOD conditions in the isokinetic case impose a linear shear profile, which takes somewhat longer to emerge for the Lowe-Andersen thermostat, and increases further for lower shear rates. Figure 2c displays a typical relaxation of the off-diagonal component of the stress tensor $\sigma_{13}$ after turning on shear. The component relaxes in less than 200 MD times as indicated by the running mean (red curve). We have verified that this holds for various shear rates and densities for the small system sizes considered in this study, and in all simulations discussed from now on, measurements start after 200 MD times to omit influences of the relaxation process.

In Figure 3 we check to which extent the Lowe-Andersen thermostat is able to thermalize the system for our given set of simulation parameters ($\Gamma\cdot\Delta t=0.001$). For high densities and shear rates, deviations in the temperature (a), trace of the stress tensor (b) and viscosity (c) are noticeable in comparison to the isokinetic thermostat, which strictly enforces temperature. As already stated in the original publication on the Lowe-Andersen thermostat [27], the viscosity of such a system is somewhat elevated (Figure 3d), and the surplus viscosity is only expected to vanish in the limit of small $\Gamma$ and large time steps, which again counteracts efforts to keep temperature and pressure fixed. On the other hand, this enables exploration of fluids with somewhat larger viscosities in the context of hybrid-continuum schemes.

In such a scenario, the model needs to be considered in conjunction with the thermostat and thermostat parameters. (Note that even though viscosity is elevated, the system still behaves like a Newtonian fluid, i.e., the viscosity changes only little as a function of the shear rate (Figure 3c).) If we are, however, interested in the viscosity of the actual model system or related properties, it is preferable to apply a thermostat, which does not enhance the latter (such as Lowe-Andersen or DPD). Even though the isokinetic thermostat with SLLOD (which we will use from now on) is somewhat unphysical, we are not interested in realistic dynamics, but properties of the steady state.

In attempt to minimize the computational effort, we investigate the dependence of the time step on pressure (Figure 4a) and viscosity (Figure 4b). While time steps larger than $10^{-3}$ (for the isokinetic thermostat) yield noticeable deviations in pressure and viscosity for dense systems and high shear rates, differences for a time step of $10^{-3}$ are probably acceptable in the context of a hybrid scheme. In this paper we use $10^{-4}$ if not noted differently.

When performing microscopic simulations in the context of a hybrid scheme, the question arises whether it is better to simulate a single large system or multiple systems of smaller size to gather statistics. From a computational point of view, the latter is typically to be preferred as MD simulations in practice often do not scale perfectly linear with the number of particles. Therefore, computational resources required to compute the contribution of ,e.g., a single particle to the stress tensor (over a given simulation period) are larger in a single large system when compared to multiple small ones. In addition, larger systems also require longer relaxation time. In Figure 4c and d, we test the lower limits of system sizes for our model by decreasing the size of the simulation box in the direction perpendicular to the applied shear. As seen in Figure 4c and d, finite-size effects start to play a role if the height of the box is smaller than four $\sigma$, which corresponds to four particle diameters. Of course, this value will change if the potential is longer-ranged or if particles become correlated, e.g., in the vicinity of a critical point.

Finally, Figure 5 shows the viscosity (inset) and the off-diagonal component of the stress-tensor $\sigma_{13}$ from which the viscosity is derived as a function of shear rate. As already indicated for large shear rates (Figure 3), the viscosity changes only little over many orders of magnitude and behaves like a Newtonian fluid. It is worth noting, however, that the computational effort to obtain meaningful values for the viscosity increases manifold when the shear rate is reduced. While fluctuations in $\sigma_{13}$ are comparable across shear rates (not shown), fluctuations in the viscosity $\eta=\sigma_{13}/\dot{\gamma}$ increase by an order of magnitude if the shear is reduced accordingly. This increase in fluctuations translates into an increase of computational effort by two orders of magnitude if we want to keep errors at the same level.

In Section 3, we will briefly outline the macroscopic simulation. In Section 4 we explain in detail the coupling of the micro and the macro level for which we need the data from Figure 5 as input.

3 Macroscale simulations

At the macroscopic level, the motion of the incompressible fluid flow is governed by the continuity and momentum equations

[TABLE]

where $\mathbf{u}$ is the velocity vector, $\boldsymbol{\sigma}$ the Cauchy stress tensor, $\mathbf{g}$ an external body force, $\rho$ the density, which is constant, and $Re$ the Reynolds number. The computational domain $\Omega$ is surrounded by the boundary $\partial\Omega=\partial\Omega_{D}\cup\partial\Omega_{P}$, where the Dirichlet and periodic boundaries are considered, respectively. In case of the Navier-Stokes equations for Newtonian fluids, $\boldsymbol{\sigma}=-p\textbf{I}+\boldsymbol{\tau}$, where $\boldsymbol{\tau}=\mu(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})$, the above momentum equations reduce to

[TABLE]

For the time integration of the continuity equation (19a) and momentum equations (19b) we apply the following multi-step projection method [13]. Using the first-order Euler method we have

[TABLE]

Here we define the average pressure $\bar{p}^{(n+1)}=1/\Delta t\int_{t^{n}}^{t^{n+1}}p^{\prime}/\rho\,dt$ with the normal derivative $\partial\bar{p}^{(n+1)}/\partial\mathbf{n}=\mathbf{n}\cdot\nabla\bar{p}^{(n+1)}$. To obtain a unique solution for (21b)-(21c), we require $\int_{\Omega}\bar{p}^{(n+1)}=0$. Note that $p^{\prime}$ is a correction to the pressure to ensure the divergence-free constraint. As one notices, the pressure is already present in the stress tensor in equation (19b).

For the second-order time integration of the velocity, we use the Adams-Bashforth method in the first step \Romannum1,

[TABLE]

with the coefficients $\beta_{0}=3/2$ and $\beta_{1}=-1/2$. However, the effective pressure, $\bar{p}^{(n+1)}$, is first-order accurate in time. If required, one can reconstruct the pressure for higher-order accuracy, see [32, 35]. In the above projection scheme, we use by construction $\nabla\cdot\mathbf{u}^{(n-j)}=0$, for all $j\geq 0$. Therefore, we can integrate the unsteady terms in the third step \Romannum3. In this way, by replacing the intermediate velocity from the first step, the right-hand side of the Poisson equation in the second step \Romannum2 is independent of the time step. This prevents the numerical instability observed in [12, 16, 37] using the dG method.

4 Hybrid multiscale method (HMM)

The particle simulations represent a bottle-neck of our hybrid method. To compute the solution of a macroscopic problem, the stress tensor is required at each quadrature point, as shown for example in Figure 6 for a two-dimensional fluid dynamics problem. In order to reduce the number of these time-consuming simulations, we do not follow the strategy of the one-to-one correspondence between the MD boxes and the quadrature points on the macroscopic level, but we split our simulations into an off-line and on-line phases to prepare approximations for the stress in advance. In an off-line training phase, hence, we collect most of the information using the greedy method (see, e.g., [39]) and build functional dependencies of the stress on shear rates based on the data approximation with Chebyshev’s orthogonal polynomials. In an on-line phase of fast multiple queries we obtain the required data using these approximations and in rare cases rebuild them if new strain rates appear which exceed the approximated intervals.

The stress data is provided by the MD method using one-dimensional shear flow (the so-called simple-shear flow) realized along one of the coordinate axes as described above. However, the macroscopic model represents a fluid flow with a two- or three-dimensional velocity field. In order to obtain one-dimensional shear rates for the particle simulations we have to rotate the MD boxes along the streamlines and change correspondingly the data basis [41].

Furthermore, the stress data provided by the MD simulations is subject to statistical and systematic errors. To reduce the noise in MD data we apply the method of Proper Orthogonal Decomposition and project our simulation data onto principal component directions. In what follows we describe our approach in more details.

4.1 Reduced-order approach for data refinement strategy

In order to find a proper approximation of the stress tensor with fewer number of particle simulations, we solve the optimization problem with a relatively small number of samples. Then we use the greedy algorithm (worst scenario search), see e.g. [39], for the data refinement to suggest the shear rate(s) for new particle simulations. If one plots the residual $(Ax-b)$ versus $\dot{\gamma}$, the proposed shear rate for a new simulation $\dot{\gamma}_{new}=0.5(\dot{\gamma}_{M}+\dot{\gamma}_{N})$, is found in the neighborhood of $\dot{\gamma}_{M}=\operatorname*{arg\,max}_{\dot{\gamma}}(Ax-b)$, where $\dot{\gamma}_{N}$ is the left or right neighbor of $\dot{\gamma}_{M}$ which corresponds to the larger residual.

4.2 Eigenvalue decomposition of strain and stress fields

Any strain rate tensor can be written as a symmetric matrix of streaming velocity gradients $\boldsymbol{S}=(1/2)(\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^{T})$. In the case of two-dimensional simulations the strain rate matrix is

[TABLE]

where $S_{12}=\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)/2$. In the case of plane flow in three-dimensional simulations the strain rate matrix is

[TABLE]

where $S_{13}=\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)/2$.

In the following part, for simplicity we focus on the derivation of the basis transform for the two-dimensional case. The plane flow in the three-dimensional case follows analogously and in the next section we also compare flow profiles of 2D Couette and 3D Poiseuille flows with analytic solutions as proof of concept.

As pointed out in [40] there exists an angle $\theta=\frac{1}{2}\arctan(-\frac{S_{11}}{S_{12}})$, where $S_{ij}$ are components of the strain rate matrix $\boldsymbol{S}$ in (21v), and a rotation matrix

[TABLE]

which transforms the strain rate tensor $\boldsymbol{S}$ to the anti-diagonal matrix

[TABLE]

This strain rate tensor $\boldsymbol{S}^{\prime}$ corresponds to a pure-shear deformation (i.e., in absence of normal stresses) with the shear rate $\dot{\gamma}/2=-S_{11}\sin(2\theta)+S_{12}\cos(2\theta)$. Eigenvalue decomposition of the pure-shear strain rate matrix $\boldsymbol{S}^{\prime}$ yields

[TABLE]

where $\lambda_{i}\in\mathbb{R}$ are the eigenvalues of $\boldsymbol{S}^{\prime}$, $\boldsymbol{P}=(\boldsymbol{p}_{1},\boldsymbol{p}_{2})$ is the matrix of the corresponding eigenvectors. By solving the eigenvalue equation $\boldsymbol{S}^{\prime}\boldsymbol{P}=\lambda\boldsymbol{P}$ one can find that $\lambda_{1,2}=\mp{\dot{\gamma}}/2$ and from $\boldsymbol{S}^{\prime}\boldsymbol{p}_{i}=\lambda\boldsymbol{p}_{i}$ it follows that $\boldsymbol{P}=\frac{1}{\sqrt{2}}\begin{pmatrix}-1&1\\ 1&1\end{pmatrix}.$

Thus, the eigenvalue decomposition of $\boldsymbol{S}^{\prime}$ reads

[TABLE]

and is same as (21y).

We are now looking for a fitting function $f(\dot{\gamma})$ approximating the MD stress data $\boldsymbol{\sigma}^{\rm MD}\equiv\{\sigma_{\alpha\beta}(\dot{\gamma})\}_{N_{\dot{\gamma}}}^{N_{\rm sets}}$ for $N_{\dot{\gamma}}$ values of the shear rate ${\dot{\gamma}}$ and ${N_{\rm sets}}$ independent data sets

[TABLE]

Thus, $f(\boldsymbol{S}^{\prime})$ should approximate $\boldsymbol{\sigma}^{\rm MD}$

[TABLE]

Substituting the eigenvalues of $\lambda_{1,2}$ into (21ac) yields

[TABLE]

This means that in order to obtain a least-square approximation $f(\dot{\gamma})$ we need to transform MD data into the eigenvector basis $\boldsymbol{P}$

[TABLE]

We can therefore assume $\tilde{\sigma}_{12}=\tilde{\sigma}_{21}=0$ and fit

[TABLE]

[TABLE]

The stress data, $\boldsymbol{\sigma}^{\rm MD}$, were obtained for various densities using the SLLOD method with a WCA potential in two- and three-dimensional MD simulations, as discussed in Section 2. In Figure 7 we show 2D MD stress data at $\rho=0.8$ and in Figure 8 the results of the orthogonal transform of this data into the eigenvector basis representation, $\tilde{\boldsymbol{\sigma}}$, via (21ae). The diagonal principal component functions $\tilde{\boldsymbol{\sigma}}_{11}(\lambda_{1})$ and $\tilde{\boldsymbol{\sigma}}_{22}(\lambda_{2})$ look very simple: $\tilde{\boldsymbol{\sigma}}_{11}$ for $\lambda_{1}\leq 0$ is nearly a straight line, $\tilde{\boldsymbol{\sigma}}_{22}$ for $\lambda_{2}\geq 0$ increases monotonously with a negative curvature, which changes with the density. These features make it easier to approximate $\tilde{\boldsymbol{\sigma}}_{11}$ and $\tilde{\boldsymbol{\sigma}}_{22}$ by low order polynomials $f(\lambda_{1})$ and $f(\lambda_{2})$, respectively. We also consider an approximation of the principal stress by $f(\lambda)$ with $\lambda=\lambda_{1}\cup\lambda_{2}$.

Note that the off-diagonal components of the matrix at the left-hand-side of (21ae) are zero. However, due to statistical errors of MD simulations (and perhaps due to a systematic error of the method) the off-diagonal components $\tilde{\sigma}_{12}$ and $\tilde{\sigma}_{21}$ of the matrix $\tilde{\boldsymbol{\sigma}}$ deviate from zero. These stochastic and systematic deviations increase with the shear rate, decrease with the density and lay in the range of few percents if compared to the values of the diagonal principal stress components. As mentioned above, we suppress these deviations by setting $\tilde{\sigma}_{12}=\tilde{\sigma}_{21}=0$ in the following processing of our data.

4.3 Proper Orthogonal Decomposition for noise reduction

In order to extract physically significant information from our MD data and to reduce the noise we apply the Proper Orthogonal Decomposition (POD) method, which is based on a Singular Value Decomposition (SVD) of a matrix. Hence, the matrix of principal stress can be represented as $\tilde{\boldsymbol{\sigma}}=\boldsymbol{V}\boldsymbol{\Sigma}\boldsymbol{Z}^{T}$, where $\boldsymbol{\Sigma}$ is a diagonal $m\times n$ matrix of rank $r\leq\min(m,n)$ with positive singular values $s_{1}\geq s_{2}\geq s_{r}>0$ stored along the main diagonal

[TABLE]

$\boldsymbol{V}$ is a $m\times m$ unitary matrix of left-singular vectors, $\boldsymbol{Z}$ is a $n\times n$ unitary matrix of right-singular vectors. $m$ is the number of $\dot{\gamma}_{i}$’s, $n$ is the number of independent simulations sampling $\boldsymbol{\sigma}^{\rm MD}(\dot{\gamma}_{i})$. The singular values $s_{i}$ are related to the eigenvalues of the correlation matrix $\boldsymbol{C}=\tilde{\boldsymbol{\sigma}}^{T}\tilde{\boldsymbol{\sigma}}$ as $s_{i}^{2}=\lambda_{i}^{(\boldsymbol{C})}$. According to the POD method one can reduce the rank of the original stress matrix $\tilde{\boldsymbol{\sigma}}$ by a low-rank approximation $\tilde{\boldsymbol{\sigma}}_{(k)}$ by keeping the first $k$ of $r$ singular values $s_{i}$ in $\boldsymbol{\Sigma}$

[TABLE]

where $k$ is a reduced rank, $\boldsymbol{\Sigma}_{(k)}$ is a reduced-rank singular-value matrix.

The SVD analysis of our principal stresses has shown that there exists one singular value, $s_{1}$, which is by two orders of magnitude greater than any other $s_{i}$, as shown in Figure 9 for 10 independent 2D MD data sets. The remaining singular values $s_{i}$ with $(i>1)$ have nearly same amplitude and correspond to the statistical noise of the MD data. Therefore, the first principal component with the largest variance along the principal direction given by the vector $\boldsymbol{V}_{1}$ approximates our principal stress data

[TABLE]

with the projection relative error

[TABLE]

For instance, in the case of singular values shown in Figure 9, the projection relative error of the rank-1 approximation, $\tilde{\boldsymbol{\sigma}}_{(1)}$, is $\epsilon_{(1)}\lessapprox 0.9\%$.

4.4 Chebyshev’s approximation with Tikhonov’s regularisation

The noise-reduced data of the principal stress, $\tilde{\boldsymbol{\sigma}}_{(1)}$, will be now approximated using the least squares method with the orthogonal polynomials of Chebyshev

[TABLE]

Here $k$ is the degree of the approximating polynomial, $a_{i}$ coefficients of the Chebyshev polynomials $T_{i}$. Our goal is to approximate simulation data by means of a simple function for $f(\dot{\gamma})$ for further use at the level of the macroscopic solver. To reduce a large number of particle simulations we use an off-line training phase and an on-line phase of fast multiple queries. For this training, we solve a least-square problem with the Tikhonov regularization for each principal component of the stress tensor. Tikhonov’s regularization improves the approximation of a badly conditioned data matrix, i.e., when the problem is ill-posed. Thus, our goal is to find a vector $\boldsymbol{x}$ minimizing an extended residuum functional

[TABLE]

Here $\boldsymbol{b}$ is the vector of $n$ ($n\geq k+1$) data points of the corresponding component of the principal stress tensor, $\tilde{\boldsymbol{\sigma}}_{(1)}$, obtained as described in the previous section. Further, $\boldsymbol{x}$ is the vector of the unknown Chebyshev coefficients $a_{i}$. $\boldsymbol{A}$ is the Vandermonde matrix of the Chebyshev polynomial basis. The penalty term $\boldsymbol{D}\boldsymbol{x}$ is added to damp the oscillations in the derivative function $\partial\boldsymbol{A}/\partial\dot{\gamma}\,\boldsymbol{x}=\boldsymbol{A}\boldsymbol{D}\boldsymbol{x}$. Parameters $\alpha_{1}$ and $\alpha_{2}$ are the regularization parameters. The equivalent problem $(\boldsymbol{A}^{T}\boldsymbol{A}+\alpha_{1}^{2}\boldsymbol{I}+\alpha_{2}^{2}\boldsymbol{D}^{T}\boldsymbol{D})\boldsymbol{x}=\boldsymbol{A}^{T}\boldsymbol{b}$ is solved by the LU factorization using the LAPACK library 111http://www.netlib.org/lapack.

The Chebyshev approximation with the Tikhonov regularization is shown in Figure 10. As it has been mentioned, the dependence $f(\lambda)$ exhibits two distinct behaviour: it is nearly linear for $\lambda<0$ and is concave for $\lambda>0$ with a density dependent curvature. Therefore, we approximate $f(\lambda)$ either by one joint function in the range $\lambda\in\mathbb{R}$ or by two separate functions $f(\lambda_{1})$ and $f(\lambda_{2})$ with $\lambda_{i}$ in ranges $\lambda_{1}<0$ and $\lambda_{2}>0$. In both cases, the joint and split approximations are indistinguishable in the plot scales, however, the degree of the optimal polynomials for the joint approximation $f(\lambda)$ is higher than that for the split approximations $f(\lambda_{i})$. From the practical point of view a low degree polynomial approximation is preferred, since it can significantly improve the numerical efficiency and stability of the automated algorithms for data processing.

4.5 Back transform of Chebyshev’s approximants

On the macroscopic level, the stress data is requested from the MD simulations level at each quadrature point for a given strain rate, $\boldsymbol{S}$. Hence, the principal stress calculated by the Chebyshev approximations, $f(\lambda_{i})$ with $\lambda_{1}=-\dot{\gamma}/2$, $\lambda_{2}=\dot{\gamma}/2$, has to be transformed back to the basis space of the macroscopic solver, $\boldsymbol{\sigma}(\boldsymbol{S})$. Using (21y), (21ad) and the property of the Galilean invariance of the rotation $\boldsymbol{\Theta}$ we have

[TABLE]

The rotation matrix $\boldsymbol{\Theta}$ is problem specific. Hence, in Figures 11 we compare the MD stress data, $\sigma_{ij}^{\rm MD}(\dot{\gamma})$, with $\boldsymbol{\sigma}^{\rm CHEB}(\dot{\gamma})$, the back transform of $f(\lambda_{i})$ to the MD basis space obtained from (21ad). One can see that the shear stress, $\sigma_{12}^{\rm MD}(\dot{\gamma})$, is described very well. Some discrepancies can be observed for the normal stresses: $\sigma_{11}^{\rm MD}$ and $\sigma_{22}^{\rm MD}$ from MD exhibit slight normal stress difference at $\rho=0.8$, whereas $\boldsymbol{\sigma}^{\rm CHEB}_{11}$ and $\boldsymbol{\sigma}^{\rm CHEB}_{22}$ are identical (in Figure 11, only $\boldsymbol{\sigma}^{\rm CHEB}_{11}$ is shown). This is a consequence of the suppression in (21ae) of the off-diagonal components, $\tilde{\sigma}_{12}$ and $\tilde{\sigma}_{21}$. For $\rho=0.8$ these terms are nearly zero and, hence, these discrepancies can hardly be recognized. Furthermore, we have checked the descriptive quality of the back transformation of the joint approximation, $f(\lambda)$ with $\lambda=\lambda_{1}\cup\lambda_{2}$, which includes higher degree Chebyshev’s polynomials. It yields very similar results, as one could expect after comparison in Figure 10.

5 Reduced-order hybrid simulations

The macroscopic hybrid simulations of the Couette flow are performed on the domain $[-1,1]\times[0,1]$. The flow is periodic in the streamwise $x$-direction. The no-slip boundary condition is applied at the walls. At the lower wall $y=0$, the velocity is zero. At the upper wall $y=1$, the velocity in $x$-direction is equal to $U$ and the velocity in $y$-direction is zero. A grid of $3\times 3$ cells is employed. A polynomial degree $k=1$ is assigned in the dG method, $Re=1$. The velocity profiles presented in Figure 12 overlap with the analytical solution $u_{x}=yU$ for the equivalent Newtonian fluid with density $\rho=0.8$ and viscosity $\mu=\left.\sigma_{12}/\dot{\gamma}\right|_{\dot{\gamma}=1}=1.96481$ calculated from the 2D MD data.

Three-dimensional hybrid simulations of the Poiseuille flow are performed for $\rho=0.8$ using the Chebyshev approximations to the $3$D MD data, analogously to the 2D case described in the previous section. Figure 13 shows the MD data and the approximations for the 3D case. In Figure 14 velocity profiles in $z$-direction are plotted for different values of the pressure gradient parameter $P$ and coincide with the analytic solution $u_{z}=x(1-x)P$. We take the shear stress $\sigma_{13}$ to estimate the viscosity of the equivalent Newtonian fluid using relation $\mu=\left.\sigma_{13}/\dot{\gamma}\right|_{\dot{\gamma}=0.5}=1.61776$. The viscosities for different densities are plotted in Figure 15, cf. [14]. The computational domain is taken as $[0,1]\times[-1,1]\times[-1,1]$. The flow is periodic in $y$- and $z$-directions. The no-slip boundary condition is applied at the walls $x=0$ and $x=1$, where the velocity is zero. The pressure gradient $f_{z}=-2\mu\,P/Re$ is applied in the streamwise $z$-direction with $Re=1$. A grid of $3\times 3\times 3$ cells is employed. The polynomial degree is $k=2$ in the dG method.

6 Conclusion

In hybrid particle continuum schemes, molecular simulations typically pose a computational bottleneck.

In the first part of this paper we discuss in detail strategies to minimize the effort of particle-based simulations. We find that it is preferable to simulate small system sizes and determine limits of the latter, which are imposed by finite-size effects such as the transition to two-dimensional behavior.

As we are interested in properties of the steady-state and not necessarily dynamics, a simple non-physical thermostat such as the isokinetic thermostat with SLLOD boundary conditions is actually preferable while thermostats such as Lowe-Andersen allow for simulations at elevated viscosity. We also investigate boundaries for the size of time steps, and discuss the mapping of time scales between micro- and the macroscale.

In the second part we propose a novel reduced-order technique, which allows us to approximate the stress-strain relation in efficient way and thereby further reduce the computational effort of the microscopic simulations. Our approach is based on a delicate combination of the eigenvalue decomposition, the least square approximation using the Chebyshev polynomials with the Tikhonov regularization and the proper orthogonal decomposition. The latter is applied to reduce statistical noise from the MD simulation data. Our hybrid dG-MD method is tested by evaluating velocity profiles of Couette and Poiseuille flow. To our knowledge this is the first time, such reduced-order techniques are applied in the context of hybrid heterogeneous atomistic-continuum simulations.

Acknowledgments

The present work is supported by German Science Foundation (DFG) under the grant TRR 146 (project C 5). We would like to thank M. Oberlack, F. Kummer and B. Müller for providing their code BoSSS and T. Raasch for fruitful discussions on the topic. The authors gratefully acknowledge the computing time granted on the HPC cluster Mogon at Johannes Gutenberg-University Mainz.

Appendix A Mapping

We assume that our particles correspond to colloidal particles. The typical size of a colloid is 1 $\mu m$ (range 100 $nm$ to 10 $\mu m$). Thus, the unit of length in the MD simulations is $\sigma=10^{-6}~{}m$.

The time scale in the simulation is given by the so-called MD time $t_{MD}$, corresponding to $\sim 10^{4}$ time steps in the simulation (for a time step of $10^{-4}$).

[TABLE]

The mass of a polystyrene bead can be determined with the density of polystyrene $\rho=1.05~{}\frac{t}{m^{3}}$:

[TABLE]

$t_{MD}$ should correspond to the structural relaxation time of a bead, i.e., a bead takes roughly $t_{MD}$ to travel to a position over the distance $\sigma$.

An alternative approach to determine relaxation times takes the diffusion into account, which leads to:

[TABLE]

in corresponding experiments [36]. Therefore, for colloidal particles simulation time scales roughly correspond to experimental time scales.

To simulate Argon, the well depth $\varepsilon=1.65\cdot 10^{-21}~{}J$. The diameter of the Argon atoms is equal to $\sigma=3.4\cdot 10^{-10}~{}m$, the density is $\rho=1.784~{}kgm^{-3}$. Inserting in the formula gives:

[TABLE]