Symulacja 6 stopni swobody: składamy samolot

0 Prolog: niebieskie pudło Edwina Linka

Pod koniec lat dwudziestych Edwin Link, syn producenta organów, zbudował z mieszków organowych i pomp kabinę, która przechylała się i pochylała w odpowiedzi na drążek. Lotnicy traktowali ją jak atrakcję wesołego miasteczka, dopóki w 1934 roku armia amerykańska, po serii wypadków pilotów wożących pocztę w chmurach, nie kupiła pierwszych trenażerów do nauki lotu według przyrządów. W czasie wojny przez „niebieskie pudła” przeszły setki tysięcy pilotów.

Trenażer Linka nie liczył żadnych równań: odtwarzał tylko, z grubsza, odczucia i wskazania przyrządów. Dzisiejszy symulator to program, który całkuje równania ruchu bryły sztywnej z siłami aerodynamicznymi i jest podstawowym narzędziem projektanta: na nim testuje się autopilota, zanim samolot oderwie się od ziemi. Ta lekcja buduje taki program od zera, z elementów, które już znamy.

Pytanie na start. Model liniowy z lekcji 4–5 i symulacja nieliniowa dają dla małego podmuchu prawie ten sam przebieg. Po co więc symulacja nieliniowa? Wymień trzy sytuacje, w których model liniowy na pewno zawiedzie. Odpowiedź w części 4.

1 Równania ruchu: 13 stanów

Samolot to bryła sztywna (lekcja 2 Bryła w 3D: tensor bezwładności, równania Eulera i oś pośrednia) o masie $m$ i tensorze bezwładności $I$. Jej stan opisuje 13 liczb:

Równania ruchu to cztery wektorowe równania pierwszego rzędu:

Pierwsze równanie to druga zasada dynamiki zapisana w obracającym się układzie ciała: człon $-\boldsymbol\omega\times\mathbf v_b$ pojawia się dlatego, że osie, w których mierzymy prędkość, same się obracają (lekcja 2 Bryła w 3D: tensor bezwładności, równania Eulera i oś pośrednia). Grawitacja jest prosta w NED, więc przenosimy ją do osi ciała macierzą $R^{\mathsf T}$. Drugie równanie to równania Eulera bryły sztywnej. Trzecie to kinematyka kwaternionu z lekcji 1, czwarte zamienia prędkość w osiach ciała na ruch po mapie.

Dlaczego kwaternion, a nie kąty Eulera? Kąty Eulera mają osobliwość przy $\theta = \pm 90^\circ$ (lekcja 1): symulacja pętli albo pionowego wznoszenia by się wysypała. Kwaternion nie ma osobliwości, kosztem jednej nadmiarowej liczby i więzu $|\mathbf q| = 1$, który po każdym kroku przywracamy, dzieląc przez normę. Kąty $\psi,\theta,\varphi$ liczymy z kwaternionu tylko do wyświetlania i dla autopilota.

2 Model aerodynamiczny: współczynniki, przeciągnięcie, opóźnienie

Całą fizykę samolotu zawierają siły $\mathbf F_a$ i momenty $\mathbf M_a$. Symulator liczy je w każdym kroku z bieżącego stanu: najpierw prędkość $V = |\mathbf v_b|$, kąt natarcia $\alpha = \operatorname{atan2}(w, u)$, kąt ślizgu $\beta = \arcsin(v/V)$ i ciśnienie dynamiczne $\bar q = \tfrac12\rho(h)V^2$ z gęstością z atmosfery wzorcowej, a potem współczynniki. Używamy dokładnie samolotu z lekcji 4 Stateczność statyczna i trym: dlaczego samolot ma ogon–7 Ruch boczny: przechylenie, spirala i holendrowanie:

Trzy rzeczy w tym modelu są nowe w porównaniu z modelami liniowymi. Po pierwsze, przeciągnięcie: powyżej $14^\circ$ siła nośna skrzydła maleje, a opór gwałtownie rośnie. To bardzo uproszczony model (prawdziwe przeciągnięcie ma histerezę, zależy od $\dot\alpha$, jest niesymetryczne i potrafi przejść w korkociąg), ale wystarcza, żeby zobaczyć granicę obszaru, w którym działa linearyzacja. Po drugie, prędkość wpływa przez $\bar q$ i $V$ w mianownikach dokładnie, a nie przez pochodne $X_u$, $Z_u$. Po trzecie, ruch podłużny i boczny są w jednych równaniach: w przechyleniu siła nośna nie jest pionowa, a to sprzęga oba ruchy.

Człon z $\dot\alpha$ wymaga uwagi. Strumień odchylony przez skrzydło dociera do usterzenia po czasie $l_H/V$ (lekcja 5 Lot w płaszczyźnie symetrii: model podłużny i pochodne stateczności: pochodna $M_{\dot\alpha}$). W symulatorze $\dot\alpha$ zależy od $\dot w$, a $\dot w$ od sił, które zależą od $\dot\alpha$: pętla algebraiczna. Najprostsze wyjście, stosowane w wielu symulatorach, to wziąć $\dot\alpha$ z poprzedniego kroku, z różnicy wstecznej $(\alpha_k - \alpha_{k-1})/\Delta t$. Przy kroku 0,01 s opóźnienie o jeden krok jest pomijalne wobec stałej czasowej $l_H/V = 0{,}08$ s.

3 Architektura symulatora: trym, całkowanie, linearyzacja

Symulator ma cztery części, które warto pisać osobno i testować osobno (tak radzi też Zipfel):

Trym daje dla naszego samolotu przy 55 m/s: $\alpha = 4{,}261^\circ$, $\delta_e = -2{,}639^\circ$, $T = 1123$ N. Model liniowy z lekcji 4 dawał $4{,}292^\circ$ i $-2{,}672^\circ$. Różnica jest prawdziwa i ma prosty powód: w symulatorze ciąg działa wzdłuż osi kadłuba, więc ma składową pionową $T\sin\alpha = 83$ N, która odciąża skrzydło (przykład 6.1).

Ostatni element to linearyzacja numeryczna: macierz $A = \partial f/\partial\mathbf x$ w punkcie trymu z różnic centralnych, kolumna po kolumnie. To pozwala porównać symulator z modelami liniowymi z lekcji 4–7. Jest jedna pułapka: $f$ zależy też od $\dot\alpha$. Linearyzacja, która to pomija, dostaje dla modu krótkookresowego $-3{,}81 \pm 5{,}05i$ ($\zeta = 0{,}60$), a lekcja 4 miała $-4{,}92 \pm 4{,}11i$ ($\zeta = 0{,}77$). Brakuje tłumienia od $M_{\dot\alpha}$. Poprawnie: $\dot{\mathbf x} = A\mathbf x + \mathbf b\,\dot\alpha$, gdzie $\dot\alpha = \mathbf c^{\mathsf T}\dot{\mathbf x}$, więc

Z tym poprawionym wzorem: mod krótkookresowy $-4{,}80 \pm 4{,}05i$, fugoida $-0{,}0166 \pm 0{,}208i$; lekcja 4: $-4{,}92 \pm 4{,}11i$ i $-0{,}0165 \pm 0{,}206i$. Fugoida zgadza się prawie dokładnie, mod krótkookresowy do 2–3%; resztę wyjaśniają uproszczenia lekcji 4 (np. opór liczony z całkowitego $C_L$, inny kierunek ciągu).

Druga pułapka linearyzacji dotyczy kwaternionu. Żeby porównać bieguny z lekcjami 4–7, sprowadzamy 13 stanów do zmiennych fizycznych $(u, w, q, \theta)$ i $(v, p, r, \varphi)$: $\mathbf x = \mathbf x_0 + T\mathbf z$, gdzie kolumny $T$ to pochodne stanu po $\theta$, $\varphi$ itd., i $A_z = T^{+}AT$. Wybranie z macierzy samych wierszy i kolumn $q_1$ albo $q_2$ jest błędem: przy $\theta_0 \ne 0$ zmiana przechylenia porusza kilka składowych kwaternionu naraz. Co więcej, w ruchu bocznym trzeba dołączyć do $T$ także kierunek $\psi$ (i potem go odrzucić, bo kurs nie wpływa na dynamikę): przy $\theta_0 \ne 0$ kierunki „zmiana $\varphi$” i „zmiana $\psi$” w przestrzeni kwaternionów nie są prostopadłe, a rzut na samo $\varphi$ myli zmianę kursu ze zmianą przechylenia. Pierwsza wersja skryptu do tej lekcji popełniła ten błąd i dostała mod spiralny trzy razy za szybki. Najpewniejszy test: wartości własne całej macierzy $13\times 13$ (pięć z nich to zera od położenia, kursu i normy kwaternionu) muszą zawierać wszystkie bieguny modeli zredukowanych.

import numpy as np

def f(x, u, ald):
    """Prawa strona: x = [u v w p q r q0 q1 q2 q3 xN yE zD], u = [de da dr T]."""
    vb, w, q, pos = x[0:3], x[3:6], x[6:10], x[10:13]
    F, M = aero(x, u, ald)                    # model z części 2, osie ciała
    R = q2R(q)                                # ciało -> NED (lekcja 1)
    vdot = (F + R.T @ [0, 0, m*g]) / m - np.cross(w, vb)
    wdot = np.linalg.solve(I, M - np.cross(w, I @ w))
    qdot = 0.5 * qmul(q, np.r_[0, w])
    return np.r_[vdot, wdot, qdot, R @ vb]

def rk4(x, u, dt, ald):
    k1 = f(x, u, ald); k2 = f(x + dt/2*k1, u, ald)
    k3 = f(x + dt/2*k2, u, ald); k4 = f(x + dt*k3, u, ald)
    x = x + dt/6*(k1 + 2*k2 + 2*k3 + k4)
    x[6:10] /= np.linalg.norm(x[6:10])        # więz |q| = 1
    return x

def simulate(x0, ctrl, TF, dt=0.01):
    x, st, out = x0.copy(), {}, []
    al_prev, ald = np.arctan2(x[2], x[0]), 0.0
    for k in range(int(round(TF/dt)) + 1):
        u = ctrl(k*dt, x, st, dt)             # autopilot: próbkowanie co krok
        out.append(np.r_[k*dt, x, u])
        x = rk4(x, u, dt, ald)
        al = np.arctan2(x[2], x[0]); ald = (al - al_prev)/dt; al_prev = al
    return np.array(out)

4 Walidacja i manewry nieliniowe

Symulator, który nie został sprawdzony, to generator przekonujących wykresów. Kolejność sprawdzania: najpierw testy bez aerodynamiki (bryła swobodna zachowuje energię i moment pędu, lekcja 2), potem trym ($f = 0$ w trymie; u nas reszta $2\cdot10^{-15}$), potem zgodność z modelem liniowym dla małych zaburzeń, a dopiero potem manewry, których model liniowy nie opisuje.

Mały podmuch. Pionowy podmuch zwiększa kąt natarcia skokowo o $2^\circ$. Symulacja nieliniowa i model liniowy z lekcji 4 różnią się najwyżej o $0{,}010^\circ$. Przy $8^\circ$ różnica rośnie do $0{,}038^\circ$, prawie dokładnie czterokrotnie. To ważna wskazówka: gdyby różnica pochodziła z nieliniowości, rosłaby jak kwadrat amplitudy (szesnastokrotnie). Rośnie liniowo, więc pochodzi z drobnych różnic między modelami (ciąg wzdłuż kadłuba, opór), a nie z nieliniowości.

Wyrwanie. Ster w górę o $\Delta\delta_e$ przez 3 s, potem z powrotem. Przy $4^\circ$: największy kąt natarcia $7{,}1^\circ$, prędkość spada do 40,9 m/s, wysokość rośnie o 73 m. Przy $8^\circ$: $9{,}9^\circ$, 27,6 m/s, 116 m. Przy $12^\circ$: $19{,}4^\circ$, czyli przeciągnięcie, prędkość 15,8 m/s. Model liniowy (z $u$ jako małym zaburzeniem) już przy $4^\circ$ rozjeżdża się w prędkości: spadek o 14 m/s to nie jest małe zaburzenie prędkości 55 m/s.

Ciasny zakręt. Autopilot z lekcji 9 Autopilot: pętle wewnętrzne i zewnętrzne zadaje przechylenie i trzyma wysokość, z kompensacją $q$ w zakręcie: w ustalonym zakręcie samolot obraca się w pochyleniu z $q = (g/V)\tan\varphi\sin\varphi$, i pętla tłumienia nie może z tym walczyć. Przy $30^\circ$ wszystko jest łagodne: samolot traci najwyżej 19,5 m. Przy $60^\circ$ przeciążenie wynosi $1/\cos 60^\circ = 2$ g; siła nośna musi się podwoić, więc opór indukowany rośnie czterokrotnie. Z automatem ciągu: ciąg chwilowo do 2458 N, strata wysokości najwyżej 39 m; całka automatu ciągu ustala się powoli i dopiero po ok. 200 s zakręt jest naprawdę ustalony: ciąg 1790 N, kąt natarcia $8{,}5^\circ$. Bez automatu ciągu ten sam autopilot ciągnie coraz mocniej, aż do kąta natarcia $20^\circ$, za przeciągnięciem: samolot traci 610 m. Modele liniowe, osobne dla ruchu podłużnego i bocznego, nie widzą tego sprzężenia w ogóle.

Puszczone stery. Samolot wytrymowany, przechylony o $20^\circ$, stery puszczone. Spirala z lekcji 7 Ruch boczny: przechylenie, spirala i holendrowanie jest stateczna, więc przechylenie maleje: w symulacji do $16{,}7^\circ$ po 10 s, $14{,}8^\circ$ po 30 s i $11{,}7^\circ$ po 60 s. Model liniowy z lekcji 7 daje $15{,}8^\circ$, $10{,}8^\circ$ i $6{,}1^\circ$. Linearyzacja symulatora wyjaśnia różnicę: mod przechylenia $-12{,}46$/s (lekcja 7: $-12{,}65$) i holendrowanie $-0{,}68 \pm 3{,}05i$ ($-0{,}58 \pm 3{,}04i$) zgadzają się nieźle, ale mod spiralny ma $-0{,}0067$/s, czyli stałą czasową 149 s zamiast 52 s. Przyczyna: lekcja 7 liczy w osiach stateczności (z zerowym $I_{xz}$ w tych osiach), a symulator przykłada te same pochodne do prędkości kątowych w osiach ciała, nachylonych o $\alpha_0 = 4{,}3^\circ$, z zerowym $I_{xz}$ w osiach ciała. To dwa trochę różne samoloty. Sprawdzenie numeryczne: gdy w symulatorze liczyć momenty boczne z prędkości kątowych w osiach stateczności i przeliczyć tensor bezwładności ($I_{xz} = -102$ kg·m² w osiach ciała), linearyzacja daje $-12{,}65$, $-0{,}58 \pm 3{,}04i$ i $-0{,}0193$, czyli dokładnie lekcję 7. Spirala, która jest różnicą dwóch małych iloczynów pochodnych (lekcja 7), wychwytuje taką niespójność pierwsza.

Odpowiedź na pytanie z prologu: model liniowy zawiedzie wszędzie tam, gdzie zaburzenie nie jest małe (prędkość zmienia się o kilkadziesiąt procent, kąt natarcia dochodzi do przeciągnięcia), gdzie ruch podłużny i boczny się sprzęgają (zakręt z dużym przechyleniem, korkociąg) i gdzie działają ograniczenia (nasycenie sterów i ciągu, logika trybów autopilota). Do tego dochodzą efekty, których nie ma w modelu w ogóle, ale to już kwestia modelu, a nie linearyzacji.

5 Monte Carlo i projekt

Jeden przebieg symulatora odpowiada jednemu samolotowi w jednych warunkach. Prawdziwy samolot lata z różnym załadunkiem, a środek ciężkości przesuwa się z paliwem i pasażerami. Analiza Monte Carlo losuje parametry z ich rozkładów, dla każdego losowania trymuje i linearyzuje (albo symuluje manewr) i patrzy na rozrzut wyników. W naszym przykładzie: położenie środka ciężkości $h \sim N(0{,}30;\ 0{,}04^2)$, masa jednostajnie od 900 do 1250 kg, bezwładności skalowane z masą, 400 losowań (ziarno 7).

Najgorszy przypadek dla tłumienia krótkookresowego to przedni środek ciężkości ($h = 0{,}20$) z dużą masą: duża sztywność $M_\alpha$ przy tym samym tłumieniu $M_q$ podnosi częstość, a nie tłumienie. Wynik jest spójny z lekcją 6 (wymagania MIL dla $\zeta_{sp}$ od 0,35 są spełnione z dużym zapasem). Monte Carlo nie dodaje fizyki; dodaje odpowiedź na pytanie „czy w całym zakresie”, którego jeden przebieg nie daje. Ma dwa ograniczenia: wynik jest tylko tak dobry, jak przyjęte rozkłady, a rzadkie zdarzenia (np. z prawdopodobieństwem $10^{-4}$) wymagają ogromnej liczby losowań.

Projekt: rozbuduj symulator

1. Test integratora. Wyłącz aerodynamikę i grawitację, nadaj bryle obrót wokół osi pośredniej (lekcja 2). Sprawdź zachowanie energii i momentu pędu dla kroków 0,01 i 0,1 s.
2. Trym w zakręcie. Rozszerz trym o $\varphi$, $\delta_a$, $\delta_r$, $\beta$ i warunek $\dot\psi =$ const. Porównaj ciąg z symulacją z autopilotem.
3. Wiatr. Dodaj prędkość wiatru w NED: aerodynamika ma liczyć siły od prędkości względem powietrza $\mathbf v_b - R^{\mathsf T}\mathbf w_{NED}$. Sprawdź, że stały wiatr tylko przesuwa tor.
4. Turbulencja. Dodaj podmuchy z filtru pierwszego rzędu sterowanego szumem (lekcja 9) i zbadaj Monte Carlo rozrzut wysokości w zakręcie $30^\circ$.
5. Tłumiki i autopilot. Dołącz tłumik odchylenia z lekcji 8 i sprawdź ślizg w zakręcie $60^\circ$.

6 Metoda i przykłady rozwiązane

7 Błędne intuicje, których trzeba się oduczyć

„Symulacja nieliniowa jest dokładniejsza niż model liniowy.”

Dokładniej całkuje równania, ale model sił może być tak samo niepewny. Dla małych zaburzeń oba modele dają to samo; różnica w spirali (część 4) pochodzi z drobnej niespójności modelu, a nie z nieliniowości.

„Skoro symulacja daje gładkie wykresy, jest poprawna.”

Błąd znaku albo zamienione osie dają gładkie, przekonujące i całkowicie fałszywe wyniki. Poprawność pokazują testy: trym, bilanse energii, zgodność z modelem liniowym.

„Mniejszy krok zawsze daje lepszy wynik.”

Po pewnym progu błąd całkowania jest dużo mniejszy niż błąd modelu i niczego już nie zmienia, a przy bardzo małym kroku rosną błędy zaokrągleń i czas obliczeń. Krok dobiera się do najszybszych modów i do częstotliwości regulatorów.

„Kwaternion nie wymaga normalizacji, bo RK4 jest dokładny.”

Każdy integrator powoli wyprowadza kwaternion z jednostkowej sfery. Normalizacja po kroku jest tania i konieczna.

„Monte Carlo znajdzie każdy groźny przypadek.”

Znajdzie typowe przypadki z przyjętych rozkładów. Zdarzenia rzadkie wymagają ogromnej liczby losowań albo analizy najgorszego przypadku; zdarzenia spoza modelu nie pojawią się wcale.

8 Zadania do przemyślenia

Te zadania wymagają więcej myślenia niż liczenia.

9 Sprawdź się

Odpowiedz w myślach, zanim rozwiniesz odpowiedź.

10 Karta wzorów i co dalej

Do przeczytania

• Obowiązkowo: B. L. Stevens, F. L. Lewis, E. N. Johnson, Aircraft Control and Simulation (Wiley, wyd. 3), rozdz. 2–3: model samolotu, trym, linearyzacja i symulacja; wzorcowy przykład w kodzie.
• Dla chętnych: P. H. Zipfel, Modeling and Simulation of Aerospace Vehicle Dynamics (AIAA): części o symulacjach 5 i 6 DOF, w tym z obracającą się, kulistą Ziemią; tegoż autora Building Aerospace Simulations in C++ o architekturze symulatorów.

Samolot jest złożony: od kwaternionu po autopilota, w jednym programie, który umie powiedzieć, kiedy prostsze modele kłamią. Kolejna część zamienia skrzydła na silnik: rakieta nie ma siły nośnej, która by ją stabilizowała, ma za to ciąg, który można obracać. Ta sama bryła sztywna i te same pętle sterowania, ale z masą, która znika w trakcie lotu.