\n", " The Navier-Stokes equation is solved for the lines of fluid flow through a 2D rectangular cavity are visualized by solving the stream function u(x), the velocity, and the vorticity:\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).\n", "$$\n", "$$\n", "\\mathbf{w} \\stackrel{def}{=} \\vec{\\nabla} \\times\n", "\\mathbf{v}(\\mathbf{x}).\n", "$$\n", "
\n", " See:\n", "Klauss A. Hoffmann, Steve T. Chiang, (2000), 4th Ed., Computational Fluid Dynamics,\n", "Vol I, Engineering Education System, Wichita." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0\n", "r1 -2.698224852071006\n", "100\n", "r1 0.05069518708248211\n", "200\n", "r1 -0.001940431667948417\n", "finished\n" ] } ], "source": [ "# CavityFlow.py: Navier-Stokes eqn for flow through cavity, hoffman vol 1 \n", "\n", "from numpy import * # needed for zero\n", "import numpy as np\n", "\n", "#Instructions: use gnuplot to plot\n", "# splot \"uhoff.dat\" w l\n", "# unset surface\n", "#unset key\n", "#set contour\n", "#set cntrparam levels 40\n", "#set view 0,0,1,1\n", "#replot\n", "\n", "Niter = 200\n", "Nx = 30\n", "h = 1\n", "Ny = 30 # grid params\n", "J1 = 15\n", "J2 = 20\n", "J3 = 6\n", "J4 = 11\n", "u = np.zeros((Nx+1, Ny+1), float) # stream \n", "w = np.zeros((Nx+1, Ny+1), float) # Vorticity\n", "V0 = 18. # velocity at inlet\n", "omega = 0.1 # relaxation param \n", "h2 = h*h\n", "nu = 25 # viscosity\n", "iter = 0 # Number of iterations\n", "R = V0*h/nu # Reynolds Number, normal units\n", "u0 = 3.\n", "\n", "def top(): \n", " for i in range(1,Nx+1): # TOP \n", " u[i-1,Ny] = u[i,Ny] # vy = -du/dx = 0\n", " w[i,Ny] = 2*(u[i,Ny]-u[i,Ny-1])/h**2 -2*u0/h\n", "\n", "def bottom():\n", " for i in range (1,Nx+1): # BOTTOM\n", " w[i,0] = 2*(u[i,0]-u[i,1])/h**2\n", " u[i,1] = u[i,0] # du/dy = 0\n", " u[i-1,0] = u[i,0] #du/dx = 0\n", " \n", "def borderleft():\n", " for j in range (1,Ny+1): # RIGHT\n", " if j < J1: \n", " w[0,j] = 2*(u[0,j]-u[1,j])/h**2 # below hole\n", " u[0,j] = u[1,j]\n", " u[0,j] = u[0,j-1]\n", " if j >= J1 and j< = J2:\n", " w[1,j] = w[0,j]\n", " u[0,j] = u[0,j-1]+V0*h\n", " if j > J2:\n", " u[0,j] = u[0,j-1] #du/dy = 0\n", " u[1,j] = u[0,j] #du/dx = 0\n", " #w[0,j] = 0.\n", " w[0,j] = 2*(u[0,j]-u[1,j])/h**2 \n", " \n", "def borderight():\n", " for j in range (1,Ny+1): # RIGHT\n", " if j < J3:\n", " w[Nx-1,j] = 2*(u[Nx-1,j]-u[Nx,j])/h**2 #below hole\n", " u[Nx-1,j] = u[Nx,j]\n", " u[Nx,j] = u[Nx,j-1]\n", " \n", " if j >= J3 and j <= J4:\n", " u[Nx,j] = u[Nx-1,j]\n", " u[Nx,j] = u[Nx,j-1]+V0*h\n", " w[Nx-1,j] = w[Nx,j]\n", " if j > J4:\n", " u[Nx,j] = u[Nx,j-1]\n", " u[Nx-1,j] = u[Nx,j]\n", " w[Nx-1,j] = 2*(u[Nx-1,j]-u[Nx,j])/h**2 \n", "\n", "def borders(iter): # Method borders: init & B.C.\n", " top()\n", " bottom() #bottom \n", " borderight() # at tight (center of hole) \n", " borderleft()\n", " \n", "def relax(iter):\n", " borders(iter) \n", " for i in range(1, Nx): \n", " for j in range (1, Ny):\n", " r1 = omega*((u[i+1,j]+u[i-1,j]+u[i,j+1]+u[i,j-1]+h*h*w[i, j])*0.25-u[i,j]) \n", " u[i,j]+ = r1\n", " if iter % 100 0:\n", " print( \"r1 \", r1)\n", " borders(iter) \n", " for i in range(1, Nx): # Relax stream function\n", " for j in range (1, Ny): \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:\n", " print (iter) # iterations counted\n", " relax(iter)\n", " iter += 1 # counter of iterations\n", "\n", "utorr = open('uhoff.dat','w') # send data to disk of u\n", "for j in range(0,Ny+1):\n", " utorr.write(\"\\n\") \n", " for i in range(0,Nx+1): \n", " utorr.write(\"%10.3e \\n\"%(u[i,j])) \n", "utorr.close()\n", "print(\"finished\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }