{ "cells": [ { "cell_type": "markdown", "id": "1d2eb3be-2943-4811-967e-c481ad04bc2b", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "# Convective Equation (Transport Equation)\n", "\n", "The convective equation also known as transport equation is a partial diferential equation (PDE) usually described by\n", "\n", "$$\n", "\\partial_t u + \\partial_{\\alpha}B_{\\alpha} = 0\n", "$$(cv-eq-lb)\n", "\n", "where $u(\\alpha,t)$ is a time and space dependent parameter that is being transporte and $B_{\\alpha}(\\alpha,t)$ is a second time and space dependent parameter that can also be depende of $u(\\alpha,t)$ in non-linear convective diffusive equations. " ] }, { "cell_type": "markdown", "id": "37b0f5cd-fa20-4a0b-b659-411688021062", "metadata": {}, "source": [ "## Numerical Discretization" ] }, { "cell_type": "markdown", "id": "c936a7d0-d22c-497b-b51a-53beae8b9663", "metadata": {}, "source": [ "Discretizing the above equation in regular grid of 1D domain, we have" ] }, { "cell_type": "markdown", "id": "3a087e51-4f2a-4877-832d-6dc8e8380489", "metadata": {}, "source": [ "```{figure} mesh-1-convective.svg\n", "---\n", "scale: 250%\n", "align: center\n", "name: mesh-1-convective\n", "---\n", "1D regular grid.\n", "```" ] }, { "cell_type": "markdown", "id": "1da843b8-89a4-4607-b5bb-4a33e77184b3", "metadata": {}, "source": [ "### FDM Discretization\n", "\n", "The Discretization through finite difference method can assume distinct forms that describe the problem with different accuracy, rate convergence and stability. In the sequence of this section, we will explore different variation that can work with the convective formulation.\n", "\n", "#### Forward-Time and Forward-Space (FTFS)\n", "Discretizing the Eq. {eq}`cv-eq-lb` through finite diference method, we have for forward-time and forward-space (FTFS) discretizations:\n", "\n", "$$\n", "\\partial_t u + \\partial_{x}B_{\\alpha} = \\frac{u(x,t + \\Delta_{t}) - u(x,t) }{\\Delta_{t}} + \\frac{B_{\\alpha}(x,t) - B_{\\alpha}(x + \\Delta_{x},t) }{\\Delta_{x}} =0,\n", "$$\n", "\n", "isolating the term $u(x,t + \\Delta_{t})$, we can describe the $u$ time evolution:\n", "\n", "$$\n", "u(x,t + \\Delta_{t}) = u(x,t) - \\frac{\\Delta_{t}}{\\Delta_{x}} \\left( B_{\\alpha}(x,t) - B_{\\alpha}(x + \\Delta_{x},t) \\right) .\n", "$$\n", "\n", "#### Forward-Time and Central-Space with Small Diffusivity (FTCS-SD)\n", "\n", "Give that discretize convective equation with pure forwar-time and central-space (FTCS) can be unstable for several cases, employ a forwar-time and central-space with small diffusivity (FTCS-SD) is a alternative to improve stability and ensure a propertie conservation.\n", "\n", "$$\n", "\\partial_t u + \\partial_{x}B_{\\alpha} - \\partial_{\\alpha}\\left( \\nu \\partial_{\\alpha}u \\right) = \\frac{u(x,t + \\Delta_{t}) - u(x,t) }{\\Delta_{t}} + \\frac{B_{\\alpha}(x+\\Delta_{x},t) - B_{\\alpha}(x - \\Delta_{x},t) }{2\\Delta_{x}} - \\nu\\frac{ u(x - \\Delta_{x},t) -2u(x,t) + u(x + \\Delta_{x},t) }{\\Delta_{x}^{2}} =0,\n", "$$\n", "\n", "isolating the term $u(x,t + \\Delta_{t})$, we can describe the $u$ time evolution:\n", "\n", "$$\n", "u(x,t + \\Delta_{t}) = u(x,t) + \\frac{\\Delta_{t}}{2\\Delta_{x}} \\left( B_{\\alpha}(x-\\Delta_{x},t) - B_{\\alpha}(x + \\Delta_{x},t) \\right) + \\frac{\\nu\\Delta_{t}}{\\Delta_{x}^{2}} \\left( u(x - \\Delta_{x},t) -2u(x,t) + u(x + \\Delta_{x},t) \\right) .\n", "$$" ] }, { "cell_type": "markdown", "id": "320b11c7-d9a3-4634-bb62-1388136ac2be", "metadata": {}, "source": [ "### Lattice Boltzmann Equation\n", "\n", "Describing the problem through the BGK lattice Boltzmann equation, we assume a equilibrium state for the present problem, where $ f_{i}= f_{i}^{eq}$ and the BGK lattice Boltzmann equation reduces to:\n", "\n", "$$\n", "f_{i}\\left(x_{\\alpha} + c_{i,\\alpha} \\Delta t, t+\\Delta t\\right)- f_{i}(x_{\\alpha}, t)= -\\left( \\frac{ f_{i} - f_{i}^{eq} }{ \\tau } \\right) \\quad \\quad \\rightarrow \\quad \\quad f_i( x_{\\alpha} + c_{i,\\alpha} \\delta t, t+\\delta t) = f_{i}^{eq} , \n", "$$(LB-cv-Eq)\n", "\n", "where the equilibrium distribution function is defined by\n", "\n", "$$\n", "f^{eq}_{i} = w_{i}\\left( u + \\frac{c_{i,\\alpha} B_{\\alpha}}{c_{s}^{2}} \\right).\n", "$$" ] }, { "cell_type": "markdown", "id": "a1f0e00d-c032-414a-8aff-b666ba180f7d", "metadata": {}, "source": [ "#### Lattice Direction Moments" ] }, { "cell_type": "code", "execution_count": 2, "id": "c8dd6283-a6ef-4060-8675-30eb45b800cf", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "import warnings\n", "warnings.filterwarnings(\"ignore\")\n", "from pylab import *\n", "from __future__ import division\n", "from sympy import *\n", "import numpy as np\n", "from sympy import S, collect, expand, factor, Wild\n", "from sympy import fraction, Rational, Symbol\n", "from sympy import symbols, sqrt, Rational\n", "import sympy as sp\n", "from IPython.display import display, Math, Latex\n", "#-------------------------------------------------Símbolos----------------------------------------------\n", "omega, u, B, w = symbols('omega, u, B_{\\\\alpha}, w')\n", "wi, cx, cy, cs = symbols('w_{i} c_{x} c_{y} c_{s}')\n", "fi, f0, f1, f2 = symbols('f_{i} f_{0} f_{1} f_{2}')\n", "#-------------------------------------------------Funções----------------------------------------------\n", "feq = Function('feq')(wi, cx, cy)\n", "fneq = Function('fneq')(wi, cx, cy)\n", "f = Function('f')(fi)\n", "#----------------------------------------------Lattice-D2Q9---Variáveis----------------------------------------------\n", "fi=np.array([f0,f1,f2])\n", "w0=Rational(4,6);w1=Rational(1,6)\n", "wi=np.array([w0,w1,w1])\n", "cx=np.array([0,1,-1])\n", "as2=Rational(3)\n", "cs2=1/as2\n", "#-------------------------------------------------Calc.Func------------------------------------------------\n", "f= fi\n", "feq=wi*(u + cx*B/cs2)" ] }, { "cell_type": "code", "execution_count": 3, "id": "2147848c-6015-4c81-8eee-6f09e6191c81", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "hide-input" ] }, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\underbrace{\\sum_{i=0} f_{i}^{eq} =\\sum_{i=0} f_{i} }_{\\textrm{Zero-Order Moment}} =u,\\quad \\quad \\underbrace{\\sum_{i=0} f_{i}^{eq}c_{i,x} }_{\\textrm{x-First-Order Moment}} =B_{\\alpha},\\quad \\quad \\underbrace{\\sum_{i=0} f_{i}^{eq} c_{i,x}c_{i,x} }_{\\textrm{xx-Second-Order Moment}} =\\frac{u}{3},\\quad \\quad \\underbrace{\\sum_{i=0} f_{i}^{eq} c_{i,x}c_{i,x}c_{i,x} }_{\\textrm{xxx-Third-Order Moment}} =B_{\\alpha}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle \\underbrace{\\sum_{i=0} f_{i}^{eq} =\\sum_{i=0} f_{i} }_{\\textrm{Zero-Order Hermite Moment}} =u,\\quad \\quad \\underbrace{\\sum_{i=0} f_{i}^{eq}c_{i,x} }_{\\textrm{x-First-Order Hermite Moment}} =B_{\\alpha},\\quad \\quad \\underbrace{\\sum_{i=0} f_{i}^{eq} \\left(c_{i,x}c_{i,x} - \\frac{1}{c_{s}^{2}}\\right) }_{\\textrm{xx-Second-Order Hermite Moment}} =0,\\quad \\quad \\underbrace{\\sum_{i=0} f_{i}^{eq} \\left(c_{i,x}c_{i,x} - \\frac{3}{c_{s}^{2}}\\right) c_{i,x} }_{\\textrm{xxx-Third-Order Hermite Moment}} =0$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "a0=simplify(sum(feq))\n", "ax=simplify(sum(feq*cx))\n", "axx=simplify(sum(feq*cx*cx))\n", "axxx=simplify(sum(feq*cx*cx*cx))\n", "display(Math(r\"\\underbrace{\\sum_{i=0} f_{i}^{eq} =\\sum_{i=0} f_{i} }_{\\textrm{Zero-Order Moment}} =\" + sp.latex(a0) \n", " +r\",\\quad \\quad \\underbrace{\\sum_{i=0} f_{i}^{eq}c_{i,x} }_{\\textrm{x-First-Order Moment}} =\" + sp.latex(ax)\n", " +r\",\\quad \\quad \\underbrace{\\sum_{i=0} f_{i}^{eq} c_{i,x}c_{i,x} }_{\\textrm{xx-Second-Order Moment}} =\" + sp.latex(axx)\n", " +r\",\\quad \\quad \\underbrace{\\sum_{i=0} f_{i}^{eq} c_{i,x}c_{i,x}c_{i,x} }_{\\textrm{xxx-Third-Order Moment}} =\" + sp.latex(axxx) ))\n", "\n", "aH0=simplify(sum(feq))\n", "aHx=simplify(sum(feq*cx))\n", "aHxx=simplify(sum(feq*(cx*cx-cs2)))\n", "aHxxx=simplify(sum(feq*(cx*cx-3*cs2)*cx))\n", "display(Math(r\"\\underbrace{\\sum_{i=0} f_{i}^{eq} =\\sum_{i=0} f_{i} }_{\\textrm{Zero-Order Hermite Moment}} =\" + sp.latex(aH0) \n", " +r\",\\quad \\quad \\underbrace{\\sum_{i=0} f_{i}^{eq}c_{i,x} }_{\\textrm{x-First-Order Hermite Moment}} =\" + sp.latex(aHx)\n", " +r\",\\quad \\quad \\underbrace{\\sum_{i=0} f_{i}^{eq} \\left(c_{i,x}c_{i,x} - \\frac{1}{c_{s}^{2}}\\right) }_{\\textrm{xx-Second-Order Hermite Moment}} =\" + sp.latex(aHxx)\n", " +r\",\\quad \\quad \\underbrace{\\sum_{i=0} f_{i}^{eq} \\left(c_{i,x}c_{i,x} - \\frac{3}{c_{s}^{2}}\\right) c_{i,x} }_{\\textrm{xxx-Third-Order Hermite Moment}} =\" + sp.latex(aHxxx) ))" ] }, { "cell_type": "markdown", "id": "16f5cc23-a6a1-4fe4-a8b6-1488556126a7", "metadata": {}, "source": [ "#### Chapmann-Enskog Analysis (Considering Equilibrium State - Consequently Disregarding Inherited Diffusive Effects)\n", "\n", "```{toggle}\n", "Applying the Chapmann-Enskog procedure to LB equation, we expand the term $f_{i}\\left(\\boldsymbol{x}+\\boldsymbol{e}_i \\Delta t, t+\\Delta t\\right)$ in a Taylor series to recover the derivative form of the equation, i.e.,\n", "\n", "$$\n", "f_{i}\\left(x_{\\alpha} + c_{i,\\alpha} \\Delta t, t+\\Delta t\\right)- f_{i}(\\boldsymbol{x}, t)=\\displaystyle\\sum_{j=1}^{\\infty}\\frac{\\Delta t^{j}}{j!}D_{t}^{j}f_{i}= -\\left( \\frac{ f_{i} - f_{i}^{eq} }{ \\tau } \\right) \\quad \\quad \\overbrace{\\rightarrow}^{\\textrm{Equilibrium State}} \\quad \\quad f_{i}\\left(x_{\\alpha} + c_{i,\\alpha} \\Delta t, t+\\Delta t\\right)=\\displaystyle\\sum_{j=0}^{\\infty}\\frac{\\Delta t^{j}}{j!}D_{t}^{j}f_{i}= f_{i}^{eq}(\\boldsymbol{x}, t).\n", "$$(EqExp-cv-Eq)\n", "\n", "Rescaling the dimensionless form of the Eq. {eq}`EqExp-cv-Eq` in terms of the Knudsen number ($Kn$), we have\n", "\n", "$$\n", "\\displaystyle \\sum^{\\infty}_{j=0} \\frac{Kn^{(j-1)}}{j!} D_{t}^{j} f_{i} = f_{i}^{eq}(x_{\\alpha}, t),\n", "$$\n", "\n", "applying the asymptotic expansion in both the distribution function ($f_{i}=f_{i}^{(0)}+Kn f_{i}^{(1)} + Kn^{2} f_{i}^{(2)}+\\cdots$) and time partial derivative ($\\partial_{t}=\\partial_{t}^{(0)}+ Kn \\partial_{t}^{(1)}+Kn^{2} \\partial_{t}^{(2)}+\\cdots$), and separating the equation in orders up to the order $Kn^{2}$:\n", "\n", "$$\n", "\\begin{array}{ll}\n", " (Kn^{(0)}):& f_{i}^{(0)}= f_{i}^{eq},\\\\\n", " (Kn^{(1)}):& \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(0)}= 0, \\\\\n", " (Kn^{(2)}):& \\partial_{t}^{(1)} f_{i}^{(0)} + \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(1)} + \\displaystyle\\frac{\\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)^{2} f_{i}^{(0)}}{2}= 0,\\\\\n", " \\textrm{ou}\\\\\n", " (Kn^{(2)}):& \\partial_{t}^{(1)} f_{i}^{(0)} + \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(1)} = 0,\\\\\n", " (Kn^{(3)}):& \\partial_{t}^{(2)} f_{i}^{(0)} +\\partial_{t}^{(1)} f_{i}^{(1)} + \\partial_{t}^{(1)}\\left( \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(0)} \\right) + \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(2)} + \\displaystyle\\frac{\\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)^{2} f_{i}^{(1)}}{2} + \\displaystyle\\frac{\\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)^{3} f_{i}^{(0)}}{6} = 0, \\\\\n", " \\textrm{ou}\\\\\n", " (Kn^{(3)}):& \\partial_{t}^{(2)} f_{i}^{(0)} +\\partial_{t}^{(1)} f_{i}^{(1)} + \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(2)} = 0. \\\\\n", " \\end{array}\n", "$$(Chap-Kn-cv-Eq)\n", "\n", "##### Zero-Order Moment Balance\n", "\n", "To retrieve the balance equation, we sum the Eq. {eq}`Chap-Kn-cv-Eq` for $Kn^{(1)}$ and $Kn^{(2)}$ over $\\sum_{i=0} $:\n", "\n", "$$\n", "\\begin{array}{l}\n", "(Kn^{(1)}): \\displaystyle\\sum_{i}\\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(0)}= 0,\\\\\n", "(Kn^{(1)}): \\partial_{t}^{(0)} u + \\partial_{\\alpha} B_{\\alpha} =0,\n", "\\end{array}\n", "$$(cv-Kn1-B0)\n", "\n", "$$\n", "\\begin{array}{l}\n", "(Kn^{(2)}):& \\displaystyle\\sum_{i}\\left[ \\partial_{t}^{(1)} f_{i}^{(0)} + \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right) f_{i}^{(1)} \\right]= 0, \\nonumber \\\\\n", "(Kn^{(2)}):& \\partial_{t}^{(1)} u + \\partial_{\\alpha} m^{(1)}_{\\alpha} = 0,\n", "\\end{array}\n", "$$(cv-Kn2-B0)\n", "\n", "and\n", "\n", "$$\n", "\\begin{array}{l}\n", "(Kn^{(3)}):& \\displaystyle\\sum_{i}\\left[ \\partial_{t}^{(2)} f_{i}^{(0)} +\\partial_{t}^{(1)} f_{i}^{(1)} + \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(2)} \\right]= 0, \\nonumber \\\\\n", "(Kn^{(3)}):& \\partial_{t}^{(2)} u + \\partial_{\\alpha} m^{(2)}_{\\alpha} = 0.\n", "\\end{array}\n", "$$(cv-Kn3-B0)\n", "\n", "By summing the Eqs. {eq}`cv-Kn1-B0`, {eq}`cv-Kn2-B0`, and {eq}`cv-Kn3-B0`, and also extending to higher orders of Knudsen, we have\n", "\n", "$$\n", "\\partial_{t} u + \\partial_{\\alpha} B_{\\alpha} + \\underbrace{\\partial_{\\alpha} \\left( m^{(1)}_{\\alpha} + m^{(2)}_{\\alpha} + ...\\right) }_{\\textrm{Numerical Errors}} =0,\n", "$$\n", "\n", "where these numerical error can be related to the numerical error from the approximation of a derivative term by an expansion in Taylor series.\n", "```" ] }, { "cell_type": "markdown", "id": "e76ff956-7520-4d95-bc10-91674c175912", "metadata": {}, "source": [ "#### Chapmann-Enskog Analysis (Full)\n", "\n", "```{toggle}\n", "Applying the Chapmann-Enskog procedure to LB equation, we expand the term $f_{i}\\left(\\boldsymbol{x}+\\boldsymbol{e}_i \\Delta t, t+\\Delta t\\right)$ in a Taylor series to recover the derivative form of the equation, i.e.,\n", "\n", "$$\n", "f_{i}\\left(x_{\\alpha} + c_{i,\\alpha} \\Delta t, t+\\Delta t\\right)- f_{i}(x_{\\alpha}, t)=\\displaystyle\\sum_{j=1}^{\\infty}\\frac{\\Delta t^{j}}{j!}D_{t}^{j}f_{i}= -\\left( \\frac{ f_{i} - f_{i}^{eq} }{ \\tau } \\right) .\n", "$$(EqExp-cv-NotEq)\n", "\n", "Rescaling the dimensionless form of the Eq. {eq}`EqExp-cv-Eq` in terms of the Knudsen number ($Kn$), we have\n", "\n", "$$\n", "\\displaystyle \\sum^{\\infty}_{j=1} \\frac{Kn^{(j-1)}}{j!} D_{t}^{j} f_{i} = - \\frac{1}{Kn}\\left( \\frac{ f_{i} - f_{eq,i} }{ \\tau } \\right) \\quad \\quad \\rightarrow \\quad \\quad \\displaystyle \\sum^{\\infty}_{j=1} \\frac{Kn^{(j)}}{j!} D_{t}^{j} f_{i} = - \\frac{ f_{i} - f_{i}^{eq} }{ \\tau },\n", "$$\n", "\n", "applying the asymptotic expansion in both the distribution function ($f_{i}=f_{i}^{(0)}+Kn f_{i}^{(1)} + Kn^{2} f_{i}^{(2)}+\\cdots$) and time partial derivative ($\\partial_{t}=\\partial_{t}^{(0)}+ Kn \\partial_{t}^{(1)}+Kn^{2} \\partial_{t}^{(2)}+\\cdots$), and separating the equation in orders up to the order $Kn^{2}$:\n", "\n", "$$\n", "\\begin{array}{ll}\n", " (Kn^{(0)}):& f_{i}^{(0)}= f_{i}^{eq},\\\\\n", " (Kn^{(1)}):& \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(0)}= - \\displaystyle\\frac{ f_{i}^{(1)} }{ \\tau } , \\\\\n", " (Kn^{(2)}):& \\partial_{t}^{(1)} f_{i}^{(0)} + \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(1)} + \\displaystyle\\frac{\\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)^{2} f_{i}^{(0)}}{2}= - \\displaystyle\\frac{ f_{i}^{(2)} }{ \\tau },\\\\\n", " \\textrm{ou}\\\\\n", " (Kn^{(2)}):& \\partial_{t}^{(1)} f_{i}^{(0)} + \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)\\displaystyle\\left(1 - \\frac{1}{2\\tau}\\right) f_{i}^{(1)} = - \\displaystyle\\frac{ f_{i}^{(2)} }{ \\tau } \\quad\\quad \\rightarrow \\quad\\quad f_{i}^{(1)}\\textrm{ Formulation} \\\\\n", " \\textrm{ou}\\\\\n", " (Kn^{(2)}):& \\partial_{t}^{(1)} f_{i}^{(0)} + \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)^{2} \\displaystyle\\left(\\frac{1}{2} - \\tau\\right) f_{i}^{(0)} = - \\displaystyle\\frac{ f_{i}^{(2)} }{ \\tau } \\quad\\quad \\rightarrow \\quad\\quad f_{i}^{(0)}\\textrm{ Formulation} \\\\\n", " (Kn^{(3)}):& \\partial_{t}^{(2)} f_{i}^{(0)} +\\partial_{t}^{(1)} f_{i}^{(1)} + \\partial_{t}^{(1)}\\left( \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(0)} \\right) + \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(2)} + \\displaystyle\\frac{\\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)^{2} f_{i}^{(1)}}{2} + \\displaystyle\\frac{\\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)^{3} f_{i}^{(0)}}{6} = - \\displaystyle\\frac{ f_{i}^{(3)} }{ \\tau }, \\\\\n", " \\textrm{ou}\\\\\n", " (Kn^{(3)}):& \\partial_{t}^{(2)} f_{i}^{(0)} + \\left(1-\\displaystyle\\frac{1}{\\tau}\\right)\\partial_{t}^{(1)} f_{i}^{(1)} + \\displaystyle\\left(\\frac{1}{2} - \\frac{1}{6\\tau}\\right) \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)^{2} f_{i}^{(1)} + \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(2)} = - \\displaystyle\\frac{f_{i}^{(3)}}{\\tau},\n", " \\end{array}\n", "$$(Chap-Kn-cv-NotEq)\n", "\n", "##### Zero-Order Moment Balance\n", "\n", "To retrieve the balance equation, we sum the Eq. {eq}`Chap-Kn-cv-Eq` for $Kn^{(1)}$ and $Kn^{(2)}$ over $\\sum_{i=0} $:\n", "\n", "$$\n", "\\begin{array}{l}\n", "(Kn^{(1)}): \\displaystyle\\sum_{i}\\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(0)}= - \\displaystyle\\sum_{i}\\left(\\displaystyle\\frac{ f_{i}^{(1)} }{ \\tau } \\right),\\\\\n", "(Kn^{(1)}): \\partial_{t}^{(0)} u + \\partial_{\\alpha} B_{\\alpha} =0,\n", "\\end{array}\n", "$$(cv-Kn1-B0-all)\n", "\n", "$$\n", "\\begin{array}{l}\n", "(Kn^{(2)}):& \\displaystyle\\sum_{i}\\left[ \\partial_{t}^{(1)} f_{i}^{(0)} + \\displaystyle\\left(1 - \\frac{1}{2\\tau}\\right)\\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right) f_{i}^{(1)} \\right]= - \\displaystyle\\sum_{i}\\left(\\displaystyle\\frac{ f_{i}^{(2)} }{ \\tau } \\right), \\nonumber \\\\\n", "(Kn^{(2)}):& \\partial_{t}^{(1)} u + \\displaystyle\\left(1 - \\frac{1}{2\\tau}\\right)\\partial_{\\alpha} m^{(1)}_{\\alpha} = 0,\n", "\\end{array}\n", "$$(cv-Kn2-B0-all)\n", "\n", "To determine the unknow term $m^{(1)}_{\\alpha}$, we sum the Eq. {eq}`Chap-Kn-cv-Eq` for $Kn^{(1)}$ over the first-order moment:\n", "\n", "$$\n", "\\begin{array}{l}\n", "(Kn^{(1)}): \\displaystyle\\sum_{i}\\left( \\partial_{t}^{(0)} + c_{i,\\beta}\\partial_{\\beta} \\right)f_{i}^{(0)}c_{i,\\alpha}= - \\displaystyle\\sum_{i}\\left(\\displaystyle\\frac{ f_{i}^{(1)} }{ \\tau } c_{i,\\alpha} \\right),\\\\\n", "(Kn^{(1)}): \\partial_{t}^{(0)} B_{\\alpha} + \\partial_{\\beta} \\left(uc_{s}^{2}\\right)\\delta_{\\alpha\\beta} = - \\displaystyle\\frac{m^{(1)}_{\\alpha}}{\\tau} \\quad \\quad \\rightarrow \\quad \\quad m^{(1)}_{\\alpha} =-\\tau\\left( \\partial_{t}^{(0)} B_{\\alpha} + \\partial_{\\beta} \\left(uc_{s}^{2}\\right)\\delta_{\\alpha\\beta} \\right).\n", "\\end{array}\n", "$$(cv-Kn1-B1-all)\n", "\n", "\n", "Extending the analysis up to $Kn^{(3)}$\n", "\n", "$$\n", "\\begin{array}{l}\n", "(Kn^{(3)}):& \\displaystyle\\sum_{i}\\left[ \\partial_{t}^{(2)} f_{i}^{(0)} + \\left(1-\\displaystyle\\frac{1}{\\tau}\\right)\\partial_{t}^{(1)} f_{i}^{(1)} + \\displaystyle\\left(\\frac{1}{2} - \\frac{1}{6\\tau}\\right) \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)^{2} f_{i}^{(1)} + \\left( \\partial_{t}^{(0)} + c_{i,\\alpha}\\partial_{\\alpha} \\right)f_{i}^{(2)} \\right] = - \\displaystyle\\sum_{i}\\left[\\displaystyle\\frac{f_{i}^{(3)}}{\\tau} \\right], \\\\\n", "(Kn^{(3)}):& \\partial_{t}^{(2)} u + \\partial_{t}^{(1)}\\left(\\partial_{\\alpha} m^{(1)}_{\\alpha} + \\partial_{\\beta} m^{(1)}_{\\beta}\\right) + \\partial_{\\alpha}\\partial_{\\beta} m^{(1)}_{\\alpha\\beta} + \\partial_{\\alpha} m^{(2)}_{\\alpha} = 0.\n", "\\end{array}\n", "$$(cv-Kn3-B0-all)\n", "\n", "By summing the Eqs. {eq}`cv-Kn1-B0`, {eq}`cv-Kn2-B0`, and {eq}`cv-Kn3-B0`, and also extending to higher orders of Knudsen, we have\n", "\n", "$$\n", "\\partial_{t} u + \\partial_{\\alpha} B_{\\alpha} - \\partial_{\\alpha}\\left(\\nu \\partial_{\\alpha} u \\right) + \\underbrace{\\partial_{\\alpha} m^{(2)}_{\\alpha} + \\partial_{\\alpha}\\partial_{\\beta} m^{(1)}_{\\alpha\\beta} + \\partial_{t}^{(1)}\\left(\\partial_{\\alpha} m^{(1)}_{\\alpha} + \\partial_{\\beta} m^{(1)}_{\\beta}\\right) - \\left(\\tau - \\frac{1}{2} \\right) \\partial_{t}^{(0)} \\partial_{\\alpha} B_{\\alpha} + ... }_{\\textrm{Numerical Errors}} =0,\n", "$$\n", "\n", "where $\\nu=c_{s}^{2}(\\tau-1/2)$. The numerical errors can be corrected using extra moments in the equilibrium distribution function and some ones can neglected due to its relevance decay exponentially with the decrease of $\\Delta x$ and $\\Delta t$.\n", "```" ] }, { "cell_type": "markdown", "id": "12d0efed-4938-41ee-90d6-e6971f1f1146", "metadata": {}, "source": [ "#### Boundary Conditions\n", "\n", "The boundary conditions for the lattices can be derived by solving a linear system of known moments." ] }, { "cell_type": "markdown", "id": "ed35c46d-0d61-4be2-9c27-68f579754ab9", "metadata": {}, "source": [ "| D1Q3: | Boundaries | | Layers | | |\n", "|------|-----------------|---|---------------------------------------------------------------------------------|---|---------------------------------------------------------------------------------|\n", "| | | | West | | East |\n", "| | Unknown $f_{i}$ | | $f_2=u_{e} - f_1 - f_0$ | | $f_1=u_{w} - f_2 - f_0$ |\n", "| | | | | | |" ] }, { "cell_type": "markdown", "id": "f18b8ebe-8315-40a1-99cb-f1a44883b1a9", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## Benchmark - 1D Non-Linear Convective formulation of the Buckley-Leverett Equation" ] }, { "cell_type": "markdown", "id": "a8cd98fa-1ac0-412c-a58d-dbbe2b4c1dcb", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "Considering a pure convective formulation of the Buckley-Leverett equation, where gravitational effects and capillary pressure are neglected, we have:\n", "\n", "$$\n", "\\partial_{t} (s_{w})+ \\partial_{\\alpha}\\left[\\displaystyle\\frac{\\left_{\\alpha}}{\\varphi} F(s_w)\\right] = 0 ,\n", "$$(EqMassMTP-BL-Sim)\n", "\n", "with the water fraction flow given by,\n", "\n", "$$\n", "F(s_w) = \\frac{ s_{w}^{2} }{ s_{w}^{2} + M_{\\mu}(1-s_{w})^{2}}.\n", "$$(EqFractionalFlow)\n", "\n", "The simulation parameters are chosen so that the average velocity component $\\langle u \\rangle_{\\alpha} / \\varphi$ equals $1\\,\\mathrm{m/s}$ and the dynamic viscosity ratio $M_{\\mu}$ equals $1$. We consider a one-dimensional porous medium of length $L_x = 20\\,\\mathrm{m}$, initially saturated with oil (see {numref}`1D-Displace-BL`). Water is injected at the inlet $(x=0)$ under a Dirichlet boundary condition $s_{w}(x=0)=1$. At the outlet $(x=L)$, a Neumann boundary condition specifying a zero gradient in water saturation is applied.\n", "\n", "```{figure} 1D-Displace-BL.svg\n", "---\n", "scale: 100%\n", "align: center\n", "name: 1D-Displace-BL\n", "---\n", "Geometry of the non-linear convective problem: initial and boundary conditions.\n", "```" ] }, { "cell_type": "markdown", "id": "f825d322-636b-402a-a91d-cb365e412b3e", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "### Transient Analytical Solution" ] }, { "cell_type": "markdown", "id": "a9a8e2b2-89ba-48a5-9962-43b8f72a0499", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "The non-linear convective problem given by the Eq. {eq}`EqMassMTP-BL-Sim` has an analytical solution dependent of the initial condition (Riemann problem). In the present case, the solution for the specific initial condition is given by\n", "\n", "$$\n", "s_{w}(x,t)= \\left\\{ \\begin{array}{ccc}\n", " t \\displaystyle F' & \\textrm{for} & x\\leq \\zeta t \\\\\n", " & & \\\\\n", " 0 & \\textrm{for} & x > \\zeta t \\\\ \n", "\\end{array} \\right.\n", "$$\n", "\n", "where $\\zeta=F'(\\tilde{s}_{w})=F(\\tilde{s}_{w})/\\tilde{s}_{w}$ is the shock wave speed, and $\\tilde{s}_{w}$ is the sharp front saturation determined analytically by\n", "\n", "$$\n", "\\tilde{s}_{w}=\\displaystyle\\frac{F(\\tilde{s}_{w})}{F'(\\tilde{s}_{w})}=\\sqrt{\\displaystyle\\frac{M_{\\mu}}{M_{\\mu}+1}}, \\quad \\quad \\textrm{where} \\quad \\quad F^{'}(s_{w})=\\frac{dF}{ds_{w}}=\\frac{2s_{w}M_{\\mu}(1-s_{w})}{(s_{w}^{2}+M_{\\mu}(1-s_{w})^{2})^{2}}.\n", "$$" ] }, { "cell_type": "markdown", "id": "9af7c174-f34d-496e-9482-2e288407b1a2", "metadata": {}, "source": [ "#### Forward Time and Forward Space - FTFS" ] }, { "cell_type": "markdown", "id": "2ce1bdfa-141e-48a6-b42a-a0e1c61a478c", "metadata": {}, "source": [ "#### FDM Solution:" ] }, { "cell_type": "markdown", "id": "17abf656-77be-48b0-997d-100c2877fc5d", "metadata": {}, "source": [ "```{toggle}\n", "\n", "```python\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from pylab import *\n", "matplotlib.rcParams['mathtext.fontset'] = 'cm'\n", "#----------------------------Simulation-Set--------------------------------------------------------\n", "Nx=125;\n", "L = 20.0 # Length of the reservoir\n", "T = 4.0 # Total simulation time\n", "dx = L / (Nx) # Spatial step size\n", "c=2**(2) # c=dx/dt\n", "dt=dx/c # Time step size\n", "nt = int(T / dt) # Time step number\n", "uo=1.0 # Constant Fluid Flux\n", "m=1.0 # Dynamic Viscosity ration\n", "print('dx=',dx,'\\t dt=',dt)\n", "#----------------------------Initilizing-Simulation------------------------------------------------\n", "sf=np.zeros((Nx),dtype=\"float64\") # Allocating Density Field\n", "sfo=np.zeros((Nx),dtype=\"float64\") # Allocating Density Buffer Field\n", "sf[:]=0.0 # Initilal Density Condition\n", "Bx=np.zeros((Nx),dtype=\"float64\") # Allocating Convective Term\n", "Bx=uo*sf**(2.0)/(sf**(2.0)+m*(1.0-sf)**(2.0)) # Initializing Convective Term\n", "for t in range(nt):\n", "# for t in range(3):\n", " sfo[1:] = sf[1:] - dt/dx * (Bx[1:] - np.roll(Bx, 1, axis=0)[1:])\n", " sfo[0] = 1.0 # Inlet boundary condition\n", " sfo[Nx-1] = 0.0 # Inlet boundary condition\n", " sf=sfo\n", " Bx=uo*sf**(2.0)/(sf**(2.0)+m*(1.0-sf)**(2.0))\n", "x = np.linspace(0, L, Nx)\n", "plt.plot(x, sf, 'k-' ,label='Water Saturation (t=4)')\n", "plt.xlabel(r\"$x$\",fontsize=18)\n", "plt.ylabel(r\"$u(x,t)$\",fontsize=18)\n", "plt.xlim(0,20)\n", "plt.ylim(0,1)\n", "plt.legend()\n", "plt.grid(True)\n", "plt.tight_layout()\n", "```" ] }, { "cell_type": "markdown", "id": "2186a0dd-9713-4390-8788-8cffd0896783", "metadata": {}, "source": [ "#### Lattice Boltzmann Framework of FTFS:" ] }, { "cell_type": "markdown", "id": "81c0d99a-c3d2-4962-bb15-39c355be317a", "metadata": {}, "source": [ "```{toggle}\n", "\n", "```python\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from pylab import *\n", "matplotlib.rcParams['mathtext.fontset'] = 'cm'\n", "#***************************************InputParameters************************************************\n", "Nx=125 #Square Domain Length\n", "#***************************************Lattice-Properties-D2Q5*************************************************\n", "cx = np.array([0, 1],dtype=\"int8\") \n", "L = 20.0 # Length of the reservoir\n", "T = 4.0 # Total simulation time\n", "dx = L / (Nx) # Spatial step size\n", "c=2**(2) # c=dx/dt\n", "dt=dx/c # Time step size\n", "nt = int(T / dt) # Time step number\n", "uo=1.0 # Constant Fluid Flux\n", "m=1.0 # Dynamic Viscosity ration\n", "print('dx=',dx,'\\t dt=',dt)\n", "#***************************************LBM-Scale************************************************\n", "ue=uo/c\n", "sk=np.zeros((Nx),dtype=\"float64\")\n", "Bx=np.zeros((Nx),dtype=\"float64\")\n", "f=np.zeros((2,Nx),dtype=\"float64\")\n", "f[0,:]=sk[:]-Bx[:]\n", "f[1,:]=Bx[:]\n", "for t in range(nt):\n", "# for t in range(3):\n", " #--------------------Collision----------------\n", " f[0,:]=sk[:]-Bx[:]\n", " f[1,:]=Bx[:]\n", " #-----------------streaming-------------------\n", " for k in range(0,2):\n", " f[k,:]=np.roll(f[k,:], cx[k], axis=0)\n", " #-----------------Boundaries-----------------------\n", " f[1,0]= 1.0 - f[0,0]\n", " f[0,Nx-1]=f[0,Nx-2]\n", " f[1,Nx-1]=f[1,Nx-2]\n", " #----------------------Macro------------------\n", " sk[:]=f[0,:]+f[1,:]\n", " Bx=ue*sk*sk/(sk*sk+(1.0-sk)**(2.0))\n", "x = np.linspace(0, L, Nx)\n", "plt.plot(x, sk, 'k-' ,label='Water Saturation $(t=4)$')\n", "plt.xlabel(r\"$x$\",fontsize=18)\n", "plt.ylabel(r\"$u(x,t)$\",fontsize=18)\n", "plt.xlim(0,20)\n", "plt.ylim(0,1)\n", "plt.legend()\n", "plt.grid(True)\n", "plt.tight_layout()\n", "```" ] }, { "cell_type": "markdown", "id": "9bea51fb-c5fa-4ee9-8a7d-bee8d61f279e", "metadata": {}, "source": [ "#### Numerical and Analytical Comparisson" ] }, { "cell_type": "markdown", "id": "71f0589c-0b60-4c83-b2b2-6ddc6e6351dd", "metadata": {}, "source": [ "```{toggle}\n", "\n", "```python\n", "plt.plot(x, sf, 'ks:' ,label='FDM (FTFS): Water Saturation $(t=4)$',fillstyle=\"none\")\n", "plt.plot(x, sk, 'ro:' ,label='LBM (FTFS): Water Saturation $(t=4)$',fillstyle=\"none\")\n", "plt.xlabel(r\"$x$\",fontsize=18)\n", "plt.ylabel(r\"$u(x,t)$\",fontsize=18)\n", "plt.xlim(0,20)\n", "plt.ylim(0,1)\n", "plt.legend()\n", "plt.grid(True)\n", "plt.tight_layout()\n", "```" ] }, { "cell_type": "markdown", "id": "bae77fa5-0db8-40b6-b1c5-f189c0efbd4c", "metadata": {}, "source": [ "```{figure} sat-FTFS.png\n", "---\n", "scale: 100%\n", "align: center\n", "name: sat-FTFS\n", "---\n", "Comparisson of saturation profile for FDM and LBM FTFS schemes.\n", "```" ] }, { "cell_type": "markdown", "id": "6931990d-e70c-44ef-8d95-f2fea281e4b3", "metadata": {}, "source": [ "### Forward-Time and Central-Space with Small Diffusivity (FTCS-SD)" ] }, { "cell_type": "markdown", "id": "20f2387b-68d3-4831-beae-ddb5625709ee", "metadata": {}, "source": [ "#### Lattice Boltzmann Method (D1Q3):" ] }, { "cell_type": "markdown", "id": "3a8f2dab-7d7c-4033-aac8-636e891b44e3", "metadata": {}, "source": [ "```{toggle}\n", "\n", "```python\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from pylab import *\n", "matplotlib.rcParams['mathtext.fontset'] = 'cm'\n", "#***************************************InputParameters************************************************\n", "Nx=126 #Square Domain Length\n", "#***************************************Lattice-Properties-D1Q3*************************************************\n", "cs=1.0/sqrt(3.0);\n", "w = np.array([4.0/6.0, 1.0/6.0, 1.0/6.0],dtype=\"float64\")\n", "cx = np.array([0, 1, -1],dtype=\"int16\") \n", "#***************************************LBM-Scale************************************************\n", "L = 20.0 # Length of the reservoir\n", "T = 10.0 # Total simulation time\n", "dx = L / (Nx-1) # Spatial step size\n", "c=2**(2) # c=dx/dt\n", "dt=dx/c # Time step size\n", "nt = int(T / dt) # Time step number\n", "uo=1.0 # Constant Fluid Flux\n", "m=1.0 # Dynamic Viscosity ration\n", "print('dx=',dx,'\\t dt=',dt)\n", "#***************************************Initialization************************************************\n", "ue=uo/c\n", "sk=np.zeros((Nx),dtype=\"float64\")\n", "Bx=np.zeros((Nx),dtype=\"float64\")\n", "f=np.zeros((3,Nx),dtype=\"float64\")\n", "for k in range(0,3):\n", " f[k,:]=w[k]*(sk[:]+cx[k]*Bx[:]*3.0)\n", "#***************************************Loop************************************************\n", "for t in range(nt):\n", " #----------------------Macro------------------\n", " sk[:]=f[0,:]+f[1,:]+f[2,:]\n", " Bx=ue*sk*sk/(sk*sk+(1.0-sk)**(2.0))\n", " #--------------------Collision----------------\n", " for k in range(0,3):\n", " f[k,:]=w[k]*(sk[:]+cx[k]*Bx[:]/cs**2)\n", " #-----------------streaming-------------------\n", " for k in range(0,3):\n", " f[k,:]=np.roll(f[k,:], cx[k], axis=0)\n", " #-----------------Boundaries-----------------------\n", " f[1,0]= 1.0 - f[0,0]-f[2,0]\n", " f[0,Nx-1]=f[0,Nx-2]\n", " f[1,Nx-1]=f[1,Nx-2]\n", " f[2,Nx-1]=f[2,Nx-2]\n", "x = np.linspace(0, L, Nx)\n", "plt.plot(x, sk, 'k-' ,label='Water Saturation $(t=10)$')\n", "plt.xlabel(r\"$x$\",fontsize=18)\n", "plt.ylabel(r\"$u(x,t)$\",fontsize=18)\n", "plt.xlim(0,20)\n", "plt.ylim(0,1)\n", "plt.legend()\n", "plt.grid(True)\n", "plt.tight_layout()\n", "```" ] }, { "cell_type": "markdown", "id": "6416111f-886b-451a-b190-0e885d8fd8f1", "metadata": {}, "source": [ "#### FDM approach of LBM Framework D1Q3:" ] }, { "cell_type": "markdown", "id": "1b091026-bd7e-4131-894e-07c110281426", "metadata": {}, "source": [ "```{toggle}\n", "\n", "```python\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from pylab import *\n", "matplotlib.rcParams['mathtext.fontset'] = 'cm'\n", "#----------------------------Simulation-Set--------------------------------------------------------\n", "Nx=126;\n", "L = 20.0 # Length of the reservoir\n", "T = 10.0 # Total simulation time\n", "dx = L / (Nx-1) # Spatial step size\n", "c=2**(2) # c=dx/dt\n", "dt=dx/c # Time step size\n", "nt = int(T / dt) # Time step number\n", "uo=1.0 # Constant Fluid Flux\n", "m=1.0 # Dynamic Viscosity ration\n", "print('dx=',dx,'\\t dt=',dt)\n", "#----------------------------Initilizing-Simulation------------------------------------------------\n", "sf=np.zeros((Nx),dtype=\"float64\") # Allocating Density Field\n", "sfo=np.zeros((Nx),dtype=\"float64\") # Allocating Density Buffer Field\n", "sf[:]=0.0 # Initilal Density Condition\n", "Bx=np.zeros((Nx),dtype=\"float64\") # Allocating Convective Term\n", "Bx=uo*sf**(2.0)/(sf**(2.0)+m*(1.0-sf)**(2.0)) # Initializing Convective Term\n", "for t in range(nt):\n", "# for t in range(3):\n", " sf=sfo\n", " Bx=uo*sf**(2.0)/(sf**(2.0)+m*(1.0-sf)**(2.0))\n", " sfo = ( (4.0*sf/6.0 + np.roll(sf, cx[1], axis=0)/6.0 + np.roll(sf, cx[2], axis=0)/6.0) \n", " + dt/dx * (np.roll(Bx, cx[1], axis=0) - np.roll(Bx, cx[2], axis=0))/2.0 )\n", " sfo[0] = 1.0 # Inlet boundary condition\n", " sfo[Nx-1] = 0.0 # Inlet boundary condition\n", "x = np.linspace(0, L, Nx)\n", "plt.plot(x, sk, 'k-' ,label='Water Saturation $(t=10)$')\n", "plt.xlabel(r\"$x$\",fontsize=18)\n", "plt.ylabel(r\"$u(x,t)$\",fontsize=18)\n", "plt.xlim(0,20)\n", "plt.ylim(0,1)\n", "plt.legend()\n", "plt.grid(True)\n", "plt.tight_layout()\n", "```" ] }, { "cell_type": "markdown", "id": "e68c1e35-4c7e-4e7c-8b13-7052e880e162", "metadata": {}, "source": [ "#### Numerical and Analytical Comparisson" ] }, { "cell_type": "markdown", "id": "6964c333-9cc8-4844-a0d8-5e2dbb954f70", "metadata": {}, "source": [ "```{toggle}\n", "\n", "```python\n", "plt.plot(x, sf, 'ks:' ,label='FDM (FTCS-SD): Water Saturation $(t=10)$',markersize=6,fillstyle=\"none\")\n", "plt.plot(x, sk, 'ro:' ,label='LBM D1Q3: Water Saturation $(t=10)$',markersize=6,fillstyle=\"none\")\n", "plt.xlabel(r\"$x$\",fontsize=18)\n", "plt.ylabel(r\"$u(x,t)$\",fontsize=18)\n", "plt.xlim(0,20)\n", "plt.ylim(0,1)\n", "plt.legend()\n", "plt.grid(True)\n", "plt.tight_layout()\n", "```" ] }, { "cell_type": "markdown", "id": "511dca43-5d8e-4d42-83dd-565ca2de6783", "metadata": {}, "source": [ "```{figure} sat-FTCS-SD.png\n", "---\n", "scale: 100%\n", "align: center\n", "name: sat-FTCS-SD\n", "---\n", "Comparisson of saturation profile for FDM and LBM FTCS-SD schemes.\n", "```" ] } ], "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.10.12" } }, "nbformat": 4, "nbformat_minor": 5 }