Spotlight #71: Symulacje z fizyki, chemii i inżynierii

Każda symulacja na tej platformie zaczyna się od konkretnego pytania: jak właściwie wygląda to równanie, gdy je uruchomimy? Niedawno wzbogacone symulacje omówione w tym przeglądzie łączy właśnie ta cecha. Obejmują one zjawiska fizyczne pojawiające się w programach studiów licencjackich i magisterskich z fizyki, chemii i inżynierii, które jednak rzadko są prezentowane interaktywnie — zwykle dlatego, że koszt obliczeniowy był historycznie zbyt wysoki dla przeglądarki. WebGL i zoptymalizowane metody numeryczne to zmieniły. Każda z sześciu opisanych tu symulacji działa z interaktywną prędkością na dowolnym nowoczesnym urządzeniu, bez żadnych wtyczek.

I. Dyfrakcja Bragga

Krystalografia rentgenowska to technika, która ujawniła strukturę DNA, kształty centrów aktywnych białek oraz rozmieszczenie atomów w każdym materiale, którego właściwości rozumiemy dogłębnie. Technika ta opiera się na formule tak eleganckiej, że wyprowadził ją William Lawrence Bragg w wieku 22 lat:

nλ = 2d sinθ

gdzie:
  n      : rząd odbicia (dodatnia liczba całkowita)
  λ      : długość fali promieniowania rentgenowskiego
  d      : odległość między płaszczyznami sieci krystalicznej (odległość międzypłaszczyznowa)
  θ      : kąt padania (poślizgowy)

Dla FCC (RSC): odbicia dozwolone, gdy h, k, l są wszystkie nieparzyste lub wszystkie parzyste
Dla BCC (RPC): odbicia dozwolone, gdy h + k + l jest parzyste
Sieć prosta regularna: wszystkie hkl dozwolone

Symulacja oblicza proszkowy wzór dyfrakcyjny dla struktur krystalicznych regularnych: ściennie centrowanej (FCC), przestrzennie centrowanej (BCC) oraz prostej regularnej, sumując wkłady wszystkich dozwolonych odbić Bragga aż do maksymalnego kąta 2θ zdefiniowanego przez użytkownika. Natężenia pików uwzględniają czynnik struktury, czynnik Lorentza-polaryzacji oraz gaussowskie poszerzenie modelujące skończony rozmiar krystalitów (równanie Scherrera: β = Kλ / L cosθ). Wynik jest wyświetlany jako diagram kreskowy nałożony na symulację pierścieni Debye'a-Scherrera.

Najbardziej pouczającą funkcją jest przełącznik struktury krystalicznej. Wystarczy przełączyć FCC na BCC, a systematyczne wygaszenia natychmiast się zmieniają: struktura FCC tłumi odbicia o mieszanej parzystości wskaźników Millera (dzięki czemu odbicia (1,0,0) oraz (1,1,0) znikają), podczas gdy BCC tłumi odbicia o nieparzystej sumie h+k+l. Dzięki temu związek między geometrią komórki elementarnej a wzorem dyfrakcyjnym staje się bezpośrednio dostrzegalny, a nie wyłącznie algebraiczny. Suwaki pozwalają regulować długość fali λ oraz parametr sieci a, umożliwiając użytkownikom odtworzenie opublikowanych wzorów dyfrakcyjnych dla rzeczywistych materiałów poprzez dopasowanie położenia pików.

II. Długość Debye'a

W plazmie swobodne ładunki przegrupowują się, aby ekranować dowolne zaburzenie elektrostatyczne na charakterystycznej odległości zwanej długością Debye'a. Wzór ma postać:

λD = √(ε0 kB T / n q²)

Ekranowany potencjał wokół ładunku próbnego:
  φ(r) = (q / 4πε0 r) exp(−r / λD)

To potencjał Yukawy: prawo Coulomba pomnożone przez
zanik wykładniczy w skali λD.

Symulacja umieszcza ładunek próbny w centrum dwuwymiarowej siatki plazmy i renderuje wynikowy potencjał elektrostatyczny jako mapę cieplną w pseudokolorach, ze skalą barwną obejmującą zakres od silnego pola centralnego po ekranowany obszar peryferyjny. Suwaki kontrolują temperaturę elektronową T oraz gęstość liczbową n. W miarę ich przesuwania mapa cieplna przerysowuje się w czasie rzeczywistym, pokazując, jak λD skaluje się jak √(T/n): wzrost temperatury zwiększa długość Debye'a (gorętsze elektrony mogą oddalić się bardziej od obszaru ekranowania), natomiast wzrost gęstości ją zmniejsza (więcej elektronów jest dostępnych, aby ekranować zaburzenie na krótszym dystansie).

Towarzyszący wykres radialny w skali logarytmicznej potwierdza zanik wykładniczy i pokazuje przejście od nieekranowanego prawa Coulomba 1/r przy krótkich odległościach do wykładniczo stłumionego ogona Yukawy przy dużych r. To podwójne przedstawienie — mapa cieplna i profil radialny jednocześnie — jest funkcją najczęściej wskazywaną przez wykładowców fizyki plazmy korzystających z symulacji.

III. Cykl Carnota

Cykl Carnota to teoretyczna górna granica sprawności silnika cieplnego. Żaden silnik pracujący między zbiornikiem gorącym o temperaturze TH a zbiornikiem zimnym o temperaturze TC nie może przekroczyć sprawności Carnota:

ηCarnota = 1 − TC / TH

Cztery etapy cyklu na diagramie P-V:
  1. Rozprężanie izotermiczne przy TH:  pochłaniane QH, gaz wykonuje pracę
  2. Rozprężanie adiabatyczne:           Q = 0, temperatura spada do TC
  3. Sprężanie izotermiczne przy TC: oddawane QC, praca wykonywana nad gazem
  4. Sprężanie adiabatyczne:         Q = 0, temperatura rośnie do TH

Praca netto = pole powierzchni zamknięte przez cykl na diagramie P-V
Sprawność η = W_netto / QH

Symulacja animuje gaz roboczy przechodzący przez wszystkie cztery etapy na żywym diagramie P-V, a punkt stanu zakreśla zamknięty cykl w trakcie trwania animacji. Suwaki kontrolują TH i TC, a towarzyszący diagram T-S przedstawia cykl w przestrzeni entropii (gdzie cykl Carnota jest prostokątem, co czyni geometryczną interpretację sprawności natychmiast zrozumiałą). Bieżące zestawienie pokazuje aktualne wartości QH, QC, Wnetto oraz η, aktualizowane z każdym cyklem. Panel entropii to funkcja, która czyni tę symulację naprawdę bardziej użyteczną niż statyczny diagram z podręcznika: studenci mogą zobaczyć, dlaczego odwracalność jest termodynamicznie wyjątkowa.

IV. Reakcja Biełousowa-Żabotyńskiego

Reakcja Biełousowa-Żabotyńskiego (BŻ) to rodzina oscylujących reakcji chemicznych, które spontanicznie generują wzory przestrzenne: koncentryczne pierścienie, obracające się spirale i biegnące fale w czymś, co wygląda na dobrze wymieszany roztwór. Mechanizm obejmuje reakcję redoks oscylującą między stanem utleniającym a redukującym, sprzężoną z autokatalitycznym sprzężeniem zwrotnym. Model Oregonatora sprowadza mechanizm do trzech substancji:

Równania różniczkowe Oregonatora (uproszczone):
  ∂u/∂t = Du∇²u + (1/ε)(u − u² − f v (u − q)/(u + q))
  ∂v/∂t = Dv∇²v + u − v

gdzie:
  u : stężenie HBrO2 (autokatalizator)
  v : stężenie Ce(IV) (utleniony katalizator)
  ε : stosunek skal czasowych (~0,04)
  f : parametr stechiometryczny (~1,4)
  q : stosunek stałych szybkości (~0,002)

Symulacja całkuje równania Oregonatora na dwuwymiarowej siatce za pomocą jawnego schematu różnic skończonych. Startując z losowych warunków początkowych, układ po krótkim etapie przejściowym spontanicznie rozwija spiralne wzory falowe — dokładnie tak, jak obserwuje się to w reakcji BŻ w szalce Petriego. Pole barwne odwzorowuje stężenie HBrO₂ na klasyczną niebiesko-czerwoną oscylację reakcji BŻ katalizowanej cerem.

Najbardziej uderzającą cechą jest interakcja fal spiralnych. Dwie przeciwbieżnie wirujące spirale, które się zbliżają, anihilują, natomiast spirala zderzająca się z granicą odbija się. Te reguły zderzeń wyłaniają się w sposób czysto emergentny z lokalnej dynamiki reakcyjno-dyfuzyjnej, bez żadnej jawnej reguły określającej, co dzieje się na granicach lub gdy fale się spotykają. Ta emergencja globalnego porządku z czysto lokalnych interakcji to własność, która czyni reakcję BŻ kanonicznym przykładem samoorganizacji w chemii dalekiej od równowagi.

V. Zasada Bernoulliego

Równanie Bernoulliego to jedno z najczęściej stosowanych i najczęściej błędnie rozumianych równań fizyki wprowadzającej. Nie jest ono samodzielnym prawem aerodynamiki; jest to stwierdzenie o zachowaniu dla przepływu nielepkiego, nieściśliwego i ustalonego wzdłuż linii prądu:

P + ½ρv² + ρgh = stała wzdłuż linii prądu

Dla przepływu poziomego (h = const):
  P1 + ½ρv1² = P2 + ½ρv2²

Równanie ciągłości (przepływ nieściśliwy):
  A1 v1 = A2 v2

Zatem węższy przekrój oznacza wyższą prędkość, co oznacza niższe ciśnienie.

Symulacja modeluje przepływ przez kanał, którego profil przekroju poprzecznego może narysować użytkownik. Pole prędkości jest rozwiązywane metodą funkcji strumienia na dwuwymiarowej siatce, a linie prądu są renderowane jako animowane ścieżki znaczników. Pole ciśnienia jest obliczane na podstawie prędkości za pomocą równania Bernoulliego i wyświetlane jako nakładka barwna, dzięki czemu studenci od razu widzą spadek ciśnienia w obszarze zwężonym. Tryb przepływomierza Venturiego dodaje wskaźniki poboru ciśnienia w dwóch punktach oraz wyświetla różnicę ciśnień i wywnioskowane natężenie przepływu — dokładnie tak, jak działałby fizyczny przepływomierz Venturiego.

Symulacja pokazuje również ograniczenia równania Bernoulliego: wystarczy przesunąć suwak liczby Reynoldsa wystarczająco wysoko, a przepływ oderwie się od ścianek kanału, tworząc strefy recyrkulacji, w których równanie Bernoulliego przestaje obowiązywać. To uczciwe potraktowanie założeń równania — ustalony, nielepki, nieściśliwy — czyni symulację bardziej użyteczną do nauczania krytycznego myślenia niż wersja, która po prostu potwierdza formułę.

VI. Ruch Browna

Ruch Browna był tematem jednej z czterech prac Einsteina z 1905 roku i dostarczył pierwszego ilościowego dowodu na słuszność hipotezy atomowej. Ziarnko pyłku obserwowane pod mikroskopem zakreśla chaotyczną, pozornie losową ścieżkę, ponieważ jest bombardowane przez pojedyncze cząsteczki wody. Einstein wykazał, że średnie kwadratowe przemieszczenie rośnie liniowo z czasem:

⟨r²(t)⟩ = 2dDt

gdzie:
  d : wymiar przestrzenny (2 dla ruchu dwuwymiarowego)
  D : współczynnik dyfuzji = kBT / (6πηR)
  η : lepkość płynu
  R : promień cząstki
  kB : stała Boltzmanna
  T : temperatura

Ścieżka jest procesem Wienera: przyrosty Δr ~ N(0, 2DΔt)

Symulacja przedstawia wiele cząstek Browna jednocześnie, kreśląc ich ścieżki jako zanikające ślady na dwuwymiarowym płótnie. Bieżący wykres średniego kwadratowego przemieszczenia w funkcji czasu rośnie liniowo, potwierdzając zależność Einsteina. Współczynnik dyfuzji jest wyznaczany poprzez dopasowanie nachylenia do krzywej MSD i porównywany z teoretyczną wartością Stokesa-Einsteina dla wybranego rozmiaru cząstki i temperatury. Suwaki kontrolują temperaturę, lepkość i promień cząstki, a współczynnik dyfuzji aktualizuje się na bieżąco, dzięki czemu zależności od 1/R i T są bezpośrednio obserwowalne.

Opcjonalny tryb pojedynczej cząstki przybliża jedną trajektorię i wyświetla rozkład przemieszczeń w wybranych przez użytkownika odstępach czasu, potwierdzając, że przemieszczenia mają rozkład Gaussa o wariancji rosnącej liniowo — definiującą statystyczną sygnaturę procesu Wienera i matematyczny fundament całego współczesnego rachunku stochastycznego.

Wspólne motywy

Patrząc na te sześć symulacji łącznie, warto jednoznacznie wskazać kilka powtarzających się zasad projektowych. Po pierwsze, każda symulacja wyświetla co najmniej dwie reprezentacje jednocześnie: wizualizację przestrzenną oraz wykres ilościowy. To podejście podwójnej reprezentacji oznacza, że użytkownicy mogą przechodzić między intuicją (wzorem) a pomiarem (krzywą), nie przełączając się między różnymi narzędziami czy stronami.

Po drugie, każda symulacja ma co najmniej jeden suwak lub interaktywny parametr, który zmienia zjawisko jakościowo, a nie tylko ilościowo. Przełączenie między FCC a BCC w symulacji Bragga nie tylko przesuwa piki — zmienia, które piki w ogóle istnieją. Podniesienie temperatury w symulacji ruchu Browna nie tylko przyspiesza cząstki — zmienia współczynnik dyfuzji w konkretny, mierzalny sposób. To właśnie w tych jakościowych przejściach zachodzi najgłębsze uczenie się.