Zakręt Macierewicza cz. 1

Raport Zespłou Parlamentarnego kierowanego przez Antoniego Macierewicza, noszący nazwę "28 miesięcy po Smoleńsku" w części technicznej nieudolnie próbuje zaprezentować alternatywny przebieg zdarzeń dnia 10.04.2010. Proponowaną sekwencję wydarzeń najlepiej prezentuje rysunek zamieszczony w raporcie dwukrotnie na stronach (wg numeracji pdf) 132 i 142 (opisany jako rysunek 21c).


Spróbujmy rozkodować sekwencję zdarzeń z tego rysunku.
    1. Pierwsza sylwetka samolotu z prawej strony, z antenami GPS mniej więcej nad
        punktem TAWS 38, z całymi skrzydłami, bez przechylenia, co zresztą potwierdza
        napis ROLL ~0.
    2. Druga sylwetka od prawej, lewe skrzydło pozbawione końcówki (przypuszczalnie po
        pierwszym "wybuchu"), narastający przechył w lewo.
    3. Trzecia sylwetka od prawej, samolot nad tak zwanym "punktem zamrożenia FMS"
        (osobiście nazywam go "FMS2 status" ponieważ jego pozycja znajduje się w statusie
        FMS drugiego pilota), przechylenie urosło do około 67 stopni (ROLL ~67).
    4. Samolot w kawałkach, przedstawienie uderzenia statecznikiem poziomym
        w miejsce, w którym leżał, sugestia odpadnięcia.
    5. Końcowy rozpad.


Na pierwszy rzut oka powyższa sekwencja nie zdradza swoich ukrytych niedostatków, dopiero zestawienie jej z analizami z mojej poprzedniej notki, nasuwa pewne wątpliwości, czy taki scenariusz jest możliwy. Problem leży w procesie zmiany toru lotu samolotu, który łatwo poddaje się szacunkom obliczeniowym.



Przechylenie

Do wykonania oszacowań zmiany toru lotu muszę stworzyć model zmiany przechylenia samolotu dostosowany do scenariusza przygotowanego przez Zespół Parlamentarny. Nie jest to trudne zadanie. W chwili zero (nad TAWS38) samolot ma mieć przechylenie 0 stopni. Nad czerwonym kwadracikiem "FMS2 status" przechylenie ma wynosić 67 stopni. Punkty te dzieli nieco ponad 140 metrów, na pokonanie których samolot potrzebuje około dwóch sekund. Przyjmiemy, że po dalszej sekundzie przechylenie wzrasta do 90 stopni. Czyli mamy taką sekwencję:

Do swoich celów użyłem funkcji:

(1)   roll( t ) = 41 * t^0,71   [wynik w stopniach]

gdzie t, to czas, a znaczek ^ oznacza potęgowanie.

Na podstawie powyższego wzoru otrzymałem następujący wykres przechylenia:




Przeciążenie pionowe

Drugim, ważnym parametrem jest przeciążenie pionowe.  Wezmę je z wykresu w raporcie MAK. Zespół z niego wywodzi dwa "potężne wstrząsy", więc z definicji musi być oryginalny i godny zaufania.

Na powyższym wykresie chwila zero (czyli moment przelotu nad TAWS38) wypada o czasie 41:02 plus minus pół sekundy. Do dalszych rozważań przyjmę dokładnie moment 41:02. Ponieważ w odpowiedzi na mój poprzedni tekst jeden z ekspertów zespołu zawarł pewne uwagi na temat szarej strefy (fragment wykresu narysowany szarą linią w okresie, gdy rejestrator nie zapisał prawidłowych danych), pójdę mu na rękę i przez pierwsze pół sekundy przyjmę przeciążenie 1,5g (dlaczego tyle, wyjaśni się w innym tekście).

Wykorzystany fragment wykresu:

Wartości w przedziale 2,5 do 3 sekundy zostały dodane w celu możliwości poprowadzenia wykresów praktycznie do miejsca upadku samolotu.



Tor lotu

Uzbrojony w powyższe dane mogę przeliczyć uproszczoną trajektorię poziomą zakrętu. Ze względu na małe kąty natarcia po 0,5 sekundy mogę zaniedbać wpływ ciągu silnika. Siłę działającą na samolot, a zanotowaną przez akcelerometr rozłożę na dwie siły: pionową i poziomą.

Zmianę toru lotu w bok powoduje składowa pozioma, której wartość wynosi:

(2)   siła_nośna * sin( przechylenie )

Obliczenia prowadzone są metodą zbliżoną do opisanej w tekście "Nie daj się robić w balona ekspertom Zespołu".

Ponieważ wartość przyśpieszenia zmienia się w czasie, do wyliczenia toru lotu nie mogę użyć wprost wzorów na ruch jednostajnie przyśpieszony. Zamiast tego dzielę czas na małe kroczki. W jednym takim kroczku (licząc się z faktem, że to tylko dość dokładne przybliżenie) stosuję znane wzory:

(3)   Δv = a * Δt ; v = v + Δv
(4)   s = s + v * Δt

W każdym kroczku, obliczam chwilowe przyśpieszenie boczne, wzory (1) oraz (2), by następnie ze wzoru (3) obliczyć zmianę prędkości, a ze wzoru (4) odległość od przedłużenia pierwotnego toru lotu. Powtórzenie tych kroczków dla wybranego przedziału czasu, umożliwia wyliczenie trajektorii poziomej.

Należy pamiętać, że z powodu nieuwzględnienia wymiaru wzdłużnego, przy oddalaniu się od pierwotnego toru lotu metoda ta powoli traci na dokładności, jednak niedokładność ta jest na tyle mała (nie osiągnie nawet jednego metra), że można ją pominąć.


Dopisek:

Ponieważ kilka osób uczepiło się tej niedokładności wynikającej z obrotu wektora siły bocznej (wraz ze zmianą kursu) tak, jakby różnica w wynikach poniżej miała jakieś istotne znaczenie, pokażę jakie modyfikacje należy wykonać, by uzyskać bardziej dokładne wyniki. Można to zrobić wprowadzając pojęcie kierunku wektora prędkości, w którym będziemy kumulować zmiany kierunku.

(5)   Δdir = asin( Δv / vgs ) ;  dir = dir +Δdir

gdzie vgs to prędkość postępowa (całkowita), można trzymać stałą lub liczyć za pomocą Twierdzenia Pitagorasa z vx i vy. Zamiast równania (3) i (4) należy zastosować

(6)   Δvx = -sin( dir ) * a * Δt ; Δvy = cos( dir ) * a *  Δt ; vx = vx + Δvx ; vy = vy + Δvy
(7)   sx = sx + vx * Δt ; sy = sy + vy * Δt

################################################
# Weryfikacja poprawności algorytmu.
# Symulacja ruchu po okręgu ze stałą prędkością 75 m/s
# i przyśpieszeniem dośrokowym 40 m/s^2
vgs = 75   # prędkość [m/s]
ay = 40    # przyśpieszenie boczne [m/s^2]
vx = vgs
vy = 0
x = 0
y = 0
dir = 0
sps = 1000

for i = 0 to 15 * sps
    t = i / sps
    dir = dir + asin( ay / sps / sqrt( vx^2 + vy^2 ) )
    vy = vy + cos( dir ) * ay / sps
    vx = vx - sin( dir ) * ay / sps
    y = y + vy / sps
    x = x + vx / sps

    if i * 10 mod sps = 0
        print( "{d};{d}{cr}", x, y )
    end if
end for
################################################


Program ten produkuje serię punktów trajektorii, które po zaimportowaniu do arkusza kalkulacyjnego pozwalają stworzyć taki wykres:

Alternatywnie można liczyć rachunkiem wektorowym - wtedy nie są potrzebne żadne funkcje trygonometryczne.

(8)  v = v + rot90( norm( v ) ) * ( a * Δt )
(9)  p = p + v * Δt

Należy pamiętać, że v i p są wielkościami wektorowymi i odpowiednio interpretować operatory.



Tor lotu na tle mapy satelitarnej

Uwaga,na wykresie oś czasu jest w pionie. Wybrałem taką orientację z powodu ograniczeń co do szerokości w blogach Salonu24.



Błękitna linia pokazuje tor lotu samolotu w płaszczyźnie horyzontalnej w oparciu o powyższe dane.



Dla zrozumienia wykresu, najlepiej wyobrazić sobie, że w chwili zero samolot znajduje się w miejscu TAWS 38 (a konkretnie jego anteny GPS), po czym przesuwa się w górę i jednocześnie powoli skręca w lewo. Proporcje rozmiarów są zachowane.



Zarówno na tym wykresie, jak i na następnych dla uproszczenia pomijam wpływ yaw (odchylenie osi podłużnej samolotu od kierunku lotu), który w niektórych przypadkach może powodować oddalenie anten GPS od toru lotu o około 5 metrów (akurat przy przechyleniu 67 stopni i więcej, yaw raczej przesuwa anteny w pionie niż w poziomie).




Cóż można powiedzieć o tym torze lotu?

W pierwszym rządzie widać jak bardzo słabo skręca w lewo, przez co nie ma szans zbliżyć się na sensowną odległość zarówno do FMS2 status jak i last GPS. Co gorsza, jego przedłużenie nie prowadzi na miejsce katastrofy, tylko około 40 metrów obok (czyżby hipoteza 2M Lite? :)

Wracając do zacytowanego rysunku z raportu Macierewicza, możemy tam zauważyć sylwetkę samolotu narysowaną dokładnie nad punktem FMS2 status. Cóż narysować można, nie ma to jednak nic wspólnego z fizyką lotu samolotu.

W następnej notce spróbuję pobawić się w adwokata diabła i pozmieniać nieco założenia Zespołu Parlamentarnego, aby zbadać przy jakich założeniach można osiągnąć rzeczywiste miejsce upadku samolotu.


Na razie tytułem podsumowania, konsekwencje zaprezentowanego scenariusza.

Konieczność przeniesienia miejsca upadku samolotu.

Jak to wytłumaczą eksperci Zespołu?


Kontynuacja Zakręt Macierewicza cz. 2.