Comparing Theory and Simulation for Thermo-osmosis

TL;DR

This paper presents a numerical comparison of three methods to calculate thermo-osmotic slip in an ideal gas model, confirming the theoretical predictions and demonstrating consistency among the approaches.

Contribution

It introduces and compares three numerical methods for calculating thermo-osmotic slip, validating the theoretical model with simulation results.

Findings

Numerical results are mutually consistent.

Results agree with theoretical predictions.

Hydrodynamics and thermodynamics are locally valid.

Abstract

We report a numerical study of thermo-osmotic slip, i.e. the particle flux induced by a thermal gradient along a solid-fluid interface. To facilitate comparison with theory, we consider a model of an ideal but viscous gas. We compare three numerical routes to obtain the slip coefficient: 1. by using the Onsager reciprocity relations 2. by using the appropriate Green-Kubo relation 3. via the excess enthalpy. The numerical results are found to be mutually consistent, and to agree with the theoretical prediction based on the assumption that hydrodynamics and thermodynamics are locally valid.

Equations

$$\sigma =\sum _i{\mathcal{F}}_i{J}_i.$$

$$J_i=\sum _j{L}_{ij}{\mathcal{F}}_j,$$

$$L_{ij}={\frac{\partial J_i}{\partial {\mathcal{F}}_j}|}_{\text{eq}},$$

$$L_{ij}=L_{ji}.$$

$$L_{ij}=\frac{1}{k_B}{\int }_0^{\infty }dt{⟨J_i(0)J_j(t)⟩}_{\text{eq}},$$

$$\sigma =-\left(\overline{Q_x}-h_B\overline{v_x}\overline{\rho }\right)\frac{\nabla T}{T^2}-\overline{v_x}\frac{\nabla P}{T},$$

$$\overline{Q_x}-h_B\overline{v_x}\overline{\rho }$$

$$\overline{v_x}$$

$$\overline{v_x}=\frac{{\sum }_{i=1}^N{v}_{x;i}}{N}=\frac{{\int }_0^{L_y}dy\rho (y)v_x(y)}{L_y\overline{\rho }},$$

$$\overline{Q_x}$$

$$u(y)=\{\begin{array}{ll}-\Delta E& y<\Delta L\\ 0& \Delta L<y<L_y-\Delta L\\ -\Delta E& y<L_y-\Delta L\end{array},$$

$$\overline{v_x^3}=\overline{{\left(v_{x;\text{th}}+v_{x;\text{drift}}\right)}^3}\approx \frac{3}{2}{k}_{B}T\overline{v_x},\overline{v_y^2{v}_x}=\frac{k_{B}T\overline{v_x}}{2}.$$

$$\overline{Q_x}\approx \overline{\rho }\left(2k_{B}T\overline{v_x}+\overline{u(y)v_x}\right).$$

$$h_B=\frac{{\sum }_i\frac{m\left(v_{x;i}^2+v_{y;i}^2\right)}{2}+P{V}_{\text{B}}}{N_{\text{B}}}=2k_{B}T,$$

$$\overline{Q_x}-h_B\overline{v_x}\overline{\rho }=\frac{{\sum }_{i}u(y_i)v_{x;i}}{L_x{L}_y}=\frac{{\int }_0^{L_y}\rho (y)u(y)v_x(y)}{L_y},$$

$$\frac{d}{dy}\left(\eta (y)\frac{d}{dy}{v}_x(y)\right)=\nabla P(y),$$

$$\rho (x,y)$$

$${\rho }_{\text{W}}$$

$${\rho }_{\text{B}}$$

$$\nabla P(x,y)$$

$$v_x(y)=\{\begin{array}{ll}{v}_{x;\text{W}}(y)& y<\Delta L,\\ v_{x;\text{B}}(y)& \Delta L\le y\le L_y-\Delta L\\ v_{x;\text{W}}(L_y-y)& y>L_y-\Delta L\end{array}$$

$$v_{x;\text{W}}(y)$$

$$v_{x;\text{B}}(y)$$

$$\begin{gathered} \overline{v_x}=-\frac{\nabla P_{\text{B}}{\rho }_B}{12{\eta }_B{L}_y\overline{\rho }}{\left(L_y-2\Delta L\right)}^3-\frac{2\nabla P_{\text{W}}{\rho }_W\Delta L^3}{3{\eta }_W{L}_y\overline{\rho }}-\frac{{\rho }_{\text{B}}}{2{\eta }_W{L}_y\overline{\rho }}\left(L_y-2\Delta L\right)\left(\Delta L(L_y-2\Delta L)\nabla P_{\text{B}}+\Delta L^2\nabla P_{\text{W}}\right) \\ -\frac{{\rho }_W}{2{\eta }_W\overline{\rho }{L}_y}\Delta L^2(L_y-2\Delta L)\nabla P_{\text{B}} \end{gathered}$$

$$\begin{gathered} \overline{Q_x}-h_B\overline{v_x}\overline{\rho }=\frac{\Delta E\Delta L^2{\rho }_W}{6{\eta }_W{L}_y}\times \\ \left(3(L_y-2\Delta L)\nabla P_{\text{B}}+4\Delta L\nabla P_{\text{W}}\right) \end{gathered}$$

$$\nabla P_{\text{W}}$$

$$\nabla P_{\text{B}}$$

$$L_{12}=-\frac{T\overline{\rho }\Delta E{L}_y\Delta L^2{e}^{\beta \Delta E}\left(3(L_y-2\Delta L)+4\Delta L{e}^{\beta \Delta E}\right)}{6{\eta }_W{\left(L_y-2\Delta L+2\Delta L{e}^{\beta \Delta E}\right)}^2}.$$

$$\nabla P_{\text{B}}=k_B\left({\rho }_{\text{B}}\nabla T+T\nabla {\rho }_{\text{B}}\right)=0,$$

$$\nabla {\rho }_{\text{B}}=-\frac{{\rho }_{\text{B}}\nabla T}{T}.$$

$$\nabla {\rho }_{\text{W}}$$

$$\nabla P_{\text{W}}$$

Full text

thermo-osmotic slip, linear irreversible thermodynamics

††preprint: AIP/123-QED

I Introduction

With the increasing importance of nano-scale transport, both in man-made devices (nano-fluidics) and in cell biology, there is a great need to improve our “microscopic” understanding of phoretic transport, i.e. transport that occurs only in the presence of interfaces.

Whilst there exists an extensive literature on this subject Anderson (1989); Levich (1962); Derjaguin et al. (1987); De Groot and Mazur (2013); Würger (2010), the focus of earlier publications has been on local continuum descriptions. Subsequently, there have been several publications on mesoscopic Yang and Ripoll (2013); Lüsebrink, Yang, and Ripoll (2012) and atomistic Han (2005); Ganti, Liu, and Frenkel (2017, 2018) simulations of phoretic flows, but there is a shortage of simulations on model systems that are simple enough to allow quantitative comparison between analytical theory and simulation.

In this paper, we present such a study, by considering thermo-osmotic flow of a simple model gas, described by Multi-Particle Collision (MPC) dynamics (also known as Stochastic Rotation Dynamics (SRD)) Gompper et al. (2009). The MPC fluid obeys the equation of state of an ideal gas, but it undergoes collisions that mediate the transport of particles, energy and momentum. The transport properties of the MPC model are (to an excellent approximation) known analytically Gompper et al. (2009). As a consequence, we can not only compute the thermo-osmotic slip coefficients for this fluid numerically but we can also predict them analytically.

We stress that the interest of the present paper does not lie in the computed numbers as such (after all, the model that we used was chosen for simplicity, not for realism), but simply in the fact that all different methods that we used to compute thermo-osmotic slip appear to be consistent. Hence, future authors may select any of the methods used, based on considerations of convenience.

In section II, we give a brief overview of the relevant framework of linear irreversible thermodynamics. This framework will then be used in section III to obtain an expression for the thermo-osmotic slip coefficient of our model system. In section IV, we will compare these theoretical predictions with numerical numerical results.

II Linear irreversible thermodynamics

In reviewing the framework of linear irreversible thermodynamics in the context of thermo-osmotic slip, we focus on those aspects that are most directly relevant (we will, for instance, ignore magnetic fields and rotating systems). For a more thorough introduction, we refer to the standard text by De Groot and Mazur De Groot and Mazur (2013).

A system can be driven out of equilibrium by a thermodynamic force, $\mathcal{F}_{i}$. This driving force could be a temperature gradient, a chemical-potential gradient or a mechanical force. Every thermodynamic force induces a thermodynamic flux, $J_{i}$, such as a heat flux or a particle flux. The entropy production rate, $\sigma$, is then given by

[TABLE]

The central assumption of linear irreversible thermodynamics is that the thermodynamic forces are small enough to justify a first-order Taylor expansion of the thermodynamic fluxes in terms of the thermodynamic forces,

[TABLE]

with

[TABLE]

the Onsager coefficients. Here, the subscript eq stands for the equilibrium condition, $\mathcal{F}_{i}=0$ for all $i$. The above expansion is typically valid if the relation between forces and fluxes is analytical. The Onsager coefficients obey an important symmetry, namely the Onsager reciprocal relations,

[TABLE]

This property has several important consequences. For example, it implies that one can determine an Onsager coefficient, even without directly measuring the associated flux.

The Onsager coefficients can also be determined from the equilibrium correlations of the system,

[TABLE]

where $k_{B}$ is the Boltzmann constant. This relation, known as the Green-Kubo relation, gives another method to determine the Onsager coefficients without directly applying a thermodynamic force.

III Model system

In this section we describe the MPC model that we used to study thermo-osmotic slip. This model was chosen for simplicity. In particular, our aim was to consider a model where simulations can be compared directly with theoretical predictions.

III.1 Model

The model system that we studied was a two-dimensional MPC gas confined between two parallel flat walls, separated by a distance $L_{y}$ (see Fig. 1.) We used periodic boundary conditions along the slit direction. In our simulations, we impose no-slip boundary conditions at the hard walls. Away from the no-slip planes, the walls interact with the gas particles via a square-well potential of width $\Delta L$ and depth $\Delta E$. In the spirit of the MPC model, we assume that the gas particles do not interact with each other directly, i.e. they behave like an ideal gas. However, the MPC particles may thermalize via collisions.

In our simulations we impose a thermal gradient $\nabla T$ and a body force inducing a pressure gradient $\nabla P$ along the slit direction (the $x$-axis). These gradients result in heat and particle fluxes, $Q_{x}$ and $v_{x}$. The entropy production rate can be written as Ganti, Liu, and Frenkel (2017)

[TABLE]

where $\sigma$ and $\overline{Q_{x}}$ are the entropy production and heat flux per unit of volume, $\overline{v_{x}}$ is the average velocity per particle, ${h}_{B}$ is the average enthalpy per particle in the bulk (i.e. away from the walls), $\overline{\rho}$ is the average density of the gas, $\nabla P$ is the pressure gradient along the $x$-direction, and $T$ is the temperature of the system. For small gradients the fluxes can be written as (Eq. (2)):

[TABLE]

The off-diagonal Onsager coefficient $L_{21}$ now determines the thermo-osmotic slip, as it gives the particle flux induced by a temperature gradient.

We can now determine explicit expressions for the fluxes $\overline{v_{x}}$ and $\overline{Q_{x}}-h_{B}\overline{v_{x}}\,\overline{\rho}$. The average velocity is by definition given by

[TABLE]

where $\rho(y)v_{x}(y)$ is defined as $\sum_{i}\delta(y_{i}-y)v_{x,i}$. The heat flux can be written as,

[TABLE]

with potential energy $u(y)$ given by

[TABLE]

cf. Fig. 1(b). This potential mimics the interaction of the particles with the wall via a potential well. To make further progress, we split the velocity in a thermal part, $v_{x/y;\textrm{th}}$ and a drift part $v_{x;\textrm{drift}}$ (only present in the x-direction), which we assume to be independent. Furthermore, we assume that $\overline{v_{x/y;\textrm{th}}}=\overline{v_{x/y;\textrm{th}}^{3}}=0$ and $\overline{v_{x;\textrm{drift}}^{2}}\ll\overline{v_{x/y;\textrm{th}}^{2}}=k_{B}T/(2m)$, where the last equality follows from the equipartition theorem. One can now write

[TABLE]

which leads to a simplification for the heat flux,

[TABLE]

Invoking the equipartition theorem once more leads to a simple formula for the enthalpy of particles in the bulk,

[TABLE]

where the sum is over all particles in the bulk. Therefore, the total flux is given by

[TABLE]

where $\rho(y)u(y)v_{x}(y)$ is defined as $\sum_{i}\delta(y_{i}-y)u(y_{i})v_{x,i}$.

III.2 Navier-Stokes equation

To determine more explicit forms for the fluxes, we first assume that the flow should satisfy a (linearized) Navier-Stokes equation,

[TABLE]

where $\eta(y)$ is the viscosity. Furthermore, we assume that the thermodynamic forces are small enough for a local equilibrium ansatz to hold:

[TABLE]

with

[TABLE]

It should be stressed that, although independent of $y$, $\rho_{\textrm{W}}$ and $\rho_{\textrm{B}}$ can vary along the $x$-axis, either via the density or via a temperature gradient. For the MCP fluid, the viscosity only depends on the density Gompper et al. (2009) and hence it should also be piecewise constant, $\eta(y)=\eta_{W}$ if $y<\Delta L$ or $y>L_{y}-\Delta L$ and $\eta(y)=\eta_{B}$ otherwise.

In the situations that we will consider, the pressure gradient is either induced by a constant body force or by a temperature gradient, therefore it will also be piecewise constant along the $y$-axis ,

[TABLE]

where the specific values of $\nabla P_{\textrm{B}}$ and $\nabla P_{\textrm{W}}$ depend on the driving forces.

Under the assumption of no-slip boundary condition and using the fact that the velocity profile should be continuous, one can now solve the Navier-Stokes equation, Eq. (16):

[TABLE]

with

[TABLE]

One can now use this result to determine the thermodynamic fluxes in terms of these quantities. The velocity flux can be derived from Eq. (9),

[TABLE]

while the heat flux, Eq. (15), is given by

[TABLE]

This completes our general analysis of the thermodynamic fluxes of a gas in a tube, as we have determined the two independent fluxes purely in terms of external parameters. To make further progress, one needs to specify the driving. This will be done in the next section to determine the Onsager coefficients.

III.3 Onsager coefficients

With these ingredients, we are now ready to determine the off-diagonal Onsager coefficients. We first calculate the heat flux induced by a pressure gradient. As the force on every particle is equal, the pressure gradients near the wall and in the bulk are given by

[TABLE]

respectively. This leads to an Onsager coefficient,

[TABLE]

On the other hand, one can try to derive the transport coefficient associated the thermo-osmotic slip directly, by studying the situation in the presence of a temperature gradient but in the absence of a bulk pressure gradient,

[TABLE]

or

[TABLE]

The gradient of the density near the wall can be calculated as

[TABLE]

[TABLE]

and therefore, Eq. (24) simplifies to

[TABLE]

which determines the thermo-osmotic slip coefficient

[TABLE]

One immediately verifies that the expressions for the off-diagonal Onsager coefficients $L_{12}$ and $L_{21}$ are equivalent and therefore, Onsager symmetry is indeed valid.

In conclusion, we have determined the transport coefficient associated with the thermo-osmotic slip fully in terms of predetermined variables. Furthermore, one can see from the above calculation that the thermo-osmotic slip can be determined even in the absence of a temperature gradient due to Onsager symmetry.

III.4 Comparison with Derjaguin’s method

The standard method to determine the thermo-osmotic slip in a fluid is by using Dejaguin’s method, which states that the slip induced by a single wall is determined by the excess enthalpy, $\Delta h(y)=h(y)-h_{\textrm{B}}$

[TABLE]

In the system under study we have two walls with limited interaction range, leading to a Derjaguin Onsager coefficient equal to Derjaguin et al. (1987); Ganti, Liu, and Frenkel (2017)

[TABLE]

This raises the question whether our result is compatible with Derjaguin’s method. One can verify that this is indeed the case if one takes into account the extra assumption $\Delta L\ll L_{y}$.

IV Numerical results

With an explicit prediction for the thermo-osmotic slip at hand, we are now ready to test the theory with numerical simulations. For realistic hydrodynamic simulations, we rely on a stochastic rotation dynamics (SRD) algorithm. We shall first discuss the implementation of the SRD algorithm. Subsequently, we will use the algorithm to determine the thermo-osmotic slip coefficient $L_{21}$.

IV.1 Stochastic rotation dynamics

In this algorithm, the tube of particles is partitioned into square simulation cells with length $a$. At every time-step, all particles in a cell get a new velocity, $\mathbf{v}_{f}$ given by

[TABLE]

with $\mathbf{v}_{\textrm{cell}}$ the average velocity of the particles inside the simulation cell, $\mathbf{v}_{i}$ the velocity of the particle before the SRD step, and $\mathbf{R}$ a rotation matrix,

[TABLE]

$\alpha$ being a rotation angle which is fixed during a simulation. Furthermore, there are also particles in the wall with fixed temperature. This leads to a thermalization of gas particles with the wall. For a throughout review of the SRD algorithm, we refer to the literatureGompper et al. (2009).

To predict a value for the thermo-osmotic transport coefficient, one needs to determine the friction coefficient $\eta$. This coefficient is known for an SRD fluid Ihle, Tüzel, and Kroll (2005); Tüzel et al. (2003),

[TABLE]

where $\Delta t$ is the duration of one SRD simulation step.

All simulations will be done with $\alpha=0.722\pi$, $a=1$ and $\Delta t=0.1$. Furthermore, we shall take periodic boundary conditions in the x-direction, and set $L_{x}=10$, $L_{y}=50$, $N=5000$ $\Delta L=10$, $\Delta E=4$, $k_{B}T=10$ and $m=1$. The theoretically predicted value for the thermo-osmotic transport coefficient is then given by

[TABLE]

IV.2 Results

To determine the thermo-osmotic transport coefficient, we explicitly use Onsager symmetry, by applying a body force to the fluid particles and calculating the induced heat flux. This gives the thermo-osmotic transport coefficient,

[TABLE]

which should hold for small pressure gradients. In agreement with our theoretical analysis, we apply an equal force $f$ in the positive $x$ direction to every particle, which induces a pressure gradient $\nabla P=-f\overline{\rho}$. We vary $f$ between [math] and $1.5\cdot 10^{-3}$ or $-\nabla P=1.5\cdot 10^{-2}$, cf. Fig. 2. After averaging over $20$ runs of duration $t=4\cdot 10^{5}$, we get

[TABLE]

which is in good agreement with the theoretically predicted value. The small discrepancy between the two results can be explained from the approximations that we made in determining the theoretical value of the transport coefficient.

An alternative numerical route to obtain the thermo-osmotic coefficient is via the Green-Kubo relation, Eq. (5). As the infinite integration cannot be done analytically, we need to introduce a cut-off time for the simulation. From Fig. 3, one can observe that $\langle(Q(0)-h^{B}v_{x}(0))v_{x}(t)\rangle\approx 0$ if $t\gtrsim 600$. Therefore, we shall calculate

[TABLE]

Averaging over $1.6\cdot 10^{5}$ runs gives

[TABLE]

which again agrees reasonably well with the theoretical predictions, although the precision is one order of magnitude worse than the Onsager method.

Whenever a temperature gradient is applied over the tube, the particles near the wall feel an effective body force related to the excess enthalpy Ganti, Liu, and Frenkel (2017),

[TABLE]

Note that this is indeed in agreement with Eqs. (29) and (32). Explicitly applying this force (in the absence of a real temperature gradient) to the system gives a third computational method to calculate the thermo-osmotic slip coefficient. If we set $-\nabla T=0.00125$ in Eq. (45), and let the system run $20$ times for a duration of $t=4\cdot 10^{5}$ we get

[TABLE]

This again is in agreement with the theoretical prediction.

V Conclusion

In this work, we have studied thermo-osmosis in a simple model consisting of an ideal gas. In particular, we have compared a theoretical approach based on hydrodynamic arguments with computational approaches based on Onsager symmetry, on the Green-Kubo relation and on the excess enthalpy. The resulting thermo-osmotic slip coefficients are in good agreement with each other, implying that mesoscopic simulations are indeed consistent with local continuum descriptions.

Acknowledgements.

KP is a postdoctoral fellow of the Research Foundation-Flanders (FWO). The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government department EWI