基于MULTI2D-Z程序的Z箍縮動態(tài)黑腔形成過程模擬?

陳忠旺 寧成

(北京應(yīng)用物理與計算數(shù)學(xué)研究所,北京 100088)

基于MULTI2D-Z程序的Z箍縮動態(tài)黑腔形成過程模擬?

陳忠旺 寧成?

(北京應(yīng)用物理與計算數(shù)學(xué)研究所,北京 100088)

(2016年12月26日收到;2017年4月12日收到修改稿)

對輻射流體力學(xué)程序MULTI-2D進行改造,增加磁場演化方程程序模塊,自洽地在運動方程模塊中增加洛倫茲力,在能量方程模塊中增加歐姆加熱,將它改造成輻射磁流體力學(xué)程序MULTI2D-Z.驗證了新增磁場程序模塊的可靠性,并發(fā)現(xiàn)溫度和密度的增大會抑制磁場的擴散,負(fù)徑向速度梯度的流體對流也會抑制磁場的擴散.利用改造好的MULTI2D-Z程序模擬了峰值為8 MA的脈沖電流驅(qū)動的鎢絲陣Z箍縮動態(tài)黑腔形成過程.得到了X光功率(約30 TW)和能量(約300 kJ)、泡沫輻射溫度(約120 eV)、箍縮軌跡等模擬結(jié)果.在動態(tài)黑腔形成過程中,磁場主要分布在鎢主體等離子體中;輻射向內(nèi)傳播,燒蝕泡沫柱而使它膨脹;輻射熱波在被撞擊的泡沫柱中傳播,其傳播速度比物質(zhì)溫度傳播得快,當(dāng)輻射熱波傳播到中心軸時泡沫柱中的輻射場變得比較均勻,并且除了沖擊波處外輻射溫度與物質(zhì)溫度基本上沒有分離.這些模擬結(jié)果可增強人們對磁場擴散和對流規(guī)律以及動態(tài)黑腔形成機制的理解,同時表明了MULTI2D-Z程序可成為Z箍縮及其應(yīng)用的新的程序模擬工具.

Z箍縮,動態(tài)黑腔,輻射磁流體,MULTI-2D程序

1 引 言

Z箍縮(Z-pinch)是由驅(qū)動電流產(chǎn)生的強磁場作用于自身載流等離子體(負(fù)載)上,使其受洛倫茲力作用而向軸心內(nèi)爆的物理過程.當(dāng)?shù)入x子體箍縮到心滯止時,產(chǎn)生高溫高密度等離子體,并輻射出很強的X光,形成等離子體強輻射源.這類強X光輻射源在慣性約束聚變、實驗室天體物理、材料輻射特性和輻射效應(yīng)等研究領(lǐng)域有著廣泛的應(yīng)用[1].

近二十年來,國內(nèi)外的Z箍縮實驗和數(shù)值模擬研究均取得了很大的進展.國內(nèi)在“強光一號”和“聚龍一號”(PTS)等脈沖功率裝置上進行了大量的絲陣Z箍縮實驗研究[2?5].在美國Sandia實驗室的Z裝置上進行的絲陣實驗中,產(chǎn)生了峰值功率200 TW、能量2 MJ、脈寬5 ns的X射線脈沖,能量轉(zhuǎn)換效率高達15%[6?8].在英國帝國理工大學(xué)的Magpie裝置上也進行了大量的絲陣Z箍縮基礎(chǔ)物理實驗研究,發(fā)展了多種先進的診斷技術(shù),對絲陣消融和內(nèi)爆的細(xì)致物理過程進行了深入的觀察和研究[9?11].

在Z箍縮的諸多應(yīng)用中,其驅(qū)動慣性約束聚變應(yīng)用研究是最令人矚目的.在Z箍縮間接輻射驅(qū)動慣性約束聚變的概念設(shè)計中,需要形成一個黑腔為氘氚的靶丸內(nèi)爆創(chuàng)造一個盡可能均勻、對稱的輻射環(huán)境.目前人們提出了三種通過Z箍縮而形成的黑腔,即雙端驅(qū)動黑腔、靜態(tài)壁黑腔和動態(tài)黑腔[12].前兩種黑腔的輻射源和靶丸相對分離,其輻射場的對稱性均比較好,大體上能滿足慣性約束聚變對輻射場對稱性的要求,但其能量耦合效率較低,需要相對較高的驅(qū)動電流才能達到聚變點火條件.在動態(tài)黑腔的形成中,其核心思想是利用單層或者雙層絲陣Z箍縮內(nèi)爆等離子體撞擊內(nèi)置的中心泡沫塑料柱,內(nèi)爆等離子體撞擊塑料泡沫柱時損失的動能以及形成的沖擊波產(chǎn)生輻射,這些輻射在泡沫柱中傳播和均勻化,并且被外面高Z絲陣等離子體俘獲,而在靶丸的周圍形成一個受約束的輻射環(huán)境.因為它是在絲陣內(nèi)爆動態(tài)過程中形成的,故稱其為動態(tài)黑腔.在動態(tài)黑腔中Z箍縮內(nèi)爆等離子體、動能-X光能量轉(zhuǎn)換體(泡沫柱)和聚變靶丸等緊密地耦合在一起,從而使它們之間的能量耦合效率比較高,進而可能在驅(qū)動電流水平較低的Z箍縮裝置上實現(xiàn)慣性約束聚變點火.目前動態(tài)黑腔成了Z箍縮間接驅(qū)動慣性約束聚變的主要技術(shù)途徑之一[13].

Z裝置上的動態(tài)黑腔實驗[14,15]表明:泡沫柱的密度由15 mg/cc被壓縮到了35 mg/cc,而電子溫度由400 eV因輻射降到了300 eV;受沖擊壓縮的泡沫等離子體產(chǎn)生了大于180 kJ的輻射能量,其中靶丸吸收了超過40 kJ的輻射能量,而此能量已是靶丸點火所需能量的1/4;沖擊波對稱性受內(nèi)爆等離子體不均性和不穩(wěn)定性的影響較小,從而有利于中心聚變靶丸實現(xiàn)球?qū)ΨQ輻射燒蝕內(nèi)爆壓縮.此外在Z裝置上還首次進行了利用Z箍縮動態(tài)黑腔驅(qū)動充氘氣靶丸出中子的實驗[16],產(chǎn)生了1×1010—5×1010個熱核中子.

美國著名的激光慣性約束聚變模擬程序LASNEX已增加了磁場模塊,已應(yīng)用于Z箍縮驅(qū)動的慣性約束聚變物理研究[17,18].原來LASNEX程序是典型的三溫(電子、離子和輻射溫度)輻射流體力學(xué)程序,與國內(nèi)大部分的激光慣性約束聚變和Z箍縮模擬程序的物理模型相似,所采用的計算方法也主要是拉格朗日流體算法.它在所增加的磁場模塊中考慮了非理想磁場演化因素,如Hall效應(yīng)和Nernst效應(yīng).在空間二維上它對動態(tài)黑腔中絲陣內(nèi)爆到泡沫柱轉(zhuǎn)換體和嵌套在泡沫柱中的靶丸內(nèi)爆,以及電極板對輻射的吸收等過程進行了整體數(shù)值模擬[18].絲陣內(nèi)爆時間在100 ns左右,而緊接著的靶丸內(nèi)爆僅有幾個納秒.模擬得到的中子產(chǎn)額比實驗測量的值大兩倍,沖擊波速度比測量值高30%,靶丸心部的氣體溫度比由光譜推測的值高.模擬結(jié)果再現(xiàn)了Z箍縮內(nèi)爆到泡沫柱和緊接著產(chǎn)生的輻射沖擊波過程,還預(yù)言了泡沫柱吸收等離子體早期因歐姆加熱和內(nèi)爆而產(chǎn)生的輻射所形成的消融沖擊波,從而幫助人們更好地理解動態(tài)黑腔的形成過程.

在國內(nèi)也開展了Z箍縮驅(qū)動動態(tài)黑腔的實驗和數(shù)值模擬研究[5,19?21].目前已發(fā)表的主要是在驅(qū)動電流比較低的實驗裝置(如強光I號)上獲得的實驗結(jié)果,得到了徑向和軸向的X光輻射功率、外邊界的內(nèi)爆軌跡、泡沫柱的發(fā)光圖像等結(jié)果.數(shù)值模擬主要是在空間一維上再現(xiàn)了動態(tài)黑腔形成過程的主要物理特征,如各物質(zhì)界面的內(nèi)爆軌跡,密度、電子溫度和輻射溫度的時空演化及其相互關(guān)系,沖擊波形成和傳播過程等.

MULTI系列程序是西班牙的Ramis教授等開發(fā)的,用于模擬高能量密度等離子體中發(fā)生的輻射和物質(zhì)相互作用過程,如激光慣性約束聚變中的激光-等離子體相互作用過程.最早的MUTLI一維程序是使用Fortran語言編寫的,1994年Ramis教授開發(fā)了一種新的程序語言r94,用于編寫MULTI程序.經(jīng)過二十多年的發(fā)展,目前MUTLI系列程序有MULTI-IFE,MULTI-fs,MULTI-1D,MULTI-2D和MULTI-3D等版本.由于它們都是開源的,因此得到了比較廣泛的應(yīng)用[22?24].相對而言,一維和二維程序已發(fā)展得比較成熟.

動態(tài)黑腔負(fù)載是由柱對稱的多層物質(zhì)嵌套組合而形成,至少需要空間二維的物理模型才能比較好地描述動態(tài)黑腔的形成過程.因此,本文對國外開源的激光慣性約束聚變模擬程序MULTI-2D進行改造,使它們能應(yīng)用于Z箍縮物理、Z箍縮驅(qū)動慣性約束聚變及其聚變能應(yīng)用方面的物理研究.在物理模型上,它最大的特點是僅有一個物質(zhì)溫度,輻射溫度是通過等離子體物質(zhì)的輻射參數(shù)和近似輻射能流計算得到,因此,具有非常高的計算效率.它采用任意歐拉-拉格朗日算法和非結(jié)構(gòu)三角形網(wǎng)格,能處理在內(nèi)爆過程中產(chǎn)生的流體大變形,較好地避免在流體模擬中因網(wǎng)格畸變而無法繼續(xù)計算的問題,以及能處理由不同幾何構(gòu)形部件組合而成的復(fù)雜的內(nèi)爆物理問題.

2 物理模型和程序框架

采用輻射磁流體力學(xué)方程對等離子體Z箍縮動力學(xué)過程進行模擬[25,26].質(zhì)量守恒方程為

這里ρ是物質(zhì)密度,u為流體速度,t為時間.

動量守恒方程為

其中p為總壓強,包括流體壓強和黏性;J×B為導(dǎo)電流體受到的單位體積洛倫茲力,J為電流密度,B為磁感應(yīng)強度.在柱對稱的Z箍縮中,可以僅考慮磁場的角向分量Bθ,此時,洛倫茲力在二維柱坐標(biāo)系(r,z坐標(biāo)系)下的表達式為

這里fr和fz分別表示f在r方向和z方向的分量.能量守恒方程為

式中,ε是流體單位質(zhì)量內(nèi)能;?·(ρεu)表示單位時間流入流體單元的內(nèi)能;p?·u是壓強做功;?·F是熱傳導(dǎo)傳輸?shù)哪芰??·S是輻射能流沉積;W是歐姆加熱,其在二維柱坐標(biāo)系下的表達式為

其中,η⊥為垂直于磁場方向的等離子體電阻率,c為光速.

在輻射輸運中,流體速度遠小于光速,在任一時刻,輻射場都可以認(rèn)為是準(zhǔn)靜態(tài)的.在局部熱力學(xué)平衡附近,忽略散射和折射,輻射輸運方程為[23]

其中,n是傳播方向;qν是單位體積輻射的沉積能;λ′ν為輻射頻率為ν的輻射平均自由程,與材料的溫度和密度有關(guān)[27];Iν是輻射強度;IPν是普朗克分布函數(shù),

這里h和k分別是普朗克常數(shù)和玻爾茲曼常數(shù),T為物質(zhì)溫度.若考慮單群方法處理輻射輸運方程,則對所有頻率積分就得到

這里平均輻射自由程λ定義為

其中I,q和IP=σT4/π分別表示Iν,qν和IPν對所有頻率的積分值;σ為Stefan-Boltzmann常數(shù).在局部熱力學(xué)平衡附近,考慮近似n·?I≈n·?IP,可得輻射強度近似表達式

從而可得輻射能流

利用輻射強度,空間某一點的輻射溫度通過以下表達式進行計算:

磁場演化方程為

只考慮θ方向的磁場,在二維柱坐標(biāo)系下(7)式變成

其中,uz,ur分別為流體沿著z和r方向的速度分量;磁擴散系數(shù)vm=η⊥c2/4π.平行磁場或無磁場下Spitzer電阻率η =1.15× 10?14z lnΛT?3/2(s),其中z為平均電離度.

本文利用實用的解析表達式來計算電離度.在Thomas-Fermi理論框架下,More提出了一個平均電離度計算公式[28].設(shè)ρ為質(zhì)量密度,單位g/cm3,Te為電子溫度,單位eV(約11605 K),Z為原子序數(shù),A為相對原子質(zhì)量.定義如下過程量R=ρ/(ZA),T0=Te/Z4/3,T1=T0/(1+T0),α =14.3139,β =0.6624,x= αQβ,f(x)=則平均電離度z=f(x)Z.其中Q的計算過程如下:

相對于原來MULTI-2D程序的物理模型,為了使它能模擬研究Z箍縮及其相關(guān)的物理問題,需要增加磁場演化方程,在動量守恒方程中添加洛倫茲力,在能量守恒方程的源項中添加歐姆加熱項.此外,為了計算電阻率而增加一個電離度的解析計算表達式.這樣就把它改造成了輻射磁流體力學(xué)程序,并仍然具有高的計算效率.

MULTI-2D是結(jié)構(gòu)化程序,由若干相對獨立的功能模塊組成.這些模塊包括流體模型建立模塊、循環(huán)模塊、數(shù)據(jù)輸出和繪圖模塊,如圖1所示.在流體模型建立模塊中,包括網(wǎng)格生成、參數(shù)輸入和數(shù)據(jù)初始化.在主體循環(huán)模塊中,有外加驅(qū)動源、流體運動、能量輸運以及狀態(tài)方程等模塊.驅(qū)動源是外加能量輸入項,如激光輻射、重離子束注入,及本文增加的電流或磁場.能量模塊包括物質(zhì)能量方程,以及輻射輸運、熱傳導(dǎo)和做功等的計算.狀態(tài)方程采用QEOS程序生成的狀態(tài)方程參數(shù)[29,30].

圖1 MULTI-2D程序框架結(jié)構(gòu)Fig.1.Programm frame of MULTI-2D.

MULTI-2D程序采用r94語言編寫.該語言于1994年開發(fā),是基于C語言的高級語言,適用于UNIX操作系統(tǒng).它通過構(gòu)建大量函數(shù)庫而進行編程,使得程序精煉,可使用X窗口進行操作.

本文在MULTI-2D程序中加入磁場演化程序模塊,在流體運動模塊中加入洛倫茲力,在物質(zhì)能量方程模塊中加入歐姆加熱,將它發(fā)展成二維輻射磁流體力學(xué)數(shù)值模擬程序MULTI2D-Z,并應(yīng)用于Z箍縮動態(tài)黑腔形成過程的模擬研究.

3 模擬結(jié)果

3.1 磁場演化模塊驗證與磁場演化分析

為了驗證所添加的磁場演化方程程序模塊的正確性,考慮如圖2所示的等離子體柱.橫坐標(biāo)軸為柱等離子體的對稱軸(z軸),縱坐標(biāo)軸為柱等離子體的半徑(r軸).等離子體柱的初始半徑為1.25 cm,半高度為0.5 cm.考慮的是鎢等離子體,密度為8.6×10?4g/cm3,初始溫度為1 eV.由于柱等離子體的對稱性,因此,實際需要的計算區(qū)域僅為實線的右半部分.在z方向上加載8 MA的恒定電流,該電流在θ方向產(chǎn)生角向磁場,該磁場將向?qū)ΨQ軸方向擴散演化.為了驗證的方便,給定等離子體速度,不求解流體運動方程.此時,磁場的演化過程由方程(8)描述.磁場的邊界條件為

其中rm為等離子體柱的徑向邊界,L為z方向邊界.

明顯地,該磁場演化過程可以由徑向(r方向)一維磁場演化方程來描述.一維方程為

其邊界條件為

圖2 磁場演化模塊驗證模型Fig.2.Plasma model for veri fi cation of magnetic fi eld evolution module.

一維磁場演化方程程序模塊在柱對稱Z箍縮等離子體的模擬中是普遍使用的,其磁場演化規(guī)律清楚,程序模塊可靠.將由新添加的二維磁場演化方程程序模塊模擬出的結(jié)果與一維磁場演化程序模擬的結(jié)果進行比較,就可以判斷基于非結(jié)構(gòu)網(wǎng)格上的新增磁場程序模塊的正確性.當(dāng)流體靜止時,僅有純磁擴散,而當(dāng)存在流體運動速度梯度時,則磁場既有擴散,也有對流.下面就這兩種情況對磁場演化模塊和相應(yīng)的磁場演化規(guī)律進行考察.

圖3(a)是在假定流體靜止時,當(dāng)8 MA的恒定電流加載到該等離子體柱時,利用一維和二維磁場演化方程對其中磁場演化過程的模擬結(jié)果進行比較.恒定電流所產(chǎn)生的角向磁場僅隨r而變化,磁場逐漸向中心軸擴散,并且很快就達到穩(wěn)定分布.可見一維和二維的模擬結(jié)果符合得很好.

圖3(b)給出的是在流體流動、并且存在速度梯度時,利用一維和二維磁場演化方程對該磁場演化過程的模擬結(jié)果.這里假定流體徑向速度隨半徑的變小而下降,為某一徑向位置,um=?5×106cm/s.磁場外邊界值Bm=2I/crm隨著流體向內(nèi)箍縮運動而增大,磁場既隨流體對流,又向內(nèi)擴散.可見,此時二維和一維模擬結(jié)果也很好地重合,說明基于非結(jié)構(gòu)網(wǎng)格上的新增二維磁場演化方程程序模塊是正確的.

下面利用我們新增的基于非結(jié)構(gòu)網(wǎng)格上的二維磁場演化方程程序模塊來研究磁擴散系數(shù)和流體對流對磁場演化的影響.在模擬磁擴散系數(shù)對磁場演化的影響時,通過改變等離子體溫度或密度而改變電阻率,使磁擴散系數(shù)發(fā)生變化.這里采用圖2所示的基準(zhǔn)等離子體溫度或密度的100倍來改變電阻率(磁擴散系數(shù)),觀察(100倍溫度,基準(zhǔn)密度),(基準(zhǔn)溫度,100倍基準(zhǔn)密度)和(基準(zhǔn)溫度,基準(zhǔn)密度)等三種情況下t=100 ns時的磁場分布情況.在考察流體對流對磁場演化的影響時,使等離子體處于三種運動狀態(tài),即零速度、恒定速度u1=um=?5×106cm/s和速度隨半徑變化,并假定磁場的邊界值為恒定的Bm.為了描述磁場擴散或向軸心演化的快慢,我們定義一個磁場擴散速度.磁場向軸心擴散而逐漸減小,因此,定義磁場為Bm/e的位置r隨時間變化率為磁場擴散速度(自然常數(shù)e=2.71828).

圖4(a)是磁場擴散的模擬結(jié)果,可見密度變化對磁場分布的影響比較小,密度增加導(dǎo)致電離度減小,從而電阻率減小,磁擴散系數(shù)變小,進而抑制磁場的擴散;溫度對磁場分布的影響比較復(fù)雜,雖然溫度升高使電離度有增大的趨勢,使得電阻率增大,但電阻率與T?3/2成正比,溫度升高又使得電阻率下降得更快.綜合起來,本模擬結(jié)果表明溫度升高抑制了磁場的擴散,磁場更集中于等離子體柱的表面.這表明等離子體的導(dǎo)電性好,則電流更趨膚,是符合物理理解的.圖5(a)是與圖4(a)相同條件下的磁擴散速度隨時間的變化情況.在電流加載的初始階段,磁擴散速度比較大,但隨著時間的增加,擴散速度快速下降,最后趨于零,即磁場分布趨于穩(wěn)定.但等離子體溫度較高時,磁擴散速度初始時相對較小,下降也緩慢,或者說磁擴散系數(shù)減小,磁場擴散速度減慢,需要更長的時間才能到達穩(wěn)定.

圖3 (網(wǎng)刊彩色)二維磁場演化程序與一維磁場演化程序模擬得到的磁場徑向分布演化的比較 (a)流體靜止;(b)存在流體速度梯度Fig.3. (color online)The comparison of magnetic fi eld pro fi les in radial direction calculated byrtwodimensional and one-dimensional magnetic inductance equations:(a)ur=uz=0;(b)ur(ri)= ium,u=0.rmz

圖4 (網(wǎng)刊彩色)磁場的分布演化 (a)流體靜止時等離子體參數(shù)變化對磁場擴散的影響;(b)流體對流對磁場分布演化的影r響,u1表示等離子體以恒定速度um=?5×106cm/s向軸心運動,u2表示等離子體以空間變化速度ur(ri)=ium向軸心運動;邊界磁場固定Bm=128 T rmFig.4.(color online)The evolutions of magnetic fi eld pro fi les:(a)The e ff ects of plasma parameters on the magnetic fi eld evolution when the the plasma is at rest;(b)the e ff ect of fl uid motion on the magnetic fi eld evolution when the plasma is moving towards the axis at a constant velocity um=?r5×106cm/s,which is denoted by u1,as well as the plasma is moving at a convective velocity ur(ri)=ium,which is denoted by u.The magnetic fi eld at boundary interface is fi xed,B=128 T.rm2m

圖5 (網(wǎng)刊彩色)磁場擴散速度隨時間變化 (a)流體靜止;(b)流體運動,uo表示流體靜止;其他與圖4中的情況相對應(yīng)Fig.5.(color online)The variation of velocity of magnetic fi eld di ff usion with time:(a)The plasma is at rest;(b)the plamsa is moving at a certain velocity.The conditions used in the simulations are the same as that in Fig.4,besides the uowhich denotes the plasma at rest.

圖4(b)是在等離子體運動的情況下磁場擴散的結(jié)果.當(dāng)?shù)入x子體以恒定的速度向軸心運動時,磁場在不同時刻的徑向分布是相似的(這是因為在拉氏坐標(biāo)下當(dāng)流體速度恒定均勻時,對流項消失,從而磁場擴散與流體靜止時的情況一樣),但隨著等離子體的壓縮,內(nèi)部的磁場逐漸增大.當(dāng)?shù)入x子體以隨r逐漸減少的速度箍縮時,磁場的擴散效應(yīng)減弱,磁場的分布更趨膚.表明此時的等離子體對流抑制了磁場的擴散,或者說在實際的Z箍縮中,在箍縮初始階段已擴散進入等離子體中的磁場,當(dāng)?shù)入x子體變得快速向心運動時,磁場將有向外反擴散而再趨膚的過程.圖5(b)是與圖4(b)相同條件下的磁擴散速度隨時間的變化情況.有意思的是,當(dāng)?shù)入x子體靜止時磁擴散速度逐漸趨于零,而當(dāng)它以恒定速度流動時磁擴散速度等于等離子體運動速度5×106cm/s,即磁場被凍結(jié)在等離子體中,隨流體一起運動;當(dāng)速度不恒定,存在負(fù)的速度梯度時,磁擴散速度略小于流體運動速度,這是因為此時的磁擴散被減弱了.

3.2 絲陣Z箍縮驅(qū)動的動態(tài)黑腔形成過程

新添加的磁場演化方程程序模塊得到驗證后,改造后MULTI-2D程序的驅(qū)動源除了原來的激光和粒子束外,又增加了磁場力,從而它就具有了模擬Z箍縮及其相關(guān)物理過程的能力.典型而簡單的動態(tài)黑腔負(fù)載構(gòu)形就是在絲陣的中心嵌入一根低密度泡沫柱(其組成成分為碳和氫,通常將它簡記為CH),當(dāng)在強脈沖電流的作用下絲陣內(nèi)爆到泡沫柱,撞擊泡沫柱而在泡沫柱中形成沖擊波,絲陣和泡沫等離子體被熱化,產(chǎn)生輻射向泡沫柱中心傳播,在絲陣等離子體的部分約束下在泡沫等離子體中形成比較勻化的輻射場.為了能利用動態(tài)黑腔來驅(qū)動慣性約束聚變,這一整體動力學(xué)過程的發(fā)展、演化過程需要通過數(shù)值模擬(當(dāng)然也需要實驗)研究、理解清楚.

下面在二維柱坐標(biāo)系下,利用改造完成的MULTI2D-Z程序?qū)Ψ逯惦娏鳛镮m=8 MA的脈沖電流驅(qū)動的嵌套了中心泡沫柱的絲陣Z箍縮內(nèi)爆過程進行模擬.

圖7是計算得到的絲陣和泡沫柱的箍縮軌跡、X光功率和能量、絲陣和泡沫柱的收縮比、電流波形等隨內(nèi)爆過程的變化.整體箍縮到心在150 ns左右,X光功率峰值接近30 TW,總能量在330 kJ左右,絲陣和泡沫的最大壓縮比分別為15和5左右.泡沫柱在被絲陣等離子體碰撞壓縮以前,略微膨脹.

圖6 動態(tài)黑腔負(fù)載模型Fig.6. The used model in simulation of dynamic hohlraum.

圖7 (網(wǎng)刊彩色)動態(tài)黑腔形成過程中一些宏觀參量隨時間的變化 (a)絲陣外半徑R、電流I、X光功率P和徑向輻射能量E;(b)絲陣和泡沫的柱外半徑(rWA,rCH),絲陣和泡沫柱的收縮比(Ratio 2,Ratio 1)Fig.7.(color online)The variations of some macroscopical quantities with time during the formation process of dynamic hohlraum:(a)Outer radius of wire-array(R),current(I),X-ray power(P)and energy(E);(b)outer radii of wire-array(rWA)and foam(rCH),convergence ratioes of wire-array(Ratio 2)and foam(Ratio 1).

圖8—圖10分別是絲陣等離子體內(nèi)爆碰撞泡沫柱后三個時刻(146,150,152 ns)密度、等離子體溫度和輻射溫度、內(nèi)爆速度、磁場、輻射場(溫度)等的空間分布情況.t=146 ns時,絲陣內(nèi)爆到了半徑為0.3 cm的地方,離泡沫柱的初始位置還有0.1 cm.此時在絲陣等離子體中形成了很強的沖擊波(見圖8(a)—(c)中的密度尖峰位置),內(nèi)爆速度達到了4×107cm/s,物質(zhì)溫度達到了370 eV左右,泡沫中的最大輻射溫度接近100 eV,邊界磁場達到約520 T.在靠近中心軸處的輻射溫度比物質(zhì)溫度高,表明輻射熱波向軸心方向傳播比物質(zhì)溫度快.在箍縮過程中磁場主要分布在主體鎢等離子體中,形成很強的磁壓,使等離子體加速向內(nèi)箍縮,在壓縮做功和歐姆加熱作用下,主體等離子體溫度較高,密度較大.溫度升高和密度增加使磁擴散系數(shù)減小,另外流體壓縮也抑制磁場的擴散.高溫高密度主體絲陣等離子體向泡沫內(nèi)部和外表面及兩端電極發(fā)射輻射,向內(nèi)輻射加熱燒蝕泡沫柱,使泡沫柱因溫度升高而向外膨脹,此時,泡沫柱的外半徑膨脹到了大約0.24 cm.向外表面(r方向)和電極(z方向)的輻射導(dǎo)致內(nèi)爆等離子體的能量漏失,使得靠近電極板區(qū)域和等離子體外則的輻射溫度均比較低.

到了150 ns(見圖9),絲陣等離子體內(nèi)爆到全部進入了原泡沫柱所在的區(qū)域,到達半徑為0.17 cm的位置,表明泡沫柱已被絲陣等離子體碰撞,泡沫柱得到了一點壓縮.此時在絲陣等離子體中的沖擊波(見圖9(a)—(c)中的密度尖峰位置)得到增強,溫度和密度增高,但內(nèi)爆速度下降了一點,最大內(nèi)爆速度下降到3.7×107cm/s,物質(zhì)溫度升高到460 eV左右,泡沫中輻射溫度在中心區(qū)域基本均勻化,其值達約120 eV,并且其中的物質(zhì)溫度和輻射溫度幾乎沒有分離,達到了平衡,邊界磁場超過了900 T.鎢等離子碰撞低密度泡沫后形成物質(zhì)沖擊波,速度減緩,動能轉(zhuǎn)換為內(nèi)能,溫度迅速升高,形成輻射沖擊波;輻射沖擊波迅速向內(nèi)傳播,在物質(zhì)沖擊波到來之前形成了較為均勻的黑腔輻射環(huán)境,黑腔中物質(zhì)溫度和輻射溫度較為重合,溫度達120 eV左右;與此同時,主體等離子體向外輻射出X光的第一個小波峰,功率在7 TW左右(見圖7(a)).在動態(tài)黑腔概念設(shè)計中,希望在絲陣碰撞壓縮泡沫柱到達靶丸之前,所形成的黑腔輻射場對靶丸的燒蝕壓縮,并使其達到聚變點火條件.因此,碰撞后主體等離子的速度大幅度減小,損失的動能轉(zhuǎn)化為內(nèi)能和輻射,大部分輻射被泡沫吸收,這是黑腔輻射場形成的重要機制.

圖8 (網(wǎng)刊彩色)146 ns時刻各參量的隨空間分布 (a)輻射溫度Tr、物質(zhì)溫度T和密度分布ρ;(b)速度ur和密度ρ分布;(c)磁場B和密度ρ分布;(d)輻射溫度Tr空間分布Fig.8.(color online)The spatial pro fi les of some physical quantities at 146 ns:(a)Radiation,plamsa temperatures,and density pro fi les;(b)velocity and density pro fi les;(c)magnetic fi eld and density pro fi les;(d)spatial distribution of radiation temperature.The part of a line with symbols of circle belongs to foam,and it with symbols of plus belongs to wire-array.

圖9 (網(wǎng)刊彩色)150 ns時各時刻參量隨空間變化 (a)輻射溫度、物質(zhì)溫度和密度分布;(b)速度和密度分布;(c)磁場和密度分布;(d)輻射溫度空間分布Fig.9.(color online)The spatial pro fi les of some physical quantities at 150 ns:(a)Radiation,plamsa temperatures,and density pro fi les;(b)velocity and density pro fi les;(c)magnetic fi eld and density pro fi les;(d)spatial distribution of radiation temperature.The part of a line with symbols of circle belongs to foam,and it with symbols of plus belongs to wire-array.

t=152 ns以后,絲陣和泡沫等離子體繼續(xù)箍縮,泡沫等離子體整體得到進一步的加速,泡沫中的溫度和密度進一步增大,在泡沫中出現(xiàn)小的密度峰,中心附近的輻射和物質(zhì)溫度稍高.由于磁壓力做功,使得泡沫的溫度和輻射溫度繼續(xù)增加,最高達到約170 eV,它們?nèi)匀换旧蠜]有分離.輻射溫度在r方向分布的均勻性變?nèi)?見圖10(a)),但對比圖9(d)可以看到,高溫和比較均勻的輻射區(qū)域在z方向有所增大.由于邊界輻射溫度的增大,使得邊界輻射能量以更高的速率損失(輻射能流約為導(dǎo)致徑向輻射溫度相對更低.因此,對于動態(tài)黑腔驅(qū)動聚變靶丸的實際應(yīng)用而言,需要根據(jù)驅(qū)動電流、靶丸大小仔細(xì)地設(shè)計動態(tài)黑腔的負(fù)載構(gòu)形,以明確能夠利用的動態(tài)黑腔中輻射場的時間窗口和空間區(qū)域,使其靶丸實現(xiàn)對稱壓縮和點火成功.在整體箍縮到心時,等離子體減速滯止,產(chǎn)生強烈的X光輻射,此時X光輻射脈沖的峰值功率達約28 TW.

在徑向一維上密度、物質(zhì)溫度(電子溫度)和輻射溫度的空間分布和隨時間的演化關(guān)系的整體特征,如絲陣等離子體的前期膨脹、輻射燒蝕泡沫柱而膨脹、輻射熱波傳播較快、輻射熱波與輻射沖擊波的空間分離等,與一維拉氏非平衡輻射磁流體力學(xué)程序CRMHA的模擬結(jié)果[20]以及國外的模擬結(jié)果[18]都是相似的.并且預(yù)期的黑腔輻射溫度與文獻[20]的結(jié)果也接近,都是120 eV左右.說明MULTI2D-Z程序模擬動態(tài)黑腔形成過程所得到的物理結(jié)果是基本正確、可靠的.但本文中泡沫柱的輻射溫度的均勻化過程要慢一些,即泡沫等離子體的輻射吸收更強烈一些,并且本文能給出輻射場在徑向和軸向方向上的均勻化情況.

4 結(jié) 論

從輻射磁流體力學(xué)基本方程出發(fā),通過在MULTI-2D程序中增加磁場演化方程程序模塊,同時自洽地在流體動量守恒方程中增加導(dǎo)電流體所受到的洛倫茲力,在能量守恒方程中增加歐姆加熱能量源項,并在相應(yīng)的程序中進行改動.通過新增程序模塊的驗證,使得MULTI-2D程序改造、發(fā)展成了MULTI2D-Z程序,成為基于非結(jié)構(gòu)網(wǎng)格的二維輻射磁流體力學(xué)程序,可以應(yīng)用于Z箍縮驅(qū)動的動態(tài)黑腔形成動力學(xué)過程模擬.

圖10 (網(wǎng)刊彩色)152 ns時各時刻參量隨空間變化 (a)輻射溫度、物質(zhì)溫度和密度分布;(b)速度和密度分布;(c)磁場和密度分布;(d)輻射溫度空間分布Fig.10.(color online)The spatial pro fi les of some physical quantities at 152 ns:(a)Radiation,plamsa temperatures,and density pro fi les;(b)velocity and density pro fi les;(c)magnetic fi eld and density pro fi les;(d)spatial distribution of radiation temperature.The part of a line with symbols of circle belongs to foam,and it with symbols of plus belongs to wire-array.

利用新研制的磁場演化方程程序模塊,對磁場的擴散和對流演化規(guī)律進行了一些研究.發(fā)現(xiàn)溫度和密度的增大會抑制磁場的擴散;負(fù)徑向速度梯度的流體(如Z箍縮等離子體的箍縮運動)也會抑制磁場的擴散.當(dāng)?shù)入x子體靜止時,磁擴散速度隨著時間的增加,迅速減少,然后趨于零,磁場分布變得穩(wěn)定;當(dāng)?shù)入x子體以恒速箍縮時,磁擴散速度也是電流加載后快速下降,然后趨同于流體運動速度,表明磁場被凍結(jié)等離子體中,隨流體一起運動;當(dāng)?shù)入x子體以隨半徑變小而減少的速度做箍縮運動時,磁場的擴散速度略小于運動速度,表明磁場擴散過程被削弱了,出現(xiàn)磁場的再趨膚現(xiàn)象.通過對磁場演化的分析,有助于理解Z箍縮過程中各種因素對磁場演化的影響趨勢.

利用MULTI2D-Z程序?qū)? MA電流水平的絲陣Z箍縮驅(qū)動的動態(tài)黑腔形成動力學(xué)過程進行了模擬,再現(xiàn)了該動態(tài)黑腔的形成過程.動態(tài)黑腔形成過程的整體特征與國內(nèi)外的模擬結(jié)果基本相似,黑腔輻射溫度與相同電流水平下一維CRMHA程序的模擬結(jié)果接近,說明了該程序用于動態(tài)黑腔形成過程模擬的可靠性.本模擬中發(fā)現(xiàn),在Z箍縮等離子體的高溫高密度和對流的共同作用下,磁場主要集中在鎢等離子體表面,使絲陣等離子體快速內(nèi)爆,并撞擊泡沫柱;鎢主體等離子體的輻射向泡沫柱傳播,使泡沫燒蝕而膨脹;絲陣等離子體碰撞泡沫柱,形成很強的物質(zhì)沖擊波和輻射,輻射向泡沫等離子體中傳播,形成輻射熱波;在泡沫等離子體中輻射熱波的傳播比物質(zhì)溫度的傳輸速度快,當(dāng)輻射熱波傳播到中心軸后,泡沫中的輻射達到了較好的均勻化(約120 eV),并且除了沖擊波附近外,輻射溫度和物質(zhì)溫度基本上重合,沒有分離現(xiàn)象;此后,隨著整體箍縮的進行,泡沫等離子體繼續(xù)被加速和壓縮,其中的輻射溫度繼續(xù)提升,輻射場在軸向方向的均勻性有所增加,但在徑向方向的均勻性有所下降.當(dāng)泡沫柱中嵌套有聚變靶丸時,需要通過程序模擬謹(jǐn)慎、精確地設(shè)計絲陣和泡沫柱的尺寸和質(zhì)量,以便能得到滿足聚變點火要求的動態(tài)黑腔輻射場時間窗口和空間區(qū)域.

動態(tài)黑腔形成過程的數(shù)值模擬對相關(guān)實驗的設(shè)計和物理過程的理解,以及Z箍縮慣性約束聚變概念設(shè)計是非常重要的.今后還需要繼續(xù)完善和發(fā)展MULTI2D-Z程序,改善其輻射模型,評估非結(jié)構(gòu)網(wǎng)格劃分對模擬結(jié)果的影響.

感謝北京應(yīng)用物理與計算數(shù)學(xué)研究所的文武老師、薛創(chuàng)老師、鄔吉明研究員和高志明老師以及中國工程物理研究院研究生部的吳福源師兄在研究工作中給予的指導(dǎo)和幫助.感謝中國工程物理研究院聚龍一號(PTS)裝置的Z箍縮實驗團隊對絲陣Z箍縮和動態(tài)黑腔實驗結(jié)果的交流和分享.

[1]Matzen M K 1997 Phys.Plasmas 4 1519

[2]Qiu A C,Kuai B,Wang L P,Wu G,Cong P T 2008 High Power Laser and Particle Beams 20 1761(in Chinese)[邱愛慈,蒯斌,王亮平,吳剛,叢培天2008強激光與粒子束20 1761]

[3]Huang X B,Zhou S T,Dan J K,Ren X D,Wang K L,Zhang S Q,Li J,Xu Q,Cai H C,Duan S C,Ouyang K,Chen G H,Ji C,Wei B,Feng S P,Wang M,Xie W P,Deng J J,Zhou X W,Yang Y 2015 Phys.Plasmas 22 072707

[4]Dan J K,Ren X D,Huang X B,Zhang S Q,Zhou S T,Duan S C,Ouyang K,Cai H C,Wei B,Ji C,He A,Xia M H,Feng S P,Wang M,Xie W P 2013 Acta Phys.Sin.62 245201(in Chinese)[但加坤,任曉東,黃顯賓,張思群,周少彤,段書超,歐陽凱,蔡紅春,衛(wèi)兵,計策,何安,夏明鶴,豐樹平,王勐,謝衛(wèi)平2013物理學(xué)報62 245201]

[5]Jiang S Q,Ning J M,Chen F X,Ye F,Xue F B,Li L B,Yang J L,Chen J C,Zhou L,Qin Y,Li Z H,Xu R K,Xu Z P 2013 Acta Phys.Sin.62 155203(in Chinese)[蔣樹慶,甯家敏,陳法新,葉繁,薛飛彪,李林波,楊建倫,陳進川,周林,秦義,李正宏,徐榮昆,許澤平 2013物理學(xué)報62 155203]

[6]Deeney C,Nash T J,Spielman R B,Seaman J F,Chandler G C,Struve K W,Porter J L,Stygar W A,McGum J S,Jobe D O,Gilliland T L,Torres J A,Vargas M F,Ruggles L E,Breeze S,Mock R C,Douglas M R,Fehl D L,McDaniel D H,Matzen M K,Peterson D L,MatuskaW,Roderick N F,MacFarlane J J 1997 Phys.Rev.E 56 5945

[7]Spielman R B,Deeney C,Chandler G A,Douglas M R,Fehl D L,Matzen M K,McDaniel D H,Nash T J,PorterJ L,Sanford T W L,Seaman J F,Stygar W A,Struve K W,Breeze S P,McGurn J S,Torres J A,Zagar D M,Gilliland T L,Jobe D O,McKenney J L,Mock R C,Vargas M,Wangoner T,Peterson D L 1998 Phys.Plasmas 5 2105

[8]Stygar W A,Ives H C,Fehl D L,Fehl D L,Cuneo M E,Mazarakis M G,Bailey J E,Bennett G R,Bliss D E,Chandler G A,Leeper R J,Matzen M K,McDaniel D H,McGurn J S,McKenney J L,Mix L P,Muron D J,Porter J L,Ramirez J S,Ruggles L E,Seamen J F,Simpson W W,Speas C S,Spielman R B,Struve K W 2004 Phys.Rev.E 69 046403

[9]Lebedev S V,Mitchell I H,Aliaga-Rossel R,Bland S N,Chittenden J P,Dangor A E,Haines M G 1998 Phys.Rev.Lett.81 4152

[10]Harvey-Thompson A J,Lebedev SV,Patankar S,Bland S N,Burdiak G,Chittenden J P,Colaitis A,de Grouchy P,Doyle H W,Hall G N,Khoory E,Hohenberger M,Pickworth L,Suzuki-Vidal F,Smith R A,Skidmore J,Suttle L,Swadling G F 2012 Phys.Rev.Lett.108 145002

[11]Swadling G F,Lebedev S V,Harvey-Thompson A J,Rozmus W,Burdiak G,Suttle L,Patankar S,Smith R A,Bennett M,Hall G N,Suzuki-Vidal F,Bland S,Yuan J 2015 Phys.Plasmas 22 072706

[12]Leeper R J,Alberts T E,Asay J R,Baca P M,Baker K L,Breeze S P,Chandler G A,Cook D L,Cooper G W,Deeney C,Derzon M S,Douglas M R,Fehl D L,Gilliland T,Hebron D E,Hurst M J,Jobe D O,Kellogg J W,Lash J S,Lazier S E,Matzen M K,McDaniel D H,McGurn J S,Mehlhorn T A,Moats A R,Mock R C,Muron D J,Nash T J,Olson R E,Porter J L,Quintenz J P,Reyes P V,Ruggles L E,Ruiz C L,Sanford T W L,Schmidlapp F A,Seamen J F,Spielman R B,Stark M A,Struve K W,Stygar W A,Tibbetts-Russell D R,Torres J A,Vargas M,Wagoner T C,Wake fi eld C,Hammer J H,Ryutov D D,Tabak M,Wilks S C,Bowers R L,McLenithan K D,Peterson D L 1999 Nucl.Fusion 39 1283

[13]Ning C,Feng Z X,Xue C 2014 Acta Phys.Sin.63 125208(in Chinese)[寧成,豐志興,薛創(chuàng) 2014物理學(xué)報63 125208]

[14]Bailey J E,Chandler G A,Slutz S A,Bennett G R,Cooper G,Lash J S,Lazier S,Lemke R,Nash T J,Nielsen D S,Moore T C,Ruiz C L,Schroen D G,Smelser R,Torres J,Vese y R A 2002 Phys.Rev.Lett.89 095004

[15]Rochau G A,Bailey J E,Maron Y,Chandler G A,Dunham G S,Fisher D V,Fisher V I,Lemke R W,Mac-Farlane J J,Peterson K J,Schroen D G,Slutz S A,Stambulchik E 2008 Phys.Rev.Lett.100 125004

[16]Ruiz C L,Cooper G W,Slutz S A,Bailey J E,Chandler G A,Nash T J,Mehlhorn T A,Leeper R J,Fehl D,Nelson A J,Franklin J,Ziegler L 2004 Phys.Rev.Lett.93 015001

[17]Zimmerman G B,Kruer W B 1975 Comments Plasma Phys.Control.Fusion 2 51

[18]Slutz S A,Peterson K J,Vesey R A,Lemke R W,Bailey J E,Varnum W,Ruiz C L,Cooper G W,Chandler GA,Rochau G A,Mehlhorn T A 2006 Phys.Plasmas 13 102701

[19]Meng S J,Huang Z C,Ning J M,Hu Q Y,Ye F,Qin Y,Xu Z P,Xu R K 2016 Acta Phys.Sin.65 075201(in Chinese)[蒙世堅,黃展常,甯家敏,胡青元,葉繁,秦義,許澤平,徐榮昆2016物理學(xué)報65 075201]

[20]Xiao D L,Sun S K,Xue C,Zhang Y,Ding N 2015 Acta Phys.Sin.64 235203(in Chinese)[肖德龍,孫順凱,薛創(chuàng),張揚,丁寧2015物理學(xué)報64 235203]

[21]Xiao D L,Ding N,Ye F,Ning J M,Hu Q Y,Chen F X,Qin Y,Xu R K,Li Z H,Sun S K 2014 Phys.Plasmas 21 042704

[22]Ramis R,Ramirez J,Meyer-ter-Vehn J 1988 Comput.Phys.Commun.49 475

[23]Ramis R,Meyer-ter-Vehn J,Ramirez J 2009 Comput.Phys.Commun.180 977

[24]Ramis R,Meyer-ter-Vehn J,Ramirez J 2016 Comput.Phys.Commun.203 226

[25]Liberman M A,De Groot J S,Toor A 2003 Physics of High-Density Z-Pinch Plasmas(Beijing:National Defence Industry Press)(in Chinese)[Liberman M A,De Groot J S,Toor A著(孫承緯譯)2003高密度Z箍縮等離子體物理學(xué)(北京:國防工業(yè)出版社)

[26]Ning C,Yang Z H,Ding N 2003 Acta Phys.Sin.52 415(in Chinese)[寧成,楊震華,丁寧 2003物理學(xué)報52 415]

[27]Murakami M,Meyer-ter-Vehn J,Ramis R 1990 J.X-Ray Sci.Technol.2 127

[28]Salamann D 1998 Atom Physics in Hot Plasmas(New York:Oxford University Press)p28

[29]Kemp A J,Meyer-ter-Vehn J 1998 Nucl.Instrum.Meth.A 415 674

[30]More R,Warren K H,Yong D A,Zimmermann G 1988 Phys.Fluids 31 3059

PACS:52.58.Lq,52.59.Qy,52.30.Cv,52.38.PhDOI:10.7498/aps.66.125202

Simulation of forming process of Z-pinch dynamic
hohlraum based on the program MULTI2D-Z?

Chen Zhong-Wang Ning Cheng?

(Institute of Applied Physics and Computational Mathematics,Beijing 100088,China)

26 December 2016;revised manuscript

12 April 2017)

The radiation hydrodynamics code MULTI-2D,which was developed by Ramis et al.in 2009(2009 Comput.Phys.Commun.180 977)and adopted the single temperature fl uid and unstructured lagrangian mesh,is modi fi ed into a radiation magnetohydrodynamics code MULTI2D-Z by adding the program module of evolution equation of magnetic fi eld,and self-consistently considering the Lorentz force in the module of motion equation and the Ohmic heating in the module of energy equation.The newly developed module for magnetic fi eld was validated to be reliable.The module is used to study the magnetic fi eld di ff usion process,and it is found that the di ff usion is weakened due to the increasing of plasma temperature and density and the fl uid convection,in which there is minus grads of velocity in radial direction.The new code MULTI2D-Z is used to simulate the formation process of dynamic hohlraum driven by tungsten wire-array Z-pinch at an 8 MA current level.The obtained results are that X-ray power and energy are,respectively,～30 TW and～300 kJ,radiation temperature in foam is～120 eV,and the implosion trajectory of wire-array is also obtained.The calculated results reveal that the magnetic fi eld is mainly distributed in the outside of tungsten plasma during the hohlraum formation.The foam expands due to the radiation heating from the shock wave created by the collision between wire-array plasma and the foam.The thermal radiation wave,which is characterized by radiation temperature,spreads towards the central axis faster than the plasma temperature.When the thermal radiation wave spreads to the central axis,the radiation temperature becomes comparatively uniform in space,and is almost equal to the plasma temperature except at the place of the shock wave.These results help the people to better understand the magnetic fi eld di ff usion and convection in Z-pinch,as well as the formation mechanism of dynamic hohlraum driven by wire-array Z-pinch.It is also indicated that the newly developed code MULTI2D-Z can be considered as a new tool for simulating Z-pinch and its applications,such as inertial con fi nement fusion and magnetically accelerated fl yer plates.

Z-pinch,dynamic hohlraum,radiation magnetohydrodynamics,MULTI-2D code

10.7498/aps.66.125202

?國家自然科學(xué)基金(批準(zhǔn)號:11675025)和國家自然科學(xué)基金重點項目(批準(zhǔn)號:11135007)資助的課題.

?通信作者.E-mail:ning_cheng@iapcm.ac.cn

?2017中國物理學(xué)會Chinese Physical Society

http://wulixb.iphy.ac.cn

*Project supported by the National Natural Science Foundation of China(Grant No.11675025)and the Key Program of the National Natural Science Foundation of China(Grant No.11135007).

?Corresponding author.E-mail:ning_cheng@iapcm.ac.cn