Kulová hvězdokupa s dostatečným množstvím hvězd je metastabilní útvar pomalu se „vypařující“ do prostoru. Simulujte hvězdokupu složenou z cca 200–400 hvězd o hmotnosti rovné hmotnosti Slunce v oblasti o velikosti řádově jednotky parseků. Verzi MD programu s opačným znaménkem interakce „nábojů“ dostanete na požádání, můžete si ho také sami napsat. Dokážete získat útvar dostatečně stabilní alespoň stovek milionů let? Přidejte výklad o aplikaci věty o viriálu na rovnováhu klastru. Literatura: V. Vanýsek: Základy astronomie a astrofyziky, Academia, Praha (1980). ------------------------------------------------------------------- Pokud si chcete simulaci napsat, doporučuji najít vhodnou numerickou metodu s adaptivním krokem pro přesný a efektivní popis blízkých průletů. Chcete-li si napsat sami i integrační metodu, doporučuji integrační metodu Runge-Kutta 4. řádu pro rovnice 2. stupně (Wilipedie, Ralston: Základy numerické matematiky). Adaptivní krok znamená, testujete poloviční a dvojnásobný časový krok a podle toho, jak výsledky souhlasí, měníte krok. Pokud budete používat MACSIMUS, dostanete na požádání „gravitační“ verzi cookstars, kdy se stejné náboje přitahují, i definiční soubor silového pole. Není příliš pravděpodobné, že se dvě hvězdy srazí, rozhodně však tento jev nebudeme modelovat. Ale aby v tomto případě simulace používající leap-frog nehavarovala, interagují „hvězdy“ odpudivým potenciálem WCA-LJ (Lennard-Jones useknutý v minimu a posunutý). Pokud bude typická vzdálenost atomů/hvězd mnohokrát větší než jejich „atomový poloměr“, budou srážky poměrně vzácné (jinak dostaneme něco jako „core collaps“ pozorovaný v některých hvězdokupách). Pokud budete jako uživatel guest na klastru, cookstars vám na požádání instaluji. Pokud jste si MACSIMUS instalovali sami, proveďte následující krok: Překlad vhodné verze (pokud nemáte cookstars z jiného zdroje): cd (kam jste MACSIMUS rozbalili)/macsimus/cook ./configure stars všechny otázky potvrdit Enter KROMĚ: ===== Select BOUNDARY CONDITIONS: f = Free (vacuum) ===== Select charge-charge interaction: g = Gravity ===== The NON-BONDED forces are: wcalj = Lennard-Jones cut at the minimum and shifted Simulační software nastavuje celkovou hybnost i moment hybnosti na nulu, zachování energie není korigováno. Obojí lze na přání ovládat (viz drift a tau.E v manuálu). ============ 0 =============== Protože cook je chemický software, musíme si převést škály. „Hvězdou“ bude atom Ar s „imaginárním nábojem“ 1 e. Doporučuji: 1 LY (světelný rok) → a=100 Å, typická vzdálenost mezi atomy/hvězdami 1 MS (hmotnost Slunce) → mAr, hmotnost atomu argonu VÁŠ PRVNÍ ÚKOL – časová škála. Vypočtěte: 1) Orbitální periodu hvězdy (o hmotnost MS) okolo pevné hvězdy o hmotnosti MS ve vzdálenosti LY. K výpočtu použijte Newtonův gravitační zákon. 2) Orbitální periodu atomu Ar o náboji e a hmotnosti mAr okolo pevného náboje e ve vzdálenosti a. K výpočtu potřebujete Coulombův zákon se změněným znaménkem interakce (stejné náboje se přitahují) a mAr. Poměr vám umožní převést „atomové“ časové údaje na „hvězdné“. Program cookstars generuje po zadání init="random" Gaussovsky rozložený mrak se σ=L/2, kde L je strana krychle takové, že hustota je daná parametrem rho (zadaným do programu v kg/m3) a rozdělení rychlostí podle Maxwellova–Boltzmannova rozdělení se zadanou teplotou T (v K). Protože jsme na začátku předpokládali, že je jedna hvězda na 1 LY, což odpovídá a=100 Å, platí L=a*N^{1/3). Gaussovský mrak NENÍ rovnovážným rozdělením hustoty hvězd v hvězdokupě, bude nutno provést relaxaci (zrovnovážnění)! Nutno tedy spočítat: HUSTOTA rho: Spočtěte z předpokladu, že je jeden atom Ar/a^3 (to odpovídá 1 hvězda/LY^3). TEPLOTA T: Odhadneme řádově z věty o viriálu. Podle té platí pro kinetickou a potenciální energii: 2 Ekin + Epot = 0. Z ekvipartičního teorému: 2 Ekin = 3*N*kB*T, kde N=počet částic a kB=Boltzmannova konstanta. Epot odhadneme v chemických jednotkách na základě úvahy, že průměrná hvězda o „hmotnosti“ (totiž náboji) e je v průměrné vzdálenosti σ od jiné hvězdy, která má „hmotnost“ také e: E1=−e^2/(4*π*ε0*σ), kde ε0=permitivita vakua. Celkem máme N^2/2 párů, tedy: Epot=E1*N^2/2, z čehož plyne odhad rovnovážné teploty: T = N*e^2/(24*π*ε0*σ*kB) (Nezapomeňte, že σ závisí na N.) Výpočet platí jen řádově, zkušenost ukazuje, že takto vypočtené T musíme ještě podělit asi 4, jinak by se nám dost hvězd vypařilo. ČASOVÝ KROK h: Hvězdy mají při přiblížení velkou rychlost. Astronomický software to řeší proměnnou délkou kroku, zde nutno pesimisticky zvolit krok velmi krátký, podle mých krátkých testů lze doporučit h=0.002/sqrt(T) kde h je v ps a teplota v K. (Rychlost je úměrná odmocnině teploty, proto dělíme sqrt(T).) V každém případě SLEDUJTE CELKOVOU ENERGII, neměla by se měnit o vic než několik %. Tyto výsledky doplníte do vstupního souboru (např.) 100stars.def. Nejste-li si jist/a výpočty, přijďte na konzultaci! ========== 1 =========== ----------- Soubor AR.mol -------------- molecule STAR ! total charge = 0.000000 parameter_set = charmm21 number_of_atoms = 1 atoms ! i Crad/atom-id a-type charge chir nb bound_atoms 0 Y999/STAR AR 1.000000 0 0 --------------------------------------- Pozn.: Y999 značí žlutou kuličku o max. velikosti 9.99 Å (aby to bylo dobře vidět). Vytvoření silového pole: blend -o star AR Vytvoří se soubor "star.ble" --------------- soubor 100stars.def ---------- N[0]=POCET ! pocet atomu=hvezd: zacnete treba 100 LJcutoff=-1 corr=4 ! NEMENIT - vyzadovano potencialem WCA rho=DOPLNIT ! dle vypoctu [kg/m3] T=DOPLNIT ! dle vypoctu [K] ! pokud se hvezdokupa rozleti, zmensit ! pokud hvezdokupa kolabuje, zvetsit Emax=T pins=0 ! pro algoritmus nahodneho vkladani h=0.003/sqrt(T)! doporucena hodnota kroku simulace [ps] noint=20 ! po noint krocich se zapisuje soubor .cp ! kratsi noint znamena jemnejsi grafy casoveho vyvoje, ale delsi soubory dt.plb=DOPLNIT ! [ps] aby odpovidalo cca statisicum let, ~300 nasobek h*noint dt.prt=dt.plb ! [ps] - jen pro kontrolu; kratsi dt.prt zpusobi velke soubory ; --------------------------------------- --------------- soubor 100stars.get ---------- init="random" ! nahodny Gaussovsky mrak no=100000 ! upravit podle rychlosti simulace a cílového času ;;;;; --------------------------------------- Start simulace: cookstars star.ble 100stars Prohlížení trajektorie (lze v průběhu simulace): show 100stars show -I%i 100stars Přerušení běžící simulace: touch 100stars.stp Veličiny v závislosti na čase (po skončení simulace): showcp -p -b1 100stars ============= 2 =============== Analýza: 1. Tlak (sloupec 5=P v 100stars.cpa) se v simulaci počítá „náhodou“ skoro stejně jako tlak. Totiž, v periodických okrajových podmínkách pro tlak platí (viz http://old.vscht.cz/fch/cz/pomucky/kolafa/simul07.pdf): P = 2Ekin/3V − ∂Epot/∂V, kde derivace podle objemu se provádí tak, že konfiguraci nafukujeme zachovávajíce relativní vzdálenosti. Protože elektrostatické a gravitační síly obsahují v Epot členy 1/r, vyjde, že ∂Epot/∂V = −Epot/3V. Celkem tedy P = 2Ekin/3V + Epot/3V. Ve kartézském prostoru ovšem nemáme definované V, program nicméně dosadí V spočtené z formální hustoty. Jinými slovy, graf „tlaku“ (sloupec 5) je úměrný viriálu a v rovnováze má kolísat okolo nuly. Zobrazte závislost tlaku na čase a okomentujte. (Pokud vám několik hvězd „uletělo“, bude viriál trochu nad nulou.) 2. Celková energie Epot+Ekin (sloupec 1) by měla být konstantní. Zobrazte závislost na čase a okomentujte. Pokud by grafy byly moc jemné, můžete vytvořit soubor zprůměrovaný po blocích délky třeba 10 takto: showcp -a -b10 100stars