14 Uvod Monte Carlo (MC) teoretične simulacije v fiziki in dru- gih vedah nam omogoča v programsko opremo računal- nika vgrajeni generator naključnih števil (GNŠ) [1, 2]. T a nam s ponavljanjem daje vedno nova enakomerno porazdeljena racionalna (praktično lahko rečemo kar re- alna) števila med 0 in 1. Čeprav je algoritem za ta števila Monte Carlo simulacija zračnega upora mag. Saša Harkai in dr. Milan Ambrožič Fakulteta za naravoslovje in matematiko Univerze v Mariboru dr. Marjan Krašna Fakulteta za naravoslovje in matematiko in Filozofska fakulteta Univerze v Mariboru Povzetek Zaradi hitrih osebnih računalnikov so postale Monte Carlo simulacije v teoretični znanosti in tudi v praksi zelo popularne. Z njimi raziskujemo naključne pojave v naravi in družbi, ki so za izključno analitično obravnavo prezap- leteni. Med zelo zahtevne fizikalne probleme spadajo tudi naloge s področja aerodinamike in hidrodinamike. Zato jih navadno rešujejo s kombiniranjem numeričnih simulacij in eksperimentov, na primer z uporabo vetrovnikov. Pri gibanju telesa skozi kapljevino ali plin se pojavi upor sredstva, za katerega pri dovolj majhnih hitrostih dobro velja linearni zakon upora (sorazmernost sile s hitrostjo telesa), pri velikih hitrostih pa kvadratni zakon (sorazmernost sile s kvadratom hitrosti telesa). Oba zakona dosledno izpeljemo iz Navier-Stokesove nelinearne parcialne diferencialne enačbe oziroma iz njenih približkov. Pri gibanju telesa skozi idealni plin pa lahko silo upora ocenimo tudi iz mikro- skopske slike. T ako smo naredili v naši Monte Carlo simulaciji in pri tem med seboj povezali različne vede: statistično termodinamiko, mehaniko in numerične metode. Ključne besede: Monte Carlo, zračni upor, idealni plin, statistična termodinamika Monte Carlo Simulation of Drag Abstract Powerful personal computers made Monte Carlo simulations very popular in theoretical science and in practice. They are used to investigate probabilistic phenomena in nature and society which are too complex to be solved using a purely analytical approach. T asks from the field of aerodynamics and hydrodynamics also belong to the class of very demanding physical problems. Therefore, they are usually solved with a combination of numerical simulations and experiments, e.g. using wind tunnels. When an object moves through liquid or gas it experiences drag, to which the linear law (force is proportional to the velocity of the object) applies for low enough velocities, while for high velocities the quadratic law applies (force is proportional to the square of the velocity of the object). Both laws can be logically derived from the Navier–Stokes nonlinear partial differential equation, i.e. from its approximations. However, when the object moves through an ideal gas, the drag force can also be estimated from the microscopic picture. This has been done in the present Monte Carlo simulation, in the process integrating different sciences: statistical thermodyna- mics, mechanics and numerical methods. Keywords: Monte Carlo, drag, ideal gas, statistical thermodynamics seveda točno določen (determinističen), imamo vseeno vtis, da so števila, ki nam jih daje GNŠ, naključna. S temi števili potem z ustrezno funkcijsko pretvorbo ge- neriramo katerokoli naključno fizikalno spremenljivko z znano verjetnostno porazdelitvijo. MC-simulacije v fiziki in v vedah, povezanih s fiziko, uporabljajo npr. pri teoretičnem preučevanju mehanskih lastnosti ma- Fizika v šoli 15 Strokovni prispevki terialov (porazdelitev trdnosti in drugo), perkolacijski teoriji, raziskavah faznih prehodov tekočih kristalov, fe- roelektričnih in feromagnetnih materialov (razni spinski modeli) itd. Zanimiv primer uporabe MC-simulacij so razni fizikalni problemi v statistični termodinamiki. T o je veda, ki se ukvarja z mikroskopskim opisom gibanja gradnikov snovi v kateremkoli agregatnem stanju in je izhodišče za termodinamične zakone v makroskopski (fenomenološki) termodinamiki [3]. Značilen zgled je obravnava gibanja posameznih molekul idealnega plina v zaprti posodi s homogenim tlakom in temperaturo, iz česar izhaja splošno znana in dokaj preprosta plinska enačba. Ker je statistična termodinamika osnovana na verjetnostnem računu, so MC-simulacije zelo primerne zanjo. V konkretnem primeru posode z ravnimi stena- mi lahko iz mikroskopske slike, to je z obravnavo trkov posameznih molekul ob steno, izračunamo tlak plina na steno – to je pravzaprav bistvo plinske enačbe. Lahko pa naredimo še korak dlje: izračunamo lahko silo plina na telo, ki se giblje skozenj. S tem lahko na mikroskopskem nivoju preverimo veljavnost linearnega in kvadratnega zakona upora, ne da bi morali reševati zapleteno Navier- -Stokesovo parcialno diferencialno enačbo [4]. Naloga je zelo poučna tudi z didaktičnega vidika, predvsem za srednješolce in študente prvega letnika naravoslovnih in tehničnih smeri, saj gre za povezavo različnih vej fizike pa še za uporabo informacijsko-komunikacijske tehno- logije [5-7]. Učitelj lahko to simulacijo opiše in prikaže npr. pri obravnavi gibalne količine in trkov ali pa pri ter- modinamiki. Linearni in kvadratni zakon upora Pri majhnih hitrostih gibanja telesa skozi tekočino (plin ali kapljevino) gre za silo upora zaradi viskoznosti sred- stva [8]. Vzrok za to silo je trenje med plastmi tekočine v okolici telesa, ki se relativno gibljejo ena glede na drugo. Tik ob površju telesa se tekočina namreč lepi na telo in se giblje skupaj z njim, dovolj daleč od njega pa tekočina miruje; vmes imajo zato različne plasti tekočine različne hitrosti. Silo upora tekočine na telo zapišemo v tem pri- meru v splošnem z enačbo: (1) Pri tem je L značilna linearna dimenzija telesa, η viskoz- nost tekočine, v hitrost telesa, k pa ustrezen številski koe- ficient, med drugim odvisen od oblike telesa. Če gre npr. za gibanje kroglice s polmerom R skozi tekočino, vzamemo v enačbi (1) L = R in k = 6p ter zapišemo Stokesovo enačbo: F u = 6pRηv. Pri velikih hitrostih pa linearni zakon upora, enačba (1), ni več uporaben. Za- radi večjih medsebojnih relativnih hitrosti plasti tekoči- ne hitro pride do turbulentnega toka, plasti tekočine se »odlepijo« od telesa itd. Za upor tekočine na telo ni več bistvena viskozna sila, temveč zastojni tlak, kar je neka- ko razlika tlačnih sil na sprednjo in zadnjo stran telesa (glede na smer gibanja) in jo približno izračunamo z Bernoullijevo enačbo. Zastojni tlak in z njim sila upora sta sorazmerna s kvadratom hitrosti telesa [8]: (2) Pri tem je S največji prečni prerez telesa (pravokotno na smer gibanja), r pa gostota tekočine. Vpeljali smo tudi koeficient upora C u , ki je odvisen od oblike telesa in katerega vrednosti za najznačilnejše aerodinamične oblike najdemo na spletu. Katera sila, (1) ali (2), prevla- duje pri neki hitrosti, najelegantneje ugotovimo z nju- nim razmerjem. Sili torej delimo med seboj, v količniku F u2 /F u1 izpustimo številske faktorje, za prečni prerez vza- memo še oceno S ≈ L 2 in dobimo Reynoldsovo število: Re = rLv/η. Če je torej število Re dovolj veliko (navad- no zaradi gotovosti vzamemo oceno Re > 1000), potem dobro velja kvadratni zakon upora (2). Za Re < 1 pa dob- ro velja linearni zakon upora (1). Vmesno območje, 1 < Re < 1000, je glede preproste ocene upora bolj ne- gotovo in si raje pomagamo z numeričnimi simulacijami in eksperimentom. Maxwellova porazdelitev hitrosti molekul v idealnem plinu Obravnavajmo enostavni plin, na primer dušik, s homo- geno temperaturo, npr. T = 300 K. Iz statistične termo- dinamike izhaja, da je povprečna hitrost molekul soraz- merna s korenom absolutne temperature. Pri 300 K je gro- ba ocena hitrosti za molekule dušika (N 2 ): v ≈ 500 m/s. Pri kisikovi molekuli O 2 ne izračunamo bistveno dru- gačne vrednosti; hitrost molekule je sicer odvisna tudi od njene mase, vendar pa je razlika med masama molekul N 2 in O 2 le okrog 14 %. Značilno je tudi to, da vse mole- kule niso enako hitre, temveč so njihove hitrosti naključ- no porazdeljene po Maxwellovi porazdelitvi [9]: (3) Simbol m v enačbi (3) označuje maso molekule, k pa je Boltzmannova konstanta: k = 1,38 · 10 -23 J/K. Funkci- ja p(v) ima pomen verjetnostne gostote: če to funkcijo integriramo med dvema vrednostma hitrosti, v 1 in v 2 , njen integral pomeni verjetnost, da bo imela naključno izbrana molekula hitrost v območju v 1  v  v 2 . Graf funkcije (3) v brezdimenzijski obliki je prikazan na sli- ki 1. Funkcija p(v) je pri majhnih hitrostih kvadratna, p µ v 2 , saj je eksponentni faktor približno enak 1. Pri ve- likih hitrostih funkcija hitro pada, pri neki hitrosti v M pa ima izrazit maksimum. T o hitrost izračunamo z zahte- vo, da je odvod dp/dv enak 0: (4 a) Vrednost funkcije pri tej hitrosti je: (4 b) 16 Simbol e v zadnji enačbi označuje osnovo naravnega logaritma. Funkcijo (3) zapišimo v brezdimenzijski obliki, primerni za grafični prikaz: (5) Omenimo še, da povprečno hitrost molekul izraču- namo takole: (6) Zaradi asimetrije funkcije p(v) sta omenjeni hitrosti različni: ≈ 1,13 v M . Slika 1: Graf verjetnostne gostote v brezdimenzijski obliki pri Maxwellovi porazdelitvi; črtkana navpična črta označuje hit- rost (v M ) z največjo vrednostjo funkcije p(v), polna navpična črta pa povprečno hitrost ( ). Poleg Maxwellove porazdelitve velikosti hitrosti mole- kul pri značilni MC-simulaciji idealnega plina upošte- vamo še, da so smeri vektorja hitrosti molekul izotrop- no porazdeljene, to pomeni, da so vse smeri enakovred- ne. V MC-simulaciji moramo poleg začetnih vektor- jev hitrosti molekul z GNŠ postaviti tudi naključne lege, enakomerno porazdeljene po simulacijski celici. Ko so tako postavljeni začetni pogoji (vektor lege in hitrosti za vsako molekulo), plin prepustimo samemu sebi. Molekule se gibljejo premo enakomerno, razen v trenutku trka s steno ali oviro. Vpliv teže lahko povsem zanemarimo. T o preverimo npr. tako, da primerjamo razliko potencialnih energij molekul na različnih viši- nah v posodi in značilno kinetično energijo. Naj bo viši- na posode 1 m, značilna hitrost molekul pa 500 m/s, kot smo ocenili zgoraj za dušik. Opazujemo torej raz- liko potencialnih energij dveh molekul pri višinski razliki 1 m in jo primerjamo s povprečno kinetično energijo: W p /W k = (mgh)/(mv 2 /2) = 2gh/v 2 ≈ 8 · 10 -5 . Potencialno energijo res lahko zanemarimo. MC-simulacija zračnega upora Namesto tridimenzionalne (3D) obravnave problema smo si zaradi nazornosti raje izbrali dvodimenzionalni (2D) prikaz in simulacijo. Kot simulacijsko celico s pli- nom smo torej vzeli kvadrat s stranico 10 m. Kot smo že omenili, smo na začetku z GNŠ vsaki molekuli priredili začetne pogoje: dve koordinati, x in y, hitrost v na osnovi porazdelitve (3), smer vektorja hitrosti pa smo podali s kotom glede na os x. T o pomeni, da smo za vsako mole- kulo plina potrebovali štiri naključna števila GNŠ. Poleg tega smo v simulacijsko celico dodali majhno ploščico (ozek pravokotnik v dvodimenzionalnem prikazu), ki se je gibala enakomerno z določeno hitrostjo v vodoravni smeri (os x). Pri izpeljavi sile plina na ploščico smo si pomagali s teorijo prožnih trkov, podobno kot pri izpe- ljavi tlaka plina na stene posode. Naj se torej ploščica (telo) giblje s hitrostjo v t = (v t , 0) v desno. V anjo naj od »zadaj« trči molekula plina s hitrostjo v = (v x , v y ), kjer je v x > v t . Ker je trk prožen in ker je molekula neprimerljivo lažja od ploščice, se molekula odbije nazaj s hitrostjo v' = (2v t – v x , v y ). T a hitrost izhaja iz zakona o ohranitvi kinetične energije in gibalne količine. Sunek sile mole- kule na ploščico v smeri osi x je torej: F t = m(v x – v x ') = 2m(v x – v t ). Molekule, ki trčijo v ploščico od »spre- daj«, v povprečju delujejo na ploščico z večjo absolutno vrednostjo sunka sile kot tiste, ki jo zadenejo od zadaj. Zanje velja enaka enačba kot za trke od »zadaj«, če prav upoštevamo predznake komponent hitrosti. V času t naj ploščico zadene M molekul, ki jih označimo z indeksom i. T edaj je sila na ploščico: (7) Čeprav je ta sila negativna, saj sila na ploščico deluje v nasprotno smer od njenega gibanja, smo na grafih silo upora prikazali kot pozitivno količino. Ker je simulacij- ska celica (posoda) končna, moramo upoštevati tudi, kaj se zgodi, ko molekula prileti do mejne ploskve. Uporabi- li smo dva različna načina obravnave takšnih primerov. Pri prvem načinu spremenimo predznak komponente hitrosti v tisti smeri, v kateri je molekula zapustila mejo celice. Če npr. molekula prekorači mejo, ki je pravokotna na os x, se komponenti v x po trku samo zamenja pred- znak, komponenta v y pa ostane nespremenjena. S takim načinom obravnavamo mejo kot steno zaprte posode. Drugi način pa je, da molekulo premaknemo na nasprot- no stran meje, ki jo je presegla: če npr. molekula preko- rači desno mejo, pravokotno na os x, potem jo premak- nemo na levo mejo pravokotno na os x na enaki višini. T ako simuliramo neskončni sistem, v katerem je simu- lacijska celica osnovna celica, ki se neskončno ponavlja v obe smeri. Po pričakovanju se je izkazalo, da dobimo v obeh primerih enake vrednosti sile upora. Uporabili smo programski paket Microsoft Visual Stu- dio 2010 v programskem jeziku C#. Izbrali smo število molekul N = 10 4 v kvadratu s stranico a = b = 10 m. Če Fizika v šoli 17 Strokovni prispevki si npr. mislimo, da je tretja stranica kvadra, ki je nismo prikazali, enaka c = 1 m, potem je prostornina kvadra V = a 2 c = 100 m 3 in je številska gostota molekul na pro- storninsko enoto: n = N/V = 100/m 3 . T o je majhna vred- nost v primerjavi z realnimi razmerami, kjer je navadno n velikostnega reda 10 25 /m 3 . Lahko bi vzeli npr. kvader s stranicami a = b = 1 mm, c = 1 mm, pa bi takoj imeli večjo prostorninsko gostoto molekul, n = 10 16 /m 3 , a še vedno precej manj od 10 25 /m 3 . Lahko izberemo renor- malizirano (precej večjo) maso molekul, kot je v resnici, pa dobimo enak učinek pri manjši prostorninski gosto- ti molekul. Izbrali smo vrednost mase m = 10 -10 kg. A tu nam ne gre za absolutne vrednosti fizikalnih količin, posebno sile na ploščico, temveč za to, kako se relativna vrednost sile spreminja s hitrostjo. Lege molekul (točk v modelu) in ploščice v nekem trenutku simulacije so pri- kazane na sliki 2, vrednosti sil pa na slikah od 3 do 5. Na sliki 3, kjer so hitrosti ploščice manjše od 200 m/s, je graf F(v t ) prikazan v navadni obliki in linearna odvisnost sile od hitrosti je očitna. T orej za v t < 200 m/s dobro velja linearni zakon upora. Na sliki 4 je prikazan graf F(v t ) za hitrosti ploščice, večje od 500 m/s. Graf ni linearen. Ali je F(v t ) potenčna funkcija, najlaže preverimo tako, da graf narišemo v logaritemskem merilu (slika 5). Linearna od- visnost v logaritemskem diagramu res dokazuje potenč- no odvisnost sile od hitrosti, smerni koeficient premice, ki podaja eksponent, pa je v našem računu k = 1,97 ± 0,02. T o je res zelo blizu pravi vrednosti 2, torej v tem območju hitrosti odlično velja kvadratni zakon upora. Slika 2: Trenutna razporeditev 10 4 molekul plina v 2D modelu in lega ploščice (tanek pravokotnik), ki se nato začne enakomerno gibati v desno. Simulacijo najdete na spletni strani: http://distan- ce.pfmb.uni-mb.si/av/fizika/idealni_plin_simulacija1.mp4. Sila, ki jo program prikazuje levo spodaj, je trenutna sila pri zgoraj določenih podatkih (npr. pri dani hitrosti ploščice), katere vred- nost zaradi statističnih fluktuacij nekoliko niha. Slika 3: Graf odvisnosti sile upora od hitrosti ploščice: sila je nor- malizirana glede na silo pri hitrosti 200 m/s. Točke pomenijo izra- čunane vrednosti sile, premica med njimi pa je najboljši približek linearne funkcije. Slika 4: Graf odvisnosti sile upora od hitrosti ploščice pri večjih hitrostih: sila je normalizirana glede na silo pri hitrosti 500 m/s. Slika 5: Graf odvisnosti sile upora od hitrosti ploščice pri večjih hitrostih v logaritemski skali. Že na pogled je videti, da je koefi- cient ustrezne premice blizu vrednosti 2. 18 Hitrosti na slikah 3 do 5 so pretirane, saj so npr. na sli- kah 4 in 5 največje hitrosti ploščice neprimerno večje od povprečne hitrosti molekul plina N 2 pri navadnih tem- peraturah. Opozoriti pa moramo še na nekaj. Čeprav smo z MC-simulacijo ugotovili, da je pri dovolj majhnih hitrostih objekta sila upora sorazmerna s hitrostjo, pa ta simulacija vseeno ni prava razlaga linearnega zakona upora (1). S simulacijo smo namreč upoštevali trke mo- lekul ob telo in npr. od stene posode, ne pa tudi medse- bojnih trkov molekul. Za viskoznost idealnega plina pa so pomembni ravno slednji, saj določajo »notranje tren- je« med molekulami zaradi različnih povprečnih hit- rosti plasti plina ob telesu. Viskoznost plina je odvisna od efektivnega premera molekul za medsebojne trke, v naši MC-simulaciji pa smo računali, kot da so molekule točkaste in se vedno »zgrešijo«. Da pri majhnih hitrostih ploščice pridemo do linearne- ga zakona upora (namesto kvadratnega zakona), ki pa nima nobene zveze z viskoznostjo, je mogoče dokazati s preprosto analitično oceno in z uporabo enačbe (7), če upoštevamo, da imajo molekule, ki zadenejo ploščico od zadaj, v povprečju manjšo relativno hitrost glede na ploš- čico kot tiste, ki ploščico zadenejo od spredaj. Mogoče je torej, da imajo direktni trki molekul ob telo tudi prispe- vek k linearnemu zakonu upora, a je ta prispevek manjši kot tisti zaradi viskoznosti. Sklep Pokazali smo, da je z Monte Carlo simulacijo, s kate- ro obravnavamo mikroskopsko gibanje idealnega plina, mogoče na poučen način preveriti linearni zakon upora pri majhnih hitrostih telesa in kvadratni zakon upora pri velikih hitrostih. Uporabili smo preprostejši dvodimen- zionalni model, čeprav bi bila razširitev na tri dimenzije relativno enostavna. Če pa bi hoteli izračunati koeficient upora v enačbi (2) za kvadratni zakon upora za telesa znanih aerodinamičnih oblik, bi seveda morali računati v treh dimenzijah. Zavedati se moramo tudi, da opisana simulacija ne daje realistične razlage za linearni zakon upora, saj viskoznost plina razložimo z medsebojnimi trki molekul, ki jih tu nismo mogli upoštevati. Molekule smo namreč obravnavali kot točkasta telesa. Naj ome- nimo še, da smo se lotili tudi simulacij v treh dimenzi- jah, npr. pri gibanju elipsoida s tremi različnimi polosmi skozi plin. T a problem je matematično zahtevnejši, npr. zato, ker je sunek sile molekule plina na elipsoid odvisen od mesta, kjer molekula zadene elipsoid. Prvi rezulta- ti potrjujejo, kar smo že ugotovili za dvodimenzional- ni sistem: pri majhnih hitrostih elipsoida velja linearna odvisnost sile upora od hitrosti, pri velikih hitrostih pa kvadratna. Prikaz dvodimenzionalne simulacije s programom, ki ga lahko prenesete na svoj računalnik, je na spletni strani: http://distance.pfmb.uni-mb.si/av/fizika/idealni_plin_ simulacija1.mp4. Program za realno simulacijo na vašem računalniku (za potrpežljivega bralca, ker ta simulacija teče počasneje) v stisnjeni obliki pa najdete na spletni strani: https://www.dropbox.com/s/tduf6l07g5gi2wx/Ideal- ni%20plin.zip?dl=0. Viri [1] https://en.wikipedia.org/wiki/Monte_Carlo_method [2] Ambrožič, M. (2013). Računalniška simulacija naključnosti. Življenje in tehnika, maj 2013, l. 64, št. 5, str. 38-44. [3] https://en.wikipedia.org/wiki/Statistical_mechanics [4] https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations [5] Gerlič, I. (2000). Sodobna informacijska tehnologija v izobraževanju. Ljubljana: DZS. [6] Svetlik, I. (2006). O kompetencah, Vzgoja in izobraževanje. 1/2006. Ljubljana: Zavod RS za šolstvo. [7] Grubelnik, V. (ur.) (2010). Opredelitev naravoslovnih kompetenc (znanstvena monografija projekta RNK). Maribor: Fakulteta za naravoslovje in matematiko. [8] Strnad, J. (2003). Mala fizika 1, Učbenik za pouk fizike v 1. in 2. letniku gimnazij in srednjih šol. Ljublja- na: DZS. [9] https://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution