ポーラプロットマグナス効果が正しいデータを表示していません

Question

回転円柱の周りの流れの速度方程式を極座標プロット上にプロットしたいと考えました。 （方程式はAndersenによる "Airodynamics of Fundamentals of Andersen"からのものです。）forループ文の中に2つの方程式を見ることができます。

私は大声で叫んで、計算されたデータを極座標プロットに表示することはできません。私はすべての考えを試したが、どこにも到着しなかった。私はデータをチェックしました、そして、これはどのようにすべきか振る舞うので、これはすべて正しいと思われます。ここで

私の最後の試みのコード：あなたが見ることができるように、私は0に今のRPMを設定している

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

流れは左から右に移動し、全体で対称でなければならないことを意味し横軸。 （。流れがシリンダの周りに両側で同じ道を行く必要があります）結果は、しかし、より次のようになります。

これは完全にナンセンスです。設定されていない場合でも、渦があるようです！私が0からpi/2までのデータしか表示しない場合でも、フローは変化します。

あなたがコードから見ることができるように、私は単位ベクトルを利用することを試みているが、明確にこれが進むべき道ではありません。私は有用な入力を感謝します。

ありがとうございます！

Answer 1

基本的な問題を自分で極性Axesオブジェクトの.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() 

ただし、ベクトル化を使用することでコードを大幅に改善することができます。必要なものを得るために各点をループする必要はありません。だからあなたのコードの同等、それでも明確：

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

少数の注目すべき変更：

• 私はthetaにendpoint=Falseを追加しました：あなたの関数は2*piで周期的であることから、あなたがドン」エンドポイントを2回プロットする必要がありません。これは現在、より多くのが表示されていることを意味します。 arrows;元の動作をしたい場合は、NumberThetasを1つ減らすことをお勧めします。
• [:,None]〜radiusを追加しました。これは2d配列になります。速度の定義における後の操作で2d配列が得られます。異なる列は異なる角度に対応し、異なる行は異なる半径に対応します。 は配列値の入力と互換性がありますので、quiverを一度呼び出して作業します。
• ベロシティが2次元配列になっているので、基本的には3次元配列のnp.linalg.normを呼び出す必要がありますが、処理する軸を指定すると期待通りに機能します。
• 極座標からデカルトのコンポーネントへの変換を行う補助関数を定義しました。これは必ずしも必要ではありませんが、私にはこのように明確に見えます。

最終的なコメント：より短い変数名とCamelCaseを持たないものを選択することをお勧めします。おそらくあなたのコーディングも速くなるでしょう。