{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# The conservative Allen-Cahn model for high Reynolds number, two phase flow with large-density and viscosity constrast" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from pystencils.session import *\n", "from lbmpy.session import *\n", "\n", "from pystencils.simp import sympy_cse\n", "from pystencils.boundaries import BoundaryHandling\n", "\n", "from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle\n", "from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel, CentralMomentMultiphaseForceModel\n", "from lbmpy.phasefield_allen_cahn.kernel_equations import *\n", "from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_parameters_rti\n", "\n", "from lbmpy.advanced_streaming import LBMPeriodicityHandling\n", "from lbmpy.boundaries import NoSlip, LatticeBoltzmannBoundaryHandling" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If `pycuda` is installed the simulation automatically runs on GPU" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "No pycuda installed\n" ] } ], "source": [ "try:\n", " import pycuda\n", "except ImportError:\n", " pycuda = None\n", " gpu = False\n", " target = ps.Target.CPU\n", " print('No pycuda installed')\n", "\n", "if pycuda:\n", " gpu = True\n", " target = ps.Target.GPU" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The conservative Allen-Cahn model (CACM) for two-phase flow is based on the work of Fakhari et al. (2017) [Improved locality of the phase-field lattice-Boltzmann model for immiscible fluids at high density ratios](http://dx.doi.org/10.1103/PhysRevE.96.053301). The model can be created for two-dimensional problems as well as three-dimensional problems, which have been described by Mitchell et al. (2018) [Development of a three-dimensional\n", "phase-field lattice Boltzmann method for the study of immiscible fluids at high density ratios](http://dx.doi.org/10.1103/PhysRevE.96.053301). Furthermore, cascaded lattice Boltzmann methods can be combined with the model which was described in [A cascaded phase-field lattice Boltzmann model for the simulation of incompressible, immiscible fluids with high density contrast](http://dx.doi.org/10.1016/j.camwa.2019.08.018)\n", "\n", "\n", "The CACM is suitable for simulating highly complex two phase flow problems with high-density ratios and high Reynolds numbers. In this tutorial, an overview is provided on how to derive the model with lbmpy. For this, the model is defined with two LBM populations. One for the interface tracking, which we call the phase-field LB step and one for recovering the hydrodynamic properties. The latter is called the hydrodynamic LB step." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Geometry Setup\n", "\n", "First of all, the stencils for the phase-field LB step as well as the stencil for the hydrodynamic LB step are defined. According to the stencils, the simulation can be performed in either 2D- or 3D-space. For 2D simulations, only the D2Q9 stencil is supported. For 3D simulations, the D3Q15, D3Q19 and the D3Q27 stencil are supported. Note here that the cascaded LBM can not be derived for D3Q15 stencils." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "stencil_phase = LBStencil(Stencil.D2Q9)\n", "stencil_hydro = LBStencil(Stencil.D2Q9)\n", "assert(len(stencil_phase[0]) == len(stencil_hydro[0]))\n", "\n", "dimensions = len(stencil_phase[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Definition of the domain size" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# domain \n", "L0 = 256\n", "domain_size = (L0, 4 * L0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Parameter definition\n", "\n", "The next step is to calculate all parameters which are needed for the simulation. In this example, a Rayleigh-Taylor instability test case is set up. The parameter calculation for this setup is already implemented in lbmpy and can be used with the dimensionless parameters which describe the problem." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# time step\n", "timesteps = 8000\n", "\n", "# reference time\n", "reference_time = 4000\n", "# density of the heavier fluid\n", "rho_H = 1.0\n", "\n", "# calculate the parameters for the RTI\n", "parameters = calculate_parameters_rti(reference_length=L0,\n", " reference_time=reference_time,\n", " density_heavy=rho_H,\n", " capillary_number=0.44,\n", " reynolds_number=3000,\n", " atwood_number=0.998,\n", " peclet_number=1000,\n", " density_ratio=1000,\n", " viscosity_ratio=100)\n", "# get the parameters\n", "rho_L = parameters.get(\"density_light\")\n", "\n", "mu_H = parameters.get(\"dynamic_viscosity_heavy\")\n", "mu_L = parameters.get(\"dynamic_viscosity_light\")\n", "\n", "tau_H = parameters.get(\"relaxation_time_heavy\")\n", "tau_L = parameters.get(\"relaxation_time_light\")\n", "\n", "sigma = parameters.get(\"surface_tension\")\n", "M = parameters.get(\"mobility\")\n", "gravitational_acceleration = parameters.get(\"gravitational_acceleration\")\n", "\n", "\n", "drho3 = (rho_H - rho_L)/3\n", "# interface thickness\n", "W = 5\n", "# coeffcient related to surface tension\n", "beta = 12.0 * (sigma/W)\n", "# coeffcient related to surface tension\n", "kappa = 1.5 * sigma*W\n", "# relaxation rate allen cahn (h)\n", "w_c = 1.0/(0.5 + (3.0 * M))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fields\n", "\n", "As a next step all fields which are needed get defined. To do so, we create a `datahandling` object. More details about it can be found in the third tutorial of the [pystencils framework]( http://pycodegen.pages.walberla.net/pystencils/). This object holds all fields and manages the kernel runs." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# create a datahandling object\n", "dh = ps.create_data_handling((domain_size), periodicity=(True, False), parallel=False, default_target=target)\n", "\n", "# pdf fields. g is used for hydrodynamics and h for the interface tracking \n", "g = dh.add_array(\"g\", values_per_cell=len(stencil_hydro))\n", "dh.fill(\"g\", 0.0, ghost_layers=True)\n", "h = dh.add_array(\"h\",values_per_cell=len(stencil_phase))\n", "dh.fill(\"h\", 0.0, ghost_layers=True)\n", "\n", "g_tmp = dh.add_array(\"g_tmp\", values_per_cell=len(stencil_hydro))\n", "dh.fill(\"g_tmp\", 0.0, ghost_layers=True)\n", "h_tmp = dh.add_array(\"h_tmp\",values_per_cell=len(stencil_phase))\n", "dh.fill(\"h_tmp\", 0.0, ghost_layers=True)\n", "\n", "# velocity field\n", "u = dh.add_array(\"u\", values_per_cell=dh.dim)\n", "dh.fill(\"u\", 0.0, ghost_layers=True)\n", "\n", "# phase-field\n", "C = dh.add_array(\"C\")\n", "dh.fill(\"C\", 0.0, ghost_layers=True)\n", "C_tmp = dh.add_array(\"C_tmp\")\n", "dh.fill(\"C_tmp\", 0.0, ghost_layers=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a next step the relaxation time is stated in a symbolic form. It is calculated via interpolation." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# relaxation time and rate\n", "tau = 0.5 + tau_L + (C.center) * (tau_H - tau_L)\n", "s8 = 1/(tau)\n", "\n", "# density for the whole domain\n", "rho = rho_L + (C.center) * (rho_H - rho_L)\n", "\n", "# body force\n", "body_force = [0, 0, 0]\n", "body_force[1] = gravitational_acceleration * rho" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Definition of the lattice Boltzmann methods" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For both LB steps, a weighted orthogonal MRT (WMRT) method is used. It is also possible to change the method to a simpler SRT scheme or a more complicated CLBM scheme. The CLBM scheme can be obtained by commenting in the python snippets in the notebook cells below. Note here that the hydrodynamic LB step is formulated as an incompressible velocity-based LBM. Thus, the velocity terms can not be removed from the equilibrium in the central moment space." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", "
MomentEq. Value Relaxation Rate
$1$$\\rho$$0$
$x$$\\rho u_{0}$$1.8208262344103532$
$y$$\\rho u_{1}$$1.8208262344103532$
$x^{2} - y^{2}$$\\rho u_{0}^{2} - \\rho u_{1}^{2}$$1$
$x y$$\\rho u_{0} u_{1}$$1$
$3 x^{2} + 3 y^{2} - 2$$3 \\rho u_{0}^{2} + 3 \\rho u_{1}^{2}$$1$
$3 x^{2} y - y$$0$$1$
$3 x y^{2} - x$$0$$1$
$9 x^{2} y^{2} - 3 x^{2} - 3 y^{2} + 1$$0$$1$
\n", " " ], "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True,\n", " weighted=True, relaxation_rates=[1, 1, 1, 1, 1])\n", "\n", "method_phase = create_lb_method(lbm_config=config_phase)\n", "method_phase.set_first_moment_relaxation_rate(w_c)\n", "\n", "# config_phase = LBMConfig(stencil=stencil_phase, method=Method.CENTRAL_MOMENT, compressible=True,\n", "# weighted=True, relaxation_rates=[0, w_c, 1, 1, 1], equilibrium_order=4)\n", "\n", "# method_phase = create_lb_method(lbm_config=config_phase)\n", "method_phase" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", "
MomentEq. Value Relaxation Rate
$1$$\\rho$$0$
$x$$u_{0}$$0$
$y$$u_{1}$$0$
$x^{2} - y^{2}$$u_{0}^{2} - u_{1}^{2}$$\\frac{1}{0.664004086170318 - 0.147603677553286 {{C}_{(0,0)}}}$
$x y$$u_{0} u_{1}$$\\frac{1}{0.664004086170318 - 0.147603677553286 {{C}_{(0,0)}}}$
$3 x^{2} + 3 y^{2} - 2$$3 u_{0}^{2} + 3 u_{1}^{2}$$1$
$3 x^{2} y - y$$0$$1$
$3 x y^{2} - x$$0$$1$
$9 x^{2} y^{2} - 3 x^{2} - 3 y^{2} + 1$$0$$1$
\n", " " ], "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, compressible=False,\n", " weighted=True, relaxation_rates=[s8, 1, 1, 1])\n", "\n", "method_hydro = create_lb_method(lbm_config=config_hydro)\n", "\n", "\n", "# config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.CENTRAL_MOMENT, compressible=False,\n", "# weighted=True, relaxation_rates=[s8, 1, 1], equilibrium_order=4)\n", "\n", "# method_hydro = create_lb_method(lbm_config=config_hydro)\n", "\n", "method_hydro" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initialization" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The probability distribution functions (pdfs) are initialised with the equilibrium distribution for the LB methods." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W)\n", "g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)\n", "\n", "h_init = ps.create_kernel(h_updates, target=dh.default_target, cpu_openmp=True).compile()\n", "g_init = ps.create_kernel(g_updates, target=dh.default_target, cpu_openmp=True).compile()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Following this, the phase field is initialised directly in python." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# initialize the domain\n", "def Initialize_distributions():\n", " Nx = domain_size[0]\n", " Ny = domain_size[1]\n", " \n", " for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):\n", " x = np.zeros_like(block.midpoint_arrays[0])\n", " x[:, :] = block.midpoint_arrays[0]\n", " \n", " y = np.zeros_like(block.midpoint_arrays[1])\n", " y[:, :] = block.midpoint_arrays[1]\n", "\n", " y -= 2 * L0\n", " tmp = 0.1 * Nx * np.cos((2 * np.pi * x) / Nx)\n", " init_values = 0.5 + 0.5 * np.tanh((y - tmp) / (W / 2))\n", " block[\"C\"][:, :] = init_values\n", " block[\"C_tmp\"][:, :] = init_values\n", " \n", " if gpu:\n", " dh.all_to_gpu() \n", " \n", " dh.run_kernel(h_init)\n", " dh.run_kernel(g_init)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA68AAAFlCAYAAAAEdN4SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbxklEQVR4nO3de5Dd5X3f8c93VxduBoQRWCAcwFZsgxOTWKbE9njskhScdAJJ7RbPxGEaOnQyOE2bzHRw25l0OkPidtKkSRq7Q2PHJE1CqRPHjBNfiGLHcXzBskOMBZaRAYMsIYmLxE1aac95+scekjVeSbC72n3Yfb1mds45z/n99jxintmzb36/89tqrQUAAAB6NrbYEwAAAICjEa8AAAB0T7wCAADQPfEKAABA98QrAAAA3ROvAAAAdG/FYk/gaE4//fR27rnnLvY0AFgKDn1tsWewNKx89WLPAIAl5Mtf/vLDrbW1R9uu+3g999xzs3nz5sWeBgBLwPChDYs9hSVh7CXelwGYP1X1reeyndOGAQAA6J54BQAAoHtHjdeq+kBV7a6qr00bO62qbquqe0a3a6Y99+6q2lZVW6vqsmnjr62qO0fP/WZV1fz/cwAAAFiKnsuR1w8mufxZY9cn2dRa25Bk0+hxquqCJFcluXC0z3urany0z/uSXJtkw+jr2d8TAAAAZnTUeG2tfSbJo88aviLJTaP7NyW5ctr4za21idbafUm2Jbm4qtYlObm19vnWWkvye9P2AQAAgCOa7Wdez2yt7UyS0e0Zo/Gzkzw4bbvto7GzR/efPT6jqrq2qjZX1eY9e/bMcooAAAAsFfN9waaZPsfajjA+o9baja21ja21jWvXHvXP/QAAALDEzTZed41OBc7odvdofHuSc6Zttz7JjtH4+hnGAQAA4KhmG6+3Jrl6dP/qJB+ZNn5VVa2uqvMydWGm20enFj9RVZeMrjL809P2AQAAgCNacbQNquqPkrw5yelVtT3JLyV5T5JbquqaJA8keXuStNa2VNUtSe5KMpnkutbaYPStfjZTVy4+PsnHRl8AAABwVEeN19baOw7z1KWH2f6GJDfMML45yauf1+wAAAAg83/BJgAAAJh34hUAAIDuiVcAAAC6J14BAADonngFAACge+IVAACA7olXAAAAuideAQAA6J54BQAAoHviFQAAgO6JVwAAALonXgEAAOieeAUAAKB74hUAAIDuiVcAAAC6J14BAADonngFAACge+IVAACA7olXAAAAuideAQAA6J54BQAAoHviFQAAgO6JVwAAALonXgEAAOieeAUAAKB74hUAAIDuiVcAAAC6J14BAADonngFAACge+IVAACA7olXAAAAuideAQAA6J54BQAAoHviFQAAgO6JVwAAALonXgEAAOieeAUAAKB74hUAAIDuiVcAAAC6J14BAADonngFAACge+IVAACA7olXAAAAuideAQAA6J54BQAAoHviFQAAgO6JVwAAALonXgEAAOieeAUAAKB74hUAAIDuiVcAAAC6J14BAADonngFAACge+IVAACA7s0pXqvq31XVlqr6WlX9UVUdV1WnVdVtVXXP6HbNtO3fXVXbqmprVV029+kDAACwHMw6Xqvq7CT/JsnG1tqrk4wnuSrJ9Uk2tdY2JNk0epyqumD0/IVJLk/y3qoan9v0AQAAWA7metrwiiTHV9WKJCck2ZHkiiQ3jZ6/KcmVo/tXJLm5tTbRWrsvybYkF8/x9QEAAFgGZh2vrbVvJ/nVJA8k2ZlkX2vtk0nObK3tHG2zM8kZo13OTvLgtG+xfTQGAAAARzSX04bXZOpo6nlJzkpyYlX91JF2mWGsHeZ7X1tVm6tq8549e2Y7RQAAAJaIuZw2/MNJ7mut7WmtHUryJ0len2RXVa1LktHt7tH225OcM23/9Zk6zfi7tNZubK1tbK1tXLt27RymCAAAwFIwl3h9IMklVXVCVVWSS5PcneTWJFePtrk6yUdG929NclVVra6q85JsSHL7HF4fAACAZWLFbHdsrX2xqj6U5CtJJpP8bZIbk5yU5JaquiZTgfv20fZbquqWJHeNtr+utTaY4/wBAABYBqq1GT922o2NGze2zZs3L/Y0AFgChg9tWOwpLAljL7lnsacAwBJSVV9urW082nZz/VM5AAAAcMyJVwAAALonXgEAAOieeAUAAKB74hUAAIDuiVcAAAC6J14BAADonngFAACge+IVAACA7olXAAAAuideAQAA6J54BQAAoHviFQAAgO6JVwAAALonXgEAAOieeAUAAKB74hUAAIDuiVcAAAC6J14BAADonngFAACge+IVAACA7olXAAAAuideAQAA6J54BQAAoHviFQAAgO6JVwAAALonXgEAAOieeAUAAKB74hUAAIDuiVcAAAC6J14BAADonngFAACge+IVAACA7olXAAAAuideAQAA6J54BQAAoHviFQAAgO6JVwAAALonXgEAAOieeAUAAKB74hUAAIDuiVcAAAC6J14BAADonngFAACge+IVAACA7olXAAAAuideAQAA6J54BQAAoHviFQAAgO6JVwAAALonXgEAAOieeAUAAKB74hUAAIDuiVcAAAC6J14BAADo3pzitapOraoPVdXXq+ruqvqhqjqtqm6rqntGt2umbf/uqtpWVVur6rK5Tx8AAIDlYK5HXn8jycdba69M8pokdye5Psmm1tqGJJtGj1NVFyS5KsmFSS5P8t6qGp/j6wMAALAMzDpeq+rkJG9K8v4kaa0dbK3tTXJFkptGm92U5MrR/SuS3Nxam2it3ZdkW5KLZ/v6AAAALB9zOfJ6fpI9SX63qv62qn6nqk5McmZrbWeSjG7PGG1/dpIHp+2/fTQGAAAARzSXeF2R5AeTvK+19gNJnsroFOHDqBnG2owbVl1bVZuravOePXvmMEUAAACWgrnE6/Yk21trXxw9/lCmYnZXVa1LktHt7mnbnzNt//VJdsz0jVtrN7bWNrbWNq5du3YOUwQAAGApmHW8ttYeSvJgVb1iNHRpkruS3Jrk6tHY1Uk+Mrp/a5Krqmp1VZ2XZEOS22f7+gAAACwfK+a4/88l+YOqWpXk3iT/MlNBfEtVXZPkgSRvT5LW2paquiVTgTuZ5LrW2mCOrw8AAMAyMKd4ba3dkWTjDE9depjtb0hyw1xeEwAAgOVnrn/nFQAAAI458QoAAED3xCsAAADdE68AAAB0T7wCAADQPfEKAABA98QrAAAA3ROvAAAAdE+8AgAA0D3xCgAAQPfEKwAAAN0TrwAAAHRvxWJPAAAWykQ7lCQZy1jGUkmS8fL/cY9k0IZJkmFahpm6f/xiTgiAZUu8ArBsfOjJl+S4OpRTx5/Ki8eeziljh3Lq2FhOGFuZFRkXsiODNsxkBnl6eCh7h8PsG67MI8MTsndwYg60lXnnYk8QgGVJvAKwbPzyB/9FhquTQye2DF58KKef8XguWvvtvOGUe3LR6gezfsVkXjS2almG7DPB+sTwYLZPrsgdE+fkb/ZtyB17zs7Du07O+KMrs/LJytjB5J03LPZsAViOxCsAy8a5v/vNZMWKZOWKDE88PpNrTsyWs74vnzv/Ndn/ygN5/YZ78xOnfyU/uHpHzhxfldW1YslH7KANM9Ems2twMF+ZOCt/vOe1+cI95+f4ratzyr3DnLpjIqc/9mTGntqfHDyUDAaJeAVgEYhXAJaNyV27p+7UWGqsMj4+nlNWrcqaE0/I8Iw1uX/DK/IfX/OqnLpxT37mvM/l0hO+kXVLNGKfidadg4O57alX5IP3/1D2bl6b0/9umFfcsy9ju7dn+NTTaQcPpg0GmRy2ZPT5VwBYDOIVgOWjtdHtYKrDJifTJiYyfPLJ1MOP5EXfXJ1Tbl+TiU+dkd+8+Mr8rzc8kus2/FV+5MRtOXN8dVbXykWd/nyZaIeyazCR2556ef7nN96cyb85LWd+6UBevO3BDB99LG1iIpODwT/89wKADohXAGgtbXIybXIyw/0HsmL3nnzP3afm4OfOyq+98Sfz+2/ZkV8475N5/XF7csrYcVlZ44s941k51AbZNzyQzx1Ym1+997I88ul1Oeuz+7Pqnm9m8NjeTB6aTIaDxZ4mAMxIvALAdMNB2sQgk7t2Z/yxvTl325o8/ZX1+cW3XJ03vvlredeZm/KKlYdyfK16wZxKPGjD7G8Hs/XQWH7jobfmC391Ydb/5aGce+e9GTz6WCYPHnSUFYDuiVcAmElrU6fPPrQrxz22N9/7jTOz9Y4L87bLvjc/97pP5W0nfzXrxo/v/ijsoTbIzsH+fOjx789vfektWffxldlw+44Md+7K5MSEaAXgBUO8AsCRtJbhgQMZfmt7Tn3ksZzy9fX54D++PH962Wvyn1720bxu9b6cPHZcd0dhB22Yx4cH8sWJNfnlb74t+z6+Lq/81N7Uvfdk8smnnB4MwAuOeAWA52I4yODxx1NbtmX9zlOzf8tL87OXX5Mr33R7/tWLP5vzV67s5oJOTw8P5luTk/mdR96YP/2ri/PST0zm7Du2ZfjIoxlOTi729ABgVsQrADwP7dDBDHbvyXF//WRecd+6fObOf5SP/egFuf7VH8+PnHB/Tl/EU4kPtUEeHuzPbU+fm1+586056WMn5RV/vSvt2w9lsH+/U4QBeEETrwDwfLWW4dNPJ9vuz9o9j2bN11+aX/nhf54/ufSb+YVzPpHXrNqfk2r1gp1KPGjDPNkm8ncHj8+vPfjP8o2/eFnO2fRUVtz99Qz2Pe4UYQCWBPEKALM1HGSwd2/Gvvx0zn/g9Dxy57m5+keuzdsu+VLeedrnc/6KHNOrEj9zFeF7J5ObHnlDPvyF12X9bS3nbf5WBnsezsBVhAFYQsQrAMzFM1cl3vFQTvrk3rzq62fl05svyYcvvSg/832fyxUn35HvWVFZXSvn7XTiQ22QiXYo35ps+cjjF+X9X319Ttt0XF75hYfTHtiRyf0HHG0FYMkRrwAwH4aDDJ96KnXPfVn70MN58d+dlf93yaX54BsvyTsu2JwfO/mOnL/iQF40tiorMv68j8YO2jCTGeSJ4cHcO7kqf/b4a/OHW16Xkz57QjZ8cV/G7r0/wyeeSHNBJgCWKPEKAPOoTU5msHdv6qtPZd23XpTh7Wfloz/wpvze616fH7pgW3789Dvyfat3ZO34MCfUeFbWeMYylrHUd3yfYVqGGeZQG+TpNsiewVjunDgrtz58UT5/18tz2pdW5GVfeSJj92/NcN8TGUwecoowAEuaeAWA+dba1FWJH3k0te+JnHnvSTnzb9bmgVd+b/7LBa/KwVfuz/efsz0bT30gLz/uobxkxb6cXBNZWcMkyYE2nqfaqjw0eVq2HXhJNu99ab764Pqsuvv4vPiuQV61dW/y0J4MH39StAKwbIhXADhWnonYRx9L7Xs8L7p/dU75wikZnLEmj6w/Lx86++XZ/5LKxGmDtJMGGVs19TnV4cR46qnxrH5kPMfvannRtwd5+fanMr57V4Z792U4MZE2GIhWAJYV8QoAx1praZOTaZOTGT79dGrX7pywdVVOPP641HHHpR2/Olm9Km3F1Odga3KYTBxM7Z9IO3Agbf+BtIMHMylYAVjGxCsALKRpIZv9+6fGaiw19g+feW1J2rAlbfj3+wDAcideAWCxPBOlbfD3nQoAzOzY/NV0AAAAmEfiFQAAgO6JVwAAALonXgEAAOieeAUAAKB74hUAAIDuiVcAAAC6J14BAADonngFAACge+IVAACA7olXAAAAuideAQAA6J54BQAAoHviFQAAgO6JVwAAALonXgEAAOieeAUAAKB74hUAAIDuiVcAAAC6J14BAADonngFAACge3OO16oar6q/raqPjh6fVlW3VdU9o9s107Z9d1Vtq6qtVXXZXF8bAACA5WE+jrz+fJK7pz2+Psmm1tqGJJtGj1NVFyS5KsmFSS5P8t6qGp+H1wcAAGCJm1O8VtX6JD+W5HemDV+R5KbR/ZuSXDlt/ObW2kRr7b4k25JcPJfXBwAAYHmY65HX/5Hk3ycZThs7s7W2M0lGt2eMxs9O8uC07baPxr5LVV1bVZuravOePXvmOEUAAABe6GYdr1X1T5Psbq19+bnuMsNYm2nD1tqNrbWNrbWNa9eune0UAQAAWCJWzGHfNyT58ar60STHJTm5qv5Pkl1Vta61trOq1iXZPdp+e5Jzpu2/PsmOObw+AAAAy8Ssj7y21t7dWlvfWjs3Uxdi+svW2k8luTXJ1aPNrk7ykdH9W5NcVVWrq+q8JBuS3D7rmQMAALBszOXI6+G8J8ktVXVNkgeSvD1JWmtbquqWJHclmUxyXWttcAxeHwAAgCVmXuK1tfbpJJ8e3X8kyaWH2e6GJDfMx2sCAACwfMzH33kFAACAY0q8AgAA0D3xCgAAQPfEKwAAAN0TrwAAAHRPvAIAANA98QoAAED3xCsAAADdE68AAAB0T7wCAADQPfEKAABA98QrAAAA3ROvAAAAdE+8AgAA0D3xCgAAQPfEKwAAAN0TrwAAAHRPvAIAANA98QoAAED3xCsAAADdE68AAAB0T7wCAADQPfEKAABA98QrAAAA3ROvAAAAdE+8AgAA0D3xCgAAQPfEKwAAAN0TrwAAAHRPvAIAANA98QoAAED3xCsAAADdE68AAAB0T7wCAADQPfEKAABA98QrAAAA3ROvAAAAdE+8AgAA0D3xCgAAQPfEKwAAAN0TrwAAAHRPvAIAANA98QoAAED3xCsAAADdE68AAAB0T7wCAADQPfEKAABA98QrAAAA3ROvAAAAdE+8AgAA0D3xCgAAQPfEKwAAAN0TrwAAAHRPvAIAANC9WcdrVZ1TVZ+qqruraktV/fxo/LSquq2q7hndrpm2z7uraltVba2qy+bjHwAAAMDSN5cjr5NJfrG19qoklyS5rqouSHJ9kk2ttQ1JNo0eZ/TcVUkuTHJ5kvdW1fhcJg8AAMDyMOt4ba3tbK19ZXT/iSR3Jzk7yRVJbhptdlOSK0f3r0hyc2ttorV2X5JtSS6e7esDAACwfMzLZ16r6twkP5Dki0nObK3tTKYCN8kZo83OTvLgtN22j8YAAADgiOYcr1V1UpI/TvJvW2uPH2nTGcbaYb7ntVW1uao279mzZ65TBAAA4AVuTvFaVSszFa5/0Fr7k9HwrqpaN3p+XZLdo/HtSc6Ztvv6JDtm+r6ttRtbaxtbaxvXrl07lykCAACwBMzlasOV5P1J7m6t/dq0p25NcvXo/tVJPjJt/KqqWl1V5yXZkOT22b4+AAAAy8eKOez7hiTvTHJnVd0xGvsPSd6T5JaquibJA0neniSttS1VdUuSuzJ1peLrWmuDObw+AAAAy8Ss47W19tnM/DnWJLn0MPvckOSG2b4mAAAAy9O8XG0YAAAAjiXxCgAAQPfEKwAAAN0TrwAAAHRPvAIAANA98QoAAED3xCsAAADdE68AAAB0T7wCAADQPfEKAABA98QrAAAA3ROvAAAAdE+8AgAA0D3xCgAAQPfEKwAAAN0TrwAAAHRPvAIAANA98QoAAED3xCsAAADdE68AAAB0T7wCAADQPfEKAABA98QrAAAA3ROvAAAAdE+8AgAA0D3xCgAAQPfEKwAAAN0TrwAAAHRPvAIAANA98QoAAED3xCsAAADdE68AAAB0T7wCAADQPfEKAABA98QrAAAA3ROvAAAAdE+8AgAA0D3xCgAAQPfEKwAAAN0TrwAAAHRPvAIAANA98QoAAED3xCsAAADdE68AAAB0T7wCAADQPfEKAABA98QrAAAA3ROvAAAAdE+8AgAA0D3xCgAAQPfEKwAAAN0TrwAAAHRPvAIAANA98QoAAED3xCsAAADdW/B4rarLq2prVW2rqusX+vUBAAB44VnQeK2q8SS/neStSS5I8o6qumAh5wAAAMALz0Ifeb04ybbW2r2ttYNJbk5yxQLPAQAAgBeYhY7Xs5M8OO3x9tEYAAAAHNaKBX69mmGsfddGVdcmuXb08Mmq2npMZ7V8nZ7k4cWeBMuStcdisfbmQdVMb+c8B9Yfi8XaY7E817X3Pc/lmy10vG5Pcs60x+uT7Hj2Rq21G5PcuFCTWq6qanNrbeNiz4Plx9pjsVh7LCbrj8Vi7bFY5nvtLfRpw19KsqGqzquqVUmuSnLrAs8BAACAF5gFPfLaWpusqncl+USS8SQfaK1tWcg5AAAA8MKz0KcNp7X250n+fKFflxk5NZvFYu2xWKw9FpP1x2Kx9lgs87r2qrXvul4SAAAAdGWhP/MKAAAAz5t4XaKq6u1VtaWqhlW18VnPvbuqtlXV1qq6bNr4a6vqztFzv1mjv4VQVaur6v+Oxr9YVecu8D+HF7Cq+s9V9e2qumP09aPTnnteaxHmoqouH621bVV1/WLPh6Wnqu4f/ey6o6o2j8ZOq6rbquqe0e2aadvP+DMQjqaqPlBVu6vqa9PGnvda837LbBxm/S3I73viden6WpKfTPKZ6YNVdUGmrvJ8YZLLk7y3qsZHT78vU39fd8Po6/LR+DVJHmutvTzJryf5r8d89iw1v95au2j09efJrNcizMpobf12krcmuSDJO0ZrEObbW0Y/6575H8fXJ9nUWtuQZNPo8dF+BsLRfDDf/d44m7Xm/ZbZ+GBmXivH/Pc98bpEtdbubq1tneGpK5Lc3FqbaK3dl2Rbkoural2Sk1trn29TH4T+vSRXTtvnptH9DyW51P+ZYx7MZi3CbF2cZFtr7d7W2sEkN2dqDcKxNv099KZ853vrd/0MXPjp8ULUWvtMkkefNfy81pr3W2brMOvvcOZ1/YnX5efsJA9Oe7x9NHb26P6zx79jn9baZJJ9SV58zGfKUvKuqvrq6DSTZ05jms1ahNk63HqD+dSSfLKqvlxV147Gzmyt7UyS0e0Zo3Frkvn2fNea91vm2zH/fU+8voBV1V9U1ddm+DrS0YSZjpi2I4wfaR9IctS1+L4kL0tyUZKdSf77M7vN8K2OthZhtqwrFsIbWms/mKnT06+rqjcdYVtrkoXi/ZaFsCC/7y3433ll/rTWfngWu21Pcs60x+uT7BiNr59hfPo+26tqRZJT8txPFWAZeK5rsar+d5KPjh7OZi3CbB1uvcG8aa3tGN3urqoPZ+o04F1Vta61tnN0mtzu0ebWJPPt+a4177fMm9barmfuH8vf9xx5XX5uTXLV6ArC52Xqw9G3j04veaKqLhl9nvWnk3xk2j5Xj+6/LclfNn8gmOdo9Ab6jJ/I1MXEktmtRZitLyXZUFXnVdWqTF084tZFnhNLSFWdWFUveuZ+kn+SqZ93099Dr853vrd+18/AhZ01S8zzWmveb5lPC/X7niOvS1RV/USS30qyNsmfVdUdrbXLWmtbquqWJHclmUxyXWttMNrtZzN19bDjk3xs9JUk70/y+1W1LVNHXK9auH8JS8B/q6qLMnUqyP1J/nWSzHItwqy01iar6l1JPpFkPMkHWmtbFnlaLC1nJvnw6HqGK5L8YWvt41X1pSS3VNU1SR5I8vbkqD8D4Yiq6o+SvDnJ6VW1PckvJXlPnv9a837L83aY9ffmhfh9rxxAAwAAoHdOGwYAAKB74hUAAIDuiVcAAAC6J14BAADonngFAACge+IVAACA7olXAAAAuideAQAA6N7/B14Dy2PIVC98AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "Initialize_distributions()\n", "plt.scalar_field(dh.gather_array(C.name))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Source Terms" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the Allen-Cahn LB step, the Allen-Cahn equation needs to be applied as a source term. Here, a simple forcing model is used which is directly applied in the moment space: \n", "$$\n", "F_i^\\phi (\\boldsymbol{x}, t) = \\Delta t \\frac{\\left[1 - 4 \\left(\\phi - \\phi_0\\right)^2\\right]}{\\xi} w_i \\boldsymbol{c}_i \\cdot \\frac{\\nabla \\phi}{|{\\nabla \\phi}|},\n", "$$\n", "where $\\phi$ is the phase-field, $\\phi_0$ is the interface location, $\\Delta t$ it the timestep size $\\xi$ is the interface width, $\\boldsymbol{c}_i$ is the discrete direction from stencil_phase and $w_i$ are the weights. Furthermore, the equilibrium needs to be shifted:\n", "\n", "$$\n", "\\bar{h}^{eq}_\\alpha = h^{eq}_\\alpha - \\frac{1}{2} F^\\phi_\\alpha\n", "$$\n", "\n", "The hydrodynamic force is given by:\n", "$$\n", "F_i (\\boldsymbol{x}, t) = \\Delta t w_i \\frac{\\boldsymbol{c}_i \\boldsymbol{F}}{\\rho c_s^2},\n", "$$\n", "where $\\rho$ is the interpolated density and $\\boldsymbol{F}$ is the source term which consists of the pressure force \n", "$$\n", "\\boldsymbol{F}_p = -p^* c_s^2 \\nabla \\rho,\n", "$$\n", "the surface tension force:\n", "$$\n", "\\boldsymbol{F}_s = \\mu_\\phi \\nabla \\phi\n", "$$\n", "and the viscous force term:\n", "$$\n", "F_{\\mu, i}^{\\mathrm{MRT}} = - \\frac{\\nu}{c_s^2 \\Delta t} \\left[\\sum_{\\beta} c_{\\beta i} c_{\\beta j} \\times \\sum_{\\alpha} \\Omega_{\\beta \\alpha}(g_\\alpha - g_\\alpha^{\\mathrm{eq}})\\right] \\frac{\\partial \\rho}{\\partial x_j}.\n", "$$\n", "\n", "In the above equations $p^*$ is the normalised pressure which can be obtained from the zeroth order moment of the hydrodynamic distribution function $g$. The lattice speed of sound is given with $c_s$ and the chemical potential is $\\mu_\\phi$. Furthermore, the viscosity is $\\nu$ and $\\Omega$ is the moment-based collision operator. Note here that the hydrodynamic equilibrium is also adjusted as shown above for the phase-field distribution functions.\n", "\n", "\n", "For CLBM methods the forcing is applied directly in the central moment space. This is done with the `CentralMomentMultiphaseForceModel`. Furthermore, the GUO force model is applied here to be consistent with [A cascaded phase-field lattice Boltzmann model for the simulation of incompressible, immiscible fluids with high density contrast](http://dx.doi.org/10.1016/j.camwa.2019.08.018). Here we refer to equation D.7 which can be derived for 3D stencils automatically with lbmpy." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W)]\n", "force_g = hydrodynamic_force(g, C, method_hydro, tau, rho_H, rho_L, kappa, beta, body_force)\n", "\n", "\n", "if isinstance(method_phase, CentralMomentBasedLbMethod):\n", " force_model_h = CentralMomentMultiphaseForceModel(force=force_h)\n", "else:\n", " force_model_h = MultiphaseForceModel(force=force_h)\n", "\n", "\n", "if isinstance(method_hydro, CentralMomentBasedLbMethod):\n", " force_model_g = CentralMomentMultiphaseForceModel(force=force_g, rho=rho)\n", "else:\n", " force_model_g = MultiphaseForceModel(force=force_g, rho=rho)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Definition of the LB update rules" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The update rule for the phase-field LB step is defined as:\n", "\n", "$$\n", "h_i (\\boldsymbol{x} + \\boldsymbol{c}_i \\Delta t, t + \\Delta t) = h_i(\\boldsymbol{x}, t) + \\Omega_{ij}^h(\\bar{h_j}^{eq} - h_j)|_{(\\boldsymbol{x}, t)} + F_i^\\phi(\\boldsymbol{x}, t).\n", "$$\n", "In our framework the pull scheme is applied as streaming step. Furthermore, the update of the phase-field is directly integrated into the kernel. As a result of this, a second temporary phase-field is needed." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "allen_cahn_lb = get_collision_assignments_phase(lb_method=method_phase,\n", " velocity_input=u,\n", " output={'density': C_tmp},\n", " force_model=force_model_h,\n", " symbolic_fields={\"symbolic_field\": h,\n", " \"symbolic_temporary_field\": h_tmp},\n", " kernel_type='stream_pull_collide')\n", "\n", "# allen_cahn_lb = sympy_cse(allen_cahn_lb)\n", "\n", "ast_allen_cahn_lb = ps.create_kernel(allen_cahn_lb, target=dh.default_target, cpu_openmp=True)\n", "kernel_allen_cahn_lb = ast_allen_cahn_lb.compile()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The update rule for the hydrodynmaic LB step is defined as:\n", "\n", "$$\n", "g_i (\\boldsymbol{x} + \\boldsymbol{c}_i \\Delta t, t + \\Delta t) = g_i(\\boldsymbol{x}, t) + \\Omega_{ij}^g(\\bar{g_j}^{eq} - g_j)|_{(\\boldsymbol{x}, t)} + F_i(\\boldsymbol{x}, t).\n", "$$\n", "\n", "Here, the push scheme is applied which is easier due to the data access required for the viscous force term. Furthermore, the velocity update is directly done in the kernel." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "scrolled": false }, "outputs": [], "source": [ "hydro_lb_update_rule = get_collision_assignments_hydro(lb_method=method_hydro,\n", " density=rho,\n", " velocity_input=u,\n", " force_model=force_model_g,\n", " sub_iterations=2,\n", " symbolic_fields={\"symbolic_field\": g,\n", " \"symbolic_temporary_field\": g_tmp},\n", " kernel_type='collide_stream_push')\n", "\n", "# hydro_lb_update_rule = sympy_cse(hydro_lb_update_rule)\n", "\n", "ast_hydro_lb = ps.create_kernel(hydro_lb_update_rule, target=dh.default_target, cpu_openmp=True)\n", "kernel_hydro_lb = ast_hydro_lb.compile()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Boundary Conditions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a last step suitable boundary conditions are applied" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "# periodic Boundarys for g, h and C\n", "periodic_BC_C = dh.synchronization_function(C.name, target=dh.default_target, optimization = {\"openmp\": True})\n", "\n", "periodic_BC_g = LBMPeriodicityHandling(stencil=stencil_hydro, data_handling=dh, pdf_field_name=g.name,\n", " streaming_pattern='push')\n", "periodic_BC_h = LBMPeriodicityHandling(stencil=stencil_phase, data_handling=dh, pdf_field_name=h.name,\n", " streaming_pattern='pull')\n", "\n", "# No slip boundary for the phasefield lbm\n", "bh_allen_cahn = LatticeBoltzmannBoundaryHandling(method_phase, dh, 'h',\n", " target=dh.default_target, name='boundary_handling_h',\n", " streaming_pattern='pull')\n", "\n", "# No slip boundary for the velocityfield lbm\n", "bh_hydro = LatticeBoltzmannBoundaryHandling(method_hydro, dh, 'g' ,\n", " target=dh.default_target, name='boundary_handling_g',\n", " streaming_pattern='push')\n", "\n", "contact_angle = BoundaryHandling(dh, C.name, stencil_hydro, target=dh.default_target)\n", "contact = ContactAngle(90, W)\n", "\n", "wall = NoSlip()\n", "if dimensions == 2:\n", " bh_allen_cahn.set_boundary(wall, make_slice[:, 0])\n", " bh_allen_cahn.set_boundary(wall, make_slice[:, -1])\n", "\n", " bh_hydro.set_boundary(wall, make_slice[:, 0])\n", " bh_hydro.set_boundary(wall, make_slice[:, -1])\n", " \n", " contact_angle.set_boundary(contact, make_slice[:, 0])\n", " contact_angle.set_boundary(contact, make_slice[:, -1])\n", "else:\n", " bh_allen_cahn.set_boundary(wall, make_slice[:, 0, :])\n", " bh_allen_cahn.set_boundary(wall, make_slice[:, -1, :])\n", "\n", " bh_hydro.set_boundary(wall, make_slice[:, 0, :])\n", " bh_hydro.set_boundary(wall, make_slice[:, -1, :])\n", " \n", " contact_angle.set_boundary(contact, make_slice[:, 0, :])\n", " contact_angle.set_boundary(contact, make_slice[:, -1, :])\n", "\n", "\n", "bh_allen_cahn.prepare()\n", "bh_hydro.prepare()\n", "contact_angle.prepare()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Full timestep" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# definition of the timestep for the immiscible fluids model\n", "def timeloop():\n", " # Solve the interface tracking LB step with boundary conditions\n", " periodic_BC_h()\n", " bh_allen_cahn() \n", " dh.run_kernel(kernel_allen_cahn_lb)\n", " dh.swap(\"C\", \"C_tmp\")\n", " \n", " # apply the three phase-phase contact angle\n", " contact_angle()\n", " # periodic BC of the phase-field\n", " periodic_BC_C()\n", " \n", " # solve the hydro LB step with boundary conditions\n", " dh.run_kernel(kernel_hydro_lb)\n", " periodic_BC_g()\n", " bh_hydro()\n", "\n", " \n", " # field swaps\n", " dh.swap(\"h\", \"h_tmp\")\n", " dh.swap(\"g\", \"g_tmp\")" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Initialize_distributions()\n", "\n", "frames = 300\n", "steps_per_frame = (timesteps//frames) + 1\n", "\n", "if 'is_test_run' not in globals():\n", " def run():\n", " for i in range(steps_per_frame):\n", " timeloop()\n", " \n", " if gpu:\n", " dh.to_cpu(\"C\")\n", " return dh.gather_array(C.name)\n", "\n", " animation = plt.scalar_field_animation(run, frames=frames, rescale=True)\n", " set_display_mode('video')\n", " res = display_animation(animation)\n", "else:\n", " timeloop()\n", " res = None\n", "res" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the video is played for 10 seconds while the simulation time is only 2 seconds!" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.7" } }, "nbformat": 4, "nbformat_minor": 2 }