2037.sagews

Circulation around a Rigid Body
Joris Vankerschaver

 

In this worksheet, we take a look at potential flow around a cylinder with circulation. We study the effect of varying the circulation on the topology of the streamlines. The standard reference for this system is Batchelor, paragraph 6.6.

x, y, U, Gamma, R = var('x, y, U, Gamma, R')

Sage has some capabilities to calculate the stream function directly from the complex potential. However, in the interest of speed it is best to define psi directly.

psi(x, y, U, Gamma, R) = U*(1 - R^2/(x^2 + y^2))*y - Gamma/(4*pi)*log(x^2 + y^2)

import matplotlib.cm
from sage.ext.fast_eval import fast_float

# We now fix U = 1.0 and R = 1.0 and let Gamma be a free parameter 
# 
# (The 'fast_float' command below optimizes the stream function for fast evaluation)

psi_fast = fast_float(psi(x, y, 1.0, Gamma, 1.0), 'x', 'y', 'Gamma')

Note that psi becomes unbounded when we approach the origin. If we were to feed psi to the contour_plot command directly, this singularity would obscure all other features of the plot. To avoid this, we define a new function which returns infinity when evaluated on a point inside the circle. These points will be left out of the contour_plot.

# The following is a helper function which does all the dirty work for us:

def plot_body_flow(Gamma):
    psi_body = lambda x, y: psi_fast(x, y, Gamma) if x^2 + y^2 > 1.0 else oo
    contour_plot(psi_body, (x, -2, 2), (y, -2, 2), cmap='jet', contours=50, aspect_ratio=1).show()

# For Gamma = 0.0, we clearly observe steady flow around a circular body 

plot_body_flow(0.0)

# We clearly see the effect of increasing the circulation:

plot_body_flow(8.0)

# When Gamma = 4*pi, the solution changes quantitatively: the two stagnation points on the boundary of the body coalesce.

plot_body_flow(4*pi)

# When Gamma is increased even further, the stagnation point moves off the boundary of the rigid body and into the flow.
# The streamlines close to the rigid body are closed curves, while far away they are unbounded.

plot_body_flow(4*pi + 0.5)