Lattice Boltzmann
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.