#!/usr/bin/env python
# coding: utf-8
# In[117]:
from IPython.display import HTML
s="""
2D Navier-Stokes for Flow in a Cavity (vs Ghia)
""";
h=HTML(s); h
# ***
# # Understand the Problem
#
# ## Question
#
# * What is the 2D velocity **and** pressure field for the **2D Navier Stokes Equation** for flow over a cavity?
#
# ## Initial Conditions
#
# * The velocity and pressure fields are zero everywhere
#
# ## Boundary Conditions
#
# Dimensionless form:
#
# * The u-boundary condition at y=1 is 1 (the lid)
#
# * The other velocity boundary conditions are zero (no slip).
#
# * The p-boundary condition at y=1 is 0 (atmospheric pressure)
#
# * The gradient of pressure at y=0, x=0 and x=1 is zero (Neumann)
#
# * Density = 1
#
# * Viscosity = 1/Re = 1/100
#
# ## Governing Equations
#
# * The Navier Stokes Momentum Equation is described as follows:
#
# $$ {\partial u \over \partial t} + u {\partial u \over \partial x} + v {\partial u \over \partial y} = -{1 \over \rho} {{\partial p} \over {\partial x}} + \nu \left ( {\partial^2 u \over \partial x^2}+ {\partial^2 u \over \partial y^2} \right ) $$
#
# $$ {\partial v \over \partial t} + u {\partial v \over \partial x} + v {\partial v \over \partial y} = -{1 \over \rho} {{\partial p} \over {\partial y}} + \nu \left ( {\partial^2 v \over \partial x^2}+ {\partial^2 v \over \partial y^2} \right ) $$
#
# * The Poisson Equation to link pressure and velocity is:
#
# $$ \nabla^2 p^{n+1} = \rho {{\nabla \cdot \mathbf{u}^n} \over {\Delta t}}-
# \rho \nabla \cdot (\mathbf{u}^n \cdot \nabla \mathbf{u}^n)+
# \mu \nabla^2 (\nabla \cdot \mathbf{u}^n)
# $$
#
# * **Assuming viscosity is small**, the Poisson Equation to link pressure and velocity is as follows:
#
# $$ {{{{\partial^2 p} \over {\partial x^2}} + {{\partial^2 p} \over {\partial y^2}} =
# {\rho \over \Delta t} \left( {{\partial u} \over {\partial x}} + {{\partial v} \over {\partial y}} \right) - \rho \left[ \left( {{\partial u} \over {\partial x}} \right)^2 +
# 2 {{\partial u} \over {\partial y}} {{\partial v} \over {\partial x}} +
# \left( {{\partial v} \over {\partial y}} \right)^2 \right]}} $$
#
# (note: this equation applies in the discrete domain)
#
# ***
#
# # Formulate the Problem
#
# ## Input Data:
#
# The Poisson Equation has **no temporal component**, so we use a number of iterations `niter`
#
# The Navier Stokes Momentum Equation **does have a temporal component**, so we use `nt`
#
# * `niter` = 51 (maximum number of iterations - for Poisson Equation)
# * `nt` = 51 (number of temporal points)
# * `nx` = 129 (number of x spatial points)
# * `ny` = 129 (number of y spatial points)
# * `tmax` = 6
# * `xmax` = 1
# * `ymax` = 1
# * `nu` = 1/100
# * `rho` = 1
#
# ## Initial Conditions:
#
# * $\forall (n, x, y) \quad n = 0 \rightarrow u = 0 \land v = 0 \land p = 0$
#
# ## Velocity Boundary Conditions:
#
# * $\forall (n, x, y) \quad y = 1 \rightarrow u = 1$
#
# * $\forall (n, x, y) \quad x = 0 \lor x = 1 \lor y = 0 \rightarrow u = 0$
#
# * $\forall (n, x, y) \quad x = 0 \lor x = 1 \lor y = 0 \lor y = 1 \rightarrow v = 0$
#
# ## Pressure Boundary Conditions:
#
# * $\forall (n, x, y) \quad y = 1 \rightarrow p = 0$
#
# * $\forall (n, x, y) \quad y = 0 \rightarrow {{\partial p} \over {\partial y}} = 0$
#
# * $\forall (n, x, y) \quad x = 0 \lor x = 1 \rightarrow {{\partial p} \over {\partial x}} = 0$
#
# ## Output Data:
#
# * $\forall (n,x,y) \quad \ p = ? \land u = ? \land v = ?$
#
# ***
#
# # Design Algorithm to Solve Problem
#
# ## Space-time discretisation:
#
# * i $\rightarrow$ index of grid in x
# * j $\rightarrow$ index of grid in y
# * n $\rightarrow$ index of time
# * m $\rightarrow$ index of iterations
#
# ## Numerical schemes for each Momentum Equation
#
# * For the **one** first derivative of velocity in time: 1st order FD in time
# * For the **two** first derivatives of velocity in space: 1st order BD in space
# * For the **one** first derivative of pressure in space: 2nd order CD in space
# * For the **two** second derivatives of velocity in space: 2nd order CD in space
#
# ## Numerical schemes for the Poisson Equation
#
# * For the **two** second derivatives of pressure in space: 2nd order CD in space
# * The the **four** first derivatives of velocity in space: 2nd order CD in space
#
# ## Discrete equation for u-Momentum Equation
#
# $$ {{u_{i,j}^{n+1} - u_{i,j}^n} \over {\Delta t}} +
# u_{i,j}^n {{u_{i,j}^n - u_{i-1,j}^n} \over \Delta x} +
# v_{i,j}^n {{u_{i,j}^n - u_{i,j-1}^n} \over \Delta y} = \\
# -{1 \over \rho} {{p_{i+1,j}^n - p_{i-1,j}^n} \over {2 \Delta x}} +
# \nu {{u_{i-1,j}^n - 2u_{i,j}^n + u_{i+1,j}^n} \over \Delta x^2} +
# \nu {{u_{i,j-1}^n - 2u_{i,j}^n + u_{i,j+1}^n} \over \Delta y^2}
# $$
#
# ## Transpose
#
# Assume $ \Delta x = \Delta y = h$
#
# $$ u_{i,j}^{n+1} = u_{i,j}^n - {{\Delta t} \over h} \left[ u_{i,j}^n(u_{i,j}^n - u_{i-1,j}^n) + v_{i,j}^n(u_{i,j}^n - u_{i,j-1}^n)\right] - \\
# {{\Delta t} \over {2 \rho h}} (p_{i+1,j}^n - p_{i-1,j}^n) + {{\Delta t \nu} \over {h^2}}
# (u_{i-1,j}^n + u_{i+1,j}^n + u_{i,j-1}^n + u_{i,j+1}^n - 4 u_{i,j}^n )
# $$
#
# ## Discrete equation for v-Momentum Equation
#
# $$ {{v_{i,j}^{n+1} - v_{i,j}^n} \over {\Delta t}} +
# u_{i,j}^n {{v_{i,j}^n - v_{i-1,j}^n} \over \Delta x} +
# v_{i,j}^n {{v_{i,j}^n - v_{i,j-1}^n} \over \Delta y} = \\
# -{1 \over \rho} {{p_{i,j+1}^n - p_{i,j-1}^n} \over {2 \Delta y}} +
# \nu {{v_{i-1,j}^n - 2v_{i,j}^n + v_{i+1,j}^n} \over \Delta x^2} +
# \nu {{v_{i,j-1}^n - 2v_{i,j}^n + v_{i,j+1}^n} \over \Delta y^2}
# $$
#
# ## Transpose
#
# Assume $ \Delta x = \Delta y = h$
#
# $$ v_{i,j}^{n+1} = v_{i,j}^n - {{\Delta t} \over h} \left[ u_{i,j}^n(v_{i,j}^n - v_{i-1,j}^n) + v_{i,j}^n(v_{i,j}^n - v_{i,j-1}^n)\right] - \\
# {{\Delta t} \over {2 \rho h}} (p_{i,j+1}^n - p_{i,j-1}^n) + {{\Delta t \nu} \over {h^2}}
# (v_{i-1,j}^n + v_{i+1,j}^n + v_{i,j-1}^n + v_{i,j+1}^n - 4 v_{i,j}^n )
# $$
#
# ## Discrete equation for Poisson Equation must be Divergence Free
#
# Since no terms have a differential temporal component - bring source term to other side and equate to forward differencing (as we did for Laplace Equation). Then take steady state after a certain number of iterations.
#
# $$ {{{{\partial p} \over {\partial t}} = {\partial^2 p \over \partial x^2} + {\partial^2 p \over \partial y^2} -
# \left[ {\rho \over \Delta t} \left( {{\partial u} \over {\partial x}} + {{\partial v} \over {\partial y}} \right) - \rho \left[ \left( {{\partial u} \over {\partial x}} \right)^2 +
# 2 {{\partial u} \over {\partial y}} {{\partial v} \over {\partial x}} +
# \left( {{\partial v} \over {\partial y}} \right)^2 \right] \right]}} $$
#
# $$ {{\partial p} \over {\partial t}} = {\partial^2 p \over \partial x^2} + {\partial^2 p \over \partial y^2} -
# b $$
#
# $$ {{p_{i,j}^{m+1}-p_{i,j}^m} \over {\Delta \tau}} = {{p_{i+1,j}^m -2p_{i,j}^m + p_{i-1,j}^m} \over \Delta x^2} + {{p_{i,j+1}^m -2p_{i,j}^m + p_{i,j-1}^m} \over \Delta y^2} - b_{i,j}^m $$
#
# Assume that $ \Delta x = \Delta y = h $
#
# $$ {{p_{i,j}^{m+1}-p_{i,j}^m} \over {\Delta \tau}} = {{p_{i+1,j}^m -2p_{i,j}^m + p_{i-1,j}^m} \over h^2} + {{p_{i,j+1}^m -2p_{i,j}^m + p_{i,j-1}^m} \over h^2} - b_{i,j}^m $$
#
#
# $$ p_{i,j}^{m+1} = p_{i,j}^m + {{\Delta \tau} \over {h^2}} \left(
# p_{i+1,j}^m + p_{i-1,j}^m + p_{i,j+1}^m + p_{i,j-1}^m - 4p_{i,j}^m \right) - b_{i,j}^m \Delta \tau $$
#
# For the fastest convergence $ r = {{\Delta \tau} \over {h^2}} = {1 \over 4} $ and $ \Delta \tau = {{h^2} \over 4} $
#
# Hence:
#
# $$ p_{i,j}^{m+1} = {1 \over 4} \left( p_{i+1,j}^m + p_{i-1,j}^m + p_{i,j+1}^m + p_{i,j-1}^m - b_{i,j}^m h^2 \right) $$
#
# So now we need to define $ b_{i,j}^m $
#
# $$
# b = {{\rho \over \Delta t} \left( {{\partial u} \over {\partial x}} + {{\partial v} \over {\partial y}} \right) - \rho \left[ \left( {{\partial u} \over {\partial x}} \right)^2 +
# 2 {{\partial u} \over {\partial y}} {{\partial v} \over {\partial x}} +
# \left( {{\partial v} \over {\partial y}} \right)^2 \right] }
# $$
#
#
#
# $$ b_{i,j}^{m} =
# {\rho \over \Delta t} \left( {{u_{i+1,j}^n - u_{i-1,j}^n} \over {2 \Delta x}} + {{v_{i,j+1}^n - v_{i,j-1}^n} \over {2 \Delta y}} \right) -
# \rho \left[ \left( {{u_{i+1,j}^n - u_{i-1,j}^n} \over {2 \Delta x}} \right)^2 + 2 \left( {{u_{i,j+1}^n - u_{i,j-1}^n} \over {2 \Delta y}} {{v_{i+1,j}^n - v_{i-1,j}^n} \over {2 \Delta x}}\right) +
# \left( {{v_{i,j+1}^n - v_{i,j-1}^n} \over {2 \Delta y}} \right)^2 \right] $$
#
# Also assume that $ \Delta x = \Delta y = h $
#
#
# $$ b_{i,j}^{m} =
# {\rho \over {2 h \Delta t}} ( {{u_{i+1,j}^n - u_{i-1,j}^n}} + {{v_{i,j+1}^n - v_{i,j-1}^n}} ) -
# {\rho \over {4h^2}} \left[ (u_{i+1,j}^n - u_{i-1,j}^n)^2 + 2 ( u_{i,j+1}^n - u_{i,j-1}^n) (v_{i+1,j}^n - v_{i-1,j}^n) + (v_{i,j+1}^n - v_{i,j-1}^n)^2 \right] $$
# ## Navier Stokes Equations - slice notation
# In[3]:
def navier_stokes_initialisation(niter, r, nx_or_ny, tmax, xmax_or_ymax):
"""
Returns the velocity field and distance for 2D linear convection
"""
# Increments:
nx = ny = nx_or_ny
xmax = ymax = xmax_or_ymax
dx = xmax/(nx-1)
dy = ymax/(ny-1)
nt = int((tmax / (r*(dx)**2))+1)
dt = tmax/(nt-1)
# Initialise data structures:
import numpy as np
p = np.zeros((nx,ny))
u = np.zeros((nx,ny))
v = np.zeros((nx,ny))
# linspace is SIMPLER than list comprehensions:
x = np.linspace(0.0,1.0,nx)
y = np.linspace(0.0,1.0,ny)
# Pressure Boundary Conditions:
p[:, ny-1] = 0.0
# Velocity Boundary Conditions:
u[:,ny-1] = 1.0
return p, x, y, u, v, nx, ny, nt, dx, dy, dt, niter, r
# In[74]:
def navier_stokes(rho, nu, niter, r, nx, tmax, xmax):
(p, x, y, u, v, nx, ny, nt,
dx, dy, dt, niter, r) = navier_stokes_initialisation(niter, r, nx, tmax, xmax)
# Increments
h = dx
error_target = 1.0e-2
import numpy as np
# Intermediate copies:
un = np.zeros((nx, ny))
vn = np.zeros((nx, ny))
pm = np.zeros((nx, ny))
bm = np.zeros((nx, ny)) # bm needs to be exactly zero at the boundaries
# Loop - use decimal points for all floating point numbers
for n in range(nt):
# We know the velocity at i=0, j=0, i=nx-1 and j=ny-1. b is zero at the boundaries.
bm[1:-1, 1:-1] = ( (rho / (2.0 * h * dt)) * ( u[2:, 1:-1] - u[0:-2, 1:-1]
+ v[1:-1, 2:] - v[1:-1, 0:-2] ) -
(rho / (4.0*h**2)) * ( (u[2:, 1:-1] - u[0:-2, 1:-1])**2.0 +
2.0*(u[1:-1, 2:] - u[1:-1, 0:-2])*(v[2:, 1:-1] - v[0:-2, 1:-1]) +
(v[1:-1, 2:] - v[1:-1, 0:-2])**2.0 ) )
# First points for p. We don't know the pressure at i=0, j=0 and i=nx-1. We DO know the pressure at j=ny-1
while True:
for m in range(niter):
pm = np.copy(p)
p[1:-1, 1:-1] = 0.25*( pm[2:, 1:-1] + pm[0:-2, 1:-1] + pm[1:-1, 2:] + pm[1:-1, 0:-2]
- bm[1:-1, 1:-1]*h**2.0 )
# Set zero gradient boundary conditions:
p[0, :] = p[1, :]
p[:, 0] = p[:, 1]
p[-1, :] = p[-2, :]
error = np.abs(np.sum(np.abs(p[1:-1, 1:-1])-np.abs(pm[1:-1, 1:-1])))
if(error < error_target):
#print "n = " + str(m) + " completed"
break
break
# First points for u and v. We know the velocity at i=0, j=0, i=nx-1 and j=ny-1.
# We are simply using the value of pressure here
un = np.copy(u)
vn = np.copy(v)
u[1:-1, 1:-1] = ( un[1:-1, 1:-1] -
(dt / h) * ( un[1:-1, 1:-1] * ( un[1:-1, 1:-1] - un[0:-2, 1:-1] ) +
vn[1:-1, 1:-1] * ( un[1:-1, 1:-1] - un[1:-1, 0:-2] ) ) -
(dt / (2.0 * rho * h)) * ( p[2:, 1:-1] - p[0:-2, 1:-1] ) +
(dt * nu / h**2.0) * ( un[0:-2, 1:-1] + un[2:, 1:-1] + un[1:-1, 0:-2] + un[1:-1, 2:] -
4.0 * un[1:-1, 1:-1] ) )
v[1:-1, 1:-1] = ( vn[1:-1, 1:-1] -
(dt / h) * ( un[1:-1, 1:-1] * ( vn[1:-1, 1:-1] - vn[0:-2, 1:-1] ) +
vn[1:-1, 1:-1] * ( vn[1:-1, 1:-1] - vn[1:-1, 0:-2] ) ) -
(dt / (2.0 * rho * h)) * ( p[1:-1, 2:] - p[1:-1, 0:-2] ) +
(dt * nu / h**2.0) * ( vn[0:-2, 1:-1] + vn[2:, 1:-1] + vn[1:-1, 0:-2] + vn[1:-1, 2:] -
4.0 * vn[1:-1, 1:-1] ) )
return u, v, p, x, y
# In[75]:
get_ipython().run_line_magic('timeit', 'navier_stokes(1, 0.01, 50, 0.5, 51, 0.5, 1.0)')
# In[76]:
u00, v00, p00, x00, y00 = navier_stokes(1, 0.01, 50, 0.5, 51, 0.5, 1.0)
# In[77]:
u01, v01, p01, x01, y01 = navier_stokes(1, 0.01, 50, 0.5, 51, 1.0, 1.0)
# In[78]:
u02, v02, p02, x02, y02 = navier_stokes(1, 0.01, 50, 0.5, 51, 1.5, 1.0)
# In[79]:
u03, v03, p03, x03, y03 = navier_stokes(1, 0.01, 50, 0.5, 51, 2.0, 1.0)
# In[80]:
u04, v04, p04, x04, y04 = navier_stokes(1, 0.01, 50, 0.5, 51, 2.5, 1.0)
# In[81]:
u05, v05, p05, x05, y05 = navier_stokes(1, 0.01, 50, 0.5, 51, 3.0, 1.0)
# In[82]:
u06, v06, p06, x06, y06 = navier_stokes(1, 0.01, 50, 0.5, 51, 3.5, 1.0)
# In[83]:
u07, v07, p07, x07, y07 = navier_stokes(1, 0.01, 50, 0.5, 51, 4.0, 1.0)
# In[84]:
u08, v08, p08, x08, y08 = navier_stokes(1, 0.01, 50, 0.5, 51, 4.5, 1.0)
# In[85]:
u09, v09, p09, x09, y09 = navier_stokes(1, 0.01, 50, 0.5, 51, 5.0, 1.0)
# In[86]:
u010, v010, p010, x010, y010 = navier_stokes(1, 0.01, 50, 0.5, 51, 5.5, 1.0)
# In[87]:
u011, v011, p011, x011, y011 = navier_stokes(1, 0.01, 50, 0.5, 51, 6.0, 1.0)
# In[88]:
u012, v012, p012, x012, y012 = navier_stokes(1, 0.01, 50, 0.5, 51, 6.5, 1.0)
# In[89]:
u013, v013, p013, x013, y013 = navier_stokes(1, 0.01, 50, 0.5, 51, 7.0, 1.0)
# In[90]:
u014, v014, p014, x014, y014 = navier_stokes(1, 0.01, 50, 0.5, 51, 7.5, 1.0)
# In[91]:
def vector_contour(u, v, p, x, y):
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
fig = plt.figure(figsize=(11,7), dpi=100)
Y,X=np.meshgrid(y,x) #note meshgrid uses y,x not x,y!!!
plt.contourf(X,Y,p[:,:],alpha=0.5) ###plotting the pressure field as a contour
plt.colorbar()
plt.contour(X,Y,p[:,:]) ###plotting the pressure field outlines
plt.quiver(X[::2,::2],Y[::2,::2],u[::2,::2],v[::2,::2]) ##plotting velocity
plt.xlabel('X')
plt.ylabel('Y')
# # Results
# ## 51 spatial intervals tmax = 0.5 seconds
# In[92]:
get_ipython().run_line_magic('matplotlib', 'inline')
vector_contour(u00, v00, p00, x00, y00)
# ## 51 spatial intervals tmax = 7.5 seconds
# In[98]:
get_ipython().run_line_magic('matplotlib', 'inline')
vector_contour(u014, v014, p014, x014, y014)
# ## How much time needs to pass before the open cavity flow settles down to a steady state?
# In[93]:
def plot_u_vs_y_2(u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, y0):
"""
Plots the 1D velocity field
"""
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
u_vel = np.zeros((51,15))
u_vel[:, 0] = u0[25, :]
u_vel[:, 1] = u1[25, :]
u_vel[:, 2] = u2[25, :]
u_vel[:, 3] = u3[25, :]
u_vel[:, 4] = u4[25, :]
u_vel[:, 5] = u5[25, :]
u_vel[:, 6] = u6[25, :]
u_vel[:, 7] = u7[25, :]
u_vel[:, 8] = u8[25, :]
u_vel[:, 9] = u9[25, :]
u_vel[:, 10] = u10[25, :]
u_vel[:, 11] = u11[25, :]
u_vel[:, 12] = u12[25, :]
u_vel[:, 13] = u13[25, :]
u_vel[:, 14] = u14[25, :]
ghia_y = [0.0000, 0.0547, 0.0625, 0.0703, 0.1016, 0.1716, 0.2813, 0.4531, 0.5000,
0.6172, 0.7344, 0.8516, 0.9531, 0.9609, 0.9688, 0.9766, 1.0000]
ghia_u = [0.0000, -0.03717, -0.04192, -0.04775, -0.06434, -0.10150, -0.15662, -0.21090,
-0.20581, -0.13641, 0.00332, 0.23151, 0.68717, 0.73722, 0.78871, 0.84123, 1.0000]
plt.figure()
ax=plt.subplot(111)
colour=iter(cm.rainbow(np.linspace(0, 1, 15)))
for n in range(15):
c=next(colour)
ax.plot(y0,u_vel[:, n],linestyle='-',c=c, label='t= '+str(0.5*n+0.5)+ ' sec')
ax.plot(ghia_y, ghia_u, linestyle='none', marker='x', c='k', label='ghia')
box=ax.get_position()
ax.set_position([box.x0, box.y0, box.width*2,box.height*2])
ax.legend( bbox_to_anchor=(1.0,1.0), loc=2)
plt.xlabel('y/L (-)')
plt.ylabel('u/U (-)')
plt.show()
# In[94]:
plot_u_vs_y_2(u00, u01, u02, u03, u04, u05, u06, u07, u08, u09, u010, u011, u012, u013, u014, y00)
# In[95]:
def plot_v_vs_x(v0, v1, v2, v3, v4, v5, v6, v7, v8, v9, v10, v11, v12, v13, v14, x0):
"""
Plots the 1D velocity field
"""
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
v_vel = np.zeros((51,15))
v_vel[:, 0] = v0[:, 25]
v_vel[:, 1] = v1[:, 25]
v_vel[:, 2] = v2[:, 25]
v_vel[:, 3] = v3[:, 25]
v_vel[:, 4] = v4[:, 25]
v_vel[:, 5] = v5[:, 25]
v_vel[:, 6] = v6[:, 25]
v_vel[:, 7] = v7[:, 25]
v_vel[:, 8] = v8[:, 25]
v_vel[:, 9] = v9[:, 25]
v_vel[:, 10] = v10[:, 25]
v_vel[:, 11] = v11[:, 25]
v_vel[:, 12] = v12[:, 25]
v_vel[:, 13] = v13[:, 25]
v_vel[:, 14] = v14[:, 25]
ghia_x = [0.0000, 0.0625, 0.0703, 0.0781, 0.0938, 0.1563, 0.2266, 0.2344, 0.5000,
0.8047, 0.8594, 0.9063, 0.9453, 0.9531, 0.9609, 0.9688, 1.0000]
ghia_v = [0.0000, 0.09233, 0.10091, 0.10890, 0.12317, 0.16077, 0.17507, 0.17527, 0.05454,
-0.24533, -0.22445, -0.16914, -0.10313, -0.08864, -0.07391, -0.05906, 0.0000]
plt.figure()
ax=plt.subplot(111)
colour=iter(cm.rainbow(np.linspace(0, 1, 15)))
for n in range(15):
c=next(colour)
ax.plot(x0,v_vel[:, n],linestyle='-',c=c, label='t= '+str(0.5*n+0.5)+ ' sec')
ax.plot(ghia_x, ghia_v, linestyle='none', c='k', marker='x', label='ghia')
box=ax.get_position()
ax.set_position([box.x0, box.y0, box.width*2,box.height*2])
ax.legend( bbox_to_anchor=(1.0,1.0), loc=2)
plt.xlabel('x/L (-)')
plt.ylabel('v/U (-)')
plt.show()
# In[96]:
plot_v_vs_x(v00, v01, v02, v03, v04, v05, v06, v07, v08, v09, v010, v011, v012, v013, v014, x00)
# Around 7.5 seconds
#