DPD Exercise 3: Transport properties of DPD fluids

Introduction

One significant advantage DPD has over molecular dynamics is its thermostat’s ability to both control system temperature and automatically ensure hydrodynamics are applied correctly [1]. If a flow field, e.g. linear shear or pressure-driven flow, is applied to a DPD simulation, the use of relative velocities between particle pairs in dissipative forces means the hydrodynamics are not affected.

The dissipative force parameter \(\gamma\) in the DPD thermostat is one way we can control the dynamic viscosity \(\mu\), an important macroscopic property of a fluid. It is limited, however, as the relationship between \(\gamma\) and \(\mu\) is complicated - even when conservative forces are omitted - and might not be completely suitable for flows of liquids.

Alternative pairwise thermostats to DPD can help us get around this problem, particularly the Lowe-Andersen [Lowe1999] and Stoyanov-Groot [Stoyanov2005] thermostats. The dissipative and random forces for DPD are replaced with changes to relative velocities \(\mathbf{v}_{ij}\) between pairs of particles to randomly-chosen values from a Maxwell-Boltzmann distribution for the required system temperature. The parameter for these thermostats, the collision frequency \(\Gamma\), replaces the dissipative force parameter \(\gamma\) for DPD and represents both the probability of a particle pair having its relative velocity changed and an analogue for fluid viscosity.

There are number of techniques to measure the viscosity of a fluid from DPD simulations based on non-equilibrium molecular dynamics (NEMD) techniques. Measuring viscosity using linear shear is one such approach: we can calculate resulting shear stresses from setting constant velocity gradients in a simulation box, with the relationship between the two providing the viscosity.

A boundary condition available in DL_MESO_DPD - the Lees-Edwards periodic boundary condition [Lees1972] - can apply a constant velocity gradient (a linear shear flow) throughout a simulation box. The products of pairwise forces and vectors between particle pairs divided by the system volume provide our stress tensor \(\sigma_{\alpha \beta}\), and values of this property averaged over a reasonably large number of timesteps gives a representative value for the flow.

Only one of the off-diagonal terms of the stress tensor is needed: for instance, if we are varying the \(x\)-component of velocity varies along \(y\), we will need to look at \(\sigma_{yx}\) [2]. This value will equal the negative of the shear stress \(\tau\) for the applied flow, i.e. \(\tau = -\sigma_{yx}\).

The viscosity of a fluid is equal to the derivative of shear stress \(\tau\) with respect to shear rate \(\dot{\gamma} = \left|\frac{\partial v_x}{\partial y}\right|\):

(1)\[\mu = \frac{d \tau}{d \dot{\gamma}},\]

which is generally a function of the shear rate. For simple fluids, often described as Newtonian, this value will be a constant value regardless of shear rate.

Aim

We are going to use the Lees-Edwards boundary conditions to produce linear shear in a simple DPD fluid. By trying out different shear rates (velocity gradients) and measuring the resulting shear stresses, we will be able to work out the fluid’s viscosity and see how it changes when we change the thermostat parameter.

We are mainly going to use the Stoyanov-Groot thermostat, but we can optionally switch to the DPD thermostat later on.

Instructions

For this exercise, you will need the main DL_MESO_DPD executable dpd.exe - either serial or parallel versions will do - and the local.exe utility to analyse and visualise some of the results. You will also need either graph plotting software (e.g. Gnuplot, Excel) or the dlmresultviewer.py Python script to look at CORREL files: the former might also be useful to plot some graphs of your own. The same Python script can be used to look at the outputs (VTK files) from the local.exe utility, although if you decide not to use this, you will also need to have Paraview installed.

Copy the CONTROL and FIELD input files from the directory dl_meso/DEMO/DPD/ShearFlow in your copy of DL_MESO to your working directory (dl_meso/WORK) and take a look at them in a text editor.

The CONTROL file is fairly similar to ones you have seen before in previous exercises, although there are three directives to pay attention to.

The first is the line used to specify writing trajectories to the HISTORY file:

trajectory 40000 10 1

which means DL_MESO_DPD will start writing to HISTORY from timestep 40,000 every 10 timesteps and include both particle positions and velocities.

The second directive is that for the ensemble:

ensemble nvt stoyanov 0.3

where we are telling DL_MESO_DPD we want to use the Stoyanov-Groot thermostat for a constant volume (NVT) system. We are also setting the Nosé-Hoover thermostat coupling parameter to 0.3, although we do not need to change this value in this exercise.

The last notable directive specifies that we wish to apply shear at the boundaries orthogonal to the \(y\)-axis:

surface shear y

The FIELD file species a single particle species (FLUID) that uses a conservative force parameter of \(A_{ij}=25\) - the frequently used value for one-molecule water beads - an interaction cutoff \(r_c = 1\) and a collision frequency \(\Gamma = 1\) [3]. The notable addition to this file is the following block:

EXTERNAL
shear 1.0 0.0 0.0

which specifies the wall velocity for the top wall (\(v_{w, x} = 1.0\)) that is undergoing shear. (The bottom wall will move with a velocity \(v_{w, x} = -1.0\).) Since the simulation box is 10 length units in the \(y\)-direction, as specified in the CONTROL file, the expected velocity gradient is \(\frac{\partial v_x}{\partial y} = \frac{\Delta v_x}{L_y} = \frac{1.0 - (-1.0)}{10} = 0.2\).

When you run DL_MESO_DPD, you may notice that the resulting system temperature (as reported in the OUTPUT file) will be higher than the value specified in the CONTROL file. This is not generally a problem for this kind of simulation: we are applying a flow field in the \(x\)-direction and DL_MESO_DPD does not directly account for this when calculating system temperature from kinetic energies. To help us check that the thermostat is working property, DL_MESO_DPD will also display partial temperatures for each direction. Since the other two directions will not include a flow field, the partial temperatures for these directions should each average out at the specified value.

Once DL_MESO_DPD has finished running, you will need to run local.exe to work out time-averaged local velocities and how those vary in the \(y\)-direction. We can do this using the following command:

./local.exe -nx 1 -ny 100 -nz 1 -av

which will divide the box in the \(y\)-direction into 100 slices - sufficient for our purposes - and write to a file called averages.vtk, which can be opened up in Paraview. After loading the file and clicking Apply in the Properties subwindow, apply the Plot Over Line filter, select the Y Axis and click Apply to display graphs of the properties that have been calculated by local.exe. Alternatively, you can open this same file using the dlmresultviewer.py script, which will plot graphs of the same properties and allow you to draw lines indicating average values and/or fit polynomial functions to the data.

Both the OUTPUT and CORREL files will display stress tensors calculated for the simulation - the former as time-averaged values separated into conservative (interaction), dissipative, random and kinetic contributions, the latter as instantaneous values.

Tasks

  1. Run DL_MESO_DPD in your working directory with the supplied input files, then run the local.exe utility to obtain time-averaged properties in slices along the \(y\)-direction. Either use the dlmresultviewer.py script or open Paraview, load in the averages.vtk file and apply the Plot Over Line filter along the \(y\)-axis to get a graph of the \(x\)-component of velocity.
    • Confirm that the time-averaged velocity profile is linear and work out the actual velocity gradient (shear rate) \(\frac{\partial v_x}{\partial y}\).
  2. Now take a look at the various stress tensor components using either plotting software or the Python script dlmresultviewer.py on the CORREL file, or the averaged values at the end of the OUTPUT file.
    • Find the average value of \(\sigma_{yx}\) and plot the absolute value of this (i.e. without positive or negative signs) against the averaged velocity gradient.
  3. Change the velocity of the shearing boundary in the FIELD file to modify the velocity gradient and repeat the calculation (steps 1 and 2). Find the averaged \(yx\)-component of the stress tensor and velocity gradient for this calculation and add this to your stress vs. stress plot.
  4. Once you have calculated shear stresses for at least five different shear rates, take a look at your plot.
    • What kind of relationship exists between the shear rate and shear stress?
    • Try fitting a straight line to the plot using regression analysis and find the viscosity of the fluid using (1) (i.e. the gradient of this plot).
  5. Change the collision frequency \(\Gamma\) in the FIELD file and repeat the above calculations of shear stress against shear rate (steps 1 to 4) to find the new fluid viscosity.
  6. Try a few different collision frequencies (at least three) and plot the viscosities you find against the collision frequencies.
    • What kind of relationship do you get? Try to find a suitable function for the plot using regression analysis.
    • Given that water’s viscosity at room temperature (298 K) is \(\mu = 8.90 \times 10^{-4}\) Pa s or approximately \(53.1\) in DPD-based units (based on one molecule of water per particle), can you work out an appropriate value of \(\Gamma\)?

The following tasks are optional but could offer more insight in how to obtain a required viscosity for a more complicated DPD system.

  1. Try switching off conservative forces by setting \(A_{ij}\) to zero and recalculate the fluid viscosity. Does this make any noticeable difference to the viscosity you measured with the original \(A_{ij}\) value for a given value of collision frequency \(\Gamma\)?
  2. Replace the line ensemble nvt stoyanov 0.3 with ensemble nvt dpdvv in the CONTROL file. This implements a scheme (known as ‘DPD Velocity Verlet’) that applies the DPD thermostat more accurately by recalculating the dissipative forces at the end of the timestep after force integration.
    • Repeat the calculations for viscosity at several shear rates for \(\gamma = 4.5\) (a frequently used value for the dissipative force parameter) and \(A_{ij} = 25\). How does the magnitude of the viscosity compare with that obtained with the Stoyanov-Groot thermostat?
    • Try varying the value of \(\gamma\) for the fluid with conservative force parameters \(A_{ij}\) of 0 and 25. How does the viscosity change with \(\gamma\) in each case? Does the higher value of \(A_{ij}\) lead to a higher viscosity? Given the required value of viscosity to represent water (shown above in DPD units), do you think the DPD thermostat can achieve it?

Footnotes

[1]Non-equilibrium molecular dynamics (NEMD) can model systems with correct hydrodynamics and constant temperatures using a standard thermostat, but in this instance the thermostat can only be applied to peculiar velocities, i.e. particle velocities without contributions from the flow field. Finding the flow field requires localised sampling of the flow field - averaging velocities in volume slices over a short period of time - which adds a significant computational cost to the simulation.
[2]Since off-diagonal stress tensor components are symmetrical, we could also look at \(\sigma_{xy}\).
[3]With the timestep \(\Delta t = 0.01\) given in the CONTROL file, this means the probability of a particle pair having its relative velocity changed is \(\Gamma \Delta t = 0.01\).

References

[Lowe1999]CP Lowe, An alternative approach to dissipative particle dynamics, EPL, 47, p. 145-151, 1999, doi: 10.1209/epl/i1999-00365-x.
[Stoyanov2005]SD Stoyanov and RD Groot, From molecular dynamics to hydrodynamics: A novel Galilean invariant thermostat, Journal of Chemical Physics, 122, 114112, 2005, doi: 10.1063/1.1870892.
[Lees1972]AW Lees and SF Edwards, The computer study of transport processes under extreme conditions, Journal of Physics C, 5, p. 1921-1928, 1972, doi: 10.1088/0022-3719/5/15/006.