Untitled

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags


def Estimate_BL_D_PFT(U_1, D, upsilon, angle_max, dx, dy):
    theta = np.arange(0, angle_max + dy, dy)
    beta = np.arccos(np.cos(theta) * np.sqrt(U_1 / D))
    delta_theta = np.diff(theta)
    delta_beta = np.diff(beta)
    delta_x = delta_theta * D / np.tan(beta)
    delta_y = delta_theta * D / np.sin(beta)
    d_x = np.cumsum(delta_x)
    U_x = U_1 * (1 - np.sqrt(np.cos(theta[1:])))
    return U_x, d_x


def implicitsolver(D, upsilon, angle_max, dx, dy, U_x, d_x, plot_interval):
    angle_steps = int(angle_max / dy) + 1
    max_steps = 10000
    epsilon = 1e-5
    u = np.zeros((angle_steps, max_steps))
    v = np.zeros((angle_steps, max_steps))
    omega = np.zeros((angle_steps, max_steps))
    A = diags([-1, 2, -1], [-1, 0, 1], shape=(angle_steps, angle_steps)).toarray() / dy ** 2
    for i in range(1, max_steps):
        u[:, i] = u[:, i - 1]
        v[:, i] = v[:, i - 1]
        omega[:, i] = omega[:, i - 1]
        for j in range(1, angle_steps - 1):
            u[j, i] = (2 * u[j, i - 1] + omega[j, i] * dx ** 2 * (v[j - 1, i] - v[j + 1, i]) +
                       upsilon * dx ** 2 * (u[j - 1, i] + u[j + 1, i]) +
                       upsilon * dy ** 2 * (u[j, i - 1] - u_x[j]) + D * omega[j, i] * dx ** 2) / (2 * (upsilon * dx ** 2 + D * dx ** 2))
            v[j, i] = (v[j - 1, i] + v[j + 1, i] +
                       upsilon * dy ** 2 * (u[j, i] - u[j, i - 1]) + D * omega[j, i] * dy ** 2) / (2 * (upsilon * dy ** 2 + D * dy ** 2))
            omega[j, i] = (omega[j - 1, i] + omega[j + 1, i]) / 2
        if np.linalg.norm(u[:, i] - u[:, i - 1]) < epsilon:
            break
    return u[:, :i + 1], v[:, :i + 1]


def makedatatip(h_plot, index):
    def on_pick(event):
        ind = event.ind[0]
        xpos = event.artist.get_xdata()[ind]
        ypos = event.artist.get_ydata()[ind]
        print(f'Index: {ind}, X: {xpos}, Y: {ypos}')

    fig, ax = plt.subplots()
    ax.plot(h_plot)
    ax.set_xlabel('Angle [°]')
    ax.set_ylabel('Derivative of tangential velocity near wall [SI]')
    ax.set_title('Point of separation (Re=10000, laminar)')
    fig.canvas.mpl_connect('pick_event', on_pick)
    plt.show()


# subtask b
angle_max = 100.2
dx = 0.1
dy = 0.1
U_1 = 20
upsilon = 1
D = upsilon * 10000 / U_1
U_x, d_x = Estimate_BL_D_PFT(U_1, D, upsilon, angle_max, dx, dy)
u, v = implicitsolver(D, upsilon, angle_max, dx, dy, U_x, d_x, 0)

# Plotting
plt.figure(figsize=(12, 10))

# Plot azimuthal velocity curves along the boundary layer
plt.subplot(2, 2, 1)
curvenum = 8
s = int(np.floor(u.shape[0] / curvenum))
x = np.linspace(0, u.shape[1] * dy, u.shape[1])
l = []
for i in range(s, s * curvenum + 1, s):
    plt.plot(x[:int(d_x[i] / dy)], u[i, :int(d_x[i] / dy)])
    l.append(f'{i * dy:.1f} [SI] distance along the wall')
plt.legend(l)
plt.xlabel('Distance from wall [SI]')
plt.ylabel('Velocity magnitude [SI]')
plt.title('Tangent velocity profiles as the boundary layer evolves (Re=10000, laminar)')

# Plot azimuthal velocity curves along the boundary layer, zoomed in
plt.subplot(2, 2, 3)
curvenum = 10
startangle = 99
endangle = 100.2
s = int((endangle / dy - startangle / dy) / curvenum)
x = np.linspace(0, u.shape[1] * dy, u.shape[1])
l = []
for i in range(startangle, startangle + s * curvenum + 1, s):
    plt.plot(x, u[i, :])
    l.append(f'{(i - 3) * dy:.1f} [SI] distance along the wall')
plt.plot(x, u[-1, :])
l.append(f'{endangle} [SI] distance along the wall')
plt.legend(l)
plt.xlim([0, 10])
plt.xlabel('Distance from wall [SI]')
plt.ylabel('Velocity magnitude [SI]')
plt.title('Tangent velocity profiles as the boundary layer evolves (Re=10000, laminar)')

# Plot normal velocity components along the boundary layer
plt.subplot(2, 2, 2)
curvenum = 8
s = int(np.floor(u.shape[0] / curvenum))
x = np.linspace(0, u.shape[1] * dy, u.shape[1])
l = []
for i in range(s, s * curvenum + 1, s):
    plt.plot(x[:int(d_x[i] / dy)], v[i, :int(d_x[i] / dy)])
    l.append(f'{i * dy:.1f} [SI] distance along the wall')
plt.legend(l)
plt.xlabel('Distance from wall [SI]')
plt.ylabel('Velocity magnitude [SI]')
plt.title('Normal velocity profiles as the boundary layer evolves (Re=10000, laminar)')

# Calculate derivative of the tangential velocity at the wall
derivative = np.zeros(u.shape[0])
for i in range(1, u.shape[0]):
    derivative[i] = (u[i, 1] - u[i, 0]) / dx

# Plot the derivative of the tangential velocity at the wall
plt.subplot(2, 2, 4)
y = np.linspace(0, u.shape[0] * dx, u.shape[0])
hplot, = plt.plot(y, derivative)
plt.xlabel('Angle [°]')
plt.ylabel('Derivative of tangential velocity near wall [SI]')
plt.title('Point of separation (Re=10000, laminar)')
plt.tight_layout()
plt.show()

# Add data tip
makedatatip(hplot, 1000)