"""Oblique-shock calculator"""
import os
from typing import Tuple, Union, Any
if os.environ.get('__USE_CUDA', '0') == '1':
import cupy as np
from cupyx.scipy.interpolate import CubicSpline
else:
import numpy as np
from scipy.interpolate import CubicSpline
from .calculator import G3DCalculator
from .normal_shock import G3DNormalShockCalculator
from g3d.defs import CFlowRegime
from g3d.defs import CFlowVariableType
from g3d.defs import AngleType
[docs]
class G3DObliqueShockCalculator(G3DCalculator):
"""Class collecting functions related to oblique shocks"""
def __init__(self, *args: Any, **kwargs: Any) -> None:
"""Initialize calculator parameters"""
super().__init__(*args, **kwargs)
self.__normal_shock = G3DNormalShockCalculator(update_code_type=False)
def _set_default(self) -> None:
super()._set_default()
self.__M = self._vars.add(name="M", long_name="Mach number")
self.__M1 = self._vars.add(name="M1", long_name="Upstream Mach number")
self.__M2 = self._vars.add(name="M2", long_name="Downstream Mach number")
self.__Mn1 = self._vars.add(name="Mn1", long_name="Upstream shock-normal Mach number")
self.__beta = self._vars.add(name="beta", long_name="Shock angle", type=CFlowVariableType.RADIANS)
self.__theta = self._vars.add(name="theta", long_name="Flow deflection angle", type=CFlowVariableType.RADIANS)
self.__PR = self._vars.add(name="PR", long_name="Pressure ratio")
def _associate(self) -> None:
super()._associate()
self.__normal_shock.gas = self._gas
[docs]
def get_Mach_angle(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
calculates Mach angle (mu [rad*]) corresponding to specified flow Mach number (M [])
* the mach angle will be given in radians unless angle_type=AngleType.DEGREES is specified
input: Mach number (M [])
example: mu = get_Mach_angle(M=M)
"""
return self._eval(["M"], self.__get_Mach_angle, *args, **kwargs)
def __get_Mach_angle(self, M: float = None, angle_type: Any = AngleType.UNDEFINED) -> float:
if M is None:
M = self.__M.value
mu = np.arcsin(1. / M)
return self._angle_export(mu, angle_type=angle_type)
[docs]
def get_theta(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
calculates flow deflection angle (theta [rad*]) from flow Mach number (M []) and shock angle (beta [rad*])
* shock angle should be given in radians unless angle_type=AngleType.DEGREES is specified
(please note that angle_type=AngleType.DEGREES will result in a calculated deflection angle in
degrees as well)
input: flow Mach number (M []), shock angle (beta [rad*])
example: theta = get_theta(M=M, beta=beta)
"""
return self._eval(["beta", "M"], self.__get_theta, *args, **kwargs)
def __get_theta(self, beta: float = None, M: float = None,
angle_type: Any = AngleType.UNDEFINED) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
if M is None:
M = self.__M.value
if beta is None:
beta = self.__beta.value
beta = self._angle_import(beta, angle_type=angle_type)
theta = np.arctan(2. * np.tan((np.pi / 2) - beta) * (np.power(M * np.sin(beta), 2) - 1.)
/ (M * M * (self._gas.gamma + np.cos(2. * beta)) + 2.))
return self._angle_export(theta, angle_type=angle_type)
[docs]
def get_normal_Mach_1(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
calculates the upstream shock-normal Mach number (the flow Mach number calculated using the velocity
component normal to the shock) (Mn1 []) using the upstream flow Mach number (M1 []) and the
shock angle (beta [rad*])
* shock angle should be given in radians unless angle_type=AngleType.DEGREES is specified
input: upstream Mach number (M1 []), shock angle (beta [rad*])
example: Mn1 = get_normal_Mach_1(M1=M1, beta=beta)
"""
return self._eval(["M1", "beta"], self.__get_normal_Mach_1, *args, **kwargs)
def __get_normal_Mach_1(self, M1: float = None, beta: float = None, angle_type: Any = AngleType.UNDEFINED) -> float:
if M1 is None:
M1 = self.__M1.value
if beta is None:
beta = self.__beta.value
beta = self._angle_import(beta, angle_type=angle_type)
return M1 * np.sin(beta)
[docs]
def get_normal_Mach_2(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
calculates the downstream shock-normal Mach number (the flow Mach number calculated using the velocity
component normal to the shock) (Mn2 []) using the upstream shock-normal Mach number (Mn1 [])
input: upstream shock-normal Mach number (Mn1 [])
example: Mn2 = get_normal_Mach_2(Mn1=Mn1)
"""
return self._eval(["Mn1"], self.__get_normal_Mach_2, *args, **kwargs)
def __get_normal_Mach_2(self, Mn1: float = None) -> float:
if Mn1 is None:
Mn1 = self.__Mn1.value
return self.__normal_shock.get_downstream_Mach(Mn1)
[docs]
def get_density_ratio(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
calculates density ratio (ROR []) over an oblique shock with the upstream flow Mach number (M1 [])
and the shock angle (beta [rad*])
* shock angle should be given in radians unless angle_type=AngleType.DEGREES is specified
input: upstream Mach number (M1 []), shock angle (beta [rad*])
example: ROR = get_density_ratio(M1=M1, beta=beta)
"""
return self._eval(["M1", "beta"], self.__get_density_ratio, *args, **kwargs)
def __get_density_ratio(self, M1: float = None, beta: float = None, angle_type: Any = AngleType.UNDEFINED) -> float:
if M1 is None:
M1 = self.__M1.value
if beta is None:
beta = self.__beta.value
beta = self._angle_import(beta, angle_type=angle_type)
Mn1 = self.__get_normal_Mach_1(M1, beta, AngleType.RADIANS)
return (self._gas.gp1 * Mn1 * Mn1) / (self._gas.gm1 * Mn1 * Mn1 + 2.0)
[docs]
def get_pressure_ratio(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
calculates pressure ratio (PR []) over an oblique shock with the upstream flow Mach number (M1 [])
and the shock angle (beta [rad*])
* shock angle should be given in radians unless angle_type=AngleType.DEGREES is specified
input: upstream Mach number (M1 []), shock angle (beta [rad*])
example: PR = get_pressure_ratio(M1=M1, beta=beta)
"""
return self._eval(["M1", "beta"], self.__get_pressure_ratio, *args, **kwargs)
def __get_pressure_ratio(self, M1: float = None, beta: float = None,
angle_type: Any = AngleType.UNDEFINED) -> float:
if M1 is None:
M1 = self.__M1.value
if beta is None:
beta = self.__beta.value
beta = self._angle_import(beta, angle_type=angle_type)
Mn1 = self.__get_normal_Mach_1(M1, beta, AngleType.RADIANS)
return 1. + ((2. * self._gas.gamma) / self._gas.gp1) * (Mn1 * Mn1 - 1.)
[docs]
def get_temperature_ratio(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
calculates temperature ratio (TR []) over an oblique shock with the upstream flow Mach number (M1 [])
and the shock angle (beta [rad*])
* shock angle should be given in radians unless angle_type=AngleType.DEGREES is specified
input: upstream Mach number (M1 []), shock angle (beta [rad*])
example: TR = get_pressure_ratio(M1=M1, beta=beta)
"""
return self._eval(["M1", "beta"], self.__get_temperature_ratio, *args, **kwargs)
def __get_temperature_ratio(self, M1: float = None, beta: float = None,
angle_type: Any = AngleType.UNDEFINED) -> float:
if M1 is None:
M1 = self.__M1.value
if beta is None:
beta = self.__beta.value
beta = self._angle_import(beta, angle_type=angle_type)
return (self.__get_pressure_ratio(M1, beta, AngleType.RADIANS) /
self.__get_density_ratio(M1, beta, AngleType.RADIANS))
[docs]
def get_downstream_Mach(self, *args: Any, **kwargs: Any) -> Union[np.float64, list[np.float64], float]:
"""
calculates downstream flow Mach number (M2 []) from upstream flow Mach number (M1 []),
shock angle (beta [rad*]) and flow deflection angle (theta [rad*])
* shock angle and flow deflection angle should be given in radians unless angle_type=AngleType.DEGREES is
specified
input: upstream flow Mach number (M1 []), shock angle (beta [rad*]), flow deflection angle (theta [rad*])
example: M2 = get_downstream_Mach(M1=M1, beta=beta, theta=theta)
"""
return self._eval(["M1", "beta", "theta"], self.__get_downstream_Mach, *args, **kwargs)
def __get_downstream_Mach(self, M1: float = None, beta: float = None, theta: float = None,
angle_type: Any = AngleType.UNDEFINED) -> float:
if M1 is None:
M1 = self.__M1.value
if beta is None:
beta = self.__beta.value
if theta is None:
theta = self.__theta.value
beta, theta = self._angle_import(beta, theta, angle_type=angle_type)
Mn1 = self.__get_normal_Mach_1(M1, beta, AngleType.RADIANS)
Mn2 = self.__get_normal_Mach_2(Mn1)
return Mn2 / np.sin(beta - theta)
def __find_theta_max(self, M: float, b0: float, db: float) -> float:
b = b0
t_old = -10
t_max = -10
while b < (0.5 * np.pi - db):
b = b + db
t = self.__get_theta(b, M, AngleType.RADIANS)
if t < t_old:
return t_max
if t > t_max:
t_max = t
t_old = t
return -1.
[docs]
def get_theta_max(self, *args: Any, **kwargs: Any) -> Union[np.float64, list[np.float64], float]:
"""
calculates maximum flow deflection angle (theta_max [rad*]) for a specified flow Mach number (M [])
* the calculated flow deflection angle will be given in radians unless angle_type=AngleType.DEGREES is specified
input: flow Mach number (M [])
example: theta_max = get_theta_max(M=M)
"""
self._reset_eval(mode=CFlowRegime.WEAK)
self._newton.reset()
return self._eval(["M"], self.__get_theta_max, *args, **kwargs)
def __get_theta_max(self, M: float = None, angle_type: Any = AngleType.UNDEFINED) -> float:
if M is None:
M = self.__M.value
b0 = self.__get_Mach_angle(M, AngleType.RADIANS)
db = (0.5 * np.pi - b0) / 10.
t_old = -10
t_max = 0.
while np.abs(t_max - t_old) > 1.e-6 * t_max:
t_old = t_max
t_max = self.__find_theta_max(M, b0, db)
b0 = self.__get_beta(t_max, M, AngleType.RADIANS)
db = db / 10.
return self._angle_export(t_max, angle_type=angle_type)
[docs]
def find_beta(self, *args, **kwargs):
"""
.. deprecated:: 2.6
redirects to :func:`get_beta`
"""
return self.get_beta(*args, **kwargs)
[docs]
def get_beta(self, *args: Any, **kwargs: Any) -> Union[np.float64, list[np.float64], float]:
"""
calculates shock angle (beta [rad*]) for specified flow Mach number (M []) and flow deflection angle
(theta [rad*])
* flow deflection angle should be given in radians unless angle_type=AngleType.DEGREES is specified
(please note that angle_type=AngleType.DEGREES will result in a calculated shock angle in degrees as well)
the oblique shock will be assumed to fall within the weak shock solution unless mode=CFlowRegime.STRONG is
specified
the shock angle is obtained using the built-in Newton-Raphson solver
input: flow Mach number (M []), flow deflection angle (theta [rad*])
example: beta = get_beta(M=M, theta=theta)
"""
self._reset_eval(mode=CFlowRegime.WEAK)
self._newton.reset()
return self._eval(["theta", "M"], self.__get_beta, *args, **kwargs)
def __get_beta(self, theta: float = None, M: float = None, angle_type: Any = AngleType.UNDEFINED) -> float:
if M is None:
M = self.__M.value
if theta is None:
theta = self.__theta.value
self.backup()
self.__M.value = M
theta = self._angle_import(theta, angle_type=angle_type)
newton = self._get_newton_kernel()
newton.target_value = theta
newton.start_value = 0.
# newton.verbose = True
if self._mode == CFlowRegime.WEAK:
newton.start_value = self.__get_Mach_angle(M, AngleType.RADIANS)
else:
newton.start_value = np.pi / 2.
# run Newton solver to find deflection angle
tmp = self._angle_type
self._angle_type = AngleType.RADIANS
beta = newton.solve(self.__get_theta)
self._angle_type = tmp
self.restore()
return self._angle_export(beta, angle_type=angle_type)
[docs]
def get_Mach(self, *args: Any, **kwargs: Any) -> Union[np.float64, list[np.float64], float]:
"""
calculate upstream Mach number (M []) corresponding to specified shock angle (beta [rad*]) and
flow deflection angle (theta [rad*])
* shock angle and flow deflection angle should be given in radians unless angle_type=AngleType.DEGREES is
specified
input: flow deflection angle (theta [rad*]), shock angle (beta [rad*])
example: M = get_Mach(theta=theta, beta=beta)
"""
return self._eval(["theta", "beta"], self.__get_Mach, *args, **kwargs)
def __get_Mach(self, theta: float = None, beta: float = None, angle_type: Any = AngleType.UNDEFINED) -> float:
if theta is None:
theta = self.__theta.value
if beta is None:
beta = self.__beta.value
beta, theta = self._angle_import(beta, theta, angle_type=angle_type)
A = np.tan(theta) / (2.0 * np.tan((np.pi / 2.0) - beta))
B = self._gas.gamma + np.cos(2.0 * beta)
C = np.power(np.sin(beta), 2)
return np.sqrt((2.0 * A + 1) / (C - A * B))
[docs]
def get_angles_from_PR_and_M1(self, *args: Any, **kwargs: Any) -> Union[np.float64, list[np.float64], float]:
"""
calculates shock angle (beta [rad*]) and flow deflection angle (theta [rad*]) corresponding to specified
pressure ratio (PR []) and upstream flow Mach number (M1 [])
* shock angle and flow deflection angle will be given in radians unless angle_type=AngleType.DEGREES is
specified
input: pressure ratio (PR []), upstream flow Mach number (M1 [])
example: beta, theta = get_angles_from_PR_and_M1(PR=PR, M1=M1)
"""
self._reset_eval(outdim=2)
return self._eval(["PR", "M1"], self.__get_angles_from_PR_and_M1, *args, **kwargs)
def __get_angles_from_PR_and_M1(self, PR: float = None, M1: float = None,
angle_type: Any = AngleType.UNDEFINED) -> float:
if PR is None:
PR = self.__PR.value
if M1 is None:
M1 = self.__M1.value
Mn1 = self.__normal_shock.get_upstream_Mach_from_pressure_ratio(PR)
beta = np.arcsin(Mn1 / M1)
theta = self.__get_theta(beta, M1, AngleType.RADIANS)
return self._angle_export(beta, theta, angle_type=angle_type)
[docs]
def find_Mach_min(self, *args: Any, **kwargs: Any) -> Union[np.float64, list[np.float64], float]:
"""
.. deprecated:: 2.6
redirects to :func:`get_Mach_min`
"""
return self.get_Mach_min(*args, **kwargs)
[docs]
def get_Mach_min(self, *args: Any, **kwargs: Any) -> Union[np.float64, list[np.float64], float]:
"""
calculates the minimum flow Mach number (M []) for a specified flow deflection angle (theta [rad*])
* flow deflection angle should be given in radians unless angle_type=AngleType.DEGREES is specified
input: flow deflection angle (theta [rad*])
example: M = get_Mach_min(theta=theta)
"""
return self._eval(["theta"], self.__get_Mach_min, *args, **kwargs)
def __get_Mach_min(self, theta: float = None, angle_type: Any = AngleType.UNDEFINED) -> float:
if theta is None:
theta = self.__theta.value
theta = self._angle_import(theta, angle_type=angle_type)
M = np.linspace(1, 20, 100)
T = np.zeros(len(M))
for i in range(1, len(M)):
T[i] = self.__get_theta_max(M[i], AngleType.RADIANS)
cs = CubicSpline(T, M)
if theta <= max(T):
return np.float64(np.round(cs(theta) + 0.5 * 10 ** (-self._precision), decimals=self._precision))
return -1.
[docs]
def analyze_oblique_shock(self, *args: Any, **kwargs: Any) -> Union[np.float64, list[np.float64], float]:
"""
calculates shock angle (beta [rad*]), downstream flow Mach number (M2 []), shock pressure ratio (PR21 []),
and shock temperature ratio (TR21 []) corresponding to the specified upstream flow Mach number (M1 []) and
flow deflection angle (theta [rad*])
* shock angle should be given in radians unless angle_type=AngleType.DEGREES is specified
(please note that angle_type=AngleType.DEGREES will result in a calculated flow deflection angle in
degrees as well)
the oblique shock will be assumed to fall within the weak shock solution unless mode=CFlowRegime.STRONG is
specified
the built-in Newton-Raphson solver is used for some of the calculations
input: upstream flow Mach number (M1 []), flow deflection angle (theta [rad*])
example: beta, M2, PR21, TR21 = analyze_oblique_shock(M1=M1, theta=theta)
"""
self._reset_eval(mode=CFlowRegime.WEAK, outdim=4)
self._newton.reset()
return self._eval(["M1", "theta"], self.__analyze_oblique_shock, *args, **kwargs)
def __analyze_oblique_shock(self, M1: float = None, theta: float = None,
angle_type: Any = AngleType.UNDEFINED) -> Tuple[float, float, float, float]:
if M1 is None:
M1 = self.__M1.value
if theta is None:
theta = self.__theta.value
self.backup()
self.__M.value = M1
theta = self._angle_import(theta, angle_type=angle_type)
beta = self.__get_beta(theta, M1, AngleType.RADIANS)
M2 = self.__get_downstream_Mach(M1, beta, theta, AngleType.RADIANS)
PR21 = self.__get_pressure_ratio(M1, beta, AngleType.RADIANS)
TR21 = self.__get_temperature_ratio(M1, beta, AngleType.RADIANS)
beta = self._angle_export(beta, angle_type=angle_type)
self.restore()
return beta, M2, PR21, TR21
[docs]
def get_downstream_Mach_from_theta_and_M1(self, *args: Any, **kwargs: Any) \
-> Union[np.float64, list[np.float64], float]:
"""
calculates flow Mach number downstream of oblique shock (M2 []) from specified upstream flow Mach number (M1 [])
and flow deflection angle (theta [rad*])
* flow deflection angle should be given in radians unless angle_type=AngleType.DEGREES is specified
the oblique shock will be assumed to fall within the weak shock solution unless mode=CFlowRegime.STRONG is
specified
the built-in Newton-Raphson solver is used for some of the calculations
input: upstream flow Mach number (M1 []), flow deflection angle (theta [rad*])
example: M2 = get_downstream_Mach_from_theta_and_M1(M1=M1, theta=theta)
"""
self._reset_eval(mode=CFlowRegime.WEAK)
self._newton.reset()
return self._eval(["M1", "theta"], self.__get_downstream_Mach_from_theta_and_M1, *args, **kwargs)
def __get_downstream_Mach_from_theta_and_M1(self, M1: float = None, theta: float = None,
angle_type: Any = AngleType.UNDEFINED) -> float:
if M1 is None:
M1 = self.__M1.value
if theta is None:
theta = self.__theta.value
self.backup()
self.__M.value = M1
theta = self._angle_import(theta, angle_type=angle_type)
beta = self.__get_beta(theta, M1, AngleType.RADIANS)
self.restore()
return self.__get_downstream_Mach(M1, beta, theta, AngleType.RADIANS)
[docs]
def find_M1_from_theta_and_M2(self, *args: Any, **kwargs: Any) -> Union[np.float64, list[np.float64], float]:
"""
.. deprecated:: 2.6
redirects to :func:`get_M1_from_theta_and_M2`
"""
return self.get_M1_from_theta_and_M2(*args, **kwargs)
[docs]
def get_M1_from_theta_and_M2(self, *args: Any, **kwargs: Any) -> Union[np.float64, list[np.float64], float]:
"""
calculates flow Mach number upstream of oblique shock (M1 []) from specified downstream flow Mach number (M2 [])
and flow deflection angle (theta [rad*])
* flow deflection angle should be given in radians unless angle_type=AngleType.DEGREES is specified
the oblique shock will be assumed to fall within the weak shock solution unless mode=CFlowRegime.STRONG is
specified
the built-in Newton-Raphson solver is used for some of the calculations
input: downstream flow Mach number (M2 []), flow deflection angle (theta [rad*])
example: M2 = get_M1_from_theta_and_M2(M2=M2, theta=theta)
"""
self._reset_eval(mode=CFlowRegime.WEAK, precision=2)
self._newton.reset()
return self._eval(["theta", "M2"], self.__get_M1_from_theta_and_M2, *args, **kwargs)
def __get_M1_from_theta_and_M2(self, theta: float = None, M2: float = None,
angle_type: Any = AngleType.UNDEFINED) -> float:
if M2 is None:
M2 = self.__M2.value
if theta is None:
theta = self.__theta.value
theta = self._angle_import(theta, angle_type=angle_type)
newton = self._get_newton_kernel()
newton.start_value = self.__get_Mach_min(theta, AngleType.RADIANS) * 1.01
newton.target_value = M2
return newton.solve(self.__get_downstream_Mach_from_theta_and_M1)
# ARRAY EXTRACTION NOT CONVERTED
# ================================================================
# ARRAY FUNCTIONS
# ================================================================
[docs]
def get_theta_beta_Mach_array(self, M1: float, N: int = 100) -> Tuple[np.ndarray, ...]:
"""
extracts shock angle vector and flow deflection angle vector for specified upstream flow Mach number (M1 [])
input: upstream flow Mach number (M1 []), vector length (N, default value: 100)
example: theta, beta = get_theta_beta_Mach_array(M1, N)
"""
beta_array = np.zeros(N)
theta_array = np.zeros(N)
b0 = self.__get_Mach_angle(M1)
db = (0.5 * np.pi - b0) / N
beta_array[0] = b0
theta_array[0] = 0.
for i in np.arange(1, N - 1):
beta_array[i] = beta_array[i - 1] + db
theta_array[i] = self.__get_theta(beta_array[i], M1)
beta_array[N - 1] = 0.5 * np.pi
theta_array[N - 1] = 0.
return theta_array, beta_array