How to plot a streamplot in python for a streamline flow arround a cylinder?

Question

I have written a program where it takes a radius range(R to 5) and theta range(0 to 2pi). It converts the polar coordinates to Cartesian coordinates. Then I have plotted the scatter graph of the x,y Cartesian coordinates.

I have a free stream velocity u(x-direction)=1 and v-velocity(y-direction)= 0. I want to plot a stream plot for the points I have obtained above. I know I have to use meshgrid function. But after trying a lot, I cannot map the Mesh grid points to the velocity values.

I want to get a plot like below picture.

The Code is given below:

import math
import numpy as np
from matplotlib import pyplot


R=1.15;  # Radius of circle
Nr=75;    # No of radial points
Ntheta=75; # No of theta points
x_start,x_end=-5,5;  
y_start,y_end=-5,5;
theta=np.linspace(0,2*np.pi,Ntheta); # Array of theta points
radial=np.linspace(R,5,Nr);          # Array of radius points
xc=-0.15; # Center of circle.
yc=0.0;   # Center of circle.

X,Y=np.meshgrid(r,theta);   # Meshgrid for stream plot

x=np.zeros((Nr*Ntheta,1),dtype=np.float64);   # To store x cordinates
y=np.zeros((Nr*Ntheta,1),dtype=np.float64);   # to store y cordinates


# Z-Plane Computation
cnt=0;
for i in range(Nr):
    for j in range(Ntheta):
        x[cnt,0]=radial[i]*np.cos(theta[j])+xc;    # Calculation of Cartesian Cordinates
        y[cnt,0]=radial[i]*np.sin(theta[j])+yc;    # Calculation of Cartesian Cordinates
        cnt+=1;
 
# Plot
fig=pyplot.figure(figsize=(10,10));
pyplot.scatter(x[:,0],y[:,0],s=1,color='k');
pyplot.xlim(-6,6);
pyplot.ylim(-6,6);
pyplot.scatter(xc,yc,s=80,color='g',marker='o');
pyplot.title('Z-Plane',fontsize=20);
pyplot.xlabel('x',fontsize=15);
pyplot.ylabel('y',fontsize=15);
pyplot.grid(color='k',which='both',axis='both',linestyle='--',linewidth=0.5);

# StreamLine Velocity
u_inf=1;
u_freestream= u_inf*np.ones((Nr,Ntheta),dtype=np.float64);
v_freestream= np.zeros((Nr,Ntheta),dtype=np.float64);

#plotting
pyplot.figure()
pyplot.streamplot(X,Y,u_freestream,v_freestream);

The error is

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-59df08985fbe> in <module>()
     46 #plotting
     47 pyplot.figure()
---> 48 pyplot.streamplot(X,Y,u_freestream,v_freestream);
     49 
     50 

3 frames
/usr/local/lib/python3.6/dist-packages/matplotlib/streamplot.py in __init__(self, x, y)
    342 
    343         if not np.allclose(np.diff(x), self.width / (self.nx - 1)):
--> 344             raise ValueError("'x' values must be equally spaced")
    345         if not np.allclose(np.diff(y), self.height / (self.ny - 1)):
    346             raise ValueError("'y' values must be equally spaced")

ValueError: 'x' values must be equally spaced

Answer 1

I am spending a lot of time with streamplots, and I thought of taking a shot at your problem. Here is what I got so far.

import numpy as np
import matplotlib.pyplot as plt

I agree with @swag2198's comment, so I started with a cartesian grid. Then I found the corresponding polar grid.

# init cartesian grid
Y, X = np.mgrid[-5:5:101j, -5:5:101j]

# get corresponding polar grid
R, T = np.sqrt(X**2 + Y**2), np.arctan(Y/X)

Then I realized that the question is missing vector field equations. So I did some googling and found the Wikipedia page for the problem. The following formula gives the vector field $$ V $$: where $$ \psi $$ is .

# stream function
U = 1 # velocity(?) constant
R_disc = 1.15 # radius of the disc
S_polar = U * (R - (R_disc**2)/R) * np.sin(T)

Then I headed to Mathematica figured out $$ \del \psi $$ . Then I wrote it in python.

# substitute and differentiate
denom = np.power(X*X + Y*Y, 3/2) * (np.sqrt(1 + ((Y*Y)/(X*X))))
Sx = 2 * (R_disc**2) * U * Y / denom
Sy = U * (R_disc**2 * (Y*Y - X*X) + (X*X + Y*Y)**2) * (1 / (X*denom))

Then, I removed the disc and the outer region from the vector fields using the fact that streamplot doesn't plot NaN values.

# remove regions
Sx[R<R_disc]=np.nan
Sy[R<R_disc]=np.nan
Sx[R>5]=np.nan
Sy[R>5]=np.nan

Finally, the last part of plotting.

# plot the figure
fig = plt.figure(figsize=(5,5))
ax = plt.subplot()
ax.streamplot(X, Y, Sy, -Sx)
disc = plt.Circle((0,0), radius=R_disc, fc='k')
plt.gca().add_patch(disc)
plt.show()

Here is the final product:

Though the vectors in $$ x < 0 $$ half point in the opposite direction, I could not figure out why this is the case. Perhaps, I made some mistake in coordinate transforms or the function.

Here is the full code for easy copy-paste:

import numpy as np
import matplotlib.pyplot as plt

# init cartesian grid
Y, X = np.mgrid[-5:5:101j, -5:5:101j]

# get corresponding polar grid
R, T = np.sqrt(X**2 + Y**2), np.arctan(Y/X)

# stream function
U = 1 # velocity(?) constant
R_disc = 1.15 # radius of the disc
S_polar = U * (R - (R_disc**2)/R) * np.sin(T)

# substitute and differentiate
denom = np.power(X*X + Y*Y, 3/2) * (np.sqrt(1 + ((Y*Y)/(X*X))))
Sx = 2 * (R_disc**2) * U * Y / denom
Sy = U * (R_disc**2 * (Y*Y - X*X) + (X*X + Y*Y)**2) * (1 / (X*denom))

# remove regions
Sx[R<R_disc]=np.nan
Sy[R<R_disc]=np.nan
Sx[R>5]=np.nan
Sy[R>5]=np.nan

# plot the figure
fig = plt.figure(figsize=(5,5))
ax = plt.subplot()
ax.streamplot(X, Y, Sy, -Sx)
disc = plt.Circle((0,0), radius=R_disc, fc='k')
plt.gca().add_patch(disc)
plt.show()