湯文輝徐彬彬 冉憲文 徐志宏

(國防科學技術大學理學院,長沙 410073)(2016年10月17日收到;2016年12月17日收到修改稿)

專題:高壓下物質(zhì)的新結(jié)構(gòu)與新性質(zhì)研究進展

高溫等離子體的狀態(tài)方程及其熱力學性質(zhì)

湯文輝?徐彬彬 冉憲文 徐志宏

(國防科學技術大學理學院,長沙 410073)(2016年10月17日收到;2016年12月17日收到修改稿)

高溫下等離子體的狀態(tài)方程及其熱力學性質(zhì)在天體物理、可控核聚變以及武器設計與破壞效應等領域有著廣泛應用.本文主要回顧了高溫等離子體在不同狀態(tài)區(qū)域的狀態(tài)方程的理論模型和處理方法.對于理想等離子體,離子之間的相互作用可以忽略,其狀態(tài)方程較簡單,已趨于完善.在超高溫下,原子完全電離,離子和電子都可以采用理想氣體狀態(tài)方程描述;當溫度不太高時,離子部分電離,可以采用Saha方程及其修正模型描述;原子在高度壓縮狀態(tài)下,其狀態(tài)方程可以采用Thomas-Fermi模型及其改進模型得到.對于非理想等離子體,離子之間存在強耦合,還沒有單一的理論模型能夠在任意密度和溫度范圍內(nèi)對離子之間的相互作用進行統(tǒng)一描述.量子分子動力學方法原則上可以在較大溫度密度范圍內(nèi)給出可靠結(jié)果,但由于計算量太大以及高溫下的計算存在收斂問題,也較難應用到溫度較高的稠密等離子體區(qū)域.半經(jīng)驗的經(jīng)典分子動力學方法雖然簡單、計算量小,但只能在一定的區(qū)域范圍內(nèi)給出較精確的狀態(tài)方程結(jié)果.在不同溫度密度區(qū)域內(nèi)采用不同的計算模型,再在空白區(qū)域進行插值從而得到全局狀態(tài)方程在目前不失為一種簡單有效的方法.

高溫等離子體,狀態(tài)方程,熱力學性質(zhì)

1 引 言

自然界中大約90%—95%的可觀測物質(zhì)處于極端高溫高壓條件下[1].太平洋底部(海平面下11 km處)壓強為0.12 GPa;地球中心處的壓強為360 GPa,溫度為0.5 eV(1 eV=11604.5 K),密度約為10—20 g/cm3;木星中心壓強約為4—6 TPa,溫度約為2 eV,密度約為30 g/cm3;太陽中心壓強為2.4×107GPa,溫度為1.6×103eV,密度約為150 g/cm3;白矮星壓強為1012—1018GPa,密度為106—109g/cm3,溫度為103eV.在慣性約束聚變(ICF)中,中心點火需要氘氚(DT)熱斑溫度達到5 keV以上,密度100 g/cm3以上,壓強達到1.5×107—2.0×107GPa.目前已知的自然界的物態(tài)范圍大致如圖1[1]所示(1 bar=105Pa).然而,在實驗室條件下只能產(chǎn)生相對較低的壓強和溫度.靜高壓加載能達到550 GPa的壓強[2],結(jié)合激光加熱技術可使實驗溫度能夠達到3000 K以上[3].動高壓加載雖然能產(chǎn)生104GPa的瞬態(tài)高壓,但持續(xù)時間短,實驗測量非常困難,測量結(jié)果誤差較大[4].事實上,人們真正了解的物質(zhì)狀態(tài)范圍只是圖1中很小的一個區(qū)域,研究更大溫度密度范圍的物質(zhì)狀態(tài)需要發(fā)展更先進的實驗手段和理論方法.

狀態(tài)方程是描述物質(zhì)系統(tǒng)中各狀態(tài)變量之間關系的函數(shù)表達式,反映物質(zhì)在一定熱力學條件下的宏觀性質(zhì)[1,4].近代自然科學和工程技術中的大量實際問題,比如天體物理、恒星以及行星內(nèi)部構(gòu)造、慣性約束聚變、武器系統(tǒng)的設計及其破壞效應等問題,都涉及到物質(zhì)處于高溫、高壓、高密度等極端條件,要對這些問題進行研究,就需要了解物質(zhì)在各種極端條件下的狀態(tài)方程及其熱力學性質(zhì)[5].

眾所周知,高溫條件下的物質(zhì)以等離子體的形式存在.隨著科學技術的發(fā)展以及國防科技的驅(qū)動,等離子體的熱力學性質(zhì)及其描述方法已成為目前的研究熱點之一.然而,等離子體的熱力學性質(zhì)及其描述方法與其所處壓強和溫度狀態(tài)密切相關,難度很大.為了便于人們對等離子體的熱力學性質(zhì)開展深入研究,本文在文獻報道的基礎上對高溫(1 eV以上)等離子體的狀態(tài)方程模型及其熱力學性質(zhì)進行了綜述.

圖1 自然界和實驗室中的物態(tài)范圍[1],圖中括號里的數(shù)字表示密度的對數(shù)值,密度單位為g/cm3;“Statics”區(qū)域?qū)o高壓實驗所能達到的溫度壓強范圍,“Dynamics”區(qū)域表示動高壓加載實驗所能達到的溫度壓強范圍Fig.1.Extreme states in natu re and in the laboratory[1],the numbers in parentheses ind icate the logarithmof density(g/cm3).The “statics” domain corresponds tostatic techniques for high-pressure production,the“dynamics” domain tothe dynamic ones.

2 等離子體狀態(tài)區(qū)域的劃分

等離子體可以根據(jù)離子間相互作用的強弱劃分為理想等離子體和非理想等離子體[6].定義ΓD為離子之間相互作用能與離子動能的比值,ΓD可用于評估離子之間耦合的強弱,稱為離子耦合系數(shù).ΓD?1表示離子之間的相互作用非常小,這時可以認為等離子體是理想等離子體;ΓD～1表示離子之間的相互作用不容忽略,等離子體是非理想等離子體;ΓD?1時,離子間存在較強的相互作用,離子間的相互作用能占主導地位.

另一方面,根據(jù)粒子量子效應的強弱可以將等離子體分為強簡并等離子體和弱簡并等離子體,粒子簡并度參數(shù)可用α表示如下:

式中,n是粒子數(shù)密度,λˉ是電子德布羅意波長,?是約化Planck常數(shù),m是粒子質(zhì)量,k是Boltzmann常數(shù),T是溫度.α=n3?1表示粒子間距遠大于粒子的德布羅意波長,粒子處于非簡并或者弱簡并態(tài),量子效應不明顯,這時粒子可以作為經(jīng)典粒子處理,遵循經(jīng)典的Boltzmann統(tǒng)計.當粒子間距較小時,粒子處于強簡并態(tài)(α?1),服從Fermi量子統(tǒng)計.相對于電子,離子的質(zhì)量較大,德布羅意波長非常小,因此離子的量子效應只有在極低溫情況下才需要考慮.對于質(zhì)量最小的氫離子,在固體密度狀態(tài)下,溫度需要降低到幾十開爾文(K),氫離子的量子效應才需要考慮.對于慣性約束聚變中心點火熱斑,氘氚燃料被壓縮到原來固體密度的1000倍以上,溫度到達5 keV以上[7],由于氘氚離子滿足α?1的條件,因此中心氘氚熱斑不需要考慮氘氚離子的量子效應,但對于電子,其量子效應不能忽略.

為了方便,我們根據(jù)離子耦合系數(shù)以及電子簡并度參數(shù)的大小,把高溫等離子體分成四個區(qū)域進行討論,如圖2所示.I區(qū)和II區(qū)分界線是原子完全電離所需溫度,II區(qū)和IV區(qū)是根據(jù)離子相互作用耦合強弱劃分的,I區(qū)和III區(qū)是根據(jù)電子簡并強弱劃分的.由于不同元素完全電離溫度不一致,圖2中的區(qū)域劃分主要適用于低Z元素,對于高Z元素,由于原子完全電離所需溫度增加,因此I區(qū)和II區(qū)的分界線需要適當上移.I區(qū)是完全電離的等離子體,在溫度非常高的情況下,原子中所有電子電離,這時粒子之間的相互作用可以忽略不計,電子也處于非簡并狀態(tài)(ΓD?1,α?1),滿足經(jīng)典Boltzmann統(tǒng)計,此區(qū)域的等離子體中電子和原子核都可以采用理想氣體狀態(tài)方程描述.II區(qū)是部分電離弱耦合等離子體,溫度不是特別高,部分電子處于自由態(tài),部分處于束縛態(tài).此區(qū)域密度較低,電子和電子之間以及電子和離子之間的相互作用較弱(ΓD? 1,α? 1),但在靠近 ΓD=1的區(qū)域,需要考慮粒子間的遠程Coulomb相互作用.III區(qū)是強簡并等離子體,物質(zhì)密度非常高,原子處于高度壓縮狀態(tài),電子分布在原子核周圍,對原子核產(chǎn)生屏蔽,原子核與原子核的相互作用可以忽略,電子處于強簡并狀態(tài)(ΓD?1,α?1).IV區(qū)是強耦合等離子體,物質(zhì)密度較高,溫度相對較低,離子之間耦合作用非常強,電子處于部分電離部分簡并狀態(tài)(ΓD?1或者ΓD～1,α～1).I區(qū)、II區(qū)以及III區(qū)的等離子體由于離子間相互作用可以忽略不計,因此可以當作理想等離子體處理,現(xiàn)有的理論模型能夠較好地對這些區(qū)域的等離子體狀態(tài)方程進行描述,而IV區(qū)是強耦合等離子體,離子之間的相互作用還沒有一個單一的理論模型能夠在任意密度和溫度范圍內(nèi)進行統(tǒng)一的描述.

原子序數(shù)為Z的離子間Coulomb作用能為

式中

其中,M為原子摩爾質(zhì)量,ρ為密度,NA為Avogadro常數(shù),e為電子電荷量,r是粒子間的平均距離,r～ n?1/3.

單個粒子的平均動能為

對于I區(qū)和II區(qū)的稀薄理想等離子體,ΓD?1.如果離子之間主要是Cou lomb相互作用,則耦合系數(shù)滿足條件

經(jīng)過變換,得到滿足稀薄理想等離子體的條件為

式中溫度單位為K.由(6)式可知,氣體在常態(tài)密度情況下,溫度只需要幾個eV就可以滿足稀薄氣體條件;對于鋰、鈹?shù)容^輕的元素,常態(tài)密度下,溫度需要達到幾十eV才能滿足稀薄氣體條件;對于重元素,則需要上百eV的溫度才能滿足稀薄氣體條件.對于完全電離的稀薄等離子體,電子和離子都滿足經(jīng)典Boltzmann條件,這時可以把電子和離子都作為理想氣體處理.對于溫度較低的稀薄理想等離子體,其狀態(tài)方程和電離度相關,而電離狀態(tài)可以采用Saha方程[8]進行描述.雖然帶電粒子間的長程Coulomb作用能相對粒子動能較小,但Cou lomb相互作用引起附加的自由能,會導致電離能減小,因此在密度稍高或者溫度較低的情況下,也就是ΓD?1條件滿足較弱的情況下,需要考慮粒子間長程Coulomb作用的影響,這時的狀態(tài)方程可以采用Debye-Hückel模型[9]計算.

圖2 溫度密度平面上等離子體區(qū)域的劃分,其中,I區(qū)為完全電離等離子體,II區(qū)為部分電離區(qū)弱耦合等離子體,III區(qū)為強簡并等離子體,IV區(qū)為強耦合等離子體;圖中粗實線是電子簡并強弱分界線,即α=1,粗虛線表示離子耦合強弱分界線,即ΓD=1,細實線表示四個區(qū)域的分界線,Sun以及ICF分別表示太陽和慣性約束聚變條件下氫的同位素所處的大致區(qū)域Fig.2.D ivision of plasma on the density and temperatu re region:Iis the fu lly ionized plasma,IIis partly ionized and week coupling plasma,IIIis the strong degenerate plasma and IV is the strong coupling plasma.The thick solid lines are the boundary ofweak and strong degenerate of electron(α=1)and the thick dashed are the boundary ofweak and strong coupling of ion(ΓD=1).The thin solid lines are the boundary of four regions.Sun and ICF mean the regions of hyd rogen’s isotopes under the condition of sun and inertial con finement fusion.

III區(qū)為強簡并等離子體,是原子被高度壓縮所形成的.原子的內(nèi)層電子被均勻擠壓,電子處于強簡并態(tài),滿足關系,電子的分布可以使用Thomas-Fermi(TF)模型[10,11]進行描述.TF模型認為原子中電子連續(xù)分布在原子核周圍[12,13],原子核被連續(xù)的電子云所包圍,電子是一團簡并氣體,并遵循Fermi量子統(tǒng)計.由于分布在原子核周圍高密度電子的屏蔽作用,原子核與原子核之間的相互作用可以忽略.當原子中電子處于強簡并狀態(tài)時,費米能為Ek～2n2/3/(2m),于是ΓD正比于n?1/3.粒子數(shù)密度n越大時,ΓD越小,因此電子強簡并條件下的等離子體在高度壓縮條件下也趨于理想等離子體.

對于IV區(qū)的強耦合等離子體,離子之間存在較強的相互作用(ΓD～1或者ΓD?1).由于電子處于部分簡并狀態(tài),這時單純采用經(jīng)典Boltzmann統(tǒng)計模型或者Fermi量子統(tǒng)計模型都不能正確地描述電子的統(tǒng)計狀態(tài).物質(zhì)處于稠密狀態(tài)時,相鄰原子之間的距離非常小,離子間存在較強的相互作用,因此必須正確描述離子間的相互作用才能準確描述物質(zhì)的熱力學性質(zhì).分子動力學方法可以較好地描述粒子之間的相互作用,但前提是需要知道粒子之間的相互作用勢.根據(jù)所使用的原子間相互作用勢的不同,分子動力學方法可以分為兩大類:一類是經(jīng)典分子動力學方法(MD),粒子之間相互作用采用經(jīng)驗勢或半經(jīng)驗勢進行描述;而另一類是量子分子動力學方法(QMD),原子間相互作用勢基于量子力學理論計算得到,比如基于密度泛函理論計算的分子動力學方法[14].兩種分子動力學方法各有優(yōu)缺點.經(jīng)典分子動力學方法由于直接使用經(jīng)驗勢,因此計算方法簡單,計算量小;但由于不同物質(zhì)的相互作用勢差別較大,一種經(jīng)驗勢一般只適用于一種原子系統(tǒng)或一部分系統(tǒng),而且經(jīng)驗勢也局限于一定的密度溫度范圍,因此使用某種經(jīng)驗勢的經(jīng)典分子動力學方法只適用于特定溫度密度范圍內(nèi)的某些物質(zhì).量子分子動力學方法不用假設粒子間的相互作用勢,而是由第一性原理計算原子相互作用勢,所以具有普適性,但由于在高溫時計算量較大,而且難以收斂,因此主要適用于溫度不太高的狀態(tài).在計算電子結(jié)構(gòu)時為了節(jié)約機時,量子分子動力學常采用贗勢的方法,因此在溫度較高時,電子躍遷和電子激發(fā)的處理也存在困難,這時需要采用一些半經(jīng)驗的分子動力學或改進的量子分子動力學模型來描述稠密物質(zhì)的狀態(tài).總之,等離子體處于溫稠密和熱稠密狀態(tài)時,難以采用單一的理論模型對其狀態(tài)方程進行精確描述,一般需要采用多種理論模型相結(jié)合,各自描述一部分密度和溫度范圍,然后通過插值得到過渡段的狀態(tài)方程.

3 高溫稀薄等離子體

對于理想氣體,不考慮原子之間的相互作用力,其狀態(tài)方程可表示為

式中P為壓強;ε為比內(nèi)能;n為粒子數(shù)密度;γ為氣體比熱比,與氣體的內(nèi)部自由度相關,如果氣體有q個自由度,則γ=1+2/q.對于單原子分子q=3,γ=5/3;對于雙原子分子q=5,γ=7/5.

高溫條件下,分子先離解成原子,隨后原子電離成各級離子.分子的離解以及原子的電離都會導致氣體的粒子數(shù)密度發(fā)生明顯變化.當溫度升高到原子完全電離(圖2中I區(qū)),等離子體的粒子(包括離子與電子)數(shù)密度由原來的n變?yōu)?1+Z)n,Z為原子的核電荷數(shù).因此,原子完全電離的高溫稀薄等離子體狀態(tài)方程可以表示為

對于部分電離的稀薄等離子體(圖2中II區(qū)),其狀態(tài)方程可采用Saha方程[8]描述.溫度不是很高時,原子并未完全電離,帶電粒子的平均動能遠遠大于粒子之間的Coulomb相互作用能,這時Saha方程可以很好地描述其電離行為.Saha方程適用的范圍是ΓD?1,對于原子序數(shù)為Z的粒子,稀薄氣體條件即是(6)式.在低密度或者高溫情況下都可以滿足Saha方程條件,但在密度較高或者溫度較低情況下(圖2中II區(qū)靠近IV區(qū)位置),電離氣體系統(tǒng)中帶電粒子間的長程Coulomb作用不可忽視.Coulomb相互作用引起附加的自由能,使電離能減小.對于這樣的系統(tǒng),必須對Saha模型添加Debye-Hückel近似進行修正.

3.1 稀薄電離氣體的Saha模型

大部分分子在1 eV以上時已經(jīng)離解為原子,有些原子已經(jīng)發(fā)生了電離,隨著溫度的進一步升高,較重原子所組成的氣體會多次電離.考察某種原子所組成的氣體,單位質(zhì)量氣體中含有N個原子,單個原子的總能量為E=Ek+Q,其中Ek為單個原子的動能,包括離子以及電子的動能;Q為原子內(nèi)電子的勢能.在氣體足夠稀薄的情況下,電子遵循Boltzmann統(tǒng)計規(guī)律,這時氣體中每一個粒子(包括中性原子、離子以及電子)的平均熱能是(3/2)kT,因此N個原子系統(tǒng)總的熱能為

其中αe為平均電離度.

用Im表示第m次激發(fā)電子所必須消耗的能量,則從原子中激發(fā)m個電子需要消耗的總能量為

在給定溫度T以及密度ρ的情況下,單位質(zhì)量氣體中含有N0個中性原子,N1個一次電離的離子,Nm個電離m次的離子,以及Ne個電子.此外電離m次的離子還具有電子激發(fā)能Hm,因此氣體總的內(nèi)能為

(11)式中,αe是氣體的平均電離度,即Ne/N;αm是電離m次離子的濃度,表示為Nm/N.由于氣體比較稀薄,在粒子之間相互作用可以忽略的條件下,所有粒子都可以當作理想氣體,因此可以得到電離氣體的壓強為

氣體的壓強與總的粒子數(shù)密度相關,總的粒子數(shù)密度等于電子數(shù)密度與離子數(shù)密度之和.等離子體電離度越大,其壓力越高.電離度可以通過原子數(shù)守恒、電荷守恒以及Saha方程聯(lián)立求得

3.2 考慮輻射超高溫等離子氣體

當?shù)入x子體氣體溫度非常高或者密度非常低時,這時輻射壓以及輻射能相對于氣體的壓強以及能量并不能被忽略,需要把輻射的能量以及壓強加入到氣體的能量和壓強中.平衡輻射條件下,輻射的比能量等于輻射能量密度除以物質(zhì)的密度,即

輻射壓強可以表示為

輻射壓和溫度的四次方成正比,而根據(jù)(12)式,物質(zhì)壓和溫度的一次方成正比,因此當溫度超過一定范圍之后,輻射壓將完全占據(jù)主導地位.令輻射壓等于物質(zhì)壓,可以得到兩種壓強平衡時溫度和密度之間的關系式為

式中密度單位取g/cm3.對標準密度的空氣來說,當空氣溫度超過100 eV時,輻射壓起主要作用;對于密度為1 g/cm3的固體來說,輻射壓要在溫度達到幾keV以上時才起主要作用.

如果輻射能流和物質(zhì)能流為同一量級,這時輻射輸運在能量輸運過程中必須考慮.輻射能量以光速進行傳播,而物質(zhì)能量以物質(zhì)的速度進行傳播,在輻射溫度不是很高的情況下,也需要對輻射能流進行考慮.同樣令物質(zhì)能流等于輻射能流,假設物質(zhì)速度為光速的千分之一,則在固體密度條件下,當溫度在幾十eV時需要考慮輻射輸運過程,這時不能單純使用流體動力學過程描述物質(zhì)的運動,而需要使用輻射流體動力學模型來考察物質(zhì)的運動.

3.3 沖擊壓縮下的等離子體氣體

在單次沖擊壓縮條件下,當沖擊波強度足夠大時,氣體的極限壓縮度可以表示為

對于單原子理想氣體,γ=5/3,h=4,即氣體密度最多只能被壓縮4倍;對于雙原子理想氣體,γ=7/5,h=6.可見,氣體比熱比γ越接近于1,其壓縮度也就越大.

在沖擊波壓縮下,波后物質(zhì)的壓強、密度與比內(nèi)能之間的關系由Hugoniot方程給出:

等離子體的內(nèi)能可以分為兩部分,即E=Ek+Q,其中Ek為粒子的平均動能,Q包含了粒子的勢能以及粒子內(nèi)部自由度的能量.由于壓強只和粒子的平動動能相關,即P=(2/3)ρEk,代入(18)式,在強沖擊波條件下,忽略沖擊波波前壓力P0和內(nèi)能E0,將ρ/ρ0用h代替,得到

從(19)式可以看出,勢能Q相對于熱運動動能忽略不計(3Q/Ek?1)時,單原子氣體的壓縮度最大值為4.在電離情況下,總的粒子數(shù)目增加,但(19)式表明,壓縮度和總的粒子數(shù)密度沒有關系,因此總粒子數(shù)目的增加并不會增加壓縮度.另一方面,電離度增加會導致用于電離的能量損耗增加,也就是勢能對總內(nèi)能的相對貢獻增加,因此原子在沖擊壓縮下發(fā)生電離時,會使氣體變得更加容易壓縮.隨著原子最外層電子開始電離,到最外層電子全部電離,這個過程用于電離的能量所占總內(nèi)能的比例增加,因此這個過程壓縮度是增加的.在原子最外層電子(K層電子)全部電離而次外層電子(M層電子)還未開始電離時,隨著沖擊波強度的增加,在相當一段范圍內(nèi)離子不再發(fā)生電離(因為次外層電子電離能比最外層電子電離能大很多),勢能Q保持不變,而溫度升高成為內(nèi)能增加的主要形式,因此Q/Ek會減小,從而導致壓縮度降低.如果沖擊波強度足夠大,能夠?qū)е麓瓮鈱与娮与婋x,這時隨著沖擊波強度的進一步增加,次外層電子會像最外層電子一樣一個個剝落,用于電離的勢能相對貢獻又增加,因此導致壓縮度再一次增大.如果是重原子氣體,原子含有多個殼層,則隨著各層電子的依次電離,沖擊壓縮條件下的壓縮度先增大再減小,然后又增大再減小.如果沖擊壓縮過后原子發(fā)生完全電離,隨著沖擊波強度進一步增加,會導致溫度急劇上升,粒子的熱運動動能急劇增加,而勢能Q保持不變,因此壓縮度最后會減小.在沖擊波強度非常大的情況下,如果不考慮熱輻射,根據(jù)(19)式,完全電離氣體的壓縮度也是趨近于4.

3.4 考慮粒子間Cou lomb作用的稀薄等離子體

對于滿足理想等離子體條件較弱的稀薄等離子體,粒子之間的Coulomb相互作用相對于熱運動能量雖然較小,但不能忽視.這時Coulomb相互作用可以采用Debye-Hückel模型進行計算.Debye-Hückel模型假定所有帶電粒子的勢場都是球?qū)ΨQ的,帶電粒子之間只考慮Coulomb作用力,但仍然假定Coulomb作用能仍遠小于粒子的熱運動動能,因此有

或者表示為

式中?(r)表示r處的電勢,ni(r)是r處第i種粒子的數(shù)密度,Zi是第i種粒子的電荷數(shù).第i種帶電粒子的勢能可以表示為

由于粒子滿足Boltzmann統(tǒng)計規(guī)律(稀薄電離氣體條件下,電子和離子都可以看作經(jīng)典氣體),粒子數(shù)密度的空間分布為

式中ni0是電荷為Zi的粒子的平均數(shù)密度.(23)式后面約等號成立的條件是粒子Coulomb作用能相對于粒子的平均動能是小量.由于滿足電荷守恒,把(23)式代入泊松方程(21)得到

在粒子中心附近有r?λD,因此有

式中第一項為中心粒子本身所產(chǎn)生的勢,第二項?′=Zie/λD是周圍其他粒子在r處所形成的勢.

2015年11月底，《中共中央國務院關于進一步推進農(nóng)墾改革發(fā)展的意見》中要求，用3年左右時間，將國有農(nóng)場承擔的社會管理和公共服務職能納入地方政府統(tǒng)一管理，基本完成農(nóng)墾國有土地使用權(quán)確權(quán)登記發(fā)證任務。

在體積V中,等離子體的Coulomb相互作用能為所有粒子相互作用勢的一半,利用(22)式有

異性帶電粒子存在吸引力作用,每個離子周圍都有大量的電子,所以Coulomb作用能以及壓強都是負值.Cou lomb作用影響氣體的狀態(tài)主要體現(xiàn)在兩個方面,一個是減小氣體的內(nèi)能和壓強,還有一個效應是降低電離能,從而使得電離更加容易發(fā)生.Coulomb相互作用中的電子并不是完全自由的,由于電子還在離子的勢場中,其勢能并不為0,而是負值,因此把電子從原子或者離子中電離出來需要消耗的電離能會下降.電離能的減小可以表示為

在較低溫度下,?In相對于kT可以忽略不計,但在上百eV高溫下,?In/In可以達到10%.在密度較高情況下,考慮粒子之間的遠程Coulomb作用,壓強和內(nèi)能會減小,也會使得電子電離能下降,從而使電離更容易發(fā)生.

4 強簡并等離子體

前面所討論的電離氣體,電子和離子都遵循經(jīng)典的Boltzmann統(tǒng)計,但只有在溫度相當高或者密度相當?shù)偷那樾蜗?才能滿足經(jīng)典Boltzmann統(tǒng)計條件.物質(zhì)在高度壓縮條件下,原子內(nèi)電子高度簡并,這時必須考慮電子的量子效應.對于物質(zhì)在高溫高壓條件下的物態(tài)性質(zhì),比較簡潔方便的處理方法是采用TF模型[15].TF模型描述的是電子對壓強和內(nèi)能的貢獻,原子核對壓強以及內(nèi)能的貢獻可以采用理想氣體狀態(tài)方程描述,

式中,na為原子核數(shù)密度.

電子對內(nèi)能的貢獻包括電子的動能、電子與其他電子以及原子核之間的相互作用能.下面主要討論電子對系統(tǒng)壓強以及內(nèi)能的貢獻.

4.1 費米簡并態(tài)自由電子氣體

自由電子在低密度或者溫度較高情況下服從經(jīng)典的Boltzmann統(tǒng)計,但在溫度較低或者密度較高的情況下,由于電子處于簡并狀態(tài),其量子態(tài)必須要考慮.零溫時,自由電子處于完全簡并狀態(tài),在dVdp相空間的量子數(shù)為4πp2dpdV/3,由于電子有兩個不同的自旋方向,因此在相空間內(nèi)量子數(shù)還需要乘以2.根據(jù)泡利不相容原理,每一個確定的自旋方向的量子狀態(tài)中只能有一個電子,假設體積V中N個電子緊湊地占據(jù)了動量從0到p0的量子態(tài),因此有關系式

定義費米溫度TF=ε0/k為

電子的簡并程度也可以根據(jù)費米溫度進行判斷,如果電子溫度T?TF,則電子處于弱簡并態(tài)或者非簡并態(tài),電子可以作為經(jīng)典氣體處理;如果T?TF,則電子處于強簡并狀態(tài),服從量子分布.比較(1)和(32)式有α∝(TF/T)3/2,可見α?1與T?TF是等價關系,都表示電子處于弱簡并狀態(tài);而α?1與T?TF也是等價關系,表明電子處于強簡并狀態(tài).因此,(1)式和(32)式都可以評判電子的簡并程度,并且評判結(jié)果是一致的.

體積V中N個電子的動能為

根據(jù)熱力學關系Td S=d E+P d V,零溫下完全簡并電子氣體的壓強為

從(34)式可以看出,在極低溫條件下,自由簡并電子的壓強正比于密度的5/3次方,而在超高溫情況下,電子服從經(jīng)典Boltzmann統(tǒng)計,這時可以采用理想氣體狀態(tài)方程描述,壓強與密度的一次方成正比,這說明電子氣體處于簡并狀態(tài)時更難壓縮.高溫壓縮下原子中的電子并不是完全自由的,而是在原子核以及其他電子所形成的勢場中運動,這時電子的行為可以使用TF模型進行描述.

4.2 TF模型

單個原子中電子在原子核周圍服從Fermi-Dirac分布,分布函數(shù)為

式中,μ是電子的化學勢,E是電子的能量.在TF模型中,電子在原子核以及其他電子所產(chǎn)生的勢場中運動,其能量主要包括動能以及在原子內(nèi)的靜電勢能,E=p2/(2me)?eU(r),U(r)是r處的電勢,p為電子的動量.利用電子的分布函數(shù)可以得到電子數(shù)密度為

原子內(nèi)電子總數(shù)等于核電荷數(shù),即

(37)式稱為電荷守恒條件,由此可得到化學勢μ.

在電子數(shù)密度已知條件下,原子內(nèi)的電勢可以通過泊松方程得到

(38)式就是TF方程.TF方程是二階微分方程,對其求解需要一定的邊界條件.在原子中心附近,也就是半徑趨向于0時,由于電子均勻連續(xù)分布在原子核周圍,因此原子中心處的電勢全部由原子核提供,這時的邊界條件可以表示為

TF方程再加上兩個邊界條件可以封閉求解,從而得到TF模型下原子中的電子分布.由于電子的分布影響勢場,而勢場又會影響電子的分布,勢場和電子分布耦合在一起,因此TF方程的解析解難以得到,一般采用數(shù)值計算或者近似方法求解.雖然TF方程沒有解析解,但TF方程的解具有普適性質(zhì),只要溫度和體積以TZ?4/3,VZ形式出現(xiàn),則壓強、內(nèi)能、比熱容以PZ?10/3,EZ?7/3,CVZ?1的形式出現(xiàn),這意味著只要知道某種元素的TF狀態(tài)方程,根據(jù)上面所得到的普適性質(zhì),可以得到任意元素的TF狀態(tài)方程.對于多種原子所組成的物質(zhì),可以利用平均原子的概念得到其TF狀態(tài)方程.平均原子序數(shù)和平均原子量可以表示為

式中xi是第i種原子數(shù)目所占總原子數(shù)目的百分比,Zi和Mi分別為第i種原子的原子序數(shù)和原子量.

TF模型不區(qū)分電子是否處于束縛狀態(tài),因此這也意味著物質(zhì)的熱力學性質(zhì)和原子的殼層結(jié)構(gòu)沒有任何關系,只和原子序數(shù)相關.從TF模型的普適性也可以看出,在相同的密度條件下,隨著溫度升高,物質(zhì)的壓強和內(nèi)能是隨著原子序數(shù)的增加而單調(diào)增加的.由于不區(qū)分電子的殼層結(jié)構(gòu),溫度的升高只是改變原子中電子密度的重新分布,因此電子勢能的變化較緩慢,電子密度在原子體積中的分布更趨于均勻.在溫度非常高的情況下,也就是電子溫度遠遠大于費米溫度(T?TF),這時電子的簡并性將被消除,等離子體從圖2中的III區(qū)過渡到I區(qū),這時的狀態(tài)方程和完全電離理想氣體狀態(tài)方程一樣,壓強和密度成正比.

4.3 改進的TF模型

簡單的TF模型只考慮了電子之間的Cou lomb相互作用,而沒有考慮交換相互作用.電子是費米子,服從泡利不相容原理,其所引起的自旋平行電子之間的相互作用在計算電子在原子核周圍的分布時需要考慮.Dirac[16]于1930年引入交換能對零溫TF模型進行了修正,由此得到了Thomas-Fermi-Dirac(TFD)模型.Cowan和Ashkin[17]進一步將TFD模型推廣到有限溫度情況.TF模型對電子的能量只考慮了電子的動能、電子與原子核以及其他電子相互作用勢能,TFD模型還考慮了電子與所有自旋平行電子的交換能.考慮交換能后電子分布函數(shù)可以表示為

式中Eex是自由電子的交換能.由Bloch[18]估算得到的交換能表示為

如果Eex～kT或者Eex? kT,則電子交換能對TF模型的修正非常明顯,而在Eex? kT時,電子交換能可以忽略,這時溫度需要滿足的條件為T? 1.37×10?15n4/3.當物質(zhì)處于固體密度時,溫度需要幾十甚至上百eV,電子交換項的影響才可以忽略不計.與TF模型結(jié)果進行比較,TFD模型的壓強在低溫下(幾十eV以下)比TF模型要低,隨著溫度的升高,兩者差別越來越小,最后會重合在一起.McCarthy[19]采用不同模型進行計算給出了鉛在不同溫度下壓強隨密度的變化,如圖3所示[19].可以看出,零溫(kT=0 eV)時,采用TFD模型得到的壓強明顯要比TF模型得到的壓強低,10 eV時兩者差別減小,在100 eV時已基本重合在一起.

圖3 鉛在100,10和0 eV下的壓強曲線:TF結(jié)果、TFD結(jié)果以及TFK結(jié)果比較[19]Fig.3.Pressu re profi le of lead at the temperatu re of 100,10 and 0 eV,the solid,dotted and dashed lines are the resu lts of TF,TFD and TFKmodels,respectively[19].

TFD模型雖然引入了交換修正,但并沒有考慮電子的量子效應.1957年,Kirzhnits[20]將電子分布函數(shù)展開,在Thomas-Fermi模型中同時引入交換修正和量子修正,形成了Thomas-Fermi-Kirzhnits(TFK)理論.TFK模型引入的修正項和TF模型一樣具有普適性,對于TF模型有

TFK模型對TF模型的修正項和TF模型的普適變化關系有所不同,其修正項可以表示為[21]

McCarthy[19]詳細地列出了較大溫度(TZ?4/3)和體積范圍內(nèi)(VZ)壓強和能量的計算結(jié)果.在計算精度要求不是很高的情況下,TFK方法的解析擬合表達式[12,22]可以為高能量密度條件下的流體模擬提供簡單的狀態(tài)方程.

從圖3可以看出,在100 eV的高溫時,三種模型(TF,TFD,TFK)基本重合在一起;10 eV時,TF模型結(jié)果稍微偏高,但在零溫且密度較低時,TF模型壓強明顯偏大,而TFD模型和TFK模型由于考慮了交換修正以及量子效應修正,比TF模型要精確.

TF模型把所有電子作為準自由電子處理,但實際上靠近原子核的電子受到原子核強烈的Coulomb作用,仍處于束縛態(tài),具有殼層結(jié)構(gòu).靠近原子核內(nèi)層電子處于束縛狀態(tài),遵從量子力學規(guī)律,而較外層的電子由于熱電離或者壓致電離較容易受激發(fā)成為自由電子,因此將原子內(nèi)的電子分成束縛電子和自由電子兩部分來處理更接近實際情況.考慮殼層結(jié)構(gòu)修正的TFS(Thomas-Fermi shell)方法[23]就是基于電子分成兩部分的設想而改進的TF模型.

將原子體積分為兩個區(qū)域:內(nèi)層靠近原子核是電子密度較高的束縛區(qū)域;外層是密度較低的自由電子區(qū)域.外層可以采用TF模型方法處理,而內(nèi)層束縛態(tài)可以采用Wenzel-Kramers-Brillouin近似處理.TFL方法把束縛電子看成準自由電子,雖然考慮了內(nèi)層束縛電子但并不求解薛定諤方程.Zink[23]分別使用TFL(Thomas-Fermi like)方法和TFS方法計算了鐵的壓強,如圖4所示.殼層修正在低溫低密度條件下起作用,在高溫或者高密度條件下,兩者差不多完全重合.在低溫低密度時,TFS結(jié)果反映了電子的凝聚現(xiàn)象,即?p/?ρ較小,因為一些束縛態(tài)能級處于空缺狀態(tài),外層自由電子較容易重新回到束縛態(tài).當密度增加到一定程度時,束縛態(tài)能級充滿電子,隨著壓縮繼續(xù)進行,將會出現(xiàn)壓致電離,這時壓強會突然升高.于是,這種電子凝聚和壓致電離交替現(xiàn)象發(fā)生在TFS結(jié)果中.Lee和Thorsos[24]對TFS方法進行了改進,認為有些電子可能處于自由態(tài)和束縛態(tài)之間的共振狀態(tài),將電子分成自由電子以及束縛電子和共振態(tài)電子三項,提出了TFSR方法.他們對Z=5到Z=30的各種元素進行了計算,結(jié)果如圖5[24]所示.可以看出,壓強隨Z呈現(xiàn)周期性變化,這和TF模型壓強隨原子序數(shù)Z單調(diào)上升有明顯區(qū)別.

圖4 鐵在7.7和770 eV下,不同密度時的壓強,圖中分別給出了TF模型、TFL模型以及TFS模型結(jié)果[23]Fig.4.Pressure profi le of iron with diff erent density at the temperature of 7.7 and 770 eV,the solid,dashed and dashed-dotted lines are the resu lts of TF,TFL and TFSmodels,respectively[23].

平均原子(average atom,AA)模型[25?27]也可用于高溫下計算等離子體中的電子結(jié)構(gòu).AA模型不區(qū)分離子的不同組態(tài),而是把等離子體中所有離子等效成一個平均原子.在這個原子中,電子按照Fermi-Dirac統(tǒng)計規(guī)律排布在各個單電子軌道上,忽略了電子之間不同的耦合方式所帶來的差異.AA模型電子軌道類似于氫原子軌道,軌道上的電子分布呈現(xiàn)出殼層結(jié)構(gòu).這是AA模型和TF模型的主要區(qū)別,也是其在低溫低密度下具有更高精度的原因.AA模型發(fā)展至今已有許多改進,比如不可分辨的躍遷系模型[28,29],超級躍遷系模型[30,31],細致譜模型[32]等,這些模型對電子殼層結(jié)構(gòu)考慮得越精細,所計算的狀態(tài)方程越精確,但所需要的計算量也越大.

圖5 不同溫度及電子數(shù)密度下TF模型和TFSR模型壓強和原子序數(shù)的關系[24]Fig.5.The relations of pressure and atomic number at d iff erent density and diff erent temperature,the dotted and solid lines are the resultsof TFSR and TFmodels,respectively[24].

5 強耦合等離子體

前面所介紹的模型都是基于單原子假設,即離子之間的相互作用能相對于離子的動能可以忽略不計,物質(zhì)的性質(zhì)可簡化為求單個原子的性質(zhì).當物質(zhì)處于密度較高,而溫度并不是特別高的溫稠密或者熱稠密狀態(tài)時,滿足ΓD～1或者ΓD?1,α～1,即離子間的耦合較強,電子處于部分電離部分簡并態(tài),這時使用單原子的性質(zhì)來推導物質(zhì)的宏觀性質(zhì)將會帶來較大誤差.

把強耦合等離子體分為電子和離子兩部分進行描述.對于電子部分的描述,可以采用TF模型、AA模型、無軌道(orbital-free,OF)模型以及量子力學理論.離子間相互作用可采用超網(wǎng)鏈(hypernetted-chain,HNC)模型[33?35]、蒙特卡羅(MC)方法以及MD方法進行描述.結(jié)合電子的描述模型以及離子相互作用模型,涌現(xiàn)了一大批計算高溫稠密物質(zhì)物態(tài)性質(zhì)的模型方法:Ofer等[36]采用TF模型來描述電子結(jié)構(gòu),離子間相互作用采用HNC近似來處理.Furukawa和Nishihara[37]采用AA模型并考慮量子修正的 HNC近似處理離子相互作用.Zérah等[38?41]采用TF模型計算電子密度空間分布,采用MD來描述離子間相互作用,提出了TFMD模型.Dharma-wardana和Francois[42]在密度泛函框架下通過Kohn-Sham方程處理電子問題,然后采用HNC近似來描述高溫稠密等離子體;Dharma-wardana和Murillo[43]采用雙溫模型解決高溫稠密物質(zhì)中的離子間相互作用問題.Zérah等[44,45]結(jié)合無軌道密度泛函理論和MD方法,發(fā)展了無軌道分子動力學方法(OFMD).侯永等[46,47]結(jié)合考慮能級展寬效應的AA模型和MD方法提出了AAMD模型.Car和Parrinello[14]結(jié)合密度泛函理論和MD方法,提出了第一性原理分子動力學方法(CPMD).Dai等[48?50]基于量子分子動力學方法,考慮了高溫下電子碰撞離子對離子位置所形成的漲落,發(fā)展了能描述從低溫凝聚態(tài)到高溫理想等離子體的量子Langevin分子動力學方法(QLMD).前面幾種基于TF模型、OF模型、AA模型的半經(jīng)驗方法,適用于特定的溫度以及密度條件,而基于第一性原理的QMD,以及基于路徑積分的蒙特卡羅(PIMC)方法[51]被認為是比較精確的模型.但它們所需計算量非常巨大,尤其是PIMC方法,計算量隨著電子個數(shù)急劇增加,其應用受限于電子數(shù)很少的低Z原子.由于計算量大,QMD方法所能給出的狀態(tài)方程數(shù)據(jù)點較少,所以一般被用來作為評判其他半經(jīng)驗方法是否適用的基準[52?54].

隨著計算方法和超級計算機的發(fā)展,分子動力學方法得到了迅速發(fā)展,分子動力學方法可以有效地描述離子之間的相互作用,但前提是需要離子之間相互作用勢.下面分別介紹使用經(jīng)驗或者半經(jīng)驗勢的經(jīng)典分子動力學方法和利用量子力學原理直接計算相互作用勢的量子分子動力學方法.

5.1 經(jīng)典分子動力學方法

MD方法是一種用來計算經(jīng)典多體體系平衡和傳遞性質(zhì)的方法,其目的是模擬系統(tǒng)的微觀細致動力學行為.在分子動力學方法中,忽略原子核運動的量子效應,原子核的運動滿足經(jīng)典力學的牛頓定律.當體系中所有粒子上的作用力已知時,便可以通過牛頓運動方程得到系統(tǒng)中所有粒子的運動軌跡,再通過統(tǒng)計平均得到該系統(tǒng)的宏觀熱力學參數(shù).由于粒子在勢場中運動,其受力等于勢函數(shù)的空間導數(shù),因此,從經(jīng)典力學出發(fā),求解原子間聯(lián)立牛頓運動方程的分子動力學方法主要包括兩個問題:一是粒子相互作用勢函數(shù)的選取,二是積分牛頓運動方程的求解.

對于一個含有N個粒子的系統(tǒng),每個粒子的坐標、速度、動量以及受力可以分別表示為ri(t),vi(t),pi(t),f(ri,t),系統(tǒng)的哈密頓量可以表示為

式中第一項是系統(tǒng)動能,第二項U(rN)為系統(tǒng)勢能函數(shù),rN=(r1,r2,...,rN)代表N個原子的坐標.哈密頓方程可以表示為

(48)式說明,離子運動所受到的力由離子所處的勢場所決定.已知離子所處的勢場,通過(48)式加上初始條件以及周期性邊界條件,可以進行數(shù)值求解,從而得到離子的運動軌跡,進而得到各種宏觀物理量.根據(jù)統(tǒng)計物理思想,系統(tǒng)經(jīng)過足夠長的時間之后會處于平衡狀態(tài),從而表征系統(tǒng)宏觀性質(zhì)的物理量A可以通過時間平均得到

在時間足夠長、系統(tǒng)達到平衡態(tài)后,有

分子動力學方法的實質(zhì)是利用計算機對系統(tǒng)中所有粒子的牛頓運動方程進行數(shù)值求解,得到粒子坐標、速度等動力學量隨時間變化的函數(shù),然后通過統(tǒng)計平均求得表征系統(tǒng)宏觀性質(zhì)的物理量.因此分子動力學計算過程可分為下面三個步驟:首先將牛頓方程進行離散化,求解各個粒子在t+?t的坐標;接下來是計算不同時刻(t+ ?t,t+2?t,...,t+l?t)的某個物理量A(r1(t+?t),r2(t+?t),...,rN(t+?t));最后對各個物理量求平均值:

經(jīng)典分子動力學采用的是經(jīng)驗勢或者半經(jīng)驗勢,相互作用勢不需要另外求解,因此計算速度較快,能夠計算的粒子數(shù)量也較大,可以達到千萬個原子量級.常用的經(jīng)驗勢函數(shù)有Lennard-Jones勢函數(shù)[55],Morse勢函數(shù)[56],Born和Mayer[57]提出的描述離子晶體中離子之間相互作用的Born-Mayer勢.至今為止,已經(jīng)有各種各樣的勢函數(shù)形式被提出,這些勢函數(shù)只考慮兩個原子之間的相互作用,僅和兩個原子位置坐標相關,在許多無機化合物中可以有效使用,但對于金屬以及共價鍵結(jié)合的有機分子,不能給出精確的結(jié)果.在20世紀80年代初,各種考慮了多體相互作用的新勢函數(shù)被相繼提出.Daw和Baskes[58]基于密度泛函理論提出了嵌入原子模型(EAM);Finnis和Sinclair[59]根據(jù)密度函數(shù)二次矩理論提出了形式上與EAM基本相同的經(jīng)驗F-S模型.Abell[60]在1985年提出了一種Morse形式的多原子間相互作用勢模型,比傳統(tǒng)的二體勢和三體勢模型更加精確.

傳統(tǒng)勢函數(shù)對處于常溫常壓下的物質(zhì)可以給出很好的描述,而對處于高溫等離子體狀態(tài)時的物質(zhì)則存在較大誤差.一些描述電子結(jié)構(gòu)的半經(jīng)驗方法:比如TF模型、AA模型、OF模型等能夠在高溫下較精確地描述離子運動的勢場.結(jié)合這些半經(jīng)驗方法和分子動力學方法便可得到描述高溫稠密物質(zhì)性質(zhì)的模型:比如TFMD模型[38,41]、AAMD模型[46,47]以及OFMD模型[44,45].但由于TF模型、AA模型以及OF模型的假設都是基于高溫的情況,因此這些模型只能在高溫下描述電子結(jié)構(gòu),在低溫情況下會帶來一定的誤差.

5.2 第一性原理分子動力學方法

第一性原理分子動力學方法是基于密度泛函理論計算原子間相互作用勢的一種量子分子動力學方法.1985年,Car和Parrinello[14]把密度泛函理論應用到分子動力學中,提出了第一原理分子動力學方法(ab initiomolecular dynamics).該方法極大地推動了計算機模擬低溫(T≤10 eV)下物態(tài)性質(zhì)的發(fā)展,已成為計算機模擬最先進和最重要的方法之一.由于不需要對離子間勢函數(shù)做任何經(jīng)驗假設,而是直接采用密度泛函理論求解離子間的相互作用勢,因此該方法具有非常強的普適性,并得到了非常廣泛的應用[61?67].

在第一性原理分子動力學中,原子核間的相互作用勢U(r)由密度泛函理論進行計算,整個分子動力學過程可分為三步:第一步,對于給定的原子核坐標ri,求解Kohn-Sham方程[68]得到粒子數(shù)密度分布,進而得到相互作用勢;第二步,計算每個原子核受到的力;第三步,求解牛頓運動方程

密度泛函理論基于絕熱近似,在求解薛定諤方程時,電子和原子核的坐標同時出現(xiàn)在電子和原子核相互作用項中,直接求解難度非常大.由于原子核質(zhì)量遠大于電子質(zhì)量,原子核的運動速度相對電子的運動速度非常小.在考慮電子的運動時,可以認為電子絕熱于原子核的運動;在考慮原子核運動時,電子可以及時調(diào)整空間分布來適應原子核的運動,這就是絕熱近似或稱Born-Oppenheimer近似[69].通過絕熱近似,求解多粒子體系薛定諤方程可以簡化為分別求解電子運動和原子核運動的問題.密度泛函理論的基本思想是將原子、分子和固體的基態(tài)物理性質(zhì)用粒子密度函數(shù)來描述.密度泛函理論的基礎是Hohenberg和Kohn關于非均勻電子氣體理論的Hohenberg-Kohn定理[70].粒子數(shù)密度函數(shù)是確定多粒子系統(tǒng)基態(tài)物理性質(zhì)的基本變量,要確定系統(tǒng)基態(tài)物理性質(zhì)需要三個確定量:粒子數(shù)密度函數(shù)、動能泛函、交換關聯(lián)能泛函.粒子數(shù)密度函數(shù)和動能泛函可以采用Kohn-Sham方程進行求解;而交換關聯(lián)能泛函一般采用局域密度近似[71]和廣義梯度近似[72]得到.

原則上,量子分子動力學方法可以描述任意條件下物質(zhì)的物態(tài)性質(zhì),但第一性原理分子動力學的應用范圍一直局限于溫度較低的區(qū)域(10eV以下),主要有以下原因:第一是在高溫條件下需要非常大的計算量;第二是大部分第一性原理程序為了減少計算量都是基于贗勢實現(xiàn)的,而在高溫高壓下對贗勢的標準和構(gòu)造要求較高;第三是高溫狀態(tài)難以達到熱平衡,計算過程中難以保證收斂性[73].當溫度較高時,許多電子成為非局域電子或自由電子,相對于離子運動的時間尺度,電子運動的時間尺度非常小,因此在每一分子動力學時間步中電子與離子發(fā)生了多次碰撞.在溫度升高時,自由電子增多,電子的運動速度增快,與離子碰撞次數(shù)也增多,這時碰撞引起離子位置的漲落是不可忽略的,Born-Oppenheimer近似也不再成立.傳統(tǒng)的第一性原理分子動力學方法基于Born-Oppenheimer近似,并沒有考慮離子位置的漲落,因此在高溫下的計算不可避免地帶來一些問題.Dai等[48?50]對電子與離子在高溫下的碰撞作用采用類似布朗運動的思想來處理,引入離子-電子碰撞項,利用Langevin方程描述離子運動,取得了很好的效果,并將第一性原理分子動力學方法拓展到了理想等離子體、高能量密度物理區(qū)域,填補了低溫凝聚態(tài)與高溫稠密狀態(tài)之間的空白.

第一性原理分子動力學方法和其改進方法雖然能夠在很大密度以及溫度范圍內(nèi)精確地得到物質(zhì)的狀態(tài)參數(shù),但由于其計算的復雜性以及所耗計算機時巨大,在精確度要求不是特別高的流體模擬中不如一些經(jīng)驗或者半經(jīng)驗的模型簡單方便.一些全局狀態(tài)方程就是在不同溫度密度區(qū)域結(jié)合各類實驗結(jié)果、理論結(jié)果再采用插值方法獲得的,比如More等[74]所提出的QEOS(quotidian equation of state),以及美國洛斯阿拉莫斯國家實驗室的SESAME數(shù)據(jù)庫[75].

5.3 全局狀態(tài)方程

高能量密度條件下的流體動力學過程涉及非常寬泛的物態(tài)范圍,比如慣性約束聚變中所涉及的溫度范圍為10?3—104eV、密度范圍為10?3—102g/cm3,因此在流體動力學模擬中需要全局狀態(tài)方程.全局狀態(tài)方程就是在非常大的密度和溫度范圍內(nèi)描述壓強P(ρ,T)和內(nèi)能E(ρ,T).盡管量子分子動力學方法原則上可以對任意狀態(tài)的物質(zhì)提供統(tǒng)一的描述,但實際上,目前還沒有哪個模型能夠覆蓋所有區(qū)域和所有物質(zhì),從而不能滿足流體動力學模擬的需求.在不同區(qū)域采用不同的近似模型,并在這些區(qū)域間進行光滑的插值,從而構(gòu)筑全局狀態(tài)方程模型是目前所采用的主要方法.目前已有一些大范圍的狀態(tài)方程(EOS)模型,比如Bushman和Fortov[76],Godval和Sikka[77],Basko[78]所提出的各種模型以及More等[74]所提出的QEOS模型.還有一些是以表格形式給出的,比如SESAME數(shù)據(jù)庫[75].SESAME數(shù)據(jù)庫的狀態(tài)方程數(shù)據(jù)是在不同溫度密度區(qū)域結(jié)合各類實驗、理論結(jié)果采用插值方法獲得的.SESAME數(shù)據(jù)庫和QEOS模型使用比較廣泛,它們也是很多文獻進行狀態(tài)方程數(shù)據(jù)分析和比對的重要參考.

QEOS把自由能分為電子項、離子項以及半經(jīng)驗的束縛修正項,表示如下:

對于電子項的貢獻Fe,可采用TF模型進行計算;對于離子項的貢獻Fi,在固體中把熱離子當作聲子氣體處理,采用玻色-愛因斯坦統(tǒng)計進行計算,流體中則采用Cowan模型處理.半經(jīng)驗的束縛項自由能不依賴于溫度,可表示為

式中E0和b是兩個經(jīng)驗參數(shù),它們滿足以下兩個關系式:第一個是在常溫和固體密度條件下,總壓強為0;第二個是需要知道物質(zhì)的體積模量以及物質(zhì)的初始密度,通過壓強和自由能的熱力學關系式再聯(lián)立(51)式可以得到壓強和兩個經(jīng)驗參數(shù)的關系,再通過B=ρ0(?P/?ρ)ρ0在已知體積模量B以及初始密度ρ0的條件下可以得到另一個關于E0和b的關系式.利用這兩個關系式可求得兩個經(jīng)驗參數(shù)E0和b的值.

QEOS模型在國際上得到了比較廣泛的使用,很多輻射流體力學程序使用QEOS模型的計算結(jié)果作為狀態(tài)方程參數(shù),QEOS有可以免費獲得的計算程序[79].Atzeni和Meyer-ter-Vehn[80]將鋁的Hugoniot曲線的理論計算結(jié)果與實驗結(jié)果進行了比較,如圖6所示.可以看出,QEOS結(jié)果和SESAME結(jié)果在104GPa以下與實驗所得到的數(shù)據(jù)[81?84]比較接近.在104GPa以上時,兩者之間有一定的差別,但總體上和實驗數(shù)據(jù)符合較好,這說明QEOS和SESAME都能給出較精確的狀態(tài)方程參數(shù).

圖6 鋁的Hugoniot曲線[80],圖中給出了QEOS結(jié)果[74],SESAME結(jié)果[75]以及一些實驗結(jié)果[81?84]Fig.6.Hugoniot curves of aluminum[80],the b lack solid line and red dashed line are the resu lts of QEOS and SESAME,respectively.The dots in figure are experimental resu lts[81?84].

圖7[85]給出了溫度為5和30 eV時不同密度下Al等離子體的狀態(tài)方程.圖中比較了TF模型[86],QEOS模型[74],SESAME數(shù)據(jù)庫[75],NPA模型[87]以及AAMD模型[88]結(jié)果.NPA(neutral pseudoatom)模型[87]采用密度泛函理論處理電子結(jié)構(gòu),采用HNC近似描述離子間相互作用.從圖7可以看出,30 eV時所有模型結(jié)果的符合度比5 eV時好.溫度較高時,一些基于高溫條件的半經(jīng)驗模型,比如TF模型以及AAMD模型,也可以得到比較精確的結(jié)果.從圖7還可以看出,在5 eV時TF模型以及AAMD模型所得到的結(jié)果在較高密度條件下明顯要比NPA模型結(jié)果、QMD結(jié)果以及SESAME數(shù)據(jù)庫結(jié)果高.

圖7 溫度為5和30 eV時不同密度下鋁的狀態(tài)方程[85],圖中分別是TF模型[86]、QEOS[74]、SESAME[75]、NPA模型[87]以及AAMD模型[88]的計算結(jié)果比較Fig.7.The EOS of aluminiumunder the temperatu re of 5 and 30 eV[85],the solid lineswith d iff erentmarks are the resu lts of TF[86],QEOS[74],SESAME[75],NPA[87]and AAMD[88]models.

圖8[89]給出了TFMD模型[90]、 OFMD模型[73]、AAMD模型[88]、QMD模型、QLMD模型[50]以及SESAME數(shù)據(jù)庫[75]在不同溫度下鐵的Hugoniot曲線.由于高溫下QMD模型不容易收斂,因此沒有給出高溫下QMD的結(jié)果.從圖8可以看出,在高溫段(100 eV以上),TFMD,AAMD,OFMD以及QLMD四個模型和SESAME數(shù)據(jù)庫符合得非常好,這說明四個MD模型在高溫下都能較精確地給出鐵的狀態(tài)方程,但在低溫下,TFMD模型和OFMD模型結(jié)果明顯要比其他幾個模型要高.相對前兩個模型,AAMD模型能在更低的溫度下給出精確結(jié)果.這三種半經(jīng)驗模型結(jié)果在低溫下比量子分子動力學結(jié)果要高,這說明了低溫下使用半經(jīng)驗模型的局限性.溫度較低時,電子交換、關聯(lián)作用需要考慮[91],而TFMD方法在計算電子壓強時沒有考慮交換、關聯(lián)作用所產(chǎn)生的負壓強對總壓強的貢獻,所以TFMD在低溫區(qū)給出的壓強偏高.AA模型雖然考慮了電子的殼層結(jié)構(gòu),但也只是簡單的結(jié)構(gòu)近似,因此低溫時AAMD的精確度也會下降,這時只有基于密度泛函理論的量子分子動力學模型才能給出較精確的結(jié)果.在0.1 eV的低溫下,電子的自旋效應不能被忽視,考慮了量子自旋效應的QLMD模型結(jié)果比未考慮量子自旋效應的結(jié)果要稍微高一些,與SESAME數(shù)據(jù)庫結(jié)果更接近.

圖8 不同溫度下鐵的Hugoniot曲線[89],圖中分別給出了考慮了交換作用的TFMD模型[90]、AAMD模型[88]、OFMD模型[73]、QMD模型以及QLMD[50]結(jié)果和SESAME數(shù)據(jù)庫[75]結(jié)果Fig.8.The Hugoniot cu rve of iron at diff erent temperatu res[89],the resu lts of TFMD[90],AAMD[88],OFMD[73],QMD,QLMD[50]and SESAME[75]models are presented,respectively.

圖9[89]針對鐵的Hugoniot曲線給出了TFD模型[90]、VAAQP模型[92](variational-averageatom-in-quantum-plasmas,基于平均原子模型近似)、QEOS模型、QLMD模型[50]以及SESAME數(shù)據(jù)庫和各種實驗結(jié)果[81,93,94]的比對.由于沖擊壓縮實驗中一些物理量的測量比較困難,實驗結(jié)果有較大的不確定性,因此各類實驗點的數(shù)據(jù)比較分散.從圖9可以看出,不同的理論模型給出的結(jié)果也存在一定差異,尤其壓強在10 Mbar以下、固體密度在2.5倍以下的溫稠密物質(zhì)區(qū)域,差別非常明顯.在此區(qū)域,TFD模型以及VAAQP模型結(jié)果偏大.QEOS模型、QLMD模型以及SESAME數(shù)據(jù)庫結(jié)果與俄羅斯的Rusbank數(shù)據(jù)庫[93]的結(jié)果十分接近.在高壓高密度情況下,各種模型雖然有一定的差異,但總體而言與實驗結(jié)果符合較好.這說明在高溫情況下,各種半經(jīng)驗的統(tǒng)計模型也能給出比較精確的結(jié)果.

圖9 鐵的Hugoniot線的實驗和理論結(jié)果的對比[89]Fig.9.The Hugoniot cu rve of iron[89],the lines are the resu lts of various theoreticalmodels and the dots are the experimental results.

6 結(jié) 語

本文根據(jù)離子耦合系數(shù)以及電子簡并度參數(shù)在溫度密度平面內(nèi)把等離子體狀態(tài)劃分為四個區(qū)域,對四個區(qū)域內(nèi)等離子體的狀態(tài)方程及其熱力學性質(zhì)進行了簡要總結(jié).在完全電離等離子體區(qū)域,溫度非常高,原子完全電離,粒子之間相互作用非常小,電子和離子都可以采用理想氣體狀態(tài)方程描述.在部分電離弱耦合等離子體區(qū)域,部分電子電離,粒子之間相互作用相對較小,狀態(tài)方程可以采用Saha方程求解,考慮粒子之間長程Coulomb相互作用后可采用Debye-Hückel方程得到.在強簡并等離子體區(qū)域,原子核之間的相互作用由于周圍電子的屏蔽可以忽略,原子核可采用理想氣體狀態(tài)方程描述,電子采用TF模型可以得到較精確的結(jié)果,一些改進型的TF模型可以把TF模型拓寬到更低溫度和更低密度.在強耦合等離子體區(qū)域,溫度不是很高,密度在物質(zhì)固體密度附近,此時物質(zhì)處于溫稠或者熱稠等離子體區(qū)域.當?shù)入x子體處于前面三個區(qū)域時,都有較準確的方法描述等離子體的狀態(tài)方程,但等離子體處于溫稠以及熱稠區(qū)域時,還沒有較完善的理論進行統(tǒng)一的描述.量子分子動力學雖然原則上能夠統(tǒng)一描述此區(qū)域物質(zhì)的性質(zhì),但由于計算量大、高溫不收斂等缺點,其使用也受到一定的限制.一些改進的量子分子動力學方法比如QLMD,雖然明顯拓寬了量子分子動力學方法的使用范圍,但其使用條件以及準確性仍需要進一步驗證確認,而且其計算量過大,需要經(jīng)過大量的計算和累計才能得到較密集的狀態(tài)方程插值表,以滿足流體動力學過程模擬的需要.一些半經(jīng)驗分子動力學方法,比如TFMD,AAMD,OFMD等,計算量雖小,但在低溫下的計算結(jié)果誤差較大.可見,發(fā)展描述溫稠和熱稠區(qū)域的先進理論模型十分必要.另一方面，流體動力學模擬所使用的狀態(tài)方程需要兼顧計算量和精度.狀態(tài)方程計算方法過于復雜,計算量太大,滿足不了流體動力學模擬的需要;而狀態(tài)方程參數(shù)精確度太低,也會對流體動力學模擬結(jié)果的可靠性產(chǎn)生較大影響.因此,采用不同的簡單模型在不同密度和溫度區(qū)域進行插值,不失為一種快速獲得等離子體狀態(tài)方程的好方法.

[1]Fortov V E 2016Extreme States of Matter(Berlin:Springer)p9

[2]Xu J A,MaoHK,Bell P M1986Science232 1404

[3]Ming L C,BassettW A1974Rev.Sci.Instrum.45 1115

[4]NellisW J 2006Rep.Prog.Phys.69 1479

[5]Q ian X S 2007Lectures on Physics(Shanghai:Shanghai Jiaotong University Press)(in Chinese)[錢學森2007物理力學講義(上海:上海交通大學出版社)]

[6]Fortov V E 2007Phys.Uspek.50 333

[7]Lind l J D,Amend t P,Berger R L,G lendinning S G,G lenzer S H,Haan SW,Kau ffman R L,Landen OL,Su ter L J 2004Phys.P lasmas11 339

[8]Saha MN,Srivastava BN 1965ATreatise on Heat(Allahabad:The India Press)

[9]Debye P,Hückel E 1923Physik Z24 185

[10]Thomas L H1927Proc.Cambridge Philos.Soc.23 542

[11]Fermi E 1928Z.Phys.48 73

[12]Xu X S,Zhang W X 1986Theory Guide of Practica l Equations of State(Beijing:Science Press)(in Chinese)[徐錫申,張萬箱 1986實用物態(tài)方程理論導引(北京:科學出版社)]

[13]Feynman R P,Metropolis N,Teller E 1949Phys.Rev.75 1561

[14]Car R,ParrinelloM1985Phys.Rev.Lett.55 2471

[15]Letter R 1955Phys.Rev.99 1854

[16]D irac P AM1930Proc.Camb.Phil.Soc.26 376

[17]Cowan R D,Ashkin J 1957Phys.Rev.105 144

[18]Bloch F 1929Zeits.F Phys.57 545

[19]McCarthy S L 1965Lawrence Livermore Laboratory ReportUCRL-14364

[20]Kirzhnits D A1957Sov.Phys.JETP5 64

[21]Tang W H,Zhang R Q 2008In troduction tothe Theory and Compution of Equations of State(Beijing:Higher Education Press)(in Chinese)[湯文輝,張若棋2008物態(tài)方程理論及計算概論(北京:高等教育出版社)]

[22]Bell AR 1980Rutherford and Appleton Laboratories ReportRL-80-091

[23]Zink JW 1968Phys.Rev.176 279

[24]Lee C M,Thorsos BI1978Phys.Rev.A17 2073

[25]Rozsnyai BF 1972Phys.Rev.A5 1137

[26]Rozsnyai BF 1982J.Quan t.Spectrosc.Radiat.Transf.27 211

[27]Rozsnyai BF,Lamou reux M1990J.Quan t.Spectrosc.Radiat.Transf.43 381

[28]Bauche-Arnou lt C,Bauche J,Klapisch 1978M.J.Opt.Soc.Am.68 1136

[29]Bauche-Arnou lt C,Bauche J,Klapisch M1979Phys.Rew.A20 2424

[30]Bar-ShalomA,Oreg J,Goldstein W H,Shvart D,Zigler A1989Phys.Rev.A40 3183

[31]Oreg J,Goldstein W H,Bar-ShalomA,Klapisch M1990J.Comput.Phys.91 460

[32]Iglesias C A,Rogers F J,W ilson BG 1987Astrophys.J.322 45

[33]Kurten KE,Ristig ML,Clark JW 1977Lett.Al NuovoCimen to20 313

[34]Kang HS,Ree F H1998Phys.Rev.E57 5988

[35]Rose D V,GenoniTC,W elch D R,C lark R E,Campbell R B,Meh lhorn TA,Flicker D G 2009Phys.P lasams16 102105

[36]D ror O,Nard i E,Rosen feld Y 1988Phys.Rev.A38 5801

[37]Fu rukawa H,N ishihara K1992Phys.Rev.A46 6596

[38]Zérah G,C lérouin J G,Pollock E L 1992Phys.Rev.Lett.69 446

[39]C lérouin J G,Pollock E L,Zérah G 1992Phys.Rev.A46 5130

[40]Flavien L,Gilles Z 2006Phys.Rev.E73 016403

[41]Danel J F,Kazand jian L,Zérah G 2006Phys.P lasmas13 092701

[42]Dharma-wardana MW C,Francois P 1982Phys.Rev.A26 2096

[43]Dharma-wardana MW C,MurilloMS 2008Phys.Rev.E77 026401

[44]Danel J F,Kazand jian L,Zérah G 2009Phys.Rev.E79 066408

[45]Horner D A,Lambert F,Kress J D,Collins L A2009Phys.Rev.B80 024305

[46]Hou Y,Jin F T,Yuan J M2006Phys.P lasmas13 093301

[47]Hou Y,Jin F T,Yuan J M2007J Phys.:Condens.Matter19 425204

[48]Dai J,Yuan J 2009Europhys.Lett.88 20001

[49]Dai J,Hou Y,Yuan J 2010Phys.Rev.Lett.104 245001

[50]Dai J,Hou Y,Yuan J 2010Astrophys.J.721 1158

[51]Militzer B2009Phys.Rev.B79 155105

[52]W unsch K,Vorberger J,Gericke D O2009Phys.Rev.E79 010201

[53]Gericke D O,W unsch K,G rinenkoA,Vorberger J 2010J.Phys.:Conf.Ser.220 012001

[54]Pey russe O,Mazevet S,Recou les V,Dorchies F,Harmand M,Levy A,Fuchs J,Mancic A,NakatsutsumiM,Renaudin P,Audebert P 2009CP,Atomic Processes in P lasmas1161 200

[55]Frenkel D,Smit B1996Understanding Molecular Simu lation(San D iego:Acedemic Press)

[56]Morse P M1929Phys.Rev.34 57

[57]Born M,Mayer J E 1931Z.Phys.75 1

[58]DawMS,Baskes MI1983Phys.Rev.Lett.50 1285

[59]Finnis MW,Sinclair J E 1984Philosophic Magazine A50 45

[60]Abell G C 1985Phys.Rev.B31 6184

[61]Collins L,Kwon I,Kress J,Trou llier N,Lynch D 1995Phys.Rev.E52 6202

[62]Desjarlais MP,Kress J D,Collins L A2002Phys.Rev.E66 025401

[63]Mazevet S,C lérouin J,Recou les V,Anglade P M,Zérah G 2005Phys.Rev.Lett.95 085002

[64]Mazevet S,Desjarlais MP,Collins L A,Kress D,Magee N H2005Phys.Rev.E71 016409

[65]C lérouin J,Noiret P 2008Phys.Rev.B78 224203

[66]Holst B,Redmer R 2008Phys.Rev.B77 184201

[67]Mazevet S,Zérah G 2008Phys.Rev.Lett.101 155001

[68]Kohn K,ShamL J 1965Phys.Rev.140 A1133

[69]Born M,Huang K1954Dynamical Theory of Crystal Lattice(Ox ford:Ox ford University Press)

[70]Hohenberg P,Kohn W 1964Phys.Rev.B136 864

[71]Hansen J P,McDonald IR 1976Theory of Simple Liquids(Lodon:Academic Press)

[72]PerdewJ P,Kieron B,Matthias E 1996Phys.Rev.Lett.77 3865

[73]Lambert F,Clérouin J,Zérah G 2006Phys.Rev.E73 016403

[74]More R M,W arren KH,Young D A,Zimmerman G B1988Phys.F luids31 3059

[75]Lyon S P,Johnson J D 1992SESAME:The Los Alamos National Laboratory Equation of State DatabaseReport No.LA-UR-92-3407

[76]Bushman AV,Fortov V E 1983Soviet Phys.Uspek.26 465

[77]Godval BK,Sikka S K1983Phys.Rep.102 121

[78]BaskoMM.Metallic 1985Soviet High Temperat.Phys.23 388

[79]KempAJ,Meyer-ter-Vehn J 1998Nucl.Instrum.Meth.Phys.Res.A415 674

[80]Atzeni S,Meyer-ter-Vehn J 2004The Physics of Inertia l Fusion(NewYork:Oxford University Press)

[81]Altshu ler L 1965Soviet Phys.Uspek.8 52

[82]Ragan C E 1982Phys.Rev.A25 3360

[83]V lad imirov AS,Voloshin N P,Nogin V N,Petrov tsev AV,SimonenkoV A1984JETP Lett.39 82

[84]Mitchell AC,Nellis W J 1981J.Appl.Phys.52 3363

[85]Hou Y 2009Ph.D.D issertation(Changsha:National University of Defence Technology)(in Chinese)[侯永2009博士學位論文(長沙:國防科學技術大學)]

[86]Letter R 1955Phys.Rev.99 1854

[87]Perrot F,Dharma-wardana MW C,Benage J 2002Phys.Rew.E65 046414

[88]Hou Y,Yuan J 2009Phys.Rev.E79 016402

[89]Dai J Y 2009Ph.D.D issertation(Changsha:National University of Defence Technology)(in Chinese)[戴佳鈺2010博士學位論文(長沙:國防科學技術大學)]

[90]Lambert F,C lérouin J,Mazevet S 2006Europhys.Lett.75 681

[91]Fromy P,Deu tsh C,Maynard G 1996Phys.P lasmas3 714

[92]Piron R,Blenski T2011Phys.Rev.E83 026403

[93]Shock W ave Database http://www.ihed.ras.ru/rusbank/gassim/[2003-07-13]

[94]Batani D,Morelli A,Tomasini M,Benuzzi-Mounaix A,Philippe F,Koenig M,Marchet B,Masclet I,Rabec M,Reverd in C,Caub le R,Celliers P,Collins G,Silva L D,Hall T,Moret M,Sacchi B,Baclet P,Cathala B2002Phys.Rev.Lett.88 235502

PACS:05.70.Ce,52.25.KnDOI:10.7498/aps.66.030505

Equations of state and thermodynamic properties of hot plasma

Tang Wen-Hui?Xu Bin-Bin Ran Xian-Wen Xu Zhi-Hong

(College of Science,National University of Defense Technology,Changsha 410073,China)(Received 17 October 2016;revised manuscript received 17 December 2016)

The equations of state(EOS)and the thermodynamics properties of plasma under high temperature are widely applied tothe fields of astrophysics,controllable fusion,weapon design and damage.In this paper wemain ly reviewthe theoreticalmodel and computing method of the EOS of hot plasma on diff erent density scales and temperature scales.For an ideal plasma,the interaction between ions can be ignored,the EOS is simple and the theories turn matured.Under the condition of extremely high temperature,ions are ionized completely and the EOSs of ions and electrons can be approximated by the EOS of ideal gas.W hen the temperature is not very high and ions are just partly ionized,the EOS can be obtained by Saha model or its modified model.W hen atoms are strongly compressed,the EOS can be calculated by Thomas-Fermimodel or its modified model.For the non-ideal plasma,there is a strong coupling between ions.Nounified theoreticalmodel can completely describe the interaction between ions at arbitrary density and arbitrary temperature.In principle,the quantummolecular dynamics(QMD)can accurately describe the EOS of plasma in large density range and large temperature range.However,due tothe enormous computation and the diffi culty in converging,it is diffi cult toapply QMD tothe plasma under high temperature.W ith simple computing method and small computation,classicalmolecu lar dynamics using semi-empirical potential can calcu late the EOS accurately at high temperature.However,it will produce great error at lower temperature.It is a simple and eff ective way toobtain a global EOS by using diff erent theoreticalmodels in diff erent density range and diff erent temperature range and by interpolating in the vacant density range and vacant temperature range.

hot plasma,equations of state,thermodynamics properties

10.7498/aps.66.030505

?通信作者.E-mail:wenhuitang@163.com

?Corresponding author.E-mail:wenhuitang@163.com