LBE Tutorial Exercise 2: Multi-Component LBE simulations¶

We are now extending our single-component LBE calculations from the previous exercise to look at systems with additional fluid components, using an approach for multi-component lattice Boltzmann (MCLB) to model how different fluids interact and affect their flows.

We will start with the custom simulation program CCP5_02.c, which has fully periodic boundary conditions and models a 2D bulk with two lattice Boltzmann fluids. It models a ‘red’ fluid drop at rest in the centre of a square lattice; periodic boundary conditions make this into an array of identical drops. It generates a single CSV (comma-separated values) file called density.csv and four other files (data_density.dat, data_rhoN.dat, data_nK.dat and data_psi.dat) that report the overall fluid density \(\rho = \rho_R + \rho_B\), the phase index \(\rho^N = \frac{\rho_R - \rho_B}{\rho}\), the interfacial curvature \(K\) and the velocity modulus (speed) \(|\mathbf{u}|\) at each grid point for the final timestep. The CSV file can be read into spreadsheet programs or used with our Python plotting scripts, while the others in fixed-length delimited format can be read into e.g. Gnuplot.

Ex1. Code familiarisation and Laplace law pressure step validation¶

Start by taking a look through CCP5_02.c, before compiling it using the following command:

[ ]:

!gcc -o LBE CCP5_02.c -lm

noting that -lm links in some required mathematics libraries. This should generate an executable called LBE that can be launched at a command-line. The code will carry out a number of timesteps (set by the value N, currently equal to 10000) and print each timestep number to the screen. To avoid producing a long stream of numbers, we have a Python script available to launch the calculation and collect those numbers in a scratch file (which we will delete afterwards), using these to keep track and display a progress bar:

[ ]:

from launchmclb import *
run_MCLB('LBE', 10000, '')

Once the calculation is complete, visualise the results using the following script:

[ ]:

from plotmclb import *
plotMCLB('density.csv')

You will see four properties:

the density \(\rho\) (pressure),
the phase field \(\rho^N\),
the interfacial curvature \(K\), and
the velocity modulus field \(|\mathbf{u}|\).

Note that the last of these will not be zero, even though the drop is at rest. The steady flow you observe is an artefact: the so-called ‘interfacial microcurrent’.

We can use the values of \(\rho^N\) to determine the location and the radius of the drop, and the values of \(\rho\) in the drop and background fluids to find the pressure drop between them. The following script takes the data from the output file, works out the centre-of-mass for the drop and then fits an ellipsoid function to points along \(\rho^N = 0\) to find the drop radius, works out the maximum and minimum densities corresponding to values inside and outside the drops (allowing us to calculate the interfacial tension), and plots the density profile through the drop.

[ ]:

from plotmclb import *
findCOMradius('density.csv')

The pressure difference between the two fluids in a circular drop should follow the two-dimensional form of the Young-Laplace equation:

\[p_R - p_B = c_s^2 \left(\rho_R - \rho_B\right) = \frac{\alpha}{R}\]

Is the Laplace pressure (the pressure difference between the fluids) consistent with the set value of the interfacial tension \(\alpha\) (DELTA2 in the code) and the (approximate) drop radius?
Try a few different values of the interfacial tension parameter, noting that you can input these directly into the code, and record characteristic values for the microcurrent velocity field, \(|\mathbf{u}_{max}|\). Is this velocity field related to the Laplace pressure?

Note that the observed interfacial microcurrents (spurious velocities) are an artefact of all numerical methods, not just LBE. The velocity should be zero at all points in the flow domain for a rest drop. These fluctuations set practical bounds on accessible values for the 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}\)).

What sets this MCLB method apart from others is that the interfacial tension \(\alpha\) is directly parameterised rather than being an emergent property, and can be much larger than those obtained from other MCLB variants, producing smaller interfacial microcurrents.

Ex2. Fluid component segregation, or re-colouring¶

A key aspect for our MCLB method is carrying out fluid collisions together and then segregating (or ‘re-colouring’) the fluids afterwards: this determines the quality of the hydrodynamics which emerges.

Continue working with CCP5_02.c, only now take a look at the segregation parameter \(\beta\), given as a symbolic constant BETA_LKR near the top of the code.

Try reducing the value of \(\beta\) and repeat the calculations and analyses for a few different values of \(\beta\). You will hopefully observe: (1) the phase field boundary is a \(\tanh\) profile as we predicted, (2) the microcurrents vary in inverse proportion to interfacial width, and (3) interfacial width is inversely proportional to \(\beta\). (Note that you may need to increase the simulation lattice size and drop radius to accommodate smaller \(\beta\) values.)
Slightly increase the value of \(\beta\). What happens? Why?

Ex3. Curvature calculations¶

The calculation of interfacial curvatures in our code is a little inefficient and we can try to simplify it. We can define our interfacial curvature as the divergence of the interfacial normal, i.e.

\[K = \frac{\partial n_x}{\partial x} + \frac{\partial n_y}{\partial y}\]

Edit the calc_obs_2() function in the code to simplify the definition of \(K\) based on the above. Compare the microcurrent activity from your new definition with the original one. (Can you still obtain a stable solution with this new curvature definition?)

To proceed with the final part of this Exercise, start this notebook.