\n", " A 2D rectangular cavity has an inlet of fluid at the left hand side ans an outlet at the right. Here the lines of fluid are visualized by solving the stream function u(x) from which the velocity is determined by the curl operator is defined as:\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", " And solve the vorticity defined as the curl of the velocity:\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.\n", " " ] }, { "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": [ "# solves Navier - Stokes equation for flow around beam, sections 19.9\n", "from numpy import * # needed for zeros\n", "#use hoffman vol 1 b c\n", "import numpy as np\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", "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 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", "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", "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", "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", " \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", "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", "def relax(iter):\n", " borders(iter) \n", " for i in range(1, Nx): \n", " for j in range (1, Ny):\n", " \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", " \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", " \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.4" } }, "nbformat": 4, "nbformat_minor": 2 }