pmocz / latticeboltzmann-python Goto Github PK

# 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

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