LBE Tutorial Exercise 2: Multi-Component LBE simulations (Part 2)¶

Code transformation¶

The real interest in MCLB applications is in their dynamics; we want to simulate e.g. immiscible fluids interacting with each other. Imagine an oil drop in water. If the ‘external’ fluid (water here) moves, it exerts stresses on the surface of the oil and causes it to move, while the interior fluid (oil) impedes the motion of the water: action and reaction and all that. The resulting dynamics are complex: given the interface (with its own physics, as you can see in Ex1 in the previous notebook) deforms, re-locates the boundary between fluids, modifies the interaction … hopefully you get the point! With MCLB, all of these physics are in place as you hopefully saw in CCP5_02.c … we will now try and convince you by switching to a more comprehensive LBE code.

DL_MESO_LBE, the LBE code in DL_MESO, is an ideal tool to observe complex, multi-component hydrodynamics. It can work in 3D but we will stick to 2D to allow fast computation and interaction/experimentation. In the last set of exercises, we shall configure and use DL_MESO_LBE to perform simulations of a neutrally buoyant drop of red fluid (oil) suspended in a blue fluid (water), which is moved at its boundaries to produce a shear flow (as seen in Figure 1 below).

Boundary conditions for sheared drop system

-:

Figure 1: Schematic of a sheared drop (Fluid1). Flow is applied to the boundaries of the external fluid (Fluid2) and the fluid motion is communicated throughout the domain, mainly by viscous forces.

We shall perform some experiments using DL_MESO_LBE - which we can view using scripts and ParaView - allowing us to observe the hydrodynamics of the system, both in terms of the change in the drop’s shape and the complex flow patterns that arise.

Ex4. Equilibration dynamics with DL_MESO_LBE¶

We have supplied a ZIP file with the source code for both of DL_MESO’s codes. Unpack this ZIP file using the command:

[ ]:

!unzip -o -q dl_meso_2.8.zip

and then go into its dl_meso/WORK directory and compile DL_MESO_LBE for a single processor core with OpenMP multithreading (to speed up the calculations) using the following commands for the provided Makefile:

[ ]:

%%bash
cd dl_meso/WORK
make -f Makefile-LBEOMP clean
make -f Makefile-LBEOMP

This provides an executable file for DL_MESO_LBE (lbe.exe) that we can use to run a simulation, provided we have the required input files.

The input files for our sheared drop calculation are in the LBE2Ex4 directory: lbin.sys provides simulation system controls (including values for \(\alpha\) and boundary velocities), lbin.spa specifies boundary conditions and lbin.init provides an initial state, which adds a drop of fluid 1 (our ‘red’ fluid) with a radius of 10 right in the middle of fluid 0 (our bulk ‘blue’ fluid). To launch the calculation, either use the following commands (ideally in a terminal window):

[ ]:

%%bash
cd LBE2Ex4
rm -f *.vts
export OMP_NUM_THREADS=6
../dl_meso/WORK/lbe.exe

or the following script to launch DL_MESO_LBE in the specified folder and keep track of its progress.

[ ]:

import os
import launchdlmeso as dlm
os.environ["OMP_NUM_THREADS"] = "6"
dlm.run_LBE('LBE2Ex4', 'dl_meso/WORK/lbe.exe', 1, True, '')

Note that we are setting the number of OpenMP threads in the commands just before launching DL_MESO_LBE. While it might be tempting to use all available threads, overheads involved in splitting sections of the calculation can actually slow things down, so a balance is required. We recommend using 6 threads if you have at least that many available on your computer.

DL_MESO_LBE will print a summary of the calculation either to the screen (if running in a terminal window) or to a file (if using the script), and generate a series of structured grid VTK files (e.g. lbout000000.vts) after an equilibration period of 5000 timesteps to allow the drop to settle first. These files are snapshots of the simulation that can be opened and visualised in ParaView. (We also have some Python scripts to look at these.)

Take a look at the simulation output files. If you are using ParaView, you should be able to generate animations of e.g. drop fluid density. The Calculation filter can be used to find properties derived from those supplied in the files, e.g. total fluid densities \(\rho\) can be calculated using:

density_0+density_1

and phase indices \(\rho^N\) using:

(density_1-density_0)/(density_0+density_1)

as expressions in the Calculation filter. The Contour filter can be used on phase indices to find the centre of the interface between the two fluids (i.e. where \(\rho^N = 0\)) and plot the drop shape.

Alternatively, we can use the following script to read in a specific output file and plot the total densities, phase indices and velocity (both the modulus and the x- and y-components), as well as show the interface between the fluids in the velocity plots as a black line.

[ ]:

from plotvtk import *
plotVTK('LBE2Ex4/lbout000200.vts')

What happens to the drop when the shearing boundaries (at the top and bottom of the domain) are applied?

The lbin.sys file consists of keywords that allow you to change operating parameters for the calculation. The most relevant keywords for this exercise are:

speed_top_0 - the x-component of velocity for the boundary at the top of the domain
speed_bot_0 - the x-component of velocity for the boundary at the bottom of the domain
interaction_0_1 - the interaction parameter (interfacial tension \(\alpha\)) between fluids 0 (blue) and 1 (red)
segregation - the segregation parameter \(\beta\) used to ‘recolour’ the fluids after collisions
interaction_type - the MCLB algorithm we want to use, currently set to Lishchuk (the name of the algorithm in DL_MESO_LBE)

Note that the value of \(\beta\) we are using in DL_MESO_LBE is around double the value we used in CCP5_02.c. This is due to a slight difference in how segregation is carried out between the two codes, although the general trends regarding interfacial width and stability will be similar.

For this system, we can define a drop Reynolds number:

\[Re = \frac{R^2 \dot{\gamma}}{\nu}\]

and capillary number (ratio of viscous forces to Laplace or interfacial tension forces):

\[Ca = \frac{\rho \nu U}{\alpha} = \frac{R \rho \nu \dot{\gamma}}{\alpha}\]

where \(R\) is the drop radius and \(\dot{\gamma}\) is the local fluid shear rate in the region of the drop (giving a characteristic velocity \(U = R \dot{\gamma}\)).

Try changing the wall velocities and/or the interfacial tension between the fluids. What happens if the shear rate decreases or the interfacial tension increases (i.e. \(Ca\) decreases)?
What happens if you increase \(Ca\) (increase shear rate or decrease interfacial tension)?
Can you find a minimum value of \(Ca\) at which the drop starts to break up due to the flow?

The very first output file (lbout000000.vts) is written immediately after equilibration and should provide a good representation of a static drop. Try visualising the drop with this file and use the following script to find a few related properties.

[ ]:

from plotvtk import *
findCOMradius('LBE2Ex4/lbout000000.vts')

Take a look at the velocity field: where are the maximum velocities in relation to the drop’s location?
What happens to these microcurrents when you vary \(Ca\)?

Ex5. Algorithm change to ‘Spencer tensor’¶

DL_MESO_LBE includes a number of different MCLB algorithms, including three more variants of the Lishchuk algorithm you have used so far in this tutorial. While this form of the algorithm is robust and provides low microcurrents at fluid-fluid interfaces, it currently requires two sets of gradients to be calculated per lattice point for each timestep using stencils: one for the interfacial normal and one for the curvature. This will require additional core-to-core MPI communications of these properties when running DL_MESO_LBE on multiple-core machines, making this method less scalable as a result.

The other variants of the Lishchuk algorithm available in DL_MESO_LBE have been devised to get around this problem. One variant, known as ‘Spencer tensor’, removes the need to calculate interfacial curvatures by replacing the interfacial force calculations with direct forcing terms applied to collisions:

\[F_i = \frac{t_i \beta \alpha \omega \rho_0 \rho_1}{c_s^4 \rho^3 \Delta t} \left(\mathbf{nn} - \mathbf{I}\right) : \left(\mathbf{c}_i \mathbf{c}_i - c_s^2 \mathbf{I} \right)\]

While it still requires the interfacial normals \(\textbf{n}\) to be calculated at each lattice point, it halves the required number of stencil-based gradient calculations and associated MPI communications when running on multiple processor cores.

We can go further and replace the stencil-based calculations of phase index gradients to directly use distribution functions at each lattice point, i.e.

\[\nabla \rho^N (\textbf{x}) \approx -\frac{1}{c_s^2 \Delta t} \sum_i t_i \hat{c}_i \left(\frac{R_i (\textbf{x}) - B_i(\textbf{x})}{R_i(\textbf{x}) + B_i(\textbf{x})}\right)\]

These are too approximate to be used directly in interfacial forces calculations, but they are sufficient for interfacial normals (unit vectors) and, along with the forcing term above, can provide a fully-local Lishchuk algorithm.

For this exercise, we will stick to the original Lishchuk algorithm and the ‘Spencer tensor’ variant with stencil-based interfacial normal calculations. To get DL_MESO_LBE to use this latter variant, change the following lines in lbin.sys:

collision_type               BGKGuo
interaction_type             Lishchuk

to:

collision_type               BGK
interaction_type             LishchukSpencerTensor

Note that we do not need Guo forcing when using the Spencer tensor variant, since we are no longer calculating interfacial forces!

Try out the input files from the previous exercise with DL_MESO_LBE after making the above changes. Do you get similar results to before? Are the calculations any faster? How do the microcurrents during equilibration compare between the two codes?

One advantage of the implementation of the Lishchuk MCLB algorithms in DL_MESO_LBE is we can extend them to more than two fluids, meaning we can have more than one drop in a background fluid. We can define different phase indices for the drops (with the background fluid and each other) and associated interfacial normals, determine curvatures between the fluids where they interact and calculate forces between each pair of fluids … or work out forcing terms for our ‘Spencer tensor’ implementation.

We have a system available with three fluids - a background fluid (0) and two immiscible drops (fluids 1 and 2) - available in the LBE2Ex5 directory: lbin.sys, lbin.spa and lbin.init. This is similar to the one from the previous exercise, in that we are subjecting the system to linear shear, but we have ‘accidentally’ (on purpose) put the two drops in at the start so they partially overlap each other. We are confident that the MCLB’s segregation step will manage to separate out the drops. (Famous last words!)

Set the calculation going for this system using DL_MESO_LBE. The input files are currently set up for the original Lishchuk algorithm, so make the above changes in lbin.sys to use the ‘Spencer tensor’ algorithm.

[ ]:

import os
import launchdlmeso as dlm
os.environ["OMP_NUM_THREADS"] = "6"
dlm.run_LBE('LBE2Ex5', 'dl_meso/WORK/lbe.exe', 1, True, '')

Once each version of the calculation has finished, visualise their results either by using ParaView or the scripts provided below. Note that we can define the phase index for the second drop using the following expression in ParaView’s Calculation filter:

(density_2-density_0)/(density_0+density_2)

although you may find for one of the cases that this expression on its own (and the similar one for the first drop) might not be enough to find the boundary between two conjoined drops. To get around this, we suggest adding a condition for each phase index to test how much the background and drop fluids contribute to the total density at a given lattice point, which we can do with an if statement in the expression for the Calculation filter, e.g.

if((density_0+density_2)/(density_0+density_1+density_2)>0.5, (density_2-density_0)/(density_0+density_2), -1.0)

[ ]:

from plotvtk import *
plotVTK('LBE2Ex5/lbout000200.vts')
findCOMradius('LBE2Ex5/lbout000200.vts')

What happens to the drops immediately after equilibration and during the simulation for both versions of the code (or both algorithms)?
Can you think why the two versions of the MCLB algorithm behave differently in this case? (Hint: how can a curvature between more than two fluids be calculated?)

Try extending the number of timesteps in the simulation by increasing the number given with total_step in the lbin.sys file, or try increasing the speeds of the top and bottom boundaries to increase the shear rate.

What happens when two moving drops collide into each other? How would or could you control that process? (Hint: Take a look in the lbin.sys file.)