极坐标图马格纳斯效应未显示正确的d

2024-03-28 23:30:30 发布

您现在位置:Python中文网/ 问答频道 /正文

我想在极坐标图上画出绕旋转圆柱体流动的速度方程。(方程式来自安徒生的《空气动力学基础》),你可以在for循环语句中看到这两个方程式。你知道吗

我无法大声地把计算出的数据表示在极坐标图上。我已经尝试了我的每一个想法,但都一无所获。我检查了数据,这一个似乎都是正确的,因为它的行为应该如何。你知道吗

这是我最后一次尝试的代码:

import numpy as np
import matplotlib.pyplot as plt


RadiusColumn = 1.0
VelocityInfinity = 10.0
RPM_Columns         = 0.0#
ColumnOmega         = (2*np.pi*RPM_Columns)/(60)#rad/s
VortexStrength      = 2*np.pi*RadiusColumn**2 * ColumnOmega#rad m^2/s

NumberRadii = 6
NumberThetas = 19

theta = np.linspace(0,2*np.pi,NumberThetas)
radius = np.linspace(RadiusColumn, 10 * RadiusColumn, NumberRadii)


f = plt.figure()
ax = f.add_subplot(111, polar=True)

for r in xrange(len(radius)):
    for t in xrange(len(theta)):


        VelocityRadius = (1.0 - (RadiusColumn**2/radius[r]**2)) * VelocityInfinity * np.cos(theta[t])
        VelocityTheta = - (1.0 + (RadiusColumn**2/radius[r]**2))* VelocityInfinity * np.sin(theta[t]) - (VortexStrength/(2*np.pi*radius[r]))
        TotalVelocity = np.linalg.norm((VelocityRadius, VelocityTheta))


        ax.quiver(theta[t], radius[r], theta[t] + VelocityTheta/TotalVelocity, radius[r] + VelocityRadius/TotalVelocity)


plt.show()

如您所见,我现在已经将RPM设置为0。这意味着水流应该从左到右,并且在横轴上是对称的。(水流应该以同样的方式在圆柱体两侧流动。)但是结果看起来更像这样:

Polar plot result - not correct!

这完全是胡说八道。似乎有一个漩涡,即使没有设置!更奇怪的是,当我只显示从0到pi/2的数据时,流就变了!你知道吗

Polar plot result first quadrant - total mayhem

正如你从代码中看到的,我尝试过使用单位向量,但显然这不是方法。如有任何有用的意见,我将不胜感激。你知道吗

谢谢!你知道吗


Tags: fornppipltrpm方程式thetaradius
1条回答
网友
1楼 · 发布于 2024-03-28 23:30:30

基本问题是,极坐标对象的.quiver方法仍然需要笛卡尔坐标中的向量分量,因此您需要自己将θ和径向分量转换为x和y:

for r in xrange(len(radius)):
    for t in xrange(len(theta)):

        VelocityRadius = (1.0 - (RadiusColumn**2/radius[r]**2)) * VelocityInfinity * np.cos(theta[t])
        VelocityTheta = - (1.0 + (RadiusColumn**2/radius[r]**2))* VelocityInfinity * np.sin(theta[t]) - (VortexStrength/(2*np.pi*radius[r]))
        TotalVelocity = np.linalg.norm((VelocityRadius, VelocityTheta))

        ax.quiver(theta[t], radius[r],
                  VelocityRadius/TotalVelocity*np.cos(theta[t])
                  - VelocityTheta/TotalVelocity*np.sin(theta[t]),
                  VelocityRadius/TotalVelocity*np.sin(theta[t])
                  + VelocityTheta/TotalVelocity*np.cos(theta[t]))

plt.show()

fixed figure

但是,通过使用矢量化,您可以大大改进代码:您不需要在每个点上循环以获得所需的内容。与你的代码相当,但更清楚:

def pol2cart(th,v_th,v_r):
    '''convert polar velocity components to Cartesian, return v_x,v_y'''

    return v_r*np.cos(th) - v_th*np.sin(th), v_r*np.sin(th) + v_th*np.cos(th)


theta = np.linspace(0,2*np.pi,NumberThetas,endpoint=False)
radius = np.linspace(RadiusColumn, 10 * RadiusColumn, NumberRadii)[:,None]

f = plt.figure()
ax = f.add_subplot(111, polar=True)

VelocityRadius = (1.0 - (RadiusColumn**2/radius**2)) * VelocityInfinity * np.cos(theta)
VelocityTheta = - (1.0 + (RadiusColumn**2/radius**2))* VelocityInfinity * np.sin(theta) - (VortexStrength/(2*np.pi*radius))
TotalVelocity = np.linalg.norm([VelocityRadius, VelocityTheta],axis=0)

VelocityX,VelocityY = pol2cart(theta,VelocityTheta,VelocityRadius)

ax.quiver(theta,radius, VelocityX/TotalVelocity, VelocityY/TotalVelocity)

plt.show()

fixed, vectorized final figure

很少有显著变化:

  • 我将endpoint=False添加到theta:因为函数在2*pi中是周期性的,所以不需要绘制两次端点。请注意,这意味着当前有更多可见的箭头;如果您想要原始行为,我建议您将NumberThetas减少1。你知道吗
  • 我在[:,None]中添加了radius:这将使它成为一个二维数组,因此后面定义速度的操作将为您提供二维数组:不同的列对应不同的角度,不同的行对应不同的半径。quiver与数组值输入兼容,因此只需调用quiver即可完成您的工作。你知道吗
  • 由于速度现在是二维数组,我们需要调用np.linalg.norm基本上是一个三维数组,但如果我们指定一个轴来处理,这是可以预期的。你知道吗
  • 我定义了pol2cart辅助函数来完成从极坐标分量到笛卡尔分量的转换;这不是必需的,但我觉得这样更清楚。你知道吗

最后一句话:我建议选择较短的变量名,以及没有大小写的变量名。这可能会让你的编码速度也更快。你知道吗

相关问题 更多 >