# Lattice Boltzmann

|   Source

This notebook is based on Python script for Coursera's Simulation and Natural Processes

This notebook simulates fluid disperse when facing obstacle. Using IPython widgets to animate image simulation.

In :
```import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

import ipywidgets as widgets
from ipywidgets import interactive
```

## Flow Definition¶

In :
```# Total number of Iteration
maxIter = 200000
#Reynolds number
# Re = 10.0
Re = 28
#Number of lattice nodes
nx, ny = 420, 180
#Height of the domain lattice units
ly = ny - 1
# Coordinates of the cylinder obstacles
cx, cy, r = nx//4, ny//2, ny//9
# Velocity in lattice units
uLB = 0.04
# Viscocity in lattice units
nulb = uLB*r/Re;
# Relaxation parameters
omega = 1 / (3 * nulb + 0.5)
```

## Lattice Constants¶

In :
```v = np.array([[1,1], [1,0], [1, -1], [0,1], [0,0],
[0,-1], [-1,1], [-1,0], [-1,-1]])
# Distribution normal (Maxwell-Botlzmann) with biggest probability at
# 1. rest/neutral point 4/9
# 2. orthogonal, 1/9
# 3. diagonal, 1/36
# t compensate different length velocities v. Diagonal slower velocities.
t = np.array([1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36])

#left
col1 = np.array([0,1,2])
#mid
col2 = np.array([3,4,5])
#right
col3 = np.array([6,7,8])
```

## Function Definition¶

In :
```def macroscopic(fin):
rho = np.sum(fin, axis=0)
u = np.zeros((2, nx, ny))
for i in range(9):
u[0, :, :] += v[i, 0] * fin[i, :, :]
u[1, :, :] += v[i, 1] * fin[i, :, :]
u /= rho
return rho, u
```
In :
```def equilibrium(rho, u):
"""Equilibrium distribution function"""
#Velocity square
usqr = 3/2 * (u ** 2 + u **2)
feq = np.zeros((9, nx, ny))
for i in range(9):
cu = 3 * (v[i,0] * u[0, :, :] + v[i, 1] * u[1, :, :])
feq[i, :, :] = rho * t[i] * ( 1 + cu + 0.5*cu**2 - usqr)
return feq
```

## Setup Cyndrical Obstacle¶

Also add velocity inlet with perturbation. Creation of a mask with 1/0 values, defining the shape of the obstacle.

In :
```def obstacle_fun(x, y):
return (x-cx)**2 + (y-cy)**2 < r**2
obstacle = np.fromfunction(obstacle_fun, (nx, ny))
```

## Initiate Velocity Profile¶

Almost Zero with a slight peturbation to trigger the instability

In :
```def inivel(d, x, y):
return (1-d) * uLB * (1 + 1e-4 * np.sin(y / ly*2*np.pi))
vel = np.fromfunction(inivel, (2, nx, ny))
```

Initialize populations at equilibrium with the given velocity

In :
```fin = equilibrium(1, vel)
```
In :
```%matplotlib inline
```
In :
```from IPython.display import display
```
In :
```from time import sleep
```
In :
```out = widgets.Output()

obstacle = np.fromfunction(obstacle_fun, (nx, ny))
vel = np.fromfunction(inivel, (2, nx, ny))
fin = equilibrium(1, vel)

def inter1(value):
sleep(0.01)
fin[col3, -1, :] = fin[col3, -2, :]

#Compute macroscopiv variables, desnity, and velocity
rho, u = macroscopic(fin)

#Left wall, inflow condition
u[:,0,:]= vel[:,0,:]
rho[0,:] = 1/ (1-u[0,0,:]) * ( np.sum(fin[col2, 0, :], axis=0) +
2*np.sum(fin[col3, 0, :], axis=0 ))

#Compute Equilibrium
feq = equilibrium(rho, u)
#Unknown population
fin[[0, 1, 2], 0, :] = feq[[0,1,2], 0, :] + fin[[8,7,6], 0, :] - feq[[8,7,6], 0, :]

#Collision Step
fout = fin - omega * (fin - feq)

#Bounce back condition for obstacle
for i in range(9):
fout[i, obstacle] = fin[8-i, obstacle]

#Streaming Step. Assign new directions in population with out population
for i in range(9):
fin[i, :, :] = np.roll(
np.roll(fout[i,:,:], v[i,0], axis=0),
v[i,1], axis=1
)
# Visualization of the velocity

plt.clf()
if value['new'] % 100 == 0:
with out:
out.clear_output()
plt.imshow(np.sqrt(u**2 + u**2).transpose(), cmap=cm.Reds)
plt.show()

```
In :
```play = widgets.Play(
interval=1,
value=0,
min=0,
max=maxIter,
step=1,
description="Press play",
disabled=False
)
slider = widgets.IntSlider(min=0, max=maxIter)