{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Torricelli Fluid Flow

\n", "\n", " A tank has a small hole in the middle of its base with water flowing out. We find the streamlines for this flow by solving the Navier-Stokes equation for the stream function $u(x)$, which determines the velocity, and the vorticity $w(x)$:\n", "

\n", " $$\n", " {\\mathbf{v} \\stackrel{def}{=} \\vec{\\nabla} \\times\\mathbf{u}(\\mathbf{x})}\n", "= \\hat{{\\epsilon}}_x\\left( \\frac{\\partial u_z}{\\partial y} - \\frac\n", "{\\partial u_y} {\\partial z}\\right) + \\hat{{\\epsilon}}_y\\left(\n", "\\frac{\\partial u_x}{\\partial z} -\\frac{\\partial u_z} {\\partial\n", "x}\\right), \\quad \\quad \n", "\\mathbf{w} \\stackrel{def}{=} \\vec{\\nabla} \\times\n", "\\mathbf{v}(\\mathbf{x}).\n", "$$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iteration 0\n", "Residual r1 5.633038634729485\n", "Residual r1 0.8853692561741199\n", "Iteration 100\n", "Residual r1 0.35031269473899446\n", "Residual r1 0.19667393154684307\n", "Iteration 200\n", "Residual r1 0.12962246620056989\n", "Residual r1 0.09354483120667396\n", "Iteration 300\n", "Residual r1 0.07156140704090319\n", "Residual r1 0.057013398963852074\n", "Iteration 400\n", "Residual r1 0.04680239465082821\n", "Residual r1 0.03931231059896448\n", "Iteration 500\n", "Residual r1 0.03362640756809015\n", "Residual r1 0.029190036488446938\n", "Iteration 600\n", "Residual r1 0.025650036011805356\n", "Residual r1 0.02277201872496164\n", "Iteration 700\n", "Residual r1 0.020394998885092264\n" ] } ], "source": [ "\"\"\" From \"COMPUTATIONAL PHYSICS\" & \"COMPUTER PROBLEMS in PHYSICS\"\n", " by RH Landau, MJ Paez, and CC Bordeianu (deceased)\n", " Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, \n", " C Bordeianu, Univ Bucharest, 2020. \n", " Please respect copyright & acknowledge our work.\"\"\"\n", "\n", "# Torricelli.py: solves Navier-Stokes equation for orifice flow\n", "\n", "import numpy as np # Need for zeros\n", " \n", "Niter = 700; Ndown = 20; Nx = 17; N2x = 2*Nx; Ny = 156 \n", "Nb = 15; h = 0.4; h2 = h*h; g = 980.; nu = 0.5; iter = 0; \n", "Vtop = 8.0e-4; omega = 0.1; R = Vtop*h/nu\n", "\n", "u = np.zeros((Nx+1, Ny+1), float); ua =np.zeros((N2x,Ny), float) \n", "w = np.zeros((Nx+1, Ny+1), float) \n", "\n", "Torri = open('Torri.dat','w'); uall = open('uall.dat','w')\n", "\n", "def BelowHole():\n", " for i in range(Nb+1,Nx+1): # Below orifice\n", " u[i,0] = u[i-1,1] # du/dy =vx=0\n", " w[i-1,0] = w[i-1,1] # Water is at floor \n", " for j in range (0,Ndown+1): \n", " if i == Nb: vy = 0\n", " if i == Nx: vy = -np.sqrt(2.0*g*h*(Ny+Nb-j))\n", " if i == Nx-1: vy = -np.sqrt(2.0*g*h*(Ny+Nb-j))/2.\n", " u[i,j] = u[i-1,j]-vy*h # du/dx=-vy \n", "\n", "def BorderRight():\n", " for j in range (1,Ny+1): # Center orifice very sensitive \n", " vy = -np.sqrt(2.0*g*h*(Ny-j)) \n", " u[Nx,j] = u[Nx-1,j]+vy*h \n", " u[Nx,j] = u[Nx,j-1] \n", " w[Nx,j] = -2*(u[Nx,j]-u[Nx,j-1])/h**2 \n", "\n", "def BottomBefore():\n", " for i in range (1,Nb+1): # Bottom, before the hole\n", " u[i,Ndown] = u[i,Ndown-1] \n", " w[i,Ndown] = -2*(u[i,0]-u[i,1]) \n", "\n", "def Top(): \n", " for i in range(1,Nx): # Top \n", " u[i,Ny] = u[i,Ny-1] \n", " w[i,Ny] = 0\n", "\n", "def Left(): \n", " for j in range (Ndown,Ny): # Left wall\n", " w[0,j] = -2*(u[0,j]-u[1,j])/h**2\n", " u[0,j] = u[1,j] # du/dx=0\n", "\n", "def Borders(iter): # Method borders: init & B.C.\n", " BelowHole()\n", " BorderRight() #right (center of hole)\n", " BottomBefore() # Bottom before the hole\n", " Top()\n", " Left()\n", " \n", "def Relax(iter):\n", " Borders(iter) \n", " for i in range(1, Nx): \n", " for j in range (1, Ny):\n", " if j <= Ndown:\n", " if i > Nb:\n", " r1 = omega*((u[i+1,j]+u[i-1,j]+u[i,j+1]+u[i,j-1]\n", " \t +h*h*w[i, j])*0.25-u[i,j]) \n", " u[i,j] += r1\n", " if j > Ndown:\n", " r1 = omega*((u[i+1,j]+u[i-1,j]+u[i,j+1]+u[i,j-1]\n", " \t +h*h*w[i, j])*0.25-u[i,j]) \n", " u[i,j] += r1 \n", " if iter%50==0:\n", " print(\"Residual r1 \", r1)\n", " Borders(iter) \n", " for i in range(1, Nx): # Relax stream function\n", " for j in range (1, Ny): \n", " if j <= Ndown:\n", " if i >= Nb:\n", " a1 = w[i+1, j]+ w[i-1,j]+w[i,j+1]+ w[i,j-1]\n", " a2 = (u[i,j+1]-u[i,j-1])*(w[i+1,j]-w[i-1,j])\n", " a3 = (u[i+1,j]-u[i-1,j])*(w[i,j+1]-w[i,j-1])\n", " r2 = omega*( (a1 + (R/4.)*(a3 - a2) )/4.0-w[i,j])\n", " w[i,j]+=r2\n", " if j > Ndown:\n", " a1 = w[i+1, j]+ w[i-1,j]+w[i,j+1]+ w[i,j-1]\n", " a2 = (u[i,j+1]-u[i,j-1])*(w[i+1,j]-w[i-1,j])\n", " a3 = (u[i+1,j]-u[i-1,j])*(w[i,j+1]-w[i,j-1])\n", " r2 = omega*( (a1 + (R/4.)*(a3 - a2) )/4.0-w[i,j])\n", " w[i,j] += r2\n", "\n", "while (iter <= Niter): \n", " if iter %100 == 0: print (\"Iteration\", iter) \n", " Relax(iter)\n", " iter += 1 # counter of iterations\n", "for j in range(0,Ny): # Send w to disk in gnuplot format\n", " for i in range(0,Nx):\n", " Torri.write(\"%8.3e \\n\"%(w[i,j]))\n", " Torri.write(\"\\n\") \n", "Torri.close()\n", "for j in range(0,Ny): # Send symmetric tank data to disk \n", " for i in range(0,N2x):\n", " if i <= Nx:\n", " ua[i,j] = u[i,j]\n", " uall.write(\"%8.3e \\n\"%(ua[i,j])) \n", " if i > Nx:\n", " ua[i,j] = u[N2x-i,j]\n", " uall.write(\"%8.3e \\n\"%(ua[i,j])) \n", " uall.write(\"\\n\") \n", "uall.close()\n", "utorr = open('Torri.dat','w') # Send u data to disk\n", "for j in range(0,Ny):\n", " utorr.write(\"\\n\") \n", " for i in range(0,Nx): \n", " utorr.write(\"%10.3e \\n\"%(u[i,j])) \n", "utorr.close()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.13" } }, "nbformat": 4, "nbformat_minor": 2 }