Untitled
unknown
python
a year ago
5.3 kB
6
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] 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()
Editor is loading...
Leave a Comment