"""Shock-tube relations calculator"""
import os
from typing import Union, Any
if os.environ.get('__USE_CUDA', '0') == '1':
import cupy as np
else:
import numpy as np
from g3d.gas import Gas
from g3d.defs import CFlowVariableType
from .calculator import G3DCalculator
[docs]
class G3DShockTubeCalculator(G3DCalculator):
"""Class collecting functions related to shock-tube problems"""
def __init__(self, *args: Any, **kwargs: Any) -> None:
"""Initialize calculator parameters"""
super().__init__(*args, **kwargs)
self._gas1 = Gas()
self._gas2 = Gas()
def _set_default(self) -> None:
super()._set_default()
self.__PR21 = self._vars.add(name="PR21", long_name="Pressure ratio P2/P1")
self.__PR41 = self._vars.add(name="PR41", long_name="Pressure ratio P4/P1")
self.__T1 = self._vars.add(name="T1", long_name="Temperature (region 1)", type=CFlowVariableType.TEMPERATURE)
self.__T4 = self._vars.add(name="T4", long_name="Temperature (region 4)", type=CFlowVariableType.TEMPERATURE)
@property
def gas1(self) -> Gas:
"""Driver gas (high-pressure chamber)"""
return self._gas1
@gas1.setter
def gas1(self, g: Gas) -> None:
self._gas1 = g
@property
def gas2(self) -> Gas:
"""Driven gas (low-pressure chamber)"""
return self._gas2
@gas2.setter
def gas2(self, g: Gas) -> None:
self._gas2 = g
[docs]
def get_PR41(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
calculates driver/driven gas pressure ratio (PR41 [], shock tube initial condition) from
incident shock pressure ratio (PR21 []), driven section temperature (T1 [k]), and
driver section temperature (T4 [K])
input: incident shock wave pressure ratio (PR21 [K]), driven section temperature (T1 [K]),
driver section temperature (T4 [K])
example: PR41 = get_PR41(PR21=PR21, T1=T1, T4=T4)
"""
return self._eval(["PR21", "T1", "T4"], self.__get_PR41, *args, **kwargs)
def __get_PR41(self, PR21: float = None, T1: float = None, T4: float = None) -> float:
if PR21 is None:
PR21 = self.__PR21.value
if T1 is None:
T1 = self.__T1.value
if T4 is None:
T4 = self.__T4.value
C1 = self._gas1.get_speed_of_sound(T1)
C4 = self._gas2.get_speed_of_sound(T4)
g1 = self._gas1.gamma
g4 = self._gas2.gamma
return PR21 * np.power((1. - (((g4 - 1.) * (C1 / C4) * (PR21 - 1.))
/ np.sqrt(2. * g1 * (2. * g1 + (g1 + 1.) * (PR21 - 1.))))),
(-2. * g4 / (g4 - 1.)))
[docs]
def find_PR21_from_PR41(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
.. deprecated:: 2.6
redirects to :func:`get_PR21_from_PR41`
"""
return self.get_PR21_from_PR41(*args, **kwargs)
[docs]
def get_PR21_from_PR41(self, *args: Any, **kwargs: Any) -> Union[Union[tuple[Any, ...], list[Any]], Any]:
"""
finds incident shock pressure ratio (PR21 []) matching specified shock-tube pressure ratio (PR41 [],
driver section pressure over driven section pressure), driver section temperature (T4 [K]), and driven
section temperature (T1 [K]).
the pressure ration is calculated using the built-in Newton-Raphson solver.
input: driver pressure of driven pressure (PR41 []), driver section temperature (T4 [K]),
driven section temperature (T1 [K])
example: PR21 = get_PR21_from_PR41(PR41=PR41, T4=T4, T1=T1)
"""
self._newton.reset()
return self._eval(["PR41", "T4", "T1"], self.__get_PR21_from_PR41, *args, **kwargs)
def __get_PR21_from_PR41(self, PR41: float = None, T4: float = None, T1: float = None) -> float:
if PR41 is None:
PR41 = self.__PR41.value
if T1 is not None:
self.__T1.value = T1
if T4 is not None:
self.__T4.value = T4
newton = self._get_newton_kernel()
newton.target_value = PR41
newton.start_value = 10.
# run Newton solver to find deflection angle
return newton.solve(self.__get_PR41)