Vorticity-Streamfunction FD Python

CH EN 6355 - Computational Fluid Dynamics

Prof. Tony Saad (www.tsaad.net)
Department of Chemical Engineering
University of Utah


Synopsis

This notebook implements a finite difference solution of the 2D vorticity streamfunction approach for the constant density incompressible Navier-Stokes equations. These are given by \begin{equation} \frac{{\partial {\omega _z}}}{{\partial t}} = - \frac{{\partial \psi }}{{\partial y}}\frac{{\partial {\omega _z}}}{{\partial x}} + \frac{{\partial \psi }}{{\partial x}}\frac{{\partial {\omega _z}}}{{\partial y}} + \nu {\nabla ^2}{\omega _z} \end{equation} \begin{equation} \nabla ^2 \psi = -\omega_z \end{equation} along with $u = \frac{\partial \psi}{\partial y}$ and $v =- \frac{\partial \psi}{\partial x}$.

The solution procedure consists of the following steps:

  1. Initialize $\omega$ in the domain and set $\omega$ and $\psi$ appropriately at the boundaries
  2. Compute the streamfunction by solving the Poisson equation for the streamfunction-vorticity relation
  3. Solve the vorticity transport equation using an appropriate discretization scheme in time and space

For this notebook, we will use an FTCS scheme (Forward in Time, Central in Space). We have \begin{equation} \frac{{\partial {\omega _z}}}{{\partial t}} \approx \frac{{\omega _{i,j}^{n + 1} - \omega _{i,j}^n}}{{\Delta t}} \end{equation} \begin{equation}

  • \frac{{\partial \psi }}{{\partial y}}\frac{{\partial {\omega z}}}{{\partial x}} = - \left( {\frac{{\psi {i,j + 1}^n - \psi {i,j - 1}^n}}{{2\Delta y}}} \right)\left( {\frac{{\omega {i + 1,j}^n - \omega _{i - 1,j}^n}}{{2\Delta x}}} \right) \end{equation} \begin{equation} \frac{{\partial \psi }}{{\partial x}}\frac{{\partial {\omega _z}}}{{\partial y}} = \left( {\frac{{\psi _{i + 1,j}^n - \psi _{i - 1,j}^n}}{{2\Delta x}}} \right)\left( {\frac{{\omega _{i,j + 1}^n - \omega _{i,j - 1}^n}}{{2\Delta y}}} \right) \end{equation} and \begin{equation} {\nabla ^2}{\omega _z} = \frac{{\omega _{i + 1,j}^n - 2\omega _{i,j}^n + \omega _{i - 1,j}^n}}{{\Delta {x^2}}} + \frac{{\omega _{i,j + 1}^n - 2\omega _{i,j}^n + \omega _{i,j - 1}^n}}{{\Delta {y^2}}} \end{equation}
In [9]:
import numpy as np
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib import cm
plt.rcParams['animation.html'] = 'html5'
In [10]:
nx = 21
ny = 21
lx = 1.0
ly = 1.0
dx = lx/(nx-1)
dy = ly/(ny-1)
In [11]:
Ut = 5.0 # u top wall
Ub = 0.0 # u bottom wall
Vl = 0.0 # V left wall
Vr = 0.0 # V right wall
In [12]:
ψ0 = np.zeros([ny,nx])
ω0 = np.zeros([ny,nx])

# apply boundary conditions
ω0[0,1:-1]  = -2.0*ψ0[1, 1:-1]/dy/dy + 2.0*Ub/dy     # vorticity on bottom wall (stationary)
ω0[-1,1:-1] = -2.0*ψ0[-2, 1:-1]/dy/dy - 2.0*Ut/dy # vorticity on top wall (moving at Uwall)
ω0[1:-1,-1] = -2.0*ψ0[1:-1,-2]/dx/dx # right wall
ω0[1:-1,0]  = -2.0*ψ0[1:-1,1]/dx/dx # left wall

# ω0[ny//2,:nx//2] = 1.0/dy
# ω0[ny//2-1,nx//2:] = 1/dy

t = 0.0
dt = 0.001

tol = 1e-3
maxIt = 30
β = 1.5
ν = 0.05
psi = []
psi.append(ψ0)
omega = []
omega.append(ω0)
# plt.contour(omega[0])
dt = min(0.25*dx*dx/ν, 4.0*ν/Ut/Ut)
tend = 4*dt
# print(ν*dt/dx/dx <= 0.25)
# print(Ut**2*dt/ν <= 4.0)
print('dt =', dt)
print('Reynods =', Ut*lx/ν)
dt = 0.008
Reynods = 100.0
In [13]:
while t < tend:

    # first find ψ by solving the Poisson equation Δψ = -ω
    it = 0
    err = 1e3
    ψ = psi[-1].copy()    
    ωn = omega[-1]    
    while err > tol and it < maxIt:
        ψk = np.zeros_like(ψ)
        ψk[1:-1,1:-1] = ψ[1:-1,1:-1] # set interior values
        for i in range(1,nx-1):
            for j in range(1,ny-1):
                rhs = (dx*dy)**2 * ωn[j,i] + dx**2 * (ψ[j+1,i] + ψ[j-1,i]) + dy**2 * (ψ[j,i+1] + ψ[j,i-1])
                rhs *= β/2.0/(dx**2 + dy**2)
                ψ[j,i] = rhs + (1-β)*ψ[j,i]
        err = np.linalg.norm(ψ - ψk)
        it += 1
    print(ψ)
    psi.append(ψ)

    # declare new array for ω
    ω = np.zeros_like(ωn)
    
    Cx  = -(ψ[2:,1:-1] - ψ[:-2,1:-1])/2.0/dy * (ωn[1:-1,2:] - ωn[1:-1,:-2])/2.0/dx
    Cy  =  (ωn[2:,1:-1] - ωn[:-2,1:-1])/2.0/dy * (ψ[1:-1,2:] - ψ[1:-1,:-2])/2.0/dx
    Dxy =  (ωn[1:-1,2:] - 2.0*ωn[1:-1,1:-1] + ωn[1:-1,:-2])/dx/dx + (ωn[2:,1:-1] -2.0*ωn[1:-1,1:-1] + ωn[:-2,1:-1])/dy/dy
    rhs = Cx + Cy + ν*Dxy
    ω[1:-1,1:-1] = ωn[1:-1,1:-1] + dt*rhs
    
    # apply boundary conditions
    ω[0,1:-1]  = -2.0*ψ[1, 1:-1]/dy/dy + 2.0*Ub/dy     # vorticity on bottom wall (stationary)
    ω[-1,1:-1] = -2.0*ψ[-2, 1:-1]/dy/dy - 2.0*Ut/dy # vorticity on top wall (moving at Uwall)
    ω[1:-1,-1] = -2.0*ψ[1:-1,-2]/dx/dx # right wall
    ω[1:-1,0]  = -2.0*ψ[1:-1,1]/dx/dx # left wall
    
    omega.append(ω)
    t += dt
[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
[[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00]
 [ 0.00000000e+00 -3.97433720e-05 -8.53052953e-05 -1.34838676e-04
  -1.86444808e-04 -2.38332260e-04 -2.88936181e-04 -3.36992444e-04
  -3.81568420e-04 -4.19374008e-04 -4.49070668e-04 -4.68926119e-04
  -4.76917208e-04 -4.71259249e-04 -4.50539329e-04 -4.13773754e-04
  -3.60547870e-04 -2.91144722e-04 -2.06605640e-04 -1.08727019e-04
   0.00000000e+00]
 [ 0.00000000e+00 -9.44345209e-05 -2.00688763e-04 -3.14579973e-04
  -4.31915598e-04 -5.48820714e-04 -6.61975715e-04 -7.68755020e-04
  -8.64027097e-04 -9.44185578e-04 -1.00546089e-03 -1.04391683e-03
  -1.05580875e-03 -1.03778777e-03 -9.87185859e-04 -9.02312584e-04
  -7.82714610e-04 -6.29381732e-04 -4.44859507e-04 -2.33236186e-04
   0.00000000e+00]
 [ 0.00000000e+00 -1.72535407e-04 -3.62847527e-04 -5.63679040e-04
  -7.67954669e-04 -9.69296649e-04 -1.16238047e-03 -1.33955121e-03
  -1.49548145e-03 -1.62364158e-03 -1.71801972e-03 -1.77283340e-03
  -1.78256978e-03 -1.74243816e-03 -1.64880857e-03 -1.49962680e-03
  -1.29481849e-03 -1.03660433e-03 -7.29654724e-04 -3.81044769e-04
   0.00000000e+00]
 [ 0.00000000e+00 -2.83530599e-04 -5.90301720e-04 -9.08985435e-04
  -1.22884715e-03 -1.54048040e-03 -1.83273593e-03 -2.09732159e-03
  -2.32507431e-03 -2.50749570e-03 -2.63642842e-03 -2.70410547e-03
  -2.70343221e-03 -2.62839635e-03 -2.47464731e-03 -2.24011355e-03
  -1.92556873e-03 -1.53506953e-03 -1.07616811e-03 -5.59824865e-04
   0.00000000e+00]
 [ 0.00000000e+00 -4.37933319e-04 -9.03394090e-04 -1.37969750e-03
  -1.85147748e-03 -2.30253451e-03 -2.71974937e-03 -3.09032221e-03
  -3.40284190e-03 -3.64638011e-03 -3.81076642e-03 -3.88651619e-03
  -3.86507849e-03 -3.73939565e-03 -3.50463304e-03 -3.15900502e-03
  -2.70457530e-03 -2.14788191e-03 -1.50024127e-03 -7.77616089e-04
   0.00000000e+00]
 [ 0.00000000e+00 -6.47287432e-04 -1.32428144e-03 -2.00728882e-03
  -2.67261570e-03 -3.29968104e-03 -3.87018884e-03 -4.36818263e-03
  -4.77945995e-03 -5.09115644e-03 -5.29137196e-03 -5.36913747e-03
  -5.31466478e-03 -5.11997565e-03 -4.77985504e-03 -4.29297657e-03
  -3.66303244e-03 -2.89965695e-03 -2.01890981e-03 -1.04312090e-03
   0.00000000e+00]
 [ 0.00000000e+00 -9.24261547e-04 -1.87708933e-03 -2.82386440e-03
  -3.73287022e-03 -4.57688713e-03 -5.33271415e-03 -5.98117653e-03
  -6.50598245e-03 -6.89279238e-03 -7.12856875e-03 -7.20125120e-03
  -7.09998452e-03 -6.81587309e-03 -6.34320723e-03 -5.68101589e-03
  -4.83470687e-03 -3.81748393e-03 -2.65117694e-03 -1.36613847e-03
   0.00000000e+00]
 [ 0.00000000e+00 -1.28296652e-03 -2.58736420e-03 -3.86538593e-03
  -5.07523251e-03 -6.18184334e-03 -7.15754020e-03 -7.98059963e-03
  -8.63369713e-03 -9.10212364e-03 -9.37253618e-03 -9.43231808e-03
  -9.26968444e-03 -8.87458301e-03 -8.24030815e-03 -7.36566212e-03
  -6.25737672e-03 -4.93235778e-03 -3.41918068e-03 -1.75822441e-03
   0.00000000e+00]
 [ 0.00000000e+00 -1.73908860e-03 -3.48476585e-03 -5.17242475e-03
  -6.74715158e-03 -8.16600283e-03 -9.39744641e-03 -1.04190086e-02
  -1.12142014e-02 -1.17696783e-02 -1.20730860e-02 -1.21119290e-02
  -1.18735236e-02 -1.13460625e-02 -1.05207352e-02 -9.39474909e-03
  -7.97492946e-03 -6.28133219e-03 -4.34999253e-03 -2.23372260e-03
   0.00000000e+00]
 [ 0.00000000e+00 -2.31210701e-03 -4.60568419e-03 -6.79301311e-03
  -8.80355042e-03 -1.05866209e-02 -1.21091033e-02 -1.33508466e-02
  -1.42993890e-02 -1.49453722e-02 -1.52792605e-02 -1.52895525e-02
  -1.49625293e-02 -1.42835093e-02 -1.32395882e-02 -1.18237960e-02
  -1.00404134e-02 -7.91079821e-03 -5.47845329e-03 -2.81142881e-03
   0.00000000e+00]
 [ 0.00000000e+00 -3.02823518e-03 -5.99740369e-03 -8.78794711e-03
  -1.13109984e-02 -1.35099734e-02 -1.53548071e-02 -1.68330632e-02
  -1.79412347e-02 -1.86777836e-02 -1.90384635e-02 -1.90138383e-02
  -1.85888057e-02 -1.77440278e-02 -1.64593272e-02 -1.47191983e-02
  -1.25204621e-02 -9.88154811e-03 -6.85176903e-03 -3.51737343e-03
   0.00000000e+00]
 [ 0.00000000e+00 -3.92497137e-03 -7.72590243e-03 -1.12387250e-02
  -1.43541931e-02 -1.70154974e-02 -1.92044875e-02 -2.09253771e-02
  -2.21909318e-02 -2.30126551e-02 -2.33947475e-02 -2.33310908e-02
  -2.28046194e-02 -2.17888351e-02 -2.02516736e-02 -1.81623163e-02
  -1.55016705e-02 -1.22767202e-02 -8.53720531e-03 -4.38969017e-03
   0.00000000e+00]
 [ 0.00000000e+00 -5.06012554e-03 -9.88942126e-03 -1.42603100e-02
  -1.80452208e-02 -2.12008595e-02 -2.37373443e-02 -2.56895809e-02
  -2.70968065e-02 -2.79901192e-02 -2.83858133e-02 -2.82824309e-02
  -2.76601991e-02 -2.64823932e-02 -2.46990506e-02 -2.22543363e-02
  -1.90996220e-02 -1.52144315e-02 -1.06354485e-02 -5.48761326e-03
   0.00000000e+00]
 [ 0.00000000e+00 -6.52982672e-03 -1.26427956e-02 -1.80216431e-02
  -2.25360840e-02 -2.61871453e-02 -2.90419628e-02 -3.11868179e-02
  -3.27002540e-02 -3.36401010e-02 -3.40382552e-02 -3.38991129e-02
  -3.31995019e-02 -3.18894485e-02 -2.98944141e-02 -2.71211301e-02
  -2.34711078e-02 -1.88679638e-02 -1.33047218e-02 -6.90923634e-03
   0.00000000e+00]
 [ 0.00000000e+00 -8.50663098e-03 -1.62428868e-02 -2.27777728e-02
  -2.80338740e-02 -3.21217406e-02 -3.52126078e-02 -3.74712434e-02
  -3.90287934e-02 -3.99754485e-02 -4.03607525e-02 -4.01955522e-02
  -3.94530686e-02 -3.80684131e-02 -3.59371999e-02 -3.29157570e-02
  -2.88289125e-02 -2.34975278e-02 -1.68059381e-02 -8.82943724e-03
   0.00000000e+00]
 [ 0.00000000e+00 -1.13284507e-02 -2.11356839e-02 -2.89171398e-02
  -3.48138202e-02 -3.91731164e-02 -4.23376277e-02 -4.45779616e-02
  -4.60854695e-02 -4.69827280e-02 -4.73353610e-02 -4.71600123e-02
  -4.64272191e-02 -4.50590068e-02 -4.29213676e-02 -3.98128525e-02
  -3.54544564e-02 -2.94971876e-02 -2.15895050e-02 -1.15886675e-02
   0.00000000e+00]
 [ 0.00000000e+00 -1.57291417e-02 -2.81229324e-02 -3.70184048e-02
  -4.32141276e-02 -4.75055592e-02 -5.04739665e-02 -5.25032198e-02
  -5.38346038e-02 -5.46115078e-02 -5.49079453e-02 -5.47437863e-02
  -5.40896024e-02 -5.28621693e-02 -5.09095686e-02 -4.79829773e-02
  -4.36922506e-02 -3.74513393e-02 -2.84616765e-02 -1.59233848e-02
   0.00000000e+00]
 [ 0.00000000e+00 -2.35052469e-02 -3.86540696e-02 -4.78689414e-02
  -5.35706584e-02 -5.72138743e-02 -5.96018146e-02 -6.11769161e-02
  -6.21858493e-02 -6.27644605e-02 -6.29801001e-02 -6.28512072e-02
  -6.23530766e-02 -6.14123184e-02 -5.98874931e-02 -5.75269848e-02
  -5.38842384e-02 -4.81528892e-02 -3.88762423e-02 -2.36326439e-02
   0.00000000e+00]
 [ 0.00000000e+00 -3.96603644e-02 -5.51420677e-02 -6.22552944e-02
  -6.60075714e-02 -6.81980828e-02 -6.95614346e-02 -7.04336784e-02
  -7.09820049e-02 -7.12924687e-02 -7.14062538e-02 -7.13346797e-02
  -7.10633779e-02 -7.05482255e-02 -6.97004458e-02 -6.83504411e-02
  -6.61601555e-02 -6.23934559e-02 -5.52501660e-02 -3.97223453e-02
   0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00]]
[[ 0.          0.          0.          0.          0.          0.
   0.          0.          0.          0.          0.          0.
   0.          0.          0.          0.          0.          0.
   0.          0.          0.        ]
 [ 0.         -0.0002371  -0.00047653 -0.00071144 -0.00093519 -0.00114168
  -0.00132557 -0.00148246 -0.00160834 -0.00169828 -0.00174937 -0.00175906
  -0.00172525 -0.00164671 -0.00152331 -0.00135615 -0.00114764 -0.00090172
  -0.00062376 -0.00032055  0.        ]
 [ 0.         -0.0004967  -0.00099786 -0.00148908 -0.00195658 -0.00238777
  -0.00277179 -0.00309908 -0.00335942 -0.00354507 -0.00364941 -0.00366696
  -0.0035937  -0.00342744 -0.0031681  -0.00281816 -0.00238291 -0.00187071
  -0.00129294 -0.00066384  0.        ]
 [ 0.         -0.0007932  -0.00159232 -0.00237435 -0.00311748 -0.00380211
  -0.00441084 -0.00492639 -0.0053352  -0.00562474 -0.00578475 -0.00580707
  -0.00568583 -0.00541792 -0.00500359 -0.00444705 -0.00375695 -0.00294675
  -0.00203475 -0.00104368  0.        ]
 [ 0.         -0.00114228 -0.00229059 -0.00341177 -0.00447486 -0.00545183
  -0.00631581 -0.00704505 -0.00761971 -0.00802318 -0.0082417  -0.0082644
  -0.0080836  -0.00769545 -0.00710068 -0.00630553 -0.00532255 -0.0041711
  -0.00287749 -0.00147445  0.        ]
 [ 0.         -0.00156114 -0.0031261  -0.00464957 -0.00608938 -0.00740571
  -0.00856509 -0.0095379  -0.01029941 -0.0108287  -0.0111087  -0.01112602
  -0.0108713  -0.01033992 -0.00953315 -0.00845942 -0.00713555 -0.00558768
  -0.0038515  -0.00197163  0.        ]
 [ 0.         -0.00206878 -0.00413563 -0.00613995 -0.00802435 -0.00973872
  -0.0112398  -0.01249146 -0.01346375 -0.01413214 -0.01447669 -0.01448178
  -0.0141364  -0.01343501 -0.01237916 -0.01097926 -0.00925659 -0.00724478
  -0.00499043 -0.00255249  0.        ]
 [ 0.         -0.00268658 -0.00535992 -0.00793893 -0.0103501  -0.01253015
  -0.01442632 -0.01599611 -0.01720551 -0.01802729 -0.01843951 -0.0184248
  -0.01797042 -0.01706948 -0.01572297 -0.01394252 -0.01175342 -0.00919719
  -0.00633298 -0.00323711  0.        ]
 [ 0.         -0.00343893 -0.00684444 -0.01011059 -0.01314407 -0.01586681
  -0.01821706 -0.02014729 -0.02162132 -0.02261112 -0.02309415 -0.02305193
  -0.02246986 -0.02133867 -0.01965719 -0.01743642 -0.01470416 -0.01150947
  -0.00792552 -0.00404975  0.        ]
 [ 0.         -0.00435442 -0.00864351 -0.01272957 -0.01649474 -0.01984561
  -0.02271269 -0.02504651 -0.02681214 -0.02798382 -0.02854062 -0.02846393
  -0.02773692 -0.02634616 -0.02428549 -0.02156172 -0.01820189 -0.01426064
  -0.00982602 -0.00502125  0.        ]
 [ 0.         -0.00546934 -0.01082501 -0.0158867  -0.02050717 -0.02457768
  -0.02802486 -0.03080278 -0.03288349 -0.03424831 -0.03488125 -0.03476497
  -0.03387975 -0.03220531 -0.02972579 -0.02643779 -0.02236115 -0.0175513
  -0.01211033 -0.00619259  0.        ]
 [ 0.         -0.0068326  -0.01347817 -0.01969827 -0.02531078 -0.03019436
  -0.0342796  -0.03753375 -0.03994491 -0.0415087  -0.04221871 -0.04206091
  -0.04101179 -0.03904029 -0.03611387 -0.0322091  -0.02732737 -0.02151466
  -0.01488228 -0.00762123  0.        ]
 [ 0.         -0.00851381 -0.01672702 -0.02432014 -0.03107115 -0.03685487
  -0.04162104 -0.04536621 -0.04810819 -0.04986702 -0.05065224 -0.05045576
  -0.04924927 -0.04698588 -0.04360665 -0.03905347 -0.03329032 -0.02633386
  -0.01829093 -0.00939205  0.        ]
 [ 0.         -0.01061846 -0.02075372 -0.0299704  -0.0380069  -0.04475593
  -0.05021474 -0.05443475 -0.05748309 -0.0594174  -0.06027097 -0.06004506
  -0.05870572 -0.05618456 -0.05238421 -0.04719149 -0.04050308 -0.03226967
  -0.02256066 -0.01163807  0.        ]
 [ 0.         -0.01331746 -0.02583992 -0.03696575 -0.04641312 -0.05414236
  -0.06024862 -0.06487698 -0.06816946 -0.07023659 -0.07114381 -0.07090546
  -0.06948197 -0.06677869 -0.06264746 -0.05689497 -0.04930723 -0.0397048
  -0.02804629 -0.01458221  0.        ]
 [ 0.         -0.01690804 -0.03244334 -0.04577925 -0.05669135 -0.06531489
  -0.07192789 -0.07682237 -0.08024398 -0.08237014 -0.083305   -0.08307936
  -0.08164974 -0.07889439 -0.07460647 -0.06848876 -0.06016185 -0.04921316
  -0.03533715 -0.01862807  0.        ]
 [ 0.         -0.02195117 -0.0413439  -0.0571284  -0.06938025 -0.07862639
  -0.08545677 -0.09037076 -0.09373988 -0.09581386 -0.09673598 -0.09655485
  -0.09522764 -0.09261272 -0.08845038 -0.08233169 -0.07366325 -0.0616595
  -0.04545833 -0.0245779   0.        ]
 [ 0.         -0.02961438 -0.05392615 -0.07209279 -0.08516379 -0.09444619
  -0.10099558 -0.10555605 -0.10861929 -0.11049208 -0.11134548 -0.11124301
  -0.11015011 -0.10792496 -0.10428636 -0.09874934 -0.0905144  -0.07830767
  -0.06025386 -0.03421189  0.        ]
 [ 0.         -0.04262301 -0.0727023  -0.09220587 -0.10479161 -0.11305571
  -0.11857951 -0.12229293 -0.12474034 -0.12623674 -0.12695283 -0.1269578
  -0.12623486 -0.12467418 -0.12203774 -0.11787542 -0.11134169 -0.10080146
  -0.08303083 -0.0520044   0.        ]
 [ 0.         -0.05787339 -0.09083797 -0.10741576 -0.11656265 -0.12201369
  -0.12544785 -0.12767942 -0.12912975 -0.13002557 -0.13048709 -0.13056777
  -0.13026981 -0.12954321 -0.12826565 -0.12618574 -0.12277052 -0.11675838
  -0.10468612 -0.07547754  0.        ]
 [ 0.          0.          0.          0.          0.          0.
   0.          0.          0.          0.          0.          0.
   0.          0.          0.          0.          0.          0.
   0.          0.          0.        ]]
[[ 0.          0.          0.          0.          0.          0.
   0.          0.          0.          0.          0.          0.
   0.          0.          0.          0.          0.          0.
   0.          0.          0.        ]
 [ 0.         -0.00040633 -0.00082319 -0.00122566 -0.00160256 -0.0019438
  -0.00224062 -0.00248565 -0.00267233 -0.00279523 -0.00285064 -0.00283615
  -0.00275071 -0.00259484 -0.00237083 -0.00208285 -0.00173704 -0.00134152
  -0.00090633 -0.00044318  0.        ]
 [ 0.         -0.00084584 -0.0017176  -0.00256035 -0.00335088 -0.00406819
  -0.00469399 -0.00521211 -0.00560833 -0.00587171 -0.00599394 -0.0059694
  -0.00579546 -0.00547277 -0.00500558 -0.00440216 -0.003675   -0.00284095
  -0.00192101 -0.00093987  0.        ]
 [ 0.         -0.00131976 -0.00268699 -0.00400825 -0.00524748 -0.00637213
  -0.00735311 -0.00816456 -0.0087855  -0.00919867 -0.00939106 -0.00935393
  -0.00908303 -0.0085791  -0.0078485  -0.00690378 -0.00576421 -0.00445601
  -0.00301213 -0.00147147  0.        ]
 [ 0.         -0.00184367 -0.00376512 -0.00562018 -0.00735856 -0.00893439
  -0.01030647 -0.01144034 -0.01230679 -0.01288251 -0.01314988 -0.01309693
  -0.01271771 -0.012013   -0.01099113 -0.00966906 -0.00807322 -0.00624001
  -0.00421564 -0.00205507  0.        ]
 [ 0.         -0.00243489 -0.00498908 -0.00745141 -0.00975479 -0.01183797
  -0.01364841 -0.01514116 -0.01627946 -0.01703398 -0.01738256 -0.01731004
  -0.01680865 -0.01587882 -0.01453056 -0.01278502 -0.010676   -0.00825093
  -0.00557094 -0.00270951  0.        ]
 [ 0.         -0.00311289 -0.00640014 -0.00956238 -0.01251239 -0.01517329
  -0.01747901 -0.01937483 -0.02081648 -0.02176921 -0.02220673 -0.0221108
  -0.02147148 -0.02028825 -0.01857193 -0.01634714 -0.01365478 -0.01055399
  -0.00712275 -0.00345642  0.        ]
 [ 0.         -0.00390002 -0.00804485 -0.01202031 -0.01571671 -0.01903874
  -0.02190676 -0.02425657 -0.02603761 -0.02721078 -0.02774656 -0.02762387
  -0.02683013 -0.02536265 -0.02323143 -0.02046282 -0.01710372 -0.0132251
  -0.00892385 -0.00432156  0.        ]
 [ 0.         -0.00482246 -0.00997699 -0.01490314 -0.01946414 -0.02354426
  -0.02705052 -0.02991084 -0.03207058 -0.03348855 -0.03413354 -0.03398213
  -0.03301838 -0.03123551 -0.02863974 -0.02525585 -0.02113393 -0.01635574
  -0.01103897 -0.00533696  0.        ]
 [ 0.         -0.00591207 -0.01226192 -0.01830354 -0.02386715 -0.02881519
  -0.03304267 -0.03647314 -0.03905205 -0.04073993 -0.0415065  -0.04132701
  -0.04018116 -0.03805503 -0.03494626 -0.03087221 -0.02588051 -0.02006023
  -0.01355074 -0.00654415  0.        ]
 [ 0.         -0.00720981 -0.01498222 -0.02233591 -0.02906118 -0.03499779
  -0.04003282 -0.044092   -0.04712802 -0.04910926 -0.05001068 -0.04980786
  -0.04847489 -0.04598634 -0.04232388 -0.0374877  -0.03151257 -0.0244867
  -0.01656927 -0.0079995   0.        ]
 [ 0.         -0.00877093 -0.01824691 -0.02714757 -0.03521466 -0.04226701
  -0.04819251 -0.05293074 -0.05645322 -0.05874572 -0.05979484 -0.05957927
  -0.05806618 -0.05521279 -0.05097411 -0.04531785 -0.03824773 -0.0298342
  -0.02024783 -0.00978339  0.        ]
 [ 0.         -0.01067372 -0.02220717 -0.03293606 -0.04254401 -0.05083677
  -0.05772043 -0.06316833 -0.06718871 -0.06979865 -0.07100553 -0.07079543
  -0.06912764 -0.06593502 -0.0611316  -0.05463007 -0.04637296 -0.03638019
  -0.02481014 -0.01201698  0.        ]
 [ 0.         -0.01303565 -0.02708334 -0.03997705 -0.05133558 -0.06097333
  -0.06884768 -0.07499828 -0.07949644 -0.08240941 -0.08377742 -0.08360005
  -0.08182887 -0.07836539 -0.07306581 -0.06575723 -0.05627392 -0.04452482
  -0.03059918 -0.01489455  0.        ]
 [ 0.         -0.01604423 -0.03321288 -0.04866961 -0.06197717 -0.07301117
  -0.08184109 -0.0886236  -0.09352913 -0.09669832 -0.0982184  -0.09811023
  -0.09632021 -0.09271432 -0.0870752  -0.07910861 -0.06847358 -0.05486282
  -0.03816885 -0.01875127  0.        ]
 [ 0.         -0.02002029 -0.04113934 -0.0596099  -0.07500098 -0.08736834
  -0.09700033 -0.10424424 -0.10941334 -0.11274571 -0.11438877 -0.1143932
  -0.11270663 -0.10916282 -0.1034649  -0.09516748 -0.08367432 -0.06829621
  -0.04846187 -0.02421827  0.        ]
 [ 0.         -0.02555839 -0.05178321 -0.07370856 -0.09113469 -0.10455162
  -0.11464031 -0.12203172 -0.12722343 -0.13056654 -0.13227519 -0.13243633
  -0.13101013 -0.1278151  -0.12249318 -0.1144504  -0.10277775 -0.08619339
  -0.06316219 -0.0326159   0.        ]
 [ 0.         -0.03386828 -0.06677832 -0.09236474 -0.11134037 -0.12512853
  -0.13504369 -0.14208404 -0.14694464 -0.15008038 -0.15176244 -0.152114
  -0.15112276 -0.14862723 -0.14426765 -0.13738092 -0.12680313 -0.11054052
  -0.08537411 -0.04707103  0.        ]
 [ 0.         -0.04647295 -0.0877149  -0.11606991 -0.13503876 -0.14778953
  -0.15646573 -0.1624101  -0.16644538 -0.16906915 -0.17057051 -0.17109347
  -0.17066402 -0.16918606 -0.16639777 -0.16175522 -0.1541533  -0.14125565
  -0.11791193 -0.07289725  0.        ]
 [ 0.         -0.05503041 -0.09883152 -0.12366617 -0.13776185 -0.14629466
  -0.1517309  -0.15531438 -0.157708   -0.1592803  -0.16023878 -0.16069357
  -0.16068589 -0.16019419 -0.15911664 -0.15720828 -0.15389204 -0.14764032
  -0.13368427 -0.09551412  0.        ]
 [ 0.          0.          0.          0.          0.          0.
   0.          0.          0.          0.          0.          0.
   0.          0.          0.          0.          0.          0.
   0.          0.          0.        ]]
In [15]:
solution = psi[-1]
plt.contour(solution,levels=20)
plt.colorbar()
Out[15]:
<matplotlib.colorbar.Colorbar at 0x11135df98>

Calculate the velocity vector field

The next step is to calculate the velocity vector field and double check that we are divergence free

In [ ]:
ψ = psi[-1]
u = np.zeros([ny,nx])
v = np.zeros([ny,nx])

u[1:-1,1:-1] = (ψ[2:,1:-1] - ψ[:-2,1:-1])/2.0/dy
v[1:-1,1:-1] = -(ψ[1:-1, 2:] - ψ[1:-1,:-2])/2.0/dx

# now will in the boundary values - they are all zero except at the top boundary where u = Ut
u[-1,1:-1] = Ut
divu = (u[1:-1,2:] - u[1:-1,:-2])/2.0/dx + (v[2:,1:-1] - v[:-2,1:-1])/2.0/dy
print(np.linalg.norm(divu.ravel()))

# note the size of the velocity vector field - it does not have points at the boundaries.

Create some beautiful plots

In [ ]:
x = np.linspace(0,1,nx)
y = np.linspace(0,1,ny)
xx,yy = np.meshgrid(x,y)
nn = 1
plt.quiver(xx[::nn,::nn],yy[::nn,::nn],u[::nn,::nn],v[::nn,::nn])
# plt.quiver(xx[1:-1:nn,1:-1:nn],yy[1:-1:nn,1:-1:nn],u[::nn,::nn],v[::nn,::nn])
# plt.axis('equal')
In [ ]:
fig = plt.figure(figsize=[6,6],dpi=600)
ax = fig.add_subplot(111)
ax.contourf(xx, yy, np.sqrt(u*u + v*v), levels = 100, cmap=plt.cm.jet)
ax.streamplot(xx,yy,u, v, color=np.sqrt(u*u + v*v),density=1.1,cmap=plt.cm.autumn,linewidth=1.5)
# ax.set_xlim([xx[0,0],xx[0,-1]])
# ax.set_ylim([0,1])
# ax.set_aspect(1)
In [ ]:
x = np.linspace(0,1,nx)
y = np.linspace(0,1,ny)
xx,yy = np.meshgrid(x,y)
nn = 3

fig = plt.figure(figsize=(6.1,5),facecolor='w',dpi=150)
ax = fig.add_subplot(111)


ims = []
# levs = np.linspace(-1,1,20)
i = 0
t = 0.0
for sol in psi:    
    if (i%10==0):      
        u = (sol[2:,1:-1] - sol[:-2,1:-1])/2.0/dy
        v = -(sol[1:-1, 2:] - sol[1:-1,:-2])/2.0/dx
        im = ax.contourf(xx[1:-1,1:-1], yy[1:-1,1:-1], np.sqrt(u*u + v*v), levels = 100, cmap=plt.cm.jet)
#         im = ax.streamplot(xx[1:-1,1:-1],yy[1:-1,1:-1],u,v, color=abs(u*u + v*v),cmap=plt.cm.autumn,linewidth=1.0)
#         im = plt.contourf(solution,cmap=cm.jet,levels=100)
        ims.append(im.collections)
    i+=1
    t += dt

ax.set_xlim([0,1])
ax.set_ylim([0,1])
ax.set_aspect(1)

# cbar = plt.colorbar()
# plt.clim(-10,10)
# cbar.set_ticks(np.linspace(0,10,5))

ani = animation.ArtistAnimation(fig, ims, interval=35, blit=True,repeat_delay=1000)

ani
In [ ]:
divu = (u[1:-1, 2:] - u[1:-1,:-2])/2.0/dx + (v[2:,1:-1] - v[:-2,1:-1])/2.0/dy
In [ ]:
plt.imshow(divu)
plt.colorbar()
In [ ]:
import urllib
import requests
from IPython.core.display import HTML
def css_styling():
    styles = requests.get("https://raw.githubusercontent.com/saadtony/NumericalMethods/master/styles/custom.css")
    return HTML(styles.text)
css_styling()