notebook.community

Edit and run



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}$

$U \rightarrow$ velocidade horizontal;
$R \rightarrow$ raio do cilindro;
$r, \theta \rightarrow$ coordenadas polar .



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)

\begin{align*} u &= + \dfrac{\partial\psi}{\partial y} \\ v &= - \dfrac{\partial\psi}{\partial x} \end{align*}



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()