蒙特卡羅方法生成三維理想氣體麥克斯韋分布

蒙特卡羅方法生成三維理想氣體麥克斯韋分布莫裕宇;李盛至【摘要】蒙特卡羅方法生成二維麥克斯韋分布的程序在一些計(jì)算物理教材中已經(jīng)做出.將三維單位球面上均勻分布的隨機(jī)向量加入到描述分子碰撞的隨機(jī)過程中，用MATLAB生成了三維麥克斯韋分布的程序.并對(duì)麥克斯韋分布的數(shù)值結(jié)果與理論預(yù)期進(jìn)行了對(duì)比.對(duì)應(yīng)地提出了一種理想氣體碰撞生成麥克斯韋分布的在統(tǒng)計(jì)上所可能遵循的散射規(guī)律.【期刊名稱】《大學(xué)物理》【年(卷)，期】2018(037)003【總頁(yè)數(shù)】6頁(yè)(P63-68)【關(guān)鍵詞】三維麥克斯韋分布數(shù)值模擬;蒙特卡羅方法;理想氣體【作者】莫裕宇;李盛至【作者單位】北京師范大學(xué)物理學(xué)系，北京100875;北京師范大學(xué)物理學(xué)系，北京100875【正文語(yǔ)種】中文【中圖分類】O552.3+1理想氣體的麥克斯韋分布是大量粒子相互作用的結(jié)果，如果能夠把這樣的實(shí)際物理過程簡(jiǎn)化為簡(jiǎn)單的隨機(jī)過程，那么當(dāng)需要用到速率分布滿足三維麥克斯韋分布的粒子進(jìn)行分子模擬時(shí)，計(jì)算量將會(huì)大大簡(jiǎn)化.目前蒙特卡羅方法生成麥克斯韋分布已有得的相關(guān)研究為文獻(xiàn)［1,2］.在文獻(xiàn)［1］中分子運(yùn)動(dòng)的速率得到了記錄，其隨機(jī)過程在于隨機(jī)抽取分子進(jìn)行碰撞，而后的分子碰撞過程中，碰撞前，分子速度方向在二維實(shí)驗(yàn)室空間的指向是各向同性的，碰撞后，各個(gè)分子速度方向在質(zhì)心坐標(biāo)系中的指向是隨機(jī)的，其速度大小由碰撞后的具體指向決定，最后記錄碰撞后的速率.本文發(fā)展了文獻(xiàn)［1］的方法，認(rèn)為通過記錄分子運(yùn)動(dòng)的速度，只需隨機(jī)碰撞后各個(gè)分子相對(duì)質(zhì)心坐標(biāo)系的速度方向，使其在三維空間內(nèi)的指向是各向同性的，即可生成三維麥克斯韋分布.文獻(xiàn)［2］中是用蒙特卡羅方法生成三維麥克斯韋分布，其隨機(jī)抽樣的思想同文獻(xiàn)［1］，但碰撞隨機(jī)過程不同，即在球坐標(biāo)系內(nèi)隨機(jī)碰撞前時(shí)兩粒子之間的速度夾角，而后交換連心線方向的速度.由于其在球坐標(biāo)系下的夾角是用均勻分布的隨機(jī)數(shù)代替的，如果是這樣，那么碰撞前各個(gè)分子的速度方向在實(shí)驗(yàn)室坐標(biāo)系中的空間指向?qū)⒉皇歉飨蛲缘?因此其正確性有待商榷.另外其做出的分布曲線并未與理論曲線進(jìn)行對(duì)比，而本文中會(huì)對(duì)數(shù)值計(jì)算結(jié)果與理論預(yù)期進(jìn)行詳細(xì)的對(duì)比.本文僅考慮經(jīng)典條件下的孤立系統(tǒng)的三維麥克斯韋分布，對(duì)于相對(duì)論情形下的麥克斯韋分布見文獻(xiàn)［3］.另外考慮到理想氣體的定義，文中在描述碰撞時(shí)，粒子將被看成是質(zhì)點(diǎn).用鋼球模型模擬二維氣體的麥克斯韋分布的見文獻(xiàn)［4］.其中，對(duì)大量粒子服從的統(tǒng)計(jì)規(guī)律性，有詳細(xì)闡述.也有對(duì)多粒子系統(tǒng)和少粒子系統(tǒng)的粒子速率和能量的分布的模擬和討論，但沒有對(duì)模擬算法的詳細(xì)說明，本文會(huì)有較大篇幅涉及對(duì)算法的說明，說明時(shí)主要與文獻(xiàn)［1］中的算法做對(duì)比，并最終證明本文與文獻(xiàn)［1］算法的等價(jià)性.模擬結(jié)果會(huì)跟對(duì)應(yīng)參數(shù)下的理論曲線進(jìn)行對(duì)比，以說明其正確性（理論曲線非調(diào)參調(diào)出來(lái)〃強(qiáng)行擬合成功的”，而是只設(shè)置了分子質(zhì)量和對(duì)應(yīng)模擬時(shí)初始化的溫度代入密度分布函數(shù)中，作圖后與統(tǒng)計(jì)直方圖進(jìn)行定性擬合）.至于實(shí)驗(yàn)驗(yàn)證，如果拿文中的用于數(shù)值模擬的模型來(lái)做實(shí)物實(shí)驗(yàn)可能難以實(shí)現(xiàn).因?yàn)榭赡芎茈y在一個(gè)箱子中統(tǒng)計(jì)兩個(gè)相互碰撞的氣體分子相對(duì)質(zhì)心散射方向.本文運(yùn)用隨機(jī)模型得到麥克斯韋分布，可以作為大量粒子碰撞形成麥克斯韋分布時(shí)所服從的統(tǒng)計(jì)規(guī)律的依據(jù)：碰撞后兩理想氣體分子相對(duì)質(zhì)心系的速度散射方向的統(tǒng)計(jì)分布是各向同性的.1基本原理1.1理想氣體模型根據(jù)熱學(xué)教材中的準(zhǔn)理想氣體的定義可知，對(duì)理想氣體分子有：分子近似看成質(zhì)點(diǎn)；氣體分子的大小忽略不計(jì).1.2分子碰撞的隨機(jī)過程把所有分子都看成是質(zhì)點(diǎn)后，在碰撞過程中決定他們運(yùn)動(dòng)末了狀態(tài)的參量的因素是:兩個(gè)粒子的初始速度，另一個(gè)是粒子碰撞后相對(duì)質(zhì)心系的速度的空間取向.由于不是運(yùn)用分子動(dòng)力學(xué)模擬的方法記錄粒子的位置和速度，而是只記錄粒子的速度，因此碰撞前后，兩粒子的連心線取向完全不知情只能假定它們碰后質(zhì)心速度在空間中指向分布各向同性（后文驗(yàn)證了這個(gè)假設(shè)能生成麥克斯韋分布）.物理圖象見圖1和圖2.總結(jié)粒子碰撞模型的出發(fā)點(diǎn)如下.首先碰撞滿足動(dòng)量能量守恒；其次，碰撞前后粒子相對(duì)質(zhì)心系的速率大小不變，但方向改變，其方向的空間指向具有隨機(jī)特性，并且其方向的隨機(jī)矢量的指向分布在空間中是各向同性的.圖1未記錄粒子位置，碰撞后兩粒子遠(yuǎn)離質(zhì)心系原點(diǎn)的方向不確定.隨機(jī)從可能的方向中選取一個(gè)，且任何方向被選中的概率均等圖2取兩邊各500個(gè)分子進(jìn)行碰撞（質(zhì)心速度為零），碰后分子遠(yuǎn)離質(zhì)心速度方向分布為各向同性.2前后速度關(guān)系式的推導(dǎo)以及能量和動(dòng)量守恒的證明2.1速度關(guān)系式基于這兩個(gè)出發(fā)點(diǎn)，這里推導(dǎo)碰撞前后的速度關(guān)系式.設(shè)碰撞前，把兩個(gè)粒子分別編號(hào)1、2.1和2的速度分別是：和(1)質(zhì)心運(yùn)動(dòng)速度是：(2)1和2相對(duì)于質(zhì)心的速度是：⑶v2cold=-v1cold⑷現(xiàn)在將這個(gè)矢量在質(zhì)心系中的方向進(jìn)行重新的分配：(程序中可以以原點(diǎn)為中心先生成一個(gè)立方體內(nèi)均勻分布的點(diǎn)，而后在把在這個(gè)立方體內(nèi)處于其內(nèi)切球內(nèi)的點(diǎn)選取出來(lái)，再把這些點(diǎn)作為矢量的失端，最后把這些矢量歸一化)設(shè)這個(gè)隨機(jī)單位矢量為⑸重新分配后1相對(duì)質(zhì)心的速度的表達(dá)式：(6)根據(jù)動(dòng)量守恒2相對(duì)質(zhì)心的速度表達(dá)式是:v2cnew=-v1cnew⑺最終1和2碰撞后的速度表達(dá)式是：=vc+v1cnew(8)=vc+v2cnew⑼2.2能量守恒動(dòng)量守恒的證明能量守恒的證明聯(lián)立式(1)一式(9)得到1和2速度的分量表達(dá)式:(10)(11)(12)(13)所以得證能量守恒.動(dòng)量守恒的證明(14)Pt=m+m=（15）3對(duì)隨機(jī)過程的說明首先有必要對(duì)數(shù)值模擬過程進(jìn)行說明.蒙特卡羅方法生成二維麥克斯韋分布的隨機(jī)過程如下：［1］1） 首先初始化粒子的速率；2） 碰撞時(shí)，隨機(jī)選取兩個(gè)粒子；3） 隨機(jī)選取兩粒子碰撞前的速度方向夾角，計(jì)算質(zhì)心速度；4） 隨機(jī)選取相對(duì)于質(zhì)心系的散射方向，碰撞，計(jì)算碰后的速率……最后生成分布.本文的模型與文獻(xiàn)［1，2］-致，都是質(zhì)點(diǎn)模型，所不同的是本文通過記錄粒子碰撞的速度矢量，從而省去文獻(xiàn)［1］中的第3）步.下面給出本文的具體數(shù)值模擬算法.1） 初始化所有粒子速度;生成各向同性隨機(jī)分布的空間單位矢量，并排列成矩陣.2） 碰撞時(shí)隨機(jī)選取兩個(gè)粒子；3） 計(jì)算兩個(gè)粒子的相對(duì)質(zhì)心的速度；4） 順序選取其中生成各向同性隨機(jī)分布的空間單位矢量，其方向用來(lái)初始化選中的兩個(gè)粒子相對(duì)質(zhì)心的速度方向，而后計(jì)算兩個(gè)粒子新的相對(duì)實(shí)驗(yàn)室坐標(biāo)系的速度;如此迭代下去，達(dá)到熱平衡的時(shí)候，也就是經(jīng)過足夠多次的蒙特卡羅計(jì)算，即可生成三維麥克斯韋分布.下面對(duì)文獻(xiàn)［1］和本文方法等價(jià)性進(jìn)行說明.本文的方法經(jīng)過足夠多次的迭代后（達(dá)到熱平衡時(shí)），粒子的速度矢量的方向在空間的分布是各向同性的，隨機(jī)選取兩個(gè)粒子i和「在選擇之前，它們的速度矢量的方向完全是隨機(jī)的（圖3）.那么從這個(gè)意義上說，它們之間的夾角也無(wú)從確認(rèn).但是有一點(diǎn)可以確認(rèn)，那就是:取定一個(gè)粒子時(shí)，它有一個(gè)速度方向，再取一個(gè)粒子時(shí)，后者的速度方向矢量與前者的相對(duì)指向完全是隨機(jī)的（因?yàn)樗粤Ｗ拥乃俣仁噶渴强臻g中均勻分布的）.這便是文獻(xiàn)［1］的蒙特卡羅方法生成二維麥克斯韋分布編程中步驟3）的思想.故迭代夠多次后，本文模型與文獻(xiàn)［1］中模型等價(jià).且把文獻(xiàn)［1］的模型推廣到了三維.圖3迭代后速度方向單位矢量的分布圖圖4單位球面上均勻分布的單位矢量的分布圖（用于與圖3進(jìn)行對(duì)比）4擬合檢驗(yàn)4.1擬合原理速度分布和速率分布的密度分布函數(shù)分別是：而正態(tài)分布的密度函數(shù)是所以麥克斯韋分布的各個(gè)方向速度分布與M=0的正態(tài)分布有相同的函數(shù)表達(dá)式，其中，（16）用normfit對(duì)x、y和z軸方向速度的計(jì)算機(jī)模擬數(shù)據(jù)進(jìn)行正態(tài)分布N（p,a2）擬合.擬合時(shí)理論上口是0,但是由于漲落，口在擬合時(shí)會(huì)略大于eps（計(jì)算時(shí)能分辨的最小精度），這不與理論矛盾.把擬合出來(lái)的3個(gè)方向上的。做平均，得最終的另外，在數(shù)值模擬初始化過程中，對(duì)所有粒子有一個(gè)統(tǒng)計(jì)上的方均速率，這也對(duì)應(yīng)這個(gè)系統(tǒng)的溫度.由于這是一個(gè)孤立系統(tǒng)，且系統(tǒng)能量守恒，方均速率不會(huì)變，那么溫度也不會(huì)變.把初始溫度，作為理論預(yù)期溫度（理論上，速率和速度概率密度分布函數(shù)中的T）.上述分布函數(shù)中m為0.029kg/N0是空氣平均分子質(zhì)量，kB為玻爾茲曼常量.把參數(shù)和kB代入密度分布函數(shù)，做出擬合曲線，與統(tǒng)計(jì)分布做對(duì)比，即得擬合結(jié)果4.2擬合結(jié)果圖5三維麥克斯韋分布數(shù)值結(jié)果與同溫度下的理論預(yù)期分布作對(duì)比（柱狀圖為數(shù)值結(jié)果，曲線為理論預(yù)期）圖6三維麥克斯韋分布中x方向速度分布的數(shù)值結(jié)果與同溫度下的理論預(yù)期分布作對(duì)比（柱狀圖為數(shù)值結(jié)果，曲線為理論預(yù)期）圖7三維麥克斯韋分布中y方向速度分布的數(shù)值結(jié)果與同溫度下的理論預(yù)期分布作對(duì)比（柱狀圖為數(shù)值結(jié)果，曲線為理論預(yù)期）圖8三維麥克斯韋分布中z方向速度分布的數(shù)值結(jié)果與同溫度下的理論預(yù)期分布作對(duì)比（柱狀圖為數(shù)值結(jié)果，曲線為理論預(yù)期）結(jié)果擬合度較好！這說明數(shù)值模擬結(jié)果從x、y和z方向的速度分布以及最終的三維的速率分布都符合理論預(yù)期.另外，我們還把初始溫度與迭代后的溫度進(jìn)行了對(duì)比：（17）某次初始化時(shí)為T〃初溫”=248.2580K（打引號(hào)因?yàn)榇藭r(shí)系統(tǒng)未達(dá)熱力學(xué)平衡，溫度直接代表了此孤立系統(tǒng)的溫度）；經(jīng)過模擬以后通過擬合3個(gè)方向的平均速度分布得出的T末溫=248.2652K;運(yùn)行十次程序得到的〃初溫”和末溫的平均相對(duì)差距是n=4.2181x10-5.所以，生成的三維麥克斯韋分布的丁是符合能量守恒定律的，這也是微觀上兩體碰撞能量守恒的必要條件.到此，證明了生成的麥克斯韋分布與理論曲線的相似度極高，并再次證明了其迭代過程是符合能量守恒的.5平均速率檢驗(yàn)麥克斯韋分布達(dá)到熱平衡時(shí)的平均速率(18)而平均速率在初始化時(shí)和足夠迭代次數(shù)以后會(huì)有很大的不同，原因是系統(tǒng)尚未達(dá)到熱平衡.但驗(yàn)證平均速率是否隨著迭代次數(shù)改變也是有必要的.驗(yàn)證：某次運(yùn)行程序，用初溫計(jì)算出的理論模擬后直接對(duì)速率平均后得到的某次運(yùn)行10次程序得到的和相對(duì)差值的平均值為q1=9.2412x10-4.且也不隨迭代次數(shù)改變.運(yùn)行10次程序，得到10個(gè)的平均值和最初的相對(duì)差值的平均值為r|2=1.2501x10-4.另外0.9212與的比值(反映了平均速率與總能量的關(guān)系)滿足某次運(yùn)行程序10次過程中，直接平均得到的與用得到的平均比值為0.9212.總之，平均速率與總能量的關(guān)系在達(dá)到熱平衡后在誤差允許的范圍內(nèi)符合理論預(yù)期.平均速率在達(dá)到熱平衡后其值在誤差允許的范圍內(nèi)保持不變.所以此粒子隨機(jī)碰撞模型的數(shù)值模擬結(jié)果與統(tǒng)計(jì)理論的預(yù)期符合較好，由此提出本文的結(jié)論:形成麥克斯韋分布的微觀粒子碰撞散射可能遵循的統(tǒng)計(jì)規(guī)律：碰撞后兩理想氣體分子相對(duì)質(zhì)心系的速度散射方向在空間中的統(tǒng)計(jì)分布是各向同性的.6總結(jié)與展望基于基本的理想氣體模型，以及能量動(dòng)量守恒定律，還有質(zhì)心系質(zhì)點(diǎn)速度散射的空間取向各向同性隨機(jī)分布，導(dǎo)出了三維蒙特卡羅方法生成麥克斯韋分布的模型.用MATLAB平臺(tái)生成了三維麥克斯韋分布，并進(jìn)行了擬合檢驗(yàn)，還把溫度和平均速度這兩個(gè)物理量進(jìn)行了理論預(yù)期檢驗(yàn).此模型在檢驗(yàn)中表現(xiàn)良好.提出了形成麥克斯韋分布的微觀粒子的碰撞散射的可能遵循的統(tǒng)計(jì)規(guī)律：碰撞后兩理想氣體分子相對(duì)質(zhì)心系的速度散射方向在空間中的統(tǒng)計(jì)分布是各向同性的.望對(duì)揭示麥克斯韋分布的物理原理有所幫助.致謝：感謝趙虎副教授和彭芳麟教授對(duì)本文發(fā)表給予的討論和支持.【相關(guān)文獻(xiàn)】彭芳麟.計(jì)算物