Untitled
unknown
python
a year ago
4.2 kB
8
Indexable
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()
Editor is loading...
Leave a Comment