import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook
import jmkfigure

Key Assignment 8

Map from circle radius $a$ to straight line using Jhukowski transformation

x = np.arange(-5, 5, 0.1)
y = x.copy()
X, Y = np.meshgrid(x, y)
a = 1

fig, axs = plt.subplots(1, 2, sharex=True, sharey=True)
theta = np.arange(0, np.pi*2, 0.01)
circle = np.exp(1j*theta) * a
axs[0].plot(np.real(circle), np.imag(circle))
axs[0].set_aspect(1)
axs[0].set_ylabel('$\eta / a$')
axs[0].set_xlabel('$\zeta / a$')
axs[0].set_title('$\chi$')

# transformed
circleZ = circle + a**2/circle
axs[1].plot(np.real(circleZ), np.imag(circleZ))
axs[1].set_aspect(1)
axs[1].set_xlim([-4.1, 4.1])

axs[1].set_ylim([-4.1, 4.1])

axs[1].set_xlabel('$x / a$')
axs[1].set_ylabel('$y / a$')
axs[1].set_title('$z$')

Text(0.5, 1.0, '$z$')

Contour the flow around both objects. Note that we need the inverse transformation to get to the new frame and that this transformation requires different branches of the quadratic on either side of $x=0$. The inverse transform is simply $\zeta = \frac{1}{2}\left( z+\sqrt{z^2-4a^2} \right)$, $x>0$ and $\zeta = \frac{1}{2}\left( z-\sqrt{z^2-4a^2} \right)$, $x<0$.

Z = X + 1j * Y
alpha = 20 * np.pi / 180
gamma = 0
U = 10
w = U * (Z * np.exp(-1j * alpha) + a**2 / Z * np.exp(1j * alpha) - 1j * gamma / 2 / np.pi * np.log(Z))
psi = np.imag(w)

# make the same thing in the new frame:

Znew = Z/2.+0.5*np.sqrt(Z**2-4*a**2)
Znew[X<0.] = Z[X<0.]/2.-0.5*np.sqrt(Z[X<0.]**2-4*a**2)

wnew = U * (Znew * np.exp(-1j * alpha) + a**2 / Znew * np.exp(1j * alpha) - 1j * gamma / 2 / np.pi * np.log(Znew))
psinew = np.imag(wnew)

fig, axs = plt.subplots(1, 2, sharex=True, sharey=True)

axs[0].plot(np.real(circle), np.imag(circle), color='k')
axs[0].set_aspect(1)
axs[0].set_ylabel('$\eta / a$')
axs[0].set_xlabel('$\zeta / a$')
axs[0].set_title('$\chi$')
axs[0].contour(x, y, psi, levels=np.arange(-U*10, U*10, U / 5))

# transformed
circleZ = circle + a**2 / circle
axs[1].plot(np.real(circleZ), np.imag(circleZ), color='k')
axs[1].set_aspect(1)
axs[1].contour(x, y, psinew, levels=np.arange(-U*10, U*10, U / 5))

axs[1].set_xlim([-4.1, 4.1])
axs[1].set_ylim([-4.1, 4.1])

axs[1].set_xlabel('$x / a$')
axs[1].set_ylabel('$y / a$')
axs[1].set_title('$z$')

Text(0.5, 1.0, '$z$')

To get the trailing edge to have the stagnation point, we need to induce sufficient circulation $-\Gamma$ so that the stagnation point at $\alpha$ is pushed down to 0 degrees. The flow along the sphere is $$u_{\theta} = -2U\sin(\theta-\alpha) + \frac{\Gamma}{2\pi a} $$, and for the separation point to be at the trailing edge we want $u_{\theta=0}=0$ or $\Gamma = - 4\pi U a \sin(\alpha)$.

Z = X + 1j * Y
alpha = 20 * np.pi / 180
U = 1

gamma = -4 * np.pi * a * U * np.sin(alpha)

w = U * (Z * np.exp(-1j * alpha) + a**2 / Z * np.exp(1j * alpha) - 1j * gamma / 2 / np.pi * np.log(Z))
psi = np.imag(w)

# make the same thing in the new frame:
Znew = 0.5*Z + np.sqrt(0.25*np.abs(Z)**2 - a**2)
# Znew =  Z + a**2 / Z

Znew = Z/2.+0.5*np.sqrt(Z**2-4*a**2)
Znew[X<0.] = Z[X<0.]/2.-0.5*np.sqrt(Z[X<0.]**2-4*a**2)

wnew = U * (Znew * np.exp(-1j * alpha) + a**2 / Znew * np.exp(1j * alpha) - 1j * gamma / 2 / np.pi * np.log(Znew))
psinew = np.imag(wnew)

fig, axs = plt.subplots(1, 2, sharex=True, sharey=True)

axs[0].plot(np.real(circle), np.imag(circle), color='k')
axs[0].set_aspect(1)
axs[0].set_ylabel('$\eta / a$')
axs[0].set_xlabel('$\zeta / a$')
axs[0].set_title('$\chi$')
axs[0].contour(x, y, psi, levels=np.arange(-U*10, U*10, U / 5))

axs[0].contour(x, y, psi, levels=[0, 0.0001], colors='0.5')


# transformed
circleZ = circle + a**2 / circle
axs[1].plot(np.real(circleZ), np.imag(circleZ), color='k')
axs[1].set_aspect(1)
axs[1].contour(x, y, psinew, levels=np.arange(-U*10, U*10, U / 5))
axs[1].contour(x, y, psinew, levels=[0, 0.0001], colors='0.5')

axs[1].set_xlim([-4.1, 4.1])
axs[1].set_ylim([-4.1, 4.1])

axs[1].set_xlabel('$x / a$')
axs[1].set_ylabel('$y / a$')
axs[1].set_title('$z$')

/var/folders/kx/4c3bx_w92t31pqyvg6ws03200000gn/T/ipykernel_25412/491758429.py:11: RuntimeWarning: invalid value encountered in sqrt
  Znew = 0.5*Z + np.sqrt(0.25*np.abs(Z)**2 - a**2)

Text(0.5, 1.0, '$z$')

Now to get the lift. First, we could just short-circuit this and say that the lift is $\rho \Gamma U = \rho 4\pi U a \sin \alpha$ and be done with it by using the lift theorem.

If we want to show that this works, then you integrate the pressure on the top and bottom of the plate:

$$ P = P_{\infty} + U^2/2 - u^2(y=0)/2$$

To do this we note that

$$ U-iV = (u+iv)/ (1-a^2/z^2)$$

If we carry this integral out we note that it is ruined by the sigularity at $\theta=\pi$, and because of this singularity, the lift is infinite.

dt=0.0000001
theta = np.arange(0, np.pi*2, dt)
z = a * np.exp(1j*theta)

Uu = (U*np.exp(-1j*alpha) - U*np.exp(1j*alpha)*a**2/z**2 - 1j*gamma/2/np.pi/z) / (1-a**2/z**2) 
fig, ax = plt.subplots()
ax.plot(theta, np.real(Uu)**2)
print(np.nansum(np.real(Uu)**2*dt))

/Users/jklymak/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in true_divide
  """

46773353.14670776