From 95bd3de45a60da51b91010718015a8fa0b2bcfea Mon Sep 17 00:00:00 2001 From: zomseffen Date: Tue, 3 Aug 2021 10:33:07 +0200 Subject: [PATCH] lattice Boltzmann example/playground --- FluidSim/LatticeBoltzmann.py | 181 +++++++++++++++++++++++++++++++++++ 1 file changed, 181 insertions(+) create mode 100644 FluidSim/LatticeBoltzmann.py diff --git a/FluidSim/LatticeBoltzmann.py b/FluidSim/LatticeBoltzmann.py new file mode 100644 index 0000000..3940a7a --- /dev/null +++ b/FluidSim/LatticeBoltzmann.py @@ -0,0 +1,181 @@ +import matplotlib.pyplot as plt +import numpy as np + +""" +Create Your Own Lattice Boltzmann Simulation (With Python) +Philip Mocz (2020) Princeton Univeristy, @PMocz +Simulate flow past cylinder +for an isothermal fluid +""" + + +def main(): + """ Finite Volume simulation """ + + # Simulation parameters + Nx = 400 # resolution x-dir + Ny = 100 # resolution y-dir + rho0 = 100 # average density + tau = 0.6 # collision timescale + Nt = 80000 # number of timesteps + plotRealTime = True # switch on for plotting as the simulation goes along + + # Lattice speeds / weights + NL = 9 + idxs = np.arange(NL) + cxs = np.array([0, 0, 1, 1, 1, 0, -1, -1, -1]) + cys = np.array([0, 1, 1, 0, -1, -1, -1, 0, 1]) + weights = np.array([4 / 9, 1 / 9, 1 / 36, 1 / 9, 1 / 36, 1 / 9, 1 / 36, 1 / 9, 1 / 36]) # sums to 1 + + # Initial Conditions + F = np.ones((Ny, Nx, NL)) # * rho0 / NL + has_fluid = np.ones((Ny, Nx), dtype=np.bool) + has_fluid[int(Ny/2):, :] = False + np.random.seed(42) + F += 0.01 * np.random.randn(Ny, Nx, NL) + X, Y = np.meshgrid(range(Nx), range(Ny)) + F[:, :, 3] += 2 * (1 + 0.2 * np.cos(2 * np.pi * X / Nx * 4)) + # F[:, :, 5] += 1 + rho = np.sum(F, 2) + for i in idxs: + F[:, :, i] *= rho0 / rho + + # Cylinder boundary + X, Y = np.meshgrid(range(Nx), range(Ny)) + cylinder = (X - Nx / 4) ** 2 + (Y - Ny / 2) ** 2 < (Ny / 4) ** 2 + inner_cylinder = (X - Nx / 4) ** 2 + (Y - Ny / 2) ** 2 < (Ny / 4 - 2) ** 2 + F[cylinder] = 0 + F[0, :] = 0 + F[Ny - 1, :] = 0 + # F[int(Ny/2):, :] = 0 + + has_fluid[cylinder] = False + has_fluid[0, :] = False + has_fluid[Ny - 1, :] = False + + # for i in idxs: + # F[:, :, i] *= has_fluid + + # Prep figure + fig = plt.figure(figsize=(4, 2), dpi=80) + + reflection_mapping = [0, 5, 6, 7, 8, 1, 2, 3, 4] + # Simulation Main Loop + for it in range(Nt): + print(it) + + # Drift + new_has_fluid = np.zeros((Ny, Nx)) + F_sum = np.sum(F, 2) + for i, cx, cy in zip(idxs, cxs, cys): + F_part = F[:, :, i] / F_sum + F_part[F_sum == 0] = 0 + to_move = F_part * (has_fluid * 1.0) + to_move = (np.roll(to_move, cx, axis=1)) + to_move = (np.roll(to_move, cy, axis=0)) + + new_has_fluid += to_move + + F[:, :, i] = np.roll(F[:, :, i], cx, axis=1) + F[:, :, i] = np.roll(F[:, :, i], cy, axis=0) + + # has_fluid = new_has_fluid > 0.5 + # new_has_fluid[F_sum == 0] += has_fluid[F_sum == 0] * 1.0 + # new_has_fluid[(np.abs(F_sum) < 0.000000001)] = 0 + + fluid_sum = np.sum(has_fluid * 1.0) + has_fluid = (new_has_fluid / np.sum(new_has_fluid * 1.0)) * fluid_sum + + print('fluid_cells: %d' % np.sum(has_fluid * 1)) + + # for i in idxs: + # F[:, :, i] *= has_fluid + + bndry = np.zeros((Ny, Nx), dtype=np.bool) + bndry[0, :] = True + bndry[Ny - 1, :] = True + # bndry[:, 0] = True + # bndry[:, Nx - 1] = True + bndry = np.logical_or(bndry, cylinder) + + # bndry = np.logical_or(bndry, has_fluid < 0.5) + + # Set reflective boundaries + bndryF = F[bndry, :] + bndryF = bndryF[:, reflection_mapping] + + sum_f = np.sum(F) + print('Sum of Forces: %f' % sum_f) + + # sum_f_cyl = np.sum(F[cylinder]) + # print('Sum of Forces in cylinder: %f' % sum_f_cyl) + + # sum_f_inner_cyl = np.sum(F[inner_cylinder]) + # print('Sum of Forces in inner cylinder: %f' % sum_f_inner_cyl) + + # if sum_f > 4000000.000000: + # test = 1 + + # F[Ny - 1, :, 5] += 0.1 + # F[0, :, 1] -= 0.1 + # F[0, :, 5] += 0.1 + # F[Ny - 1, :, 1] -= 0.1 + + # Calculate fluid variables + rho = np.sum(F, 2) + ux = np.sum(F * cxs, 2) / rho + uy = np.sum(F * cys, 2) / rho + + ux[(np.abs(rho) < 0.000000001)] = 0 + uy[(np.abs(rho) < 0.000000001)] = 0 + + # print('minimum rho: %f' % np.min(np.abs(rho))) + # print('Maximum F: %f' % np.max(F)) + # print('Minimum F: %f' % np.min(F)) + + # Apply Collision + Feq = np.zeros(F.shape) + for i, cx, cy, w in zip(idxs, cxs, cys, weights): + Feq[:, :, i] = rho * w * ( + 1 + 3 * (cx * ux + cy * uy) + 9 * (cx * ux + cy * uy) ** 2 / 2 - 3 * (ux ** 2 + uy ** 2) / 2) + + F += -(1.0 / tau) * (F - Feq) + + # Apply boundary + F[bndry, :] = bndryF + + # plot in real time - color 1/2 particles blue, other half red + if (plotRealTime and (it % 10) == 0) or (it == Nt - 1): + plt.cla() + ux[cylinder] = 0 + uy[cylinder] = 0 + vorticity = (np.roll(ux, -1, axis=0) - np.roll(ux, 1, axis=0)) - ( + np.roll(uy, -1, axis=1) - np.roll(uy, 1, axis=1)) + vorticity[cylinder] = np.nan + + # vorticity *= has_fluid + + cmap = plt.cm.bwr + cmap.set_bad('black') + + # plt.imshow(vorticity, cmap='bwr') + plt.imshow(has_fluid * 2.0 - 1.0, cmap='bwr') + # plt.imshow(bndry * 2.0 - 1.0, cmap='bwr') + + plt.clim(-.1, .1) + ax = plt.gca() + ax.invert_yaxis() + ax.get_xaxis().set_visible(False) + ax.get_yaxis().set_visible(False) + ax.set_aspect('equal') + plt.pause(0.001) + + # Save figure + # plt.savefig('latticeboltzmann.png', dpi=240) + plt.show() + + return 0 + + +if __name__ == "__main__": + main() \ No newline at end of file