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 get_y_hh(dy, l, Re_l):
    return dy * np.sqrt(l * Re_l)


# subtask c
angle_max = 45
dx = 0.1
dy = 0.1
U_1 = 20
upsilon = 1
D = upsilon * 10000 / U_1
R = D / 2
U_2 = 25000 / 10000 * U_1
endangle = 45
U_x10, d_x10 = Estimate_BL_D_PFT(U_1, D, upsilon, endangle, 0.1, 0.1)
u10, v10 = implicitsolver(D, upsilon, endangle, dx, dy, U_x10, d_x10, 0)
U_x25, d_x25 = Estimate_BL_D_PFT(U_2, D, upsilon, endangle, 0.1, 0.1)
u25, v25 = implicitsolver(D, upsilon, endangle, dx, dy, U_x25, d_x25, 0)

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

# Plot azimuthal velocity profiles at 45°
plt.subplot(2, 2, 1)
x = np.linspace(0, u10.shape[1] * dy, u10.shape[1])
plt.plot(x[:int(d_x10[int(endangle / dx)] / dy)], u10[int(endangle / dx), :int(d_x10[int(endangle / dx)] / dy)])
x = np.linspace(0, u25.shape[1] * dy, u25.shape[1])
plt.plot(x[:int(d_x25[int(endangle / dx)] / dy)], u25[int(endangle / dx), :int(d_x25[int(endangle / dx)] / dy)])
plt.legend(["Re=10000", "Re=25000"])
plt.xlabel("Distance from wall [SI]")
plt.ylabel("Velocity magnitude [SI]")
plt.title("Tangent velocity profiles in the boundary layer around a cylinder at 45° angle (laminar)")

# Calculate u', y'', v'' for RE=10000 and RE=25000
u_hip10 = u10[int(45 / dx), :] / U_1
u_hip25 = u25[int(45 / dx), :] / U_2
v_hip10 = v10[int(45 / dx), :] / U_1
v_hip25 = v25[int(45 / dx), :] / U_2
l = R * np.deg2rad(45)
Re_l_10 = U_1 * l / upsilon
Re_l_25 = U_2 * l / upsilon
v_hiphip10 = v_hip10 * np.sqrt(Re_l_10)
v_hiphip25 = v_hip25 * np.sqrt(Re_l_25)
y_hh10 = np.arange(0, (u10.shape[1] - 1) * get_y_hh(dy, l, Re_l_10), get_y_hh(dy, l, Re_l_10))
y_hh25 = np.arange(0, (u25.shape[1] - 1) * get_y_hh(dy, l, Re_l_25), get_y_hh(dy, l, Re_l_25))

# Plot u' in function of y''
plt.subplot(2, 2, 2)
plt.plot(y_hh10[:int(d_x10[int(45 / dx)] / dy)], u_hip10[:int(d_x10[int(45 / dx)] / dy)])
plt.plot(y_hh25[:int(d_x25[int(45 / dx)] / dy)], u_hip25[:int(d_x25[int(45 / dx)] / dy)])
plt.legend(["Re=10000", "Re=25000"])
plt.xlabel("Dimensionless distance from wall y'' [SI]")
plt.ylabel("Dimensionless velocity magnitude u' [SI]")
plt.title("u' velocity magnitudes in function of y'' (laminar)")

# Plot v'' in function of y''
plt.subplot(2, 2, 3)
plt.plot(y_hh10[:int(d_x10[int(45 / dx)] / dy)], v_hiphip10[:int(d_x10[int(45 / dx)] / dy)])
plt.plot(y_hh25[:int(d_x25[int(45 / dx)] / dy)], v_hiphip25[:int(d_x25[int(45 / dx)] / dy)])
plt.legend(["Re=10000", "Re=25000"])
plt.xlabel("Dimensionless distance from wall y'' [SI]")
plt.ylabel("Dimensionless velocity magnitude v'' [SI]")
plt.title("v'' velocity magnitudes in function of y'' (laminar)")

# Plot u' scaled in function of y''
plt.subplot(2, 2, 4)
plt.plot(y_hh10[:int(d_x10[int(45 / dx)] / dy)], u_hip10[:int(d_x10[int(45 / dx)] / dy)] / u_hip10[-1] * u_hip25[-1])
plt.plot(y_hh25[:int(d_x25[int(45 / dx)] / dy)], u_hip25[:int(d_x25[int(45 / dx)] / dy)] / u_hip25[-1])
plt.plot(y_hh10[:int(d_x10[int(45 / dx)] / dy)], v_hiphip10[:int(d_x10[int(45 / dx)] / dy)] / v_hiphip10[-1] * v_hiphip25[-1])
plt.plot(y_hh25[:int(d_x25[int(45 / dx)] / dy)], v_hiphip25[:int(d_x25[int(45 / dx)] / dy)] / v_hiphip25[-1])
plt.legend(["u' Re=10000", "u' Re=25000", "v'' Re=10000", "v'' Re=25000"])
plt.xlabel("Dimensionless distance from wall y'' [SI]")
plt.ylabel("Dimensionless velocity magnitude u' v'' [SI]")
plt.title("u' and v'' velocity magnitudes in function of y'' scaled to same U(x) boundary layer velocities (laminar)")

plt.tight_layout()
plt.show()