In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from sympy.abc import x, y
from sympy import sin, atan2, sqrt, lambdify, latex
from IPython.display import display
from sympy.interactive import printing
printing.init_printing()
$\psi = U\left( r - \dfrac{R^2}{r} \right) \sin{\theta}$
In [ ]:
def cylinder_stream_function(U=1, R=1):
r = sqrt(x**2 + y**2)
theta = atan2(y, x)
return U * (r - R**2 / r) * sin(theta)
In [ ]:
def velocity_field(psi):
u = lambdify((x, y), psi.diff(y), 'numpy')
v = lambdify((x, y), -psi.diff(x), 'numpy')
ut = lambdify((x, y), psi.diff(y).diff(y), 'numpy')
vt = lambdify((x, y), -psi.diff(x).diff(y), 'numpy')
return u, v, ut, vt
In [ ]:
def plot_streamlines(ax, u, v, ut, vt, xlim=(-1, 1), ylim=(-1, 1)):
x0, x1 = xlim
y0, y1 = ylim
Y, X = np.ogrid[y0:y1:100j, x0:x1:100j]
ut, vt = ut(X, Y), vt(X, Y)
acc = np.sqrt(ut**2, vt**2)
kw = dict(cmap=plt.cm.rainbow, alpha=0.2, zorder=1, shading='gouraud')
cs = ax.pcolormesh(X, Y, acc, vmax=1, **kw)
fig.colorbar(cs, orientation='horizontal')
ax.streamplot(X, Y, u(X, Y), v(X, Y), color='black', zorder=0)
def format_axes(ax):
ax.set_aspect('equal')
ax.figure.subplots_adjust(bottom=0, top=1, left=0, right=1)
ax.xaxis.set_ticks([])
ax.yaxis.set_ticks([])
for spine in ax.spines.values():
spine.set_visible(False)
In [ ]:
radius = 1
psi = cylinder_stream_function(U=1, R=radius)
u, v, ut, vt = velocity_field(psi)
xlim = ylim = (-4, 4)
fig, ax = plt.subplots(figsize=(6, 6))
acc = plot_streamlines(ax, u, v, ut, vt, xlim, ylim)
c = plt.Circle((0, 0), radius=radius, facecolor='cornflowerblue')
ax.add_patch(c)
format_axes(ax)
In [ ]:
psi.simplify()
In [ ]:
u = psi.diff(y)
u.simplify()
In [ ]:
v = -psi.diff(x)
v.simplify()