{
"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": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"vel =np.sqrt(uh**2 + vh**2)\n",
"plt.contourf(vel,cmap='hsv')\n",
"plt.quiver(u,v,color='k') "
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"divu = div(u,v,dx,dy)\n",
"plt.imshow(abs(divu), origin='bottom')\n",
"plt.colorbar()"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"divu = div(uh,vh,dx,dy) - dt * laplacian(p,dx,dy)\n",
"plt.imshow(abs(divu), origin='bottom')\n",
"plt.colorbar()"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"CSS style adapted from https://github.com/barbagroup/CFDPython. Copyright (c) Barba group\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import urllib\n",
"import requests\n",
"from IPython.core.display import HTML\n",
"def css_styling():\n",
" styles = requests.get(\"https://raw.githubusercontent.com/saadtony/NumericalMethods/master/styles/custom.css\")\n",
" return HTML(styles.text)\n",
"css_styling()\n",
"\n"
]
}
],
"metadata": {
"anaconda-cloud": {},
"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.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}