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]


# subtask b
angle_max = 95
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, 0.1, 0.1)
u, v = implicitsolver(D, upsilon, angle_max, dx, dy, U_x, d_x, 1)

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

# Plot azimuthal velocity curves along the boundary layer
plt.subplot(2, 2, 1)
curvenum = 8
s = np.floor(u.shape[0] / curvenum).astype(int)
x = np.linspace(0, u.shape[1] * dy, u.shape[1])
l = [f"{i * dy}[SI] distance along the wall" for i in range(s, s * curvenum, s)]
for i in range(s, s * curvenum, s):
    plt.plot(x[:int(d_x[i] / dy)], u[i, :int(d_x[i] / dy)])
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, turbulent)")

# Plot azimuthal velocity curves along the boundary layer (zoomed in)
plt.subplot(2, 2, 3)
curvenum = 10
startangle = 90
endangle = 95
s = int((endangle / dy - startangle / dy) / curvenum)
x = np.linspace(0, u.shape[1] * dy, u.shape[1])
l = [f"{(i - 3) * dy}[SI] distance along the wall" for i in range(startangle, endangle + s, s)]
for i in range(startangle, endangle + s, s):
    plt.plot(x, u[i, :])
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, turbulent)")

# Plot normal velocity components along the boundary layer
plt.subplot(2, 2, 2)
curvenum = 8
s = np.floor(u.shape[0] / curvenum).astype(int)
x = np.linspace(0, u.shape[1] * dy, u.shape[1])
l = [f"{i * dy}[SI] distance along the wall" for i in range(s, s * curvenum, s)]
for i in range(s, s * curvenum, s):
    plt.plot(x[:int(d_x[i] / dy)], v[i, :int(d_x[i] / dy)])
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, turbulent)")

# 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])
derivative = np.zeros(u.shape[0])
for i in range(1, u.shape[0]):
    derivative[i] = (u[i, 1] - u[i, 0]) / dx
plt.plot(y, derivative)
plt.xlabel("Angle [°]")
plt.ylabel("Derivative of tangential velocity near wall [SI]")
plt.title("Point of separation (Re=10000, turbulent)")

plt.show()