From 8a5de47da3f02147c5128097e8bc5dcb13602637 Mon Sep 17 00:00:00 2001 From: zomseffen Date: Fri, 15 Oct 2021 19:21:17 +0200 Subject: [PATCH] lava lamp lattice boltzmann --- FluidSim/FluidSimParameters.py | 24 ++++ FluidSim/LatticeBoltzmann.py | 209 ++++++++++++++++++++++++++------- 2 files changed, 188 insertions(+), 45 deletions(-) create mode 100644 FluidSim/FluidSimParameters.py diff --git a/FluidSim/FluidSimParameters.py b/FluidSim/FluidSimParameters.py new file mode 100644 index 0000000..2170be4 --- /dev/null +++ b/FluidSim/FluidSimParameters.py @@ -0,0 +1,24 @@ +class FluidSimParameter: + viscosity = 0.1 / 3.0 + # Pr = 1.0 + Pr = 100.0 + # vc = 1.0 + vc = 0.5 + + def __init__(self, height: int): + self.t1 = 3 * self.viscosity + 0.5 + self.t2 = (2 * self.t1 - 1) / (2 * self.Pr) + 0.5 + self.g = (self.vc ** 2) / height + + self.R = self.Pr * self.g * (height ** 3) / (self.viscosity ** 2) + + +class MagmaParameter(FluidSimParameter): + viscosity = 10 ** 19 + Pr = 10 ** 25 + + +class WaterParameter(FluidSimParameter): + viscosity = 8.9 * 10 ** -4 + Pr = 7.56 + vc = 0.05 diff --git a/FluidSim/LatticeBoltzmann.py b/FluidSim/LatticeBoltzmann.py index 3940a7a..b39a8b3 100644 --- a/FluidSim/LatticeBoltzmann.py +++ b/FluidSim/LatticeBoltzmann.py @@ -1,5 +1,6 @@ import matplotlib.pyplot as plt import numpy as np +from FluidSim.FluidSimParameters import * """ Create Your Own Lattice Boltzmann Simulation (With Python) @@ -13,48 +14,63 @@ def main(): """ Finite Volume simulation """ # Simulation parameters + epsilon = 0.000000001 Nx = 400 # resolution x-dir Ny = 100 # resolution y-dir - rho0 = 100 # average density + rho0 = 1 # average density tau = 0.6 # collision timescale Nt = 80000 # number of timesteps plotRealTime = True # switch on for plotting as the simulation goes along + params = FluidSimParameter(Ny) + # params = WaterParameter(Ny) + # params = MagmaParameter(Ny) + # 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 + xx, yy = np.meshgrid(range(Nx), range(Ny)) # Initial Conditions - F = np.ones((Ny, Nx, NL)) # * rho0 / NL + N = np.ones((Ny, Nx, NL)) # * rho0 / NL + temperature = np.ones((Ny, Nx, NL), np.float) # * 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) + N += 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) + N[:, :, 3] += 2 * (1 + 0.2 * np.cos(2 * np.pi * X / Nx * 4)) + # N[:, :, 5] += 1 + rho = np.sum(N, 2) + temperature_rho = np.sum(temperature, 2) for i in idxs: - F[:, :, i] *= rho0 / rho + N[:, :, i] *= rho0 / rho + temperature[:, :, i] *= 1 / temperature_rho + # N[50:, :] = 0 + temperature[:, :] = 0 + # temperature += 0.01 * np.random.randn(Ny, Nx, NL) + # 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 + N[cylinder] = 0 + N[0, :] = 0 + N[Ny - 1, :] = 0 + + temperature[cylinder] = 0 + # N[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 + # N[:, :, i] *= has_fluid # Prep figure fig = plt.figure(figsize=(4, 2), dpi=80) @@ -66,9 +82,9 @@ def main(): # Drift new_has_fluid = np.zeros((Ny, Nx)) - F_sum = np.sum(F, 2) + F_sum = np.sum(N, 2) for i, cx, cy in zip(idxs, cxs, cys): - F_part = F[:, :, i] / F_sum + F_part = N[:, :, 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)) @@ -76,8 +92,11 @@ def main(): new_has_fluid += to_move - F[:, :, i] = np.roll(F[:, :, i], cx, axis=1) - F[:, :, i] = np.roll(F[:, :, i], cy, axis=0) + N[:, :, i] = np.roll(N[:, :, i], cx, axis=1) + N[:, :, i] = np.roll(N[:, :, i], cy, axis=0) + + temperature[:, :, i] = np.roll(temperature[:, :, i], cx, axis=1) + temperature[:, :, i] = np.roll(temperature[:, :, i], cy, axis=0) # has_fluid = new_has_fluid > 0.5 # new_has_fluid[F_sum == 0] += has_fluid[F_sum == 0] * 1.0 @@ -89,7 +108,7 @@ def main(): print('fluid_cells: %d' % np.sum(has_fluid * 1)) # for i in idxs: - # F[:, :, i] *= has_fluid + # N[:, :, i] *= has_fluid bndry = np.zeros((Ny, Nx), dtype=np.bool) bndry[0, :] = True @@ -101,51 +120,138 @@ def main(): # bndry = np.logical_or(bndry, has_fluid < 0.5) # Set reflective boundaries - bndryF = F[bndry, :] - bndryF = bndryF[:, reflection_mapping] + bndryN = N[bndry, :] + bndryN = bndryN[:, reflection_mapping] - sum_f = np.sum(F) - print('Sum of Forces: %f' % sum_f) + bndryTemp = temperature[bndry, :] + bndryTemp = bndryTemp[:, reflection_mapping] - # sum_f_cyl = np.sum(F[cylinder]) + sum_f = np.sum(N) + print('Sum of Particles: %f' % sum_f) + print('Sum of Temperature: %f' % np.sum(temperature)) + + # sum_f_cyl = np.sum(N[cylinder]) # print('Sum of Forces in cylinder: %f' % sum_f_cyl) - # sum_f_inner_cyl = np.sum(F[inner_cylinder]) + # sum_f_inner_cyl = np.sum(N[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 + # N[Ny - 1, :, 5] += 0.1 + # N[0, :, 1] -= 0.1 + # N[0, :, 5] += 0.1 + # N[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 + rho = np.sum(N, 2) + temp_rho = np.sum(temperature, 2) + ux = np.sum(N * cxs, 2) / rho + uy = np.sum(N * cys, 2) / rho - ux[(np.abs(rho) < 0.000000001)] = 0 - uy[(np.abs(rho) < 0.000000001)] = 0 + ux[(np.abs(rho) < epsilon)] = 0 + uy[(np.abs(rho) < epsilon)] = 0 + + g = -params.g * (temp_rho - yy / Ny) + # uy[np.abs(rho) >= epsilon] += g[np.abs(rho) >= epsilon] / 2.0 + uy += g / 2.0 + + # u_length = np.maximum(np.abs(ux), np.abs(uy)) + u_length = np.sqrt(np.square(ux) + np.square(uy)) + + u_max_length = np.max(u_length) + if u_max_length > np.sqrt(2): + ux = (ux / u_max_length) * np.sqrt(2) + uy = (uy / u_max_length) * np.sqrt(2) + + print('max vector part: %f' % u_max_length) + # ux /= u_max_length + # uy /= u_max_length + + # scale = abs(np.max(np.maximum(np.abs(ux), np.abs(uy))) - 1.0) < epsilon + # if scale: + # g = 0.01 * (temp_rho - yy / Ny) + # + # # F = np.zeros((Ny, Nx), dtype=np.bool) + # # F = -0.1 * rho + # + # # uy[np.abs(rho) >= epsilon] += tau * F[np.abs(rho) >= epsilon] / rho[np.abs(rho) >= epsilon] + # uy[np.abs(rho) >= epsilon] += g[np.abs(rho) >= epsilon] / 2.0 + # + # u_length = np.maximum(np.abs(ux), np.abs(uy)) + # u_max_length = np.max(u_length) + # + # print('max vector part: %f' % u_max_length) + # + # ux /= u_max_length + # uy /= u_max_length # print('minimum rho: %f' % np.min(np.abs(rho))) - # print('Maximum F: %f' % np.max(F)) - # print('Minimum F: %f' % np.min(F)) + # print('Maximum N: %f' % np.max(N)) + # print('Minimum N: %f' % np.min(N)) # Apply Collision - Feq = np.zeros(F.shape) + temperature_eq = np.zeros(temperature.shape) + Neq = np.zeros(N.shape) for i, cx, cy, w in zip(idxs, cxs, cys, weights): - Feq[:, :, i] = rho * w * ( + Neq[:, :, 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) + temperature_eq[:, :, i] = temp_rho * w * ( + 1 + 3 * (cx * ux + cy * uy) + 9 * (cx * ux + cy * uy) ** 2 / 2 - 3 * (ux ** 2 + uy ** 2) / 2) + # test1 = np.sum(Neq) + test2 = np.sum(N-Neq) + if abs(test2) > 0.0001: + test = '' + + print('Overall change: %f' % test2) + + n_pre_sum = np.sum(N[np.logical_not(bndry)]) + temperature_pre_sum = np.sum(temperature[np.logical_not(bndry)]) + + N += -(1.0 / params.t1) * (N - Neq) + temperature += -(1.0 / params.t2) * (temperature - temperature_eq) # Apply boundary - F[bndry, :] = bndryF + N[bndry, :] = bndryN + temperature[bndry, :] = bndryTemp + + # temperature[0, :, 0] = 0 + # temperature[1, :, 0] = 0 + + temperature[0, :, 0] /= 2 + temperature[1, :, 0] /= 2 + + temperature[Ny - 1, :, 0] = 1 + temperature[Ny - 2, :, 0] = 1 + + # n_sum = np.sum(N, 2) + # n_sum_min = np.min(n_sum) + # if n_sum_min < 0: + # N[np.logical_not(bndry)] += abs(n_sum_min) + # N[np.logical_not(bndry)] /= np.sum(N[np.logical_not(bndry)]) + # N[np.logical_not(bndry)] *= n_pre_sum + # print('Sum of Forces: %f' % np.sum(N)) + + # temperature_sum = np.sum(temperature, 2) + # temperature_sum_min = np.min(temperature_sum) + # if temperature_sum_min < 0: + # temperature[np.logical_not(bndry)] += abs(temperature_sum_min) + # temperature[np.logical_not(bndry)] /= np.sum(temperature[np.logical_not(bndry)]) + # temperature[np.logical_not(bndry)] *= temperature_pre_sum + # print('Sum of Temperature: %f' % np.sum(temperature)) + + no_cylinder_mask = np.sum(N, 2) != 0 + print('min N: %f' % np.min(np.sum(N, 2)[no_cylinder_mask])) + print('max N: %f' % np.max(np.sum(N, 2))) + + print('min Temp: %f' % np.min(np.sum(temperature, 2)[no_cylinder_mask])) + print('max Temp: %f' % np.max(np.sum(temperature, 2))) # plot in real time - color 1/2 particles blue, other half red if (plotRealTime and (it % 10) == 0) or (it == Nt - 1): + fig.clear() plt.cla() ux[cylinder] = 0 uy[cylinder] = 0 @@ -158,16 +264,29 @@ def main(): cmap = plt.cm.bwr cmap.set_bad('black') - # plt.imshow(vorticity, cmap='bwr') - plt.imshow(has_fluid * 2.0 - 1.0, cmap='bwr') + plt.subplot(2, 2, 1) + plt.imshow(vorticity, cmap='bwr') + plt.clim(-.1, .1) + # plt.imshow(has_fluid * 2.0 - 1.0, cmap='bwr') # plt.imshow(bndry * 2.0 - 1.0, cmap='bwr') + # plt.imshow(np.sum(N, 2) * 2.0 - 1.0, cmap='bwr') + # plt.imshow((np.sum(temperature, 2) / np.max(np.sum(temperature, 2))) * 2.0 - 1.0, cmap='bwr') + plt.subplot(2, 2, 2) + max_temp = np.max(np.sum(temperature, 2)) + # plt.imshow(np.sum(temperature, 2) / max_temp * 2.0 - 1.0, cmap='bwr') + plt.imshow(np.sum(temperature, 2) * 2.0 - 1.0, cmap='bwr') + plt.clim(-.1, .1) + + plt.subplot(2, 2, 3) + max_N = np.max(np.sum(N, 2)) + plt.imshow(np.sum(N, 2) / max_N * 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') + # 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