3D_LBM

[Armin268]

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()

[pmocz]

I tried out the code. If I set V0 = 0, it gives me the correct result that z velocity is 0. If I set it to 2 like you have, then I find the z velocity is 0.07. Maybe you were expecting it to be 2 in that case? It is not 2 because you are adding terms to an un-normalized distribution function F, so you either need to add the appropriate normalization where you add V0 to F, or multiply V0 to get the desired velocity in the simulation

[Armin268]

Thank you for the prompt answer. Oh, really? I tried to check the velocity in z direction by "print (uz)" after "# Apply boundary" section. But even after setting V0 = 0, I got something around 0.07267267! Could you please explain how you achieved a zero velocity in the z-direction by V0 = 0? Thank you.

[Armin268]

I have one more inquiry related to the obstacles (cylinder) in your code. I would like to incorporate moving spherical obstacles in this 3D code, where the particles move in each iteration. My ultimate goal is to couple LBM with DEM ! I must say, your code has been incredibly helpful for me. I am just curious whether it is feasible to implement the particle scenario I described within your code. Thank you for your valuable assistance.