Peemka Peemka
18717
BLOG

Analiza "symulacji dynamicznej" prof. Pawła Artymowicza

Peemka Peemka Rozmaitości Obserwuj notkę 477

Podczas odbytej pod koniec maja w Kaziemierzu Dolnym konferencji zorganizowanej pod patronatem Polskich Linii Lotniczych LOT XV Konferencji Mechanika w Lotnictwie astrofizyk z Uniwersytetu w Toronto, specjalista od zagadnień formowania się pozasłonecznych planet, prof. Paweł Artymowicz przedstawił prezentację dotyczącą problemu półbeczki, jaką wykonać miał tupolew po utracie na pancernej brzozie części lewego skrzydła. Prezentacja Artymowicza bazowała na wyliczeniach dokonanych przy pomocy udostępnionego przez niego skryptu, pracującego w środowsku IDL. Była ona już krytykowana m.in. przez Marka Dąbrowskiego, który wykazał że prof. Artymowicz aby dopasować swoje wyniki do oficjalnej wersji przebiegu katastrofy poprzesuwać musiał geograficzne koordynaty szeregu drzew i zignorować dane z radiowysokościomierza. W tym tekście przeanalizuję niektóre aspekty programu Artymowicza.

Złożoność kodu

Udostępniony przez Artymowicza skrypt liczy 703 linijki kodu. Z tego około 150 to komentarze, około 45 to początkowe deklaracje zmiennych i stałych, większość reszty kodu odpowiedzialna jest za generowanie grafik i opisów. Numeryczne wartości do tych grafik i wykresów przeliczane są w petli o długości 28 linijek (nie licząc komentarzy). To w tej pętli, a ściślej w jej 4700 iteracjach przeliczane są podstawowe wartości parametrów, którymi Artymowicz podpiera swoje tezy. Zatem funcjonalne serce "zredukowanej symulacji dynamicznej" Artymowicza to 28 linijek, czyli około 4 procent kodu.

"Lokalna bezwładność brzozy"

Co zwraca uwagę w początkowych deklaracjach zmiennych to parametr oznaczony jako impact_dk. W sekcji komentarzy Artymowicz lakonicznie definiuje ten parametr jako "kąt o który obróciła się trajektoria tupolewa", zatem nie wiadomo o jaką trajektorię chodzi, ale w dalszych anglojęzycznych komentarzach jest już bardziej dokładny. Otóż parametr ten to "bounds from the local inertia of birch" czyli skokowa zmiana kursu spowodowana lokalną bezwładnością brzozy. Wartość tego parametru to 0.005 radiana czyli 0.29 stopnia. Skąd Artymowicz wziął tę „lokalną bezwładność brzozy” i dlaczego akurat 0.29 stopnia? Na poparcie owej „lokalnej bezwładności brzozy” nie podaje on nawet śladowych wyliczeń. Stała ta użyta jest w programie Artymowicza tylko raz, gdzie dodaje on 0.29 stopnia do początkowej wartości kursu samolotu. Pytanie po co, skoro nie zmienia to w żaden sposób trajektorii kursu samolotu. Poniżej grafika wygenerowana na podstawie programu Artymowicza w środowisku GDL (kliknij, by powiększyć).

 

 

Czarna krzywa to oryginalna zmiana kursu samolotu, którą znaleźć można w grafice nr 5 na 33 stronie prezentacji Artymowicza. Do początkowej wartości kursu dodane jest właśnie owe 0.29 stopnia wynikające z „bezwładności brzozy”. Czerwona krzywa to trajektoria bez dodanego początkowego parametru impact_dk. Jak widać, różnice są śladowe. Jako, że sam Artymowicz w konkluzjach swojej prezentacji (punkt 5, s. 39) stwierdza, że odtwarzał odchył kursowy samolotu z dokładnością do 1 stopnia zatem różnica 0.29 stopnia spokojnie mieści sie w tym marginiesie. Skąd więc ta „lokalna bezwładność brzozy”? Najprawdopodobniej to danina płacona przez Artymowicza komisji Millera. Według jej członków uderzenie lewym skrzydłem w pancerną brzozę spowodować miało natychmiastową zmianę kursu o 3.5 stopnia:

Zapoczątkowało to obrót samolotu w lewo względem osi podłużnej z jednoczesną zmianą kierunku lotu około 3,5°. Zmiana kierunku lotu została spowodowana reakcją samolotu na uderzenie w jego konstrukcję w odległości 10,8 m od osi pionowej samolotu. (Załącznik do raportu nr 4, Konfiguracja samolotu w chwili wypadku, s. 8)

Zatem mamy odpowiedź na pytanie skąd koncepcja „lokalnej bezwładności brzozy”. Problem w tym, że raport komisji Millera też nie podaje na jakiej podstawie stwierdzono lub wyliczono wartość „lokalnej bezwładności brzozy”. Co gorsza 3.5 stopnia z raportu Millera nie równa się 0.29 stopnia z programu Artymowicza. Zatem Artymowicz płaci ekspertom komisji Millera daninę, ale płaci im fałszywą monetą, dwunastokrotnie zmniejszając podaną przez nich wartość zmiany kursu w wyniku uderzenia w brzozę, tak aby pasowała do reszty jego założeń. Jak zatem wygladałaby trajektoria odchylenia kursu według Artymowicza gdybyśmy użyli oryginalnej wartości „lokalnej bezwładności brzozy” z załączników do raportu Millera, czyli 3.5 stopnia? Poniżej oryginalna krzywa bazująca na "lokalnej bezwładności brzozy" rzędu 2.9 stopnia z prezentacji Artymowicza strona 33 (kliknij, by powiększyć).

 

 

A poniżej jak to wygląda z wartością "lokalnej bezwładności brzozy" rzędu 3.5 stopnia z raportu PKBWL - wszystkie inne parametry bez zmian (kliknij, by powiększyć).

 

 

Jak widać, zmiana jest znacząca i kompletnie nie pasuje do trajektorii wyznaczonej za pomocą ścietych gałązek. Teraz można zrozumieć dlaczego Artymowicz musiał dwunastokrotnie zmniejszyć wartość skokowej zmiany kursu samolotu w wyniku uderzenia w brzozę podaną przez komisję Millera – oryginalna wartość z raportu poważnie zniekształcałaby jego wyliczenia. W swojej prezentacji taktownie to jednak przemilczał.

Zmiana kursu w wyniku "lokalnej bezwładności brzozy" o 3.5 stopnia znacząco wpływa także na odchylenie boczne trajektorii ( zob. rys. 7 z prezentacji Artymowicza, s. 25). Poniżej grafika wygenerowana przez program Artymowicza, ale z wartością "lokalnej bezwładności brzozy" zgodnej z raportem PKBWL, czyli 3.5 stopnia - wszystkie inne parametry bez zmian (kliknij, by powiększyć):

 

 

Trajektoria samolotu "precyzyjnie zmierzona na drzewach" rozmija się zarówno z alertem TAWS-a nr 38, jak i ostatnim zapisem FMS.

Kurs samolotu

Z problemem „lokalnej bezwładności brzozy” bezpośrednio związany jest inny – zmiany kursu samolotu. Kurs ten Artymowicz wylicza w modelu dwuwymiarowym, poprzez dzielenie wartości prędkości postępowej samolotu przez przyspieszenie boczne, które z kolei jest sinusem kąta odchylającej się od pionu siły nośnej. W konkluzjach swojej prezentacji Artymowicz stwierdził, że jego skrypt odtwarza trajektorię samolotu „dając bardzo dobrą zgodność z zapisami parametrycznymi” (s. 38). Twierdzenie to jest fałszywe. Omawiany powyżej rysunek nr 5 przedstawia zmianę kursu samolotu po zderzeniu z brzozą do momentu upadku. Jak to wygląda w zestawieniu z danymi parametrycznymi z czarnych skrzynek i z logów TAWS-a? Kiepsko. Zacznijmy od danych z QAR-ATM przedstawionych w załącznikach do raportu Millera (kliknij, by powiększyć):

 

 

Grafika pochodzi z załącznika do raportu Millera nr 2, rys. 17. Przebieg wybranych parametrów na ścieżce podejścia do lądowania – kanał poprzeczny – sterowanie lotkami, s. 20. Czerwona pionowa linia pochodzi z oryginalnej grafiki z załączników i opisana jest tam jako “Moment uderzenia w drzewo i oderwania części lewego skrzydła”. Na powyższym wykresie kurs samolotu od momentu uderzenia w brzozę przez blisko dwie sekundy jest równiusieńki jak spod linijki. Tymczasem skrypt Artymowicza przedstawia zupełnie inny obrazek: patrząc na grafikę już 1.7 sekundy po uderzeniu skrzydłem w brzozę odchylenie kursowe wynosi około 5 stopni. Można to sprawdzić dokładnie nie polegając tylko na wykresie, ale generując z programu Artymowicza dane numeryczne: dokładna wartość zmiany kursu w tym momencie (1.7 sekundy po uderzeniu w brzozę) to 5.45 stopnia (marker "t = 1.7s" na grafice Artymowicza jest minimalnie przesunięty w stosunku do faktycznego kroku czasowego).  Taką wartość zmiany kursu zdecydowanie wykluczają dane z QAR-ATM.

Jeszcze gorzej sprawa wygląda z logiem TAWS nr 38. Został on zapisany około 140 metrów po rzekomym uderzeniem skrzydłem w pancerną brzozę. Znając dosyć dokładnie prędkość postępową samolotu (75 m/s) znaczy to, że log ten został zapisany około 1.7-1.9 sekundy za brzozą. TAWS zapisuje zarówno kurs samolotu (geograficzny) oraz bezpośrednio prędkość kątową zmiany kursu. Parametry zapisane w logu TAWS nr 38 pokazują kurs geograficzny (true track) prawie -93 stopnie (dokładnie -92,988281). Konwencja zapisu jest nieco inna niż standardowa przy róży wiatrów, mianowicie dla kierunków NW oraz SW system loguje dane o wartości negatywnej przeciwnie do ruchu wskazówek zegara w stosunku do geograficznej północy. Czyli w tej konwencji kurs geograficzny -93 stopnie to 267 stopni w standardowym zapisie, czyli 260 według kursu magnetycznego, co widać na wykresach z rejestratorów (deklinacja magnetyczna dla obszaru Smoleńska to około 7 stopni). Kurs zapisany w tym logu doskonale koresponduje z magnetycznym kursem zapisanym w QAR-ATM. Jest on dokładnie taki sam, jaki miał samolot podczas stabilnej fazy lotu: 267 stopni. Żadnych pięciu stopni odchylenia Artymowicza.

TAWS zapisywał także bezpośrednio wartość prędkości odchylania się samolotu (track rate) od kursu, to znaczy jak szybko samolot skręca. Dla alertu nr 38 wynosi ona 0,064736 stopnia na sekundę. Setna część stopnia na sekundę. W samym środku rzekomej smoleńskiej półbeczki samolot utrzymuje dokładnie poprzedni kurs a szybkość skrętu jest znikoma. Jak to wygląda według skryptu Artymowicza? Nie zakodował on w swoim programie mechanizmu bezpośrednio przeliczającego prędkość kątową zmiany kursu, ale znając inne parametry łatwo taki mechanizm implementować. Na grafice poniżej biała krzywa to oryginalna zmiana kursu samolotu na osi odległości od brzozy; czerwona to prędkość kątowa tejże zmiany kursu (kliknij, by powiększyć).

 

 

W momencie zapisania sie logu TAWS nr 38, około 1.8 sekundy po hipotetycznym uderzeniu w brzozę prędkość kątowa wynosiła około 8 st/s. Znowu można to dokładnie sprawdzić korzystając z wygenerowanych wartości numerycznych: dla przedziału od 1.6 do 2 sekund po uderzeniu w brzozę prędkość kątowa według skryptu Artymowicza wzrastała od 7.46 do 8.64 st/s. Tymczasem wedle danych TAWS-a wartość ta wynosiła setne części stopnia na sekundę. I znowu dane parametryczne wykluczają wyniki generowane przez prościutki algorytm Artymowicza.

Przyspieszenie pionowe

Przyspieszenie pionowe w skrypcie (w układzie samolotu) nie jest wyliczane ze stanów układu samolotu, w rodzaju zmian kąta natarcia, prędkości, itd. Artymowicz zakodował je niezależnie, oświadczając, że pochodzi ono z zapisów parametrów lotu. Na stronie 21 swojej prezentacji (rys. 1) opisał on krzywą przyspieszenia pionowego jako "teoria aerodynamiczna zgodna z KBWL i MAK". To twierdzenie również jest fałszywe - ani teoria, ani zgodna z danymi parametrycznymi. Wyrażenie odpowiedzialne za przeliczanie przyspieszenia pionowego ma w skrypcie postać:

nnz = n0 + (n1-n0)*exp(-(tt/n_time)^7)

Gdzie:

n1 = 1.32 początkowe przyspieszenie pionowe

n0 = 0.36 końcowe przyspieszenie pionowe

n_time = 2.9 Artymowicz komentując tę wartość tajemniczo stwierdza, że jest to jakiś parametr obecny na rysunku 15 w załączniku do raportu Millera nr 2. Rysunek ten przedstawia podstawowe parametry lotu w ostatnich kilkudziesięciu sekundach lotu. Artymowicz nie wskazał o jaką wartość chodzi.

Podstawiając zdefiniowane wcześniej wartości liczbowe do wyrażnenia zarządzającego wartościami przyspieszenia pionowego otrzymamy:

nnz = 0.36 + 0.96*exp(-(tt/2.9)^7)

Wyrażenie tt zdefiniowane jest w skrypcie jako:

tt = (findgen(ii)+0.5)/ii*t2

Gdzie ii = 4700 oraz t2 = 4.7. Wyrażenie to zwraca w istocie prosty ciąg arytmetyczny, który wyrazić można jako:

Kluczowym elementem kontrolującym zmianę przyspieszenia pionowego w wyrażeniu nnz jest oczywiście funkcja wykładnicza. Jej wykres, odtworzony w programie MathCad jest następujący:

 

Gwarantuje ona, że druga stała w wyrażeniu nnz (0.96) najpierw utrzymuje swoją wartość przez co Artymowicz zagwarantował ponad dwusekundowe maksymalną wartość przyspieszenia pionowego, by następnie zacząć wykładniczo maleć, a tym samym cała funkcja nnz przeliczająca przyspieszenie pionowe wykładniczo maleje od wartości 1.32 do 0.36. Funkcja przyspieszenia pionowego odtworzona w pakiecie MathCad wygląda następująco:

Jak zatem widać, "teoria aeorodynamiczna" Artymowicza to raptem prosta funkcja malejąca wykładniczo w postaci:

Gdzie tt jest zdefiniowanym powyżej ciągiem artymetycznym.

No dobrze, powie ktoś. Nazywanie funkcji wykładniczej teorią aerodynamiczną to zwykłe napinanie się, ale co z drugą częścią twierdzenia Artymowicza, jakoby wartość przyspieszenia pionowego pasowała do danych MAK i PKBWL? Tu jest jeszcze gorzej.

Komentując parametr nnz Artymowicz stwierdza:

;    this model of vertical accleration fits Miller data (Zal.2) as well as
;     MAK report data, with noisy wave-like behavior omitted

Czyli opisane powyżej modelowanie przyspieszenia pionowego ma pasować do danych z raportu Millera i MAK tyle że z pominiętym "szumem" czyli falowymi zmianami zapisu. Jest tylko jeden problem. Zapis tego parametru od momentu uderzenia w brzozę to nic innego niż "wave-like behavior". Co więcej, to właśnie na podstawie skoków przyspieszenia PKBWL synchronizowała zapisy parametryczne i głosowe z kabiny pilotów, była to więc dla nich kluczowa charakterystyka. W przeciwieństwie zatem do PKBWL Artymowicz zdecydował się "pominąć" gwałtowne zmiany przeciążenia. Poniżej odpowiednie wykresy przyspieszenia pionowego z raportu MAK (górny) i PKBWL (dolny):

 

 

Jak widać, wartość przyspieszenia pionowego najpierw w okolicach pancernej brzozy wyraźnie spada, by powrócić do wartości około 1.3 g, następnie znowu gwałtownie spada, podnosi się by ostatecznie spaść do wartości 0.32 g. Wszystko to mało przypomina gładki spadek przyspieszenia pionowego wymuszony przez skrypt Artymowicza. Jak wyglądałoby modelowanie Artymowicza, gdyby choć trochę bardziej realistycznie odtworzyć widoczne w zapisach parametrycznych przyspieszenie pionowe? Funkcja cos((0.5*tt)^1.8)^2 nieco lepiej symuluje zmiany wartości przyspieszenia pionowego:

 

Poniżej krzywa przyspieszenia pionowego wygenerowana przez MathCad, gdzie zamiast funkcji wykładniczej Artymowicza zastosowana jest bardziej odpowiadająca realiom funkcja cosinusa (kliknij, by powiększyć):

 

Na korzyść modelu Artymowicza zignorowałem zupełnie pierwszy spadek przyspieszenia. Jak to zmieni jego wyliczenia? Zastępując w jego programie oryginalne wyrażenie nnz wyrażniem:

nnz = n0 + (n1-n0)*cos((0.5*tt)^1.8)^2

otrzymamy wyraźnie różne od oryginalnych trajektorie (kliknij, by powiększyć):

 

 

Końcówka lewego statecznika uderza w ziemię kilkadziesiąt metrów wcześniej w dodatku w tym samym miejscu co kikut lewego skrzydła. W tym samym momencie zapisuje się log TAWS nr 38. Dopiero chwilę później mają szanse zapisać się ostatnie dane FMS-a, które zresztą pokazują wysokość barometryczną około 15 metrów. Odchylenie boczne także jest wyraźnie mniejsze, gdzie końcówka statecznika nie miałaby pasującego do śladów odchylenia w kierunku południowym (kliknij, by powiększyć):

 

 

Oczywiście, dokładniejsze ujęcie zmian przyspieszenia pionowego wprowadziłoby jeszcze większe różnice. Jak z powyższego wynika parametr ten jest jednym z kluczowych elementów w modelowaniu hipotetycznej półbeczki wykonanej przez samolot. Modelowanie to jest bardzo wrażliwe nawet na jego niewielkie zmiany - to dlatego Artymowicz tak starannie dobrał przebieg przyspieszenia pionowego w swoim skrypcie.

Założony początkowy moment obrotowy

Innym kluczowym założeniem zakodowanym w skrypcie jest moment obrotowy, a ściślej początkowa wartość przyspieszenie kątowego (parametr dwdt). Jej wartość to 80 st/s^2. Dlaczego akurat taka? Dobre pytanie - Artymowicz nie przedstawił wyliczeń skąd się ta wartość wzięła. Pochodzić ona ma z innych wyliczeń dokonanych skryptem ykwroll.pro. I znowu pozostałe parametry są bardzo wrażliwe na zmiany tej wartości. Przykładowo, zmniejszenie początkowego przyspieszenia kątowego z 80 do 60 st/s^2 powoduje, że samolot przelatuje powyżej miejsca upadku (kliknij, by powiększyć):

 

 

Na moment pisania tego tekstu skrypt ykwroll.pro mający obliczać rozkład sił nośnych na skrzydłach tupolewa pozostaje nieudostępniony. "Nie chowamy niczego i nie zamiatamy pod dywan" - odgrażał się Artymowicz w swojej prezentacji (s. 10). Obietnicy nie dotrzymał.

Wnioski

Bardzo uproszczone, dwuwymiarowe modelowanie Artymowicza za główny cel wydaje się mieć nie próbę odtworzenia zachowania się samolotu po utracie części lewego skrzydła, ale możliwe maksymalne wpasowanie zachowania się samolotu do oficjalnej wersji. Pomimo tych wysiłków szereg generowanych przez jego program paramentrów wciąż nie pasuje do uzyskanych z czarnych skrzynek danych parametrycznych, jak dane kursowe czy przeciążenia pionowego, bazuje na wątpliwych - i starannie ukrytych - założeniach, jak początkowy moment obrotowy oraz początkowa prędkość kątowa. Artymowicz obficie korzysta z taktyki zwanej cherry picking, czyli wybieranie rodzynek: starannie dobiera on parametry początowe oraz przebieg kluczowych funkcji tak, aby jakoś wpasować swój model do oficjalnej narracji. Równie starannie ignoruje on dane, które nie pasują do tego modelu, jak skokowe zmiany przyspieszenia pionowego i dane kursowe oraz barometryczne zapisane w parametrach lotu i logu TAWS. Co można dedykować prof. Artymowiczowi:

Praktyka wybiórczego dobierania z konkurujących dowodów w taki sposób, aby pokreślić tylko te rezultaty, które wspierają określone stanowisko, przy jednoczesnym ignorowaniu lub odrzucaniu jakichkolwiek danych, które go nie wspierają znana jest jako cherry picking i stanowi znak firmowy kiepskiej nauki lub pseudonauki.
 
(Richard Somerville, Testimony before the U.S House of Representatives Committee on Energy and Commerce Subcommittee on Energy and Power, 8 March 2011.)

Właściwe podejście jest przeciwieństwem podejścia Artymowicza: mając szereg danych początkowych zapisanych w parametrycznych rejestratorach, jak zachowałby się samolot po hipotetycznej utracie części skrzydła? To pytanie jest jak najbardziej zasadne. Najwyraźniej jednak fizyka smoleńsla prof. Artymowicza nie dorosła jeszcze do tego etapu.

Peemka
O mnie Peemka

Nowości od blogera

Komentarze

Inne tematy w dziale Rozmaitości