Animation of fluid flow over a cylinder

Question

I am trying to write a Python code to animate fluid flow over a cylinder. I can plot the streamlines fine but I am trying to draw velocity vectors on the plot. I am having problems with some of the vectors, as they appear very large for some reason.

The code:

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

# Function to calculate the stream function for flow over a cylinder
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

# Function to calculate the velocity components (u, v) from the stream function
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

# Function to plot streamlines
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')
    
    # Solid streamlines below the cylinder
    ax.contour(X, Y, psi, levels=[0], colors='b')

# Function to plot velocity vectors
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)

# Create a figure and axis
fig, ax = plt.subplots()

# Cylinder properties
D = 2  # Diameter of the cylinder in inches
R = D / 2  # Radius of the cylinder
U = 0.01  # Fluid velocity in feet/second

# Set plot limits
ax.set_xlim(-D, D)
ax.set_ylim(-D, D)

# Set aspect ratio to be equal
ax.set_aspect('equal')

# Create a grid for the streamlines
x = np.linspace(-D, D, 200)
y = np.linspace(-D, D, 200)
X, Y = np.meshgrid(x, y)

# Function to update the plot for each frame
def update(frame):
    ax.clear()
    
    # # Plot streamlines
    # plot_streamlines(ax, X, Y, stream_function, U, R)
    
    # Plot velocity vectors
    plot_velocity_vectors(ax, X, Y, velocity_components, U, R, frame)
    
    # Set plot limits
    ax.set_xlim(-D, D)
    ax.set_ylim(-D, D)
    
    # Set aspect ratio to be equal
    ax.set_aspect('equal')
    
    # Set title
    ax.set_title(f'Fluid Flow Over a Cylinder - Frame {frame}')

# Create the animation
animation = FuncAnimation(fig, update, frames=np.arange(0, 4, 1), interval=100)  # 10 frames per second

# Display the initial plot
plt.show()

# Display the animation
animation.save('fluid_flow_over_cylinder_with_vectors.gif', writer='imagemagick', fps=10)  # Save the animation as a gif

Below is what the plot looks like. I would also like it to be animated, possibly by showing the velocity vectors move along the streamlines?

enter image description here

Answer 1

Your streamfunction is slightly wrong - it doesn't need the free-stream velocity added, as the free-stream velocity is what it already tends to at large r.

Your velocity functions are both wrong.

Since you ultimately want streamfunction and velocity at specific (x,y) rather than (r,theta) I think you would be better re-expressing your stream function in terms of x and y: that is, psi = Uy(1-R^2/r^2), where r^2=x^2+y^2. Then it makes it easier to find velocity components u=d.psi/d.y and v=-d.psi/dx.

I've removed the animation from the code below, but you can easily reinstate it when you are happy with everything else. Actually, I think you would probably get a better animation by moving the cylinder than having a steady onset flow, but that's your choice.

import numpy as np
import matplotlib.pyplot as plt

# Function to calculate the stream function for flow over a cylinder
def stream_function(x, y, U, R):
    rsq = x ** 2 + y ** 2
    psi = U * y * ( 1.0 - R ** 2 / rsq )
    return np.where( rsq >= R * R, psi, 0.0 )

# Function to calculate the velocity components (u, v) from the stream function
def velocity_components(x, y, U, R):
    rsq = x ** 2 + y ** 2
    Rsq = R ** 2
    u = U * ( 1 - Rsq / rsq + 2 * Rsq * y ** 2 / ( rsq * rsq ) )
    v = -2 * U * Rsq * x * y / ( rsq * rsq )
    u = np.where( rsq >= Rsq, u, 0.0 )
    v = np.where( rsq >= Rsq, v, 0.0 )
    return u, v

# Function to plot streamlines
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), 50)
    ax.contour(X, Y, psi, levels=levels, colors='b' )

# Function to plot velocity vectors
def plot_velocity_vectors(ax, X, Y, velocity_components, U, R ):
    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='r', scale=0.3 )


# Create a figure and axis
fig, ax = plt.subplots()

# Cylinder properties
D = 2  # Diameter of the cylinder
R = D / 2  # Radius of the cylinder
U = 0.01  # Fluid velocity

# Set plot limits
ax.set_xlim(-D, D)
ax.set_ylim(-D, D)

# Set aspect ratio to be equal
ax.set_aspect('equal')

# Create a grid for the streamlines
x = np.linspace( -D, D, 200 )
y = np.linspace( -D, D, 200 )
X, Y = np.meshgrid(x, y)

#plot_streamlines( ax, X, Y, stream_function, U, R )
plot_velocity_vectors( ax, X, Y, velocity_components, U, R )
    
ax.set_xlim(-D, D)
ax.set_ylim(-D, D)
ax.set_aspect('equal')
ax.set_title(f'Fluid Flow Over a Cylinder')

plt.show()

Streamfunction:

Velocity vectors: