From c59e452d7c55dda46cb5bed7f5d93a2638719ce7 Mon Sep 17 00:00:00 2001 From: ratnania Date: Fri, 23 Feb 2024 12:57:11 +0100 Subject: [PATCH] adding applications and sphinx documentation of psydac and sympde --- .gitignore | 2 +- _config.yml | 4 + _toc.yml | 34 ++- chapter2/elliptic-general-form.ipynb | 278 ++++++++++++++++++ chapter5/bingham_plastic_flow_in_a_pipe.ipynb | 78 +++++ .../buoyancy-driven_natural_convection.ipynb | 81 +++++ chapter5/casson_fluid_flow_in_a_channel.ipynb | 74 +++++ chapter5/cfd.md | 1 + chapter5/compressible_flow_in_a_nozzle.ipynb | 80 +++++ ...tion-diffusion_equation_in_a_channel.ipynb | 78 +++++ ...ss_power_law_fluid_flow_in_a_channel.ipynb | 77 +++++ chapter5/electromagnetics.md | 56 ++++ chapter5/free_surface_flow.ipynb | 80 +++++ chapter5/fsi.md | 75 +++++ chapter5/heat_conduction_in_a_solid.ipynb | 74 +++++ ...erschel-bulkley_fluid_flow_in_a_pipe.ipynb | 79 +++++ .../incompressible_flow_past_a_cylinder.ipynb | 75 +++++ chapter5/magnetohydrodynamics_flow.ipynb | 81 +++++ chapter5/material-science.md | 1 + chapter5/mhd.md | 55 ++++ chapter5/multiphysics.md | 1 + chapter5/non-newtonian-fluids.md | 1 + .../oldroyd-b_fluid_flow_in_a_channel.ipynb | 79 +++++ ...apanastasiou_fluid_flow_in_a_channel.ipynb | 76 +++++ chapter5/particle-laden_flow.ipynb | 79 +++++ .../power_law_fluid_flow_in_a_channel.ipynb | 76 +++++ .../stokes_flow_in_a_lid-driven_cavity.ipynb | 78 +++++ chapter5/two-phase_flow.ipynb | 78 +++++ documentation/psydac.rst | 19 ++ documentation/sympde.rst | 14 + intro.md | 10 + logo.png | Bin 7812 -> 10614 bytes 32 files changed, 1872 insertions(+), 2 deletions(-) create mode 100644 chapter2/elliptic-general-form.ipynb create mode 100644 chapter5/bingham_plastic_flow_in_a_pipe.ipynb create mode 100644 chapter5/buoyancy-driven_natural_convection.ipynb create mode 100644 chapter5/casson_fluid_flow_in_a_channel.ipynb create mode 100644 chapter5/cfd.md create mode 100644 chapter5/compressible_flow_in_a_nozzle.ipynb create mode 100644 chapter5/convection-diffusion_equation_in_a_channel.ipynb create mode 100644 chapter5/cross_power_law_fluid_flow_in_a_channel.ipynb create mode 100644 chapter5/electromagnetics.md create mode 100644 chapter5/free_surface_flow.ipynb create mode 100644 chapter5/fsi.md create mode 100644 chapter5/heat_conduction_in_a_solid.ipynb create mode 100644 chapter5/herschel-bulkley_fluid_flow_in_a_pipe.ipynb create mode 100644 chapter5/incompressible_flow_past_a_cylinder.ipynb create mode 100644 chapter5/magnetohydrodynamics_flow.ipynb create mode 100644 chapter5/material-science.md create mode 100644 chapter5/mhd.md create mode 100644 chapter5/multiphysics.md create mode 100644 chapter5/non-newtonian-fluids.md create mode 100644 chapter5/oldroyd-b_fluid_flow_in_a_channel.ipynb create mode 100644 chapter5/papanastasiou_fluid_flow_in_a_channel.ipynb create mode 100644 chapter5/particle-laden_flow.ipynb create mode 100644 chapter5/power_law_fluid_flow_in_a_channel.ipynb create mode 100644 chapter5/stokes_flow_in_a_lid-driven_cavity.ipynb create mode 100644 chapter5/two-phase_flow.ipynb create mode 100644 documentation/psydac.rst create mode 100644 documentation/sympde.rst diff --git a/.gitignore b/.gitignore index 6700fd7..17ee1da 100644 --- a/.gitignore +++ b/.gitignore @@ -18,4 +18,4 @@ usr .env .iga-python - +*_autosummary* diff --git a/_config.yml b/_config.yml index 9de863f..eaf1248 100644 --- a/_config.yml +++ b/_config.yml @@ -38,5 +38,9 @@ html: sphinx: extra_extensions: - sphinx_proof + - sphinx.ext.autodoc + - sphinx.ext.autosummary + config: + autosummary_generate: True exclude_patterns: [README.md, .iga-python] diff --git a/_toc.yml b/_toc.yml index c0a3846..2c0b69e 100644 --- a/_toc.yml +++ b/_toc.yml @@ -37,6 +37,7 @@ parts: chapters: - file: chapter2/poisson - file: chapter2/poisson-nitsche + - file: chapter2/elliptic-general-form - file: chapter2/biharmonic - file: chapter2/vector-poisson - file: chapter2/advection-diffusion @@ -65,4 +66,35 @@ parts: - file: chapter4/subdomains sections: - file: chapter4/poisson-nitsche - + - caption: Applications + chapters: + - file: chapter5/cfd + sections: + - file: chapter5/buoyancy-driven_natural_convection + - file: chapter5/compressible_flow_in_a_nozzle + - file: chapter5/convection-diffusion_equation_in_a_channel + - file: chapter5/free_surface_flow + - file: chapter5/heat_conduction_in_a_solid + - file: chapter5/incompressible_flow_past_a_cylinder + - file: chapter5/magnetohydrodynamics_flow + - file: chapter5/particle-laden_flow + - file: chapter5/stokes_flow_in_a_lid-driven_cavity + - file: chapter5/two-phase_flow + - file: chapter5/non-newtonian-fluids + sections: + - file: chapter5/power_law_fluid_flow_in_a_channel + - file: chapter5/bingham_plastic_flow_in_a_pipe + - file: chapter5/oldroyd-b_fluid_flow_in_a_channel + - file: chapter5/herschel-bulkley_fluid_flow_in_a_pipe + - file: chapter5/cross_power_law_fluid_flow_in_a_channel + - file: chapter5/casson_fluid_flow_in_a_channel + - file: chapter5/papanastasiou_fluid_flow_in_a_channel + - file: chapter5/fsi + - file: chapter5/electromagnetics + - file: chapter5/mhd + - file: chapter5/material-science + - file: chapter5/multiphysics + - caption: Documentation + chapters: + - file: documentation/sympde + - file: documentation/psydac diff --git a/chapter2/elliptic-general-form.ipynb b/chapter2/elliptic-general-form.ipynb new file mode 100644 index 0000000..34a98c7 --- /dev/null +++ b/chapter2/elliptic-general-form.ipynb @@ -0,0 +1,278 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b5e0474c", + "metadata": {}, + "source": [ + "# Elliptic equation in the general form\n", + "\n", + "We consider here, the following general form of an elliptic partial differential equation,\n", + "\n", + "$$\n", + "\\begin{align}\n", + "- \\nabla \\cdot \\left( A \\nabla u \\right) + \\mathbf{b} \\cdot \\nabla u + c u &= f, \\quad \\Omega \\\\\n", + "u &= 0, \\quad \\partial \\Omega\n", + "\\end{align}\n", + "$$\n", + "\n", + "## Variational Formulation\n", + "\n", + "An $H^1$-conforming variational formulation of reads\n", + "\n", + "$$\n", + "\\begin{align}\n", + " \\text{find $u \\in V$ such that} \\quad a(u,v) = l(v) \\quad \\forall v \\in V,\n", + "\\end{align}\n", + "$$\n", + "\n", + "where \n", + "\n", + "- $V \\subset H^1(\\Omega)$, \n", + "- $a(u,v) := \\int_{\\Omega} \\left( \\left( A \\nabla u \\right) \\cdot \\nabla v + \\left( \\mathbf{b} \\cdot \\nabla u \\right) v + cuv \\right) ~ d\\Omega$,\n", + "- $l(v) := \\int_{\\Omega} f v ~ d\\Omega$.\n" + ] + }, + { + "cell_type": "markdown", + "id": "c5c944a1", + "metadata": {}, + "source": [ + "## Formal Model" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "c2cee2d9", + "metadata": {}, + "outputs": [], + "source": [ + "from sympde.expr import BilinearForm, LinearForm, integral\n", + "from sympde.expr import find, EssentialBC, Norm, SemiNorm\n", + "from sympde.topology import ScalarFunctionSpace, Square, element_of\n", + "from sympde.calculus import grad, dot, div\n", + "from sympde.core import Vector, Matrix\n", + "\n", + "from sympy import pi, sin\n", + "\n", + "from psydac.api.discretization import discretize\n", + "\n", + "domain = Square()\n", + "\n", + "V = ScalarFunctionSpace('V', domain)\n", + "\n", + "x,y = domain.coordinates\n", + "\n", + "u,v = [element_of(V, name=i) for i in ['u', 'v']]\n", + "\n", + "c = x*y\n", + "b = Vector([1e-2, 1e-1], name='b')\n", + "A = Matrix([[1,1], [0,1]], name='A')\n", + "\n", + "# bilinear form\n", + "expr = dot(grad(v), A * grad(u)) + dot(b, grad(u))*v + c*u*v\n", + "a = BilinearForm((u,v), integral(domain, expr))" + ] + }, + { + "cell_type": "markdown", + "id": "366c2d67", + "metadata": {}, + "source": [ + "### Manifactured solution" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "2b05c79f", + "metadata": {}, + "outputs": [], + "source": [ + "# the analytical solution and its rhs\n", + "ue = sin(pi * x) * sin(pi * y)\n", + "\n", + "L = lambda u: - div(A*grad(u)) + dot(b,grad(u)) + c*u\n", + "f = L(ue)" + ] + }, + { + "cell_type": "markdown", + "id": "66b50430", + "metadata": {}, + "source": [ + "### Formal Equation" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "d0d386a2", + "metadata": {}, + "outputs": [], + "source": [ + "l = LinearForm(v, integral(domain, f*v))\n", + "\n", + "# Dirichlet boundary conditions\n", + "bc = [EssentialBC(u, 0, domain.boundary)]\n", + "\n", + "# Variational problem\n", + "equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc)" + ] + }, + { + "cell_type": "markdown", + "id": "0614c740", + "metadata": {}, + "source": [ + "## Discretization" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "96c29a27", + "metadata": {}, + "outputs": [], + "source": [ + "degree = [2,2]\n", + "ncells = [8,8]" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "ef6a9709", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/ranania/PYCCEL/IGA-Python/.iga-python/lib/python3.10/site-packages/sympy/matrices/repmatrix.py:90: SymPyDeprecationWarning: \n", + "\n", + "non-Expr objects in a Matrix has been deprecated since SymPy 1.9. Use\n", + "list of lists, TableForm or some other data structure instead. See\n", + "https://github.com/sympy/sympy/issues/21497 for more info.\n", + "\n", + " ).warn()\n" + ] + } + ], + "source": [ + "# Create computational domain from topological domain\n", + "domain_h = discretize(domain, ncells=ncells, comm=None)\n", + "\n", + "# Create discrete spline space\n", + "Vh = discretize(V, domain_h, degree=degree)\n", + "\n", + "# Discretize equation\n", + "equation_h = discretize(equation, domain_h, [Vh, Vh])" + ] + }, + { + "cell_type": "markdown", + "id": "c286e572", + "metadata": {}, + "source": [ + "### Solving the PDE" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "3bafe9f5", + "metadata": {}, + "outputs": [], + "source": [ + "equation_h.set_solver('gmres', info=False, tol=1e-8)" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "b88daf50", + "metadata": {}, + "outputs": [], + "source": [ + "uh = equation_h.solve()" + ] + }, + { + "cell_type": "markdown", + "id": "d865a17f", + "metadata": {}, + "source": [ + "## Computing the error norm" + ] + }, + { + "cell_type": "markdown", + "id": "be88e26c", + "metadata": {}, + "source": [ + "### Computing the $L^2$ norm" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "3440e74a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.00021948770944141364\n" + ] + } + ], + "source": [ + "u = element_of(V, name='u')\n", + "\n", + "# create the formal Norm object\n", + "l2norm = Norm(u - ue, domain, kind='l2')\n", + "\n", + "# discretize the norm\n", + "l2norm_h = discretize(l2norm, domain_h, Vh)\n", + "\n", + "# assemble the norm\n", + "l2_error = l2norm_h.assemble(u=uh)\n", + "\n", + "# print the result\n", + "print(l2_error)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5d7d9180", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/bingham_plastic_flow_in_a_pipe.ipynb b/chapter5/bingham_plastic_flow_in_a_pipe.ipynb new file mode 100644 index 0000000..3d1717f --- /dev/null +++ b/chapter5/bingham_plastic_flow_in_a_pipe.ipynb @@ -0,0 +1,78 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "1f0b7959", + "metadata": {}, + "source": [ + "# Bingham Plastic Flow in a Pipe\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the Bingham plastic model for non-Newtonian fluids with yield stress:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\frac{\\partial u}{\\partial t} + u\\frac{\\partial u}{\\partial x} &= -\\frac{1}{\\rho}\\frac{\\partial p}{\\partial x} + \\frac{\\tau_0}{\\rho}\\frac{\\partial^2 u}{\\partial y^2} \\\\\n", + "\\frac{\\partial^2 u}{\\partial y^2} &\\geq 0 \\quad \\text{(Yield stress condition)}\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $u$ is the velocity, $p$ is the pressure, $\\rho$ is the density, $\\tau_0$ is the yield stress, and $y$ is the radial coordinate.\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Inlet: Specify velocity or pressure conditions.\n", + "- Outlet: Specify outflow conditions.\n", + "- Pipe walls: No-slip conditions.\n", + "\n", + "### Initial Conditions\n", + "\n", + "Specify initial conditions for velocity and pressure.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $u \\in H_0^1(\\Omega)$ and $p \\in L^2(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\frac{\\partial u}{\\partial t}v \\, d\\Omega + \\int_{\\Omega} u\\frac{\\partial u}{\\partial x}v \\, d\\Omega + \\int_{\\Omega} \\frac{1}{\\rho}\\frac{\\partial p}{\\partial x}v \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\frac{\\tau_0}{\\rho}\\frac{\\partial^2 u}{\\partial y^2}\\frac{\\partial v}{\\partial y} \\, d\\Omega \\\\\n", + "&= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $v \\in H_0^1(\\Omega)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3104c066", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/buoyancy-driven_natural_convection.ipynb b/chapter5/buoyancy-driven_natural_convection.ipynb new file mode 100644 index 0000000..767c6a4 --- /dev/null +++ b/chapter5/buoyancy-driven_natural_convection.ipynb @@ -0,0 +1,81 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "315d2f4f", + "metadata": {}, + "source": [ + "# Buoyancy-Driven Natural Convection\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the Boussinesq approximation for buoyancy-driven natural convection:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\frac{\\partial \\mathbf{v}}{\\partial t} + (\\mathbf{v} \\cdot \\nabla)\\mathbf{v} &= -\\frac{1}{\\rho_0}\\nabla p + \\nu \\nabla^2 \\mathbf{v} + \\beta g (T - T_0)\\mathbf{k} \\\\\n", + "\\frac{\\partial T}{\\partial t} + \\mathbf{v} \\cdot \\nabla T &= \\alpha \\nabla^2 T\n", + "\\end{align}\n", + "$$\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Vertical walls: $\\mathbf{v} \\cdot \\mathbf{t} = 0$ and $\\frac{\\partial T}{\\partial n} = 0$.\n", + "\n", + "- Inlet and outlet: Specify velocity or temperature conditions.\n", + "\n", + "- Top and bottom: No-slip and adiabatic conditions.\n", + "\n", + "### Initial Conditions\n", + "\n", + "Start with an initial velocity and temperature distribution.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $\\mathbf{v} \\in H_0^1(\\Omega)$ and $T \\in H^1(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\left(\\frac{\\partial \\mathbf{v}}{\\partial t} + (\\mathbf{v} \\cdot \\nabla)\\mathbf{v}\\right) \\cdot \\mathbf{w} \\, d\\Omega + \\int_{\\Omega} \\frac{1}{\\rho_0}\\nabla p \\cdot \\mathbf{w} \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\nu \\nabla \\mathbf{v} : \\nabla \\mathbf{w} \\, d\\Omega - \\beta g (T - T_0)\\mathbf{k} \\cdot \\mathbf{w} \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\left(\\frac{\\partial T}{\\partial t} + \\mathbf{v} \\cdot \\nabla T - \\alpha \\nabla^2 T\\right) \\psi \\, d\\Omega \\\\\n", + "&= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $\\mathbf{w}, \\psi \\in H_0^1(\\Omega)$.\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc53cb64", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/casson_fluid_flow_in_a_channel.ipynb b/chapter5/casson_fluid_flow_in_a_channel.ipynb new file mode 100644 index 0000000..9c7b3b3 --- /dev/null +++ b/chapter5/casson_fluid_flow_in_a_channel.ipynb @@ -0,0 +1,74 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e9f1cb0e", + "metadata": {}, + "source": [ + "# Casson Fluid Flow in a Channel\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the Casson model for non-Newtonian fluids:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\sigma_{ij} &= -p\\delta_{ij} + \\sqrt{\\frac{\\tau_0}{\\rho}}\\sqrt{\\frac{\\tau_0}{\\rho} + \\frac{1}{2}\\rho\\left(\\frac{\\partial u_i}{\\partial x_j} + \\frac{\\partial u_j}{\\partial x_i}\\right)^2} \\\\\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $\\sigma_{ij}$ is the stress tensor, $p$ is the pressure, $\\rho$ is the density, $u_i$ is the velocity, and $\\tau_0$ is the yield stress.\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Inlet and outlet: Specify velocity or pressure conditions.\n", + "- Channel walls: No-slip conditions.\n", + "\n", + "### Initial Conditions\n", + "Start with an initial guess for velocity and pressure.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $u_i \\in H_0^1(\\Omega)$ and $p \\in L^2(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\sigma_{ij} \\left(\\frac{\\partial u_i}{\\partial x_j} + \\frac{\\partial u_j}{\\partial x_i}\\right) \\, d\\Omega + \\int_{\\Omega} \\frac{\\partial p}{\\partial x_i}u_i \\, d\\Omega \\\\\n", + "&= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $u_i \\in H_0^1(\\Omega)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "52879066", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/cfd.md b/chapter5/cfd.md new file mode 100644 index 0000000..c156b86 --- /dev/null +++ b/chapter5/cfd.md @@ -0,0 +1 @@ +# Computational Fluid Dynamics diff --git a/chapter5/compressible_flow_in_a_nozzle.ipynb b/chapter5/compressible_flow_in_a_nozzle.ipynb new file mode 100644 index 0000000..d5b4e52 --- /dev/null +++ b/chapter5/compressible_flow_in_a_nozzle.ipynb @@ -0,0 +1,80 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "896cf22a", + "metadata": {}, + "source": [ + "# Compressible Flow in a Nozzle\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the compressible Euler equations for flow in a converging-diverging nozzle:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\frac{\\partial \\rho}{\\partial t} + \\nabla \\cdot (\\rho \\mathbf{v}) &= 0 \\\\\n", + "\\frac{\\partial (\\rho \\mathbf{v})}{\\partial t} + \\nabla \\cdot (\\rho \\mathbf{v} \\otimes \\mathbf{v} + p\\mathbf{I}) &= 0 \\\\\n", + "\\frac{\\partial E}{\\partial t} + \\nabla \\cdot \\left[(E + p) \\mathbf{v}\\right] &= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Inlet: Specify inflow conditions ($\\rho$, $\\mathbf{v}$, $p$).\n", + "\n", + "- Outlet: Specify outflow conditions or use characteristic-based boundary conditions.\n", + "\n", + "- Nozzle walls: $\\mathbf{v} \\cdot \\mathbf{n} = 0$ (no penetration).\n", + "\n", + "### Initial Conditions\n", + "\n", + "Specify initial conditions for density, velocity, and pressure.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $\\rho, \\mathbf{v}, p \\in H^1(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\left(\\frac{\\partial \\rho}{\\partial t} + \\nabla \\cdot (\\rho \\mathbf{v})\\right) \\phi \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\left(\\frac{\\partial (\\rho \\mathbf{v})}{\\partial t} + \\nabla \\cdot (\\rho \\mathbf{v} \\otimes \\mathbf{v} + p\\mathbf{I})\\right) \\cdot \\mathbf{w} \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\left(\\frac{\\partial E}{\\partial t} + \\nabla \\cdot \\left[(E + p) \\mathbf{v}\\right]\\right) \\psi \\, d\\Omega \\\\\n", + "&= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $\\phi, \\mathbf{w}, \\psi \\in H^1(\\Omega)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8a3bde01", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/convection-diffusion_equation_in_a_channel.ipynb b/chapter5/convection-diffusion_equation_in_a_channel.ipynb new file mode 100644 index 0000000..d283e61 --- /dev/null +++ b/chapter5/convection-diffusion_equation_in_a_channel.ipynb @@ -0,0 +1,78 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "6889a63e", + "metadata": {}, + "source": [ + "# Convection-Diffusion Equation in a Channel\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the 2D convection-diffusion equation in a channel:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\rho u \\frac{\\partial c}{\\partial x} + \\rho v \\frac{\\partial c}{\\partial y} - \\nabla \\cdot (D \\nabla c) = f\n", + "\\end{align}\n", + "$$\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Inlet: $c = c_{\\text{inlet}}$ with \\$\\mathbf{v} \\cdot \\mathbf{n} < 0$ (inward flow).\n", + "\n", + "- Outlet: $\\frac{\\partial c}{\\partial x} = \\frac{\\partial c}{\\partial y} = 0$ (zero gradient).\n", + "\n", + "- Lateral faces: $-D \\nabla c \\cdot \\mathbf{n} = 0$ (insulated).\n", + "\n", + "### Initial Conditions\n", + "\n", + "Assume an initial concentration distribution $c_0(x, y)$ at $t = 0$.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $c \\in H_0^1(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\left(\\rho u \\frac{\\partial c}{\\partial x} + \\rho v \\frac{\\partial c}{\\partial y}\\right) \\phi \\, d\\Omega - \\int_{\\Omega} D \\nabla c \\cdot \\nabla \\phi \\, d\\Omega \\\\\n", + "&= \\int_{\\Omega} f \\phi \\, d\\Omega\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $\\phi \\in H_0^1(\\Omega)$.\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d3c6c22", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/cross_power_law_fluid_flow_in_a_channel.ipynb b/chapter5/cross_power_law_fluid_flow_in_a_channel.ipynb new file mode 100644 index 0000000..13fb81b --- /dev/null +++ b/chapter5/cross_power_law_fluid_flow_in_a_channel.ipynb @@ -0,0 +1,77 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "45606393", + "metadata": {}, + "source": [ + "# Cross Power Law Fluid Flow in a Channel\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the cross power-law model for non-Newtonian fluids:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\sigma_{ij} &= -p\\delta_{ij} + \\lambda \\varepsilon_{kk}\\delta_{ij} + 2\\mu\\left(\\varepsilon_{ij} + \\alpha\\varepsilon_{ij}^2\\right) \\\\\n", + "\\varepsilon_{ij} &= \\frac{1}{2}\\left(\\frac{\\partial v_i}{\\partial x_j} + \\frac{\\partial v_j}{\\partial x_i}\\right) \\\\\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $\\sigma_{ij}$ is the stress tensor, $p$ is the pressure, $\\lambda$ and $\\mu$ are material constants, $v_i$ is the velocity, $\\varepsilon_{ij}$ is the strain rate, and $\\alpha$ is a power-law index.\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Inlet and outlet: Specify velocity or pressure conditions.\n", + "\n", + "- Channel walls: No-slip conditions.\n", + "\n", + "### Initial Conditions\n", + "\n", + "Start with an initial guess for velocity and pressure.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $v_i \\in H_0^1(\\Omega)$ and $p \\in L^2(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\sigma_{ij} \\varepsilon_{ij} \\, d\\Omega + \\int_{\\Omega} \\frac{\\partial p}{\\partial x_i}v_i \\, d\\Omega \\\\\n", + "&= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $v_i \\in H_0^1(\\Omega)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "85946710", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/electromagnetics.md b/chapter5/electromagnetics.md new file mode 100644 index 0000000..9f793c7 --- /dev/null +++ b/chapter5/electromagnetics.md @@ -0,0 +1,56 @@ +# Electromagnetics + +Electromagnetic problems are commonly described by Maxwell's equations, which govern the behavior of electric and magnetic fields. The Finite Element Method (FEM) provides a powerful numerical approach for solving these equations in complex geometries. This section provides a concise overview of the mathematical formulation for electromagnetic problems using finite elements. + +## Maxwell's Equations + +The time-harmonic Maxwell's equations in free space are given by: + +$$ +\begin{align} + \nabla \times \mathbf{E} &= -j\omega \mu \mathbf{H}, \label{eq:maxwell1} \\ + \nabla \times \mathbf{H} &= j\omega \varepsilon \mathbf{E}, \label{eq:maxwell2} \\ + \nabla \cdot \mathbf{B} &= 0, \label{eq:maxwell3} \\ + \nabla \cdot \mathbf{D} &= 0, \label{eq:maxwell4} +\end{align} +$$ + +where $\mathbf{E}$ is the electric field, $\mathbf{H}$ is the magnetic field, $\mathbf{B}$ is the magnetic flux density, $\mathbf{D}$ is the electric displacement field, $\omega$ is the angular frequency, $\mu$ is the permeability, and $\varepsilon$ is the permittivity. + +## Weak Formulation + +The weak form of Maxwell's equations is obtained by multiplying each equation with suitable test functions and integrating over the domain $\Omega$. For the electric field $\mathbf{E}$, the weak form is given by: + +$$ +\begin{equation} + \int_{\Omega} \nabla \times \mathbf{E} \cdot \nabla \times \mathbf{v} \, d\Omega + \omega^2 \mu \varepsilon \mathbf{E} \cdot \mathbf{v} \, d\Omega = 0, +\end{equation} +$$ + +where $\mathbf{v}$ is a test function belonging to the function space $H(\text{curl}; \Omega)$. + +Similarly, for the magnetic field $\mathbf{H}$, the weak form is given by: + +$$ +\begin{equation} + \int_{\Omega} \nabla \times \mathbf{H} \cdot \nabla \times \mathbf{u} \, d\Omega - \omega^2 \mu \varepsilon \mathbf{H} \cdot \mathbf{u} \, d\Omega = 0, +\end{equation} +$$ + +where $\mathbf{u}$ is a test function belonging to the function space $H(\text{curl}; \Omega)$. + +## Finite Element Discretization + +To apply the finite element method, the domain $\Omega$ is discretized into elements, and the electromagnetic fields are approximated using piecewise basis functions. The discretized weak form leads to a system of linear equations, which can be solved numerically to obtain the finite element solution. + +## References + +The following references provide in-depth coverage of the mathematical background for electromagnetism using finite elements: + +\begin{enumerate} + \item \cite{jin2015finite} + \item \cite{bossavit1998whitney} + \item \cite{monk2003finite} +\end{enumerate} + +This mathematical background serves as the foundation for implementing finite element simulations of electromagnetic problems. Researchers interested in detailed mathematical derivations and computational techniques are encouraged to explore the referenced works. diff --git a/chapter5/free_surface_flow.ipynb b/chapter5/free_surface_flow.ipynb new file mode 100644 index 0000000..d482cab --- /dev/null +++ b/chapter5/free_surface_flow.ipynb @@ -0,0 +1,80 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "fea980ee", + "metadata": {}, + "source": [ + "# Free Surface Flow (Navier-Stokes with Free Surface)\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the Navier-Stokes equations coupled with a level set or volume-of-fluid method for tracking a free surface:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\rho \\frac{\\partial \\mathbf{v}}{\\partial t} + \\rho (\\mathbf{v} \\cdot \\nabla) \\mathbf{v} &= -\\nabla p + \\mu \\nabla^2 \\mathbf{v} + \\rho \\mathbf{g} \\\\\n", + "\\nabla \\cdot \\mathbf{v} &= 0 \\\\\n", + "\\frac{\\partial \\phi}{\\partial t} + \\mathbf{v} \\cdot \\nabla \\phi &= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $\\phi$ is the level set or volume fraction representing the free surface.\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Inlet and outlet: Specify velocity or pressure conditions.\n", + "- No-slip conditions on solid boundaries.\n", + "- Implement the level set or volume fraction boundary conditions.\n", + "\n", + "### Initial Conditions\n", + "\n", + "Specify initial conditions for velocity, pressure, and the level set or volume fraction.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $\\mathbf{v} \\in H_0^1(\\Omega)$, $p \\in L^2(\\Omega)$, and $\\phi \\in H^1(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\rho \\frac{\\partial \\mathbf{v}}{\\partial t} \\cdot \\mathbf{w} \\, d\\Omega + \\int_{\\Omega} \\rho (\\mathbf{v} \\cdot \\nabla) \\mathbf{v} \\cdot \\mathbf{w} \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\mu \\nabla \\mathbf{v} : \\nabla \\mathbf{w} \\, d\\Omega - \\int_{\\Omega} p \\nabla \\cdot \\mathbf{w} \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\nabla \\cdot \\mathbf{v} q \\, d\\Omega + \\int_{\\Omega} \\frac{\\partial \\phi}{\\partial t} q \\, d\\Omega + \\int_{\\Omega} \\mathbf{v} \\cdot \\nabla \\phi q \\, d\\Omega \\\\\n", + "&= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $\\mathbf{w} \\in H_0^1(\\Omega)$, $q \\in L^2(\\Omega)$, and $r \\in H^1(\\Omega)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b3ac645a", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/fsi.md b/chapter5/fsi.md new file mode 100644 index 0000000..fb4f1c7 --- /dev/null +++ b/chapter5/fsi.md @@ -0,0 +1,75 @@ +# Fluid-Structure Interaction + +Fluid-Structure Interaction (FSI) involves the coupled interaction between a fluid and a structure, where the motion of one influences the behavior of the other. The Finite Element Method (FEM) is a powerful tool for simulating FSI problems. This section provides an overview of the mathematical formulation for fluid-structure interaction using finite elements. + +## Coupled Fluid-Structure Equations + +The governing equations for fluid-structure interaction involve the Navier-Stokes equations for the fluid and the equations of motion for the structure. In a partitioned approach, the coupled system is given by: + +\begin{align} + \text{Fluid Domain:} \quad \rho_f \left(\frac{\partial \mathbf{u}_f}{\partial t} + (\mathbf{u}_f \cdot \nabla)\mathbf{u}_f\right) &= -\nabla p_f + \mu_f \nabla^2 \mathbf{u}_f + \mathbf{f}_f + \mathbf{F}_s, \label{eq:fsi_fluid_momentum} \\ + \nabla \cdot \mathbf{u}_f &= 0, \label{eq:fsi_fluid_continuity} \\ + \text{Structure Domain:} \quad \rho_s \frac{\partial^2 \mathbf{u}_s}{\partial t^2} &= \nabla \cdot \mathbf{P}_s + \mathbf{f}_s, \label{eq:fsi_structure_motion} +\end{align} + +where: + +- $\mathbf{u}_f$ is the fluid velocity, +- $p_f$ is the fluid pressure, +- $\rho_f$ is the fluid density, +- $\mu_f$ is the fluid dynamic viscosity, +- $\mathbf{f}_f$ is the fluid body force, +- $\mathbf{F}_s$ is the fluid-structure interaction force on the fluid by the structure, +- $\mathbf{u}_s$ is the structure displacement, +- $\rho_s$ is the structure density, +- $\mathbf{P}_s$ is the structure internal force, +- $\mathbf{f}_s$ is the structure external force. + +## Coupling Conditions + +For a fluid-structure interaction problem, coupling conditions at the fluid-structure interface are essential. These conditions ensure continuity of velocities, pressures, and forces between the fluid and structure. A common approach is to enforce kinematic and dynamic conditions: + +$$ +\begin{align} + \text{Kinematic Condition:} \quad \mathbf{u}_f &= \mathbf{u}_s, \label{eq:fsi_kinematic_condition} \\ + \text{Dynamic Condition:} \quad \mathbf{P}_f \cdot \mathbf{n} &= \mathbf{P}_s \cdot \mathbf{n}, \label{eq:fsi_dynamic_condition} +\end{align} +$$ + +where $\mathbf{n}$ is the unit outward normal vector at the fluid-structure interface, and $\mathbf{P}_f$ and $\mathbf{P}_s$ are the fluid and structure stresses, respectively. + +## Weak Formulation + +The weak form of the coupled fluid-structure problem involves the variational formulation of the fluid and structure equations, incorporating the coupling conditions. For the fluid domain, the weak form is given by: + +$$ +\begin{equation} + \int_{\Omega_f} \rho_f \left(\frac{\partial \mathbf{u}_f}{\partial t} + (\mathbf{u}_f \cdot \nabla)\mathbf{u}_f\right) \cdot \mathbf{v}_f \, d\Omega_f + \int_{\Omega_f} \mu_f \nabla \mathbf{u}_f : \nabla \mathbf{v}_f \, d\Omega_f = \int_{\Omega_f} (\mathbf{f}_f + \mathbf{F}_s) \cdot \mathbf{v}_f \, d\Omega_f, +\end{equation} +$$ + +where $\mathbf{v}_f$ is a test function belonging to the function space $H^1(\Omega_f)$, and $\Omega_f$ is the fluid domain. + +For the structure domain, the weak form is given by: + +$$ +\begin{equation} + \int_{\Omega_s} \rho_s \frac{\partial^2 \mathbf{u}_s}{\partial t^2} \cdot \mathbf{v}_s \, d\Omega_s + \int_{\Omega_s} \nabla \cdot \mathbf{P}_s \cdot \mathbf{v}_s \, d\Omega_s = \int_{\Omega_s} \mathbf{f}_s \cdot \mathbf{v}_s \, d\Omega_s, +\end{equation} +$$ + +where $\mathbf{v}_s$ is a test function belonging to the function space $H^1(\Omega_s)$, and $\Omega_s$ is the structure domain. + +The coupling conditions, such as \eqref{eq:fsi_kinematic_condition} and \eqref{eq:fsi_dynamic_condition}, should be enforced during the assembly of the global matrices. + +## References + +The following references provide comprehensive coverage of the mathematical background for fluid-structure interaction using finite elements: + +\begin{enumerate} + \item \cite{hron2006fluid} + \item \cite{bathe2014finite} + \item \cite{quarteroni2017fluid} +\end{enumerate} + +Researchers interested in detailed mathematical derivations and computational techniques are encouraged to explore the referenced works. diff --git a/chapter5/heat_conduction_in_a_solid.ipynb b/chapter5/heat_conduction_in_a_solid.ipynb new file mode 100644 index 0000000..dfe455d --- /dev/null +++ b/chapter5/heat_conduction_in_a_solid.ipynb @@ -0,0 +1,74 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b3466488", + "metadata": {}, + "source": [ + "# Heat Conduction in a Solid\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the heat conduction equation in a solid domain \\(\\Omega\\):\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\rho C_p \\frac{\\partial T}{\\partial t} - \\nabla \\cdot (k \\nabla T) = Q\n", + "\\end{align}\n", + "$$\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Dirichlet: $T = T_{\\text{inlet}}$ on the inlet face.\n", + "- Neumann: $-k \\frac{\\partial T}{\\partial n} = 0$ on the outlet face.\n", + "- Insulated: $-k \\frac{\\partial T}{\\partial n} = 0$ on the lateral faces.\n", + "\n", + "### Initial Conditions\n", + "\n", + "Assume an initial temperature distribution $T_0(x, y, z)$ at $t = 0$.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $T \\in L^2(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\int_{\\Omega} \\rho C_p \\frac{\\partial T}{\\partial t} \\phi \\, d\\Omega + \\int_{\\Omega} k \\nabla T \\cdot \\nabla \\phi \\, d\\Omega = \\int_{\\Omega} Q \\phi \\, d\\Omega\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $\\phi \\in H_0^1(\\Omega)$.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "55ba709f", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/herschel-bulkley_fluid_flow_in_a_pipe.ipynb b/chapter5/herschel-bulkley_fluid_flow_in_a_pipe.ipynb new file mode 100644 index 0000000..0313108 --- /dev/null +++ b/chapter5/herschel-bulkley_fluid_flow_in_a_pipe.ipynb @@ -0,0 +1,79 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a0c1f909", + "metadata": {}, + "source": [ + "# Herschel-Bulkley Fluid Flow in a Pipe\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the Herschel-Bulkley model for non-Newtonian fluids with yield stress:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\frac{\\partial u}{\\partial t} + u\\frac{\\partial u}{\\partial x} &= -\\frac{1}{\\rho}\\frac{\\partial p}{\\partial x} + \\frac{\\tau_0}{\\rho}\\left(\\frac{\\partial u}{\\partial x}\\right)^{n-1} \\\\\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $u$ is the velocity, $p$ is the pressure, $\\rho$ is the density, $\\tau_0$ is the yield stress, and $n$ is the flow behavior index.\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Inlet: Specify velocity or pressure conditions.\n", + "\n", + "- Outlet: Specify outflow conditions.\n", + "\n", + "- Pipe walls: No-slip conditions.\n", + "\n", + "### Initial Conditions\n", + "\n", + "Specify initial conditions for velocity and pressure.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $u \\in H_0^1(\\Omega)$ and $p \\in L^2(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\frac{\\partial u}{\\partial t}v \\, d\\Omega + \\int_{\\Omega} u\\frac{\\partial u}{\\partial x}v \\, d\\Omega + \\int_{\\Omega} \\frac{1}{\\rho}\\frac{\\partial p}{\\partial x}v \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\frac{\\tau_0}{\\rho}\\left(\\frac{\\partial u}{\\partial x}\\right)^{n-1}v \\, d\\Omega \\\\\n", + "&= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $v \\in H_0^1(\\Omega)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1f098463", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/incompressible_flow_past_a_cylinder.ipynb b/chapter5/incompressible_flow_past_a_cylinder.ipynb new file mode 100644 index 0000000..71188b8 --- /dev/null +++ b/chapter5/incompressible_flow_past_a_cylinder.ipynb @@ -0,0 +1,75 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0b948ce3", + "metadata": {}, + "source": [ + "# Incompressible Flow Past a Cylinder\n", + "\n", + "## Mathematical Model:\n", + "\n", + "Consider the steady, incompressible Navier-Stokes equations for fluid flow around a cylinder:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\rho (\\mathbf{v} \\cdot \\nabla) \\mathbf{v} &= -\\nabla p + \\mu \\nabla^2 \\mathbf{v} + \\rho \\mathbf{g} \\\\\n", + "\\nabla \\cdot \\mathbf{v} &= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Inlet: $\\mathbf{v} = (U, 0)$ (uniform velocity).\n", + "- Outlet: $\\frac{\\partial \\mathbf{v}}{\\partial x} = 0$ (zero gradient).\n", + "- Cylinder surface: $\\mathbf{v} \\cdot \\mathbf{n} = 0$ (no-slip condition).\n", + "\n", + "### Initial Conditions\n", + "\n", + "Assume a quiescent fluid, i.e., $\\mathbf{v} = (0, 0)$ at $t = 0$.\n", + "\n", + "### Weak Formulation\n", + "\n", + "Find $\\mathbf{v} \\in H_0^1(\\Omega)$ and $p \\in L^2(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\int_{\\Omega} \\rho (\\mathbf{v} \\cdot \\nabla) \\mathbf{v} \\cdot \\mathbf{w} \\, d\\Omega & = -\\int_{\\Omega} p \\nabla \\cdot \\mathbf{w} \\, d\\Omega + \\mu \\int_{\\Omega} \\nabla \\mathbf{v} : \\nabla \\mathbf{w} \\, d\\Omega + \\int_{\\Omega} \\rho \\mathbf{g} \\cdot \\mathbf{w} \\, d\\Omega \\\\\n", + "\\int_{\\Omega} \\nabla \\cdot \\mathbf{v} q \\, d\\Omega & = 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $\\mathbf{w} \\in H_0^1(\\Omega)$ and $q \\in L^2(\\Omega)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ed14fcc5", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/magnetohydrodynamics_flow.ipynb b/chapter5/magnetohydrodynamics_flow.ipynb new file mode 100644 index 0000000..6b0f078 --- /dev/null +++ b/chapter5/magnetohydrodynamics_flow.ipynb @@ -0,0 +1,81 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "6c15c8b4", + "metadata": {}, + "source": [ + "# Magnetohydrodynamics (MHD) Flow\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the MHD equations coupling the Navier-Stokes equations with Maxwell's equations:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\rho \\frac{\\partial \\mathbf{v}}{\\partial t} + \\rho (\\mathbf{v} \\cdot \\nabla) \\mathbf{v} &= -\\nabla p + \\mu \\nabla^2 \\mathbf{v} + \\mathbf{J} \\times \\mathbf{B} \\\\\n", + "\\frac{\\partial \\mathbf{B}}{\\partial t} &= \\nabla \\times (\\mathbf{v} \\times \\mathbf{B} - \\eta \\nabla \\times \\mathbf{B}) \\\\\n", + "\\nabla \\cdot \\mathbf{v} &= 0 \\\\\n", + "\\nabla \\cdot \\mathbf{B} &= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $\\mathbf{J} = \\nabla \\times \\mathbf{B}$ is the current density.\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Inlet and outlet: Specify velocity or pressure conditions.\n", + "- Magnetic field conditions depending on the problem.\n", + "\n", + "### Initial Conditions\n", + "\n", + "Specify initial conditions for velocity, pressure, and the magnetic field.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $\\mathbf{v} \\in H_0^1(\\Omega)$, $p \\in L^2(\\Omega)$, and $\\mathbf{B} \\in H(\\text{curl}; \\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\rho \\frac{\\partial \\mathbf{v}}{\\partial t} \\cdot \\mathbf{w} \\, d\\Omega + \\int_{\\Omega} \\rho (\\mathbf{v} \\cdot \\nabla) \\mathbf{v} \\cdot \\mathbf{w} \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\mu \\nabla \\mathbf{v} : \\nabla \\mathbf{w} \\, d\\Omega - \\int_{\\Omega} p \\nabla \\cdot \\mathbf{w} \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\frac{\\partial \\mathbf{B}}{\\partial t} \\cdot \\boldsymbol{\\psi} \\, d\\Omega - \\int_{\\Omega} \\nabla \\times (\\mathbf{v} \\times \\mathbf{B} - \\eta \\nabla \\times \\mathbf{B}) \\cdot \\boldsymbol{\\psi} \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\nabla \\cdot \\mathbf{v} q \\, d\\Omega + \\int_{\\Omega} \\nabla \\cdot \\mathbf{B} r \\, d\\Omega \\\\\n", + "&= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $\\mathbf{w} \\in H_0^1(\\Omega)$, $q \\in L^2(\\Omega)$, $\\boldsymbol{\\psi} \\in H(\\text{curl}; \\Omega)$, and $r \\in L^2(\\Omega)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "850e45f1", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/material-science.md b/chapter5/material-science.md new file mode 100644 index 0000000..1d2e5f8 --- /dev/null +++ b/chapter5/material-science.md @@ -0,0 +1 @@ +# Material Science diff --git a/chapter5/mhd.md b/chapter5/mhd.md new file mode 100644 index 0000000..cf5f17c --- /dev/null +++ b/chapter5/mhd.md @@ -0,0 +1,55 @@ +# MHD + +In magnetohydrodynamics (MHD), finite element analysis can be applied to solve a variety of problems related to the behavior of electrically conducting fluids (plasmas or liquid metals) in the presence of magnetic fields. Here are some common problems in magnetohydrodynamics that can be addressed using finite element methods. + +Finite element analysis provides a powerful and flexible framework for addressing these complex problems in magnetohydrodynamics by discretizing the governing equations and solving them numerically. The specific application will determine the details of the problem setup and the required numerical techniques. + +## MHD Equations + +Magnetohydrodynamics (MHD) is a branch of physics and fluid dynamics that studies the magnetic properties and behavior of electrically conducting fluids, such as plasmas, liquid metals, and saltwater. The MHD equations describe the coupled interactions between the magnetic field, fluid flow, and electric current. Here are the mathematical models for the MHD equations: + +1. MHD Continuity Equation: + $$ + \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0 + $$ + where: + - $\rho$ is the fluid density, + - $\mathbf{v}$ is the fluid velocity. + +2. MHD Momentum Equation: + $$ + \rho \frac{D\mathbf{v}}{Dt} = -\nabla p + \rho \mathbf{g} + \nabla \cdot \boldsymbol{\tau} + \mathbf{J} \times \mathbf{B} + $$ + where: + - $\frac{D\mathbf{v}}{Dt}$ is the material derivative of velocity, + - $p$ is the pressure, + - $\mathbf{g}$ is the gravitational acceleration, + - $\boldsymbol{\tau}$ is the stress tensor, + - $\mathbf{J}$ is the current density, + - $\mathbf{B}$ is the magnetic field. + +3. MHD Induction Equation: + $$ + \frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{v} \times \mathbf{B} - \eta \nabla \times \mathbf{B}) + $$ + where: + - $\eta$ is the magnetic diffusivity. + +4. MHD Ohm's Law: + $$ + \mathbf{J} = \sigma (\mathbf{E} + \mathbf{v} \times \mathbf{B}) + $$ + where: + - $\sigma$ is the electrical conductivity, + - $\mathbf{E}$ is the electric field. + +5. MHD Energy Equation: + $$ + \rho C_p \frac{DT}{Dt} = -p \nabla \cdot \mathbf{v} + \nabla \cdot (k \nabla T) + \frac{1}{\sigma}(\mathbf{J} \cdot \mathbf{E}) + $$ + where: + - $C_p$ is the specific heat at constant pressure, + - $T$ is the temperature, + - $k$ is the thermal conductivity. + +These equations, when solved together with appropriate boundary conditions, describe the behavior of magnetized fluids in the presence of electric and magnetic fields. Keep in mind that the specific form of these equations may vary based on the assumptions made and the type of fluid being considered (e.g., ideal MHD vs. resistive MHD). diff --git a/chapter5/multiphysics.md b/chapter5/multiphysics.md new file mode 100644 index 0000000..3748ec9 --- /dev/null +++ b/chapter5/multiphysics.md @@ -0,0 +1 @@ +# Multiphysics diff --git a/chapter5/non-newtonian-fluids.md b/chapter5/non-newtonian-fluids.md new file mode 100644 index 0000000..9ae636c --- /dev/null +++ b/chapter5/non-newtonian-fluids.md @@ -0,0 +1 @@ +# Non-Newtonian Fluids diff --git a/chapter5/oldroyd-b_fluid_flow_in_a_channel.ipynb b/chapter5/oldroyd-b_fluid_flow_in_a_channel.ipynb new file mode 100644 index 0000000..b7a8440 --- /dev/null +++ b/chapter5/oldroyd-b_fluid_flow_in_a_channel.ipynb @@ -0,0 +1,79 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "aea3782b", + "metadata": {}, + "source": [ + "# Oldroyd-B Fluid Flow in a Channel\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the Oldroyd-B model for viscoelastic fluids:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\frac{\\partial \\mathbf{u}}{\\partial t} + \\mathbf{u}\\cdot\\nabla\\mathbf{u} &= -\\frac{1}{\\rho}\\nabla p + \\nabla\\cdot\\boldsymbol{\\tau} \\\\\n", + "\\frac{\\partial \\boldsymbol{\\tau}}{\\partial t} + \\mathbf{u}\\cdot\\nabla\\boldsymbol{\\tau} &= \\boldsymbol{\\tau}\\cdot\\nabla\\mathbf{u} - \\frac{\\boldsymbol{\\tau}}{\\lambda} + 2\\mu\\left(\\nabla\\mathbf{u} - \\frac{\\boldsymbol{\\tau}}{\\lambda}\\right) \\\\\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $\\mathbf{u}$ is the velocity, $p$ is the pressure, $\\rho$ is the density, $\\boldsymbol{\\tau}$ is the extra stress tensor, $\\lambda$ is the relaxation time, and $\\mu$ is the viscosity.\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Inlet and outlet: Specify velocity or pressure conditions.\n", + "- Channel walls: No-slip conditions.\n", + "\n", + "### Initial Conditions\n", + "\n", + "Specify initial conditions for velocity and pressure.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $\\mathbf{u} \\in H_0^1(\\Omega)$ and $p \\in L^2(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\frac{\\partial \\mathbf{u}}{\\partial t}\\cdot\\mathbf{w} \\, d\\Omega + \\int_{\\Omega} \\mathbf{u}\\cdot\\nabla\\mathbf{u}\\cdot\\mathbf{w} \\, d\\Omega + \\int_{\\Omega} \\frac{1}{\\rho}\\nabla p\\cdot\\mathbf{w} \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\frac{\\partial \\boldsymbol{\\tau}}{\\partial t}:\\mathbf{v} \\, d\\Omega + \\int_{\\Omega} \\mathbf{u}\\cdot\\nabla\\boldsymbol{\\tau}:\\mathbf{v} \\, d\\Omega \\\\\n", + "&- \\int_{\\Omega} \\boldsymbol{\\tau}\\cdot\\nabla\\mathbf{u}:\\mathbf{v} \\, d\\Omega + \\int_{\\Omega} \\frac{\\boldsymbol{\\tau}}{\\lambda}:\\mathbf{v} \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} 2\\mu\\left(\\nabla\\mathbf{u} - \\frac{\\boldsymbol{\\tau}}{\\lambda}\\right):\\nabla\\mathbf{w} \\, d\\Omega \\\\\n", + "&= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $\\mathbf{w} \\in H_0^1(\\Omega)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b949a16a", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/papanastasiou_fluid_flow_in_a_channel.ipynb b/chapter5/papanastasiou_fluid_flow_in_a_channel.ipynb new file mode 100644 index 0000000..c0a4bd3 --- /dev/null +++ b/chapter5/papanastasiou_fluid_flow_in_a_channel.ipynb @@ -0,0 +1,76 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "d8b8b5bd", + "metadata": {}, + "source": [ + "# Papanastasiou Fluid Flow in a Channel\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the Papanastasiou model for non-Newtonian fluids:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\sigma_{ij} &= -p\\delta_{ij} + 2\\mu\\left(\\varepsilon_{ij} + \\frac{\\tau_0}{2\\mu}\\right)\\left(\\frac{\\partial u_i}{\\partial x_j} + \\frac{\\partial u_j}{\\partial x_i}\\right) \\\\\n", + "\\varepsilon_{ij} &= \\frac{1}{2}\\left(\\frac{\\partial u_i}{\\partial x_j} + \\frac{\\partial u_j}{\\partial x_i}\\right) \\\\\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $\\sigma_{ij}$ is the stress tensor, $p$ is the pressure, $\\mu$ is the viscosity, $u_i$ is the velocity, $\\varepsilon_{ij}$ is the strain rate, and $\\tau_0$ is the yield stress.\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Inlet and outlet: Specify velocity or pressure conditions.\n", + "- Channel walls: No-slip conditions.\n", + "\n", + "### Initial Conditions\n", + "\n", + "Start with an initial guess for velocity and pressure.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $u_i \\in H_0^1(\\Omega)$ and $p \\in L^2(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\sigma_{ij} \\left(\\frac{\\partial u_i}{\\partial x_j} + \\frac{\\partial u_j}{\\partial x_i}\\right) \\, d\\Omega + \\int_{\\Omega} \\frac{\\partial p}{\\partial x_i}u_i \\, d\\Omega \\\\\n", + "&= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $u_i \\in H_0^1(\\Omega)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "02de07e4", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/particle-laden_flow.ipynb b/chapter5/particle-laden_flow.ipynb new file mode 100644 index 0000000..eaf2574 --- /dev/null +++ b/chapter5/particle-laden_flow.ipynb @@ -0,0 +1,79 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b4a03bc4", + "metadata": {}, + "source": [ + "# Particle-Laden Flow (Lagrangian-Eulerian Approach)\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the incompressible Navier-Stokes equations coupled with the Lagrangian motion of particles:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\rho \\frac{\\partial \\mathbf{v}}{\\partial t} + \\rho (\\mathbf{v} \\cdot \\nabla) \\mathbf{v} &= -\\nabla p + \\mu \\nabla^2 \\mathbf{v} + \\mathbf{F}_{\\text{drag}} \\\\\n", + "\\nabla \\cdot \\mathbf{v} &= 0 \\\\\n", + "\\frac{d \\mathbf{x}}{dt} &= \\mathbf{v}\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $\\mathbf{F}_{\\text{drag}}$ represents the drag force on the particles.\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Inlet and outlet: Specify velocity or pressure conditions.\n", + "- Solid particle boundaries: Implement no-slip conditions.\n", + "\n", + "### Initial Conditions\n", + "\n", + "Specify initial conditions for velocity, pressure, and particle positions.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $\\mathbf{v} \\in H_0^1(\\Omega)$ and $p \\in L^2(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\rho \\frac{\\partial \\mathbf{v}}{\\partial t} \\cdot \\mathbf{w} \\, d\\Omega + \\int_{\\Omega} \\rho (\\mathbf{v} \\cdot \\nabla) \\mathbf{v} \\cdot \\mathbf{w} \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\mu \\nabla \\mathbf{v} : \\nabla \\mathbf{w} \\, d\\Omega - \\int_{\\Omega} p \\nabla \\cdot \\mathbf{w} \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\nabla \\cdot \\mathbf{v} q \\, d\\Omega \\\\\n", + "&= -\\int_{\\Omega} \\mathbf{F}_{\\text{drag}} \\cdot \\mathbf{w} \\, d\\Omega\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $\\mathbf{w} \\in H_0^1(\\Omega)$ and $q \\in L^2(\\Omega)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe99e28f", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/power_law_fluid_flow_in_a_channel.ipynb b/chapter5/power_law_fluid_flow_in_a_channel.ipynb new file mode 100644 index 0000000..5170b3e --- /dev/null +++ b/chapter5/power_law_fluid_flow_in_a_channel.ipynb @@ -0,0 +1,76 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "33aa7908", + "metadata": {}, + "source": [ + "# Power Law Fluid Flow in a Channel\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the steady-state power-law fluid model for non-Newtonian fluids:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\sigma_{ij} &= -p\\delta_{ij} + \\lambda \\varepsilon_{kk}\\delta_{ij} + 2\\mu \\varepsilon_{ij} \\\\\n", + "\\varepsilon_{ij} &= \\frac{1}{2}\\left(\\frac{\\partial v_i}{\\partial x_j} + \\frac{\\partial v_j}{\\partial x_i}\\right) \\\\\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $\\sigma_{ij}$ is the stress tensor, $p$ is the pressure, $\\lambda$ and $\\mu$ are material constants, $v_i$ is the velocity, and $\\varepsilon_{ij}$ is the strain rate.\n", + "\n", + "### Boundary Conditions:\n", + "\n", + "- Inlet and outlet: Specify velocity or pressure conditions.\n", + "- Channel walls: No-slip conditions.\n", + "\n", + "### Initial Conditions\n", + "\n", + "Start with an initial guess for velocity and pressure.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $v_i \\in H_0^1(\\Omega)$ and $p \\in L^2(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\sigma_{ij} \\varepsilon_{ij} \\, d\\Omega + \\int_{\\Omega} \\frac{\\partial p}{\\partial x_i}v_i \\, d\\Omega \\\\\n", + "&= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $v_i \\in H_0^1(\\Omega)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f5664e1e", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/stokes_flow_in_a_lid-driven_cavity.ipynb b/chapter5/stokes_flow_in_a_lid-driven_cavity.ipynb new file mode 100644 index 0000000..13f8d51 --- /dev/null +++ b/chapter5/stokes_flow_in_a_lid-driven_cavity.ipynb @@ -0,0 +1,78 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "fbe877e5", + "metadata": {}, + "source": [ + "# Stokes Flow in a Lid-Driven Cavity\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the steady-state Stokes equations for flow in a 2D lid-driven cavity:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "- \\nabla p + \\mu \\nabla^2 \\mathbf{v} &= 0 \\\\\n", + "\\nabla \\cdot \\mathbf{v} &= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Lid: $\\mathbf{v} = (1, 0)$ (lid-driven at constant velocity).\n", + "\n", + "- Other boundaries: $\\mathbf{v} = (0, 0)$ (no-slip condition).\n", + "\n", + "- Outlet: $\\frac{\\partial p}{\\partial n} = 0$ (zero-pressure gradient).\n", + "\n", + "### Initial Conditions\n", + "\n", + "Start with an initial guess for velocity and pressure.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $\\mathbf{v} \\in H_0^1(\\Omega)$ and $p \\in L^2(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\mu \\nabla \\mathbf{v} : \\nabla \\mathbf{w} \\, d\\Omega - \\int_{\\Omega} p \\nabla \\cdot \\mathbf{w} \\, d\\Omega \\\\\n", + "&= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $\\mathbf{w} \\in H_0^1(\\Omega)$.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ce2a125a", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/chapter5/two-phase_flow.ipynb b/chapter5/two-phase_flow.ipynb new file mode 100644 index 0000000..a425d84 --- /dev/null +++ b/chapter5/two-phase_flow.ipynb @@ -0,0 +1,78 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "94e72217", + "metadata": {}, + "source": [ + "# Two-Phase Flow (Immersed Boundary Method)\n", + "\n", + "## Mathematical Model\n", + "\n", + "Consider the incompressible Navier-Stokes equations with the Immersed Boundary Method for simulating two-phase flows:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\rho \\frac{\\partial \\mathbf{v}}{\\partial t} + \\rho (\\mathbf{v} \\cdot \\nabla) \\mathbf{v} &= -\\nabla p + \\mu \\nabla^2 \\mathbf{v} + \\mathbf{f}_{\\text{IB}} \\\\\n", + "\\nabla \\cdot \\mathbf{v} &= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $\\mathbf{f}_{\\text{IB}}$ is the immersed boundary force.\n", + "\n", + "### Boundary Conditions\n", + "\n", + "- Inlet and outlet: Specify velocity or pressure conditions.\n", + "- Immersed boundary conditions to model the interface.\n", + "\n", + "### Initial Conditions\n", + "\n", + "Specify initial conditions for velocity, pressure, and the immersed boundary.\n", + "\n", + "## Weak Formulation\n", + "\n", + "Find $\\mathbf{v} \\in H_0^1(\\Omega)$ and $p \\in L^2(\\Omega)$ such that:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "&\\int_{\\Omega} \\rho \\frac{\\partial \\mathbf{v}}{\\partial t} \\cdot \\mathbf{w} \\, d\\Omega + \\int_{\\Omega} \\rho (\\mathbf{v} \\cdot \\nabla) \\mathbf{v} \\cdot \\mathbf{w} \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\mu \\nabla \\mathbf{v} : \\nabla \\mathbf{w} \\, d\\Omega - \\int_{\\Omega} p \\nabla \\cdot \\mathbf{w} \\, d\\Omega + \\int_{\\Omega} \\mathbf{f}_{\\text{IB}} \\cdot \\mathbf{w} \\, d\\Omega \\\\\n", + "&+ \\int_{\\Omega} \\nabla \\cdot \\mathbf{v} q \\, d\\Omega \\\\\n", + "&= 0\n", + "\\end{align}\n", + "$$\n", + "\n", + "for all test functions $\\mathbf{w} \\in H_0^1(\\Omega)$ and $q \\in L^2(\\Omega)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2542c254", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".iga-python", + "language": "python", + "name": ".iga-python" + }, + "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 +} diff --git a/documentation/psydac.rst b/documentation/psydac.rst new file mode 100644 index 0000000..dd90723 --- /dev/null +++ b/documentation/psydac.rst @@ -0,0 +1,19 @@ +Psydac API Reference +==================== + +.. autosummary:: + :toctree: _autosummary + :recursive: + + psydac.api + psydac.cad + psydac.cmd + psydac.core + psydac.ddm + psydac.feec + psydac.fem + psydac.linalg + psydac.mapping + psydac.polar + psydac.pyccel + psydac.utilities diff --git a/documentation/sympde.rst b/documentation/sympde.rst new file mode 100644 index 0000000..3fc651f --- /dev/null +++ b/documentation/sympde.rst @@ -0,0 +1,14 @@ +SymPDE API Reference +==================== + +.. autosummary:: + :toctree: _autosummary + :recursive: + + sympde.core + sympde.calculus + sympde.topology + sympde.expr + sympde.printing + sympde.exterior + sympde.utilities diff --git a/intro.md b/intro.md index 02df7d4..f24f164 100644 --- a/intro.md +++ b/intro.md @@ -2,5 +2,15 @@ Welcome to IGA-Python, a tutorial for Isogeometric Analysis using Python. +## TODO + +- Subdomains +- Geometry and Meshes +- More models +- Visualization +- Documentation using docstrings from sympde and psydac +- GLT +- DOcumentation using docstrings from gelato + ```{tableofcontents} ``` diff --git a/logo.png b/logo.png index fc529c5267f0dfe1373be6b7a62999611728f731..de817f9ff91c15bec557a279fa697dfaa870d4d0 100644 GIT binary patch literal 10614 zcmeHt2T)UOw=Sa6dzIb7>PlJ?*kq8F|hg9?K9U~l^ zO9YZFVyfpWZc$+uh!>1o><+i6C=d!2bN~UJ zJ%Mgsg6@tmEQF=IK^#C%AV-ID9#LUY@ms>8w?xHEL`1m7l*P|}M1=*VL`5ym`GL-k z?iT`eApsyT7|1QEE+Ht46~$|QwnqiKGPUz|hk2i`bwKXkept*;TGCch{2cT2!HN|T z6%iEX!%{RH?Oh>QDWbw6!q`g$3#o%V9i1T9ZR)?d&Q<$3dU}B%?iZGd35p0xo)a8l zK<9H*%h~6iho})$#VkP031A}P=BIP-oCoF)bv);BhCrObj_0VCJrKM69UiROvv}A+ zJRKZ8&mpW+qTFI?+``IO!r3LJdJ#bEjXTioVy6K%CburOgDwK;ZlbPi>jafDxeJr< zRR(}PRc`+d;9sTvk?H2A;R%Gg=t3ME!C1i#eisIbi;17xN%>9^&rqXn&E_3yDBaPl)f|*_FLfneqRC28k`zy6-JhpHp^EZBeIc%cVdaR29f{}Bz(f3(K`MjW;Z&kE58 z2==~+F1M&880(XY69n6OqQYJ<|BIrP^zc4&_BY{{*YEtFm8ygU6k9{sM*Ks^VopH& zb7_AtRfxAI2wUl*!g`Lre`B?WfFb7{@Q0UMR7^zjHqZfEkN<=uozGPNvCR&<_zh@k zoS9(Y4|Bn`>V@Ihu>OPo7wp`hzsvc*>f-|upn-z~)X7gv3k(bp1G#`;Agun2sA0|g z-?{q76+3Jc|L(>A$lP0%-A@0}=s$eHP z_}@W5Fv$6AEZSr1?`(FdU?(LEgzZD+3xXTS!QrfARD7`!a5RRV*Mcv0>R>TXh_}1L z?*{;Oz+vwbwWs8Faj?olk-1hpIQAs>0T_)Mf3k>_KT@vdH$0r@|*UzdHq)@ zzo~zo5a)!84QG$#^TERRm5mw)haE@rj(bKFsl%d1nO}$6L^jv4DzMjm%1a=_ zd-qJJ*e{)36f}a2mk6KW5I%WmdWF}@%#@y3`!}BwI+yee#vwFUdx&y0^{%Qsmp~&p z`Wyb&hO>bFD@y-Al!j6}!I!$=uVs9#FBO~ufM&(;LrXHhGP|oYv1QmCzIiih?m|_W zq@+96GqC?&@sujUa;4m%uHeh&) zOYix2yUn*8Zq{&sY}a@3>awns&YkMQ5bKgnW2Bba=KD$BU%4k6-Avcv znlZ6ld~4J-!Dr?1-r>1xGESsG?kySttEI2OdVKG%kgAH#2zVvdsstBT+&-Dun0pG} zwVTghaL7`fJ*sQ2PhML4E_5Q}QN~QD!_<1EfDR#&E&WjM<2`v@QmOX!x3T<22;=#Z zM_X(&(WSg`I#wLGG>GV&SV}>iPnqF-wE~=tCDS5+PtO~lzsnkW99S?a+IhG57o%i& znZVTdsTaRgN47n*0P5(jARhJnskZ3tdg&mqD#v2CmAc-iGe$&!s=2gaP67iYrxCe; z)qCv^ddUg!j1xLDdJolPehl+mAJiL8D+ z>_5HrRXWI5#8IT2(5uaqM6Ix#Ce%o>Fooai^{W2u?fR;Yp|Rx2%Ul{R zY1RPl2;w_Y29TZCNg*R@_=vLMay?x(*1$_F4EvRDOU({~P5{&>DB~9Eq`n zOoG)zxUKkv#>%nkdM8MiJKYs&R#cS=x`chgN4V`o8 zQHFiK$8EHw0!?cm_dr0EW5nLb{I$MJ+w}~p@CR4oy;F{Gb?g^Wf!}I28ogqIrtTDx zb65d@uC0qFUFpl_-u18vrj6L*9nG2o4b$Q1=d~3&+fwvf!34veBN$B&I)wFe$C9GD z+8?(dR9zElcNf*rb@Jc&W*S3{nNCf=>_ur>!hhI{wSWBh%BWc8i)c>wH-w*iQfbDe zPRbJ(Hdc&lafpJBO@YqrC*^4Gd$;EsLkoklZ_R(0$)Jz?Y#=*Jix{oeiRGdl|JooF zJ;d|$3I3j-mLf>}30ql~g_IkGaXhPWz&_U|w4J6FzJ?Ny;5ojEUXFr({TzbJ6{Tq3 zwhp<*_lpi?${y`ErP$lt#CQC1=II+6R6$dIVZ(~B#^a-g!a?;W@sJ+N-sGK1E#}Hb z$6)8G{1RESwWkK*J!wOF=$6CO$>$w4i;FKGtkseB1m=(IT`d5p*u<7I#)kD@{V}Ve zK>W0aH18mPR19D9m0-&rGi38)dt5tL7cKqHd&y(YrIX;jnr_&~9t@@Arrgbsm#=}c zlw>#X=aGP2DGA0_bN`07)p@R$z?h<_O}UAI@yg*5I9mGkqlgtNgRUUoS@gYbXFe-O zDna8}#Y36s<~!c$ssj>4{ zKS+NTFL?P`J4PQ~i?L+4DpEo|U(r{>UujzBNImqX@tH2X1Mk_ww9nXK=-R#fowO#@ zHyr2doR!PAKbE`}7!IYG*`>TmTn5opq+URe6WG-H`WJl5v>0jo(g$z%RKnlxnkp=BDPQX`dBw`7%Q8+-XB znbpZCTCmhI{xwWS)yd16QIuC@B%?l=mfXNV#$@S;I-cdtpRHmV42+HRaJ3!g3ck-s z;^C!e+WBL=)$Q?IH;v-8X5|iFYkzQ-7-SWZE}7PFc`d-TZ;jtishr6piH+zjS%yX3 zvRl8$DLBBnrwnGK;rzJk(o-A3M`GWu&Yl+fEhwQy2S5825|yU>>Zf+x0pQ4 z_8;_yHM)*?vyzX!-1R34wztsxt2Ha7cyU~YY?m*BbkA=645C&y&p z+t+BDrp*=5Vgv-%OTH$GM3&w}&99LXq;ADR7oIG_tdWHR7W~`CO%YkAH4Rzn z;t6JVW?jYvt$t({apJ7Mw*hU0S3Y>_#ur&yE$>CS02bQ7B>U-^0`uWg>UA#3@)e<3 z_dFA%Z_CpZI&ox;JY!E3hbfum_9^rx>+g7H!uh4Df>>@8(^SX5&Q9-l{qT_|e{D%w zF;=V^uQKzp2=hq=`yz>4(2ZZ8%Q@`kA{@E=nUE9l9)rXqVMejk)CrFkZmwJJ87J`)Q0Mk&w(PMw!4D%U!nhB^kqk==$-Fa6EF-C)?xigz~C1(LE?o)dCwm7 z&N&_91SdPj*4sWH!xk#Lkn$>AhStVV7$G=A%7T@Y*`H&~flo%#WL= zE1%H%{OJ}G@9~-pgOZ;VD=35>?mxx^>sl^EN0@}PobW}7TS}Ie$m_cHa$E_ttX}XI z^JC17H>}y$m^akOwukROUz*s>NBXl30AwV>{huhKV;_UzxyVh$M8k;jd|cO)Zqo>bL{ zxoUDh>-B-A+=|*)ersRs;#_e1jgVdJR#4NQZR3IAeQC*F1U?gM(1<%}ll9-F!^s2M zmN0wQ$G^goSkkAQw3vi08QgCw%%!_w^9h%1X zWMAG4-#-P0X442?Evu1SO&-ZDvXw7N?HxEJ2{`PKw^bg9Fy!3F_pdPxXiKr|DfZkY z7&U2%FX+Kc<2_0(*JC(YQu?xs(Ek>D0{IX*k@c|NjH!W2#B(*B^buL&`_Js?E`JT4 zBsg9<(~#BH{?Es4mHKTYZy$kji-d!2zwA*N39 z9%kUYL{bdAHvqQ?Fjn@juoWwJccJdics8)LhQV|4@*K_3E&4geWPj96d+RPz%P=^O zinYxrj4NDEJ<4sIWNbg}Ro^j7z|xu|ab=JAD+AZw=|&N7xa?uV3&-?UYQ*6OOlAtb z;06E|JxN!n*Q*@Z{#nwiCjhM78-gJ7j4@aijN6^0UpSG84+$c=R1wqKFoeO4)Q*O? zpwpT?m7-Th#a_n4)k&Fgx8nVAi2%1FO3PO$%RCt5KfKXm{XX?<@>LaY5wzgB6CEed zjJHj8Ox8qp@8{toH94VU;Sf@;KCX%|;>oAF>j_qLUQ9wCpvv!tN*%r#F0rhvTQveo zxVFrs>6X+^!USC2p%Q?+3`TBEHwj8(7~X#kZoob6uYk^uQulfBf^omCQAk~)SxXKeiAJTngzJ6l z>~!ZM*pjC}v{MM`>kF~;6K96DZf$Lhh2vKp(YAH_a+McbDI41GC?Qp{f4a#Tld-un zf~ITy1%qF97_-HA%Y?Z}pee*=L32DkC<>IR85wX@yc@)8HBW)4Mv&>y9^UCPTg;)a zZdIQit?s5aVCV?q3T{!{;(C^NlHW|BHv3L*n?x~z^$lOwj$FsBCTxk6!?9fEv$6tnSh`N+NBtx0gqTohR>G`J;&xPnb24GNP%# zlHR2%0yx>NYcD+RQ#f)6s%p%kjt>dw5SyL^u52E<8fgAU(5u_a$9A96Nv@KTWZbCl z{)ocSOKTuH>9*jHG_;F=$_A&^vi$OB=-{xDc^0YgG~iaijQ+G={1N5yzSKd2gkd8S z*XvRd$9V>Q;gE(0Lp!ccrlgSNRtqB3*t|7C{Wuri9>uLkWEy-L6K@cCK&440@12kB zDV-^sX>pI;E43A7h@=|{^f(teedcLt5WGX6;3$+S9-1iTSq}M=2-P;_b`!_z=rF+j z2TSvl4|~PiXn)%V7`(JrbI4WnenGOk*W{~gH46d*oSJlt(S~1e7w}z+{p$6vtd0&pe zTW1dRF|ahM0aspK8|d4v%&95QW_!byBmdjWa#+<@187U??LoNBcxpYuy0b>eT#MDA z#oy}P;^WD3(j{Y3pV=zie9snD-pGrJ1U+RkY{@1%N?5`SCRREDLig?n0IecF3vKtB z-?%hjK_;Vw-aXvEPb@?K;0+-_S7Y%w;#p~2(`QZ=>-7&2Nh!xq_!Om{W-Pwf4MWlX zsF6FQi7eGYMP_BLaaJnqe(*J`Q*dcrcG-xb6TCuIZkgHl^lC|`GiQr_YpI)R1!Gxk zT{qQn)J8#cDU>j1^Jl=GtquICu=T#H&J8wvjf4}*_;Ig0jzmkAB;}*KKWk(;;fb{z zsh#$=0~76Ym{cDW(`q@B4?Ku$ceS-8%Nv;KF8shI)^Llt2V#~>zGBH&0jKVMWkM^X ze1PrY>sF(}ESMz%M)i(%Fyk#V8243WLw(v=q>SF<>;g@>=->kP@(>=nYBiqXA z6Zdeb6bbt_uvnXCT#ul{FiZ5BK*p2h4@qzVO$-|XY?}waWUJCg#cI9@^Sy z{v;>{e7XC6SEf|N0sJywh|X7@cy=(RL5%$rMY8_yoPv2?ui73N!mG!6%K)8>Ja%*$HM( z)BEio8Vx4<^ulL*@8cu9(A$oFU0F2Sx0KOY>&F@#RyFW>YJ4)~Df-2uBs1hY;9YL# zv}t?eAGqUUQP2A*2U1GoF6mR=LZ^fdz4Jg=H38TdI81#ar(SCpwSI7KSz0QY{Qf~M zd8^wMURlB;a&o`Z^hT~>OK$!6*Pdvi2x1MF-~X=J@5|^v<`ZQdb77^J8GE|8Yy8aj zi)H!By)j{c)lJHI9veVLq;?YSux!gtTi}DVum&ohs@~56g92c4nq@JMz!TW*ReH!~ zqD!EBxR{2&?M+5747EJuoYiLa(WVe09QlAzx5AC6X(jtb$XGN>Z8ciAa7y*2^qu17 zG&wO)q`*>h)5a$GN`Jyra#X{|=ap*%#!7c(x$drZEEMH?sX&ZRDr{)eaOj;2x|Nt~ z!Z1ZSqNo()Rd^y1YB)||5$u&@D@r_>l3S$IoMe7Kcr6&?)mG){AkOv$&D7J^57NofQ=TVK2;mKUIGE#` zvKfqh%S(yQlN2rGFoKY5O%S*T(uFaFDC&#~gL*%?wl1~$eL0qtZ4GVX57#qKR9Yg8 zyW%{vv~RuM^TKeE>Ue`geCW|(gP0zvWZ|K^>g9b(TlGT&Z4TDUY%D`MJ%OuCgT;3a zn~OW7JqYZ%JhWBOG3H#rNNwkASq+f+c&yEjEK1wOoe`YT4`#^0)0?L}jJXmWAPVlR z(4wS8zItF+6sY`0Z}ZQUbhF1rPxdB(o5c<}}9`_X)$`aJFQqc&aQ%Mpi_ zVw+Tklh#TA^D_pXDbKV`2x2;3LB)e0#)n~*2d12v zq^N@1_vN>FblQXZQ#BcCSe6L3jgYpFPE6E_>6x(dGe@@40&-WZj(K_{;a?A_`rfil zEl*VZ8nZ82J;};}s#R42vR^^f=IAw%)eQnthg7QN{oeHV16;BnMlEOG-SR3E7^mPt z3PhALk@m<{$nAJctKO^enO!6^XWMuislcOled3bwY*gQUDXYgTc?-4)i(pDb&fBWy z1x(g^Oi&l5_A%+FBcST$cfVAdBf<~g=XgK{)2fy5>?_@(mg*^T zN`IBf z>4b#l>exfhdy!Shq$Hdn9{o;}L8|NuxDcoLv6E-8d{fMo7%n9=e*(R6lF}!_FYYVu zI$f&``y$=jSv)3*^UOAtZ6=$kVyaY?PjDhXqqOXw-XIS(=`5LzT(NPD5h!kg)mvWU zA-@$&NC@V}9vX5#QMyX7#Jg^4#HOEXqlCO3@rhaj+9xaRB6~&dkhRTa1)5OJZ=e|% zDSj+TsSirEU?+ zVc|QvpLtV^&?F2evXm)E$_VQ1e0{ObVu|J}C*zt~jv8RCLm TQwFf#Q{rf<>D@u7+z+eXUjrLKmb6?l;8D^motP^SkH%&bj~0IXUruKhN`ip5=X>=S?E=yqTe( zz##!19v(rXvuA$g;n_k3*17!Kf%kc#JT)F3zD$Dt1p z$g8WW*l*~=JP_VoLQ`xY5{-sQK@6bsV1Vk_1w}b;1C-+9$2gg8nJ=;y1{_u zdU^n}94=QZ9*)3oKtNCsDJ6X=unr*Lyp;5~2?QRzVHobt^FU7IxZ9E3MC$7M_`B#? zdb*pKnH#D4qX;mQ3;)Wt$%!Eu;$b*XQ!E^T28iJ#u0zU7${Rkx12zy>R-IzonFe+9Kfu&kbEQDAzp>einiB@Z7Wv9REgq7}%yV7x49j!Lj}u``4kr zp2|iZH?xVyVgbIt9NJ&Sn@Sjx9s1<{d!E0G(Mu6_9uCF1lZ=hgus|iGCz5~!sQzC=-3-MQNd5Ow{96HY z#?7CJhY%1XPC`8i7&K7moJs0Sz++JeJ?tJV+JJDd zKg)Z);{WF!`?J#w5m*d@fDZs@xr`7+Zhvs+nxYc-+5gWep}dKC{&hmAaus3R$>i|| zH^=LLr+Wb%Zu9@0|4Mp*&fiRtb=i0G`mZeO(#;!UL%`j@xt2Eui`>JoV>~><_l?fz zSdpBWgS%rLtzk`{mxqq-IeUfO+xy87^qa)X8T@5~JwfMnO00JA>42hsejq$xzL&ABUEIVflp4wq3$*ZP>8|wC6O>mRq2$xB84N_X11M9xzBm zbjwzLiOZL5w}KLE2ld_7aen!-4dnt5zkESuRO6A$JX`tqIHI=!sythrf7s*bOcaUM z=TJyg^lazwf~xmPR8&vOABqB0B~Igt1UahwB&pN5fKDf6VGiNt5m(W5yPu3Nk9^6k z@;u40%|m;x!cmnNiQ}P-S^K}auOT(uli~kx zNvPF>4rx}UDM??jv1iukkh(5Ty-GWLac4fY3FwWc6hr#nNBWmXa`mHUTco#yqg3V% z8Zj^Z3Bw0LG+ug6*CAiq49X(=Wo5ORibeC;kf6f`ZI*q%n5d4VoJW)CcQt@k?gPY-8mxutR zRXk7So>O1i?6ejX>_@MK@XL#wtDQ;#y3#43vw%k2L7l?%kgR6wauDM-bPcgAc42!T zfo8QWxHgZ~QQtRfb)qe4I)lrLGRNr}ua?zl z*3rVvEl?nf)cFbJE)S#bkE)JR4P9KjDl`BbH%3cJ`L@TV?uL&Ykj5M1kxHDh^t z`9reFAdqG7*};B6>5XuTwz%LzPuwX*rTnyF=hy@!|cg~#KcUsziDmMb!Z-= zU+7E{Hs!N@D0!mvoeoq&Sy}^9NIU4&Ts1$$RumVC?ny`BOz9t!PEfv>Xqfg;MzGQq zxgfV9bzp9_=FAWr64|XaGTZN1kKC*=K$p0tuf!V4$OUA!j``1d(~0?#P4Pv$=j7}O zsZXZw*&5G|X$GozCHl*}TFjv~-?u9o5H0)ZzDr?Mdh}Q8OdUimc};uYy^+VG!(G)E z)I4EBj+fZC`wy9~FbTNmo=>A;3vuN#b)_YFtbuj|?O*%xSu9DR;Ky|o20dDp#g3;= z4Z7@tQX)#9xGnYDdKXxIeW&`%<8SGXpV;3dcD|M`QQhC%UYdjTY#OcUq!tT?w@toS zQ0QnT7x*v^k00EIY2E#?N~Nx$$CF~Flj?;oJUi3&sBPweKt$VWAdx;|)qlF-b?fXK zmb^_GGNpT(?#chYq;mXkrpKZZS zk}~7rxF}si%c9l>>B3}z=45Ty>mL@fMF$G0<4?iL!wPj+CLi&w>b~t|l6!ZF96#I` zpA$flS#yQbtarx7q^5cH$X&yIxqd$4i}Z@2pVxyyaf_<8>`|Vf&P$Y-QEz#^l+s(m zeGWFUtEo97eW>t{0ga8CylIc#%Ja5oO70N3uTC2K;zHI|@K%z1gNe+nyOBQ0VkX`7 zMY=dNZ`&Vk)h!!wC5yEz3~faX^oo{exdzRJ6i|I@e~y~EuT*Gf0T%=DQ#B-B#~tw! zpi+m-y>>4zJ-b+C8o@-JNMtuIQ5HoC?3{bj5@KF^ze!zP-bW*jraF4h>R8M*q9@O^Jw!&0}=OGyiijKZzo+&fRBn&aW>0;}QSbs5z|VNZ^fd`!0V zON55@c2mW~o%tT}(1@q~g2v*gQWeUcq2>_`smQP=r#!pQ1eK>GgyxgY$<4wE+jRuz zn%{>`O(2R@dPetNkMZ&oo}MwF*BWp8&{Cbdn>ust?GU}m9@foxF}m$}9aIy7-A;Ig zAzT`@4gO83r|WGpns|8f+O$|->VXMXwd=U}ith5%5*DWIYVh# zh*8P(P}7iM0~+)S)Kt`X7MZBP7#5dbp%9mG%0A)N%MUmnjJ81vbl>ZSF04rktZ+5kQHE_!MA#G%eA2@BT6;Rg*aDfjO^@Jmad z^%hIMd@=3#ss7;nC(gbateP(=gRBWp*;-`Pw~jl_0ozbXt2G@OPy%MLwLRvjxu6EG zEj%K5sXzG_H7fPZ50HY`zQdvJ4Hp9#(o^ipp*-k{>q?J5I-{U8z%I+-aEW(+iFbj# zT;5WjwafSi8bW58eR8?4O{6blpdyI%J2S_B*m7uESfaw+e+}MO*%zL-tRv<_Il&iQ zsycD@-EN>=zl9V^29=>cp$2S+ZtGq&QU|G=dXnzb7L}AEeAu2n&{tdgywxzVIsDtl z;mW$~An^pV2;{l4B;^jgN2Q_$WkHWOcrjA)4s?8>TyAY@P5h^SUzq_&|wd)#lI3iic&J@AtC#c$CTLmC78w+3?Mp5sMmB-70WD4!HR9r{mYs3_ABY zGIWe0SQL%Xz_p3Xr#`PzKO|+DCqfC{0mamU0ZaSrj`vw5tW|}+pio2PL8rR`CZYW| zf;@E5_B2Fv>K9m-=jpDuE0G0PV`g?dl{m(iw0%kN(8T*ps8iKcjhP`<^(Mvt(|9Fq z$@eaW-tjuM8#S|HKDFOb<52M3?AGA~s`(|+0(wWsas#aAvD^HsRM%v>kmt?BISa>1 zl|?5(msdB3Aa zQ)m-Qsz64#y_Djyysn$%ByUMFJs*(j1zeUE$@->hBN#oE$ zZEA83R$N1k{c*>vgCZEY_?#Ry=UL{23$DHyjBNE4zpIBQTT$y zY3Bk`VwNf^^-k1uXt!*IbdUx+me$!@th2JZdPcOipuP4{#}YYBJ#i%F8stmus$kjk!k8Lx4XS3isM)>&S__tZpd`af4?#W>3pFZi{_n`mDk zIYLb|1KV17RTh2${RmnnpsCY zdGlPTi-YC+-^BQtw_kjI@{HI_8hts79F0uf3R%iE49a!#k53i_PIzFVsWM+*j%IzB zq){o}UQe*W?2^QrLL#~u#R_KQ&XL<_q2PS#FoE*MdK+59=9C(lv7mOb@b zIc4Np2FoJaKc@PRRkz2_Zj(**Lu2i%m`mM~$Kr?PqB0u>Qpvxo9=HXfIq=yhuS)jr zV&>n^6fJC?&7lM_s~4NRG}y%;+fuDYJaS9;6wi?KhvIfQDHN&>N zR%%FwP~@ZbyN=hmPLD5{?mfLekjU!+9^bCv9YURyG!4^@scz%5M=8Z z89YjMV0-60{T6w!@)!Rw>E=_yb>$C>%a40tg16yE4PeW}!ExJKX$^ZbV%>`tD0KyM zkni<4s5^|_E!0)wt&M_cT)$+zW;riO8fDNW--WXpRco@-KpAaUOKocj9mCT%Li1z~ z+?3IcolWZiMbedp%x{V34x^;Q^VH6hQ!#2ucG4YlGw<)J4w{up3sFU*<3InF%6Fmt zc=eSpFJ!?nkLGTF2;|lOj1+IsOExxVDp0P)Tz#J$HIBlmD?ZuvdGFY5^@+JT`svBe z)TQ}TQmq>|Z2Qi#f>c^5_H>MS_wza-;-?R-WxOeU%Z_WK*`mfQ7shp_Qs6*))4re? z)`aewy5+8c$VAffE%bYB6Pj7WFwtI%F{i;Pkz!P=J-iC7yN>~ zuF!*AS(;--RLU(^T^UcRQ|c|Psx0@=xPEdf3m8X;@;9DoywYs5gpkQB+^=-SxRBB3 z&UD=CdGFF#eimy-V$y{}y|r7TY?`eho=@pALJnf&9Q zR9mwbI)p*)Re+)Xpz8|(3gFMz_OU3xLsJ4r^{ZI2F9alR7xB>I6xBM#NLw>mm|l0S zm5Du+m8d0_wA)ir;&5-{HuDImLA}%E5xnX{5a&0bu$_R-okckU>%Rm2BlssJ{}jW2 fv_D4Hc&sg2e*9p1=^pUkWga7avoloPi{bwT>szXh