Molecular Simulation of Thermo-osmotic slip

TL;DR

This paper compares three molecular simulation methods to predict thermo-osmotic slip at solid-fluid interfaces, demonstrating their practical feasibility despite minor differences under various conditions.

Contribution

It introduces and evaluates three different molecular simulation techniques for thermo-osmotic slip prediction, highlighting their consistency and practical applicability.

Findings

Different simulation approaches show minimal differences in results.

Practical molecular simulations of thermo-osmotic slip are feasible.

The methods are robust across various physical conditions.

Abstract

Thermo-osmotic slip -- the flow induced by a thermal gradient along a surface -- is a well-known phenomenon, but curiously there is a lack of robust molecular-simulation techniques to predict its magnitude. Here, we compare three different molecular simulation techniques to compute the thermo-osmotic slip at a simple solid-fluid interface. Although we do not expect the different approaches to be in perfect agreement, we find that the differences are barely significant for a range of different physical conditions, suggesting that practical molecular simulations of thermo-osmotic slip are feasible.

Equations

$$v_s=-\frac{2}{\eta }{\int }_0^{\infty }\mathrm{d}zz\Delta h(z)\frac{\nabla T}{T},$$

$$VdP=\sum _{i=1}^n{N}_{i}d{\mu }_i+SdT.$$

$$\frac{\partial P}{\partial x}=\left(\sum _{i=1}^n{\rho }_i\frac{\partial {\mu }_i}{\partial T}+\frac{S}{V}\right)\frac{\partial T}{\partial x},$$

$${\left(\frac{\partial {\mu }_i}{\partial T}\right)}_P=-s_i^B$$

$$\frac{\partial P(z)}{\partial x}=\left(\sum _{i=1}^n{\rho }_i(z)\frac{\partial {\mu }_i}{\partial T}+\frac{S(z)}{V}\right)\frac{\partial T}{\partial x}.$$

$$\frac{\partial P(z)}{\partial x}=\left(-\sum _{i=1}^n{\rho }_i(z)s_i^B+\sum _{i=1}^n{\rho }_i(z)s_i(z)\right)\frac{\partial T}{\partial x}.$$

$$\frac{\partial P(z)}{\partial x}$$

$$=\left(\frac{\Delta h(z)}{T}\right)\frac{\partial T}{\partial x},$$

$$\eta \left(\frac{{\partial }^2{v}_x(z)}{\partial z^2}\right)=\left(\frac{\partial P(x,z)}{\partial x}\right).$$

$$v_s=-\frac{1}{\eta }{\int }_0^{\infty }\mathrm{d}zz\left(\frac{\Delta h(z)}{T}\right)\frac{\partial T}{\partial x}.$$

$$T\rho \dot{s}=-v_x\nabla P-Q_x\frac{\nabla T}{T},$$

$$v_x=-{\beta }_{11}\nabla P-{\beta }_{12}\frac{\nabla T}{T}$$

$$Q_x=-{\beta }_{21}\nabla P-{\beta }_{22}\frac{\nabla T}{T}.$$

$$Q_x=h^B{v}_x+\frac{1}{2d}{\int }_{-d}^{+d}\Delta h(z)v_x(z)\mathrm{d}z.$$

$${\beta }_{21}=-{\left(\frac{Q_x-h^B{v}_x}{\nabla P}\right)}_T.$$

$${\beta }_{21}=\frac{1}{\eta }{\int }_0^{\infty }\mathrm{d}zz\Delta h(z).$$

$${\beta }_{12}=-{\left(\frac{v_x}{(\nabla T/T)}\right)}_P.$$

$$f_x(z)=-\left(\frac{P_{xx}^{eq,T_2}(x,z)-P_{xx}^{eq,T_1}(x,z)}{T_2-T_1}\right)\left(\frac{\partial T}{\partial x}\right),$$

$$P_{xx}^{Vir}(z)=⟨\rho (z)⟩k_{B}T-\frac{1}{V(z)}⟨\frac{1}{2}\sum _i^{N(z)}\sum _{j\ne i}\frac{x_{ij}^2}{r_{ij}}{\varphi }^{\prime }(r_{ij})⟩.$$

$$\begin{array}{rl}{P}_{xx}^{IK}(z)& =⟨\rho (z)⟩k_{B}T\\ & -\frac{1}{2A}⟨\sum _i^N\sum _{j\ne i}\frac{x_{ij}^2}{r_{ij}}\frac{{\varphi }^{\prime }(r_{ij})}{|z_{ij}|}\Theta \left(\frac{z-z_i}{z_{ij}}\right)\Theta \left(\frac{z_j-z}{z_{ij}}\right)⟩.\end{array}$$

$$h(z)=u(z)+\frac{P_{xx}^{Vir}(z)}{\rho (z)},$$

$$f_x(z)=-\frac{\rho (z)(h(z)-h^B)}{T}\left(\frac{\partial T}{\partial x}\right).$$

$$V_{\mathrm{trunc}}(r)=\{\begin{array}{ll}4ϵ\left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right]-V(r_c)& r\le r_c\\ 0& r>r_c.\end{array}$$

Full text

Thermo-osmosis and thermophoresis are phenomena of great practical interest in the context of non-isothermal hydrodynamics Morthomas and Würger (2009); Ajdari and Bocquet (2006), non-equilibrium thermodynamics Dhont et al. (2007), thermophoresis Anderson (1989); Piazza and Parola (2008); Würger (2010), and the propulsion of active matter Golestanian et al. (2007). Thermo-osmosis is usually described phenomenologically as the induced slippage of fluid along an interface, due to an applied temperature gradient. Phoretic motion is driven by the interfacial stresses induced by a temperature gradient in a microscopic boundary region of thickness $\lambda$, where the properties of the solvent are influenced by the interactions with the surface (or interface) Anderson (1989); Piazza and Parola (2008); Würger (2010).

Clearly, it would be useful to be able to predict thermophoretic slip on the basis of a molecular description of the solid-liquid interface. However, in practice this is not simple because much of the existing theoretical framework is couched in terms that assume the validity of a local continuum theory (e.g. Debye-Hückel plus the (Navier-)Stokes equation). Yet, crucially, near an interface, a continuum description of the structure or dynamics of a liquid is not allowed. More ominously, the definition of the stress in a liquid is not unique. This non-uniqueness has no effect on the computed value of, say, the liquid-liquid surface tension Schofield and Henderson (1982), but it could affect the prediction of phoretic flows, where the local value of the stress gradient is what drives the flow. In this paper, we consider this problem and explore various ‘microscopic’ methods to predict thermo-osmotic slippage in a simple model system.

The ‘classical’ approach to predict thermophoretic slippage is based on Onsager’s reciprocity relations (see Ref. de Groot and Mazur (1962)). Derjaguin Derjaguin et al. (1987) used Onsager’s Linear Non-Equilibrium Thermodynamics (LNET) approach to derive an expression for thermo-osmotic slip. Derjaguin’s approach exploits the relation between the flow caused by a temperature gradient and the excess heat flux due to hydrodynamic flow, resulting in the following equation for thermo-osmotic slip:

[TABLE]

where $\Delta h(z)$ is the excess enthalpy density at a height $z$ above the surface and $\eta$ is the viscosity. The difficulty with this expression is that there is some ambiguity in the microscopic definition of the local excess enthalpy $\Delta h(z)$, a quantity that is also not easy to probe in experiments Anderson (1989). Moreover, Eq (1) assumes that the (Navier-)Stokes equation with constant viscosity holds very close to a surface. A correct, microscopic description would not make such a continuum assumption. For these reasons, molecular simulations should be the method of choice to predict and study thermo-osmotic flows.

To place the various descriptions of thermo-osmotic slip in a broader context, we first consider the classical thermodynamic approach to the problem, based on the assumption of Local Thermal Equilibrium (LTE). We note that neither temperature gradients nor, for that matter, chemical potential gradients can exert a net force on a fluid element in a bulk liquid. Mechanical forces in liquids can only be caused by body forces such as gravity, or by pressure gradients. If temperature gradients cause flow near a surface, it is only because a local pressure gradient is induced. To clarify this, we first consider the thermodynamics of the problem. Consider a temperature gradient $+x$ direction parallel to a hard wall; the $z$ coordinate measures distance perpendicular to the wall. Starting from the Gibbs-Duhem relation for an $n$-component mixture, we write

[TABLE]

The Gibbs-Duhem equation makes use of the fact that the system is homogeneous. A stratified system in equilibrium, is homogeneous in the directions parallel to the stratification, but not perpendicular to it. Hence, here and in what follows, the ‘pressure’ $P$ refers to a component of the pressure tensor parallel to the surface (e.g. $P_{xx}$). Dividing Eq (2) by $V$ and differentiating with respect to $x$ gives the following expression:

[TABLE]

where $\rho_{i}$ is the number density of species $i$. In the bulk, the pressure is equalised quickly and the fluid reaches hydrostatic equilibrium. Since the bulk pressure is constant, Eq (3) reduces to

[TABLE]

using $S=\sum_{i}N_{i}s_{i}$, where $s_{i}$ is the specific entropy of species $i$ and the superscript $B$ denotes a bulk quantity. We now write a similar expression for the pressure gradient at a position $z$ above the surface

[TABLE]

We assume that there are no gradients of $\mu_{i}$ and $T$ perpendicular to the surface. Using Eq (4), Eq (5) can be rewritten as

[TABLE]

Eq (6) can be simplified by noting that the expression in brackets is the difference between the specific entropy at position $z$ and the bulk specific entropy. Since $\mu_{i}$ and $T$ do not depend on $z$, $\mu_{i}=h_{i}-Ts_{i}$ can be used to rewrite Eq (6) as

[TABLE]

where $\Delta{h(z)}$ is the excess enthalpy density at a distance $z$ from the surface. To compute the flow velocity, we combine the expression for the pressure gradient with the linearised (Navier-)Stokes equation

[TABLE]

Assuming no surface slip, this equation can be integrated twice to yield

[TABLE]

This is equivalent to Derjaguin’s expression in Eq (1) in the case of a single surface. The advantage of the LTE approach is that Eq (7) demonstrates the relationship between the external temperature gradient and the local pressure gradient that induces thermo-osmotic slip. The magnitude of this coupling is determined exclusively by the excess surface enthalpy.

To relate our expression to Derjaguin’s LNET approach, consider the slit pore as depicted in Fig. 1. Following Derjaguin, a pressure and temperature gradient is maintained across the slit. Fluid flows in the $-x$ direction as depicted by the arrows. For a one-component fluid, the rate of entropy production can be written as

[TABLE]

where $v_{x}$ is the fluid velocity ($m/s$) and $Q_{x}$ is the heat flux ($J/(m^{2}\cdot s)$). Eq (11) implies the following phenomenological equations:

[TABLE]

Consider the isothermal heat flux across the pore in Fig. 1

[TABLE]

The first term on the right-hand side of Eq (14) is the advective term and equivalent to the bulk heat content of the fluid that is transported across the pore. The second term is the transport of excess heat in the boundary layers. The second term is what determines the strength of thermo-osmotic slip de Groot and Mazur (1962). Therefore, we can write

[TABLE]

$\beta_{21}$ defined here is conventionally known as the ‘mechano-caloric’ coefficient. Using Derjaguin’s approximation that the velocity profile is linear in the boundary layer $v_{x}(z)=-2dz\nabla{P}/\eta$ and considering only the bottom wall, $\beta_{21}$ can also be expressed as

[TABLE]

By considering the isobaric mass flux in Eq (12), we can write

[TABLE]

Substituting Eq (10) for $v_{x}$ in Eq (17) immediately shows $\beta_{12}=\beta_{21}$ as expected. This provides the link between our LTE and Derjaguin’s LNET approach.

The usual definition of the ‘slip’ velocity is the extrapolated velocity at the interface, where the fluid density approaches zero. For a thin boundary layer, the slip velocity is equal to the fluid velocity in the bulk just outside the boundary layer. $\beta_{12}$ in Eq (17) is the ‘thermo-osmosis coefficient.’

An alternative approach to the LTE route is to compute the thermo-osmotic slip, using a mechanical route, i.e. by computing the force on a volume element directly from the gradient of the microscopic stress. Such an approach has been used by Han Han (2005), but it could be problematic due to the non-uniqueness of the definition of the microscopic stress. We start with the relation between the stress gradient and $f_{x}(z)$ the force per unit volume on a fluid element at a distance $z$ from the surface. Rather than computing the stress gradient in a non-equilibrium simulation, we use the fact that $P_{xx}$ depends on $x$, only through $T$. Hence,

[TABLE]

where the superscript $eq$ denotes equilibrium calculations that are both carried out at the same bulk pressure. With this method, $\Delta{P}/\Delta{T}$ is determined, and for any $\partial T/\partial x$, $f_{x}(z)$ can be computed via Eq (18). The thermo-osmotic force per particle $f^{P}_{x}(z)=f_{x}(z)/\rho_{ave}(z)$ where $\rho_{ave}(z)=(\rho_{T_{1}}(z)+\rho_{T_{2}}(z))/2$. To compute the thermo-osmotic flow, we now carry out a second simulation at the average temperature $T_{ave}=(T_{1}+T_{2})/2$, where we apply the local body force $f_{x}(z)$. The resulting slip velocity and therefore $\beta_{12}$ can then be computed. Note that in this approach, we make no continuum assumptions.

The calculation as described above is complicated by the fact that the pressure tensor in an inhomogeneous fluid is not unique Hafskjold and Ikeshoji (2002); Rowlinson and Widom (2003). Irving and Kirkwood (IK) Irving and Kirkwood (1950) proposed an expression by integrating the total momentum flux acting across a virtual surface element. This approach gives the appropriate mechanical force balance normal to the interface. However, as argued by Schofield and Henderson Schofield and Henderson (1982), the definition of the pressure tensor is not unique since any term with a vanishing divergence can be added without changing the momentum flux. All common definitions do, however, yield the correct surface tension.

In a simulation, we need to know the thermo-osmotic force acting on atoms, as opposed to the force on the fictitious surface of a volume element. This would suggest that the atom-based virial (V) expression for pressure might be preferable.

In order to determine if the choice of the pressure affects the computed thermo-osmosis coefficient, we computed $P_{xx}$ in Eq (18) using both the V and IK expressions. The V pressure is given by Hansen and McDonald (1990)

[TABLE]

The IK pressure is computed using Walton et al. (1983)

[TABLE]

In addition to these ‘mechanical’ expressions for the thermo-osmotic slip, consider the right-hand side of Eq (7). We express the local specific enthalpy as

[TABLE]

where $u$ is the specific internal energy. In Eq (21), we have made explicit that the pressure that enters into the expression for the local enthalpy must be the component that is parallel to the surface, as argued below Eq (2). The body force on a fluid element at a height $z$ above the surface is then given by

[TABLE]

Eq (22) can be computed in a simulation thermostatted at $T_{ave}$ and applied as a body force in the same vein as Eq (18).

We compare the above calculations of the slip coefficient with the result for $\beta_{12}$ that follows from Derjaguin’s approach based on the Onsager reciprocity relations. In this case, a uniform pressure gradient represented by a body force is applied during a simulation at $T_{ave}$. $\beta_{21}$ is computed via Eq (15) (see Supplemental Material S2). The LTE approaches for computing $\beta_{12}$ and the ‘Derjaguin’ method for computing $\beta_{21}$ should be equivalent if the temperature and pressure gradients are small enough to ensure that the resulting response is linear.

All Molecular Dynamics simulations reported here were performed using the LAMMPS package Plimpton (1995). The simulation setup is depicted in Fig. 1. The system consists of a Lennard-Jones fluid of $N=2640$ atoms interacting with a solid atomic wall.

[TABLE]

Interaction parameters are set so that $\sigma_{\text{fluid-fluid}}=\sigma_{\text{wall-fluid}}=\sigma$ and $\epsilon_{\text{fluid-fluid}}=\epsilon$ with $r_{c}=4\sigma$. Two different wall-fluid interactions were considered: a less attractive Lennard-Jones where $\epsilon_{\text{wall-fluid}}=0.55\epsilon$ and a purely repulsive Weeks-Chandler-Andersen potential Weeks et al. (1971) such that $r_{c}=2^{1/6}\sigma$. Wall atoms are fixed via harmonic springs. All computed quantities are expressed in Lennard-Jones units.

NVT dynamics with a time step $\Delta{t}=0.001$ were run to equilibrate the system. This was accomplished using a Nosé-Hoover thermostat for $100,000$ MD steps. For an additional $100,000$ steps, the system was barostatted at P$\approx 0.122$ by applying a downward force to the top wall atoms.

The solid surface introduces local anisotropy in the pressure that vanishes in the bulk. The response of this anisotropy to a change in temperature drives thermo-osmotic slip. Using the pressure profiles (see Supplemental Material: Figs. S1(a,b)), $\Delta{P_{xx}(z)}/\Delta{T}$ was computed for the three temperatures shown in Figs. 2(a) and 2(b). Encouragingly, the choice of the pressure tensor makes no significant difference to the measured response.

At constant temperature, $T_{ave}$, the specific kinetic energy is uniform everywhere and therefore computed by dividing the total average kinetic energy by the number of atoms. For the same temperatures, the specific potential energy profiles were spatially averaged in $z$. Using the profiles of $P_{xx}$ (Fig. S1(a)) and density profiles (Fig. S1(c)), $\Delta{h(z)}/T$ was computed via Eqs (21) and (22) and shown in Fig. 2(c).

We note that the V and IK expressions (Figs. 2(a), 2(b)) and the LTE quantity Fig. 2(c) show similar qualitative behavior. The body force vanishes in the bulk as it should. In both cases, the profiles flatten and shift outward as the temperature is increased.

The body force per particle $f^{P}_{x}(z)$ can be computed by dividing the profiles in Figs. 2(a-c) by $\rho(z)$ (Fig. S1(c)) and multiplying by a sufficiently small gradient e.g. $\partial T/\partial x=0.0005$ for WCA wall-fluid interactions. To compute slip, non-equilibrium simulations were carried out by applying these forces to the equilibrated systems at the appropriate temperatures. To obtain reasonable statistics, forces were applied for $10^{8}$ steps until the fluid approached a steady velocity. The slip was then computed for an additional $2\times 10^{8}$ steps. Figs. 3(a,c) show calculations of the slip velocity. There appears to be reasonably good agreement between the three approaches and the choice of the pressure tensor appears to make no difference. The flow profile computed at $T=0.9$ (Fig. 3(b)) shows that for WCA wall interactions the velocity decreases monotonically indicating that the viscosity close to the surface is constant. For less attractive Lennard-Jones (Fig. 3(d)), the viscosity is clearly not constant showing significant departure from (Navier-)Stokes and Derjaguin’s result (Eq (1)).

To compare our LTE approaches and Derjaguin’s method (see Supplemental Material S2), $\beta_{12}$ was computed via Eq (17) using the slip calculations shown in Figs. 3(a,c). Fig. 4 shows a comparison of all three methods. For the range of temperatures $T\sim 0.8-1.0$, there appears to be reasonable agreement. There is some discrepancy between $\beta_{12}$ computed via Eq (22) and $\beta_{21}$ in Fig. 4(a). This may be due to noise in the force profile or fluctuations in the barostat. As expected, $\beta_{12}$ and $\beta_{21}$ for solely repulsive wall-fluid interactions are considerably larger than those for interactions with an attractive component. This is entirely due to the amount of slip at the surface in response to either a temperature or pressure gradient. Furthermore, both cases demonstrate an approximately linear dependence of the thermo-osmosis coefficient with respect to temperature.

In addition to the methods described above, we also attempted to compute $\beta_{21}$ using linear-response theory Luttinger (1964). However, no reliable results were obtained as the statistical noise overwhelmed the signal.

In summary, we have considered four different methods to compute thermo-osmotic slip on the basis of molecular simulations. The first approach is based on a computation of the thermally-induced stress gradient method, computed using equilibrium simulations and then represented as a body force in non-equilibrium simulations. We find no evidence that different choices for the pressure tensor lead to different results. In the second approach, we compute the excess enthalpy density near the wall and use a local-thermodynamics formalism to derive the body force acting on the fluid. These methods do not assume that macroscopic thermodynamics or hydrodynamics holds close to an interface. The final approach is based on Onsager’s reciprocal relations, which allow us to derive thermo-osmotic slip from the excess heat flux due to a pressure gradient.

Our results are encouraging and surprising: we find that all methods yield results for the thermo-osmotic slip that do not differ significantly. Hence, the choice of the method to compute thermo-osmotic slip seems to be a matter of taste or convenience.

We gratefully acknowledge numerous discussions with Lydéric Bocquet, Mike Cates, Patrick Warren, Ignacio Pagonabarraga and Benjamin Rotenberg. Additionally, we are grateful to Peter Wirnsberger for his help with the error analysis. RG gratefully acknowledges a PhD Grant from the Sackler Fund.