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 [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

import ipywidgets as widgets
from ipywidgets import interactive

Flow Definition

In [14]:
# 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 [15]:
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 [16]:
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 [17]:
def equilibrium(rho, u):
    """Equilibrium distribution function"""
    #Velocity square
    usqr = 3/2 * (u[0] ** 2 + u[1] **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 [18]:
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 [19]:
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 [20]:
fin = equilibrium(1, vel)
In [21]:
%matplotlib inline
In [22]:
from IPython.display import display
In [23]:
from time import sleep
In [24]:
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[0]**2 + u[1]**2).transpose(), cmap=cm.Reds)
            plt.show()
            
In [25]:
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)
widgets.jslink((play, 'value'), (slider, 'value'))
# interactive_plot = interactive(inter1, play)
play.observe(inter1, 'value')
# output = interactive_plot.children[-1]
# output.layout.height = '350px'
widgets.VBox([play, slider, out])

And here the oscilation pattern will be called by Karman Vortex Street.