OpenEngine/ode_functions.py at main · WiktorBoltrukiewicz/OpenEngine · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
"""
ode_functions.py — Prawe strony równań różniczkowych (ODE) silnika rakietowego.

Zawiera dwie wersje układu ODE:
  1. myNozzleODE      — uproszczona (tylko geometria dyszy)
  2. myNozzleODEfull2  — pełna (z tarciem i wymianą ciepła)

Zmienne stanu Y = [N, P, T]:
  N — kwadrat liczby Macha (M^2)
  P — ciśnienie statyczne [Pa]
  T — temperatura statyczna [K]

Odpowiada funkcjom myNozzleODE, myNozzleODEfull2 i eventN1 z MATLAB-a.
"""

import numpy as np


def my_nozzle_ode(x, Y, params):
    """
    Uproszczony ODE: przepływ izoentropowy w dyszy zmiennego przekroju.

    Równania zależą tylko od A(x) i dA/dx — brak tarcia i ciepła.
    Gamma jest interpolowana lokalnie w punkcie x (PCHIP).
    """
    N, P, T = Y
    gamma = float(params['gamma_interp'](x))

    A = params['A_func'](x)
    dA = params['dA_func'](x)

    # Równanie na kwadrat liczby Macha
    dNdx = -(N / (1 - N)) * ((2 + (gamma - 1) * N) / A) * dA
    # Równanie na ciśnienie
    dPdx = (P / (1 - N)) * ((gamma * N) / A) * dA
    # Równanie na temperaturę
    dTdx = (T / (1 - N)) * (((gamma - 1) * N) / A) * dA

    return [dNdx, dPdx, dTdx]


def my_nozzle_ode_full2(x, Y, params, dQdx_fun, dFdx_fun):
    """
    Pełny ODE: przepływ z uwzględnieniem wymiany ciepła (dQ/dx) i tarcia (dF/dx).

    dQdx_fun, dFdx_fun — interpolanty obliczone wcześniej z profilu chłodzenia.
    Gamma, Cp i Rs są interpolowane lokalnie w punkcie x (PCHIP).
    """
    N_2, P_2, T_2 = Y
    gamma = float(params['gamma_interp'](x))
    Cp    = float(params['Cpcg_interp'](x))
    Rs    = float(params['Rs_interp'](x))

    A = params['A_func'](x)
    dA = params['dA_func'](x)
    dq = -dQdx_fun(x)       # minus, bo ciepło oddawane z gazu
    df = dFdx_fun(x)

    tolN = 1e-10
    if not (np.isfinite(dq) and np.isfinite(df) and A > 0
            and T_2 > 0 and abs(1 - N_2) > tolN):
        raise ValueError(
            f"Bad domain at x={x}: A={A}, T={T_2}, |1-N2|={abs(1-N_2)}, dq={dq}, df={df}"
        )

    dNdx = (N_2 / (1 - N_2)) * (
        ((1 + gamma * N_2) / (Cp * T_2)) * dq
        + ((2 + (gamma - 1) * N_2) / (Rs * T_2)) * df
        - ((2 + (gamma - 1) * N_2) / A) * dA
    )

    dPdx = (
        -(P_2 / (1 - N_2)) * ((gamma * N_2) / (Cp * T_2)) * dq
        - (P_2 / (1 - N_2)) * ((1 + (gamma - 1) * N_2) / (Rs * T_2)) * df
        + (P_2 / (1 - N_2)) * ((gamma * N_2) / A) * dA
    )

    dTdx = (
        (T_2 / (1 - N_2)) * ((1 - gamma * N_2) / (Cp * T_2)) * dq
        - (T_2 / (1 - N_2)) * (((gamma - 1) * N_2) / (Rs * T_2)) * df
        + (T_2 / (1 - N_2)) * (((gamma - 1) * N_2) / A) * dA
    )

    return [dNdx, dPdx, dTdx]


def event_N1(x, Y, *args):
    """
    Zdarzenie: N = 1  ⟹  Mach = 1  (bariera dźwięku).

    Solver ODE zatrzyma się, gdy ta funkcja zwróci 0.
    """
    return Y[0] - 1.0


# Atrybuty wymagane przez scipy.integrate.solve_ivp
event_N1.terminal = True    # zatrzymaj solver po wykryciu
event_N1.direction = 1      # wykryj tylko gdy N rośnie przez 1