Geometria tupolewa jest bardzo dobrze znana. Można też oszacować z dokładnościa 10-15% jego najważniejszy moment bezwładności.
Geometria tupolewa jest bardzo dobrze znana. Można też oszacować z dokładnościa 10-15% jego najważniejszy moment bezwładności.
you-know-who you-know-who
34856
BLOG

30. Aerodynamika beczki smoleńskiej najprościej, jak można

you-know-who you-know-who Katastrofa smoleńska Obserwuj temat Obserwuj notkę 1030

[Wersja z poprawionymi obliczeniami, oryginalny tekst zawierał pomyłki i literówki. Styczeń 2018.]

"Things should be as simple as possible, but not simpler" (A. Einstein o edukacji)


W rozdziale 28-ym bloga skomentowałem obliczenia różnicy sily nośnej przed i po urwaniu końcowki lewego skrzydla PLF 101 ktore zostaly zaprezentowane w tej animacji, i porównałem wyniki mojego programu aerodynamicznego liczacego siły i momenty sił medotą panelową VLM z oszacowaniami robionymi przez inne osoby. Tutaj pokażę analityczną, uproszczoną, metodę obliczeń aerodynamicznych, która daje wyniki zbliżone do pełniejszej teorii pół-numerycznej zaprezentowanej wcześniej. Nawet trójwymiarowa trajektoria przewidywana przez obliczenia analityczne nie jest najgorsza! Ten rozdział to kropka nad "i" w modelowaniu beczki. Pokazuje, ze nawet uproszczone teorie sa zgodne z zasadniczymi sladami na ziemi i w rejestratorach.

Pisząc artykuł do żurnala mechaniki teor. i stosowanej JTAM zauważyłem, że istnieje proste heurystyczne wyprowadzenie zasadniczych dwóch wielkości dawanych przez wspomniane wyliczenia robione uogólniona metodą panelową Prandtla-Weissingera. Daje to bezpośredni wgląd w to, skąd biorą się konkretne wartości przyspieszenia i tempa obrotu, bez potrzeby obliczeń numerycznych, jeśli naturalnie nie żadamy wysokiej dokładności. (Te rzeczy leżą u podstawy dowodu, że tupolew obrócił się niemal na plecy w czasie 4.7s, że przelecial tuż nad terenem ponad 300m po dokładnie tej trajektorii, o której mówią nam teraz pozostałe pod Smoleńskiem, połamane gałęzie i drzewa, że spadł dokładnie tam gdzie spadł, no i że nie było do tego potrzeba żadnych cudów, maskirowek ani zamachów dwubombowych. Dlatego są tak ważne.)  

Opisze analityczne oszacowania nie dlatego, ze metody pół-numeryczne przeze mnie stosowane (trzeba numerycznie rozwiazac  uklad rownan liniowych) są  aż tak skomplikowane, że niemożliwe do powtórzenia. W istocie, każdy student wyższych lat fizyki lub inżynierii lotniczej potrafi takie obliczenia przeprowadzić pisząc samodzielnie program w Matlabie czy innym języku. Wiem to stąd, ze jeden z nich, David, korzysta z mojej pomocy w studiach nad takimi rzeczami na zaliczenie w ramach kursu PSCD01, który jest końcowym projektem indywidualnym w zakresie nauk fizycznych (David chce być doktorantem w naszym znanym w świecie instytucie aeronautyki na UofT). Opiszę je, bo tak sie robi w nauce: porównanie metod analitycznych i pół-numerycznych wykazuje, że te drugie nie zawieraja żadnego grubego błędu, że można zaufać wynikom. Bez wyliczeń analitycznych, w ekstremalnym przykładzie czysto numerycznych obliczeń przy użyciu programu LD-DYNA, wiedzieliśmy do czego może dojść. Wyniki mogą być wręcz niefizyczne. Ale jak mowię, to przypadek ekstremalnie jednotorowego podejscia do zagadnienia.

W przypadku moich obliczeń macierzowych, droga naszkicowana jest na tyle dokładnie, ze wiele osób może ją przebyć ponownie i potwierdzić wyniki. Jednak nie wszystkim będzie się chciało wnikać w szczegóły moich obliczeń, ani też nie każdy usiądzie by rozwiązać 400 równań z 400 niewiadomymi, żeby policzyć (jak to tu opisałem) :

(i) że tutka zaraz po przycięciu skrzydła do 13m rozpiętości rozkręca przechył w tempie dω/dt ~ 83o/s2(bez uwzgl. lotki;  z uwzgl. lotki to niecałe 80o/s2; już to samo pokazuje ze pozostała prawa lotka nie miała najmniejszej szansy na zahamowanie beczki),

(ii) że tutka przy takim obciążeniu (n=1.32g) jakie miała przy brzozie Bodina nie może obracać się szybciej, niż około ω ~ 48o/s.

To są właśnie te zasadnicze dwie liczby, dla uzyskania ktorych liczyłem dokładnie aerodynamikę samolotu. Liczyłem całą siatke modeli, ale okazało się, że wszystkie je można opisać jednym czy dwoma wzorami analitycznymi z dokładnością paru procent; wzorami które zawierają dwie stałe podane powyżej (bądź bliskie im wartości, jeśli chcemy uwzględnić  szczegóły takie, jak odchylenie kierownicy wolantu w prawo).

Reszta moich obliczeń to osobna symulacja dynamiczna opisana w uprzednich rozdziałach, oparta częściowo na danych ze skrzynek (n(t), pochylenie(t), V(t)), a częściowo na rekonstrukcji dynamicznej przechyłu i odchyleń w pionie i poziomie od orginalnej trajektorii sprzed brzozy Bodina. Ale te obliczenia podałem juz w całości, także jako dobrze udokumentowaną treść programu w jezyku IDL, zresztą programu nieco prostszego niż ten liczący dwie opisane powyżej wielkości.

 

MOMENT BEZWŁADNOŚCI

Ponieważ do policzenia dω/dt potrzeba znajomości momentu bezwładności Ixx wokół osi podłużnej, pozwolę sobie naszkicować jak się wyznacza Ixx. Dzielimy myślowo samolot na jak najwięcej części składowych: kadłub, bagaże, pasażerowie, paliwo w osobnych kilku zbiornikach, skrzydła, stateczniki i tak dalej. Następnie ze znanej geometrii samolotu dość dokładnie możemy ocenić średnio-kwadratową odległość od osi tych części składowych. Jest to jakby promień równoważnego walca mającego taki sam moment bezwładności, jak dana część.

Czasem można użyć do wyznaczenia tego promienia prostego całkowania, np. trzeba to zrobić w przypadku rozciągłych obiektów - skrzydeł. Innym razem ta odległość jest od razu oczywista, jak np. w przypadku silników. W takim razie potrzeba już "tylko" wyznaczyć masy poszczególnych składowych. I to nie jest łatwe do wyszukania na sieci! Dokładne plany konstrukcyjne i dokładne masy to szczegóły techniczne znane producentowi, ale na ogol strzeżone przed okiem postronnych. Jednak te masy trzeba znać, by projektowac sensownie nowe samoloty i nie zaczynać zawsze od Adama i Ewy. Dlatego w pracach z bibliografii bloga takich jak: Sadraey (2009) i Kroo (2006), oraz Cleveland (1970), rozwijane sa prawa skalowania. Pomagają one zrozumieć w oparciu o wytrzymałość materiałów i zasady budowy skrzydła kesonowo - dźwigarowego, dlaczego masa skrzydeł skaluje się tak a nie inaczej wraz z całkowitą (maksymalną) masą startową samolotu transportowego, a po tym jak sie wyprowadzi odpowiednie zależności (Kroo 2006), jak się ta teoria ma do rzeczywistości - w oparciu o znane dawne konstrukcje samolotów. I te zgadzają się z prawami skalowania, jeśli chodzi o skrzydła. Dlatego mówię tu tyle o skrzydłach, ponieważ nawet jeśli pomylimy się trochę w ocenie masy kadłuba, nie będzie to bardzo ważne dla momentu bezwładności wokół osi podłużnej (kwadrat odległości jest stosunkowo mały), najważniejsze jest natomiast ocenić dobrze mase skrzydła, a także masę zespołu podwozia, bo ciężkie wózki koła głównego są oddalone aż o 6m od osi.

Wstawiłem do wzoru Kroo wartości odpowiadające TU-154M i otrzymałem oszacowanie z prawa skalowania równe 11.4 t dla masy strukturalnej skrzydła o kształcie skrzydła tupolewa, zaprojektowanego na wytrzymywanie przeciążeń do 2.5g*1.5 (gdzie 1.5 to lotniczy margines bezpieczeństwa). Po dodaniu masy obiektów siedzących w skrzydle ale nie nadających mu żadnej wytrzymałości, co stanowi typowo 25% dodatkowej masy, otrzymałem masę nieuszkodzonego skrzydła (tj. obu skrzydeł) równą 14.2 t. Masa strukturalna 11.4 t zgadza sie wspaniale z całym szeregiem mas w historycznych konstrukcjach takich jak: B 707, DC-8 DC-10, B 737, B 727 i in. Wszędzie tam, nota bene, jako reguły kciuka moża użyć zamiast skomplikowanego wzoru Kroo po prostu 11% TOW (max take-off weight) i będzie ok. Nota bene, dr Berczyński z zespołem twierdzi, ze skrzydło tupolewa było niesłychanie masywne, co jest wyssane z palca. Masa uzasadniona tutaj zgadza się z moimi wcześniej opublikowanymi w blogu ocenami pola przekrojów poprzecznych aluminum D16 w skrzydłach.

Masy i równoważne promienie sa w tupolewie następujące:

utracona końcowka skrzydła: 0.7t w odległości 16m,

kadłub wraz pasażerami i cargo: 40t w odl. 1.5m,

paliwo: 6t w odl. 2m oraz 5t w odl. 3.7m (dwa rodzaje zbiorników),

silniki: 6t w odl. 3.7m,

nieuszkodzone skrzydła: 14t w odl. 8m,

stabilizatory: 4t w odl. 5m,

podwozie wraz konsolami: 4t w odl. 6m.

________

Masa całkowita 78.6t przez zderzeniem z brzoza, 77.9 po nim.

Moment bezwładności Ixx = 1405 t*m2 przed utratą, oraz

 Ixx =1223 t*m2

po utracie końcówki skrzydła. (ta wartość zawiera w sobie malutka poprawkę na te okoliczność, ze oś bezwładności przesuwa sie o 0.125m w prawo po utracie końcówki skrzydła). Wynik ten odpowiada masie tupolewa (po urwaniu części skrzydła równej M=78 ton=78 tys. kg) rozłożonej na walcu o promieniu prawie rx=4 m, tak że

Ixx = M rx2

 

 

WYPROWADZENIE MOMENTU SIŁY NA POCZĄTKU BECZKI

Teraz możemy już przejść do zasadniczego oszacowania heurystycznego momentu siły T (ang.: torque) i  ściśle z nim związanego dω/dt, oraz maksymalnej pr. kątowej ω w czasie beczki.  To będzie proste oszacowanie, ale nie prostsze niż konieczne (Einstein mówił: trzeba wszystko wyjaśniać tak prosto jak sie da, ale tylko tak i nie prościej.) Przyjmiemy do oszacowań zasadniczy wynik teorii skrzydła o skończonej długości Ludwika Prandtla: siła jest tym "lepiej" rozłożona wzdłuż rozpiętości skrzydła, im bardziej jej wykres zbliżony jest do elipsy. W jakim sensie "lepiej"? W sensie minimalnego oporu indukowanego (zob. polskie wiki pt . Opór aerodynamiczny a lepiej te angielska stronę).

Skrzydła tupolewa z ich dużym skosem zaprojektowane zostały tak, że rozkład siły (wyliczony np. przez mój program aerodynamiczny) nieobracającego się skrzydła jest w istocie dość podobny do tej idealnej ze względu na opór indukowany funkcji

dF/dy =  f0 (1 - y2/S2)1/2

gdzie  f0=const (oczywiście ta stała jest proporcjonalna do ciśnienia dynamicznego ρ V2/2, ale ta informacja nie będzie nam tu potrzebna, wyliczymy fze znanego przeciążenia). To podobieństwo nie jest oczywiste, bo kształt skrzydła tupolewa nie jest wcale zbliżony do eliptycznego. 

S = długość czy dokładniej pół-rozpiętość skrzydła. W naszym problemie rozpiętosci lewego i prawego skrzydła równają się

L = 13 m

R = 18.77 m

i takie wartości  stałej S należy przyjąć  na odpowiednich skrzydłach. Na obu skrzydłach przyjmijmy wiec rozkład eliptyczny siły o tej samej wartości siły maksymalnej, ale o rożnej rozpiętości. W naszym oszacowaniu dla uproszczenia zaniedbujemy wpływ wąskiego kadłuba, gdyż wyniki są znacznie bardziej czułe na to, co dzieje się na końcach skrzydeł, niż przy kadłubie.

Całkując iloczyn ramienia działania siły y oraz siły dF(y), dostajemy wkład dT do momentu siły T0 (indeks 0 wskazuje, że mowa o zerowej szybkości obrotu skrzydła). Całkowanie zaczynamy na kikucie lewego skrzydła przy y=-L, a kończymy na końcu prawego przy y=+R.  Zatem

T0 = ∫-LR (dF/dy) y dy.

Używając powyższego wzoru na dF/dy, osobno dla skrzydła lewego i prawego, dostajemy

T0 = (f0/2) (R2-L2)  ∫01  (1 - x)1/2 dx = (f0/3) (R2-L2)

Potrzebną nam wartość f0 wyznaczamy ze znanej siły nośnej F, równej F = n M, gdzie M to masa samolotu, zaś n jest przyspieszeniem pionowym środka masy (1.32 g ~13 m/s2):

F = ∫-LR (dF/dy) dy = f0(R+L)  ∫01  (1 - x2)1/2 dx = (π f0/4) (R+L) = n M,

f0 = 4nM/[π(L+R)].

Wstawiając f0 do wyrażenia na moment siły, mamy

 T0 =  (4/3π)(R-L) n M

czyli moment siły to siła  F=n M razy ramię działania równe (4/3π) (R-L).

Pozostaje nam zapisać przyspieszenie kątowe krótką chwilę po urwaniu skrzydła jako

 dω/dt =  T0/ Ixx = (4/3π) (R-L) n M/Ixx= (4/3π) (R-L) n/rx2

Wstawiamy wartości opisujące tupolewa (n = 1.32 g = 12.95 m/s2, rx = 4 m, R = 18.77 m, L = 13 m)  i dostajemy oszacowanie

dω/dt ~  114o/s2.

To oszacowanie jest o ok.jedną trzecią większe niż dokładniejsza wartość dω/dt ~ 83o/s2, wyliczona dla skrzydła z opuszczonymi jak w tupolewie klapami i niewychyloną lotką, metodą macierzową, z uogólnionej teorii Prandtla-Weissingera.

 

KOŃCOWA (MAKSYMALNA) PRĘDKOŚĆ KĄTOWA PRZECHYŁU

Teraz chcemy uwzględnić wpływ obrotu wokół osi podłużnej. Ponieważ siła nośna zależy od lokalnego kąta natarcia strugi powietrza na linie średniego profilu skrzydła, zapiszmy ten kąt (ang.: AOA) jako α:

 α = αkadl - αwzn + αskr - αL0 - αobr

(kadl = kąt kadłuba czyli pochylenie do poziomu; wzn = kąt wznoszenia ; skr = kat skręcenia (washout) skrzydła; L0 = kąt zerowej siły nośnej, na przykład w tupolewie z klapami opuszczonymi 36o wynosi on średnio dla całego skrzydła tuz nad ziemia, αL0= -6.6o ;  obr = kąt wywołany obrotem skrzydła robiącego beczkę). Zauważmy, ze α to taki kąt natarcia, który znika gdy siła nośna znika; nie jest to kąt natarcia w stosunku do jakichś innych, bardziej typowych kierunków, jak np. oś samolotu. Nie jest to więc kat AOA (angle of attack) taki, jakim się go zazwyczaj widuje w podręcznikach i instrukcjach obsługi samolotu TU-154. Przy AOA=0 samolot produkuje bowiem dodatnią siłę nośną, a my musimy ustalić kąt do którego siła nośna będzie, w przybliżeniu małych kątów, wprost proporcjonalna.

Kąt αobr to lokalna prędkość pionowa danego fragmentu skrzydła podzielona przez prędkość postępową V = 75 m/s), co daje ten mały kąt w radianach. Pomijając obrót i przyjmując średni kat skręcenia skrzydła równy 1o, co jest średnia arytmetyczna z ekstremalnych kątów skręcenia +3 i -1 stopni na początku i końcu prawego skrzydła, przy brzozie tupolew miał

 α0 =  αkadl - αwzn + αskr - αL0 

 α0 ~ 13.5o - 4.8o   +1o   +6.6o  ~   16.3o

Uwzględniając obrót wokół osi podłużnej z prędkością kątową ω,  mamy

 α = α0 - ω y/V  = α0 (1 - ω y/Vα0).

Zatem T = ∫-LR f0 (1 - y2/S2)1/2  y (1 - ω y/Vα0) dy

gdzie S=L  dla  y<0,  zaś  S=R  dla  y>0.

W wyrażeniu podcałkowym (1-ωy/Vα0), 1 produkuje w wyniku końcowym T0, natomiast (-ωy/Vα0) produkuje moment siły zależny od rotacji, nazwijmy go  Trot

  Trot = -f0 ω /(Vα0)  ∫-LR  (1 - y2/S2)1/2 y2 dy.

Można w całce zamienić zmienne z  y na t, gdzie t=(y/S)2, a wtedy całka przybiera postać typu (1/2) ∫0 [t(1 - t)]1/2 dt =  π/16.  Asymptotyczna prędkość kątowa obrotu to takie ωmax, dla którego  T = 0.

T =  T0 [1  - (3π/16) (R3+L3)/(R2-L2) (ωmax/ Vα0)] = 0, wiec

ωmax = (16/3π) (R2-L2)/(R3+L3) Vα0

Wstawiając teraz L=13 m i R=18.77 m, V=75 m/s, oraz  α0 ~16.3o, dostajemy następujące oszacowanie asymptotycznego tempa obrotu:

ωmax  ~ 43.2o/s

co jest dość dokładnym oszacowaniem tej wielkości, równym 90% dokładniejszej wartości  ω ~ 48o/s, otrzymanej z uwzględnieniem klap, lotki, steru kierunku, stopniowego skręcenia kształtu skrzydła, skosu i zmiennej szerokości skrzydeł oraz wzajemnego oddziaływania różnych części skrzydeł na lokalne kąty ataku strug w wyniku produkcji wirów, a także wpływ kadłuba  (uogólniona teoria Prandtla-Weissingera).

 

WYNIKI ALTERNATYWNEJ SYMULACJI 

Jeśli użyjemy wyprowadzonych w przybliżeniu wartości 43o/s zamiast 48o/s oraz 114o/s2 zamiast 80o/s2, dostaniemy następującą trajektorię: rys. 1-4, a tutaj rys. 4-8, która nie pasuje do danych smoleńskich aż tak dokładnie, jak wyliczona poprzednio, ale można zauważyć, że także pasuje "nieźle". Jest tak zapewne dlatego, że o ile przyspieszenie obrotu wzrosło w mniej dokładnym obliczeniu, to jednak końcowy kąt obrotu przy uderzeniu w ziemię nie zwiększył się tylko zmalał do 140 stopni przechyłu, gdyż asymptotyczne szybkości obrotu w mniej dokładnym oszacowaniu są zauważalnie mniejsze, niż w teorii pełniejszej. Błędy się częściowo skasowały, a trajektoria --nawet ta niedokładna-- ma sens: kikut lewego skrzydła nadal leci nad a nie pod ziemią (4-9m zamiast poprzednich 5-10m nad terenem), samolot schodzi ok. 20 stopni z kursu i tak dalej. Niewielkie zmiany omawianych wyników aerodynamicznych nie są więc ważne dla scenariusza katastrofy, który został tu po raz kolejny potwierdzony. 

 

SUMMA SUMMARUM  

Uzyskaliśmy proste oszacowania dwóch najważniejszych wielkości cechujących samolot TU-154M po utracie 1/3 lewego skrzydła: przyspieszenia kątowego i maksymalnej prędkości kątowej przy znanym przeciążeniu. Te oszacowania sa zgodne z dużo uważniej uprzednio wyliczonymi wartościami,  z dokładnością, odpowiednio, ok. 37% i 10%. Podobne wyliczenia wielu z was, np. andrzejmat, zrobiło już dawno. Mi było po prostu łatwiej pokazać na 3-wymiarowej trajektorii co z tych wartości wynika. 

 

Na koniec dwa postscriptum:

1.  O ile przyspieszenie kątowe tupolewa miało rzeczywiście początkowo przybliżoną wartość 80 stopni na sekundę kwadrat, to maksymalna prędkosc kątowa 48o/s nie została jednak przez samolot osiągnięta, ponieważ po czasie potrzebnym do zbliżenia się do niej (~2 s) zmieniły się warunki lotu: wzrosło a zaraz potem spadło przyspieszenie n(t) i zmieniły się wychylenia sterów. To wszystko uwzględniono w symulacji dynamicznej Artymowicza (2014, JTAM, in prep.) posługującej się wyliczonymi tam parametrami i prawami skalowania.

2. Nie wspominałem jeszcze o dwuczęściowym wywiadzie  prasowym dla polonijnego dziennika kanadyjskiego: dziennik Dziennik. Zachęcam do kliknięcia na marginesie po prawej stronie w dziale 'Polecane strony'.

Dodatkowo, dla części agnostycznej i antyreligijnej Polaków, ktorzy czytają czasopismo Fakty i Mity, tez udzielilem tam wywiadu ktory ukazał się wczoraj i został zalinkowany na "Polecanych stronach".

Informacji udzielam (bez wynagrodzenia!) wszystkim, od FiM do GP i ND. Powiedziałem to dziennikarzom tych ostatnich gazet, którzy jako jedyni (prasy tzw. lemingowej nie było) znaleźli sie w Kazimierzu nad Wisłą w końcu maja br. Jeszcze do tej pory nie skontaktowali się ze mną. Ja to rozumiem i czekam cierpliwie, nie obiecywali mi nic przecież. Mają już moje zdjęcia stamtąd, robiliśmy w pięknej kamienicy, gdzie odbywała się konferencja.

 

Zobacz galerię zdjęć:

Symulacja dynamiczna oparta o uproszczone wyprowadzenie parametrow - rysunki 5-8. Link do obrazu w lepszej rozdz. w tekście.
Symulacja dynamiczna oparta o uproszczone wyprowadzenie parametrow - rysunki 5-8. Link do obrazu w lepszej rozdz. w tekście. Alternatywna symulacja beczki z parametrami uzyskanymi w sposob uproszczony. Link do lepszej rozdzielcz. obrazu w tekście. Wszystkie beczki opisane sa takimi samymi prawami aerodynamiki. Beczka smoleńska była szybką beczką lotkową (nie autorotacyjną).

Nazywam się Paweł Artymowicz, ale wolę tu występować jako YKW. Moje wyniki zatwierdził w 2018 r. i podał za wzór W. Biniendzie jako wiarygodne wódz J. Kaczyński (naprawdę! oto link). Latam wzdłuż i wszerz kontynentu amerykańskiego (link do mapki), w 2019 r. 40 godz. za sterami, ok. 10 tys. km; Jestem niezłym (link), szeroko cytowanym profesorem fizyki i astrofizyki [link] (zestawienie ze znanymi osobami poniżej). Kilka krajów nadało mi najwyższe stopnie naukowe. Ale cóż, że byłem stypendystą Hubble'a (prestiżowa pozycja fundowana przez NASA) jeśli nie umiałbym nic policzyć i rozwikłać części "zagadki smoleńskiej". To co mówię i liczę wybroni się samo. Nie mieszam się do polityki, ale gdy polityka zaczyna gwałcić fizykę, a na dodatek moje ulubione hobby - latanie, to bronię tych drugich, obnażając różne obrażające je teorie z zakresu "fizyki smoleńskiej". Zwracam się do was per "drogi nicku" lub per pan/pani jeśli się podpisujecie nazwiskiem. Zapraszam do obejrzenia wywiadów i felietonów w artykule biograficznym wiki. Uzupełnienie o wskaźnikach naukowych w 2014 (za Google Scholar): Mam wysoki indeks Hirscha h=30, i10=41, oraz ponad 4 razy więcej cytowań na pracę niż średnia w mojej dziedzinie - fizyce. Moja liczba cytowań to ponad 4100 [obecnie 7500+, h=35]. Dla porównania, prof. Binienda miał wtedy dużo niższy wskaźnik h=14,  900 cytowań oraz 1.2 razy średnią liczbę cytowań na pracę w dziedzinie inżynierii. Inni zamachiści (Nowaczyk, Berczyński, Szuladzinski, Rońda i in. 'profesorowie') są kompletnie nieznaczący w nauce/inż. Częściowe  archiwum: http://fizyka-smolenska.blogspot.com. Prowadziłem też blog http://pawelartymowicz.natemat.pl. 

Nowości od blogera

Komentarze

Inne tematy w dziale Polityka