## Theory to LBE Exercise 1: Pressure-driven flows and Kármán vortex streets¶

The Reynolds number ($$Re$$) is an important dimensionless quantity in fluid mechanics that can be used to predict flow patterns in different flow situations. Two flow systems with the same Reynolds number are regarded as behaving similarly, even if their system dimensions and fluid properties are different. This allows experimentation with smaller scale systems in preparation of larger systems with confidence that the latter will behave in an expected manner.

The definition of the Reynolds number is:

$Re = \frac{u L}{\nu} = \frac{\rho u L}{\mu}$

where $$u$$ is the (average) fluid speed with respect to a given object (system or obstacle size), $$\rho$$ is the fluid density, $$L$$ is a characteristic length scale (e.g. pipe diameter), $$\mu$$ and $$\nu$$ are respectively the dynamic and kinematic viscosities of the fluid.

One major use of the Reynolds number - which can also be regarded as the ratio of inertial to viscous forces - is to characterise flow types. For instance, for fully-developed flows in cylindrical pipes, $$Re$$ values of less than 2300 typically indicate laminar (sheet-like, viscosity-dominated) flows, while values greater than 2900 indicate turbulent (chaotic, inertially-driven) flows.

In LBE simulations, the kinematic viscosity depends upon the relaxation time:

$\nu = \frac{1}{3} \left(\tau - \frac{1}{2}\right) \frac{\Delta x^2}{\Delta t}$

We have not actively set physical length or time scales for our simulations, but we can use the system size and velocity in internal units to find the Reynolds number as these will be consistent.

For a system of pressure-driven flow between two parallel plates of distance $$2H$$ apart, we can use the Navier-Stokes equations or a balance of viscous and inertial forces to determine its velocity profile when the flow is laminar. If the flow is in the $$x$$-direction and varies with vertical position $$y$$ from the centre between the plates, the $$x$$-component of velocity will be:

$u_x (y) = \frac{1}{2 \mu} \left(-\frac{\partial p}{\partial x} \right) \left(H^2 - y^2\right)$

where $$-\frac{\partial p}{\partial x} = \Delta p / L$$ is the pressure drop along the length of plates $$L$$. The maximum velocity exists at $$y=0$$ and equals:

$u_{max} = \frac{H^2 \Delta p}{2 \mu L},$

while the mean velocity (found by integrating the profile over $$y$$ and $$z$$ and then dividing by the cross-sectional area) is equal to half of the maximum value.

In a two-dimensional lattice Boltzmann simulation, the pressure drop can be applied either by using constant density boundary conditions at each end of the simulation domain or by applying a constant body force to the fluid:

$\Delta p = c_s^2 \left(\rho_{in} - \rho_{out}\right) = \frac{F_x}{2H}$

Return to LBE Exercise 1: Pressure-driven flows and Kármán vortex streets to complete the exercise.

## Theory to LBE Exercise 2: Phase separation and equations of state¶

The Shan/Chen pseudopotential method of applying mesoscale fluid interactions is based on calculating forces that depend on gradients of a pseudopotential $$\psi$$, generally expressed as a function of density and temperature. For a one-fluid system, the interaction force acting on the fluid can be expressed by:

$\mathbf{F} \left(\mathbf{x}\right) = -\psi g \sum_i w_i \psi \left(\mathbf{x} + \mathbf{e}_i \Delta t \right) \mathbf{e}_i$

with the gradient of the pseudopotential approximated here (rather well, as it happens) by a stencil using pseudopotentials at all neighbouring grid points. The standard approach to applying this force to the LBE fluid is to modify the velocity used to calculate the local equilibrium distribution function $$f_i^{eq}$$ used in the collision operator [Martys1996]:

$\mathbf{v}_i = \mathbf{u}_i + \frac{\tau \mathbf{F}}{\rho}$

although it is also possible to use alternative schemes, e.g. the Equal Difference Method (EDM) [Kupershtokh2006] or the Guo scheme [Guo2002], to apply this force.

The above interaction force results in the following equation of state:

$p = \rho c_s^2 + \frac{1}{2} g c_s^2 \psi^2$

where the pseudopotential and the interaction strength $$g$$ control the pressure contributions that deviate from standard lattice fluid behaviour. If you wish to apply a particular equation of state, you can rearrange the above equation to find the form of the required pseudopotential:

$\psi = \sqrt{\frac{2\left(p \left(\rho, T\right) - \rho c_s^2\right)}{g c_s^2}}$

where the interaction strength $$g$$ is set to $$\pm 1$$ to ensure the square root in the above expression is evaluated with a positive value.

The Swift free-energy approach does not require interaction forces to be calculated, but uses a modified form of the local equilibrium distribution function to apply a bulk pressure and interfacial tension terms. For a two-dimensional simulation, this can be expressed as follows:

$f_i^{eq} = w_i^{00} \rho + w_i \left[ \rho \left\{ \left(\mathbf{e}_i \cdot \mathbf{u} \right) + \frac{3}{2} \left(\mathbf{e}_i \cdot \mathbf{u}\right)^2 - \frac{1}{2}u^2 \right\} + \lambda \left\{3 \left(\mathbf{e}_i \cdot \mathbf{u} \right) \left(\mathbf{e}_i \cdot \nabla \rho\right) + \left[ \gamma_i \left(\mathbf{e}_i \cdot \mathbf{e}_i \right) + \delta_i \right] \left(\mathbf{u} \cdot \nabla \rho \right) \right\} \right] + w_i^{p} P_0 - w_i^{t} \kappa \rho \nabla^2 \rho + w_i^{xx} \kappa \left(\partial_x \rho \right)^2 + w_i^{yy} \kappa \left(\partial_y \rho \right)^2 + w_i^{xy} \kappa \partial_x \rho \partial_y \rho$

where $$P_0$$ is the bulk pressure (equal to the required pressure for the equation of state $$p$$), $$\kappa$$ is the interfacial tension parameter and $$\lambda$$ is a parameter dependent on the equation of state and density that is used to ensure Galilean invariance. The required first and second order gradients of density can be calculated using stencils in a similar manner to those used for Shan/Chen interaction forces. Apart from being able to specify the interfacial tension (which is not directly possible with Shan/Chen pseudopotential interactions), the free-energy approach also guarantees correct thermodynamic behaviour.

For this exercise, the Peng-Robinson cubic equation of state [Peng1976] has been selected. This takes the form:

$p = \frac{\rho R T}{1 - b \rho} - \frac{a \alpha \left(T_r, \omega\right) \rho^2}{1 + 2 b \rho - b^2 \rho^2}$

where $$R$$ is the universal gas constant, $$a$$ is a temperature-independent attraction parameter, $$b$$ is a finite volume parameter and $$\alpha$$ is a function of reduced temperature $$T_r = \frac{T}{T_c}$$ and acentric factor $$\omega$$, usually given as:

$\alpha \left( T_r, \omega \right) = \left[1 + \left( 0.37464 + 1.54226\omega - 0.26992 \omega^2 \right)\left(1 - \sqrt{T_r}\right) \right]^2$

For a given fluid, the parameters $$a$$ and $$b$$ can be calculated from the critical temperature and pressure ($$T_c$$ and $$p_c$$):

$a \approx 0.45724 \frac{R^2 T_c^2}{p_c}$
$b \approx 0.07780 \frac{R T_c}{p_c}$

Return to LBE Exercise 2: Phase separation and equations of state to complete the exercise.

## Theory to LBE Exercise 3: Drop flows¶

The Lishchuk algorithm for modelling immiscible fluids is based upon applying interfacial stresses between them at a continuum level. This can be achieved by defining a phase index between two fluids ($$a$$ and $$b$$):

$\rho_{ab}^{N} = \frac{\rho^a - \rho^b}{\rho^a + \rho^b}$

which will equal +1 for pure fluid $$a$$ and -1 for pure fluid $$b$$. An interfacial normal can be determined using a first-order spatial differential of the phase index:

$\widehat{n}_{ab} = \frac{\nabla \rho_{ab}^{N}}{| \rho_{ab}^{N} |}$

and the local curvature at the interface can be calculated from a surface differential of the interfacial normal:

$K_{ab} = -\nabla_{S} \cdot \widehat{n}_{ab}$

These can be used to calculate the interfacial tension force acting between the two fluids:

$\mathbf{F}^{ab} = \frac{1}{2} g_{ab} K_{ab} \nabla \rho_{ab}^{N}.$

To apply this interfacial force, the individual distribution functions for the fluids are combined together into what is often described as an ‘achromatic’ distribution function (a reference to earlier forms of this algorithm that identified different fluids as colours):

$f_i = f^{a}_i + f^{b}_i$

This ‘achromatic’ distribution function is collided with a forcing term (normally the Guo forcing term [Guo2002]) and is then segregated to produce post-collisional distribution functions for each fluid using the following function:

$f_i^{a} \left(\mathbf{x}, t^{+}\right) = \frac{\rho^{a} \left(\mathbf{x}\right)}{\rho \left(\mathbf{x}\right)} f_i \left(\mathbf{x}, t^{+}\right) + \beta w_i \frac{\rho^a \left(\mathbf{x}\right) \rho^b \left(\mathbf{x}\right)}{\rho^2 \left(\mathbf{x}\right)} \mathbf{e}_i \cdot \widehat{n}_{ab} \left(\mathbf{x}\right)$

where $$\beta$$ is a segregation parameter which results in a diffuse interface between the two fluids but reduces non-physical effects such as pinning of drops to a lattice, spatial anisotropy in interfacial tension and spurious microcurrents.

It can be shown that the interfacial force parameter $$g_{ab}$$ can be directly related to the interfacial tension between the two fluids:

$\sigma_{ab} = \frac{4 g_{ab} \nu^2 \langle \rho \rangle}{c_s^4 \left(2 \tau - 1\right)^2 \delta x}$

We can define the dimensionless capillary number $$Ca$$, which describes the relative effects of viscosity and surface tension on a fluid drop surrounded by a continuous fluid. Based on the ratio of viscous to surface tension forces, it is given by the following formula:

$Ca = \frac{\mu u}{\sigma} = \frac{\rho \nu u}{\sigma}$

where $$\rho$$ is the density of the fluid drop and $$\sigma$$ is the interfacial tension between the two fluids. For drops with small capillary numbers (less than $$10^{-5}$$), we can assume flows in porous materials or capillaries will be dominated by interfacial tension.