"""ODE-based q1d calculator : area changes, friction and heat addition"""
import os
from typing import Tuple, Any
from collections.abc import Callable
if os.environ.get('__USE_CUDA', '0') == '1':
import cupy as np
else:
import numpy as np
from .calculator import G3DCalculator
[docs]
class Q1DODESolver(G3DCalculator):
"""
quasi-one-dimensional ODE solver handles cross-section area changes,
heat addition, and friction
"""
def __init__(self, *args: Any, **kwargs: Any) -> None:
"""Initialize calculator parameters"""
# default
self.__default_r_fun = None
self.__default_q_fun = None
self.__default_f_fun = None
self.__default_r = 1.
self.__default_x1 = 0.
self.__default_rho1 = 1.
self.__default_u1 = 0.
self.__default_p1 = 1.e5
self.__default_dksi = 0.
self.__default_L = 1.e20
self.__default_h = 1.e-4
# initialize Q1D-solver parameters
self._r_fun = self.__default_r_fun
self._q_fun = self.__default_q_fun
self._f_fun = self.__default_f_fun
self._r = self.__default_r
self._x1 = self.__default_x1
self._rho1 = self.__default_rho1
self._u1 = self.__default_u1
self._p1 = self.__default_p1
self._dksi = self.__default_dksi
self._L = self.__default_L
self._h = self.__default_h
# run parent class constructor
super().__init__(*args, **kwargs)
def _set_default(self) -> None:
super()._set_default()
self.__nullify()
def __nullify(self) -> None:
self._r_fun = self.__default_r_fun
self._q_fun = self.__default_q_fun
self._f_fun = self.__default_f_fun
self._r = self.__default_r
self._x1 = self.__default_x1
self._rho1 = self.__default_rho1
self._u1 = self.__default_u1
self._p1 = self.__default_p1
self._dksi = self.__default_dksi
self._L = self.__default_L
self._h = self.__default_h
[docs]
def reset(self, *args: Any, **kwargs: Any) -> None:
"""Resets solver parameters to their default values"""
self.__nullify()
self.set(*args, **kwargs)
[docs]
def set_from_key_value(self, key: str, value: Any) -> bool:
"""Set parameters from provided key => value input"""
if super().set_from_key_value(key, value):
return True
return self.__set_from_key_value(key, value)
def __set_from_key_value(self, key: str, value: Any) -> bool:
if key == 'r_fun':
self._r_fun = value
return True
elif key == 'q_fun':
self._q_fun = value
return True
elif key == 'f_fun':
self._f_fun = value
return True
elif key == 'r':
self._r = value
return True
elif key == 'x1':
self._x1 = value
return True
elif key == 'u1':
self._u1 = value
return True
elif key == 'rho1':
self._rho1 = value
return True
elif key == 'p1':
self._p1 = value
return True
elif key == 'dksi':
self._dksi = value
return True
elif key == 'L':
self._L = value
return True
elif key == 'h':
self._h = value
return True
return False
@property
def r_fun(self) -> Callable[[float], float]:
"""Cross-section radius as a function of axial coordinate r(x) [m]"""
return self._r_fun
@r_fun.setter
def r_fun(self, f: Callable[[float], float]) -> None:
self._r_fun = f
@property
def q_fun(self) -> Callable[[float], float]:
"""Added heat as a function of axial coordinate q(x) [J/kg]"""
return self._q_fun
@q_fun.setter
def q_fun(self, f: Callable[[float], float]) -> None:
self._q_fun = f
@property
def f_fun(self) -> Callable[[float], float]:
"""Friction factor as a function of axial coordinate f(x) []"""
return self._f_fun
@f_fun.setter
def f_fun(self, f: Callable[[float], float]) -> None:
self._f_fun = f
@property
def r(self) -> float:
"""Cross-section radius (r [m], used unless function r_func is specified)"""
return self._r
@r.setter
def r(self, value: float) -> None:
self._r = value
@property
def x1(self) -> float:
"""Start coordinate (x1 [m], default value: 0.)"""
return self._x1
@x1.setter
def x1(self, value: float) -> None:
self._x1 = value
@property
def u1(self) -> float:
"""Initial velocity (u1 [m/s])"""
return self._u1
@u1.setter
def u1(self, value: float) -> None:
self._u1 = value
@property
def p1(self) -> float:
"""Initial pressure (p1 [Pa])"""
return self._p1
@p1.setter
def p1(self, value: float) -> None:
self._p1 = value
@property
def rho1(self) -> float:
"""Initial density (rho1 [kg/m^3])"""
return self._rho1
@rho1.setter
def rho1(self, value: float) -> None:
self._rho1 = value
@property
def dksi(self) -> float:
"""Axial spacing between function evaluations"""
return self._dksi
@dksi.setter
def dksi(self, value: float) -> None:
self._dksi = value
@property
def L(self) -> float:
"""Pipe length (L [m])"""
return self._L
@L.setter
def L(self, value: float) -> None:
self._L = value
@property
def h(self) -> float:
"""Axial spacing used for numerical differentiation (h [m], default value: 1.e-4)"""
return self._h
@h.setter
def h(self, value: float) -> None:
self._h = value
def __get_derivatives(self, x: float, rho: float, u: float, p: float) -> Tuple[float, float, float, float]:
x2 = self._x1 + self._L
r = self._r
if self._r_fun is not None:
r = self._r_fun(x)
A = np.pi * r ** 2
Lp = 2. * np.pi * r
qw = 0.
if self._q_fun is not None:
qw = self._q_fun(x)
tauw = 0.
if self._f_fun is not None:
tauw = self._f_fun(x) * rho * u ** 2 / 8.
dAdx = 0.
if self._r_fun is not None:
if (x - self._h) < self._x1:
dAdx = (np.pi * self._r_fun(x + self._h) ** 2 - A) / self._h
elif (x + self._h) > x2:
dAdx = (A - np.pi * self._r_fun(x - self._h) ** 2) / self._h
else:
dAdx = np.pi * (self._r_fun(x + self._h) ** 2 - self._r_fun(x - self._h) ** 2) / (2. * self._h)
drho_dksi = (1. / A) * (-rho * u ** 2 * dAdx + tauw * Lp + self._gas.gm1 * (qw + tauw * u) * Lp / u)
du_dksi = (1. / (rho * A)) * (self._gas.gamma * p * u * dAdx - u * tauw * Lp - self._gas.gm1 *
(qw + tauw * u) * Lp)
dP_dksi = (1. / A) * (-self._gas.gamma * p * u ** 2 * dAdx + self._gas.gamma * p /
rho * tauw * Lp + self._gas.gm1 * (qw + tauw * u) * Lp * u)
dx_dksi = u ** 2 - self._gas.gamma * p / rho
return drho_dksi, du_dksi, dP_dksi, dx_dksi
[docs]
def solve(self, *args: Any, **kwargs: Any) -> Tuple[np.ndarray, ...]:
"""
Q1D ODE solver for investigations of combined affects of area changes, heat addition, and friction
input: functions describing area distribution (r_fun - r(x) [m]), heat addition (q(x) [J/kg]), and
friction factor (f_fun - f(x) [], default value: None),
initial conditions: axial position (x1 [m], default value: 0.), density (rho1 [kg/m^3], default value: 1.),
velocity (u1 [m/s], default value: 0.), pressure (p1 [Pa], default value: 1.e5),
axial increment per step (dksi [m], default value: 0.), pipe length (L [m], default value: 1.e20),
delta (h [m], default value: 1.e-4)
example: x, rho, u, p, T, M, A = solve(r_fun=r_fun, q_fun=q_fun, f_fun=f_fun, x1=x1, rho1=rho1,
u1=u1, p1=p1, dksi=dksi, L=L, h=h)
"""
self.set(*args, **kwargs)
T1 = self._gas.get_temperature(self._p1, self._rho1)
M1 = self._u1 / self._gas.get_speed_of_sound(T1)
dksi = abs(self._dksi)
if M1 < 1.:
dksi *= -1.
x = self._x1
r = self._r
if self._r_fun is not None:
r = self._r_fun(x)
rho = self._rho1
u = self._u1
p = self._p1
N = 1
while True:
last_entry = False
drho_dksi, du_dksi, dp_dksi, dx_dksi = self.__get_derivatives(x, rho, u, p)
if dx_dksi >= 0.:
break
dk = dksi
if (x + dx_dksi * dk - self._x1) > self._L:
dk = (self._L - self._x1 - x) / dx_dksi
last_entry = True
rho += drho_dksi * dk
u += du_dksi * dk
p += dp_dksi * dk
x += dx_dksi * dk
N += 1
if last_entry:
break
x = np.zeros(N)
rho = np.zeros(N)
u = np.zeros(N)
p = np.zeros(N)
T = np.zeros(N)
M = np.zeros(N)
A = np.zeros(N)
x[0] = self._x1
rho[0] = self._rho1
u[0] = self._u1
p[0] = self._p1
T[0] = T1
M[0] = M1
A[0] = np.pi * r ** 2
for i in range(N - 1):
last_entry = False
drho_dksi, du_dksi, dp_dksi, dx_dksi = self.__get_derivatives(x[i], rho[i], u[i], p[i])
dk = dksi
if (x[i] + dx_dksi * dk - self._x1) > self._L:
dk = (self._L - self._x1 - x[i]) / dx_dksi
last_entry = True
rho[i + 1] = rho[i] + drho_dksi * dk
u[i + 1] = u[i] + du_dksi * dk
p[i + 1] = p[i] + dp_dksi * dk
x[i + 1] = x[i] + dx_dksi * dk
if self._r_fun is not None:
r = self._r_fun(x[i + 1])
A[i + 1] = np.pi * r ** 2
T[i + 1] = self._gas.get_temperature(p[i + 1], rho[i + 1])
M[i + 1] = u[i + 1] / self._gas.get_speed_of_sound(T[i + 1])
if last_entry:
break
return x, rho, u, p, T, M, A