{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Relax and hold steady" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Stokes flow" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook presents the coding assignment for **Module 5** of the course [*\"Practical Numerical Methods with Python.*](https://github.com/numerical-mooc/numerical-mooc) Your mission is to solve Stokes flow in a square cavity, using the vorticity-streamfunction formulation.\n", "\n", "Stokes flow, also known as *creeping flow*, refers to flows that are dominated by viscous forces and not by the advective/convective forces. The Stokes-flow assumption works well for flows that have very low Reynolds number, much smaller than 1: very slow, highly viscous or flows at microscopic length scales.\n", "\n", "Stokes flow allows us to simplify the Navier-Stokes equations, eliminating the non-linearity. Let's run through a quick derivation of the vorticity-transport equation with Stokes-flow assumptions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Vorticity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start with the Navier-Stokes equations for incompressible flow:\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial u}{\\partial t} + u \\cdot \\nabla u = -\\frac{1}{\\rho}\\nabla p + \\nu\\nabla^2 u\n", "\\end{equation}\n", "$$\n", "\n", "If we scale Equation $(1)$ to make it non-dimensional, we can rewrite it as\n", "\n", "$$\n", "\\begin{equation}\n", "Re \\left(\\frac{\\partial u^*}{\\partial t} + u^* \\cdot \\nabla u^* \\right) = -\\nabla p^* + \\nabla^2 u^*\n", "\\end{equation}\n", "$$\n", "\n", "Where $u^*$ and $p^*$ are the non-dimensional velocity and pressure, respectively. \n", "\n", "To obtain Stokes flow, we assume that the Reynolds number approaches zero. Applying that assumption to Equation $(2)$ and dropping the stars, yields\n", "\n", "$$\n", "\\begin{equation}\n", "0 = - \\nabla p + \\nabla^2 u\n", "\\end{equation}\n", "$$\n", "\n", "That simplified things! Now, we apply the curl operator on both sides of the equation:\n", "\n", "$$\n", "\\begin{equation}\n", "\\nabla \\times 0 = \\nabla \\times \\left( - \\nabla p + \\nabla^2 u\\right)\n", "\\end{equation}\n", "$$\n", "\n", "The left-hand side remains zero, while the first term on the right-hand side is\n", "\n", "$$\n", "\\begin{equation}\n", "\\nabla \\times - \\nabla p = 0\n", "\\end{equation}\n", "$$\n", "\n", "because $\\nabla \\times \\nabla \\phi = 0$ where $\\phi$ is a scalar.\n", "\n", "Finally,\n", "\n", "$$\n", "\\begin{equation}\n", "\\nabla \\times \\nabla^2 u =\\nabla^2\\omega\n", "\\end{equation}\n", "$$\n", "\n", "where $\\nabla \\times u = \\omega$ is the vorticity. \n", "\n", "Combining all of these equations, we arrive at the simplified vorticity transport equation for Stokes flow:\n", "\n", "$$\n", "\\begin{equation}\n", "\\nabla ^2 \\omega = 0\n", "\\end{equation}\n", "$$\n", "\n", "Look familiar?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Stream function" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the stream function $\\psi$, such that\n", "\n", "$$\n", "\\begin{equation}\n", "u = \\frac{\\partial \\psi}{\\partial y} \\text{ and } v = - \\frac{\\partial \\psi}{\\partial x}\n", "\\end{equation}\n", "$$\n", "\n", "In 2D, we can write out the vorticity as\n", "\n", "$$\n", "\\begin{equation}\n", "\\omega = \\frac{\\partial v}{\\partial x} - \\frac{\\partial u}{\\partial y}\n", "\\end{equation}\n", "$$\n", "\n", "which, combined with the previous equation yields another familiar looking equation:\n", "\n", "$$\n", "\\begin{equation}\n", "\\nabla^2 \\psi = -\\omega\n", "\\end{equation}\n", "$$\n", "\n", "We have a system of two coupled equations that can describe the fluid flow in a lid-driven cavity at very low Reynolds numbers. \n", "\n", "$$\n", "\\begin{equation}\n", "\\nabla^2 \\omega = 0\n", "\\end{equation}\n", "$$\n", "\n", "$$\n", "\\begin{equation}\n", "\\nabla^2 \\psi = -\\omega\n", "\\end{equation}\n", "$$\n", "\n", "Note that by substituting Equation $(12)$ into $(11)$, we arrive at a new equation: the *biharmonic equation*:\n", "\n", "$$\n", "\\begin{equation}\n", "\\nabla^4 \\psi= 0\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solving the biharmonic equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Is it possible to discretize a 4th-order partial differential equation? Of course! Are we going to? No!\n", "\n", "There's nothing wrong with a 4th-order equation, but in this course module we learned about the Poisson equation and that's what we're going to use. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cavity flow" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You will solve a problem called *lid-driven cavity flow*. This is a common test problem for Navier-Stokes solvers—we'll be using it to examine Stokes flow. \n", "\n", "Assume that the lid of a square cavity moves at a constant velocity of $u=1$, with no fluid leaking out past the moving lid. We we want to visualize what the flow field inside the cavity looks like at steady state. \n", "\n", "All of the surfaces, including the lid, are assumed to have no-slip boundary conditions. The boundary conditions are all specified in terms of the streamfunction $\\psi$, as shown below in Figure $(1)$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Figure 1. Lid-driven Cavity Flow" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Boundary conditions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One of the major hurdles with the vorticity-streamfunction formulation is the treatment of boundary conditions. \n", "\n", "The boundary conditions are all specified in terms of $\\psi$ and its derivatives, but the Laplace equation\n", "\n", "$$\n", "\\nabla \\omega^2 = 0\n", "$$\n", "\n", "has no $\\psi$ value. Instead, we need a way to represent the boundary conditions for $\\omega$ in terms of $\\psi$. \n", "\n", "Consider the equation $\\nabla ^2 \\psi = -\\omega$ along the top surface of the cavity (the moving surface). The streamfunction $\\psi$ along $x$-direction is a zero constant, so $\\frac{\\partial ^2 \\psi}{\\partial x^2}$ goes to zero and the equation simplifies to\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial ^2 \\psi}{\\partial y^2} = -\\omega\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A 2nd-order central difference discretization gives\n", "\n", "$$\n", "\\begin{equation}\n", "\\omega_j = - \\left(\\frac{\\psi_{j+1} - 2\\psi_j + \\psi_{j-1}}{\\Delta y^2}\\right)\n", "\\end{equation}\n", "$$\n", "\n", "but the value $\\psi_{j+1}$ is outside of the domain. Now take a 3rd-order discretization of $\\frac{\\partial \\psi}{\\partial y}$ evaluated along the top edge.\n", "\n", "$$\n", "\\begin{equation}\n", "\\left.\\frac{\\partial \\psi}{\\partial y}\\right|_j = \\frac{2\\psi_{j+1} + 3\\psi_j - 6\\psi_{j-1} + \\psi_{j-2}}{6 \\Delta y}\n", "\\end{equation}\n", "$$\n", "\n", "$\\frac{\\partial \\psi}{\\partial y}$ is a given boundary value in the problem along the top edge\n", "\n", "$$\n", "\\begin{equation}\n", "\\left.\\frac{\\partial \\psi}{\\partial y}\\right|_j = u_j\n", "\\end{equation}\n", "$$\n", "\n", "which leaves us with a value for $\\psi_{j+1}$ that consists only of points within the domain. \n", "\n", "$$\n", "\\begin{equation}\n", "\\psi_{j+1} = \\frac{6\\Delta y u_j - 3\\psi_j + 6 \\psi_{j-1} - \\psi_{j-2}}{2}\n", "\\end{equation}\n", "$$\n", "\n", "Plug in that result into the initial discretization from Equation $(15)$ and we have a boundary condition for $\\omega$ along the top surface in terms of $\\psi$:\n", "\n", "$$\n", "\\begin{equation}\n", "\\omega_{i,j} = -\\frac{1}{2 \\Delta y^2} (8\\psi_{i, j-1} - \\psi_{i, j-2}) - \\frac{3u_j}{\\Delta y} + \\mathcal{O}(\\Delta y^2)\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Coding assignment" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve Stokes flow in a lid-driven cavity using the parameters given below. \n", "\n", "You should iteratively solve for both $\\omega$ and $\\psi$ until the L1 norm of the difference between successive iterations is less than $1$$\\tt{E}$$^-6$ for **both** quantities. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Set parameters.\n", "nx, ny = 41, 41 # number of points in each direction\n", "L = 1.0 # length of the square cavity\n", "dx = L / (nx - 1) # grid spacing in the x direction\n", "dy = L / (ny - 1) # grid spacing in the y direction" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def l1_norm(u, u_ref):\n", " \"\"\"\n", " Computes and returns the L1-norm of the difference\n", " between a solution u and a reference solution u_ref.\n", "\n", " Parameters\n", " ----------\n", " u : numpy.ndarray\n", " The solution as an array of floats.\n", " u_ref : numpy.ndarray\n", " The reference solution as an array of floats.\n", "\n", " Returns\n", " -------\n", " diff : float\n", " The L2-norm of the difference.\n", " \"\"\"\n", " diff = numpy.sum(numpy.abs(u - u_ref))\n", " return diff" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The final result should resemble the plot shown in Figure $(2)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Hint" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The boundary conditions for $\\omega$ depend upon the current value of $\\psi$. The two equations are *coupled*. If you try to solve them in a *uncoupled* way, things will go poorly." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Figure 2. Contour plot of streamfunction at steady state" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### References" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Fletcher, C. A. (1988). Computational Techniques for Fluid Dynamics: Volume 2: Specific Techniques for Different Flow Categories.\n", "\n", "* Ghia, U. K. N. G., Ghia, K. N., & Shin, C. T. (1982). High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. Journal of computational physics, 48(3), 387-411.\n", "\n", "* Greenspan, D. (1974). Discrete numerical methods in physics and engineering (Vol. 312). New York: Academic Press.\n", "\n", "* Heil, Matthias (2007). [Viscous Fluid Flow Vorticity Handout (pdf)](http://www.maths.manchester.ac.uk/~mheil/Lectures/Fluids/Material_2007/Vorticity.pdf)\n", "\n", "* Non-dimensionalization and scaling of the Navier Stokes equations. (n.d.). In *Wikipedia*. Retrieved January 30, 2015 [http://en.wikipedia.org/w/index.php?title=Non-dimensionalization_and_scaling_of_the_Navier-Stokes_equations](http://en.wikipedia.org/w/index.php?title=Non-dimensionalization_and_scaling_of_the_Navier%E2%80%93Stokes_equations&oldid=641860920)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "###### The cell below loads the style of this notebook." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.core.display import HTML\n", "css_file = '../../styles/numericalmoocstyle.css'\n", "HTML(open(css_file, 'r').read())" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3 (MOOC)", "language": "python", "name": "py36-mooc" }, "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.6.5" } }, "nbformat": 4, "nbformat_minor": 1 }