{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Navier-Stokes Solver using Finite Differences\n", "## CH EN 6355 - Computational Fluid Dynamics\n", "**Prof. Tony Saad (www.tsaad.net)
Department of Chemical Engineering
University of Utah**\n", "
" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "%matplotlib inline\n", "%config InlineBackend.figure_format = 'svg'\n", "import matplotlib.pyplot as plt\n", "import matplotlib.animation as animation\n", "from matplotlib import cm\n", "plt.rcParams['animation.html'] = 'html5'" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "'''\n", "Compute first derivative using central differencing on f\n", "'''\n", "def ddx(f, dx):\n", " result = np.zeros_like(f)\n", " result[1:-1,1:-1] = (f[1:-1,2:] - f[1:-1,:-2])/2.0/dx\n", " return result\n", "\n", "def ddy(f, dy):\n", " result = np.zeros_like(f)\n", " result[1:-1,1:-1] = (f[2:,1:-1] - f[:-2,1:-1])/2.0/dy\n", " return result\n", " \n", "def laplacian(f, dx, dy):\n", " result = np.zeros_like(f)\n", " result[1:-1,1:-1] = (f[1:-1,2:] - 2.0*f[1:-1,1:-1] + f[1:-1,:-2])/dx/dx \\\n", " + (f[2:,1:-1] -2.0*f[1:-1,1:-1] + f[:-2,1:-1])/dy/dy\n", " return result\n", "\n", "def div(u,v,dx,dy):\n", " return ddx(u,dx) + ddy(v,dy)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Re = 20.0\n" ] } ], "source": [ "nx = 32\n", "ny = 32\n", "lx = 1.0\n", "ly = 1.0\n", "dx = lx/(nx-1)\n", "dy = ly/(ny-1)\n", "\n", "ν = 0.05\n", "Ut = 1.0 # m/s\n", "\n", "dt = min(0.25*dx*dx/ν, 4.0*ν/Ut/Ut)\n", "print('Re =', Ut*lx/ν)\n", "u = np.zeros([ny,nx])\n", "v = np.zeros([ny,nx])\n", "uh = np.zeros([ny,nx])\n", "vh = np.zeros([ny,nx])\n", "p = np.zeros([ny,nx])" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "def pressure_poisson(p, dx, dy, b):\n", " pn = np.empty_like(p)\n", " it = 0\n", " err = 1e5\n", " tol = 1e-3\n", " maxit = 50\n", " while it < maxit and err > tol:\n", " pn = p.copy()\n", " p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 + \n", " (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2) /\n", " (2 * (dx**2 + dy**2)) -\n", " dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * \n", " b[1:-1,1:-1])\n", "\n", " p[:, -1] = p[:, -2] # dp/dy = 0 at right wall\n", " p[0, :] = p[1, :] # dp/dy = 0 at y = 0\n", " p[:, 0] = p[:, 1] # dp/dx = 0 at x = 0\n", " p[-1, :] = 0 # p = 0 at the top wall\n", " err = np.linalg.norm(p - pn, 2)\n", " it += 1\n", " \n", " return p, err" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "t = 0.0\n", "tend = 4000*dt\n", "\n", "while t < tend:\n", " # set boundary conditions\n", " # bottom wall\n", " u[0,:] = 0.0 \n", " v[0,:] = 0.0 \n", " # top wall \n", " u[-1,:] = Ut\n", " v[-1,:] = 0.0\n", " # left wall\n", " u[:,0] = 0.0\n", " v[:,0] = 0.0\n", " # right wall\n", " u[:,-1] = 0.0\n", " v[:,-1] = 0.0\n", " \n", " # do the x-momentum RHS\n", " # u rhs: - d(uu)/dx - d(vu)/dy + ν d2(u)\n", " uRHS = - ddx(u*u,dx) - ddy(v*u,dy) + ν*laplacian(u,dx,dy)\n", " # v rhs: - d(uv)/dx - d(vv)/dy + ν d2(v)\n", " vRHS = - ddx(u*v,dx) - ddy(v*v,dy) + ν*laplacian(v,dx,dy)\n", " \n", " uh = u + dt*uRHS\n", " vh = v + dt*vRHS\n", " \n", " # next compute the pressure RHS: prhs = div(un)/dt + div( [urhs, vrhs])\n", " prhs = div(uh,vh,dx,dy)/dt\n", " p,err = pressure_poisson(p,dx,dy,prhs)\n", " \n", " # finally compute the true velocities\n", " # u_{n+1} = uh - dt*dpdx\n", " u = uh - dt*ddx(p,dx)\n", " v = vh - dt*ddy(p,dy)\n", " t += dt" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "