"""Quasi-one-dimensional flow calculator (nozzles and diffusers)"""
import os
from typing import Union, Any, Tuple
if os.environ.get('__USE_CUDA', '0') == '1':
import cupy as np
else:
import numpy as np
from .calculator import G3DCalculator
from .expansion_wave import G3DExpansionWaveCalculator
from .normal_shock import G3DNormalShockCalculator
from .oblique_shock import G3DObliqueShockCalculator
from g3d.defs import CFlowRegime
from g3d.defs import CFlowVariableType
[docs]
class G3DNozzleCalculator(G3DCalculator):
"""
class collecting functions related to nozzle flows
(Q1D one-dimensional flows with smooth cross-section area changes)
"""
def __init__(self, *args: Any, **kwargs: Any) -> None:
"""Initialize calculator parameters"""
super().__init__(*args, **kwargs)
def _set_default(self) -> None:
super()._set_default()
self.__expansion = G3DExpansionWaveCalculator(update_code_type=False)
self.__normal_shock = G3DNormalShockCalculator(update_code_type=False)
self.__oblique_shock = G3DObliqueShockCalculator(update_code_type=False)
self.__M = self._vars.add(name="M", long_name="Mach number")
self.__AR = self._vars.add(name="AR", long_name="Area ratio")
self.__Po = self._vars.add(name="Po", long_name="Total pressure", type=CFlowVariableType.PRESSURE)
self.__To = self._vars.add(name="To", long_name="Total temperature", type=CFlowVariableType.TEMPERATURE)
self.__At = self._vars.add(name="At", long_name="Throat area", type=CFlowVariableType.AREA)
self.__mp = self._vars.add(name="mp", long_name="Massflow", type=CFlowVariableType.MASSFLOW)
self.__Pb = self._vars.add(name="Pb", long_name="Back pressure", type=CFlowVariableType.PRESSURE)
self.__Ae = self._vars.add(name="Ae", long_name="Nozzle-exit area", type=CFlowVariableType.AREA)
self.__M1 = self._vars.add(name="M1", long_name="Mach number upstream of shock")
self.__Pe = self._vars.add(name="Pe", long_name="Nozzle-exit pressure", type=CFlowVariableType.PRESSURE)
def _associate(self) -> None:
super()._associate()
self.__normal_shock.gas = self._gas
[docs]
def area_ratio(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
.. deprecated:: 2.6
redirects to :func:`get_area_ratio`
"""
return self.get_area_ratio(*args, **kwargs)
[docs]
def get_area_ratio(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
calculates area ratio (AR []) for specified Mach number (M []) using the area-Mach-number
relation
input: Mach number (M [])
example: AR = get_area_ratio(M=M)
"""
return self._eval(["M"], self.__get_area_ratio, *args, **kwargs)
def __get_area_ratio(self, M: float = None) -> float:
if M is None:
M = self.__M.value
return (np.sqrt(np.power(((2. / self._gas.gp1) * self._stagnation.get_total_temperature(M)),
(self._gas.gp1 / self._gas.gm1)) / (M * M)))
[docs]
def find_Mach_from_area_ratio(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
.. deprecated:: 2.6
redirects to :func:`get_Mach_from_area_ratio`
"""
return self.get_Mach_from_area_ratio(*args, **kwargs)
[docs]
def get_Mach_from_area_ratio(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
calculates Mach number (M []) corresponding to the specified area ratio (AR [])
the Mach number is calculated using the built-in Newton-Raphson solver
the flow is assumed to be subsonic unless mode='supersonic' is specified
input: area ratio (AR [])
additional input: mode = 'subsonic' (default if not specified) or mode = 'supersonic'
example: M = get_Mach_from_area_ratio(AR=AR)
"""
self._reset_eval(mode=CFlowRegime.SUBSONIC)
self._newton.reset()
return self._eval(["AR"], self.__get_Mach_from_area_ratio, *args, **kwargs)
def __get_Mach_from_area_ratio(self, AR: float = None) -> float:
if AR is None:
AR = self.__AR.value
newton = self._get_newton_kernel()
if self._mode == CFlowRegime.SUPERSONIC:
newton.start_value = 1.5
newton.lower_bound = 1.
else:
newton.start_value = 0.02
newton.upper_bound = 1.
# run Newton solver to find Mach number
return newton.solve(self.__get_area_ratio, target_value=AR)
[docs]
def get_choked_mass_flow(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
calculates choked massflow (mp [kg/s]) for specified upstream total pressure (Po [Pa]),
total temperature (To [K]), and throat area (At [m^2])
input: total pressure (Po [Pa]), total temperature (To [K]), throat area (At [m^2])
example: me = get_choked_mass_flow(Po=Po, To=To, At=At)
"""
return self._eval(["Po", "At", "To"], self.__get_choked_mass_flow, *args, **kwargs)
def __get_choked_mass_flow(self, Po: float = None, At: float = None, To: float = None) -> float:
if Po is None:
Po = self.__Po.value
if To is None:
To = self.__To.value
if At is None:
At = self.__At.value
return (Po * At / np.sqrt(To) * np.sqrt((self._gas.gamma / self._gas.R)
* np.power((2. / self._gas.gp1),
self._gas.gp1 / self._gas.gm1)))
[docs]
def get_choked_area(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
calculates choked cross-section area (At [m^2]) for specified choked massflow (mp [kg/s]),
total pressure (Po [Pa]), and total temperature (To [K])
input: massflow (mp [kg/s]), total pressure (Po [Pa]), total temperature (To [K])
example: As = get_choked_area(mp=mp, Po=Po, To=To)
"""
return self._eval(["mp", "Po", "To"], self.__get_choked_area, *args, **kwargs)
def __get_choked_area(self, mp: float = None, Po: float = None, To: float = None) -> float:
if mp is None:
mp = self.__mp.value
if Po is None:
Po = self.__Po.value
if To is None:
To = self.__To.value
return mp / self.__get_choked_mass_flow(Po, 1.0, To)
[docs]
def internal_shock_get_exit_Mach_number(self, *args: Any, **kwargs: Any) \
-> Union[np.float64, list[np.float64], float]:
"""
calculates exit Mach number (Me []) for nozzle flow with internal shock
input: total pressure upstream of shock (Po [Pa], same as nozzle inlet total pressure),
static pressure downstream of nozzle exit (Pb [Pa]), nozzle exit-to-throat area ratio (AR=Ae/At [])
example: Me = internal_shock_get_exit_Mach_number(Po=Po, Pb=Pb, AR=AR)
"""
return self._eval(["Po", "Pb", "AR"], self.__internal_shock_get_exit_Mach_number, *args, **kwargs)
def __internal_shock_get_exit_Mach_number(self, Po: float = None, Pb: float = None, AR: float = None) -> float:
if Po is None:
Po = self.__Po.value
if Pb is None:
Pb = self.__Pb.value
if AR is None:
AR = self.__AR.value
x = Po / Pb / AR
y = (2. / self._gas.gp1) ** (self._gas.gp1 / self._gas.gm1)
return np.sqrt(-(1. / self._gas.gm1) + np.sqrt((1. / self._gas.gm1) ** 2 + (2. / self._gas.gm1) * y * x ** 2))
[docs]
def internal_shock_get_exit_conditions(self, *args: Any, **kwargs: Any) \
-> Union[Tuple[np.float64, ...], list[Tuple[np.float64, ...]]]:
"""
calculates exit Mach number and pressure (Me [], Pe [Pa]) for nozzle flow with internal shock
some values are calculated using the built-in Newton-Raphson solver
input: nozzle inlet Mach number (M1 []), nozzle inlet total pressure (Po [Pa]), throat area (At [m^2]),
nozzle exit area (Ae [m^2])
example: Me, Pe = internal_shock_get_exit_conditions(M1=M1, Po=Po, At=At, Ae=Ae)
"""
self._reset_eval(mode=CFlowRegime.SUBSONIC)
self._newton.reset()
return self._eval(["M1", "Po", "At", "Ae"],
self.__internal_shock_get_exit_conditions, *args, **kwargs, outdim=2)
def __internal_shock_get_exit_conditions(self, M1: float = None, Po: float = None,
At: float = None, Ae: float = None) -> tuple[float, float]:
if M1 is None:
M1 = self.__M1.value
if Po is None:
Po = self.__Po.value
if At is None:
At = self.__At.value
if Ae is None:
Ae = self.__Ae.value
As = self.__normal_shock.get_downstream_Astar(M1=M1, As1=At)
Po2 = self.__normal_shock.get_downstream_total_pressure(M1=M1, Po1=Po)
AR = Ae / As
Me = self.__get_Mach_from_area_ratio(AR)
Pe = self._stagnation.get_static_pressure(M=Me, Po=Po2)
return Me, Pe
[docs]
def internal_shock_get_exit_pressure(self, *args: Any, **kwargs: Any) \
-> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
calculates exit pressure (Pe [Pa]) for nozzle flow with internal shock
calculations are done using the built-in Newton-Raphson solver
input: nozzle inlet Mach number (M1 []), nozzle inlet total pressure (Po [Pa]), throat area (At [m^2]),
nozzle exit area (Ae [m^2])
example: Pe = internal_shock_get_exit_pressure(M1=M1, Po=Po, At=At, Ae=Ae)
"""
self._reset_eval(mode=CFlowRegime.SUBSONIC)
self._newton.reset()
return self._eval(["M1", "Po", "At", "Ae"], self.__internal_shock_get_exit_pressure, *args, **kwargs)
def __internal_shock_get_exit_pressure(self, M1: float = None) -> float:
Me, Pe = self.__internal_shock_get_exit_conditions(M1)
return Pe
[docs]
def find_Mach_internal_shock(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
.. deprecated:: 2.6
redirects to :func:`get_Mach_internal_shock`
"""
return self.get_Mach_internal_shock(*args, **kwargs)
[docs]
def get_Mach_internal_shock(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
calculate Mach number (M1 []) just upstream of shock for choked nozzle flow with internal shock
calculations are made using the built-in Newton-Raphson solver
input: inlet total pressure (Po [Pa], total pressure upstream of shock), throat area (At [m^2]),
exit area (Ae [m^2]), nozzle exit pressure (Pe [Pa])
example: M1 = get_Mach_internal_shock(Po=Po, At=At, Ae=Ae, Pe=Pe)
"""
self._reset_eval(mode=CFlowRegime.SUPERSONIC)
self._newton.reset()
return self._eval(["Po", "At", "Ae", "Pe"], self.__get_Mach_internal_shock, *args, **kwargs)
def __get_Mach_internal_shock(self, At: float = None, Ae: float = None, Pe: float = None) -> float:
if At is None:
At = self.__At.value
if Ae is None:
Ae = self.__Ae.value
if Pe is None:
Pe = self.__Pe.value
AR = Ae / At
# supersonic flow ahead of internal shock
self._mode = CFlowRegime.SUPERSONIC
M = 0.5 * (1.0 + self.__get_Mach_from_area_ratio(AR))
newton = self._get_newton_kernel()
# subsonic flow at exit downstream of internal shock
self._mode = CFlowRegime.SUBSONIC
return newton.solve(self.__internal_shock_get_exit_pressure, target_value=Pe, start_value=M)
# ================================================================
# METHOD OF CHARACTERISTICS
# ================================================================
[docs]
def moc(self, Me: float, n: int = 10, t0: float = 0.01) -> Tuple[np.ndarray, ...]:
"""
simple Methods-Of-Characteristics (MOC) deign of 2D planar divergent nozzle shape
input: exit Mach number (design condition/supercritical flow) (Me []),
number of characteristic lines (n, default value: 10),
start angle (t0 [rad], default value: 0.01)
example: Km, Kp, t, v, M, mu, x, y = moc(Me=Me, n=n, t0=t0)
"""
t_max = self.__expansion.get_Prandtl_Meyer_function(Me) / 2.
t0 *= np.pi / 180.
dt = (t_max - t0) / (n - 1)
NN = int(0.5 * n * (3 + n))
t = np.zeros(NN)
v = np.zeros(NN)
Kp = np.zeros(NN)
Km = np.zeros(NN)
M = np.zeros(NN)
mu = np.zeros(NN)
x = np.zeros(NN)
y = np.zeros(NN)
for i in range(n):
t[i] = t0 + dt * i
v[i] = t[i]
Km[i] = t[i] + v[i]
Kp[i] = t[i] - v[i]
M[i] = self.__expansion.get_Mach(v[i])
mu[i] = self.__oblique_shock.get_Mach_angle(M[i])
i0 = 0
for k in range(n - 1):
i0 += n - k + 1
Km[i0] = Km[i0 - n + k]
t[i0] = 0.
v[i0] = Km[i0]
Kp[i0] = t[i0] - v[i0]
M[i0] = self.__expansion.get_Mach(v[i0])
mu[i0] = self.__oblique_shock.get_Mach_angle(M[i0])
for j in range(n - 2 - k):
i = i0 + 1 + j
Km[i] = Km[i - n + k]
Kp[i] = Kp[i - 1]
t[i] = 0.5 * (Km[i] + Kp[i])
v[i] = 0.5 * (Km[i] - Kp[i])
M[i] = self.__expansion.get_Mach(v[i])
mu[i] = self.__oblique_shock.get_Mach_angle(M[i])
i0 = 0
for k in range(n):
i0 += (n - k)
Kp[i0] = Kp[i0 - 1]
Km[i0] = Km[i0 - 1]
t[i0] = t[i0 - 1]
v[i0] = v[i0 - 1]
M[i0] = M[i0 - 1]
mu[i0] = mu[i0 - 1]
i0 += 1
# calculate node coordinates
Rt = 1.0
x[0] = -Rt / np.tan(t[0] - mu[0])
y[0] = 0.
for i in range(1, n):
dxdy_av = 0.5 * (mu[i - 1] + t[i - 1] + mu[i] + t[i])
dxdy = t[i] - mu[i]
x[i] = (Rt - y[i - 1] + x[i - 1] * np.tan(dxdy_av)) / (np.tan(dxdy_av) - np.tan(dxdy))
y[i] = np.tan(dxdy) * x[i] + Rt
dxdy_av = 0.5 * (t[n - 1] + t[n] + mu[n - 1] + mu[n])
t_av = 0.5 * (t_max + t[n])
x[n] = (y[n - 1] - Rt - x[n - 1] * np.tan(dxdy_av)) / (np.tan(t_av) - np.tan(dxdy_av))
y[n] = Rt + x[n] * np.tan(t_av)
i0 = n + 1
for k in range(1, n):
j = i0 - (n - k + 1)
dxdy_av = 0.5 * (t[j] - mu[j] + t[i0] - mu[i0])
x[i0] = x[j] - y[j] / np.tan(dxdy_av)
y[i0] = 0.
for kp in range(1, n - k):
i = i0 + kp
j = i - (n - k + 1)
s1 = np.tan(0.5 * (t[i] + t[i - 1] + mu[i] + mu[i - 1]))
s2 = np.tan(0.5 * (t[i] + t[j] - mu[i] - mu[j]))
x[i] = (y[j] - y[i - 1] + s1 * x[i - 1] - s2 * x[j]) / (s1 - s2)
y[i] = y[i - 1] + (x[i] - x[i - 1]) * s1
i = i0 + n - k
j = i0 - 1
dxdy_av = 0.5 * (t[i - 1] + t[i] + mu[i - 1] + mu[i])
t_av = 0.5 * (t[j] + t[i])
x[i] = ((x[j] * np.tan(t_av) - y[j] + y[i - 1] - x[i - 1] * np.tan(dxdy_av)) /
(np.tan(t_av) - np.tan(dxdy_av)))
y[i] = y[j] + (x[i] - x[j]) * np.tan(t_av)
i0 += (n - k + 1)
print("Node K_ K+ theta v M mu x y")
for i in range(NN):
print(('{0:4d} '
'{1:7.3f} '
'{2:7.3f} '
'{3:7.3f} '
'{4:7.3f} '
'{5:7.3f} '
'{6:7.3f} '
'{7:7.3f} '
'{8:7.3f}').format(i + 1,
Km[i] * 180 / np.pi,
Kp[i] * 180 / np.pi,
t[i] * 180 / np.pi,
v[i] * 180 / np.pi,
M[i],
mu[i] * 180 / np.pi,
x[i],
y[i]))
return Km, Kp, t, v, M, mu, x, y