# TLBfind: a Thermal Lattice Boltzmann code for concentrated emulsions with
FINite-size Droplets

Francesca Pelusi  Matteo Lulli  Mauro Sbragaglia  Massimo Bernaschi  Helmholtz
Institute Erlangen-Nürnberg for Renewable Energy, Forschungszentrum Jülich,
91058 Erlangen, Germany  Department of Mechanics and Aerospace Engineering,
Southern University of Science and Technology, Shenzhen, 518055 Guangdong,
China  Department of Physics & INFN, University of Rome “Tor Vergata”, 00133
Rome, Italy  Istituto per le Applicazioni del Calcolo (IAC) - CNR, 00185 Rome,
Italy

###### Abstract

In this paper, we present TLBfind, a GPU code for simulating the hydrodynamics
of droplets along with a dynamic temperature field. TLBfind hinges on a two-
dimensional multi-component lattice Boltzmann (LB) model simulating a
concentrated emulsion with finite-size droplets evolving in a thermal
convective state, just above the transition from conduction to convection. The
droplet concentration of the emulsion system is tunable and at the core of the
code lies the possibility to measure a large number of physical observables
characterising the flow and droplets. Furthermore, TLBfind includes a parallel
implementation on GPU of the Delaunay triangulation useful for the detection
of droplets’ plastic rearrangements, and several types of boundary conditions,
supporting simulations of channels with structured rough walls.

###### keywords:

lattice Boltzmann; soft suspensions; finite-size droplets; thermal convection;
rough channels.

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: TLBfind
CPC Library link to program files: (to be added by Technical Editor)
Developer’s repository link: https://github.com/FrancescaPelusi/TLBfind
Code Ocean capsule: (to be added by Technical Editor)
Licensing provisions: MIT
Programming language: CUDA-C
Nature of problem: Hydrodynamics of concentrated emulsions with finite-size
droplets in a thermal convective state.
Solution method: Single relaxation time Lattice Boltzmann (LB) method to solve
Navier-Stokes equations for fluids, coupled with the temperature field
dynamics . The output describes the dynamics of finite-size droplets of
concentrated emulsions in presence of a temperature field. The temperature
field obeys the advection-diffusion equation.
Additional comments including restrictions and unusual features: Plastic
rearrangements of droplets are detected via the parallel implementation of the
Delaunay triangulation, and boundary conditions are tunable.

##  1 Introduction

Understanding the hydrodynamic behaviour of concentrated emulsions – and in
general soft particles suspensions – represents an intriguing subject of study
in the context of fluid dynamics, with a wide range of applications, from
everyday life situations to modern technologies [1, 2, 3, 4, 5, 6, 7].
Although many questions have been answered in the last decades [8, 9, 10, 11,
12], many other aspects still deserve further scrutiny, such as the precise
understanding of the heat transfer properties in these materials when evolving
in a thermal convective state. This may be relevant in a variety of situations
including the convective motion of magma [13] composed of a melt with crystal
suspensions where convection may be limited to the hotter (less crystalline)
portions, the oil recovery industrial processes [14] aiming at using emulsions
as packer fluids in order to significantly reduce the heat transfer rate
during extraction, the production of slurry ice [15], i.e., ice crystals
distributed in water or an aqueous solution typically used to replace ice and
salt for food cooling. The need for numerical simulations allowing to study
the convective heat transfer in model emulsions motivated the present work.
Emulsions are structurally characterised by a collection of “soft domains”
(i.e., droplets) of a dispersed phase in another continuous phase. The
emulsion droplets concentration influences the rheological properties of the
system: dilute emulsions behave as Newtonian fluids, with the viscosity being
constant regardless of the applied shear rate; more concentrated emulsions
display a non-Newtonian mechanical response and exhibit a viscosity that
depends non-linearly on the applied shear rate [16]. Furthermore, there exists
a critical concentration above which emulsions can be categorised as yield-
stress fluids [10, 17, 18], wherein the system shows a stress threshold (the
so-called yield stress) below which no flow is observed and above which the
emulsion flows displaying shear-thinning. Emulsions have been characterised in
the literature via many experimental studies (see [19, 20, 21] for reviews).
However, experimental insights may be arduous to capture at scales comparable
with the droplet size. This issue warrants the use of numerical simulations to
investigate the response of emulsions and, more generally, soft particle
suspensions. For the latter, different numerical methods have been exploited,
such as the boundary integral method [22], the discrete element method [23,
24], molecular dynamics simulations [25, 26, 27], and lattice Boltzmann (LB)
models [28, 29, 30] just to cite some examples. The focus of this paper is on
LB models, which have become very popular and attractive in the last two
decades, due to their simplicity, efficiency and applicability in different
contexts [31, 32]. Several open-access LB implementations are available, e.g.,
Ludwig [33], LB3D [34], LBsoft [35], Palabos [36], and LBfoam [37]. However,
an LB software that simulates the dynamics of finite-size droplets stabilised
against coalescence and subjected to thermal convection is not available to
the best of our knowledge in an open-access version. We aim at filling this
gap by presenting TLBfind, a code based on a two-dimensional LB scheme devised
to simulate emulsions with finite-size droplets where the temperature field
dynamics is coupled with the emulsion momentum equation. We consider an
emulsion confined between two parallel walls, heated from below and cooled
from above, i.e., in the paradigmatic set-up of the Rayleigh-Bénard convection
[38, 39, 40, 41]. Studies on the thermal response of this kind of emulsions
have been presented in a recent work [42]. TLBfind was tested to be optimal to
study the heat transfer properties in the regime where the system sustains a
convective state, just above the transition from conduction to convection. In
TLBfind it is possible to tune the droplet concentration, allowing to
systematically transition from diluted to concentrated emulsions. Further
flexibility of TLBfind lies in a variety of different boundary conditions that
can be explored, and the high degree of parallelisation obtained on GPU using
the CUDA-C language entails a huge efficiency, by saving computational costs
and time.
The paper is organised as follows: in Section 2 we briefly explain the LB
scheme implemented in TLBfind. In Section 3 we give a test case of the
Rayleigh-Bénard convection in concentrated emulsion when the walls of the
channel are flat; a test case with rough walls will be analysed in Section 4,
then we summarise the potentialities of TLBfind in Section 5.

##  2 Method

We leverage an LB method [31, 32] for non-ideal multi-component systems that
allows the mesoscopic simulation of a collection of finite-size droplets
evolving in a thermal convective state. The isothermal counterpart of the
model has already been presented and validated in many publications [43, 44,
45, 28, 46, 47, 48]. TLBfind considers the advection-diffusion dynamics of a
scalar temperature together with the multi-component model. We summarise the
main essential features of the methodology, by treating separately the
isothermal multi-component fluid and the temperature dynamics.

###  2.1 Multi-component lattice Boltzmann

We consider a two-dimensional system consisting of two components (labelled by
$\ell=1,2$). The LB dynamics considers mesoscopic probability density
functions $f_{\ell,i}({\bf x},t)$, representing the fluid particle (mass)
density of the component $\ell$ at the space-time location $({\bf x},t)$ and
velocity ${\bf c}_{i}$, where the position $\mathbf{x}$ takes on the integer
values of the nodes coordinates in a squared lattice. The $f_{\ell,i}({\bf
x},t)$’s evolve via the discrete lattice Boltzmann equations:

$f_{\ell,i}({\bf x}+{\bf c}_{i},t+1)-f_{\ell,i}({\bf
x},t)=\Omega_{\ell,i}({\bf x},t),$ (1)

where both the lattice spacing $\Delta x$ and time step $\Delta t$ are
considered equal to the unity. The r.h.s. of Eq. (1) encodes the physical
effects of collisions, i.e., to redistribute particles among the populations
$f_{\ell,i}({\bf x},t)$ at each node. Collisions are implemented via the
Bhatnagar–Gross–Krook (BGK) operator [49], approximating the relaxation of
$f_{\ell,i}({\bf x},t)$ towards the local equilibrium distribution
$f_{\ell,i}^{(eq)}$ with a relaxation time $\tau$:

$\Omega_{\ell,i}({\bf x},t)=-\frac{1}{\tau}\left[f_{\ell,i}({\bf
x},t)-f_{\ell,i}^{(eq)}\left(\rho_{\ell}({\bf x},t),{\bf u}_{\ell}({\bf
x},t)\right)\right].$ (2)

The local equilibrium $f_{\ell,i}^{(eq)}$ is the Maxwellian distribution,
given in the form

$f_{\ell,i}^{(eq)}\left(\rho_{\ell},{\bf
u}_{\ell}\right)=w_{i}\rho_{\ell}\left[1+\frac{u_{\ell,k}c_{ik}}{c_{s}^{2}}+\frac{u_{\ell,k}u_{\ell,p}(c_{ik}c_{ip}-c_{s}^{2}\delta_{kp})}{2c_{s}^{4}}\right],$
(3)

where the $w_{i}$ are the lattice weights, $c_{s}^{2}=\sum_{i}w_{i}|{\bf
c}_{i}|^{2}/d=1/3$ is the squared sound velocity (with $d=2$ the space
dimension), and repeated indices are summed upon. The index $i$ runs on a
finite number of values depending on the LB model used. In TLBfind, we
leverage the widely used D2Q9 scheme (see Fig. 1(a), orange arrows), with $9$
lattice velocities (i.e., $i=0,\ldots,8$) in a two-dimensional domain. The
associated weights $w_{i}$ are reported in Fig. 1(b). After the collision, the
l.h.s. of Eq. (1) performs a streaming of each fluid particle (mass) density
on the lattice following the set of lattice velocities ${\bf c}_{i}$.
Information on coarse-grained fields of densities of each component
($\rho_{\ell}$) and emulsion velocity (${\bf u}$) can be extracted by by
taking moments in the discrete velocity space as

$\rho_{\ell}({\bf x},t)=\sum_{i}f_{\ell,i}({\bf x},t),$ (4) ${\bf u}({\bf
x},t)=\frac{1}{\rho}\sum_{\ell,i}{\bf c}_{i}f_{\ell,i}({\bf x},t),$ (5)

where we used the total density defined as $\rho=\sum_{\ell}\rho_{\ell}$.

---
(a)
| $w_{i}=w(\bf{c}_{i})$ | $p_{i}=p(\bf{c}_{i})$ | $|\bf{c}_{i}|$ | $\bf{c}_{i}$ | $i$
---|---|---|---|---
4/9 | 247/420 | 0 | (0,0) | 0
1/9 | 4/63 | 1 |  $(\pm 1,0)$;$(0,\pm 1)$ | 1-4
1/36 | 4/135 | 2 | $(\pm 1,\pm 1)$ | 5-8
0 | 1/180 | 4 |  $(\pm 2,0)$ ; $(0,\pm 2)$ | 9-12
0 | 2/945 | 5 |  $(\pm 1,\pm 2)$; $(\pm 2;\pm 1)$ | 13-20
0 | 1/15120 | 8 | $(\pm 2,\pm 2)$ | 21-24
(b)
Figure 1: Panel (a): a scheme of the two neighbours sets of discrete
velocities used for the evaluation of the interaction forces: a “nearest-
neighbours” (NN) and a “next-to-nearest neighbours” (NNN) zone. The NN zone
corresponds to the set of velocities ${\bf c}_{i}$ of the D2Q9 LB method used,
with index $i=0...8$. NNN are links with index $i=9...24$. Concerning the
interactions, each fluid component self-interacts via competing interactions,
i.e., attractive forces involving the NN (Eq. (10)), and repulsive ones (Eq.
(11)) acting on both zones. On the other hand, the fluid-fluid cohesion
interactions Eq. (7) are given by summing on sites in the NN zone. Panel (b)
shows the list of weights appearing in Eqs. (7), (10) and (11).

The effects of interaction forces, ${\bf F}^{\mbox{\tiny int}}_{\ell}({\bf
x},t)$, and external volume forces, ${\bf F}^{\mbox{\tiny ext}}_{\ell}({\bf
x},t)$ (the latter will be discussed in Section 2.3), are encompassed in a
source term ${\bf F}_{\ell}({\bf x},t)={\bf F}^{\mbox{\tiny int}}_{\ell}({\bf
x},t)+{\bf F}^{\mbox{\tiny ext}}_{\ell}({\bf x},t)$ that enters in Eq. (2) as
a shift in the hydrodynamic velocity (5):

${\bf u}_{\ell}({\bf x},t)={\bf u}({\bf x},t)+\dfrac{\tau{\bf F}_{\ell}({\bf
x},t)}{\rho_{\ell}}.$ (6)

The interaction forces ${\bf F}^{\mbox{\tiny int}}_{\ell}({\bf x},t)$ hinge on
the multi-range Shan-Chen methodologies [50, 51, 52, 53] and include three
different contributions for the case presented in this paper [43, 54]: ${\bf
F}^{x}_{\ell}({\bf x},t)$, ${\bf F}^{a}_{\ell}({\bf x},t)$, and ${\bf
F}^{r}_{\ell}({\bf x},t)$. The term ${\bf F}^{x}_{\ell}({\bf x},t)$ represents
the interactions between fluid elements of the two different species; it is
introduced with the aim of creating phase segregation and the formation of
stable diffuse interfaces, i.e., interfaces characterised by a finite width.
In formulae, it reads:

${\bf F}^{x}_{\ell}({\bf x},t)=-{\cal G}_{12}\psi_{\ell}({\bf
x},t)\sum_{\ell^{\prime},\ell^{\prime}\neq\ell}\sum_{i\in\text{NN}}w_{i}\psi_{\ell^{\prime}}({\bf
x}+{\bf c}_{i},t){\bf c}_{i},$ (7)

where $\psi_{\ell}({\bf x},t)=\psi(\rho_{\ell}({\bf x},t))$ is the pseudo-
potential function of the Shan-Chen formulation [50], and ${\cal G}_{12}$ is a
positive coupling constant dictating the strength of these interactions
111${\cal G}_{12}$ also controls the width of the interface.. In the specific
case, we use the simplest pseudo-potential, i.e.,

$\psi_{\ell}({\bf x},t)=\dfrac{\rho_{\ell}}{\rho_{0}}$ (8)

where $\rho_{0}$ is a reference density. The set NN refers to the set of
“nearest neighbours” nodes of ${\bf x}$ on the lattice, which coincide with
the set of the D2Q9 directions used for the streaming of the LB populations
(see Fig. 1(a)). The effect of ${\bf F}^{x}_{\ell}({\bf x},t)$ is a change in
the bulk pressure, which now features ideal contributions summed to the non-
ideal ones:

$P_{b}(\rho_{1},\rho_{2})=\underbrace{c_{s}^{2}\rho_{1}+c_{s}^{2}\rho_{2}}_{\text{ideal}}+\underbrace{c_{s}^{2}{\cal
G}_{12}\rho_{1}\rho_{2}}_{\text{non-ideal}}.$ (9)

The two other contributions to ${\bf F}_{\ell}^{\mbox{\tiny int}}$, i.e.,
${\bf F}^{a}_{\ell}({\bf x},t)$ and ${\bf F}^{r}_{\ell}({\bf x},t)$, represent
short-range attractive ($a$) and long-range repulsive ($r$) competing
interactions, respectively [43, 46]. They are introduced to simulate the
action of a positive disjoining pressure at the interface that inhibits the
coalescence of droplets [44] (see Fig. 3). In formulae, they read [55]:

$\displaystyle{\bf F}^{a}_{\ell}({\bf x},t)=$ $\displaystyle-{\cal
G}^{a}_{\ell\ell}\psi_{\ell}({\bf
x},t)\sum_{i\in\text{NN}}w_{i}\psi_{\ell}({\bf x}+{\bf c}_{i},t){\bf c}_{i}$
(10) $\displaystyle{\bf F}^{r}_{\ell}({\bf x},t)=$ $\displaystyle-{\cal
G}^{r}_{\ell\ell}\psi_{\ell}({\bf
x},t)\left[\sum_{i\in\text{NN}}p_{i}\psi_{\ell}({\bf x}+{\bf c}_{i},t){\bf
c}_{i}+\sum_{i\in\text{NNN}}p_{i}\psi_{\ell}({\bf x}+{\bf c}_{i},t){\bf
c}_{i}\right],$ (11)

with ${\cal G}_{\ell\ell}^{a}<0$ and ${\cal G}_{\ell\ell}^{r}>0$. In Eqs.
(10)-(11), the pseudo-potential is used in the original Shan-Chen form [50]:

$\psi_{\ell}({\bf x},t)=\rho_{0}\left[1-\exp(-\rho_{\ell}({\bf
x},t)/\rho_{0})\right].$ (12)

The set of nodes NNN in (11) refers to “next-to-nearest neighbours”, i.e., an
additional layer of 16 lattice velocities beyond the D2Q9 links (see Fig.
1(a), light-blue arrows). The values of the weights $p_{i}$ are reported in
Fig. 1(b).
By summing over the components $\ell$, the total force ${\bf F}({\bf
x},t)=\sum_{\ell}{\bf F}_{\ell}$ acting on the momentum density (Eq. (6)) is
obtained.
The reference hydrodynamic equations at large scales are the diffuse-interface
Navier-Stokes equations:

$\rho\left(\partial_{t}+u^{\mbox{\tiny(H)}}_{k}\partial_{k}\right)u^{\mbox{\tiny(H)}}_{i}=-\partial_{j}P_{ij}+\eta_{0}\partial_{j}\left(\partial_{i}u^{\mbox{\tiny(H)}}_{j}+\partial_{j}u^{\mbox{\tiny(H)}}_{i}\right)+f^{\mbox{\tiny
ext}}_{i},$ (13)

where $\rho{\bf u}^{\mbox{\tiny(H)}}=\rho{\bf u}+{\bf F}/2$ is the
hydrodynamical momentum density, and $i=x,y$. Notice that the hydrodynamic
velocity ${\bf u}^{\mbox{\tiny(H)}}$ differs from Eq. (5) because of forcing
renormalisations. The bare viscosity $\eta_{0}$ in the hydrodynamic equations
is related to the relaxation time $\tau$ of the LB equation (1) as:

$\eta_{0}=\rho c_{s}^{2}\left(\tau-\frac{1}{2}\right),$ (14)

hence it can be tuned by changing the $\tau$ in the LB dynamics. The viscosity
$\eta_{0}$ represents the “bare” viscosity, thus the bare viscous stress
tensor
$\eta_{0}\left(\partial_{i}u^{\mbox{\tiny(H)}}_{j}+\partial_{j}u^{\mbox{\tiny(H)}}_{i}\right)$
is supplemented by the pressure tensor $P_{ij}$ that depends on the density
heterogeneities [28]. These two contributions of stress sum up to give the
total stress that is used in the rheological characterisations of the
emulsions (see Section 3.2.2). Finally, the term $f^{\mbox{\tiny ext}}_{i}$ in
Eq. (13) appears on behalf of the density of external forces (the buoyancy
force density will be discussed in Section 2.3).

###  2.2 Boundary conditions and rough walls

TLBfind allows for the exploration of both fully periodic and confined
systems. Rough-wall-flags – one associated with each wall – are introduced
with the aim of switching on/off the presence of roughness at the walls (these
features will be discussed in detail in Sections 3 and 4). If the rough-wall-
flags are switched off, we are in the presence of flat walls and the code
implements a modified mass-conserving bounce-back rule [32, 31], that is a
bounce-back dynamics designed to assign the desired input value to the
hydrodynamical velocity at the walls

${\bf u}^{\mbox{\tiny(H)}}_{x}({\bf x}_{b},t)=u_{\mbox{\tiny w}}{\bf e}_{x}$
(15)

where ${\bf x}_{b}$ are the boundary nodes. This set-up in TLBfind is used
when performing rheological experiments (see Section 3.2.2) in a Couette
channel, where each wall moves with its own velocity $u_{\mbox{\tiny w}}{\bf
e}_{x}$ along the stream-flow direction $x$. In the presence of structured
rough walls, rough-wall-flags are switched on and TLBfind operates with simple
half-way bounce-back boundary conditions, with $u_{\mbox{\tiny w}}{\bf
e}_{x}=0$.

####  2.2.1 Wetting conditions

In the presence of walls, both flat and rough, TLBfind allows handling the
fluid-wall interactions, i.e., the wetting conditions. The implementation is
very simple, regarding just the definition of density values of the two
components at the wall: the code sets these densities (and the pseudo-
potential accordingly) in the ghost nodes along the $y$ direction ($y=0,ny+1$)
222The ghost nodes are introduced along both $x$ and $y$ directions as part of
the lattice to facilitate the collision and streaming step close to the
boundaries with the aim to supply the boundary nodes with otherwise missing
populations. and in the wall nodes as equal to the values $\rho_{\mbox{\tiny
w, max}}$ and $\rho_{\mbox{\tiny w, min}}$, respectively. More details on how
to perform simulations with different wetting conditions are provided in
Section 3.

###  2.3 Thermal lattice Boltzmann

An auxiliary probability distribution function $g_{i}({\bf x},t)$ is
introduced with the purpose of simulating the dynamics of a scalar temperature
field $T({\bf x},t)$. Similarly to the fluid populations $f_{\ell,i}$,
$g_{i}({\bf x},t)$ are governed by a discrete lattice Boltzmann equation [56,
32, 31]

$g_{i}({\bf x}+{\bf c}_{i},t+1)-g_{i}({\bf
x},t)=-\frac{1}{\tau_{g}}\left[g_{i}({\bf x},t)-g_{i}^{(eq)}\left(T({\bf
x},t),{\bf u}^{\mbox{\tiny(H)}}({\bf x},t)\right)\right],$ (16)

and the temperature field is obtained by taking the moment of order zero of
the distribution functions:

$T({\bf x},t)=\sum_{i}g_{i}({\bf x},t).$ (17)

The local equilibrium $g_{i}^{(eq)}$ for the temperature field takes the form
333The temperature field $T({\bf x},t)$ has to be interpreted as the relative
temperature with respect to some reference value.

$g_{i}^{(eq)}(T,{\bf
u}^{\mbox{\tiny(H)}})=w_{i}T\left[1+\frac{u^{\mbox{\tiny(H)}}_{k}c_{ik}}{c_{s}^{2}}+\frac{u^{\mbox{\tiny(H)}}_{k}u^{\mbox{\tiny(H)}}_{p}(c_{ik}c_{ip}-c_{s}^{2}\delta_{kp})}{2c_{s}^{4}}\right].$
(18)

The long-wavelength limit of (16) approximates the advection-diffusion
equation for the temperature field

$\partial_{t}T+u^{\mbox{\tiny(H)}}_{k}\partial_{k}T=\kappa\partial_{kk}T,$
(19)

where the thermal diffusivity $\kappa$ is related to the thermal relaxation
time $\tau_{g}$ as follows:

$\kappa=c_{s}^{2}\left(\tau_{g}-\frac{1}{2}\right).$ (20)

The buoyancy term in Eq. (13) is defined as:

$f^{(\mbox{\tiny ext})}_{i}=\rho\alpha gT\delta_{iy},$ (21)

with $\alpha$ being the thermal expansion coefficient, $g$ the gravity
acceleration, and ${\bf e}_{y}$ the unit vector in the wall-to-wall direction.
With the aim to obtain the buoyancy force in (21), we include in Eq. (6) an
external volume force that reads

${\bf F}_{\ell}^{\mbox{\tiny ext}}({\bf x},t)=\rho_{\ell}({\bf
x},t)\,\alpha\,g\,T({\bf x},t)\,{\bf e}_{y}.$ (22)

###  2.4 Delaunay triangulation

---
(a) (b)
Figure 2: Domains topology analysis of a two-dimensional concentrated
emulsion. Dark-yellow regions refer to emulsion droplets and blue regions are
occupied by the continuous phase/component. Panel (a): construction of
Delaunay triangulation (black lines) and Voronoi tassellation (orange lines)
based on the centres-of-mass positions (black dots). Panel (b): a sketch of
detection of a T1 plastic rearrangement. The old (dashed red) link between
droplets $a$ and $b$ goes to zero and the new (blue continuous) link is
created.

TLBfind performs a run-time analysis of Delaunay triangulation [57] on the
centres of mass of the droplets. The main reason for executing the Delaunay
triangulation on GPU is the performance improvement we obtain. In a very
preliminary version of TLBfind we resorted to an existing CPU library to carry
out the triangulation. However, that choice had two disadvantages: $i)$ we had
to copy back data from the GPU to CPU and it is well know that the bandwidth
between CPU and GPU is limited so that the copy is slow; $ii)$ the CPU library
we used and other open-source alternatives we found were implemented without
considering high-performance requirements. So we decided to implement our own,
high-performance parallel implementation. The Delaunay triangulation can be
seen as an adjacency matrix of the droplets, storing information about the
topology of the droplets arrangement. In particular, the triangulation is
obtained in such a way that the circumference circumscribed at any of the
triangles does not contain any other point of the set. A reference depiction
is shown in Fig. 2. The algorithm implemented in TLBfind was described in
detail in [58], and here we report the main ingredients. The first analysis
step consists in identifying the bulk of the droplets as those regions of the
lattice where the density field of the first component is larger than a
certain threshold $B_{th}$. Then, these sets of lattice points are identified
by means of a GPU implementation of a clustering algorithm adapted from those
used in accelerated Ising model simulations at the critical point [59, 60]. We
remark that, while in the case of the Ising model the construction of the
cluster is probabilistic in nature, in our present case we use the algorithm
in a _geometric_ way in order to identify all the nodes related to the bulk of
a droplet. This allows us assigning unique labels to all droplets and compute
their centre of mass (black dots in Fig. 2$(a)$). The labels associated with
each droplet are stored and employed to define a _digital Voronoi_
tessellation [61] where each point of the lattice is “coloured” by the label
of the nearest centre of mass (orange lines). In the continuum construction, a
Voronoi tessellation is dual to a Delaunay triangulation, thus encoding the
same information but in a metric language. The vertices at which adjacent
digital Voronoi cells converge are used to identify the Delaunay triangulation
(black links in Fig. 2$(a)$). Subsequent Delaunay triangulations are compared
in order to detect T1 plastic rearrangements of the droplets, which have the
effect of stably changing the topology of the droplets arrangement [58]. These
topological events only involve 4 droplets, which can be thought of as being
at the four vertices of a quadrilateral labelled $a$, $b$, $c$ and $d$, and
they are associated with a ”flip” of the diagonal link of the triangulation.
For example, in Fig. 2$(b)$, the diagonal link between $a$ and $b$ (dashed red
line), is replaced by the link between $c$ and $d$ (solid blue line). From the
Voronoi tesselation perspective, at the degenerate configuration, i.e., when
the four droplets all lie on a common circle, the Voronoi edge orthogonal to
the Delaunay side shrinks to zero. This is an equivalent way, though less
robust, to implement the T1 event detection. Avoiding dealing with such a
metric structure makes the detection of plastic events via Delaunay
triangulation more robust and easier to implement.

###  2.5 GPU Implementation

One of the crucial requirements to achieve a good performance on the NVIDIA
GPU is that global memory accesses (both read and write) should be coalesced.
This means that memory access needs to be aligned and coordinated within a
group of threads. The basic rule is that the thread with id
$n\in\\{0,\ldots,N-1\\}$ should access element $n$ at byte address
$StartingAddress+sizeof(type)\cdot n$ where $sizeof(type)$ is equal to either
$4$, $8$ or $16$ and $StartingAddress$ is a multiple of $16\cdot
sizeof(type)$.

With this aim, the fluid populations of a lattice site are not contiguous in
the GPU global memory so that they are ordered following the structure-of-
arrays data layout. All data not modified during the simulation, such as
coefficients, are pre-computed during the initialisation phase and stored in
the GPU constant memory, which has performances analogous to those of
registers if, as in our case, all the threads running on the same
multiprocessor access the same constant memory locations.

After the initialisation phase, all the computations required for the LB
update are performed on the GPU. A single step of the simulation is
implemented through a sequence of CUDA kernels guaranteeing the correct
sequential order of the sub-steps. Each CUDA kernel implements a sub-step of
the update procedure (e.g., collision, streaming) by splitting the work among
a configurable number of threads and blocks, which may be fine-tuned to
achieve optimal performance on different CUDA devices. Each thread works
sequentially on a group of lattice nodes assigned to it. For each lattice
node, the thread copies data from the global memory into registers, performs
the computation and writes the results back in the global memory. In order to
manage the parallelisation of the streaming phase without causing conflicts
among multiple threads, fluid populations are stored in the global memory
using a “double buffer” policy. At the end of the simulation, the final
results are copied back to the CPU main memory in order to be saved on a
secondary storage device. Through the input file it is possible to require
also the saving of partial results of the simulation at regular intervals
(e.g., for check-pointing purposes).

Most of the global memory read-and-write operations are coalesced, with the
exception of a few reads relative to the computation of the interaction forces
and a few writes relative to the streaming phase. In the first case, the
calculation of the force for a lattice node depends on values related to other
lattice nodes, which must be loaded from global memory even if alignment
requirements for coalesced accesses are not satisfied. In the second case,
target locations of the streaming phase are defined by the lattice topology,
and in general, they don’t comply with the memory alignment requirements.

For the function that computes the value of the hydrodynamic variables all
memory operations are “local”, meaning that only the fluid populations of a
lattice site are required and that the resulting hydrodynamic variables belong
to the same lattice site. As a consequence, there is not a single uncoalesced
memory access.

Finally, as already mentioned, the fluid populations once uploaded on the GPU
memory do not need to be copied back to the main memory unless a dump of the
whole configuration is required. However, hydrodynamic variables or other
observables derived from them might be written back to the main memory much
more frequently since they represent the main physical output of the
simulation. Although the number of hydrodynamic variables per lattice site is
small compared to the number of fluid populations (there are 4 hydrodynamic
variables vs. 9 fluid populations), so that the run-time overhead of the copy
from the GPU-memory to the CPU-memory is small compared to the initialisation
overhead, better performance may be obtained by reducing the number of these
copy-back operations. Although the main goal of the present paper is to make
TLBfind available and explain how to use it, we report some basic data about
its computational performance. It is common to measure the performance of an
LB code by using the LUPS (Lattice Updates per Second) metrics. TLBfind has,
at least, two sets of populations for each fluid node and this means that both
the amount of data to be moved from/to GPU global memory and the number of
operations should be (at last) doubled. The number of MLUPS (Million of LUPS)
that TLBfind is able to perform is $\sim 525$ on a Titan-V a GPU featuring
5120 CUDA cores (Volta architecture). This figure is perfectly in line with
the performance (515 MLUPS) of sailfish [62] an LB open-source code
(https://github.com/sailfish-team/sailfish.git) on the same hardware. By using
Nsight, a user-friendly Nvidia visual profiler, we collected several
additional information about the performance of TLBfind. In Table 1 we report
two of them: the warp execution efficiency that shows a very high degree of
utilisation and the throughput of the load global memory operations for three
representatives kernels. For the Titan V the peak memory bandwidth is $\sim
650$ GB/s, and TLBfind achieves $>80\%$ of that value. It is well known that
LB codes are memory bandwidth bound meaning that the number of floating-point
operations they execute is limited with respect to the amount of data that
need to be loaded and stored from/to memory. This is true for both CPU and GPU
codes and TLBfind is not an exception. As mentioned above, although TLBfind
has, potentially, up to three sets of populations per node, the ratio between
the number of arithmetic operations required by the execution of each LB time
step and the amount of data moved from/to memory is fundamentally the same of
any other “simple” 2D LB code. As to the Delaunay triangulation, we did not
carry out specific performance measurements. However, we noticed that the
Delaunay procedure takes a time that is, roughly, equivalent to two LB time
steps. Further details about the GPU implementation can be found in [58] and
[63].

CUDA performance metrics
---
kernel |  metrics |  average value
AVERAGE_1 |  Warp Execution Efficiency |  99.74%
AVERAGE_1 |  Global Load Throughput |  548 GB/s
moveplusforcingconstructWW |  Warp Execution Efficiency |  95.6% GB/s
moveplusforcingconstructWW |  Global Load Throughput |  524 GB/s
forcingconstructWW |  Warp Execution Efficiency |  99.85%GB/s
forcingconstructWW |  Global Load Throughput |  538 GB/s
Table 1: Some CUDA performance metrics obtained by using Nsight.

##  3 Test case with flat walls

As a first example, we show a simulation of Rayleigh-Bénard convection of a
model emulsion, by explaining in detail which are the parameters involved.
With this aim, hereafter in the text, all words in Computer Modern Typewriter
font are referred to as input parameters, and what will be defined as Boolean
is considered turned off (false = 0) or on (true = 1). The same font will be
used also for output files names. All tunable simulation parameters are listed
in the input file tlbfind.inp.
We remark that, before starting a simulation, TLBfind requires a preparation
step during which the structural properties of the system under study (e.g.,
droplet concentration) are defined. It is possible to reproduce all the
results from the examples discussed below by following the instructions
contained in README_howToPreparation and README_howToRun for the preparation
and simulations steps, respectively.

###  3.1 Output files options

In this paragraph, we describe the parameters that can be set in the input
file tlbfind.inp controlling the output files properties, such as dump
frequency and format type. All ASCII files including two-dimentional fields
can be plotted by using the pm3d command of gnuplot. The VTK files can be read
with programs like Paraview [64]. Table 2 reports typical values of the output
parameters.
Independently from the output options that are set, the program will output a
few files at the beginning of the preparation step, containing the information
about the initial state of the various fields describing the system. These
files are init_rho1.dat, init_rho2.dat, and init_temperature.dat, including
the initial values of the density fields of the two components and the
temperature field, respectively.

First of all we describe the options related to the density and velocity
fields of the two components. The main option is nout density which is an
integer setting the time interval for the ASCII dumps firstdensity.#.dat and
seconddensity.#.dat. The program also writes the files firstdensity.dat and
seconddensity.dat containing the last dumped values and which are overwritten
at each dump. If nout density is set to a value $\leq$ 0 no dump is performed.
The variable write vtk file (Boolean), acts as a primary switch changing the
dumpfile format from ASCII to binary VTK 2.0 files. In order to assure the VTK
dump one needs to toggle the (Boolean) variables write vtk file rho1 and write
vtk file rho2 which will make the program write on disk firstdensity.#.vtk and
seconddensity.#.vtk, respectively, with firstdensity.vtk and seconddensity.vtk
containing the last dumped configurations. The variable nout velocity sets the
time interval for the ASCII dump number of steps for the vector velocity
field, veloconf.#.dat, with the last configuration dumped in veloconf.dat. If
write vtk file is set to 1 the binary VTK 2.0 files veloconf.#.vtk are also
dumped.

The variable nout temperature is an integer setting the time interval for the
ASCII dumps temperature.#.dat, with temperature.dat containing the last
configuration. The variable write vtk file temperature (Boolean), toggles the
dumps in VTK 2.0 binary file temperature.#.vtk.
The variable write energy file (Boolean) is the main switch for writing the
values of the global kinetic energy $E=\langle|{\bf
u}|^{2}\rangle_{x,y}/(\mbox{{nx}}\cdot\mbox{{ny}})$ that are appended every
nout energy steps to the file timeEnergy.dat. The variable nout tensor sets
the number of steps for dumping the stress tensor ASCII files Ptot_xy.#.dat,
with Ptot_xy.dat containing the last configuration. The variable nout average
is an integer for the number of steps for dumping the average of the
$x$-component in the $x$-direction of (i) the velocity field in u_av.#.dat
files, and (ii) the stress tensor in Pxy_av.#.dat files. The last
configurations are dumped in u_av.dat and Pxy_av.dat files, respectively. More
details on the structure of the above-mentioned output files can be found in
the README_Output file.

Concerning the variables used for dumping the state of the system in order to
allow the simulation restart, noutconfig is an integer setting the dumping
number of steps for binary files conf1_#.in, conf2_#.in, and confG_#.in,
containing populations for the first component, second component, and
temperature, respectively. Furthermore, the number of successive populations
files is written in a sequence, modulo noutconfigMod, for fail-safe purposes.

Some variables in tlbfind.inp are referred to the detection of plastic events.
The integer ncheckdelaunay sets the number of steps at which the program
checks for T1 plastic events (see Section 2.4). The output files are written
in the sub-directories ./delaunayTriggerDir, collecting the dumps related to
plastic events, and ./delaunayNoTriggerDir, when no plastic event occurs.
Within these folders, the Delaunay binary files delaunayNowTime#,
delaunayPastTime#, and delaunayIsoTriggerTime# are dumped, containing
information about the position of the centre-of-mass of all droplets, and the
links connecting nearest droplets. By inspecting the functions contained in
delaunayCuda.cu, one can find the details about the structure of those binary
files. Furthermore, the trigger for the detection of plastic events also dumps
a few ASCII files in the Delaunay directories that can be used for a
preliminary run-time analysis of the simulation itself. In the files
arisingLinks and arisingLinksBoundary the user can find data of the
triangulation links that are formed after a plastic event. On the other hand,
breakingLinks and breakingLinksBoundary files store data about the links
vanishing after a plastic event is detected. Finally, the time sequence of the
total number of droplets and the total number of links in the Delaunay
triangulation are contained in nBubblesOut and nLinksOut, respectively. A
detailed description of these files is contained in the README_Output file.

###  3.2 Preparation

All simulation parameters in tlbfind.inp used in the preparation step are
listed in Table 3. At this step, no driving force or coupling with the
temperature field has to be considered. A step-by-step description of commands
to run a preparation-step simulation can be found in the
README_howToPreparation file.

####  3.2.1 Two-component fluid

The two-dimensional system is resolved with a regular grid with size nx along
the stream-flow direction $x$ and ny along the vertical $y$-direction. As
discussed in Section 2.5, TLBfind is parallelised on GPU threads via the CUDA
toolkit thus, once the simulation domain is fixed, it is extremely important
to declare in the input file Which GPU to use 444It depends on the machine
used., the number of used threads (nthreads) per block, and the number of
thread blocks (nblocks). Make sure that the product of nthreads with nblocks
is equal to the whole simulation domain nx$\cdot$ny 555We remark that in the
CUDA framework the size of a block of threads should be assigned a multiple of
32, which is the size of the workload “unit” called _warp_ , with the optimal
performances that are typically obtained setting the minimum value at 128..

Figure 3: An emulsion is prepared in an initial honeycomb structure (top
panel), then it settles in a structured configuration (bottom panel). This
situation is taken as the initial condition of a convection simulation
experiment (see Section 3.3). At the preparation step, the temperature field
is fixed to be linear between the walls, but it does not influence the fluid
flow (alphaG = 0). For details on how to produce these plots see the
README_howToPreparation file.

Each component comprises bulk regions with initial densities assigned via the
parameters rhoMax and rhoMin. The bulk densities will change during the
simulation to match the values compatible with the interactions chosen. Thus,
it is preferable to choose the values of rhoMax and rhoMin once the phase-
separation diagram is obtained from a dedicated experiment. For an easier
understanding, in Fig. 8(a) we report the phase-separation diagram, showing
the densities of each component as a function of ${\cal G}_{12}$ (for fixed
competing interactions). The corresponding density values are given in 8(b).
More details on the phase-separation experiment are given in B.
The bare viscosity $\eta_{0}$ in Eq. (14) is fixed by the relaxation time tau
in the LB equation (1). Moreover, the strength of the phase segregation
interaction in Eq. (7) can be tuned with G12 and rho0, and the competing
interactions in Eqs. (10) and (11) are handled via G11a, G22a, G11r, and G22r
666In the nomenclature of these constants, the numbers are referred to the
species (1 and 2) and the $a$ and $r$ indicate the nature of the interaction
(attractive or repulsive).. The droplets are created via the balance of phase
segregation and competing interactions and can be initialised in many
different ways. In the present version of TLBfind, we initialise the soft
domains in an honeycomb structure (droplet initialisation = 1), where we
establish the total number of droplets $N_{\mbox{\tiny droplets}}$ by
calibrating the relative number of droplets along $x$ and $y$
777$N_{\mbox{\tiny droplets}}=$ number of droplets x$\cdot$ number of droplets
y, the droplet diameter, and the (spacing) between two consecutive droplets
(see Fig. 3, top panel). Then the system settles in a structured
configuration, which will be the starting point of the convection simulation
step (see Fig. 3, bottom panel).
In this test case, we sketch out a system confined between two walls in the
$y$-direction 888It means that periodic boundary condition along y = 0 with
periodic boundary condition along x as shown in Fig. 3 999In TLBfind no wall
is implemented with the normal parallel to the $x$-direction.. By considering
flat walls, the Boolean rough-wall-flags, roughWallUp and roughWallDown, must
be fixed to zero 101010Note that, if roughWallUp and roughWallDown are true
and periodic boundary condition along y = 1, simultaneously, there will be an
inconsistency and the code will not work correctly., and the wetting
properties have to be properly calibrated by changing the densities of the two
species at wall with rhoWallMax and rhoWallMin, respectively. A flag.dat file
containing the wall (binary) density field is necessary to perform every
simulation: it is the output file of the buildObstacle program, which in turn
reads input parameters from inputflag.inp. The parameters used for the flat
walls case are listed in Table 4, and they will be discussed in detail in the
next Section 111111We wish to remark that the wall parameters shared in
inputflag.inp and tlbfind.inp files should be equal..
The above-listed information is sufficient to start from scratch (Boolean,
true) a simulation of a static fluid for a number of nsteps of simulation
times.

####  3.2.2 Rheological characterisation

Before performing a simulation of an emulsion under thermal convection, we
need to learn about the system response to an applied external driving: a
shear rheology experiment is of fundamental importance for this purpose. The
rheological characterisation is carried out in a Couette channel [17, 65],
where a shear flow is applied through walls moving along $x$ with opposite
directions, with velocities uWallUp and uWallDown. Fig. 4 shows the flow curve
for an emulsion initialised in Fig. 3: in this experiment, the shear rate
$\dot{\gamma}$ can be calculated as

$\dot{\gamma}=\dfrac{\mbox{{uWallUp}}-\mbox{{uWallDown}}}{\mbox{{ny}}},$ (23)

while the stress $\Sigma$ is computed by averaging in space and time the total
stress, given by the sum of all components contributions 121212This
information is dumped in Pxy_av.#.dat files (for more details on the output
files structure see the README_Output file and Section 2.1 as a reference)..

Figure 4: Rheological characterisation of the emulsion shown in Fig. 3: the
stress $\Sigma$ is reported as a function of the shear rate $\dot{\gamma}$.
All dimensional quantities are given in simulation units.

####  3.2.3 Temperature field

One of the key point of TLBfind lies in the possibility to simulate the
dynamics of the temperature field $T({\bf x},t)$ in Eq. (16), and couple it
with the momentum dynamics. One can decide to activate the dynamics of $T({\bf
x},t)$ by setting to 1 the Boolean parameter THERMAL. The thermal relaxation
time tauG ($\tau_{g}$ in Eq. (16)) impacts the thermal diffusivity $\kappa$
(see Eq. (20)).
In the setup of the Rayleigh-Bénard thermal convection, the system is confined
between two walls at a different temperature: a hot bottom wall at temperature
Tdown and a top cold one at temperature Tup (see Fig. 5). The temperature
profile is here initialised with a linear profile, corresponding to the
temperature initialization=2. If needed (e.g. study of situations at constant
temperature) the temperature can be initialised with a constant profile
(temperature initialization=0).

###  3.3 Simulation in convection

All values of simulation parameters in tlbfind.inp that need to be changed to
start a convective simulation are listed in Table 5. A step-by-step
description of the commands needed to run a simulation in convection can be
found in the README_howToRun file.

####  3.3.1 The onset of convection

At this point, we are ready to perform a simulation of an emulsion in
convection. The onset of convection is known to be dependent on the magnitude
of the buoyancy term in Eq. (22): at fixed walls temperatures, system size,
and fluid properties, alphaG is the order parameter that triggers the
transition from a conductive to a convective state. Note that, a finite
perturbation is necessary to destabilise the conductive state; this initial
velocity perturbation depends on the emulsion rheology [66]: the more the
system is non-Newtonian, the greater the required perturbation is.

With the aim to restart a simulation from the last preparation configuration
(i.e., situation in bottom Panel of Fig. 3), TLBfind requires four files: the
three final populations files of the components 1,2 and temperature from the
preparation step (i.e., conf1.in_0, conf2.in_0, and confG.in_0, respectively);
the dumpcount.in file, which counts the enumeration of time steps and dumping
number of output files 131313The dumpcount.in file must be manually created..
At $t=0$ it is equal to a series of seven zeros 0 0 0 0 0 0 0; it is very
useful especially in the case of a restart of a finished or suspended
simulations 141414For details on how to restart a finished or suspended
simulation see the README_howToRestart file.. As highlighted in Table 5, at
this step start from scratch has to be switched off while post preparation
temp has to be switched on 151515It is helpful to distinguish the case of
simulation after preparation from a simple restart: in the latter case no
initial velocity perturbation has to be applied, thus post preparation temp is
false..

Figure 5: Rayleigh-Bénard convection experiment: the convection simulation
step. Black arrows refer to droplet displacements: the typical convective
rolls are visible. For details on how to produce this plot see the
README_howToRun file.

Fig. 5 shows the emulsion prepared in Fig. 3 when subjected to a sufficiently
large buoyancy force: the convective state, with the characteristic circular
convective rolls, is highlighted thanks to the overlap to the density map of
Lagrangian droplet displacement vector fields $\mbox{\bf
d}(\mathbf{X}_{i}(t),t)$ (black arrows), computed for all droplet centre-of-
mass $\mathbf{X}_{i}(t)$ ($i=1...N_{\mbox{\tiny droplets}}$). Both centre-of-
mass positions and displacement of all droplets at any time can be computed as
part of the Delaunay analysis (see Section 2.4). The threshold used in order
to identify the droplets region can be tuned via the bubbleThreshold (i.e.,
$B_{th}$) parameter: it defines the floating-point value used as a threshold
on the first density field to detect the high-density compact regions
constituting the droplets. Then, nmindelaunay sets the starting time for the
Delaunay analysis. This can be useful for discarding the initial transient of
a flow from the analysis. From the knowledge of centre-of-mass positions, the
computation of droplet displacement is an immediate consequence. As mentioned
in Section 3.1, Delaunay binary files can be analysed in order to extract
droplets information by compiling the source codes of the following two
programs:

  1. 1.

deltaAnalysis, which executes the analysis for the calculation of the
displacement of the centre of mass of the droplets yielding further
information as well; the main output is contained in the ASCII file
DeltaField, see the README_Analysis file for further details;

  2. 2.

dropletStats, which appends in the ASCII files dropletsHostNow and
dropletsHostPast some relevant information about the droplets, such as centre-
of-mass coordinates and size (see the README_Analysis file for details).

####  3.3.2 Heat transfer analysis

---
(a)
(b)
Figure 6: Panel (a): Time behaviour of the Nusselt number defined in Eq. (24)
(expressed in simulation units) for emulsion of Fig. 5 in the statistically
steady state. The average Nusselt number is indicated with the dashed line.
Panel (b): a snapshot of the droplet Nusselt number defined in Eq. (25). For
further details on how producing these plots see the README_howToRun file.

The heat-flux properties of the system can be studied in TLBfind. A measure of
the average heat flux can be performed in terms of the Nusselt number (Nu), a
dimensionless observable defined as [67, 68, 69, 70]

$\mbox{Nu}(t)=\frac{\langle
u_{y}(x,y,t)T(x,y,t)\rangle_{x,y}-\kappa\langle\partial_{y}T(x,y,t)\rangle_{x,y}}{\kappa\frac{\Delta
T}{\mbox{{ny}}}},$ (24)

where $\langle(\dots)\rangle_{x,y}$ denotes a space average, and $\Delta
T=\mbox{{Tdown}}-\mbox{{Tup}}$. Nu gives an estimate of the balance between
convective and conductive transport at macroscopic scales. The program
nusseltNumber allows computing Nu as a function of time (see the
README_Analysis file for details), as Fig. 6(a) shows for the emulsion
displayed in Fig. 5. Because of the finite-size effects at mesoscales [42], in
some cases there may be the need to deal with the heat flux contributions of
the single droplets. For this reason, a Nusselt number associated with each
single droplet can be defined as [42]:

$\mbox{Nu}^{(\mbox{\tiny
drop})}_{i}(t)=\frac{u^{(i)}_{y}(t)T^{(i)}(t)-\kappa(\partial_{y}T)^{(i)}(t)}{\kappa\frac{\Delta
T}{\mbox{{ny}}}},$ (25)

where $u^{(i)}_{y}(t)=u_{y}(\mathbf{X}_{i}(t),t)$,
$T^{(i)}(t)=T(\mathbf{X}_{i}(t),t)$ and
$(\partial_{y}T)^{(i)}(t)=\partial_{y}T(\mathbf{X}_{i}(t),t)$ are the fluid
velocity, temperature, and temperature gradient evaluated for the $i$-th
droplet, respectively.

The program nusseltNumberDroplet computes Eq. (25) for all droplets at any
simulation time step by reading velocity and temperature fields and combining
them with the centres-of-mass positions extracted from Delaunay analysis (see
the README_Analysis file for further details).

##  4 Test case with rough walls

In the previous section, we have discussed a test case of Rayleigh-Bénard
numerical experiment in the presence of flat walls.Further complexity may be
introduced by the presence of roughness along the walls. Here we report the
case of rough walls with an asymmetric trapezoidal shape, as shown in Fig.
7(a). The height, the width, the periodicity (via lambda), and the angle (via
alpha) can be varied. Note that the rough-wall-flags roughWallUp and
roughWallDown have to be switched on in order to design the roughness at the
walls. Table 6 lists roughness parameters in the inputflag.dat file used in
the present test case. We build here a very general roughness shape by
accessing all input parameters involved. As a result of this code, a symmetric
trapezoidal-shape roughness was used in [48], where simulations with a
symmetric roughness (asymmetric = 0) built only on the bottom wall
(roughWallUp = 0) were performed. Moreover, in a very similar context [29,
30], the dynamics of emulsions in a channel with a simple rectangular-shape
roughness (alpha = 0.0) was analysed in detail.
In order to perform simulation of Rayleigh-Bénard convection, the steps to
follow are the same as shown for the flat-walls case: after the system
preparation step, we run the simulation with the emulsion in convection. For
both these steps, we modified the input files as well. Input parameters
changed with respect to the flat case are shown in Table 7. A snapshot of an
emulsion in convection with these boundary conditions is shown in Fig. 7(b).

|
---|---
(a) | (b)
Figure 7: Panel (a): sketch of a roughness on the bottom wall (roughWallDown =
1): the obstacles have asymmetric trapezoidal-shape with inclination given by
the angle alpha, specific height, width, and repeated every lambda lattice
nodes. Panel (b): a snapshot of Rayleigh-Bénard cell in the proximity to the
wall during a convection experiment. Black arrows indicate the droplet
displacement field.

##  5 Conclusions

We have presented TLBfind, a high-performance CUDA code for the simulation of
thermal flows of finite-size droplet suspensions with non-trivial boundary
conditions. TLBfind implements a multi-component Lattice Boltzmann (LB) method
where the non-ideal interactions are tuned in such a way to promote inhibition
of coalescence of adjacent droplets. With the LB method, it is possible to
simulate advection and diffusion of a temperature field coupled to the
droplets dynamics, and easily implement boundary conditions with complex
geometry. We demonstrated various code features by means of two examples,
corresponding to the case with flat and rough walls. We expect to work in the
future on the extension of the present code to three dimensions. A multi-GPU
implementation is also part of our future plans.

##  6 Acknowledgements

We thank Roberto Benzi, Andrea Scagliarini and Sauro Succi. FP thanks Fabio
Bonaccorso for technical support in accessing and operating GPUs at the
University of Rome Tor Vergata. ML acknowledges the support of the National
Science Foundation of China Grant 12050410244.

##  Appendix A Input parameters lists

Output files parameters
---
nout density |  10000 |  write energy file |  1
write vtk file |  0 |  nout energy |  1000
write vtk file rho1 |  1 |  nout tensor |  10000000
write vtk file rho2 |  0 |  nout average |  1000
nout velocity |  1000 |  noutconfig |  1000
nout temperature |  1000 |  noutconfigMod |  2
write vtk file temperature |  0 |  ncheckdelaunay |  1000
Table 2: Output files parameters in tlbfind.inp. All ASCII files of the
various two-dimensional fields can be plotted with pm3d on gnuplot. All VTK
files are VTK 2.0 binary files. Input preparation parameters
---
nx |  947 |  number of droplet x |  31
ny |  431 |  number of droplet y |  16
rhoMax |  1.4 |  diameter |  23.92
rhoMin |  0.1 |  spacing |  6.4
tau |  1.0 |  WD |  0.6
bubbleThreshold |  0.8 |  threhold_WD |  0.4
rho0 |  0.83 |  THERMAL |  1
G12 |  0.405 |  temperature initialisation |  2
G11a |  -9.0 |  Tup |  -0.5
G22a |  -8.0 |  Tdown |  0.5
G11r |  -8.1 |  tauG |  1.0
G22a |  -7.1 |  alphaG |  0.0
rhoWallMax |  0.612 |  periodic boundary condition along x |  1
rhoWallMin |  0.612 |  periodic boundary condition along y |  0
start from scratch |  1 |  nsteps |  $10^{5}$
droplet initialisation |  1 |  initial velocity perturbation |  0.0
delaunayDebug |  0 |  |
Table 3: Preparation parameters in tlbfind.inp. WD and threshold_WD are used
to introduce randomization in the initial condition; the used values produce
the initial condition displayed in Fig. 3. Input flat wall parameters
---
roughWallUp |  0 |  width |  0
roughWallDown |  0 |  alpha |  0.0
lambda |  0 |  asymmetric |  0
height |  0 |  |
Table 4: Input parameters in inputflag.inp for a flat wall. Input simulation
parameters
---
start from scratch |  0 |  alphaG |  $1.24\ 10^{-5}$
nsteps |  $10^{6}$ |  initial velocity perturbation |  $10^{-4}$
post preparation temp |  1 |  |
Table 5: Input parameters varied in tlbfind.inp to start a convective
simulation after the preparation step. Input rough wall parameters
---
roughWallUp |  1 |  width |  10
roughWallDown |  1 |  alpha |  0.25
lambda |  150 |  asymmetric |  1
height |  12 |  |
Table 6: Input parameters in inputflag.inp for a rough wall as in Fig. 7.
Input simulation parameters
---
roughWallUp |  1 |  rhoWallMax |  6.0
roughWallDown |  1 |  rhoWallMin |  0.1
Table 7: Input parameters varied (compared to the flat wall case) in
tlbfind.inp to run the simulation with rough walls discussed in the text (see
Tables 3 and 5).

##  Appendix B Dependence of densities $\rho_{1}$ and $\rho_{2}$ on the
coupling parameter ${\cal G}_{12}$

As discussed in Section 3.2.1, initial bulk densities are assigned via the
parameters rhoMax and rhoMin. A dedicated phase-separation experiment is
needed in order to set these values in such a way that they will not change
sensibly during the simulation. Here we provide an experiment example: we
performed some dedicated simulations with parameters in Table 8, and we
measured the maximum and minimum density values from the last configuration of
the first component (see Fig. 8(a)) by varying only the value of G12 161616The
same is valid if the measure is done with the second component, but with an
inversion of the maximum and minimum values.. In this kind of experiment, as
well as for the test cases presented in the text, the preparation-step
features must be followed, as explained in Section 3.2, with droplet
initialization = 3. With this initialisation, we prepared a simple system of
two layers: one half of the domain is occupied with the majority of the first
component, and the second half with the second component. Then, we apply no
force, but we wait a sufficient number of time steps to observe a stabilised
phase separation, and we look at the last density configuration of the first
component. The resulting phase-separation diagram is shown in Fig. 8(a), and
the corresponding measured values of bulk densities are listed in Table 8(b).
Note that, in this specific case, for values of ${\cal G}_{12}$ below
$\approx$ 0.25, the minimum density value coincides with the maximum one: in
such conditions, there is no phase separation and the two components mix.

Phase-segregation experiment parameters
---
nx |  128 |  THERMAL |  0
ny |  2 |  start from scratch |  1
rhoMax |  1.4 |  nsteps |  $5\ 10^{4}$
rhoMin |  0.1 |  periodic boundary condition along x |  1
tau |  1.0 |  periodic boundary condition along y |  1
rho0 |  0.83 |  droplet initialisation |  3
G12 |  0.405 |  nout density |  10000
G11a |  -9.0 |  write vtk file |  0
G22a |  -8.0 |  write vtk rho1 |  1
G11r |  -8.1 |  noutconfig |  10000
G22a |  -7.1 |  noutconfigMod |  2
Table 8: Phase-segregation experiment parameters. Regarding the output files
and preparation parameters present in Tables 2 and 3, and not listed here,
they can be all switched off because not necessary.
---
(a)
| ${\cal G}_{12}$ | $\rho_{\mbox{\tiny max}}$ | $\rho_{\mbox{\tiny min}}$ | ${\cal G}_{12}$ | $\rho_{\mbox{\tiny max}}$ | $\rho_{\mbox{\tiny min}}$
---|---|---|---|---|---
0.05 | 0.7601748 | 0.7601377 | 0.50 | 1.445688 | 0.08590994
0.10 | 0.7603156 | 0.7599969 | 0.55 | 1.470143 | 0.06602583
0.15 | 0.7614925 | 0.7588201 | 0.60 | 1.487723 | 0.05179497
0.20 | 0.7708472 | 0.7494691 | 0.65 | 1.500535 | 0.04135837
0.25 | 0.8346988 | 0.6855015 | 0.70 | 1.509980 | 0.03360042
0.30 | 1.108622 | 0.3624054 | 0.75 | 1.517030 | 0.02778616
0.35 | 1.267322 | 0.2425093 | 0.80 | 1.522366 | 0.02340400
0.40 | 1.364685 | 0.1549724 | 0.85 | 1.526470 | 0.02008613
0.45 | 1.410851 | 0.1149639 | 0.90 | 1.529679 | 0.01756290
(b)
Figure 8: Panel (a): phase-separation diagram for a system with initial
densities rhoMax=1.4 and rhoMin=0.1 and fixed competing interactions given in
Table 8. Panel (b) shows values of densities as a function of ${\cal G}_{12}$
producing the upper plot.

##  Appendix C Output Files

Here we schematically report the list of the output files generated by the
simulation program itself (not the analysis programs) along with a brief
description and the associated flags in the input file

Simulation output files
---
init_rho1.dat
init_rho2.dat
init_temperature.dat |  Initial conditions of the components densities and the temperature fields
firstdensity.#.dat
seconddensity.#.dat |  Components density fields at the # step. Input file option: nout density $>0$

firstdensity.dat
seconddensity.dat |  Components density fields at the last iteration

firstdensity.#.vtk |  Input file options: write vtk file and write vtk file rho1 set to 1

seconddensity.#.vtk |  Input file options: write vtk file and write vtk file rho2 set to 1

temperature.#.dat |  Temperature field at the # step: input option nout temperature $>0$

temperature.dat |  same as above with last iteration data

temperature.#.vtk |  Input file options: nout temperature $>0$ and write vtk file temperature set to 1
Table 9: On the left, list of simulation output files for densities and
temperature fields. On the right, the relevant parameters of the input file
tlbfind.inp. Simulation output files
---
veloconf.#.dat |  Velocity field at the # step. Input file option nout velocity$>0$

veloconf.dat |  Velocity field at the last iteration

veloconf.#.vtk |  Input file options: nout velocity $>0$ and write vtk file set to 1

u_av.#.dat |  x-wise x-component average of the velocity field. Input file options: nout velocity $>0$, nout average $>0$ and write vtk file set to 0

u_av.dat |  same as above with last iteration data
Ptot_xy.#.dat |  Stress tensor field at the # step. Input file option: nout tensor $>0$

Ptot_xy.dat |  same as above with last iteration data

Pxy_av.#.dat |  x-wise average of the stress tensor field at the # step. Input file options: nout tensor $>0$ and nout average $>0$
Pxy_av.dat |  same as above with last iteration data
Table 10: On the left, list of simulation output files for the velocity and
stress tensor fields. On the right, the relevant parameters of the input file
tlbfind.inp. Simulation output files
---
timeEnergy.dat |  Input file option: write energy file set to 1 and nout energy $>0$
conf1_#.in |  First component populations field dumped as a binary file. Input file options: noutconfig $>0$, noutconfigMod $>1$

conf2_#.in |  Second component populations field dumped as a binary file. Input file options: noutconfig $>0$, noutconfigMod $>1$

confG_#.in |  Temperature field populations dumped as a binary file. Input file oprions: noutconfig $>0$, noutconfigMod $>1$
Table 11: On the left, list of simulation output files for the system energy
as a function of time and the binary dumps for the populations. On the right,
the relevant parameters of the input file tlbfind.inp. Simulation Delaunay
output files
---
./delaunayTriggerDir |  Directory created if ncheckdelaunay$>0$ containing a plethora of files which are described below

arisingLinks |  File containing triangulation links created after a plastic event
arisingLinksBoundary |  File containing links created after a plastic event one the boundary of the triangulation

breakingLinks |  File containing triangulation links vanishing after a plastic event
breakingLinksBoundary |  File containing links vanishing after a plastic event on the boundary of the triangulation

nBubblesOut |  Time sequence of the number of droplets. It can be used to monitor the coalescence events

nLinksOut |  Time sequence of the number of links in the Delaunay triangulation

delaunayNowTime#
delaunayPastTime#
delaunayIsoTriggerTime# |  Binary files needed for the analysis of plastic events using the programs that can be compiled from the source files deltaAnalysis.cu and dropletStats.cu

./delaunayNoTriggerDir |  Directory containing only the binary files described above related to the frames for which no plastic event is detected
Table 12: On the left, list of simulation directories and output files for
data related to the Delaunay triangulation and the detection of plastic
events. On the right, the relevant parameters of the input file.

## References

  * [1] G. Hetsroni, Handbook of multiphase systems, McGraw-Hill Book Co., New York, NY, 1982.
  * [2] C. Gallegos, J. Franco, Rheology of food, cosmetics and pharmaceuticals, Current opinion in colloid & interface science 4 (4) (1999) 288–293.
  * [3] B. A. Khan, N. Akhtar, H. M. S. Khan, K. Waseem, T. Mahmood, A. Rasul, M. Iqbal, H. Khan, Basics of pharmaceutical emulsions: A review, African Journal of Pharmacy and Pharmacology 5 (25) (2011) 2715–2725.
  * [4] J. Shao, J. Darkwa, G. Kokogiannakis, Review of phase change emulsions (pcmes) and their applications in hvac systems, Energy and buildings 94 (2015) 200–217.
  * [5] D. J. McClements, Food emulsions: principles, practices, and techniques, CRC press, 2015.
  * [6] M. Yukuyama, D. Ghisleni, T. Pinto, N. Bou-Chacra, Nanoemulsion: process selection and application in cosmetics–a review, International journal of cosmetic science 38 (1) (2016) 13–24.
  * [7] F. Wang, W. Lin, Z. Ling, X. Fang, A comprehensive review on phase change material emulsions: Fabrication, characteristics, and heat transfer performance, Solar Energy Materials and Solar Cells 191 (2019) 218–234.
  * [8] H. Princen, A. Kiss, Rheology of foams and highly concentrated emulsions: Iv. an experimental study of the shear viscosity and yield stress of concentrated emulsions, J. Colloid Interface Sci. 128 (1) (1989) 176 – 187.  doi:https://doi.org/10.1016/0021-9797(89)90396-2.
URL http://www.sciencedirect.com/science/article/pii/0021979789903962

  * [9] L. Schramm, Emulsions, American Chemical Society (ACS), 1992.
  * [10] H. A. Barnes, Rheology of emulsions—a review, Colloids and Surfaces A: Physicochemical and Engineering Aspects 91 (1994) 89–95.
  * [11] P. Coussot, Rheometry of Pastes, Suspensions, and Granular Materials, Wiley-Interscience, 2005.
  * [12] J.-L. Barrat, Elasticity and plasticity of disordered systems, a statistical physics perspective, Physica A: Statistical Mechanics and its Applications 504 (2017) 20–30.  doi:10.1016/j.physa.2017.11.146.
  * [13] D. J. Stein, F. J. Spera, Rheology and microstructure of magmatic emulsions: theory and experiments, Journal of Volcanology and Geothermal Research 49 (1-2) (1992) 157–174.
  * [14] Y. Zhou, D. Yin, W. Chen, B. Liu, X. Zhang, A comprehensive review of emulsion and its field application for enhanced oil recovery, Energy Science & Engineering 7 (4) (2019) 1046–1058.
  * [15] P. W. Egolf, M. Kauffeld, From physical properties of ice slurries to industrial ice slurry applications, International journal of refrigeration 28 (1) (2005) 4–12.
  * [16] R. Pal, Shear viscosity behavior of emulsions of two immiscible liquids, Journal of colloid and interface science 225 (2) (2000) 359–366.
  * [17] R. G. Larson, The Structure and Rheology of Complex Fluids, Oxford University Press, 1999.
  * [18] D. Bonn, M. M. Denn, L. Berthier, T. Divoux, S. Manneville, Yield stress materials in soft condensed matter, Rev. Mod. Phys. 89 (2017) 035005.  doi:10.1103/RevModPhys.89.035005.
URL https://link.aps.org/doi/10.1103/RevModPhys.89.035005

  * [19] P. K. Kilpatrick, Water-in-crude oil emulsion stabilization: review and unanswered questions, Energy & Fuels 26 (7) (2012) 4017–4026.
  * [20] S. N. Kale, S. L. Deore, Emulsion micro emulsion and nano emulsion: a review, Systematic Reviews in Pharmacy 8 (1) (2017) 39.
  * [21] D. J. McClements, S. M. Jafari, Improving emulsion formation, stability and performance using mixed emulsifiers: A review, Advances in colloid and interface science 251 (2018) 55–79.
  * [22] A. Rahimian, S. K. Veerapaneni, G. Biros, Dynamic simulation of locally inextensible vesicles suspended in an arbitrary two-dimensional domain, a boundary integral method, Journal of Computational Physics 229 (18) (2010) 6466–6484.
  * [23] W. Yang, Z. Zhou, A. Yu, D. Pinson, Particle scale simulation of softening–melting behaviour of multiple layers of particles in a blast furnace cohesive zone, Powder Technology 279 (2015) 134–145.
  * [24] M. Kroupa, M. Vonka, M. Soos, J. Kosek, Utilizing the discrete element method for the modeling of viscosity in concentrated suspensions, Langmuir 32 (33) (2016) 8451–8460.
  * [25] J. R. Seth, L. Mohan, C. Locatelli-Champagne, M. Cloitre, R. T. Bonnecaze, A micromechanical model to predict the flow of soft particle glasses, Nature materials 10 (11) (2011) 838–843.
  * [26] V. Mansard, A. Colin, P. Chaudhuri, L. Bocquet, A molecular dynamics study of non-local effects in the flow of soft jammed particles, Soft Matter 9 (2013) 7489–7500.  doi:10.1039/C3SM50847A.
URL http://dx.doi.org/10.1039/C3SM50847A

  * [27] G. Jung, S. M. Fielding, Wall slip and bulk yielding in soft particle suspensions, Journal of Rheology 65 (2) (2021) 199–212.
  * [28] B. Dollet, A. Scagliarini, M. Sbragaglia, Two-dimensional plastic flow of foams and emulsions in a channel: experiments and lattice boltzmann simulations, J. Fluid Mech. 766 (2015) 556–589.  doi:10.1017/jfm.2015.28.
URL https://www.cambridge.org/core/article/two-dimensional-plastic-flow-of-
foams-and-emulsions-in-a-channel-experiments-and-lattice-boltzmann-
simulations/A8E06679EA25FD323CEDD1695246EF8B

  * [29] L. Derzsi, D. Filippi, G. Mistura, M. Pierno, M. Lulli, M. Sbragaglia, M. Bernaschi, P. Garstecki, Fluidization and wall slip of soft glassy materials by controlled surface roughness, Phys. Rev. E 95 (5) (2017) 052602.
  * [30] L. Derzsi, D. Filippi, M. Lulli, G. Mistura, M. Bernaschi, P. Garstecki, M. Sbragaglia, M. Pierno, Wall fluidization in two acts: from stiff to soft roughness, Soft matter 14 (7) (2018) 1088–1093.
  * [31] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E. M. Viggen, The lattice boltzmann method, Springer International Publishing 10 (978-3) (2017) 4–15.
  * [32] S. Succi, The lattice Boltzmann Equation, Oxford University Press, 2018.
  * [33] J.-C. Desplat, I. Pagonabarraga, P. Bladon, Ludwig: A parallel lattice-boltzmann code for complex fluids, Computer Physics Communications 134 (3) (2001) 273–290.
  * [34] S. Schmieschek, L. Shamardin, S. Frijters, T. Krüger, U. D. Schiller, J. Harting, P. V. Coveney, Lb3d: A parallel implementation of the lattice-boltzmann method for simulation of interacting amphiphilic fluids, Computer Physics Communications 217 (2017) 149–161.
  * [35] F. Bonaccorso, A. Montessori, A. Tiribocchi, G. Amati, M. Bernaschi, M. Lauricella, S. Succi, Lbsoft: A parallel open-source software for simulation of colloidal systems, Computer Physics Communications 256 (2020) 107455\.
  * [36] J. Latt, O. Malaspinas, D. Kontaxakis, A. Parmigiani, D. Lagrava, F. Brogi, M. B. Belgacem, Y. Thorimbert, S. Leclaire, S. Li, et al., Palabos: parallel lattice boltzmann solver, Computers & Mathematics with Applications 81 (2021) 334–350.
  * [37] M. Ataei, V. Shaayegan, F. Costa, S. Han, C. B. Park, M. Bussmann, Lbfoam: An open-source software package for the simulation of foaming using the lattice boltzmann method, Computer Physics Communications 259 (2021) 107698.
  * [38] D. Moore, N. Weiss, Two-dimensional rayleigh-bénard convection, Journal of Fluid Mechanics 58 (2) (1973) 289–312.
  * [39] H. Bénard, Les tourbillons cellulaires dans une nappe liquide, Rev. Gen. Sci. Pures Appl. 11 (1900) 1261–1271.
  * [40] L. Rayleigh, Lix. on convection currents in a horizontal layer of fluid, when the higher temperature is on the under side, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 32 (192) (1916) 529–546.
  * [41] D. Lohse, K.-Q. Xia, Small-scale properties of turbulent rayleigh-bénard convection, Annual Review of Fluid Mechanics 42 (2010) 335–364.
  * [42] F. Pelusi, M. Sbragaglia, R. Benzi, A. Scagliarini, M. Bernaschi, S. Succi, Rayleigh–bénard convection of a model emulsion: anomalous heat-flux fluctuations and finite-size droplet effects, Soft Matter 17 (2021) 3709–3721.  doi:10.1039/D0SM01777A.
URL http://dx.doi.org/10.1039/D0SM01777A

  * [43] R. Benzi, M. Sbragaglia, S. Succi, M. Bernaschi, S. Chibbaro, Mesoscopic lattice boltzmann modeling of soft-glassy systems: Theory and simulations, J. Chem. Phys. 131 (10) (2009).  doi:http://dx.doi.org/10.1063/1.3216105.
URL http://scitation.aip.org/content/aip/journal/jcp/131/10/10.1063/1.3216105

  * [44] M. Sbragaglia, R. Benzi, M. Bernaschi, S. Succi, The emergence of supramolecular forces from lattice kinetic models of non-ideal fluids: applications to the rheology of soft glassy materials, Soft Matter 8 (2012) 10773–10782.  doi:10.1039/C2SM26167G.
URL http://dx.doi.org/10.1039/C2SM26167G

  * [45] R. Benzi, M. Sbragaglia, P. Perlekar, M. Bernaschi, S. Succi, F. Toschi, Direct evidence of plastic events and dynamic heterogeneities in soft-glasses, Soft Matter 10 (2014) 4615–4624.  doi:10.1039/C4SM00348A.
URL http://dx.doi.org/10.1039/C4SM00348A

  * [46] L. Fei, A. Scagliarini, A. Montessori, M. Lauricella, S. Succi, K. H. Luo, Mesoscopic model for soft flowing systems with tunable viscosity ratio, Physical Review Fluids 3 (10) (2018) 104304.
  * [47] F. Pelusi, M. Sbragaglia, R. Benzi, Avalanche statistics during coarsening dynamics, Soft Matter 15 (22) (2019) 4518–4524.
  * [48] F. Pelusi, M. Sbragaglia, A. Scagliarini, M. Lulli, M. Bernaschi, S. Succi, On the impact of controlled wall roughness shape on the flow of a soft material, EPL (Europhysics Letters) 127 (3) (2019) 34005.  doi:10.1209/0295-5075/127/34005.
  * [49] P. L. Bhatnagar, E. P. Gross, M. Krook, A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems, Physical review 94 (3) (1954) 511.
  * [50] X. Shan, H. Chen, Lattice boltzmann model for simulating flows with multiple phases and components, Physical review E 47 (3) (1993) 1815.
  * [51] X. Shan, H. Chen, Simulation of nonideal gases and liquid-gas phase transitions by the lattice boltzmann equation, Physical Review E 49 (4) (1994) 2941.
  * [52] G. Falcucci, G. Bella, G. Chiatti, S. Chibbaro, M. Sbragaglia, S. Succi, et al., Lattice boltzmann models with mid-range interactions, Communications in computational physics 2 (6) (2007) 1071–1084.
  * [53] G. Falcucci, S. Ubertini, S. Succi, Lattice boltzmann simulations of phase-separating flows at large density ratios: the case of doubly-attractive pseudo-potentials, Soft Matter 6 (18) (2010) 4357–4365.
  * [54] M. Sbragaglia, R. Benzi, L. Biferale, S. Succi, K. Sugiyama, F. Toschi, Generalized lattice boltzmann method with multirange pseudopotential, Physical Review E 75 (2) (2007) 026702.
  * [55] R. Benzi, M. Sbragaglia, A. Scagliarini, P. Perlekar, M. Bernaschi, S. Succi, F. Toschi, Internal dynamics and activated processes in soft-glassy materials, Soft Matter 11 (2015) 1271–1280.  doi:10.1039/C4SM02341B.
URL http://dx.doi.org/10.1039/C4SM02341B

  * [56] P. Ripesi, L. Biferale, M. Sbragaglia, A. Wirth, Natural convection with mixed insulating and conducting boundary conditions: low-and high-rayleigh-number regimes, Journal of fluid mechanics 742 (2014) 636–663.
  * [57] B. Delaunay, Sur la sphère vide, Bulletin de l’Académie des Sciences de l’URSS, Classe des sciences mathématiques et naturelles 6 (1934).
URL http://mi.mathnet.ru/eng/izv4937

  * [58] M. Bernaschi, M. Lulli, M. Sbragaglia, GPU based detection of topological changes in voronoi diagrams, Comput. Phys. Commun. 213 (2017) 19 – 28.  doi:http://dx.doi.org/10.1016/j.cpc.2016.11.005.
URL http://www.sciencedirect.com/science/article/pii/S0010465516303599

  * [59] R. H. Swendsen, J.-S. Wang, Nonuniversal critical dynamics in Monte Carlo simulations, Physical Review Letters 58 (2) (1987) 86–88.  doi:10.1103/PhysRevLett.58.86.
URL
http://link.aps.org/doi/10.1103/PhysRevLett.58.86%5Cnhttp://kobus.ca/teaching/cs645/spring09/ua_only/swendsen-
wang-87.pdfhttps://link.aps.org/doi/10.1103/PhysRevLett.58.86

  * [60] Y. Komura, Y. Okabe, GPU-based Swendsen–Wang multi-cluster algorithm for the simulation of two-dimensional classical spin systems, Computer Physics Communications 183 (6) (2012) 1155–1161.  arXiv:1202.0635, doi:10.1016/j.cpc.2012.01.017.
URL
http://dx.doi.org/10.1016/j.cpc.2012.01.017https://linkinghub.elsevier.com/retrieve/pii/S001046551200032X

  * [61] G. Voronoi, Nouvelles applications des paramètres continus à la théorie des formes quadratiques. Deuxième mémoire. Recherches sur les parallélloèdres primitifs., Journal für die reine und angewandte Mathematik (Crelles Journal) 1908 (134) (1908) 198–287.  doi:10.1515/crll.1908.134.198.
URL https://www.degruyter.com/document/doi/10.1515/crll.1908.134.198/html

  * [62] M. Januszewski, M. Kostur, Sailfish: A flexible multi-GPU implementation of the lattice boltzmann method, Computer Physics Communications 185 (9) (2014) 2350–2368.  doi:10.1016/j.cpc.2014.04.018.
URL https://doi.org/10.1016%2Fj.cpc.2014.04.018

  * [63] M. Bernaschi, L. Rossi, R. Benzi, M. Sbragaglia, S. Succi, Graphics processing unit implementation of lattice boltzmann models for flowing soft systems, Phys. Rev. E 80 (2009) 066707.  doi:10.1103/PhysRevE.80.066707.
URL http://link.aps.org/doi/10.1103/PhysRevE.80.066707

  * [64] J. Ahrens, S. Jourdain, P. O'Leary, J. Patchett, D. H. Rogers, M. Petersen, An image-based approach to extreme scale in situ visualization and analysis, in: SC14: International Conference for High Performance Computing, Networking, Storage and Analysis, IEEE, 2014.  doi:10.1109/sc.2014.40.
URL https://doi.org/10.1109%2Fsc.2014.40

  * [65] R. Pal, Effect of droplet size on the rheology of emulsions, AIChE Journal 42 (11) (1996) 3181–3190.
  * [66] J. Zhang, D. Vola, I. A. Frigaard, Yield stress effects on rayleigh–benard convection, J. Fluid Mech. 566 (2006) 389–419.  doi:10.1017/S002211200600200X.
  * [67] B. I. Shraiman, E. D. Siggia, Heat transport in high-rayleigh-number convection, Physical Review A 42 (6) (1990) 3650.
  * [68] G. Ahlers, S. Grossmann, D. Lohse, Heat transfer and large scale dynamics in turbulent rayleigh-bénard convection, Reviews of modern physics 81 (2) (2009) 503.
  * [69] R. J. Stevens, R. Verzicco, D. Lohse, Radial boundary layer structure and nusselt number in rayleigh-bénard convection, Journal of Fluid Mechanics 643 (2010) 495–507.
  * [70] F. Chillà, J. Schumacher, New perspectives in turbulent rayleigh-bénard convection, The European Physical Journal E 35 (7) (2012) 1–25.