OpenEngine/geometry.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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
"""
geometry.py — Budowa geometrii konturu dyszy de Lavala.

Tworzy profil dyszy złożony z odcinków kołowych, prostych i krzywej Béziera.
Zwraca siatkę punktów (x, r), pola przekrojów A(x) i ich pochodne dA/dx,
oraz interpolanty PCHIP gotowe do użycia w solverze ODE.

Odpowiada funkcji buildNozzleGeometry() z MATLAB-a.
"""

import numpy as np
from scipy.interpolate import PchipInterpolator


def build_nozzle_geometry(R_param=0.01878, E_r=5, n_grid=100):
    """
    Buduje geometrię konturu dyszy rakietowej.

    Parameters
    ----------
    R_param : float — promień gardzieli [m] (domyślnie 0.01878)
    E_r : float — stosunek ekspansji A_exit/A_throat (domyślnie 5)
    n_grid : int — liczba punktów siatki (domyślnie 100)

    Returns
    -------
    x_grid : np.ndarray  — pozycje osiowe [m]
    r_grid : np.ndarray  — promień konturu w każdym punkcie [m]
    A_grid : np.ndarray  — pole przekroju A = pi * r^2 [m^2]
    dA_grid_dx : np.ndarray — pochodna dA/dx [m]
    A_interp : PchipInterpolator — interpolant A(x)
    dA_interp : PchipInterpolator — interpolant dA/dx(x)
    """

    # --- Punkty charakterystyczne konturu (sekcja zbieżna) ---
    # Obliczamy współrzędne kluczowych punktów
    # na podstawie promieni zaokrągleń i kątów (norma TRD / Rao).
    X1 = 1.5 * R_param * np.cos(np.radians(-120))
    Y1 = 1.5 * R_param * np.sin(np.radians(-120)) + 2.5 * R_param

    # Nachylenie stycznej
    m = X1 / np.sqrt((1.5 * R_param) ** 2 - X1 ** 2)

    Y3 = -0.0306 + 0.07265 / np.sqrt(1 + m ** 2)
    X3 = X1 + (Y3 - Y1) / m
    Xc = X3 - (-m * 0.07265) / np.sqrt(1 + m ** 2)

    # --- Punkt N (początek krzywej Béziera – sekcja rozbieżna) ---
    X2 = 0.382 * R_param * np.cos(np.radians(-68))
    Y2 = 0.382 * R_param * np.sin(np.radians(-68)) + 1.382 * R_param

    Nx, Ny = X2, Y2
    m1 = np.tan(np.radians(22))
    m2 = np.tan(np.radians(12))
    c1 = Ny - m1 * Nx

    # Punkt wylotowy E
    Ey = np.sqrt(E_r) * R_param
    Ex = 0.8 * (((np.sqrt(E_r) - 1) * R_param) / np.tan(np.radians(15)))
    c2 = Ey - m2 * Ex

    # Punkt kontrolny Q krzywej Béziera
    Qx = (c2 - c1) / (m1 - m2)
    Qy = (m1 * c2 - m2 * c1) / (m1 - m2)

    # --- Budowa 6 segmentów konturu ---
    N_pts = 500  # punktów na segment (duża gęstość dla dokładności)

    # Segment 1: komora spalania (stały promień)
    x1 = np.linspace(-0.14262, Xc, N_pts)
    y1 = 0.04205 * np.ones(N_pts)

    # Segment 2: łuk okrągły (zbieżna, duży R)
    x2 = np.linspace(Xc, X3, N_pts)
    y2 = -0.0306 + np.sqrt(0.07265 ** 2 - (x2 - Xc) ** 2)

    # Segment 3: odcinek prosty (zbieżna)
    x3 = np.linspace(X3, X1, N_pts)
    y3 = m * (x3 - X1) + Y1

    # Segment 4: łuk okrągły tuż przed gardzielą
    x4 = np.linspace(X1, 0, N_pts)
    y4 = -np.sqrt((1.5 * R_param) ** 2 - x4 ** 2) + 2.5 * R_param

    # Segment 5: łuk okrągły tuż za gardzielą
    x5 = np.linspace(0, X2, N_pts)
    y5 = -np.sqrt((0.382 * R_param) ** 2 - x5 ** 2) + 1.382 * R_param

    # Segment 6: krzywa Béziera (rozbieżna, profil Rao)
    t = np.linspace(0, 1, N_pts)
    x6 = (1 - t) ** 2 * Nx + 2 * (1 - t) * t * Qx + t ** 2 * Ex
    y6 = (1 - t) ** 2 * Ny + 2 * (1 - t) * t * Qy + t ** 2 * Ey

    # --- Sklejenie segmentów (bez duplikatów na złączeniach) ---
    X_raw = np.concatenate([x1, x2[1:], x3[1:], x4[1:], x5[1:], x6[1:]])
    Y_raw = np.concatenate([y1, y2[1:], y3[1:], y4[1:], y5[1:], y6[1:]])

    # Przesuwamy tak, żeby x zaczynało się od 0
    X_raw = X_raw - X_raw[0]

    # --- Parametryzacja po długości łuku → równomierna siatka 100 punktów ---
    S = np.zeros(len(X_raw))
    S[1:] = np.cumsum(np.sqrt(np.diff(X_raw) ** 2 + np.diff(Y_raw) ** 2))
    S_100 = np.linspace(0, S[-1], n_grid)

    x_grid = np.interp(S_100, S, X_raw)
    r_grid = np.interp(S_100, S, Y_raw)

    # --- Obliczenie pola przekroju i jego pochodnej ---
    A_grid = np.pi * r_grid ** 2
    dr_dx = np.gradient(r_grid, x_grid)          # dr/dx (numerycznie)
    dA_grid_dx = 2 * np.pi * r_grid * dr_dx      # dA/dx = 2*pi*r * dr/dx

    # --- Interpolanty PCHIP (gładkie, monotoniczne) ---
    A_interp = PchipInterpolator(x_grid, A_grid, extrapolate=True)
    dA_interp = PchipInterpolator(x_grid, dA_grid_dx, extrapolate=True)

    return x_grid, r_grid, A_grid, dA_grid_dx, A_interp, dA_interp