2

我在圆柱表面的一部分上有一些方位角循环数据，我想复制这些数据以创建一个封闭的圆柱表面，然后在 (theta,z) 平面上绘制数据。

我的输入网格和压力场在 (x,y,z) 笛卡尔坐标系中如下所示：

输入数据和网格

以下是在 (theta,z) 平面中展开的数据：

压力场展开

所以我对此很满意。当我通过沿 z 轴旋转主题来复制数据时，我遇到了不连续 theta 值的问题（请参见下面的代码）。例如，以下是仅 2 次重复获得的结果：

(theta,z) 平面中的重复数据 重复数据

这是我用来获得此结果的代码：

# coding: utf-8

import matplotlib.pyplot as plt
import numpy as np

# Define some useful functions
def cartesianToPolar(x, y):
    '''Transform cartesian data (x, y) into polar coordinates (r, theta)'''
    r = np.sqrt(x**2 + y**2)
    theta = np.arctan2(y, x)
    return r, theta 

def rotate(x,y,phi):
    '''Rotate the points by an angle phi in degree'''
    phiRad = np.deg2rad(phi)
    xRotated = x*np.cos(phiRad) - y*np.sin(phiRad)
    yRotated = x*np.sin(phiRad) + y*np.cos(phiRad)
    return xRotated,yRotated

# Load the 3d data
x,y,z,p = np.loadtxt('data.dat',unpack=False)

# Compute the cylindrical coordinates
r, theta = cartesianToPolar(x, y)

# Define the mesh caracteristics
nbTheta = 21 # Number of points in the azimuthal direction
nbBlades = 3 # Number of blades for the duplication

# Generate structured arrays
ix = np.lexsort((y,x,z))
ny = nbTheta
nx = np.size(x)/ny

X = np.reshape(x[ix],(nx,ny))
Y = np.reshape(y[ix],(nx,ny))
Z = np.reshape(z[ix],(nx,ny))
P = np.reshape(p[ix],(nx,ny))

# 3D plot of the input wall surface
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_aspect('equal', 'datalim')
ax.set_xlabel('Z')
ax.set_ylabel('X')
ax.set_zlabel('Y')
ax.plot_surface(Z,X,Y,rstride=1, cstride=1,color='b',facecolors=cm.jet(P/P.max()))
plt.show()

# Duplicate the data 
xunroll = np.copy(x)
yunroll = np.copy(y)
zunroll = np.copy(z)
punroll = np.copy(p)

for i in range(1,nbBlades):
    x, y = rotate(x,y,45)
    xunroll = np.append(xunroll,x)
    yunroll = np.append(yunroll,y)
    punroll = np.append(punroll,p)
    zunroll = np.append(zunroll,z)

runroll, tunroll = cartesianToPolar(xunroll, yunroll) # Compute the radial and azimuthal coordinates

# Generate the structured arrays after duplication
ix = np.lexsort((tunroll,zunroll))
ny = nbTheta*nbBlades
nx = np.size(xunroll)/ny

Xunroll = np.reshape(xunroll[ix],(nx,ny))
Yunroll = np.reshape(yunroll[ix],(nx,ny))
Zunroll = np.reshape(zunroll[ix],(nx,ny))
Tunroll = np.reshape(tunroll[ix],(nx,ny))
Punroll = np.reshape(punroll[ix],(nx,ny))

# Plot the unrolled duplicated data in the (theta,z) plane
fig = plt.figure()
plt.gca().set_aspect('equal', 'datalim')
plt.pcolormesh(Zunroll/0.018,Tunroll,Punroll,edgecolors='gray',lw=0.05)
plt.xlabel('z/R')
plt.ylabel(r'$\theta$ [rad]')
cbar = plt.colorbar()
cbar.set_label('Pressure',labelpad=10)
plt.grid()

# 3d plot of the rolled duplicated data
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_aspect('equal', 'datalim')
ax.set_xlabel('Z')
ax.set_ylabel('X')
ax.set_zlabel('Y')
ax.plot_surface(Zunroll,Xunroll,Yunroll,rstride=1, cstride=1,alpha=0.5,facecolors=cm.jet(Punroll/Punroll.max()))
plt.show()

请注意，我使用了 numpy 中的 arctan2 函数来计算方位角 theta 角，并且我的网格是结构化的。

最后，我想让我的数据完全旋转 -pi <= theta <= pi (8 x input data) 并在 (theta,z) 平面上展开。

如果你想和他们一起玩，这里是数据。

任何帮助都将受到热烈欢迎:)

此致，

弗朗索瓦

4

0 回答 0