1

我试图在一个身体周围的流线上放置实验室，该身体的对称轮廓是由涡流和均匀流动产生的，到目前为止我必须得到这样的东西（带有标签）

涡流加均匀流

我用下一个代码得到它：

import numpy as np
import matplotlib.pyplot as plt

vortex_height = 18.0
h = vortex_height
vortex_intensity = 55.0
cv = vortex_intensity
permanent_speed = 10.0
U1 = permanent_speed

Y, X = np.mgrid[-25:25:100j, -25:25:100j]
U = 5.0 +  37.0857 * (Y - 18.326581) / (X ** 2 + (Y - 18.326581) ** 2) +- 37.0857 * (Y + 18.326581) / (X ** 2 + (Y + 18.326581) ** 2)
V = - 37.0857 * (X) / (X ** 2 + (Y - 18.326581) ** 2) + 37.0857 * (X) / (X ** 2 + (Y + 18.326581) ** 2)
speed = np.sqrt(U*U + V*V)

plt.streamplot(X, Y, U, V, color=U, linewidth=2, cmap=plt.cm.autumn)
plt.colorbar()
#CS = plt.contour(U, v, speed)
#plt.clabel(CS, inline=1, fontsize=10)
#f, (ax1, ax2) = plt.subplots(ncols=2)
#ax1.streamplot(X, Y, U, V, density=[0.5, 1])

#lw = 5*speed/speed.max()
#ax2.streamplot(X, Y, U, V, density=0.6, color='k', linewidth=lw)
plt.savefig("stream_plot5.png")
plt.show()

所以我正在更改下一个示例代码（使用 pylab）：

from numpy import exp,arange
from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show

# the function that I'm going to plot
def z_func(x,y):
    return (1-(x**2+y**3))*exp(-(x**2+y**2)/2)

x = arange(-3.0,3.0,0.1)
y = arange(-3.0,3.0,0.1)
X,Y = meshgrid(x, y) # grid of point
Z = z_func(X, Y) # evaluation of the function on the grid

im = imshow(Z,cmap=cm.RdBu) # drawing the function
# adding the Contour lines with labels
cset = contour(Z,arange(-1,1.5,0.2),linewidths=2,cmap=cm.Set2)
clabel(cset,inline=True,fmt='%1.1f',fontsize=10)
colorbar(im) # adding the colobar on the right
# latex fashion title
title('$z=(1-x^2+y^3) e^{-(x^2+y^2)/2}$')
show()

有了这个情节：

样本图

最后我得到它是这样的：

import numpy as np
from numpy import exp,arange,log
from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show

# PSI = streamline
def streamLine(x, y, U = 5, hv = 18.326581, cv = 37.0857):
    x2 = x ** 2
    y2plus = (y + hv) ** 2
    y2minus = (y - hv) ** 2
    PSI_1 = U * y
    PSI_2 =   0.5 * cv * log(x2 + y2minus)
    PSI_3 = - 0.5 * cv * log(x2 + y2plus)
    psi = PSI_1 + PSI_2 + PSI_3
    return psi
"""
def streamLine(x, y):
    return  0.5 * 37.0857 * log(x ** 2 + (y - 18.326581) ** 2)
# (5.0 * y + 0.5 * 37.0857 * math.log(x ** 2 + (y - 18.326581) ** 2) - 0.5 * 37.0857 * math.log(x ** 2 + (y + 18.326581) ** 2))
"""
x = arange(-20.0,20.0,0.1)
y = arange(-20.0,20.0,0.1)
X,Y = meshgrid(x, y) # grid of point
#Z = z_func(X, Y) # evaluation of the function on the grid
Z= streamLine(X, Y)

im = imshow(Z,cmap=cm.RdBu) # drawing the function
# adding the Contour lines with labels
cset = contour(Z,arange(-5,6.5,0.2),linewidths=2,cmap=cm.Set2)
clabel(cset,inline=True,fmt='%1.1f',fontsize=9)
colorbar(im) # adding the colobar on the right
# latex fashion title
title('$phi= 5.0 y + (1/2)* 37.0857 log(x^2 + (y - 18.326581)^2)-(1/2)* 37.085...$')
show()
#print type(Z)
#print len(Z)

但后来我得到了下一个情节：

result_of_vortex_plus_uniform_flow

这让我想知道出了什么问题，因为轴不在它们应该在的位置。

4

1 回答 1

3

contour()绘制标量场的等高线，并且streamplot()是矢量场的绘制，可以使用梯度算子从标量场构造矢量场。

这是一个例子：

import numpy as np
from numpy import exp,arange,log
from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show,streamplot

# PSI = streamline
def f(x, y, U = 5, hv = 18.326581, cv = 37.0857):
    x2 = x ** 2
    y2plus = (y + hv) ** 2
    y2minus = (y - hv) ** 2
    PSI_1 = U * y
    PSI_2 =   0.5 * cv * log(x2 + y2minus)
    PSI_3 = - 0.5 * cv * log(x2 + y2plus)
    psi = PSI_1 + PSI_2 + PSI_3
    return psi

x = arange(-20.0,20.0,0.1)
y = arange(-20.0,20.0,0.1)
X,Y = meshgrid(x, y) # grid of point
#Z = z_func(X, Y) # evaluation of the function on the grid
Z= f(X, Y)

dx, dy = 1e-6, 1e-6
U = (f(X+dx, Y) - f(X, Y))/dx
V = (f(X, Y+dy) - f(X, Y))/dy
streamplot(X, Y, U, V, linewidth=1, color=(0, 0, 1, 0.3))

cset = contour(X, Y, Z,arange(-20,20,2.0),linewidths=2,cmap=cm.Set2)
clabel(cset,inline=True,fmt='%1.1f',fontsize=9)
colorbar(im) # adding the colobar on the right


# latex fashion title
title('$phi= 5.0 y + (1/2)* 37.0857 log(x^2 + (y - 18.326581)^2)-(1/2)* 37.085...$')

show()

输出：

在此处输入图像描述

于 2013-10-28T05:01:21.447 回答