**2. Computational domains**

The flow within a cavity of height **h** and wide **w** where the bottom wall is moving at constant velocity *U*<sup>0</sup> Fig.1 is going to be model. The cavity is completely filled by an incompresible fluid with constant density *ρ* and cinematic viscosity *ν*.

Fig. 1. Cavity

### **3. Flow modelling by LBM with vorticity stream-function variables**

Is important to introduce the equations that govern the vorticity transport and a few definitions that will be used during the present study.

**Definition 0.1.** *A vortex is a set of fluid particles that moves around a common center*

The vorticity vector is defined as *ω* = ∇ × *v* and its transport equation is given by

$$\frac{\partial \omega}{\partial t} + [\nabla \omega] v = [\nabla v] \omega + \nu \nabla^2 \omega. \tag{1}$$

which is obtained by calculating the *curl* of the NS equation. For a 2D flow Eq.(1) is simplified to obtain

$$\frac{\partial \omega}{\partial t} + [\nabla \omega] v = \nu \nabla^2 \omega. \tag{2}$$

In order to recover the velocity field from the vorticity field the Poisson equation for the stream function needs to be solved. The Poisson equation wich involves the stream function is stated as

$$
\nabla^2 \psi = -\omega \tag{3}
$$

where *ψ* is the stream function who carries the velocity field information as

$$
\mu = \frac{\partial \psi}{\partial y} \; \; \; v = -\frac{\partial \psi}{\partial x} \; \; \; \; \tag{4}
$$

and ensures the mass conservation. The motivation for adopting vorticity as the primitive variables lies in the fact that every potential, as the pressure, is eliminated which is physicaly desirable because being the vorticity an angular velocity, the pressure, which is always normal to the fluid can not affect the angular momentum of a fluid element.

### **3.1 Numerical method**

2 Will-be-set-by-IN-TECH

The flow within a cavity of height **h** and wide **w** where the bottom wall is moving at constant velocity *U*<sup>0</sup> Fig.1 is going to be model. The cavity is completely filled by an incompresible

**2. Computational domains**

Fig. 1. Cavity

to obtain

as

fluid with constant density *ρ* and cinematic viscosity *ν*.

definitions that will be used during the present study.

*∂ω*

**3. Flow modelling by LBM with vorticity stream-function variables**

*∂ω*

where *ψ* is the stream function who carries the velocity field information as

to the fluid can not affect the angular momentum of a fluid element.

*<sup>u</sup>* <sup>=</sup> *∂ψ*

**Definition 0.1.** *A vortex is a set of fluid particles that moves around a common center*

The vorticity vector is defined as *ω* = ∇ × *v* and its transport equation is given by

Is important to introduce the equations that govern the vorticity transport and a few

which is obtained by calculating the *curl* of the NS equation. For a 2D flow Eq.(1) is simplified

In order to recover the velocity field from the vorticity field the Poisson equation for the stream function needs to be solved. The Poisson equation wich involves the stream function is stated

*<sup>∂</sup><sup>y</sup>* , *<sup>v</sup>* <sup>=</sup> <sup>−</sup>*∂ψ*

and ensures the mass conservation. The motivation for adopting vorticity as the primitive variables lies in the fact that every potential, as the pressure, is eliminated which is physicaly desirable because being the vorticity an angular velocity, the pressure, which is always normal

*<sup>∂</sup><sup>t</sup>* + [∇*ω*]*<sup>v</sup>* = [∇*v*]*<sup>ω</sup>* <sup>+</sup> *<sup>ν</sup>*∇2*ω*. (1)

*<sup>∂</sup><sup>t</sup>* + [∇*ω*]*<sup>v</sup>* <sup>=</sup> *<sup>ν</sup>*∇2*ω*. (2)

<sup>∇</sup>2*<sup>ψ</sup>* <sup>=</sup> <sup>−</sup>*<sup>ω</sup>* (3)

*<sup>∂</sup><sup>x</sup>* . (4)

Consider a set of particles that moves in a bidimensional lattice and each particle with a finite number of movements. Now a vorticity distribution function *gi*(*x*, *t*) will be asigned to each particle with unitary velocity *ei* giving to it a dynamic consistent with two principles:


Fig. 2. D2Q5 Model.2 dimensions and 5 possible directions of moving

**Observation 0.2.** *The method only considers binary particle collisions.*

The evolution equation is discribed by

$$g\_k(\vec{\mathfrak{x}} + c\vec{e}\_k \Delta t, t + \Delta t) - g\_k(\vec{\mathfrak{x}}, t) = -\frac{1}{\tau} [g\_k(\vec{\mathfrak{x}}, t) - g\_k^{\epsilon q}(\vec{\mathfrak{x}}, t)]^1 \tag{5}$$

where *ek* are the posible directions where the vorticity can be transported as shown in Fig.2. *c* = Δ*x*/Δ*t* is the fluid particle speed, Δ*x* and Δ*t* the lattice grid spacing and the time step respectively and *τ* the dimensionless relaxation time. Clearly Eq.(5) is divided in two parts, the first one emulates the advective term of (1) and the collision term, which is in square brackets, emulates the diffusive term of equation (1).

The equilibrium function is calculed by

$$g\_k^{eq} = \frac{w}{5} \left[ 1 + 2.5 \frac{\vec{e\_k} \cdot \vec{u}}{c} \right]. \tag{6}$$

The vorticity is calculed as

$$w = \sum\_{k\geq 0} \mathbf{g}\_k \tag{7}$$

and *τ*, the dimensionless relaxation time, is determined by Re number

$$Re = \frac{5}{2c^2(\tau - 0.5)}.\tag{8}$$

<sup>1</sup> The evolution equations were taken from (Chen et al., 2008) and (Chen, 2009). Is strongly recomended to consult the latter references for a deeper understanding of the evolution equations and parameter calculations.

of Lid-Driven Cavities 5

Flow Evolution Mechanisms of Lid-Driven Cavities 415

8. Solution of Poisson equation: In order to solve Poisson equation the evolution equation Eq.(9) for the stream-function distribution was implemented within a loop wishing to

While the simulations were ran, it was found that the algorithm was demanding finer meshes for higher Re numbers, i.e. 700x700 nodes mesh for Re 6,000, increasing the computational cost and most of the times ending in overflows own to the fluid regime. To overcome this

The principal characteristic of a turbulent flow is that its velocity field is of random nature. Considering this, the velocity field can be split in a deterministic term and in a random term i.e. *U*(*x*, *t*) = *U*¯ (*x*, *t*) + *u*(*x*, *t*), being the deterministic and random term respectively. In order to solve the velocity field, the NS equations are recalculated in deterministic variables adding to the set a closure equation own to the loss of information undertaken by solving only the deterministic term. At introducing a turbulent model there exist three different approaches: algebraic models, closure models and Large Eddy Simulations (LES) being the latter used in the present study. LES were introduce by James Deardorff on 1960 (Durbin & Petersson-Rief, 2010). Such simulations are based in the fact that the bulk of the system energy is contained in the large eddys of the flow making not neccesary to calculate all the vortex disipative range which would imply a high computational cost (Durbin & Petersson-Rief, 2010). If small scales are ommited, for example by increasing the spacing by a factor of 5, the number of grid points is substantially reduced by a factor of 125 (Durbin & Petersson-Rief, 2010). In LES context the elimination of these small scales is called filtering. But this filtering or omission of small scales is determined as follows: the dissipative phenomenon is replaced by an alternative that produces correct dissipation levels without requiring small scale simulations. The Smagorinsky model was introduced where another flow viscosity (usually known as subgrid viscosity) is considered which is calculated based on the fluid deformation stress.

*<sup>k</sup>* <sup>−</sup> *fk*<sup>|</sup> <sup>&</sup>lt; <sup>10</sup><sup>−</sup>3.

*Dt* <sup>=</sup> <sup>∇</sup>2*<sup>ψ</sup>* <sup>+</sup> *<sup>ω</sup>* <sup>=</sup> 0. For the latter

3. Velocity field calculation using Eq.(4)

5. Colission term calculation using Eq.(12) 6. Probability transport using Eq.(13) 7. Vorticity field calculation using Eq.(7)

loop the process terminated when

**4. Introduction of turbulence in LBM**

4. Equilibrium probability calculation using Eq.(6)

compare *fk*'s values (i.e. *<sup>ψ</sup>*) aiming to achive that *<sup>D</sup><sup>ψ</sup>*

Specifically it is model as *<sup>ν</sup><sup>t</sup>* = (*C*Δ)2|*S*|Chen et al. (2008) where

*∂ω*

*Sij* <sup>=</sup> <sup>1</sup> 2

Assuming this new **subgrid viscosity** *ν<sup>t</sup>* the momentum equation is given by

*∂x νe ∂ω ∂x* + *∂ ∂y νe ∂ω ∂y* 

*<sup>∂</sup><sup>t</sup>* + [∇*ω*]*<sup>v</sup>* <sup>=</sup> *<sup>∂</sup>*

 *∂U*¯*<sup>i</sup> ∂xi* + *∂U*¯ *<sup>j</sup> ∂xj*

Δ is the filter width and *C* the Smagorinsky constant. In the present study *C* = 0.1 and Δ = Δ*x*.

 ,

∑*x*,*y* | *f* +

situations a turbulence model was introduced to the LBM proposed by (Chen, 2009).

In order to calculate the velocity field Poisson equation must be solved (3). In order to do this (Chen et al., 2008) introduces another evolution equation.

$$f\_k(\vec{x} + c\vec{e}\_k \Delta t, t + \Delta t) - f\_k(\vec{x}, t) = \Omega\_k + \hat{\Omega}\_k. \tag{9}$$

Where

$$
\Omega\_{\vec{k}} = \frac{-1}{\tau\_{\Psi}} [f\_{\vec{k}}(\vec{x}, t) - f\_{\vec{k}}^{\varepsilon} q(\vec{x}, t)]\_{\prime} \,\, \widehat{\Omega}\_{\vec{k}} = \Delta t \xi\_{\vec{k}}^{\varepsilon} \theta D \tag{10}
$$

and *D* = *<sup>c</sup>*<sup>2</sup> <sup>2</sup> (0.5 − *τψ*). *τψ* is the dimensionless relaxation time of the latter evolution equation wich can be chosen arbitrarly. For the sake of understanding the evolution equations, the equation (9) consist on calculating *<sup>D</sup><sup>ψ</sup> Dt* <sup>=</sup> <sup>∇</sup>2*<sup>ψ</sup>* <sup>+</sup> *<sup>ω</sup>* until *<sup>D</sup><sup>ψ</sup> Dt* = 0, having found a solution *ψ* for the Poisson equation.

By last, the equlibrium distribution function is defined as

$$f\_k^{eq} = \begin{cases} \zeta\_k \psi \ k = 1, 2, 3, 4 \\ -\psi \quad k = 0 \end{cases} \tag{11}$$

where *ξ<sup>k</sup>* and *ζ<sup>k</sup>* are weight parameters of the equation.

#### **3.2 Algorithm implementation**

In order to implement the evolution equation Eq.(5) two main calculations are considered. First, the collision term is calculated as

$$\mathcal{g}\_{k}^{\rm int} = -\frac{1}{\tau} [\mathcal{g}\_{k}(\vec{x}, t) - \mathcal{g}\_{k}^{\rm eq}(\vec{x}, t)] \tag{12}$$

and next the vorticity distributions is transported as

$$\mathbf{g}\_k(\vec{\mathbf{x}} + c\vec{e}\_k \Delta t, t + \Delta t) = \mathbf{g}\_k^{\text{int}} + \mathbf{g}\_k(\vec{\mathbf{x}}, t) \tag{13}$$

which is, as mentioned, the basic concept that governs the LBM, collisions and transportation of determined distribution in our case a vorticity distibution.

#### **3.2.1 Algorithm and boundary conditions**

1. Paramater Inicialization


$$
\omega|\_{\partial\Omega} = \frac{7\psi\_w - 8\psi\_{w-1} + \psi\_{w-2}}{2\Delta n^2} \tag{14}
$$

$$
\omega|\_{\partial\Omega} = \frac{7\psi\_w - 8\psi\_{w-1} + \psi\_{w-2}}{2\Delta n^2} - \frac{3\mathcal{U}\_0}{\Delta n} \tag{15}
$$

Both equations came from solving Poisson equation Eq.(3) on the walls by a second order Taylor approximation. Eq.(15) is used on the moving wall nodes.

<sup>2</sup> For the sake of clarity Re number is imposed in the method by the user which intrinsically is imposing different flow viscosities.

3. Velocity field calculation using Eq.(4)

4 Will-be-set-by-IN-TECH

In order to calculate the velocity field Poisson equation must be solved (3). In order to do this

wich can be chosen arbitrarly. For the sake of understanding the evolution equations, the

In order to implement the evolution equation Eq.(5) two main calculations are considered.

*<sup>τ</sup>* [*gk*(*<sup>x</sup>*, *<sup>t</sup>*) <sup>−</sup> *<sup>g</sup>*

which is, as mentioned, the basic concept that governs the LBM, collisions and transportation

*<sup>ω</sup>*|*∂*<sup>Ω</sup> <sup>=</sup> <sup>7</sup>*ψ<sup>w</sup>* <sup>−</sup> <sup>8</sup>*ψw*−<sup>1</sup> <sup>+</sup> *<sup>ψ</sup>w*−<sup>2</sup>

*<sup>ω</sup>*|*∂*<sup>Ω</sup> <sup>=</sup> <sup>7</sup>*ψ<sup>w</sup>* <sup>−</sup> <sup>8</sup>*ψw*−<sup>1</sup> <sup>+</sup> *<sup>ψ</sup>w*−<sup>2</sup>

Both equations came from solving Poisson equation Eq.(3) on the walls by a second order

<sup>2</sup> For the sake of clarity Re number is imposed in the method by the user which intrinsically is imposing

<sup>2</sup>Δ*n*<sup>2</sup> <sup>−</sup> <sup>3</sup>*U*<sup>0</sup>

*gk*(*x* + *cek*Δ*t*, *t* + Δ*t*) = *gint*

*eq*

*Dt* <sup>=</sup> <sup>∇</sup>2*<sup>ψ</sup>* <sup>+</sup> *<sup>ω</sup>* until *<sup>D</sup><sup>ψ</sup>*

*<sup>ζ</sup>k<sup>ψ</sup> <sup>k</sup>* = 1, 2, 3, 4

<sup>2</sup> (0.5 − *τψ*). *τψ* is the dimensionless relaxation time of the latter evolution equation

[ *fk*(*<sup>x</sup>*, *<sup>t</sup>*) <sup>−</sup> *<sup>f</sup> <sup>e</sup>*

*fk*(*<sup>x</sup>* <sup>+</sup> *<sup>c</sup>ek*Δ*t*, *<sup>t</sup>* <sup>+</sup> <sup>Δ</sup>*t*) <sup>−</sup> *fk*(*<sup>x</sup>*, *<sup>t</sup>*) = <sup>Ω</sup>*<sup>k</sup>* <sup>+</sup> <sup>Ω</sup><sup>ˆ</sup> *<sup>k</sup>*. (9)

*<sup>k</sup> q*(*x*, *t*)], Ω *<sup>k</sup>* = Δ*tξkθD* (10)

<sup>−</sup>*<sup>ψ</sup> <sup>k</sup>* <sup>=</sup> <sup>0</sup> (11)

*Dt* = 0, having found a solution *ψ*

*<sup>k</sup>* (*x*, *t*)] (12)

*<sup>k</sup>* + *gk*(*x*, *t*) (13)

<sup>2</sup>Δ*n*<sup>2</sup> (14)

<sup>Δ</sup>*<sup>n</sup>* (15)

(Chen et al., 2008) introduces another evolution equation.

<sup>Ω</sup>*<sup>k</sup>* <sup>=</sup> <sup>−</sup><sup>1</sup> *τψ*

By last, the equlibrium distribution function is defined as

where *ξ<sup>k</sup>* and *ζ<sup>k</sup>* are weight parameters of the equation.

and next the vorticity distributions is transported as

*f eq <sup>k</sup>* =

*gint <sup>k</sup>* <sup>=</sup> <sup>−</sup> <sup>1</sup>

of determined distribution in our case a vorticity distibution.

• *ψ*|*∂*<sup>Ω</sup> = 0, own to the fact that no particle is crossing the walls. • *u* = *v* = 0 in the whole cavity excepting the moving wall.

Taylor approximation. Eq.(15) is used on the moving wall nodes.

equation (9) consist on calculating *<sup>D</sup><sup>ψ</sup>*

for the Poisson equation.

**3.2 Algorithm implementation**

First, the collision term is calculated as

**3.2.1 Algorithm and boundary conditions**

• Moving wall velocity: *U*<sup>0</sup> = 1.

1. Paramater Inicialization

• Re number definition<sup>2</sup> 2. Wall vorticity calculation

different flow viscosities.

Where

and *D* = *<sup>c</sup>*<sup>2</sup>


$$\sum\_{\mathbf{x}, \mathbf{y}} |f\_{\mathbf{k}}^{+} - f\_{\mathbf{k}}| < 10^{-3}.$$

While the simulations were ran, it was found that the algorithm was demanding finer meshes for higher Re numbers, i.e. 700x700 nodes mesh for Re 6,000, increasing the computational cost and most of the times ending in overflows own to the fluid regime. To overcome this situations a turbulence model was introduced to the LBM proposed by (Chen, 2009).

#### **4. Introduction of turbulence in LBM**

The principal characteristic of a turbulent flow is that its velocity field is of random nature. Considering this, the velocity field can be split in a deterministic term and in a random term i.e. *U*(*x*, *t*) = *U*¯ (*x*, *t*) + *u*(*x*, *t*), being the deterministic and random term respectively. In order to solve the velocity field, the NS equations are recalculated in deterministic variables adding to the set a closure equation own to the loss of information undertaken by solving only the deterministic term. At introducing a turbulent model there exist three different approaches: algebraic models, closure models and Large Eddy Simulations (LES) being the latter used in the present study. LES were introduce by James Deardorff on 1960 (Durbin & Petersson-Rief, 2010). Such simulations are based in the fact that the bulk of the system energy is contained in the large eddys of the flow making not neccesary to calculate all the vortex disipative range which would imply a high computational cost (Durbin & Petersson-Rief, 2010). If small scales are ommited, for example by increasing the spacing by a factor of 5, the number of grid points is substantially reduced by a factor of 125 (Durbin & Petersson-Rief, 2010). In LES context the elimination of these small scales is called filtering. But this filtering or omission of small scales is determined as follows: the dissipative phenomenon is replaced by an alternative that produces correct dissipation levels without requiring small scale simulations. The Smagorinsky model was introduced where another flow viscosity (usually known as subgrid viscosity) is considered which is calculated based on the fluid deformation stress. Specifically it is model as *<sup>ν</sup><sup>t</sup>* = (*C*Δ)2|*S*|Chen et al. (2008) where

$$S\_{i\bar{j}} = \frac{1}{2} \left( \frac{\partial \mathcal{U}\_i}{\partial x\_{\bar{i}}} + \frac{\partial \bar{\mathcal{U}}\_{\bar{j}}}{\partial x\_{\bar{j}}} \right) \mathcal{I}$$

Δ is the filter width and *C* the Smagorinsky constant. In the present study *C* = 0.1 and Δ = Δ*x*. Assuming this new **subgrid viscosity** *ν<sup>t</sup>* the momentum equation is given by

$$\frac{\partial \omega}{\partial t} + [\nabla \omega] v = \frac{\partial}{\partial x} \left( \nu\_\varepsilon \frac{\partial \omega}{\partial x} \right) + \frac{\partial}{\partial y} \left( \nu\_\varepsilon \frac{\partial \omega}{\partial y} \right).$$

of Lid-Driven Cavities 7

Flow Evolution Mechanisms of Lid-Driven Cavities 417

Several studies have proposed to study the deep cavity geometry (Gustafson, 1991; Patil et al., 2006) but none has reached to simulate high Re numbers possibly because the mesh sizes. Due to the LBM low computational cost it was decided to present the study of a deep cavity with

A general description is presented emphasizing the most important configurations through

• **Step 1 Fig.4(a)** The positive vortex creates a negative vortex that arises from the right wall

• **Step 2 Fig.4 (b)** The negative vortex that arises from the right wall has taken the whole

• **Step 3 Fig.4(c,d)** Positive vortices have joined in one by an interesting process discribed in

• **Step 4 Fig.4(e)** The positive vortex expands into the cavity moving upward the negative vortex until the steady state is reached in which both vortices occupy the same space of the cavity. Is worth to notice that this vortex distribution is not achieved in the square cavity

During the evolution it was observed that after positive vortices joined (Fig.4(c, d)) the new big positive vortex acted as a moving wall for the negative vortex injecting vorticity to it. Reproducing the behavior seen in the square cavity, now by the negatie vortex. Ergo a *quasi square cavity* was created in the top of the cavity but instead having a moving wall it had a vortex. The phenomenon is shown in Fig.5 where it is clear that the top of the deep cavity is a "reflection" of the square cavity with respect to an imaginary vertical axis drawn between

A particular process for Re 10,000 in square and deep cavities was found to take place through evolution. This process occurs several times throughout evolution, named Vortex Binding. In this process isolated vortices get connected forming a "massive" vortex which eventually will configure the steady state vortices distribution. A binding process that occured through evolution is shown in Fig.6 binding a positive vortex that appeared in the upper right corner

In order to explain the binding process, which is illustrated in Fig.6, recall the vorticity transport equation Eq.(1). The transport equation is divided in two terms that dictate the transport of vorticity, the diffusive term *<sup>ν</sup>*∇2*<sup>ω</sup>* and the advective term [∇*ω*]*v*. For a high Re number flow the diffusive term can be neglected, turning the attention in the advective term. As the flow evolved it was seen that the vorticity and stream-function contour lines tended to align as shown in Fig.7(a) making the vorticity gradient vector and velocity vector orthogonal

As shown in Fig.7(b) vorticity contour lines started to curve, due to its own vorticity, crossing with the stream-function contour lines and making [∇*ω*]*v* �= 0. In Fig.7(b)can be seen that

with the positive vortex that came from the movement of the bottom wall.

at different places (Fig.7(a)) causing [∇*ω*]*v* = 0, i.e. no vorticity transport.

triggering an interaction since the begining of the evolution.

Sec6. This union creates a "mirror" phenomenon inside the cavity.

cavity confining the positive vortex to the bottom.

**5.1 Deep cavities**

**5.1.1 Vortex dynamics**

steady state.

these two.

**6. Vortex binding**

**5.1.2 Mirror phenomenon**

evolution to steady state:

an aspect ratio (AR) of 1.5 for Re 8,000.

where

$$\nu\_{\ell} = \nu\_{\ell} + \nu\_{\ast}$$

As the transport equation has changed, the LBM evolution equation has also changed

$$\mathcal{g}\_k(\vec{\mathfrak{x}} + c\vec{e}\_k \Delta t, t + \Delta t) - \mathcal{g}\_k(\vec{\mathfrak{x}}, t) = -\frac{1}{\tau\_\mathcal{\mathfrak{x}}} [\mathcal{g}\_k(\vec{\mathfrak{x}}, t) - \mathcal{g}\_k^{eq}(\vec{\mathfrak{x}}, t)] \tag{16}$$

where

$$\tau\_{\mathfrak{e}} = \tau + \frac{5(\mathbf{C}\Lambda)^2|\mathbf{S}|}{2c^2\Delta t} \text{ and } |\mathbf{S}| = |\omega|^3.$$

Having a new evolution equation Eq.(16) the algorithm has to be modified adding a new step where *τ<sup>e</sup>* is calculated based on the vorticity field. After making this improvement to the method, the algorithm began to work eficiently allowing to achive higher Re numbers without compromising the computer cost, justifing the use of a LBM.

### **5. Steady state study for different Re numbers**

It is said that the flow has reached steady state when collisions and transport do not affect each node probability. Concerning the algorithm it was considered that the flow had reached the steady state when its energy had stabilized and when the maps of vorticity and stream function showed no changes through time.

Steady state vortex configuration for Re 1,000 and Re 10,000 is shown in Fig.3. It worth to notice that both are very similar, a positive vortex that fills the cavity and two negative vortices at the corners of the cavity. This configuration was observed from Re 1,000 to Re 10,000 being a prime characteristic of cavity flows. It is also important to clarify that for Re 10,000 the steady state presents a periodicity which is located in the upper left vortex that we shall see later, indeed Fig.3(b) is a "snapshot" of the flow.

(a) Stream-function map in steady state for Re 1,000. (b) Stream-function map in steady state for Re 10,000.

Fig. 3. Steady states. Maps were taken at 100,000 and 110,000 iterations respectively.

<sup>3</sup> Is strongly recomended to consult (Chen, 2009) for a deeper understanding of the evolution equations and parameter calculations.
