.. _DPD_Exercise3:
DPD Exercise 3: Transport properties of DPD fluids
==================================================
Introduction
------------
One important aspect of the DPD thermostat is its ability to control temperature while ensuring the hydrodynamics of the fluid represented in the simulation are correct. This makes DPD suitable for studying systems with applied flow fields such as pressure-driven flow or linear shear.
The dissipative force parameter :math:`\gamma` in the DPD thermostat is one way we can control the dynamic viscosity :math:`\mu`, an important macroscopic property of a fluid. It is limited, however, as the relationship between :math:`\gamma` and :math:`\mu` is complicated - even when conservative forces are omitted - and might not be completely suitable for flows of liquids.
To get around this problem, we can make use of other pairwise thermostats. Two interesting ones for us are the Lowe-Andersen [Lowe1999]_ and Stoyanov-Groot [Stoyanov2005]_ thermostats. These randomly select particle pairs whose relative velocities are replaced with values chosen from a Maxwell-Boltzmann distribution for the required temperature. The dissipative force parameter :math:`\gamma` is replaced with a collision frequency :math:`\Gamma` that determines the probability of a particle pair having its relative velocity changed.
To determine the viscosity of a fluid from simulations, we need to find the relationship between shear stress and shear rate (velocity gradient). We can make use of a special boundary condition available in DL_MESO_DPD. The Lees-Edwards periodic boundary condition [Lees1972]_ can apply a linear shear flow field to a particle simulation by changing the velocity and tangential position of a particle when it moves through the boundary. This applies a constant velocity gradient throughout the box, giving a constant shear rate for the system.
We can calculate a stress tensor for the system at a given timestep by summing products of components of pairwise forces and vectors between particle pairs, i.e.
.. math :: \sigma_{\alpha \beta} = \sum_i \left(m_i v_{i,\alpha} v_{i,\beta} + \sum_{j>i} F_{ij,\alpha} r_{ij,\beta} \right), ~~~~~ (\alpha, \beta = x, y, z)
which can be averaged over a reasonably large number of timesteps to obtain representatives values. Only one of the off-diagonal terms is needed, e.g. if the :math:`x`-component of velocity varies along :math:`y`, we will need to look at :math:`\sigma_{yx}`. This value is often divided by system volume to give the pressure tensor :math:`P_{yx} = \sigma_{yx}/V`.
The viscosity of a fluid is equal to the derivative of shear stress :math:`\tau = |\sigma_{yx}|` with respect to shear rate :math:`\dot{\gamma} = \left|\frac{\partial v_x}{\partial y}\right|`:
.. math :: \mu = \frac{d \tau}{d \dot{\gamma}},
:label: viscosity
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.
Please refer to :ref:`DPD_Ex3Theory` if you would like more details about the theoretical background to this Exercise.
Aim
---
We are going to use the Lees-Edwards boundary conditions to produce linear shear in a simple DPD fluid and work out its viscosity. We are going to use the Stoyanov-Groot thermostat to begin with, but we will switch to the DPD thermostat later on.
Instructions
------------
Copy the CONTROL and FIELD input files from the directory ``DEMO/DPD/ShearFlow`` in your copy of DL_MESO to your working directory.
These simulation input files will model a cubic box (with dimensions :math:`10 \times 10 \times 10`) with Lees-Edwards shearing applied to the boundaries orthogonal to the :math:`y`-axis, containing 3000 particles of a fluid with a conservative force parameter :math:`A_{ij}=25` (i.e. the standard value for water). The simulation will last for 50,000 timesteps: we will only be collecting trajectory data in the HISTORY file after 40,000 timesteps have passed, but the data will include particle velocities as well as positions.
We have specified in the CONTROL file that we are using the Stoyanov-Groot thermostat, setting the collision frequency :math:`\Gamma = 1.0` in the FIELD file (last number specified for an interaction between species pairs). [#]_
.. [#] With the timestep :math:`\Delta t = 0.01` given in the CONTROL file, this means the probability of a particle pair having its relative velocity changed is :math:`\Gamma \Delta t = 0.01`. The other parameter for the Stoyanov-Groot thermostat :math:`\alpha` is set to a default value of 0.3: we are not going to worry about this one!
The FIELD file currently specifies that one of the shearing boundaries will produce a velocity of :math:`v_x = 1.0`, while the opposing boundary will produce an equal and opposite velocity (:math:`v_x = -1.0`). This means the current shear rate is the difference in velocities divided by the length of the box in the :math:`y`-dimension, i.e. :math:`\frac{\partial v_x}{\partial y} = \frac{\Delta v_x}{L_y} = \frac{1.0 - (-1.0)}{10} = 0.2`.
Tasks
-----
#. Run DL_MESO_DPD in your working directory with the supplied input files.
- If you take a look at the OUTPUT file, you may notice that the average system temperature will be significantly higher than the value requested in the CONTROL file. Do not panic: we are applying a velocity field in a particular direction (:math:`x`) and basing our temperature calculations on kinetic energy, so the measured temperature will naturally be higher than usual. If you look at the partial temperatures in the directions other than the flow field (:math:`y` and :math:`z`), you should find the averages of these should be more or less correct.
- Run the ``local.exe`` utility and split the box into several slices along the :math:`y`-direction (e.g. 100) but only one in the :math:`x`- and :math:`z`- directions: you will only need the time-averaged results for this. Open Paraview, load in the averages.vtk file and apply the Plot Over Line filter along the :math:`y`-axis to get a graph of the :math:`x`-component of velocity. Confirm that the time-averaged velocity profile is linear and work out the actual velocity gradient.
- We now want to take a look at the various pressure tensor components, particularly the :math:`yx`-component as this relates to the dimensions of the velocity and its change in value. Either look at the CORREL file created by DL_MESO_DPD using either your favourite plotting software - Excel, Gnuplot, :download:`correl.py ` etc., or take a look at the values at the end of the OUTPUT file (which are separated out into conservative, random, dissipative and kinetic contributions). **Don't forget to multiple these pressure tensor values by the system volume to get the stress tensor!**
- Find the average value of :math:`\sigma_{yx}` and plot the *absolute* value of this against the averaged velocity gradient.
#. Change the velocity of the shearing boundary in the FIELD file to modify the velocity gradient and repeat the calculation. Find the averaged stress tensor and velocity gradient for this calculation and add this to the plot.
#. 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 (velocity gradient) and shear stress?
- Try fitting a straight line to the plot using regression analysis and find the viscosity of the fluid by using equation :eq:`viscosity` (i.e. the gradient of this plot).
#. Change the collision frequency :math:`\Gamma` in the FIELD file and repeat the above calculations of shear stress against shear rate (steps 1 to 3) to find the new fluid viscosity.
#. Try a few different collision frequencies 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 its viscosity at room temperature (298 K) is :math:`\mu = 8.90 \times 10^{-4}` Pa s or approximately :math:`53.1` in DPD-based units (based on one molecule of water per particle), can you work out an appropriate value of :math:`\Gamma` for water?
Given the limited time available, you might like to either share the following tasks among you or select just one of these to concentrate on for this practical session.
#. Try changing the conservative force parameter :math:`A_{ij}` for the fluid and recalculate the viscosity. Can you see if this makes any noticeable difference to the measured viscosity if you use the same value for the collision frequency :math:`\Gamma`?
#. Try using :math:`A_{ij} = 0` and find the relationship between viscosity and collision frequency. How similar is this to the relationship for the original conservative force parameter?
#. Replace the line ``ensemble nvt stoyanov 0.3`` with ``ensemble nvt dpdvv`` in the CONTROL file. This will implement a scheme that applies the DPD thermostat more accurately by recalculating the dissipative forces at the end of the timestep.
- Repeat the calculations for viscosity at several shear rates for :math:`\gamma = 4.5` (a frequently used value for the dissipative force parameter). How does the magnitude of the viscosity compare with that obtained with the Stoyanov-Groot thermostat?
- Try varying the value of :math:`\gamma` for the fluid with conservative force parameters :math:`A_{ij}` of 0 and 25. How does the viscosity change with :math:`\gamma` in each case? Does a higher value of :math:`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?
.. [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.