Lattice Boltzmann fluid simulation
Philip Mocz (2020) Princeton Univeristy, @PMocz
Simulate a flow past cylinder with the Lattice Boltzmann method
python latticeboltzmann.py
Lattice Boltzmann fluid simulation
Home Page: https://medium.com/swlh/create-your-own-lattice-boltzmann-simulation-with-python-8759e8b53b1c
License: GNU General Public License v3.0
Lattice Boltzmann fluid simulation
Simulate a flow past cylinder with the Lattice Boltzmann method
python latticeboltzmann.py
Hi @pmocz ,
I wanted to change the Reynolds number of the simulation to create an ensemble with e.g. Re from 50 to 250. I know the Reynolds number equation but am not sure which variables to change in the code, e.g. dynamic viscosity, fluid velocity, or fluid density. I tried to change the average density (rho0
) but it didn't affect the flow. What is the best way to change the Reynolds number in this simulation? For now, I only managed to modify the flow by changing the cylinder boundary parameters.
Thank you for providing the code. I attempted to transform the 2D code into a 3D (D3Q15) simulation. Although the new code runs, I encountered an issue with the computed velocities not matching the desired results accurately. For instance, when attempting to simulate still water (no velocity) in the Z direction, the computed value was 0.07, which significantly deviates from the expected value of 0.
I would appreciate your assistance in resolving this matter. Please let me know if you need any further information or if there are specific areas of the code that require attention to achieve the desired simulation outcome.
import numpy as np
import os
import sys
import math
import re
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def main ():
#######################################################################
# Model Parameters
Nx = 15 # Number of lattices in x-direction
Ny = 15 # Number of lattices in Y-direction
Nz = 30 # Number of lattices in Z-direction
rho0 = 1000 # Fluid density
mu = 1e-3 # Fluid dynamic viscosity
tau = 0.6 #single-relaxation-time (Equation 3)
V0 = 2 # initial velocity
r0 = 0.5 # Particle Radius
n = 4 # n th Lattice node
cs_square = 1 / 3 # Speed of sound square (lattice-dependent)_ Check Table_3
dx = 1.0 # Lattice spacing in the x-direction
dy = 1.0 # Lattice spacing in the y-direction
dz = 1.0 # Lattice spacing in the z-direction
Nt = 4000
#######################################################################
# Lattice speeds & weights _ Check Table_3
NL = 15
idxs = np.arange(NL)
lx = np.array([0, 1, -1, 0, 0, 0, 0, 1, -1, 1, -1, 1, -1, -1, 1])
ly = np.array([0, 0, 0, 0, 0, 1, -1, 1, -1, -1, 1, 1, -1, 1, -1])
lz = np.array([0, 0, 0, 1, -1, 0, 0, 1, -1, 1, -1, -1, 1, 1, -1])
weights = np.array([2/9, 1/9, 1/9, 1/9, 1/9, 1/9, 1/9, 1/72, 1/72, 1/72, 1/72, 1/72, 1/72, 1/72, 1/72])
#######################################################################
# Initial Conditions
F = np.ones((Nx, Ny, Nz, NL)) # * rho0 / NL # Make the calculation more stable by "np.ones"
np.random.seed(42)
F += 0.01 * np.random.randn(Nx, Ny, Nz, NL)
X, Y, Z = np.meshgrid(range(Nx), range(Ny), range(Nz))
F[:, :, :, n] = V0 * (1 + 0.2 * np.cos(2 * np.pi * X / Nx * 4)) * np.sin(2 * np.pi * Y / Ny * 2) * np.cos(2 * np.pi * Z / Nz * 4) # velocity in the n th lattice node direction
rho = np.sum(F,3)
for i in idxs:
F[:, :, :,i] *= rho0 / rho
for it in range(Nt):
#print (it)
# Read ball coordinates from "ball_coor.txt"
with open("ball_coor.txt", "r") as f:
ball_x = float(f.readline().strip())
ball_y = float(f.readline().strip())
ball_z = float(f.readline().strip())
# Convert the coordinates to three decimal places
ball_x = round(ball_x, 3)
ball_y = round(ball_y, 3)
ball_z = round(ball_z, 3)
X, Y, Z = np.meshgrid(range(Nx), range(Ny), range(Nz))
particle = (X - ball_x)**2 + (Y - ball_y)**2+ (Z - ball_z)**2 < (r0)**2
#Streaming
for i, cx, cy, cz in zip(range(NL), lx, ly, lz):
F[:, :, :, i] = np.roll(F[:, :, :, i], cx, axis=0)
F[:, :, :, i] = np.roll(F[:, :, :, i], cy, axis=1)
F[:, :, :, i] = np.roll(F[:, :, :, i], cz, axis=2)
# Set reflective boundaries
bndryF = F[particle,:]
bndryF = bndryF[:, [0, 2, 1, 4, 3, 6, 5, 8, 7, 10, 9, 12, 11, 14, 13]]
# Fluid variables
rho = np.sum(F, 3)
print (rho)
#Velocity components
ux = np.sum(F * lx, 3) / rho
uy = np.sum(F * ly, 3) / rho
uz = np.sum(F * lz, 3) / rho
#Pressure
p = cs_square * rho
#Gradient of pressure components
px = (np.roll(p, -1, axis=0) - np.roll(p, 1, axis=0)) / (2 * dx)
py = (np.roll(p, -1, axis=1) - np.roll(p, 1, axis=1)) / (2 * dy)
pz = (np.roll(p, -1, axis=2) - np.roll(p, 1, axis=2)) / (2 * dz)
# Collision _ Equation 15
Feq = np.zeros(F.shape)
for i, cx, cy, cz, w in zip(range(NL), lx, ly, lz, weights):
Feq[:, :, :, i] = rho * w * (1 + 3 * (cx * ux + cy * uy + cz * uz) + 9 * (cx * ux + cy * uy + cz * uz)**2/2 - 3 * (ux**2 + uy**2 +
uz**2)/2)
F = F + -(1/tau)*(F - Feq)
# Apply boundary
F[particle,:] = bndryF
ux[particle] = 0
uy[particle] = 0
uz[particle] = 0
if name== "main":
main()
Thank you
How to get all points with colors and velocity values?
Get the following warning
latticeboltzmann.py:90: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("bwr").copy()
cmap.set_bad('black')
Using Python 3.10.5 with Matplotlib 3.5.2
In the collision step it looks like the python code uses ux * ux + uy * uy = dot(u, u)
dot(u,u)^2
:On a different note, thanks for the great work on the small, self-contained simulation examples, they are very educational.
Hi,
I noticed that while the simulation continues the values of the flow field are getting smaller and smaller. That is the case if I use e.g. 10K time steps of the simulation. Why is that so? Is it meant to be this way? I wanted to use this data in my project related to flow learning via neural networks and values which are becoming very small is somewhat problematic.
Is it possible to add a way to simulate the diffusion of a hot flow of air in a colder one?
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.