李樹(shù) 陳耀樺 姬志成 章明宇 任國(guó)利 霍文義 閆威華韓小英 李志超劉杰3)4) 藍(lán)可3)4)

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

2)(中國(guó)工程物理研究院,激光聚變中心,綿陽(yáng) 621900)

3)(北京大學(xué),應(yīng)用物理與技術(shù)中心,北京 100871)

4)(上海交通大學(xué),聚變科學(xué)與應(yīng)用協(xié)同創(chuàng)新中心,上海 200240)

1 引 言

激光間接驅(qū)動(dòng)慣性約束聚變(ICF)的原理是:將激光束入射到一個(gè)內(nèi)裝靶丸的高Z元素構(gòu)成的黑腔內(nèi)壁上,內(nèi)壁吸收激光能量并以很大的份額轉(zhuǎn)換為X射線(xiàn),然后通過(guò)輸運(yùn)使黑腔內(nèi)部的X射線(xiàn)盡可能均勻燒蝕靶丸表面,驅(qū)動(dòng)靶丸實(shí)現(xiàn)球形內(nèi)爆,進(jìn)而聚變點(diǎn)火、燃燒放能[1].因此,間接驅(qū)動(dòng)研究中的一項(xiàng)重要課題是黑腔輻射場(chǎng)研究.對(duì)于靶丸內(nèi)爆而言,理想的條件是黑腔能夠提供溫度足夠高且分布均勻的輻射場(chǎng),因此,設(shè)計(jì)出激光-X光耦合效率足夠高、均勻性足夠好的黑腔是黑腔研究及激光聚變研究的重點(diǎn)之一[2,3].

間接驅(qū)動(dòng)慣性約束聚變概念提出至今已有四十多年,其間美國(guó)主要致力于柱型黑腔的理論與實(shí)驗(yàn)研究[3].為開(kāi)展激光聚變研究,美國(guó)建造了迄今世界上最大的激光裝置——國(guó)家點(diǎn)火裝置(NIF)[4],當(dāng)時(shí)計(jì)劃于2012年實(shí)現(xiàn)點(diǎn)火,但實(shí)驗(yàn)未能獲得成功[5],其中一個(gè)重要原因是柱腔中難以獲得足夠均勻的輻射場(chǎng)[6].近年,我國(guó)ICF黑腔團(tuán)隊(duì)提出一種六孔球型黑腔并開(kāi)展了系統(tǒng)的理論研究,給出了六孔球腔的概念、設(shè)計(jì)、創(chuàng)新構(gòu)型和皮實(shí)度分析[7?11].初步理論研究表明,六孔球腔具有高輻照均勻度優(yōu)勢(shì),在無(wú)需任何輔助技術(shù)手段的情況下,該型黑腔擁有全維度高輻照均勻度、高能量耦合效率、高皮實(shí)度三大特性[12].

六孔球腔具有顯著的三維特性,因此需要三維程序才能開(kāi)展合適的數(shù)值模擬研究.近年,我國(guó)研究人員發(fā)展了基于隱式蒙特卡羅(implicit Monte Carlo,IMC)方法[13]的三維輻射輸運(yùn)數(shù)值模擬程序IMC3D[14],目前可以初步應(yīng)用于黑腔物理研究.六孔球腔的實(shí)驗(yàn)尚未開(kāi)展,但先期已在我國(guó)的神光系列激光裝置上開(kāi)展了一些兩孔球型黑腔的相關(guān)實(shí)驗(yàn)[15].本文針對(duì)2015年在神光III激光裝置上實(shí)施的一次實(shí)驗(yàn)開(kāi)展數(shù)值模擬研究,比對(duì)分析數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果,檢驗(yàn)黑腔輻射場(chǎng)數(shù)值模擬能力.

2 兩孔球型黑腔實(shí)驗(yàn)概述

神光III上兩孔球腔如圖1(a)所示.黑腔形狀為球形,半徑1.8 mm,腔壁材料為金,壁厚50μm.黑腔上下端(北極、南極)開(kāi)有激光注入孔(LEH),孔半徑0.6 mm,激光由兩個(gè)LEH分別注入并在腔內(nèi)壁形成若干光斑,光斑在腔內(nèi)壁的分布(縱剖面圖)如圖1(b)所示.實(shí)驗(yàn)中共有32路(上24下8)激光在3 ns時(shí)間范圍內(nèi)均勻注入,激光總能量86.4 kJ.黑腔中心分為兩種情況:布置靶丸和不布置靶丸.靶丸為半徑0.48 mm的CH小球.黑腔放置在半徑為3 m的球形靶室中心,靶室上布置若干平板響應(yīng)式X射線(xiàn)探測(cè)器( fl at response X-ray detector,FXRD),如圖1(c)中藍(lán)點(diǎn)所示.

圖1 (a)激光注入黑腔示意圖;(b)光斑分布示意圖;(c)FXRD位置示意圖Fig.1.(a)Laser beams in hohlraum;(b)spots on the wall of hohlraum;(c)locations of FXRDs.

3 數(shù)值模擬方法與程序

兩孔球腔的結(jié)構(gòu)為二維旋轉(zhuǎn)對(duì)稱(chēng),但是考慮激光光斑的分布時(shí),輻射場(chǎng)分布將不具備二維旋轉(zhuǎn)對(duì)稱(chēng)性,用三維結(jié)構(gòu)來(lái)描述整個(gè)系統(tǒng)才更加合理.本文利用三維輻射輸運(yùn)數(shù)值模擬程序IMC3D進(jìn)行模擬.IMC3D采用隱式蒙特卡羅方法求解輻射輸運(yùn)問(wèn)題.單一離散時(shí)間步[tn,tn+1]內(nèi),求解如下IMC輻射輸運(yùn)方程[14]及物質(zhì)溫度方程:

式中I=chνn(r,?,ν,t)為輻射強(qiáng)度,n(r,?,ν,t)為相空間(r,?,ν,t)的光子數(shù)密度,c為光速,h為普朗克常數(shù),r為光子位置,?為光子飛行方向,ν為光子頻率,t為時(shí)間;其他所有帶有下標(biāo)n的物理量均表示離散時(shí)間步[tn,tn+1]起始時(shí)刻tn(或初始時(shí)刻t0)的物理量,且均為已知量:Tn為物質(zhì)溫度,bn為歸一輻射普朗克函數(shù),Qn為獨(dú)立輻射外源,cv,n為物質(zhì)比熱,fn為 fl eck因子,σa,n為吸收不透明度系數(shù),σea,n=fnσa,n為有效吸收系數(shù),σea,n=(1?fn)σa,n為有效散射系數(shù),σp,n為普朗克平均吸收截面,Ur,n為輻射能量密度,ζn為局域再發(fā)射譜.(1)式和(2)式中輻射強(qiáng)度I和物質(zhì)溫度Tn+1為待求未知變量.

對(duì)兩孔球腔問(wèn)題進(jìn)行模擬的步驟如下:

1)對(duì)球腔作空間離散,制成三維網(wǎng)格模型,作時(shí)間離散;

2)按離散時(shí)間步,將與時(shí)間相關(guān)的激光能量分別加在腔壁光斑區(qū)對(duì)應(yīng)網(wǎng)格上,求解能量(溫度)方程,得到網(wǎng)格溫度(升高),各網(wǎng)格按照溫度及發(fā)射(吸收)系數(shù)產(chǎn)生相應(yīng)的輻射光子源(Planck譜);

3)利用蒙特卡羅方法模擬跟蹤光子輸運(yùn)過(guò)程,包括飛行、碰撞吸收、散射等物理過(guò)程;

4)統(tǒng)計(jì)計(jì)數(shù)光子能量沉積(交給物質(zhì))、現(xiàn)存及漏失光子的能譜和角度譜分布等,并按照能量方程計(jì)算新的物質(zhì)溫度,開(kāi)始下一離散時(shí)間步計(jì)算;

5)完成所有離散時(shí)間步計(jì)算后,根據(jù)統(tǒng)計(jì)計(jì)數(shù)結(jié)果處理其他物理量,如黑腔、激光注入孔平均輻射溫度,靶丸表面輻射流等效溫度,靶室內(nèi)表面等效輻射溫度,探測(cè)器等效輻射溫度等.

黑腔平均輻射溫度計(jì)算公式為

式中a為輻射常數(shù),V為黑腔體積,Er,H為黑腔現(xiàn)存光子總能量.如果Er,H為黑腔中某網(wǎng)格的現(xiàn)存光子總能量,V為黑腔中某網(wǎng)格的體積,則計(jì)算得到該網(wǎng)格體的平均溫度.

LEH平均輻射溫度為

式中c為光速,Δt為時(shí)間步長(zhǎng),SLEH為L(zhǎng)EH面積,Er,L為當(dāng)前時(shí)間步從LEH漏失的光子能量.

靶丸表面輻射流等效溫度為

式中SCap為靶丸表面積,Er,C為當(dāng)前時(shí)間步從靶丸外表面流入、流出的光子能量.

FXRD等效輻射溫度為

式中SD為FXRD探頭面積(接收光子面積),θD為FXRD探頭瞄準(zhǔn)方向與LEH外法向夾角,L為FXRD與LEH之間的距離,Er,D為當(dāng)前時(shí)間步進(jìn)入FXRD的光子能量.(6)式中將FXRD處的實(shí)際輻射溫度與LEH處的平均輻射溫度進(jìn)行了可比等效轉(zhuǎn)換.

4 數(shù)值模擬結(jié)果及與實(shí)驗(yàn)結(jié)果的對(duì)比

利用IMC3D輻射輸運(yùn)程序模擬兩孔球腔實(shí)驗(yàn),由于IMC3D程序目前暫時(shí)未實(shí)現(xiàn)與流體程序的耦合,因此激光能量只能轉(zhuǎn)化成輻射能和物質(zhì)的內(nèi)能兩部分.然而實(shí)際過(guò)程中有一部分(約25%)激光能量會(huì)轉(zhuǎn)化為物質(zhì)的動(dòng)能,另外還有約5%的背反激光能量損失.將這兩個(gè)因素考慮在內(nèi),因而在實(shí)際模擬過(guò)程中將加在網(wǎng)格上的激光能量減小了30%.下面分別給出帶靶丸模型和不帶靶丸模型的數(shù)值模擬結(jié)果并作定性分析,同時(shí)將模擬結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行比較.

對(duì)于帶靶丸的模型,至激光注入結(jié)束時(shí)刻(3 ns時(shí)刻),由IMC3D模擬得到的能量分配情況如下:注入激光能量為3.2×1017MeV,物質(zhì)內(nèi)能為2.549×1017MeV,現(xiàn)存輻射能為5.482×1017MeV,總漏失輻射能為5.959×1017MeV.

激光能量大約有20%轉(zhuǎn)化為輻射能,其他約80%被物質(zhì)吸收后轉(zhuǎn)化為物質(zhì)內(nèi)能.另外,從計(jì)算結(jié)果可知,至3 ns時(shí)刻,靶丸凈吸收能量為4.18×1016MeV,約占注入激光總能量的13%,這與已有的認(rèn)知基本一致.

圖2(a)所示為黑腔、注入孔輻射溫度及靶丸入流、出流等效輻射溫度隨時(shí)間的演化情況,黑腔平均輻射溫度在3 ns時(shí)刻達(dá)到最高值,為230 eV.圖2(b)所示為靶丸表面出射輻射流與入射輻射流之比(反照率)隨時(shí)間的演化情況.圖3所示為3和4 ns時(shí)刻的黑腔輻射溫度場(chǎng)(縱剖面).從圖2和圖3可以看出,在激光注入階段,黑腔溫度與靶丸入流溫度基本一致,靶丸出流溫度低于黑腔溫度,原因在于開(kāi)始階段靶丸溫度較低,進(jìn)入靶丸的輻射能被CH材料吸收后主要轉(zhuǎn)化為物質(zhì)內(nèi)能,靶丸向外輻射的能量低于吸收能量;同時(shí),靶丸質(zhì)量較大,需要加熱的時(shí)間較長(zhǎng),到3 ns時(shí)物質(zhì)溫度與輻射溫度未能達(dá)到平衡,入射輻射流大于出射輻射流(反照率約為0.83).激光停止注入后,黑腔溫度先有短暫的快速下降過(guò)程,然后下降速度減緩,靶丸出流溫度高于黑腔溫度及入流溫度,原因在于金壁被激光加熱的質(zhì)量較小,且其吸收發(fā)射能力很強(qiáng),在沒(méi)有激光源的情況下很快與黑腔輻射場(chǎng)達(dá)到平衡,而黑腔有開(kāi)口,輻射溫度快速下降使得金壁溫度也下降較快,當(dāng)溫度下降到低于靶丸物質(zhì)溫度后,靶丸成為輻射源的凈提供體,出射輻射流大于入射輻射流(反照率超過(guò)1.0),而靶丸被加熱的質(zhì)量較大且吸收發(fā)射能力比金弱得多,因此黑腔輻射溫度緩慢降低.整個(gè)過(guò)程中注入口等效溫度略低于黑腔平均溫度,原因在于注入口不是X光轉(zhuǎn)換區(qū),而是黑腔輻射場(chǎng)的壑,從圖3可清楚看到.另外,圖3(a)中的部分腔壁及附近溫度明顯更高,原因是該腔壁位置存在激光光斑;圖3(b)中由于激光注入已經(jīng)結(jié)束,腔壁溫度十分均勻.

圖2 (a)黑腔、注入孔及靶丸表面出入流輻射溫度時(shí)間演化;(b)靶丸反照率時(shí)間演化Fig.2.(a)Temporal evolution of average radiation temperatures of hohlraum,LEH,and capsule;(b)temporal evolution of capsule re fl ection ratio.

圖3 (a)黑腔輻射溫度場(chǎng)分布(3 ns時(shí)刻);(b)黑腔輻射溫度場(chǎng)分布(4 ns時(shí)刻)Fig.3. (a)Radiation temperature distribution of hohlraum at 3 ns;(b)radiation temperature distribution of hohlraum at 4 ns.

圖4 輻射溫度隨時(shí)間的變化 (a)上0°;(b)上16°;(c)上42°;(d)下20°Fig.4.Radiation temperature versus time:(a)Up 0°;(b)up 16°;(c)up 42°;(d)down 20°.

圖4所示為從激光注入口外4個(gè)角度觀測(cè)到的腔內(nèi)輻射溫度隨時(shí)間演化的模擬與實(shí)驗(yàn)結(jié)果的對(duì)比.圖中黑色實(shí)線(xiàn)為數(shù)值模擬結(jié)果,紅色實(shí)線(xiàn)為實(shí)驗(yàn)中FXRD探測(cè)到的結(jié)果.圖4(a)為探測(cè)器布置在上部且與LEH外法向夾角為0°(正對(duì)上部LEH)時(shí)的結(jié)果.圖4(b)和圖4(c)分別為探測(cè)器與上LEH外法向夾角為16°,42°的結(jié)果,圖4(d)是探測(cè)器與下LEH外法向夾角為20°的結(jié)果.從圖中結(jié)果可以看出:數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果總體上比較接近,其中0°角方向的結(jié)果差別較大,其他三個(gè)方向的結(jié)果較為符合.造成0°角輻射溫度的差異,尤其是溫度較低階段(0—1.5 ns),初步分析主要原因如下.數(shù)值模擬沒(méi)有考慮流體運(yùn)動(dòng),故上0°角FXRD所能接收的輻射光子主要由兩部分組成,一部分由中心靶丸發(fā)射,另一部分由下腔口附近腔壁發(fā)射.但是在實(shí)際情況中,除了這兩部分之外,還有腔壁(尤其是光斑區(qū))向腔內(nèi)噴射的稀薄金等離子體所發(fā)射的光子.在靶丸溫度上升之前(1.5 ns時(shí)刻),噴射入腔內(nèi)的金等離子體可能起主導(dǎo)作用,但當(dāng)靶丸溫度足夠高后,靶丸輻射光子占主導(dǎo)地位.對(duì)于其他3個(gè)角度方向,FXRD能接收的輻射光子絕大部分來(lái)自腔壁,噴射入腔內(nèi)的金等離子體發(fā)射光子所占份額較小,而且金等離子體還對(duì)腔壁光子存在一定的阻擋、彌散作用,因此數(shù)值模擬結(jié)果略小于實(shí)驗(yàn)結(jié)果,但總體符合較好.

對(duì)于不帶靶丸的模型,至激光注入結(jié)束時(shí)刻(3 ns時(shí)刻),由IMC3D模擬得到的能量分配情況如下:注入激光能量為3.2×1017MeV,物質(zhì)內(nèi)能為2.397×1017MeV,現(xiàn)存輻射能為6.414×1015MeV,總漏失輻射能為7.390×1016MeV.與帶靶丸模型的計(jì)算結(jié)果相比,轉(zhuǎn)化為輻射能的激光能量多5%,轉(zhuǎn)化為物質(zhì)內(nèi)能的激光能量少5%.

圖5所示為不帶靶丸和帶靶丸的黑腔輻射溫度、LEH輻射溫度比較,不帶靶丸的黑腔平均輻射溫度在3 ns時(shí)刻達(dá)到最高值,為238 eV,略高于帶靶丸情況下的計(jì)算結(jié)果(230 eV).這是因?yàn)闆](méi)有靶丸吸能,更多的激光能量轉(zhuǎn)化為輻射能.圖6所示為3和4 ns時(shí)刻的黑腔輻射溫度場(chǎng)(縱剖面),可以看出除光斑區(qū)和注入口附近外,黑腔輻射場(chǎng)比較均勻.

圖5 黑腔和激光注入孔輻射溫度比較Fig.5.Temporal distributions of average radiation temperatures of hohlraum and LEH.

圖7所示為不帶靶丸情況下FXRD處輻射溫度隨時(shí)間的演化.與帶靶丸情況相似,0°角方向的結(jié)果差別較大,其他三個(gè)方向的結(jié)果符合較好.0°角方向的差異(尤其是后期)比帶靶丸情況更大,這是因?yàn)閿?shù)值模擬沒(méi)有考慮噴射入腔內(nèi)的金等離子體對(duì)FXRD的輻射貢獻(xiàn),對(duì)于不帶靶丸的情況,數(shù)值模擬中進(jìn)入0°角FXRD的只有下腔口附近腔壁發(fā)射的光子,由于沒(méi)有在后期起主導(dǎo)作用的靶丸輻射光子,因此忽略噴射入腔內(nèi)的金等離子體輻射光子的貢獻(xiàn)而引入的相對(duì)誤差必然更大.

圖7 輻射溫度隨時(shí)間的變化 (a)上0°;(b)上16°;(c)上42°;(d)下20°Fig.7.(a)Radiation temperature versus time:(a)Up 0°;(b)up 16°;(c)up 42°;(d)down 20°.

圖8所示為帶靶丸和不帶靶丸兩種情況下0°角FXRD輻射溫度的比較.在大部分時(shí)間范圍內(nèi),帶靶丸時(shí)輻射溫度比不帶靶丸時(shí)高得多,這是靶丸輻射光子進(jìn)入FXRD所致.但是,在前期(0—1.2 ns)不帶靶丸的輻射溫度更高,原因在于靶丸對(duì)0°角FXRD的輻射流(溫度)貢獻(xiàn)可以分為兩部分(圖9),其一為正作用,靶丸上半部輻射的光子進(jìn)入FXRD,其二為負(fù)作用,靶丸阻擋下腔口附近腔壁輻射的光子進(jìn)入FXRD.在靶丸溫度上升之前(約1.2 ns時(shí)刻),靶丸發(fā)射光子能力弱,吸收光子能力強(qiáng),故以負(fù)作用為主,計(jì)算得到的帶靶丸情況下FXRD溫度低于不帶靶丸情況;之后,靶丸上半部輻射的光子貢獻(xiàn)大于其阻擋作用,且隨著靶丸溫度的上升FXRD溫度差距愈加顯著.

圖8 FXRD輻射溫度比較(上0°)Fig.8.Temporal distribution of calculated FXRD radiation temperature for up 0°.

圖9 黑腔、靶丸與0°角FXRD位置示意圖Fig.9.Location of hohlraum,capsule and FXRD for up 0°.

5 結(jié) 論

利用三維IMC輻射輸運(yùn)數(shù)值模擬程序模擬神光III實(shí)驗(yàn)中兩孔球腔輻射場(chǎng),模擬結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行比較發(fā)現(xiàn):帶靶丸的黑腔平均輻射溫度最高可達(dá)230 eV,靶丸吸能效率約為13%;不帶靶丸的黑腔平均輻射溫度最高可達(dá)238 eV.在激光加源時(shí)段,靶丸表面入射輻射流高于出射輻射流,靶丸反照率逐漸上升至0.83左右;激光加源結(jié)束后,靶丸表面出射輻射流高于入射輻射流,靶丸反照率大于1.0,是輻射源的主要供體.在靶丸溫度上升之前,帶靶丸情況下0°角方向的FXRD輻射溫度低于不帶靶丸的情況,但在靶丸溫度上升之后,帶靶丸情況下0°角方向的FXRD輻射溫度遠(yuǎn)高于不帶靶丸的情況.FXRD處輻射溫度隨時(shí)間的演化與實(shí)驗(yàn)結(jié)果總體上較為一致,其中處于16°,20°和42°三個(gè)方位的結(jié)果與實(shí)驗(yàn)結(jié)果符合較好,而處于0°角方位的結(jié)果與實(shí)驗(yàn)結(jié)果差異較大,這可能是模擬中未能考慮腔壁(主要是光斑區(qū))向腔內(nèi)噴射稀薄金等離子體的影響所致.綜合數(shù)值模擬結(jié)果及其與實(shí)驗(yàn)結(jié)果的對(duì)比可知,IMC3D程序具備較好的黑腔三維輻射輸運(yùn)數(shù)值模擬能力,能夠給出較為豐富的物理量,尤其是能夠靈活方便地得出實(shí)驗(yàn)診斷類(lèi)物理量,模擬結(jié)果基本可靠.下一步需要開(kāi)展的主要工作是與三維輻射流體力學(xué)數(shù)值模擬程序相耦合,更加真實(shí)、全面地模擬黑腔輻射場(chǎng)的動(dòng)態(tài)演化過(guò)程;還要繼續(xù)開(kāi)展計(jì)算方法研究,進(jìn)一步擴(kuò)大模擬樣本規(guī)模,提高數(shù)值模擬精度以增強(qiáng)六孔球腔三維數(shù)值模擬置信度,為黑腔理論研究及實(shí)驗(yàn)設(shè)計(jì)提供參考.

[1]Zhang J,Chang T Q 2004Fundaments of the Target Physics for Laser Fusion(Beijing:National Defense Industry Press)(in Chinese)[張均,常鐵強(qiáng) 2004激光核聚變靶物理基礎(chǔ)(北京:國(guó)防工業(yè)出版社)]

[2]Atzeni S,Meyer-ter-Vehn J(Shen B F,Transl.)2008The Physics of Inertial Fusion(Beijing:Science Press)(in Chinese)[阿蔡塞,邁耶特費(fèi) (沈百飛 譯)2008慣性聚變物理(北京:科學(xué)出版社)]

[3]Lindl J D 1995Phys.Plasmas2 3933

[4]Moses E I,Boyd R N,Remington B A,Keane C J,Al-Ayat R 2009Phys.Plasmas16 041006

[5]Moses E I,Lindl J D,Spaeth M L,Patterson R W,Sawicki R H,Atherton L J,Baisden P A,Lagin L J,Larson D W,Magowan B J,Miller G H,Rardin D C,Roberts V S,van Wonterghem B M,Wegner P J 2016Fusion Sci.Technol.69 1

[6]Lindl J D 2014Phys.Plasmas21 020501

[7]Lan K,Liu J,Lai D X,Zheng W D,He X T 2014Phys.Plasmas21 010704

[8]Lan K,He X T,Liu J,Zheng W D,Lai D X 2014Phys.Plasmas21 052704

[9]Lan K,Zheng W D 2014Phys.Plasmas21 090704

[10]Huo W Y,Liu J,Zhao Y Q,Zheng W D,Lan K 2014Phys.Plasmas21 114503

[11]Li S,Lan K,Liu J 2015Laser Part.Beams15 263

[12]Lan K,Liu J,Li Z C,Xie X F,Huo W Y,Chen Y H,Ren G L,Zheng C Y,Yang D,Li S W,Yang Z W,Guo L,Li S,Zhang M Y,Han X Y,Zhai C L,Hou L F,Li Y K,Deng K L,Yuan Z,Zhan X Y,Wang F,Yuan G H,Zhang H J,Jiang B B,Huang L Z,Zhang W,Du K,Zhao R C,Li P,Wang W,Su J Q,Deng X W,Hu D X,Zhou W,Jia H T,Ding Y K,Zheng W G,He X T 2016Matter Radiat.Extremes1 8

[13]Fleck J A,Cummings J D 1971J.Comput.Phys.8 313

[14]Li S,Li G,Tian D F,Deng L 2013Acta Phys.Sin.62 249501(in Chinese)[李樹(shù),李剛,田東風(fēng),鄧力 2013物理學(xué)報(bào)62 249501]

[15]Huo W Y,Li Z C,Yang D,Lan K,Liu J,Ren G L,Li S W,Yang Z W,Guo L,Hou L F,Xie X F,Li Y K,Deng K L,Yuan Z,Zhan X Y,Yuan G H,Zhang H J,Jiang B B,Huang L Z,Du K,Zhao R C,Li P,Wang W,Su J Q,Ding Y K,He X T,Zhang W Y 2016Matter Radiat.Extremes1 2