| 56 |
|
(RNEMD) obtains transport coefficients (thermal conductivity and shear |
| 57 |
|
viscosity) in a fluid by imposing an artificial momentum flux between |
| 58 |
|
two thin parallel slabs of material that are spatially separated in |
| 59 |
< |
the simulation cell.\cite{MullerPlathe:1997xw,Muller-Plathe:1999ek} The |
| 59 |
> |
the simulation cell.\cite{MullerPlathe:1997xw,ISI:000080382700030} The |
| 60 |
|
artificial flux is typically created by periodically ``swapping'' either |
| 61 |
|
the entire momentum vector $\vec{p}$ or single components of this |
| 62 |
|
vector ($p_x$) between molecules in each of the two slabs. If the two |
| 201 |
|
hot and cold slabs are separated along the z-axis, the energy flux is |
| 202 |
|
given simply by the decrease in kinetic energy of the cold bin: |
| 203 |
|
\begin{equation} |
| 204 |
< |
(1-x^2) K_c^x + (1-y^2) K_c^y + (1-z^2) K_c^z = J_z |
| 204 |
> |
(1-x^2) K_c^x + (1-y^2) K_c^y + (1-z^2) K_c^z = J_z\Delta t |
| 205 |
|
\end{equation} |
| 206 |
|
The expression for the energy flux can be re-written as another |
| 207 |
|
ellipsoid centered on $(x,y,z) = 0$: |
| 208 |
|
\begin{equation} |
| 209 |
< |
x^2 K_c^x + y^2 K_c^y + z^2 K_c^z = (K_c^x + K_c^y + K_c^z - J_z) |
| 209 |
> |
x^2 K_c^x + y^2 K_c^y + z^2 K_c^z = K_c^x + K_c^y + K_c^z - J_z\Delta t |
| 210 |
|
\label{eq:fluxEllipsoid} |
| 211 |
|
\end{equation} |
| 212 |
|
The spatial extent of the {\it flux ellipsoid equation} is governed |
| 221 |
|
|
| 222 |
|
One may also define momentum flux (say along the x-direction) as: |
| 223 |
|
\begin{equation} |
| 224 |
< |
(1-x) P_c^x = j_z(p_x) |
| 224 |
> |
(1-x) P_c^x = j_z(p_x)\Delta t |
| 225 |
|
\label{eq:fluxPlane} |
| 226 |
|
\end{equation} |
| 227 |
|
The above {\it flux equation} is essentially a plane which is |
| 273 |
|
an approach similar to the above would be sufficient for this as well. |
| 274 |
|
|
| 275 |
|
\section{Computational Details} |
| 276 |
< |
Our simulation consists of a series of systems. |
| 276 |
> |
Our simulation consists of a series of systems. All of these |
| 277 |
> |
simulations were run with the OOPSE simulation software |
| 278 |
> |
package\cite{Meineke:2005gd} integrated with RNEMD methods. |
| 279 |
|
|
| 280 |
|
A Lennard-Jones fluid system was built and tested first. In order to |
| 281 |
|
compare our method with swapping RNEMD, a series of simulations were |
| 282 |
|
performed to calculate the shear viscosity and thermal conductivity of |
| 283 |
< |
argon. 2592 atoms were in a orthorhombic cell, which was ${10.06 \sigma |
| 284 |
< |
\times 10.06 \sigma \times 30.18 \sigma}$ by size. The reduced density |
| 285 |
< |
${\rho^* = \rho\sigma^3}$ was thus 0.849, which enabled direct |
| 286 |
< |
comparison between our results and others. |
| 283 |
> |
argon. 2592 atoms were in a orthorhombic cell, which was ${10.06\sigma |
| 284 |
> |
\times 10.06\sigma \times 30.18\sigma}$ by size. The reduced density |
| 285 |
> |
${\rho^* = \rho\sigma^3}$ was thus 0.85, which enabled direct |
| 286 |
> |
comparison between our results and others. These simulations used |
| 287 |
> |
Verlet time-stepping algorithm with reduced timestep ${\tau^* = |
| 288 |
> |
4.6\times10^{-4}}$. |
| 289 |
|
|
| 290 |
|
For shear viscosity calculation, the reduced temperature was ${T^* = |
| 291 |
< |
k_B T / \epsilon = 0.72}$. Simulations were run in microcanonical |
| 291 |
> |
k_B T/\varepsilon = 0.72}$. Simulations were run in microcanonical |
| 292 |
|
ensemble (NVE). For the swapping part, Muller-Plathe's algorithm was |
| 293 |
|
adopted.\cite{ISI:000080382700030} The simulation box was under |
| 294 |
< |
periodic boundary condition, and devided into 20 slabs. In each swap, |
| 295 |
< |
the top slab ${(n = 0)}$ exchange the most negative $x$ momentum with the |
| 296 |
< |
most positive $x$ momentum in the center slab ${(n = N/2)}$. Referring |
| 294 |
> |
periodic boundary condition, and devided into ${N = 20}$ slabs. In each swap, |
| 295 |
> |
the top slab ${(n = 1)}$ exchange the most negative $x$ momentum with the |
| 296 |
> |
most positive $x$ momentum in the center slab ${(n = N/2 + 1)}$. Referred |
| 297 |
|
to Tenney {\it et al.}\cite{tenneyANDmaginn}, a series of swapping |
| 298 |
< |
frequency were chosen. Corresponding to each result from swapping |
| 298 |
> |
frequency were chosen. According to each result from swapping |
| 299 |
|
RNEMD, scaling RNEMD simulations were run with the target momentum |
| 300 |
< |
flux parameter set to produce a similar momentum flux and shear |
| 301 |
< |
rate. Furthermore, various scaling frequencies and corresponding flux |
| 302 |
< |
can be tested for one swapping rate. |
| 300 |
> |
flux set to produce a similar momentum flux and shear |
| 301 |
> |
rate. Furthermore, various scaling frequencies can be tested for one |
| 302 |
> |
single swapping rate. To compare the performance between swapping and |
| 303 |
> |
scaling method, temperatures of different dimensions in all the slabs |
| 304 |
> |
were observed. |
| 305 |
|
|
| 306 |
< |
After each simulation, the shear viscosities were calculated in |
| 307 |
< |
reduced unit. |
| 306 |
> |
After each simulation, the shear viscosity was calculated in reduced |
| 307 |
> |
unit. The momentum flux was calculated with total unphysical |
| 308 |
> |
transferred momentum ${P_x}$ and simulation time $t$: |
| 309 |
> |
\begin{equation} |
| 310 |
> |
j_z(p_x) = \frac{P_x}{2 t L_x L_y} |
| 311 |
> |
\end{equation} |
| 312 |
> |
And the velocity gradient ${\langle \partial v_x /\partial z \rangle}$ |
| 313 |
> |
can be obtained by a linear regression of the velocity profile. From |
| 314 |
> |
the shear viscosity $\eta$ calculated with the above parameters, one |
| 315 |
> |
can further convert it into reduced unit ${\eta^* = \eta \sigma^2 |
| 316 |
> |
(\varepsilon m)^{-1/2}}$. |
| 317 |
|
|
| 318 |
+ |
For thermal conductivity calculation, simulations were first run under |
| 319 |
+ |
reduced temperature ${T^* = 0.72}$ in NVE ensemble. Muller-Plathe's |
| 320 |
+ |
algorithm was adopted in the swapping method. Under identical |
| 321 |
+ |
simulation box, in each swap, the top slab exchange the molecule with |
| 322 |
+ |
least kinetic energy with the molecule in the center slab with most |
| 323 |
+ |
kinetic energy, unless this ``coldest'' molecule in the ``hot'' slab |
| 324 |
+ |
is still ``hotter'' than the ``hottest'' molecule in the ``cold'' |
| 325 |
+ |
slab. According to swapping RNEMD results, target energy flux for |
| 326 |
+ |
scaling RNEMD simulations can be set. Also, various scaling |
| 327 |
+ |
frequencies can be tested for one target energy flux. To compare the |
| 328 |
+ |
performance between swapping and scaling method, distributions of |
| 329 |
+ |
velocity and speed in different slabs were observed. |
| 330 |
+ |
|
| 331 |
+ |
For each swapping rate, thermal conductivity was calculated in reduced |
| 332 |
+ |
unit. The energy flux was calculated similarly to the momentum flux, |
| 333 |
+ |
with total unphysical transferred energy ${E_{total}}$ and simulation |
| 334 |
+ |
time $t$: |
| 335 |
+ |
\begin{equation} |
| 336 |
+ |
J_z = \frac{E_{total}}{2 t L_x L_y} |
| 337 |
+ |
\end{equation} |
| 338 |
+ |
And the temperature gradient ${\langle\partial T/\partial z\rangle}$ |
| 339 |
+ |
can be obtained by a linear regression of the temperature |
| 340 |
+ |
profile. From the thermal conductivity $\lambda$ calculated, one can |
| 341 |
+ |
further convert it into reduced unit ${\lambda^*=\lambda \sigma^2 |
| 342 |
+ |
m^{1/2} k_B^{-1}\varepsilon^{-1/2}}$. |
| 343 |
+ |
|
| 344 |
+ |
\section{Results And Discussion} |
| 345 |
+ |
\subsection{Shear Viscosity} |
| 346 |
+ |
|
| 347 |
|
\section{Acknowledgments} |
| 348 |
|
Support for this project was provided by the National Science |
| 349 |
|
Foundation under grant CHE-0848243. Computational time was provided by |
