python - 如何在Python中绘制圆柱体周围流线流的流线图？

Question

我编写了一个程序，它需要半径范围(R 到 5)和 theta 范围(0 到 2pi)。它将极坐标转换为笛卡尔坐标。然后我绘制了 x,y 笛卡尔坐标的散点图。

我有一个自由流速度 u(x 方向)=1 和 v 速度(y 方向)= 0。我想为上面获得的点绘制流图。我知道我必须使用网格函数。但经过多次尝试，我无法将网格网格点映射到速度值。

我想要得到如下图所示的情节。

代码如下:

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

错误是

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

我花了很多时间在streamplots上，我想尝试一下你的问题。这是我到目前为止所得到的。

import numpy as np
import matplotlib.pyplot as plt

我同意@swag2198的评论，所以我从笛卡尔网格开始。然后我找到了对应的极坐标网格。

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

然后我意识到问题缺少向量场方程。所以我做了一些谷歌搜索并发现the Wikipedia page for the problem 。下面的公式给出了向量场$$ V $$:其中 $$\psi $$ 是 .

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

然后我前往 Mathematica 计算出 $$\del\psi $$ 。然后我用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))

然后，我利用 streamplot 不绘制 NaN 值的事实，从矢量场中删除了圆盘和外部区域。

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

这是最终产品:

虽然 $$ x < 0 $$ half 中的向量指向相反的方向，但我无法弄清楚为什么会出现这种情况。也许，我在坐标变换或函数上犯了一些错误。

以下是轻松复制粘贴的完整代码:

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