diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..ed8ebf5 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +__pycache__ \ No newline at end of file diff --git a/src/turborocket/meanline/lossmodels.py b/src/turborocket/meanline/lossmodels.py new file mode 100644 index 0000000..379575a --- /dev/null +++ b/src/turborocket/meanline/lossmodels.py @@ -0,0 +1,150 @@ +""" + +This file contains a series of loss-model utilities functions for mean-line design of turbines + +Each Loss model is encapsulated into an object with a series of functions calculating loss contributions. + +Loss-Models currently implemented: + + - craigcox: Loss Model by Craig and Cox + - kackerOkapuu: Loss Model by Kacker and Okapuu (1982) + - Aungier: Loss Modely by Aungier + +""" + +import numpy as np + +class Aungier: + """ + + The Aungier Loss Model + + Variables: + Y_p: Profile Loss parameters + K_mod: Experience factor (derived from the kackerOkapuu method) + K_inc: Off-design incidence factor + K_m: Mach Number Factor + K_p: Compressibility Factor + K_RE: Reynolds Number Factor + + Y_p1: Profile Loss coefficient for nozzle blades (beta_1 = 0) + Y_p2: Profile Loss coefficent for rotor blades (beta_1 = alpha_2) + + Y_s: Secondary Flow Losses for low aspect ratio + + """ + + def __init__(): + # Intialising core parameters of the loss model + return + + def Y(self): + """ + Total Pressure Loss Coefficient (Y) + + This function solves for the profile loss coefficient of the turbine stage, characterising the increase in total pressure of the turbine exit stage + + Y = (P_t1 - P_t2) / (P_t2 - P_2) + + Variables: + Y_p: The Profile Loss Coefficient + Y_s: The Secondary Flow loss coefficient + Y_tcl: The blade clearnace loss coefficient + Y_te: The trailing edge loss coefficient + Y_ex: The supersonic expansion loss coefficient + Y_sh: The shock loss coefficient + T_lw: The lashing wire loss coefficient (rotors) + """ + + return Y_p + Y_s + Y_tcl + T_te + Y_ex + Y_sh + Y_lw + + def delta_h_par(self): + """ + Parasitic Losses (delta_h_par) + + The parasitic losses accounts for waster work due to losses in disk friction, shear forces and partial admission. + + These are reflected in an increase in total enthalph of the discharge flow + relative to the value produced when the flow passes through the blade row. + + In contrast to the pressure loss coefficient, these losses do not affect total pressure, but do influence stage efficiency. + + Variables: + delta_h_adm: Partial admission losses (rotors) + delta_h_df: Disk friction work (diaphragm-disk rotors) + delta_h_seal: Leakage bypass loss (shourded rotors and rotor balance holes) + delta_h_gap: Clearnace gap windage loss (shrouded blades and nozzles of diaphragm-disk rotors) + delta_h_q: moisture work loss (rotors) + """ + + return delta_h_adm + delta_h_df + delta_h_seal + delta_h_gap + delta_h_q + + def yp(self): + """ + The Profile Loss Coefficient (Y_p) + + This pressure loss coefficient characterises the profile losses of the turbine stage + + Variables: + K_mod: kackerOkapuu Experience factor + K_inc: Correction for off-design incidence effects + K_m: Correction for Mach Number effects + K_p: correction for Compressibility effects + K_re: correction for Reynolds Number effects + Y_p1: profile loss coefficient for nozzle blades (beta_1 = 90) + Y_p2: profile loss coefficients for impulse blades (alpha_2 = beta_1) + + Returns: + _type_: _description_ + """ + + prt = (y_p1 + ((beta_1/alpha_2)**2 ) * (y_p2 - y_p1))*((5*t/c)**(beta_1/alpha_2)) - delta_y_te + + y_p = k_mod * k_inc * k_m * k_p * k_re * prt; + + return y_p + + def F_ar(self): + + if h_c < 2: + f_ar = 0.5 * (2*c/h)**(0.7) + else: + f_ar = (c/h) + + return f_ar + + def Y_s(self): + y_bar_s = 0.0334 * F_ar * (C_l / (s_c))**2 * (np.cos(alpha_2)/np.cos(beta_1))*(np.cos(alpha_2)**2/np.cos(alpha_m)***3) + + y_s = k_s*k_re*(y_bar_s**2 / (1 + 7.5*y_bar_s))**(1/2) + + return y_s + + def Y_sh(self): + + y_bar_sh = 0.8*x_1**2 + x_2**2 + + y_sh = (y_bar_sh**2 / (1 + y_bar_sh**2))**(1/2) + + return y_sh + + def Y_ex(self): + + y_ex = ((M2 - 1)/M2)**2; + + return y_ex + + def Y_te(self): + + y_te = ((t_2/s) * np.sin(beta_g) - t_2)**2 + + return y_te + + def Y_tcl(self): + + if TIP_TYPE = "shrouded": + B = 0.37 + else: + B = 0.47 + + Y_tcl = B*(c/h)*(t_cl / c)**(0.78) * (C_l/(s/c))**2 * (np.cos(alpha_2)**2 /np.cos(alpha_m)**3) \ No newline at end of file diff --git a/src/turborocket/profiling/Supersonic/circular.py b/src/turborocket/profiling/Supersonic/circular.py new file mode 100644 index 0000000..9b479eb --- /dev/null +++ b/src/turborocket/profiling/Supersonic/circular.py @@ -0,0 +1,127 @@ +""" +Circular.py + +This source file encapsulates all the functions required for sizing the circular sections of the tubrine blade, following the free vortex assumption. + +The following equations follows the method for supersonic turbine design presented in NASA TN D-4421. + +Free vortex flow is assumed for the blade flow. + +# noqa: E501 +""" + +import numpy as np + + +def M_star(gamma, M): + """ + Critical Velocity Ratio + + This function computes the critical velocity ratio of the flow as a function of Mach number and specific heat ratio. + + Variables: + gamma: Specific Heat Ratio + M: Mach Number + + # noqa: E501 + """ + + crit_vel_rat = ( + (((gamma + 1) / 2) * M**2) / (1 + ((gamma - 1) / 2) * M**2) + ) ** (1 / 2) + + return crit_vel_rat + + +def prandtl_meyer(gamma, crit_vel_rat): + """ + Prandtl Meyer Angle (v) + + This function computes the Prandtl-Meyer angle v based on the mach number of a flow. + + # noqa: E501 + """ + + v = np.pi / 4 * (((gamma + 1) / (gamma - 1)) ** (1 / 2) - 1) + (1 / 2) * ( + (((gamma + 1) / (gamma - 1)) ** (1 / 2)) + * np.arcsin((gamma - 1) * crit_vel_rat**2 - gamma) + + np.arcsin((gamma + 1) / crit_vel_rat**2 - gamma) + ) + + return v + + +def arc_angles_upper(beta_o, beta_i, v_i, v_o, v_u): + """ + Upper Circular Arc Angles Alpha + + This function computes the circular arc angles based on the inlet and outlet angles, + coupled with their prantl-meyer angles + + Variables: + alpha_l_i: Lower Arc Inlet Angle + alpha_l_o: Lower Arc Outlet Angle + alpha_u_i: Upper Arc Inlet Angle + alpha_u_o: Upper Arc Outlet Angle + beta_i: Inlet relative flow angle + beta_o: Exit relative + + # noqa: E501 + """ + + alpha_u_i = beta_i - (v_u - v_i) + + alpha_u_o = beta_o + (v_u - v_o) + + return [alpha_u_i, alpha_u_o] + + +def arc_angles_lower(beta_o, beta_i, v_i, v_o, v_l): + """ + Lower Circular Arc Angles Alpha + + This function computes the arc angles based on the inlet and outlet angles, coupled with their prantl-meyer angles + + Variables: + alpha_l_i: Lower Arc Inlet Angle + alpha_l_o: Lower Arc Outlet Angle + alpha_u_i: Upper Arc Inlet Angle + alpha_u_o: Upper Arc Outlet Angle + beta_i: Inlet relative flow angle + beta_o: Exit relative + + # noqa: E501 + """ + + alpha_l_i = beta_i - (v_i - v_l) + + alpha_l_o = beta_o + (v_o - v_l) + + return [alpha_l_i, alpha_l_o] + + +def beta_o(M_i, M_o, gamma, beta_i): + """ """ + + exit_o = -np.arccos( + ( + (M_i / M_o) + * ((1 + (gamma - 1) / 2 * M_o**2) / (1 + (gamma - 1) / 2 * M_i**2)) + ** ((gamma + 1) / (2 * (gamma - 1))) + ) + * np.cos(beta_i) + ) + + return exit_o + + +def A_rat(beta_i, beta_o): + """ + The Area Ratio (A_i / A_o) + + This equation computes the area ratio of the supersonic turbine blade, assuming that the spacings are equal on the inlet and exit sides of the turbine (Axial turbine) + + # noqa: E501 + """ + + return np.cos(beta_i) / np.cos(beta_o) diff --git a/src/turborocket/profiling/Supersonic/transition.py b/src/turborocket/profiling/Supersonic/transition.py new file mode 100644 index 0000000..255b8ea --- /dev/null +++ b/src/turborocket/profiling/Supersonic/transition.py @@ -0,0 +1,225 @@ +""" +transition.py + +This source file encapsulates all the functions required for sizing the transition arc sections of the turbine blade, utilising the method of characteristics. + +The following equations follows the method for supersonic turbine design presented in NASA TN D-4421. + + + +# noqa: E501 + +""" + +from src.turborocket.solvers.solver import adjoint + +import numpy as np + + +def func_R_star(R_star, gamma): + """ + Function of R* ( f(R*) ) + + This equation describes the generalised function of the non-dimentionalised radius + defining the expansion line angle. + + Variables: + gamma: Specific Heat Ratio + R_star: Non-dimentionalised Radius + + """ + + f_r_star = ((gamma + 1) / (gamma - 1)) ** (1 / 2) * np.arcsin( + (gamma - 1) / R_star**2 - gamma + ) + np.arcsin((gamma + 1) * R_star**2 - gamma) + + return f_r_star + + +def func_R_star_k(gamma, v_i, dv, k): + f_r_star_k = ( + 2 * v_i + - np.pi / 2 * (((gamma + 1) / (gamma - 1)) ** (1 / 2) - 1) + - 2 * (k - 1) * dv + ) + + return f_r_star_k + + +def R_star(gamma, v_i, dv, k, guess): + # We use the ajoint equation + + target = func_R_star_k(gamma, v_i, dv, k) + + r_star = adjoint(func_R_star, guess, -0.01, 50, 0.1, target, params=[gamma]) + + return r_star + + +def phi_k(v_i, v_l, k, dv): + phi_k_i = v_i - v_l - (k - 1) * dv + + return phi_k_i + + +def vortex_coords(R_star_k, phi_k): + x_star = -R_star_k * np.sin(phi_k) + y_star = R_star_k * np.cos(phi_k) + + return [x_star, y_star] + + +def mue_k(gamma, R_star_k): + mach_angle = -np.arcsin( + (((gamma + 1) / 2) * R_star_k**2 - ((gamma - 1) / 2)) ** (1 / 2) + ) + + return mach_angle + + +def mach_slope(phi_k_1, phi_k_2, mue_k_1, mue_k_2): + # This function computes the gradient of the mach line + + m_k = np.tan(((phi_k_1 + phi_k_2) / 2) + ((mue_k_1 + mue_k_2) / 2)) + + return m_k + + +def wall_slope(phi_k_1): + # This function computes the slope the wall at a given segment + + m_bar_k = np.tan(phi_k_1) + + return m_bar_k + + +def wall_coords(x_star_l_k_1, y_star_l_k_1, y_star_k, x_star_k, m_bar_k, m_k): + x_star_l_k = ( + (y_star_l_k_1 - m_bar_k * x_star_l_k_1) - (y_star_k - m_k * x_star_k) + ) / (m_k - m_bar_k) + y_star_l_k = ( + m_k * (y_star_l_k_1 - m_bar_k * x_star_l_k_1) + - m_bar_k * (y_star_k - m_k * x_star_k) + ) / (m_k - m_bar_k) + + return [x_star_l_k, y_star_l_k] + + +def transform_coord(x_star_l_k, y_star_l_k, alpha_l_i): + # This function does the co-ordinate transformation for the transition co-ords + + x_star_l_k_t = x_star_l_k * np.cos(alpha_l_i) - y_star_l_k * np.sin(alpha_l_i) + + y_star_l_k_t = x_star_l_k * np.sin(alpha_l_i) + y_star_l_k * np.cos(alpha_l_i) + + return [x_star_l_k_t, y_star_l_k_t] + + +def moc(k_max, v_i, v_l, gamma, alpha_l_i): + # This function encapsulates the overall method of characteristics procedure for solving + # for the shape of the transition arcs + + # First we need to define the k_max value, this needs to be an integer, + # so what we do is we can invert this + + delta_v = (v_i - v_l) / k_max + + # We setup our history arrays accordingly + phi_hist = [] + r_star_hist = [] + xk_hist = [] + yk_hist = [] + mue_hist = [] + mk_hist = [] + m_bar_k_hist = [] + xlk_hist = [] + ylk_hist = [] + xlkt_hist = [] + ylkt_hist = [] + + # We can now do the intial look to get the intial values of the loop. + phi = phi_k(v_i, v_l, k_max + 1, delta_v) + + phi_hist.append(phi) + + # We now must solver R* using the adjoint method + r_star = R_star(gamma, v_i, delta_v, k_max + 1, 1) + + r_star_hist.append(r_star) + + # We can now compute the co-ordinates of the expansion line + [xk, yk] = vortex_coords(r_star, phi) + + xk_hist.append(xk) + yk_hist.append(yk) + + # We need to now compute the mach angle + mue = mue_k(gamma, r_star) + + mue_hist.append(mue) + + # Gradients need to be 0 at the intial point which is vertical + mk_hist.append(0) + m_bar_k_hist.append(0) + + # Wall co-ordinates are the same as the intial co-ordinates + xlk_hist.append(xk) + ylk_hist.append(yk) + + # Finally, we do a co-ordinate transformation + + [xlkt, ylkt] = transform_coord(xk, yk, alpha_l_i) + + xlkt_hist.append(xlkt) + ylkt_hist.append(ylkt) + + # Having established the delta_v, we can then loop through k values + + for k in np.linspace(k_max, 1, num=k_max): + # We evaluate our phi angle + phi = phi_k(v_i, v_l, k, delta_v) + + phi_hist.append(phi) + + # We now must solver R* using the adjoint method + r_star = R_star(gamma, v_i, delta_v, k, 1) + + r_star_hist.append(r_star) + + # We can now compute the co-ordinates of the expansion line + [xk, yk] = vortex_coords(r_star, phi) + + xk_hist.append(xk) + yk_hist.append(yk) + + # We need to now compute the mach angle + mue = mue_k(gamma, r_star) + + mue_hist.append(mue) + + # We must thus compute the slope of the mach line, based on the average of the current point and previous (k+1) + m_k = mach_slope(phi_hist[k_max - k], phi, mue_hist[k_max - k], mue) + + mk_hist.append(m_k) + + # We can now compute the wall segment slope + m_bar_k = wall_slope(phi_hist[k_max - k]) + + m_bar_k_hist.append(m_bar_k) + + # Wall Co-ordinates + [xlk, ylk] = wall_coords(xlk_hist, ylk_hist, yk_hist, xk_hist, m_bar_k, m_k) + + xlk_hist.append(xlk) + ylk_hist.append(ylk) + + # Finally, we do a co-ordinate transformation + + [xlkt, ylkt] = transform_coord(xlk, ylk, alpha_l_i) + + xlkt_hist.append(xlkt) + ylkt_hist.append(ylkt) + + # Finally, we return the transformed co-ordinates + + return [xlkt_hist, ylkt_hist] diff --git a/src/turborocket/profiling/supersonic_profile.py b/src/turborocket/profiling/supersonic_profile.py new file mode 100644 index 0000000..e912984 --- /dev/null +++ b/src/turborocket/profiling/supersonic_profile.py @@ -0,0 +1,7 @@ +# This file encapsulates the main function for computing the supersonic profile + + +def supersonic_profile(beta_I, beta_o, M_i, M_o): + + # The first process required for the design of the supersonic blade profile is the definition of the upper and lower circular arcs. + diff --git a/src/turborocket/solvers/solver.py b/src/turborocket/solvers/solver.py new file mode 100644 index 0000000..85304e4 --- /dev/null +++ b/src/turborocket/solvers/solver.py @@ -0,0 +1,62 @@ +import numpy as np +import matplotlib.pyplot as plt + + +def adjoint(func, x_guess, dx, n, relax, target, params=[], RECORD_HIST=False): + # We first evaluate the error as a function of x + + x_hist = [] + guess_hist = [] + error_hist = [] + error_grad_hist = [] + dx_hist = [] + + # We first evaluate the error of the guess + guess = func(x_guess, *params) + + error = guess - target + + guess_hist.append(guess) + error_hist.append(error) + x_hist.append(x_guess) + error_grad_hist.append(0) + + for k in np.linspace(1, n, num=n, dtype=int): + # We now timestep in the x space + guess = func(x_hist[k - 1] + dx, *params) + + error = guess - target + + guess_hist.append(guess) + error_hist.append(error) + x_hist.append(x_hist[k - 1] + dx) + dx_hist.append(dx) + # Now we computed the gradient + + error_grad = (error_hist[k] - error_hist[k - 1]) / (x_hist[k] - x_hist[k - 1]) + + error_grad_hist.append(error_grad) + if abs(error) < 1e-5: + break + + # based on this gradient, we can thus find the new x-step based on this gradient # noqa: E501 + + dx = -(1 / error_grad) * error_hist[k] * relax + + if RECORD_HIST is True: + plt.plot(x_hist, error_hist) + + plt.show() + + print("Guess History") + print(guess_hist) + print("Error History") + print(error_hist) + print("x History") + print(x_hist) + print("Error Grad History") + print(error_grad_hist) + print("Dx hist") + print(dx_hist) + + return x_hist[k] diff --git a/tests/profiling/supersonic/test_circular.py b/tests/profiling/supersonic/test_circular.py new file mode 100644 index 0000000..21c7325 --- /dev/null +++ b/tests/profiling/supersonic/test_circular.py @@ -0,0 +1,46 @@ +from src.turborocket.profiling.Supersonic.circular import ( + M_star, + prandtl_meyer, + arc_angles_lower, + arc_angles_upper, + beta_o, + A_rat, +) + +import numpy as np + + +def test_M_star(): + assert np.isclose(M_star(1.4, 2), 1.63299316, rtol=1e-3) + + +def test_prandtl_meyer(): + assert np.isclose(prandtl_meyer(1.4, 1.63), 0.457243, rtol=1e-3) + + +def test_arc_angles_lower(): + alpha_l_i, alpha_l_o = arc_angles_lower( + -np.pi / 3, np.pi / 3, np.pi / 3, -np.pi / 3, np.pi / 2 + ) + + assert np.isclose(alpha_l_i, 1.570796327, rtol=1e-3) + + assert np.isclose(alpha_l_o, -3.665191429, rtol=1e-3) + + +def test_arc_angles_upper(): + alpha_u_i, alpha_u_o = arc_angles_upper( + -np.pi / 3, np.pi / 3, np.pi / 3, -np.pi / 3, np.pi / 2 + ) + + assert np.isclose(alpha_u_i, 0.5235987759, rtol=1e-3) + + assert np.isclose(alpha_u_o, 1.570796327, rtol=1e-3) + + +def test_beta_o(): + assert np.isclose(beta_o(2, 1.5, 1.4, np.pi / 4), -1.055440078, rtol=1e-3) + + +def test_A_rat(): + assert np.isclose(A_rat(np.pi / 4, -np.pi / 4), 1, rtol=1e-3) diff --git a/tests/profiling/supersonic/test_transition.py b/tests/profiling/supersonic/test_transition.py new file mode 100644 index 0000000..eecfc72 --- /dev/null +++ b/tests/profiling/supersonic/test_transition.py @@ -0,0 +1,69 @@ +from src.turborocket.profiling.Supersonic.transition import ( + func_R_star, + func_R_star_k, + R_star, + phi_k, + vortex_coords, + mue_k, + mach_slope, + wall_slope, + wall_coords, + transform_coord, + moc, +) + +import numpy as np + + +def test_func_R_star(): + assert np.isclose(func_R_star(0.9967412602371031, 1.4), -2.276853164, rtol=1e-3) + + +def test_func_R_star_k(): + assert np.isclose(func_R_star_k(1.4, np.pi, 2, 3), -3.993667857, rtol=1e-3) + + +def test_R_star(): + assert np.isclose(R_star(1.4, np.pi, 1.5708, 3, 1), 1, rtol=1e-3) + + +def test_phi_k(): + assert np.isclose(phi_k(np.pi / 2, np.pi / 4, 3, 1.57), -2.354601837, rtol=1e-3) + + +def test_vortex_coords(): + xk, yk = vortex_coords(0.8, np.pi / 2) + + assert np.isclose(xk, -4 / 5, rtol=1e-3) + + assert np.isclose(yk, 0, rtol=1e-3) + + +def test_mue_k(): + assert np.isclose(mue_k(1.4, 0.8), -0.8536095489, rtol=1e-3) + + +def test_mach_slope(): + assert np.isclose( + mach_slope(np.pi / 2, np.pi / 4, np.pi / 3, np.pi / 2), -0.767326988, rtol=1e-3 + ) + + +def test_wall_slope(): + assert np.isclose(wall_slope(np.pi / 4), 1, rtol=1e-3) + + +def test_wall_coords(): + xkl, ykl = wall_coords(1, 1, 0.8, 0.8, -0.76, 1) + + assert np.isclose(xkl, 1, rtol=1e-3) + + assert np.isclose(ykl, 1, rtol=1e-3) + + +def test_transform_coord(): + xklt, yklt = transform_coord(2, 3, np.pi / 3) + + assert np.isclose(xklt, -1.598076, rtol=1e-3) + + assert np.isclose(yklt, 3.232050808, rtol=1e-3) diff --git a/tests/solvers/test_solver.py b/tests/solvers/test_solver.py new file mode 100644 index 0000000..f87003b --- /dev/null +++ b/tests/solvers/test_solver.py @@ -0,0 +1,11 @@ +from src.turborocket.solvers.solver import adjoint + +import pytest + + +def dummy_function(x): + return x**2 - 4 + + +def test_adjoint(): + assert adjoint(dummy_function, 4, 0.1, 100, 0.6, 0) == pytest.approx(2, 1e-3)