Python实现圆柱绕流动画：速度向量异常及动画优化求助

我正在编写Python代码实现圆柱绕流的动画，目前流线绘制正常，但绘制速度向量时部分向量异常偏大。同时希望实现速度向量沿流线运动的动画效果，相关代码如下：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# 计算圆柱绕流的流函数
def stream_function(x, y, U, R):
    r = np.sqrt(x**2 + y**2)
    theta = np.arctan2(y, x)
    psi_cylinder = U * (r - R**2 / r) * np.sin(theta)
    psi_freestream = U * y
    psi = np.where(r < R, 0, psi_freestream + psi_cylinder)
    return psi

# 从流函数计算速度分量(u, v)
def velocity_components(x, y, U, R):
    psi_x = -U * (1 - R**2 / x**2) * np.sin(np.arctan2(y, x))
    psi_y = U * (1 + R**2 / y**2) * np.cos(np.arctan2(y, x))
    psi_x = psi_x * ( x ** 2 + y ** 2 >= R **2 )
    psi_y = psi_y * ( x ** 2 + y ** 2 >= R **2 )
    return psi_x, psi_y

# 绘制流线
def plot_streamlines(ax, X, Y, stream_function, U, R):
    psi = stream_function(X, Y, U, R)
    levels = np.linspace(np.min(psi), np.max(psi), 100)
    ax.contour(X, Y, psi, levels=levels, colors='b', linestyles='dashed')
    
    # 绘制圆柱下方的实流线
    ax.contour(X, Y, psi, levels=[0], colors='b')

# 绘制速度向量
def plot_velocity_vectors(ax, X, Y, velocity_components, U, R, frame):
    spacing = 5
    x_points = X[::spacing, ::spacing]
    y_points = Y[::spacing, ::spacing]

    psi_x, psi_y = velocity_components(x_points, y_points, U, R)
    ax.quiver(x_points, y_points, psi_x, psi_y, color='g', scale=0.5, width=0.01)

# 创建画布和坐标轴
fig, ax = plt.subplots()

# 圆柱参数
D = 2  # 圆柱直径（英寸）
R = D / 2  # 圆柱半径
U = 0.01  # 来流速度（英尺/秒）

# 设置绘图范围
ax.set_xlim(-D, D)
ax.set_ylim(-D, D)

# 设置等比例显示
ax.set_aspect('equal')

# 创建流线网格
x = np.linspace(-D, D, 200)
y = np.linspace(-D, D, 200)
X, Y = np.meshgrid(x, y)

# 帧更新函数
def update(frame):
    ax.clear()
    
    # # 绘制流线
    # plot_streamlines(ax, X, Y, stream_function, U, R)
    
    # 绘制速度向量
    plot_velocity_vectors(ax, X, Y, velocity_components, U, R, frame)
    
    # 重置绘图范围
    ax.set_xlim(-D, D)
    ax.set_ylim(-D, D)
    
    # 重置等比例显示
    ax.set_aspect('equal')
    
    # 设置标题
    ax.set_title(f'圆柱绕流 - 帧 {frame}')

# 创建动画
animation = FuncAnimation(fig, update, frames=np.arange(0, 4, 1), interval=100)  # 每秒10帧

# 显示初始图
plt.show()

# 保存动画
animation.save('fluid_flow_over_cylinder_with_vectors.gif', writer='imagemagick', fps=10)

当前绘制效果显示部分速度向量异常偏大，期望优化向量显示并实现速度向量沿流线运动的动画。

1. 修复速度向量异常偏大问题 #

向量异常的核心原因是velocity_components函数的推导错误，误用直角坐标的x²、y²替代了极坐标下的r²。圆柱绕流的速度分量需基于极坐标公式推导：

• 径向速度：(v_r = U\left(1 - \frac{R^2}{r2}\right)\cos\theta)
• 切向速度：(v_\theta = -U\left(1 + \frac{R^2}{r2}\right)\sin\theta)
  再转换为直角坐标分量：
• (u = v_r\cos\theta - v_\theta\sin\theta)
• (v = v_r\sin\theta + v_\theta\cos\theta)

修正后的速度分量函数：

def velocity_components(x, y, U, R):
    r = np.sqrt(x**2 + y**2)
    theta = np.arctan2(y, x)
    # 极坐标速度分量
    v_r = U * (1 - R**2 / r**2) * np.cos(theta)
    v_theta = -U * (1 + R**2 / r**2) * np.sin(theta)
    # 转换为直角坐标
    u = v_r * np.cos(theta) - v_theta * np.sin(theta)
    v = v_r * np.sin(theta) + v_theta * np.cos(theta)
    # 圆柱内部速度置0
    mask = r >= R
    u = np.where(mask, u, 0)
    v = np.where(mask, v, 0)
    return u, v

同时调整quiver的scale参数适配速度量级，比如改为scale=0.02，避免向量过长：

ax.quiver(x_points, y_points, psi_x, psi_y, color='g', scale=0.02, width=0.005)

2. 实现速度向量沿流线运动的动画 #

通过预定义追踪点、每帧更新点位置的方式实现运动效果，完整代码如下：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

def stream_function(x, y, U, R):
    r = np.sqrt(x**2 + y**2)
    theta = np.arctan2(y, x)
    psi_cylinder = U * (r - R**2 / r) * np.sin(theta)
    psi_freestream = U * y
    psi = np.where(r < R, 0, psi_freestream + psi_cylinder)
    return psi

def velocity_components(x, y, U, R):
    r = np.sqrt(x**2 + y**2)
    theta = np.arctan2(y, x)
    v_r = U * (1 - R**2 / r**2) * np.cos(theta)
    v_theta = -U * (1 + R**2 / r**2) * np.sin(theta)
    u = v_r * np.cos(theta) - v_theta * np.sin(theta)
    v = v_r * np.sin(theta) + v_theta * np.cos(theta)
    mask = r >= R
    u = np.where(mask, u, 0)
    v = np.where(mask, v, 0)
    return u, v

def plot_streamlines(ax, X, Y, stream_function, U, R):
    psi = stream_function(X, Y, U, R)
    levels = np.linspace(np.min(psi), np.max(psi), 100)
    ax.contour(X, Y, psi, levels=levels, colors='b', linestyles='dashed')
    ax.contour(X, Y, psi, levels=[0], colors='b')
    # 绘制圆柱
    circle = plt.Circle((0, 0), R, color='gray')
    ax.add_artist(circle)

# 初始化追踪点
def init_trace_points(num_points, D, R):
    points = []
    # 在圆柱上游和两侧生成初始点（避免落在圆柱内部）
    for i in range(num_points):
        x = np.random.uniform(-2*D, D)
        y = np.random.uniform(-D, D)
        while x**2 + y**2 < R**2:
            x = np.random.uniform(-2*D, D)
            y = np.random.uniform(-D, D)
        points.append([x, y])
    return np.array(points)

# 更新追踪点位置
def update_trace_points(points, U, R, dt=0.1):
    x, y = points[:, 0], points[:, 1]
    u, v = velocity_components(x, y, U, R)
    # 按速度更新位置
    x_new = x + u * dt
    y_new = y + v * dt
    # 超出边界的点重置到上游
    reset_mask = (np.abs(x_new) > 2*D) | (np.abs(y_new) > D)
    x_new[reset_mask] = np.random.uniform(-2*D, -D, size=np.sum(reset_mask))
    y_new[reset_mask] = np.random.uniform(-D, D, size=np.sum(reset_mask))
    # 避免进入圆柱内部
    inside_mask = (x_new**2 + y_new**2) < R**2
    x_new[inside_mask] = np.random.uniform(-2*D, -D, size=np.sum(inside_mask))
    y_new[inside_mask] = np.random.uniform(-D, D, size=np.sum(inside_mask))
    return np.column_stack((x_new, y_new))

fig, ax = plt.subplots()

D = 2
R = D / 2
U = 0.5  # 增大速度让动画效果更明显

ax.set_xlim(-2*D, 2*D)
ax.set_ylim(-D, D)
ax.set_aspect('equal')

x = np.linspace(-2*D, 2*D, 200)
y = np.linspace(-D, D, 200)
X, Y = np.meshgrid(x, y)

# 初始化30个追踪点
trace_points = init_trace_points(30, D, R)

def update(frame):
    ax.clear()
    # 绘制流线和圆柱
    plot_streamlines(ax, X, Y, stream_function, U, R)
    
    # 更新追踪点位置
    global trace_points
    trace_points = update_trace_points(trace_points, U, R)
    
    # 绘制追踪点的速度向量和点本身
    u, v = velocity_components(trace_points[:,0], trace_points[:,1], U, R)
    ax.quiver(trace_points[:,0], trace_points[:,1], u, v, color='g', scale=2, width=0.005)
    ax.scatter(trace_points[:,0], trace_points[:,1], color='r', s=10)
    
    ax.set_title(f'圆柱绕流动画 - 帧 {frame}')
    ax.set_xlim(-2*D, 2*D)
    ax.set_ylim(-D, D)
    ax.set_aspect('equal')

animation = FuncAnimation(fig, update, frames=np.arange(0, 50, 1), interval=100)

plt.show()
animation.save('cylinder_flow_animation.gif', writer='pillow', fps=10)

代码说明 #

• 修正了速度分量计算逻辑，彻底解决向量异常问题
• 添加追踪点的初始化与更新逻辑，实现点沿流线运动的效果
• 增大来流速度并调整向量缩放参数，让动画更直观
• 加入圆柱可视化，完善整体画面

内容的提问来源于stack exchange，提问作者remusconnor