čtvrtek 13. srpna 2009

Srážka rotujících n-rozměrných koulí a termodynamika

Srážka rotujících n-rozměrných koulí je něco, co mi vždycky vrtalo hlavou. Ačkoliv to někteří z vás možná budou považovat za triviální problém, našel jsem v něm celkem netriviální aspekty, o které se v tomto příspěvku podělím. Jaké je přesně zadání?


Mějme soubor dokonale tuhých koulí o daných poloměrech. Koule mohou také obecně rotovat. Sráží se tak, že se zachovává energie, hybnost i moment hybnosti. Najděme takový přepis působících sil, který bude při dostatečně velkém souboru takových koulí imitovat „správnou“ termodynamiku.
Co myslím tou „správnou“ termodynamikou? V tom se právě ukáže být háček, ale původně jsem měl namysli splnění ekvipartičního principu, který se většinou považuje za hodně obecný a v klasické mechanice všeobecně platný.

Dynamika
Jak srážet nerotující koule je poměrně jednoduché.
  1. (1) Převedeme do těžišťové soustavy.
  2. (2) Rozložíme hybnost na směr na střed a na něj kolmý směr.
  3. (3) U složky na střed otočíme znaménko.
  4. (4) Převedeme do původní soustavy.
Takováto simulace automaticky ekvipartiční princip splňuje, i pokud máme různé hmoty koulí. To nás příliš nepřekvapí (i když by mělo..).

Jak je to ale s rotací? Jak ji vlastně reprezentovat? Většinou se k popisu používá (pseudo)vektor úhlové rychlosti. Pro zobecnění do n-rozměrů ale není příliš vhodný - to je vidět už když zkoumáme rotaci dvojrozměrného disku - ten má přece jen jednu úhlovou rychlost, ne dvě, jak bychom od vektoru očekávali! To je dané tím, že vektor úhlové rychlosti není vektor, nýbrž pseudovektor (jako všechny objekty vytvořené vektorovým součinem) - jde o tzv. Hodgeův duál „správně“ invariantního objektu - antisymetrického tenzoru druhého řádu.

To ostatně můžeme ukázat jednoduchou analogií - co je rychlost? Snadno nahlédneme z Taylorova rozvoje funkce polohy s(t), že jde o první derivaci, vektor, který stojí před členem dt:

s(t) = s(0) + v(0) dt + ...

Pokud je rotace daná maticí otočení A(t), která vzniká obecným součinem n vhodně zvolených generátorů rotace kolem jednotlivých os:


,




,



 kde v osách jedné roviny jsou vždy na příslušných místech sinus a kosinus a zbytek vypadá jako jednotková matice. Pokud úhly pootočení závisí na čase jako φ(t), ψ(t) a rozvineme vše do Taylorova rozvoje do prvního řádu v dt, rozvité generátory vypadají jako


,




,



Pokud tedy rozvineme A(t) do prvního řádu v čase, můžeme získat jedině to, co bychom získali nějakým součinem rozvitých generátorů Gij. Všimněme si, že nediagonální členy těchto matic jsou antisymetrické. Tuto vlastnost bude mít do prvního řádu v dt i každý jejich součin (protože na diagonálu se dostanou jen násobky dvou mimodiagonálních elementů, alespoň dt2 , které zanedbáme + to chce trochu rozmýšlení). Proto se dá napsat Taylorův rozvoj matice otočení jako:

A(t) = A(0) + ω(0) dt,

kde A(0) je původní pootočení (sem jsme posbírali všechny členy bez dt) a ω(0) je (nutně antisymetrická) matice úhlové rychlosti (sem jsme posbírali všechny členy úměrné dt). To je jen náhled. Analogicky se moment hybnosti nedefinuje v n-rozměrech jako (pseudo)vektor, ale jako tenzor druhého řádu vytvořený jako

Lij = xi pj - xj pi.

Tento tvar je podezřele známý z kvantové mechaniky. Důležité je, že při srážce se zachovává. Nyní tedy můžeme shrnout to, co děláme, do bodů:
  1. (1) Převedeme do těžišťové soustavy.
  2. (2) Rozložíme hybnost na směr na střed a na něj kolmý směr.
  3. (3) U složky na střed otočíme znaménko.
  4. (4) Z kolmé složky hybnosti a poloh těžišť spočteme moment hybnosti Lij = xi pj - xj pi
  5. (5) Libovolně vyměníme moment hybnosti mezi momenty hybnosti obou míčků (1lij , 2lij) a momentem hybnosti Lij daném pohybem těžišť v těžišťové soustavě.
  6. (6) Vhodnou inverzí soustavy z bodu 4. spočteme z Lij novou hybnost kolmou na střed pro každý z míčků (každý z nich dostane moment hybnosti nepřímo úměrný své hmotě, aby se zachovala i hybnost - v tom, jak rozdělovat moment hybnosti mezi ně volnost nemáme, narozdíl od nezávislých rotací míčků).
  7. (7) Přidáme dvojici sil ve směru spojnice středů mířící od nich. A to takové velikosti, aby výsledná energie byla stejná jako před srážkou. (Touto dodatečnou silou se zjevně nezmění moment hybnosti vůbec ničeho.)
  8. (8) Převedeme do původní soustavy.
Termodynamika
Podobnou simulaci jsem si provedl - výměnu v bodě (4) jsem zkusmo specifikoval tak, že momenty hybnosti se mění tak, aby se vyrovnaly. (To mi přišlo přirozené - termodynamická rovnováha systému koneckonců k takovému výsledku vede.) Byl jsem ale silně překvapený, zjevně nebyl splněn ekvipartiční princip - rotační stupně volnosti měly podstatně méně energie než by jim příslušelo! Docela dlouho jsem bádal nad tím, čím je to způsobené (ekvipartiční princip jsem si zafixoval jako obecné pravidlo platné skoro pro všechny systémy). Přišlo mi hodně zajímavé, že kromě jednoduché kinematiky modelu a integrálů pohybu je potřeba splnit další požadavky, abychom dostali „korektní termodynamiku“. Jedniné řešení, které jsem ale vymyslel je takové, že síly budou zvoleny tak, aby energie každé srážky byla po odrazu vybrána podle Boltzmannova rozdělení - což, jak jistě uznáte, je trochu podvod .

Ve skutečnosti jsem se ale mýlil - se spolužáky jsme prošli předpoklady odvození ekvipartičního principu a zjistili jsme, že zdaleka není tak obecný, jak jsem myslel. Obecně se odvozuje ze střední hodnoty derivace Hamiltoniánu násobené obecnou souřadnicí, podle které derivujeme, ale známá formulace „každý stupeň volnosti má stejnou střední energii“ platí jen pro Hamiltoniány kvadratické v potenciálu či hybnosti. (Sem mj. spadá v dobrém přiblížení v podstatě celá chemie (pokud nejsou stupně volnosti kvantově zamrzlé), takže není divu, že je tato formulace často popisována jako tolik obecná.)

Popsaný model navíc z žádného Hamiltoniánu nevychází - je to jenom přepis pro síly. Takže bychom vlastně neměli očekávat nic - je spíš zajímavou otázkou, proč dokud nezahrneme rotace, ekvipartiční princip vůbec platí?

P.S.
Při přemýtání nad tímto problémem a následnou diskusí jsem narazil na zajímavá skripta z kinetiky, zde, zde, zde, zde a zde.

10 komentářů:

Unknown řekl(a)...

Tedy, v tom formalismu s rotacemi se takhle hned neorientuju, nicméně k těm nerotujícím koulím:

Vůbec mi nepřijde zvláštní, že nerotující koule splňují ekvipartiční teorém. Ekvipartiční teorém je přece splňen tehdy, pokud má každá složka hybnosti každé koule normální rozdělení s rozptylem kT. No ale žejo, to zní docela rozumně, že složky hybnosti budou mít normální rozdělení, nebo ne? :) Hlavně tedy pro ten mikrokanonický soubor, je spíš nutnost, aby se ty koule srážely, a potom teprve bude platit ekvipartiční teorém. Protože bych tak intuitivně řekl, že právě normální rozdělení rychlostí bude v tom vzniklém srážkovém chaosu to nejpravděpodobnější (tj. rovnovážné), a nějaký statistik jistě dodá hezký důkázek proč.

Dokonce bych řekl, že to odvození na wiki, ač se tak netváří, tímto statistickým důkazem je.

irigi řekl(a)...

> Ekvipartiční teorém je přece splňen tehdy, pokud má každá složka hybnosti každé koule normální rozdělení s rozptylem kT.

No, ono rozdeleni nebude normalni, v lepsim pripade bude Boltzmannovo. (Mohlo by byt libovolne, dokud je jeho stredni hodnota 1/2 k T, ekviparticni teorem bude splnen, ale termodynamicky spravne rozdeleni je Boltzmannovo.) Ale uz to, ze je rozdeleni Boltzmannovo, mi prijde prekvapive. Normalne bych takovou vec pro mikrokanonicky soubor odvozoval z jeho Hamiltoninu. (Pro takovou formulaci plati to odvozeni na wiki.)

Ale system micku s temi silami nevychazi z Hamiltonianu, jen z prepisu pro sily. Proto mi prijde prekvapive, ze je rozdeleni Boltzmannovo a e. princip plati, i kdyz (mozna) neni splneny tento predpoklad, ze ktereho jsou odvozeny.

Unknown řekl(a)...

Nerozumím termínu "přepis pro síly", vysvětlit prosím :)

Jinač, mluvím o jednotlivých složkách hybnosti, ne o hybnosti jako celku.

> Mohlo by byt libovolne, dokud je jeho stredni hodnota 1/2 kT, ekviparticni teorem bude splnen

Střední hodnota rozdělení složek hybnosti musí být 0! jinak celé ty koule potečou v nějakém směru a to asi není dobře.

V těch názvech je pak taky trochu guláš. Říďme se wiki. Zaprvé tu máme boltzmannovo rozdělení (dále B), které chápu jako rozdělení, určující pravděpodobnost, že kanonický systém bude mít nějakou energii E, při kontaktu (a rovnováze) s rezervoárem o teplotě T. To se nás podle mě vůbec netýká, protože máme mikrokanonický soubor.

Pak tu je Maxwell-Boltzmannovo rozdělení (dále MB), které je právě uváděno u plynů apod. a popisuje rozdělení _velikostí_ rychlostí částic. A důležitý fakt: není to nic jiného, než rozdělení velikostí náhodného vektoru, jehož všechny složky mají normální rozdělení s nulovou střední hodnotou a stejným rozptylem.

Tedy, aby skupina těch koulí měla MB rozdělení rychlostí, stačí aby měly jejich hybnosti v každé složce normální rozdělení se stejnými rozptyly (které nazveme kT). A taková intuice různých centrálních limitních vět mi říká, že náhodně se srážející koule nejspíše složky hybnosti normálně rozdělené mít budou (a proto mi to nepřijde překvapivé). A nepochybuju o tom, že existuje i nějaký pěkný rigorózní důkaz. A nakonec ani žádné MB rozdělení zavádět nepotřebujeme. To je jen další vlastnost, že normálně rozdělené složky hybnosti generují rozdělení velikostí hybnosti jako MB. Ale už jen to normální rozdělení po složkách, má za důsledek splnění ekvipartičního teorému.

Rozptyl rychlostí jsme nazvali kT/m, tudíž 1/2 m < v^2 > = 1/2 kT pro každou složku.

A svými slovy o MB jsem si téměř jist minimálně v 3D :)

Unknown řekl(a)...

Aha, už možná chápu co jsi měl na mysli. I v tom případě co popisuju já, tj. složky hybnosti rozdělené normálně, velikost hybnosti podle MB, bude mít energie částic rozdělení B, tj navíc splňující ekvipartiční teorém. Nicméně původ to má v té normalitě složek hybnosti, kterou skromně považuji za očekávanou :)

irigi řekl(a)...

Teď jsem si uvědomil, že jsem Tvoji úvahu o normálním rozdělení odmítl neprávem. Jsou to totiž IMHO vlastně dvě rovnocenné věci. Mohu se na to dívat tak, že:

1a) Celkově mám mikrokanonický soubor, každá kulička v něm se ale chová jako kanonický soubor v kontaktu s rezervoárem. Na každé své složce hybnosti tedy bude mít B rozdělení e^(-E/kT) = e^(-p_x^2/2mkT).
2a) MB rozdělení vzniká tak, že místo rozdělení složek zavedeme rozdělení velikostí. Důvod, proč nemůžeme rovnou MB rozdělení počítat z B rozdělení je, že větší rychlosti mají větší míru degenerace (je více stavů s vyšší rychlostí, protože povrch koule v prostoru rychlostí roste).

Tak jsem se na to díval já. Ve skutečnosti ale v případě jednosložkového rozdělení mohu použít to, co zjevně používáš Ty:

1b) Rozdělení hybností je normální e^(-p_x^2/2mkT) (krát norma), protože vzniká složením mnoha vlivů. (To je tady shodou okolností stejné jako B.)
2b) To samé co 1a)

Takže v tomhle už si snad rozumíme. Teď si myslím, že ani jedno ze dvou odvození nevysvětluje, proč se plyn kuliček s různou velikostí a hmotností (bez rotace) řídí ekvipartičním principem, protože:

a) Pokud je každá kulička kanonický soubor, vyplývá odvození B-rozdělení z maximalizace entropie přes energie (tady OK), a z faktu, že mikrokanonický soubor systém + rezervoár je popsaný nějakým Hamiltoniánem a že z Hamiltonových pohybových rovnic a dalších věcí plyne, že střední hodnota {x_j partial H / partial x_i} = delta_ij k T . Pokud ale moje simulace běží tak, že při nárazu akorát vygeneruji nějakou sílu, pak nemůžu tvrdit, že vycházím z nějakého Hamiltoniánu, nebo že je ten systém vůbec nějaký Hamiltonián má (jinak bych musel síly počítat z něj, ne je prostě předepsat).

b) IMHO taky nefunguje, protože nenajdu žádný důvod, proč by rozptyl normálního rozdělení pro různě těžké míčky měl být stejný. (Pokud by byly stejné, celkem věřím, že rozdělení bude normální a že rozptyl nějaký být musí, takže proč ne třeba k T / 2 m (nebo takový, který rozdělení výš má). Ale nenacházím žádný důvod, proč by k nemělo záviset na hmotnosti a velikosti míčku (kromě fyzikálního, že je to Boltzmannova konstanta, ale to by byl důkaz kruhem..))

Proto si myslím, že by nás mělo překvapovat, že ekvipartiční princip tady platí. (To může být dáno třeba tím, že přecejenom nějaký Hamiltonián simulace má, takže projde a), ale nevidím to.)

Nebo se pletu?

irigi řekl(a)...

P.S. ale jinak bych řekl, že obecnější bude to odvození typu "podsystém je kanonický soubor" - to, že něco má normální rozdělení je sice běžné (protože to rozdělení maximalizuje entropii přes mnoho vlivů), ale to ještě neznamená, že je to nějak fundamentální, nebo že by mělo být všude. Zatímco u B rozdělení mám pocit, že se alespoň odrážím o Hamiltonův formalismus.. Řekl bych třeba, že pro nějaký plyn fotonů nebo relativistických částic normální rozdělení platit nebude.

Unknown řekl(a)...

Ale přece, nějaké k nás vůbec nezajímá :) Musí být konstantní, je to jen nějaké škálování které nám dělá stupně Kelvina stejně velké jako Celsia nebo tak něco.

Předkládám dle svého pohledu až na statistické detaily, plné odvození:

Mám systém koulí, každá má jinou hmotnost, na počátku libovolně náhodně rozdělené rychlosti. Nakonec mě ani nezajímá nějaká jejich poloha nebo tak něco. Klíčové vlastnosti tohohle systému, díky kterým dojde ke splnění ekvipartičního teorému jsou dvě:
1) Zachovává se celková hybnost (tedy i po složkách). Bez újmy na obecnosti je celková hybnost nulová.
2) Srážka dvou koulí probíhá tak, že jedna koule předá část svojí hybnosti jiné kouli.

Ke srážkám koulí dochází naprosto náhodně, a nezajímá mě jak konkrétně, stačí že dochází k náhodnému přečerpávání hybností mezi koulemi. Hybnost se neztrácí.

Vezmu teď libovolnou složku hybnosti. Protože se zachovává, má každá koule nějakou svoji "přihrádku" na tuhle složku hybnosti. Celková hybnost je poté rozdělena mezi jednotlivé přihrádky (v přihrádce jde však s hybností jít i do mínusu). Každá srážka dvou koulí, přesype náhodné množství hybnosti z jedné přihrádky do druhé.

Teď přijde nějaký statistik, a řekne, že nejpravděpodobnější náhodné rozsypání hybnosti do jednotlivých přihrádek je takové, že v každé přihrádce je stejné rozdělení té hybnosti. Protože máme celkovou hybnost búno nulovou, bude střední hodnota patrně nulová, a navíc rozptyl všude stejný, označme ho V. Navíc to rozdělení hybnosti bude všude normální, to statistik taky nějak vysvětlí (použitím CLV? aktuální hodnota v každé přihrádce je součet nějakých přisypání nebo odebrání hybnosti, se střední hodnotou nula).

A nyní již mám co jsem chtěl. Fůru přihrádek na hybnost a v každé té přihrádce je rozdělení té hybnosti stejné a navíc normální, a může za to ta fyzika, a sice zachování hybnosti a mechanismus srážky, tj. neustálé náhodné přelévání hybností mezi přihrádkami.

Rozptyl hybností ve všech přihrádkách je tedy nějaké V.

Rychlosti koulí mají potom rozdělení s rozptylem V/m.

Střední hodnota kinetické energie každé koule potom bude N*1/2*V, kde N je počet dimenzí.

Nyní přijde teoretický fyzik, a definuje teplotu systému tak, aby jí byla úměrná ta střední hodnota energie. A nejlepší co ho napadne, je přejmenovat to V na T, tedy pojmenovat ten rozptyl hybnosti, který se ustálí všude stejný. A protože má za sebou nějaké historické důvody, a chce aby jednotka jeho teploty byla nějaká rozumně velká, tak to ještě naškáluje konstantou k. A je hotovo.

Unknown řekl(a)...

Ještě doplněk k tomu k. Když se zavádí teplota T statisticky (třeba ve Waldramovi), tak se zavede 1/kT = \parc{ln g}/\parc{E} kde g je hustota stavů. A o k se tam právě mluví jen jako o škálování, aby mi ty stupně Kelvina vyšly nějaké pěkné.

Pak taky, pokud ty koule budu mít v nějaké idealizované krabici, tak se asi hybnost zachovávat nebude (při odrazu koule od krabice). Nicméně, střední hodnota toho nezachování bude asi stejně nulová, protože koule narážejí na všechny stěny přibližně stejně často. Pokud ani to nebude stačit, tak zahrnu do systému i krabici, dam jí třeba nějakou velice vysokou hmotnost, aby se nehýbala, a jedu.

Když už o tom tak přemýšlím, tak vidím že vůbec nemusím uvažovat pouze koule. Jediná podmínka je, aby docházelo k těm náhodným přeléváním zachovávající se hybnosti (a od té byl také jediný příspěvek do energie)

irigi řekl(a)...

Omlouvám se za pozdní odpověď:

Tohle odvození jsem už v minulém příspěvku odsouhlasil, jen s malou výhradou. A to sice:

1) Fakt, že se hybnost rozhazuje podle normálního rozdělení je poměrně dobře uvěřitelné. Není to však IMHO obecně nutné. Třeba v relativitě platí Maxwell-Juttnerovo rozdělení rychlostí, kde rozdělení hybností obsahuje myslím Besselovy funkce. Apriori mi třeba není jasné, proč by se podle normálního rozdělení neměla rozhazovat energie (pokud namítneš že proto, že je zdola omezená, tak třeba hybnost normovaná velikostí na energii, nebo nějaká podobná vymyšlenost, která už zdola omezená není a přitom se chová k oběma míčkům symetricky.)
2.) Dejme tomu, že si bodu 1 nebudu všímat - koneckonců předpokládat normální rozdělení je celkem rozumné. Pak mi pořád není jasné, proč by se míčky s různými hmotnostmi měly řídit normálním rozdělením se stejným rozptylem. Proč nemůže mít každá hmotnost míčků normální rozdělení, ale jejich rozptyl bude různý?

Teď když nad tím přemýšlím se mi zdá, že:

1) Za to, že se normálně rozdělí zrovna hybnosti asi bude moci ekvivalence vztažných soustav - fakt, že při změně vztažné soustavy by měly hybnosti pořád normální rozdělení, jen posunuté, mi připadne silný. Třeba pro energie, nebo jiné vymyšlené veličiny, to platit nebude - pak bychom měli privilegovanou soustavu. V relativitě by roli normálně rozdělené veličiny možná mohla hrát rapidita násobená klidovou hmotou, ale to jenom hádám.
2) Za vyhovění této námitce asi může fakt, že právě hybnost nerozlišuje hmotnosti míčků (narozdíl třeba od rychlosti). Tím, že hybnosti těžkých a lehkých míčků vystupují při srážkách symetricky, se asi námitka 2 řeší.

Takže Ti dávám zapravdu úplně. Ale připadne mi, že tyhle dvě námitky jsou vcelku důležité a nejsou (alespoň mně) zřejmé na první pohled. Jakmile se začneme zabývat srážkami s rotačními stupni volnosti (byť se zachováním všech integrálů pohybu), už nás asi žádná symetrie nezachrání. Proto se z tohohle pohledu nedivím nesplnění ekvipartičního principu v tomto případě.

P.S. Že B konstanta představuje jenom škálování vím, já jsem to myslel jinak, asi ne úplně trefně, tak bych to nechal být.

Jinak díky moc za zajímavou diskusi!

irigi řekl(a)...

P.S. To s tou rapiditou beru zpátky - byl to první nástřel co mne napadl a zrovna jsem si vybral veličinu, která se nezachovává a možná ani správně netransformuje. Ale zdá se mi, že by měla nějaká podobná veličina, která se v relativitě rozdělí normálně a invariantně vzhledem k pozorovateli být. To bych ale musel rozmyslet víc..