Marcus N. Bannerman School of Engineering, University of Aberdeen m.campbellbannerman@abdn.ac.uk

• Aberdeen and the University
• Event-Driven Molecular Dynamics (EDMD)
• Thermal conductivity in nanofluids

Aberdeen City

University of Aberdeen

Event-Driven Molecular Dynamics

I distribute an Event-Driven Molecular Dynamics package, called DynamO (dynamomd.org).

EDMD is a general molecular dynamic technique and it is used to study a wide range of fundamental model systems:

A binary hard-sphere fluid, used as an exploratory model for transport in nanofluids (more later).

A coarse-grained simulation of the binding of two heteropolymers. Event-driven simulations are popular in the protein folding/structure literature.

Asymmetric-dumbells or "Snowmen" molecules, used as an model for anisotropic colloids, or Na-Cl crystals, and have a rich phase behavior.

What is event-driven dynamics?

Although it is not as "popular" as time-stepping dynamics, event-driven dynamics is a natural approach to dynamics.

Consider how you would solve this problem:

How are event times calculated?

Assuming no drag, the exact solution to the projectile motion is a parabola: \begin{align*} r_y(t+\Delta t) = \color{aqua}{r_y(t)} + \color{lime}{v_y(t)}\,\Delta t+g_y\,\frac{\Delta t^2}{2} \end{align*}

\begin{align*} r_y(t+\Delta t) = \color{aqua}{r_y(t)} + \color{lime}{v_y(t)}\,\Delta t+g_y\,\frac{\Delta t^2}{2} \end{align*} When it hits the ground at a time $t+\Delta t$ we have: \begin{align*} r_y(t+\Delta t)=r_{ground} \end{align*}

To find the time of this "event", we are looking for the roots to the equation: \begin{align*} r_y - r_{ground} + v_y\,\Delta t+g_y\,\frac{\Delta t^2}{2} = 0 \end{align*}

In this case, this is a simple quadratic with two roots: \begin{align*} \Delta t = \frac{-v_y \pm \sqrt{v_y^2 - 4\,g_y(r_y-r_{ground})}}{2\,g_y} \end{align*}

Event-detection is driven by root finding.

Generalizing for multiple particles

Consider two dueling catapults:

We calculate $\Delta t$ for all possible events: \begin{align*} \vec{\Delta t}=\left\{\Delta t_{1-ground},\,\Delta t_{2-ground},\,\Delta t_{1-2}\right\} \end{align*} The earliest event is the next event!

Event processing is driven by sorting.

Event execution

What actually happens at an event?

Consider a two-body event between particle $i$ and particle $j$, which has an energy change of $\Delta E$.

Conservation of momentum requires: \begin{align*} m_i\,\Delta \vec{v}_i = - m_i\,\Delta \vec{v}_i = \Delta \vec{P} \end{align*}

Where the impulse, $\Delta \vec{P}$, is a solution to energy balance: \begin{align*}\scriptsize \frac{1}{2}m_i\, {v}_i^2 + \frac{1}{2}m_j\, v_j^2 + \Delta E = \frac{1}{2}m_i\, \left(\vec{v}_i +\frac{\Delta \vec{P}}{m_i}\right)^2 + \frac{1}{2}m_j\, \left(\vec{v}_j-\frac{\Delta \vec{P}}{m_i}\right)^2 \end{align*}

Event execution is driven by root finding for an impulse.

The EDMD algorithm

Note, we do not decide the time of the simulation! We decide how many events to run it for and the time is calculated.

The EDMD algorithm

Modern EDMD algorithms are extremely complex and include:

• Neighbor lists and particle region tracking.
• Bounded calendar priority queues for sorting.
• Lazy deletion methods for optimal event list updates.

Perhaps the most unusual optimization is the time warp algorithm.

As events only involve small, localized particles, the system can run out of sync. i.e., each particle is left at different points in time:

Initially, the simulation is run in synchronous mode. The time-warp algorithm is then enabled (with the same simulation speed) and only the required updates to the particles are shown.

Precision

What is unique to DynamO is its stability.

EDMD is an exact (to machine precision) solution of the dynamics of the model used.

All approximation is explicit in the model used. (no worries about the time step)

The high-precision of the method conserves detailed balance. This is ideal for mixed MC/MD techniques (e.g., replica exchange, multicanonical).

Consider repeated bounces of the projectile: Using a constant coefficient of restitution $e<1$ we have: \begin{align*} \Delta P_{i,y} = -m_i\,(1+e) v_{i,y} \end{align*}

The velocity on impact/event decays to zero, and an infinite number of events happens in a finite time.

This "inelastic collapse" is not an instability, but a phenomena of the model.

Trying to simulate an "inelastic collapse" in a naive way highlights a numerical issue.

Although $r_y$ should equal $r_{ground}$ after an event, limitations of the precision of the computer leave it at $\approx \pm 10^{-14}$.

As $v_y\to0$, numerical precision is lost and the argument of the square-root can turn negative! \begin{align*} \Delta t = \frac{-v_y \pm \sqrt{v_y^2 - 4\,g_y(r_y-r_{ground})}}{2\,g_y} \end{align*}

This gives no real-roots/events and the projectile falls through the ground!

Although this problem is acute in dissipative systems with gravity, it also plagues elastic/molecular systems for certain sequences of events:

Machine precision is not enough to stop particles entering "forbidden" states, which may get rapidly worse.

DynamO has a (small) compile-time computer algebra library and event searching routines:

Variable t;
double deltat = nextroot(rij + t * (vij  + aij * t / 2), t);
	      

Additional "stabilizing" events are added if particles are in invalid conditions and this situation is deteriorating.

The library can find all roots of $N$-order polynomials and certain classes of other functions in a stable way, allowing a wide range of models/events to be detected.

Can we simulate "normal" molecular models using EDMD?

Between discontinuities/events/impulses there is no force: Exact solution to the motion is possible!

Thomson, C., Lue, L. & Bannerman, MN. . 'Mapping continuous potentials to discrete forms'. J. Chem. Phys, 140, 034105 (2014)

These stepped approximations can closely reproduce the phase behavior of the LJ system:

Thomson, C., Lue, L. & Bannerman, MN. . 'Mapping continuous potentials to discrete forms'. J. Chem. Phys, 140, 034105 (2014)

And the dynamics, for example, the shear viscosity:

Thomson, C., Lue, L. & Bannerman, MN. . 'Mapping continuous potentials to discrete forms'. J. Chem. Phys, 140, 034105 (2014)

But we should only use this approach for gases as for liquids it is not as efficient as standard approaches (as expected).

Thomson, C., Lue, L. & Bannerman, MN. . 'Mapping continuous potentials to discrete forms'. J. Chem. Phys, 140, 034105 (2014)

• EDMD and standard MD are equivalent.
• EDMD performs exponentially better for low/gas density systems, but is slower for liquids.
• EDMD is more efficient for:
  • Stiff systems (hard sphere, WCA, spring).
  • Coarse-grained potentials (hard sphere, square-well).
  • Multi-timescale systems (granular systems).

Thermal conductivity of nanofluids

Nanofluids are suspensions of nano-sized particles (solid or liquid) in a base fluid.

They are used to attempt to "alloy" the properties of two phases without the macroscopic problems (abrasion, stability).

In particular, metallic nanoparticles can be used to attempt to increase the thermal conductivity of water.

Thermal conductivity of nanofluids

Early attempts to explore the thermal properties of nanofluids were extremely encouraging.

Dispersion of copper nanoparticles in ethylene glycol exhibited enormous thermal conductivity enhancements for low volume fractions of nanoparticle.

These enhancements were well beyond initial predictions, earning the title of an "anomalous" enhancement.

Thermal conductivity of nanofluids

A massive effort, from 30 different organizations, tested identical nanofluid samples with a variety of techniques as part of the International Nanofluid Property Benchmark Exercise.

Finally, in 2010, Eapen et. al. published a paper which highlighted that most results lay within the classical limits of series and parallel conductors (right).

Eapen concluded that there are no "anomalous" results.

An equilibrium simulation of a hard-sphere nanofluid at a size ratio of 10:1 and mass ratio of 1000:1.

Thermal conductivity, $L_{\lambda\lambda}$, of a 10:1 nanoparticle fluid as a function of nanoparticle mole fraction, $x_A$, at (a) fixed packing fractions, $\phi$, or (b) fixed pressures, $p$. Symbols are simulation results, lines are predictions from kinetic theory.1

Bannerman, M. N. & Lue, L., "Transport properties of highly asymmetric hard-sphere mixtures" J. Chem. Phys., 130, 164507 (2009)

Analyzing the kinetic theory predictions, the dominant contribution has this form: \begin{align*} J_q &\approx \left(C_v^\text{Nano} - C_v^\text{Base}\right)k_B\,T\, L_{A\lambda}\frac{\nabla T}{T} \end{align*}

• Negative thermophoresis causes the nanoparticles to move towards the hot plate.
• Conservation of momentum causes an equal mass of base fluid to move towards the cold plate:

\begin{align*} J_q &\approx \left(C_v^\text{Nano} - C_v^\text{Base}\right)k_B\,T\, L_{A\lambda}\frac{\nabla T}{T} \end{align*}

• Nanoparticles have $10^3$ times the mass of the base particles, causing a large ``pumping'' effect.
• As the mass of the nanoparticle is comparatively so high, its mass specific heat capacity is low. \begin{align*} C_v^\text{Nano} &= \frac{5}{2\,m_\text{Nano}} & C_v^\text{Base} &= \frac{5}{2\,m_\text{Base}} \end{align*}
• Therefore, thermophoresis also causes a ``pumping'' of heat iff nanoparticles and base fluids have different masses/Cv!
• This effect can oppose or enhance thermal conductivity, depending on the sign of $L_{A\lambda}$.

A non-equilibrium thermal conductivity simulation of a hard-sphere nanofluid at a size ratio of 10:1 and mass ratio of 1000:1. Thermophoresis causes the larger species to diffuse to the hot plate.

Isotopic 10:1 mass ratio non-equilibrium molecular dynamics measurements of the thermal conductivity. For this system, the equilibrium thermal conductivity is $L_{\lambda\lambda}=2.35$. Packing fraction is 0.1, mole fraction of the heavier component is 0.1.

Summary

• Nanofluids are an exciting development which can revolutionize heat transfer.
• Anomalous (en/de)hancements of $L_{\lambda\lambda}$ are possible!
• It arise from a transport of heat due to thermophoresis, provided there is a significant heat capacity difference between the nanoparticle and base fluid.
• The mechanism is transient which may cause the apparent difficulties in experimental measurements of the thermal conductivity.
• More data is needed on the thermal diffusivity of nanofluids, to verify this effect.