D2Q9 - Multiphase flow: Rayleigh-Taylor instability

Before start, it is assumed that the reader is already familiar with:

Basic model creation in the TCLB code, see first tutorial.
Basics of the LB physics, see obligatory SRT LBM tutorial.
Basics of running TCLB

We will cover following topics:

Governing equations for multiphase flow
One of the phase field LBMs available in TCLB
Using RTemplate to write multiple input files
Using globals to track interface locations

This tutorial re-creates the Rayleigh-Taylor instabilities in:

Fakhari, A., Mitchell, T., (2017). Improved locality of the phase-field lattice-Boltzmann model for immiscible fluids at high density ratios, Physical Review E 96(5), 053301.

Governing equations for multiphase flow

In order to resolve the governing dynamics of multiphase (two-phase in this case) flow, this tutorial takes a 'single-fluid' approach. Here, we solve a single set of Navier-Stokes equations where the density and viscosity are determined based on a phase parameter that tracks the various fluids in our systems. The phase parameter evolves based on the Allen-Cahn equation. This is in contrast to 'multi-fluid' approaches in which phases are treated as interpenetrating media with their own respective set of Navier-Stokes equations weighted by a void fraction. The system of equations that we are resolving here are:

$\partial_t \rho + \nabla \cdot \rho \mathbf{u} = 0$ $\rho \left( \partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla p + \nabla \cdot (\mu [\nabla \mathbf{u} + (\nabla \mathbf{u})^T]) + \mathbf{F}_s + \mathbf{F}_b$ $\partial_t \phi + \nabla \cdot \phi \mathbf{u} = \nabla \cdot M \left[ \left(\nabla \phi - \frac{\nabla \phi}{|\nabla \phi|}\frac{1-4(\phi-\phi_{avg})^2}{W} \right) \right]$

Here, we have introduced a volumetric force that applies in the diffuse interface between the phases to describe the surface effects,

$\mathbf{F}_s = \mu_{\phi} \nabla \phi,$

where $ \mu_{\phi} $ is the chemical potential described as,

$\mu_{\phi} = 4\beta (\phi - \phi_L)(\phi-\phi_H)(\phi-\phi_{avg}) - \kappa \nabla^2 \phi.$

In this, $ \beta, \kappa $ are related to the surface tension, $\sigma$, and diffuse interface width, $W$, by,

$\beta = 12\sigma/ W,$ $\kappa = 3\sigma W/2.$

Other parameters of interest in the Allen-Cahn equation include the mobility, $M$, which we relate to the relaxation time for one of our lattice populations.

These equations are recovered through a double-distribution lattice Boltzmann method described in detail in the paper listed above. Here, we use a population, $g_i$, to recover the hydrodynamics and $h_i$ for the interfacial dynamics. The lattice Boltzmann equations for these can be written as,

$g_i(x+e_i\delta t, t + \delta t) = g_i(x,t) + \Omega_i + F_i$ $h_i(x+e_i\delta t, t + \delta t) = h_i(x,t) - \frac{ h_i - h_u^{eq} }{ \tau_{\phi} + 0.5 } + F_i$

The details of these equations are beyond the scope of this tutorial, but for those interested please read Fakhari et al. (2017) and Mitchell et al. (2018). In TCLB this model in 2D is the 'd2q9_pf_velocity' model and in 3D is the 'd3q27_pf_velocity' model.

The Rayleigh-Taylor instability (RTI)

The RTI is an instability of an interface between two-fluids caused when the heavier is situated above the lighter within a gravitational field. Due to the buoyancy forces created, the light fluid rises while the heavy fluid propagates downwards in a characteristic fashion. Some nice examples from Wikipedia include the behaviour of water suspended above oil, mushroom clouds from volcanic eruptions, supernova explosions where we see expanding gas accelerate into a denser shell layer of gas, as well as instabilities in plasma fusion reactors and inertial confinement fusion.

This instability is commonly used to benchmark multiphase lattice Boltzmann solvers due to the interesting dynamics associated in addition to the range of applications where it finds relevance. In this tutorial we will look at two cases; (i) low density ratio, (ii) high density ratio. In order to create this file, we are going to learn to use the RTemplate tool. In doing this, we are going to write one .xml file that will generate our two run cases using the file in your TCLB directory in tools/RT.

The domain that we are going to use for the RTI simlation will be $L\times 4L$, with $L=256$. The interface between the two fluids is initially set at half way, but then perturbed by a cosine wave with an applitude of $0.1L$. The following parameters are used to specify the case, with the values given representing the {low density, high density} ratio scenarios:

$At = \frac{\rho_H - \rho_L}{\rho_H+\rho_L} = {0.5, 0.998}$ $t_0 = \sqrt{L/g At} = {16000,8000}$ $U_0 = \sqrt{gL}$ $Re = \frac{\rho_H*U_0*L}{\mu_H} = {3000, 3000}$ $Ca = \frac{\mu_H*U_0}{\sigma} = {0.26, 0.44}$ $Pe = \frac{U_0 L}{M}$

These are the Atwood number, reference velocity, Reynolds number, Capillary number, Peclet number and the reference time, respectively. With these, we can rearrange to find all the required simulation parameters for the simulation. This fact is highlighted as they are used as inputs to our template file.

Creating the template run file

Create a new file called 'RTI_template.xml.Rt' and to start with, we are going to define the parameters above (for simplicity, we have specified the density and viscosity ratios in place of using the Atwood number):

<?R
L=c(256, 256)

rhostar = c(3,1000)
mustar  = c(3,100)

At = c(0.5,    0.998)
Re = c(3000,   3000)
Ca = c(0.26,   0.44) 
Pe = c(1000,   1000)
to = c(16000,  8000)

With this, we can now start calculating the simulation parameters:

W = c(5,       5)

rho = c(3,     1)
rhol= rho / rhostar

g  = L/(to^2*At)
Uo = sqrt(L*g)
muh= rho*Uo*L/Re
mul= muh/mustar
nuh=muh/rho
nul=mul/(rho/rhostar)

sigma = muh*Uo/Ca
M = Uo*L / Pe

The final parameters that we will set before building the actual TCLB run file are the names of the files we want and their dimensionless runtime:

names = c('RTI_LDR.xml', 'RTI_HDR.xml')
time = c(3, 2)

This effectively formulates our preamble, and we now can loop through the file names creating the necessary files:

for (i in c(1,2)){
    print(names[i])
    sink(names[i]) ?>
<?xml version="1.0"?>
<!--Model:  d2q9_pf_velocity 
    Created:    17-01-2020 
    By:     T.Mitchell -->
<CLBConfig version="2.0" output="output/" permissive="true">
    <Geometry nx="<?%f L[i] ?>" ny="<?%f 4*L[i] ?>">
        <MRT>
            <Box/>
        </MRT>
        <SpikeTrack>
            <Box dx="<?%f 0.5*L[i] ?>" nx="1" dy="10" fy="<?%f 4*L[i] - 10 ?>"/>
        </SpikeTrack>
        <BubbleTrack>
            <Box nx="1" dy="10" fy="<?%f 4*L[i] - 10 ?>"/>
        </BubbleTrack>

        <Wall mask="ALL">
            <Box dy="-1"/>
        </Wall>
        <Wall mask="ALL">
            <Box ny="1"/>
        </Wall>
    </Geometry>
    <Model>
        <Param name="MidPoint" value="<?%f 2*L[i] ?>"/>
        <Param name="Perturbation" value="0.1"/>
        <Param name="Period" value="<?%f L[i] ?>"/>
        <Param name="sigma" value="<?%f sigma[i] ?>"/>
        <Param name="M" value="<?%f M[i] ?>"/>
        <Param name="W" value="<?%f W[i] ?>"/>
        <Param name="GravitationY" value="<?%.8f -1*g[i] ?>"/>
        <Param name="Viscosity_h" value="<?%.8f nuh[i] ?>"/>
        <Param name="Viscosity_l" value="<?%.8f nul[i] ?>"/>
        <Param name="PhaseField_init" value="1.0"/>
        <Param name="Density_l" value="<?%f rhol[i] ?>"/>
        <Param name="Density_h" value="<?%f rho[i]  ?>"/>
    </Model>
    <VTK/>
    <Solve Iterations="<?%f time[i] * to[i] ?>" output="output/">
        <VTK Iterations="<?%f  0.5*to[i] ?>"/>
        <Log Iterations="<?%f 0.01*to[i] ?>" />
    </Solve>
</CLBConfig>
<?R
    sink()
}
?>

Here we can see that everything is being sunk into the file names that we originally created and the characters <? ?> are being used to reference the calculated values for both the low and high density ratio cases. It is also evident that there are built in parameters for perturbing the interface of the two fluids, one could also create this shape in geometry and assign zones for the fluid parameters. Another important point in the geometry of this model is that the ADDITIONAL node types SpikeTrack and BubbleTrack have been specified. Here, two vertical lines are created through which one tracks the front of the rising low density fluid and the other tracks the downward propagation of the heavy fluid. This is using the global option 'op="MAX"' in the Dynamics.R file for the model. As such, to ensure the lowest point of the high density ratio fluid is tracked, the location is stored as $4\times L - Y_{interface}.$

Running the simulation

With our template file now created, run the command:

tools/RT -f RTI_template.xml.Rt

This will generate two new .xml files ready for use with the TCLB solver.

As mentioned, we will be using the d2q9_pf_velocity model as such run from your TCLB directory:

make -j 8 d2q9_pf_velocity

Now we can run the two simulations:

CLB/d2q9_pf_velocity/main RTI_LDR.xml
CLB/d2q9_pf_velocity/main RTI_HDR.xml

Post-processing

The run files here will create a series of vtk files in the output folder which you can view and enjoy the complex dynamics of the Rayleigh-Taylor instability. If you would like to take this tutorial a step further, comparing the interface location of the low density bubble and high density spike with the referenced paper is a worthwhile activity. To visualise the progression of the interface in the simulations you have just run, the Log files in the output can be used. These can be opened in Excel or analysed with your favourite language, for example if we use R:

myData = read.csv("output/RTI_LDR_Log_00000000.csv")
t0 = 16000
L  = 256
plot(myData$Iteration/t0, myData$RTIBubble/L, ylim=c(0,4))
points(myData$Iterations/t0, 4-myData$RTISpike/L)

This has hopefully given you a nice introduction to the multiphase solver in TCLB as well as using RTemplate to generate multiple run files, this is particularly useful for parameter sweeps and keeping track of changes in simulations. If you have any troubles with this tutorial, feel free to contact the development team!

Tip

A good starting point explaining the 'zoology' of various LBM models is a modern (2017) book 'The Lattice Boltzmann Method: Principles and Practice' written by T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E.M. Viggen .

LINK TO THE GIT PAGE FOR THE GEOMETRY ↩