# Untitled

unknown
python
14 days ago
4.9 kB
1
Indexable
Never
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()

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()