孫貝貝 葉文華 張維巖

1) (中國工程物理研究院研究生院,北京 100088)

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

1 引言

最近,美國國家點火裝置(NIF)上的慣性約束聚變(inertial confinement fusion,ICF)實 驗實現(xiàn)了能量凈增益,標志著ICF 研究的巨大進展.充分了解流體不穩(wěn)定性對內(nèi)爆性能的影響對ICF 點火是十分必要的[1?5].目前研究較多的是中心點火方式,通過內(nèi)爆壓縮在靶丸中心產(chǎn)生點火熱斑.為了形成點熱斑,需要產(chǎn)生極高的殼層內(nèi)爆速度(不小于350 km/s)[6],加速階段燒蝕面上會產(chǎn)生Rayleigh-Taylor (RT)不穩(wěn)定性[7?10].擾動種子在加速階段會被RT 不穩(wěn)定性指數(shù)放大,同時擾動通過Feed-in 過程耦合到內(nèi)界面(DT 冰和DT 氣體界面)引起內(nèi)界面擾動增長,形成減速階段內(nèi)界面RT 不穩(wěn)定性的擾動種子.流體不穩(wěn)定性發(fā)展會破壞靶丸的均勻性和對稱性,導(dǎo)致能量增益降低甚至點火失敗.內(nèi)爆中流體不穩(wěn)定性是不可避免的,ICF 研究工作的重點之一是確定并減輕不穩(wěn)定性來源,通過減小擾動種子來抑制不穩(wěn)定性增長[11?13].

燒蝕面的擾動種子來源于靶丸缺陷、輻照能量不均勻性以及靶丸定位偏差等因素.材料內(nèi)部的密度擾動(同種材料內(nèi)部的不均勻性)是靶丸缺陷的一種.在制靶過程中燒蝕材料和DT 冰內(nèi)部會產(chǎn)生密度擾動.另外,DT 燃料中的氚衰變會將能量沉積到燒蝕材料和DT 冰殼層中,在DT 冰中產(chǎn)生空泡,并引起燒蝕材料的局部膨脹[14].

類似界面的Richtmyer-Meshkov (RM)不穩(wěn)定性[15,16],沖擊波穿過密度擾動時,由于斜壓效應(yīng)而產(chǎn)生渦量的沉積,引起不穩(wěn)定性增長,因此將這種不穩(wěn)定性稱為類RM 不穩(wěn)定性[17?21].類RM 不穩(wěn)定性耦合到初始無擾動界面,會引起界面的不穩(wěn)定性增長.最新的研究工作表明,在ICF 中考慮初始密度擾動對內(nèi)爆性能的影響十分重要[13,22?26].Haines 等[24,25]工作表明,間接驅(qū)動內(nèi)爆中密度擾動產(chǎn)生的流體力學(xué)不穩(wěn)定性會降低燃料的壓縮性.文獻[13,26]表明了材料內(nèi)部缺陷(包括密度擾動和孤立缺陷)會耦合到燒蝕面,形成加速階段RT 不穩(wěn)定性的種子,經(jīng)加速階段的RT 不穩(wěn)定性放大降低內(nèi)爆性能.殼層厚度和擾動放置位置會影響燒蝕面擾動種子的大小,燒蝕面附近的擾動會產(chǎn)生最大的燒蝕面擾動種子.沖擊波與密度擾動相互作用產(chǎn)生的熵渦波凍結(jié)在波后流體中,由于質(zhì)量燒蝕效應(yīng),渦量會向燒蝕面輸運,影響燒蝕面擾動幅值.

目前,針對密度擾動引起的流體不穩(wěn)定性及其對內(nèi)爆性能影響的研究工作還比較少,對相關(guān)的物理認識仍然存在不足.本文采用數(shù)值模擬的方法,研究了線性階段類RM 不穩(wěn)定性的增長規(guī)律及不穩(wěn)定性到界面的耦合規(guī)律.第2 節(jié)介紹了數(shù)值模擬的初始設(shè)置.第3 節(jié)研究了沖擊波與密度擾動相互作用后,線性階段密度擾動的類RM 不穩(wěn)定性增長規(guī)律.分析了不同擾動波長λy、沖擊波馬赫數(shù)Ma和密度擾動大小η對不穩(wěn)定性增長的影響.第4 節(jié)研究了密度擾動的不穩(wěn)定性增長與界面的耦合規(guī)律.分析了不同密度擾動位置、界面Atwood數(shù)對密度擾動與界面耦合的影響.第5 節(jié)為總結(jié)與討論.

2 數(shù)值模擬初始設(shè)置

本文所有算例計算使用組內(nèi)開發(fā)的二維歐拉程序,已經(jīng)過多個算例驗證和考核[27].控制方程為無黏、無傳熱的歐拉方程,狀態(tài)方程為理想氣體狀態(tài)方程.在本文算例中,各部分氣體絕熱指數(shù)均為γ=5/3.

圖1 給出了數(shù)值模擬初始設(shè)置示意圖.Ⅲ區(qū)為沖擊波上游區(qū)域,存在密度擾動,密度擾動的中心位置為x0,密度擾動在y方向上以余弦函數(shù)的形式周期分布,在x方向上由擾動中心位置向兩側(cè)衰減.Ⅲ區(qū)中的密度分布為[13]

圖1 數(shù)值模擬初始設(shè)置示意圖Fig.1.Schematics of the initial configuration.

式中,ρ3=1 g/cm3(物理量下標1,2,3 分別代表Ⅰ,Ⅱ,Ⅲ區(qū))為CH 燒蝕材料的密度,η為密度擾動的相對大小,η的取值區(qū)間為[0,1),本文中η的取值較小以保障擾動的不穩(wěn)定性增長處于線性階段.kx,y=2π/λx,y,λy為擾動波長,λx表征密度擾動在x方向上的寬度.λx越大,密度擾動在x方向上衰減越慢,密度擾動區(qū)域越寬.Ⅲ區(qū)設(shè)置為等壓區(qū),p3=0.017 Mbar (1 Mbar=1011Pa)為CH材料溫度為300 K 時的壓強.沖擊波沿x正向傳播,Ⅱ區(qū)為沖擊波后區(qū)域.Ⅱ區(qū)中的密度、速度和壓強通過求解朗金-雨貢紐關(guān)系式得到.Ⅱ區(qū)和Ⅰ區(qū)之間為接觸間斷,界面兩側(cè)速度和壓強相等,密度不同.界面上Atwood 數(shù)At=(ρ2-ρ1)/(ρ2+ρ1).令界面兩側(cè)密度相等時,Ⅱ區(qū)和Ⅰ區(qū)之間的界面不存在,研究密度擾動的類RM 增長規(guī)律.程序中空間項的離散使用五階WENO 格式,時間步的離散使用兩步龍格-庫塔方法.x方向計算寬度為200 μm,y方向計算寬度為λy,使用128 個網(wǎng)格來分辨不同波長的擾動,網(wǎng)格寬度 Δx=Δy=λy/128.x方向上采用出流邊界條件,y方向上采用周期邊界條件.

3 線性階段密度擾動的類RM 不穩(wěn)定性

在ICF 內(nèi)爆中,驅(qū)動能量輻照在靶丸表面,燒蝕出噴射等離子體產(chǎn)生燒蝕壓,形成向靶丸內(nèi)部傳播的沖擊波.沖擊波經(jīng)過存在密度擾動的流體時,密度梯度與壓強梯度不平行,斜壓效應(yīng)會在密度擾動區(qū)域沉積渦量,產(chǎn)生不穩(wěn)定性增長.這種不穩(wěn)定性增長機制與沖擊波和擾動界面相互作用后產(chǎn)生的RM 不穩(wěn)定性是相同的.將這種不穩(wěn)定性稱為類RM 不穩(wěn)定性.為了研究密度擾動的類RM 不穩(wěn)定性到界面的耦合,首先研究類RM 不穩(wěn)定性的增長規(guī)律.

本節(jié)中算例選取λx=5 μm.圖2 給出了算例SPI1(算例設(shè)置見表1)在1 ns 和3 ns 時的渦量ω云圖.藍色圈起部分為密度擾動區(qū)域,紅色虛線標注了沖擊波位置.沖擊波與密度擾動相互作用之后產(chǎn)生了旋轉(zhuǎn)方向相反的渦對,流場中渦量最大的地方位于密度擾動區(qū)域.沖擊波與密度擾動相互作用之后,沖擊波陣面上產(chǎn)生擾動,在沖擊波和密度擾動區(qū)域之間產(chǎn)生凍結(jié)在流體上的熵渦波.平面沖擊波擾動幅值是振蕩衰減的,擾動沖擊波產(chǎn)生的渦量隨著到密度擾動區(qū)域的距離增大而衰減.密度擾動與沖擊波之間的區(qū)域近似等壓.沖擊波與密度擾動作用之后同樣會產(chǎn)生向左傳播的反射波,在η較大時密度擾動左側(cè)區(qū)域也幾乎是沒有熵增的,反射波可以做聲學(xué)近似.

表1 不同算例的參數(shù)設(shè)置Table 1.Initial physical parameters in different cases.

圖2 算例SPI1 在1 ns (a)和3 ns (b)時的渦量場.藍色圈起部分為密度擾動區(qū)域,紅色虛線標注了沖擊波位置Fig.2.Contour of vorticity at 1 ns (a) and 3 ns (b) of case SPI1.Blue circled part is the density perturbation region,and the red dashed line indicates the position of the shock.

密度擾動的類RM 不穩(wěn)定性并不存在明確的界面,無法像RM 不穩(wěn)定性一樣找到確定的界面來統(tǒng)計擾動幅值.圖2 的渦量場分布表明可以使用流場中最大渦量的渦對寬度來表征密度擾動的寬度.選取ωc=ηωωmax為臨界值(ωc應(yīng)大于由擾動沖擊波產(chǎn)生的渦量最大值),統(tǒng)計渦量大于臨界值的區(qū)域?qū)挾茸鳛闇u對的寬度.沖擊波過后,密度擾動區(qū)域獲得擾動速度,沖擊波作用產(chǎn)生的渦對發(fā)生旋轉(zhuǎn),密度擾動區(qū)域變寬.顯然,擾動速度越大密度擾動的寬度增大越快,渦對的寬度增加越快.因此,可以使用密度擾動區(qū)域渦對的寬度增長來表征類RM 不穩(wěn)定性的增長.另外,沖擊波與密度擾動作用產(chǎn)生的渦量與y方向的切向速度密切相關(guān)[17?21,28],也可以使用y方向擾動速度的最大值來表征不穩(wěn)定性的增長.圖3 給出了算例SPI1 的渦對寬度D以及y方向最大擾動速度隨時間的演化.本文主要研究線性階段的類RM 不穩(wěn)定性,可以看到擾動寬度隨時間線性增大,擾動速度隨時間的演化與RM 不穩(wěn)定性類似,在線性值附近振蕩.圖3(a)中實線由1.5 ns 后的數(shù)據(jù)進行線性擬合得到的.盡管ηω分別為0.2,0.4,0.6 和0.8 時擬合出的增長速度不同,但均可反應(yīng)線性增長規(guī)律.本文后續(xù)的算例選取ηω=0.2 來統(tǒng)計密度擾動的寬度.

圖3 算例SPI1 的渦對寬度(a)和y 方向最大擾動速度(b)隨時間的變化Fig.3.Time histories of the width of the vortex pair (a) and the maximum tangential velocity (b) of case SPI1.

對不同密度擾動大小η,沖擊波馬赫數(shù)Ma,以及不同擾動波長λy的沖擊波與密度擾動作用問題進行數(shù)值模擬,研究這些物理參數(shù)對類RM 不穩(wěn)定性增長的影響,不同算例的參數(shù)設(shè)置在表1 列出.算例SPI1-SPI4 擾動波長不同,算例SPI5-SPI8的密度擾動大小不同,算例SPI9-SPI12 的沖擊波馬赫數(shù)不同.其中M0=9.5,為程序計算出的激光功率密度為10 TW/cm2時的平臺脈沖產(chǎn)生的沖擊波馬赫數(shù).

通過對數(shù)值模擬結(jié)果的分析,發(fā)現(xiàn)類RM 不穩(wěn)定性的增長速度隨ky,Δu以及η線性增大,即δv ∝kyΔuη,如圖4 所示.其中Δu為沖擊波后的流體速度減沖擊波前的流體速度.圖4 中的點是由數(shù)值模擬結(jié)果統(tǒng)計出的密度擾動寬度的增大速度,直線為數(shù)值模擬結(jié)果的線性擬合.密度擾動的類RM 增長是由激波與密度擾動作用時累積的渦量決定的.擾動波長越長,在沖擊波與密度擾動相互作用時密度和壓強的斜壓角度越小,沉積的渦量越小.當η減小時,密度梯度變小,同樣的會減少密度擾動區(qū)域的渦量沉積.沖擊波馬赫數(shù)減小時,減小了斜壓效應(yīng)中的壓強梯度,也會減少渦量沉積.在強沖擊的情況下,波后流體速度正比于沖擊波馬赫數(shù),因此也可以寫為δv∝kyvsη,vs為沖擊波速度.類RM 不穩(wěn)定性的增長速度與RM 不穩(wěn)定性沖擊模型的形式類似[15],只是沖擊模型中的Atwood 數(shù)替換為η.

圖4 密度擾動寬度增長速度δv 隨kyΔuη 的變化Fig.4.Curve of density perturbation width growth rate δv versus kyΔuη.

4 密度擾動與無擾動界面的耦合

在ICF 內(nèi)爆中早期階段,靶內(nèi)密度擾動引起的類RM 不穩(wěn)定性耦合到燒蝕面,產(chǎn)生加速階段的擾動種子.密度擾動存在于靶丸內(nèi)部材料(高密度物質(zhì))中,燒蝕面外為燒蝕出的等離子體(低密度物質(zhì)),燒蝕面上的Atwood 數(shù)是大于0 的.擾動種子的形成與流體不穩(wěn)定性耦合、渦量輸運以及質(zhì)量燒蝕效應(yīng)等因素相關(guān)[13,29,30].本節(jié)研究了無質(zhì)量燒蝕、熱傳導(dǎo)的情況下類RM 不穩(wěn)定性與界面的耦合.沖擊波后有一個無擾動界面,界面兩側(cè)流體速度相同,壓強相同,密度不同ρ2＞ρ1,與沖擊波后流體一起運動.本節(jié)選取λx=2 μm.圖5 展示了密度擾動到界面耦合的密度云圖,計算參數(shù)Ma=9.5,η=0.3,λy=20 μm,At=0.77.沖擊波壓縮后界面到密度擾動區(qū)域中心的距離L=2.35 μm.左側(cè)界面為初始無擾動界面,右側(cè)界面為沖擊波陣面.密度擾動的類RM 不穩(wěn)定性耦合到無擾動界面引起界面的不穩(wěn)定性增長.

圖5 類RM 不穩(wěn)定性與界面耦合問題的密度云圖,此算例中的物理參數(shù)為Ma=9.5,η=0.3,λy =20 μm,At=0.77Fig.5.Contour of density of the RM-like instability and interface coupling problem,the physical parameters in this case are Ma=9.5,η=0.3,λy =20 μm,At=0.77.

擾動增長速度可以用y方向的擾動速度來表征,為了說明密度擾動的類RM 不穩(wěn)定性到界面的耦合機制,圖6 給出了不同L,同一時刻的y方向擾動速度分布.界面位置在圖中用黑色虛線標出.圖6(a)界面右側(cè)第一個渦對區(qū)域為密度擾動區(qū)域,圖中界面距密度擾動區(qū)域較遠,熵渦波并未傳播到界面位置,這種情況下不穩(wěn)定性通過聲波耦合到界面,沖擊波與密度擾動作用產(chǎn)生的聲波會在界面上沉積渦量.聲波擾動的壓強和密度滿足,其中cs為聲速.因此在均勻流體中聲波擾動產(chǎn)生的密度梯度和壓強梯度是同向的.當聲波傳播到界面位置,密度梯度主要由界面處的密度分布決定(擾動聲波帶來的密度擾動是小量),此時密度梯度和壓強梯度不再同向,產(chǎn)生斜壓效應(yīng)從而累積渦量.圖6(b)界面上擾動速度大于圖6(a),表明界面距離密度擾動越近,聲波在界面上引起的不穩(wěn)定性增長越大.另外,聲波在界面上發(fā)生反射向右傳播,還會增強密度擾動和沖擊波之間的渦量場,這點可以從圖6 中密度擾動區(qū)域右側(cè)的擾動速度對比得出.圖6(c)中界面與密度擾動距離較近,密度擾動很快傳播到界面位置,密度擾動上的渦和界面上的渦會發(fā)生渦合并.渦合并是密度擾動與界面耦合的第2 個機制.界面上的渦量和密度擾動的渦量方向相同,渦合并時渦量增強,會增大界面上的擾動速度.

圖6 y方向擾動速度云圖(a) L=5.53 μm;(b) L=3.53 μm;(c) L=1.53 μm.物理參數(shù)Ma=9.5,η=0.05,λy =20 μm,At =0.77Fig.6.Contour of the y-component of the perturbation velocity: (a) L=5.53 μm;(b) L=3.53 μm;(c) L=1.53 μm.Physical parameters: Ma=9.5,η=0.05,λy =20 μm,At =0.77.

數(shù)值模擬結(jié)果表明,密度擾動的類RM 不穩(wěn)定性到界面的耦合與距離L是相關(guān)的.圖7 為不同L時,界面擾動幅值隨時間的演化.考慮20 μm和40 μm 兩個波長,其余物理參數(shù)與圖6 相同.界面早期的擾動增長是由于氣泡和尖釘處受到壓強相反的聲波影響,運動方向相反,會產(chǎn)生一個很陡峭的增長.界面上擾動開始增長的時間與聲波傳播到該距離的時間是一致的.L越小,界面擾動增長越快,隨著L增大界面擾動增長速度逐漸減小,kyL=3.3 時,20 μm 波長的界面擾動幾乎無增長.將圖7 中算例界面擾動線性增長的一段進行擬合,得到界面增長速度δvi.前面的研究發(fā)現(xiàn),密度擾動的增長速度正比于kyΔuη,使用該值對界面增長速度進行歸一化.RM 不穩(wěn)定界面到無擾動界面的耦合系數(shù)與界面間距離L的關(guān)系通常以的形式出現(xiàn)[31,32],在解析上目前還沒有較好的方法來分析密度擾動到界面的耦合,在分析L對密度擾動與界面耦合的影響時,也假定其是以的形式存在的.圖8 給出了歸一化后的界面增長速度δvi/(kyΔuη)隨的變化.數(shù)值模擬結(jié)果表明,當kyL比較大時,δvi/(kyΔuη)∝,即密度擾動到界面的耦合隨著kyL增加以e 指數(shù)的形式衰減.當密度擾動與界面距離比較近的時候,擾動速度得到增強,界面擾動增長速度偏離與的線性關(guān)系.

圖7 不同L 時,界面擾動幅值隨時間的變化 (a) λy=20 μm;(b) λy=40 μm.Fig.7.Variation of interface perturbation amplitudes with time for different L: (a) λy=20 μm;(b) λy=40 μm.

圖8 界面擾動增長速度δvi/(kyΔuη)隨 的變化Fig.8.Curve of the interface disturbance growth rate δvi/(kyΔuη) versus .

界面的Atwood 數(shù)同樣會影響界面的擾動速度,圖9 給出了不同Atwood 數(shù)時界面擾動幅值的演化.具體物理參數(shù)為Ma=9.5,L=5.5 μm,λy=20 μm,η=0.05.數(shù)值模擬結(jié)果表明,界面Atwood 數(shù)越大,界面的擾動速度越大.Atwood 數(shù)增大會增大界面上的密度梯度,從而增加反射聲波在界面上的渦量沉積,Atwood 數(shù)大于0 時,界面上的擾動速度與密度擾動的擾動速度方向相同,因而耦合在一起時是互相增強的,會影響界面上密度梯度及界面的過渡層.考慮e 指數(shù)形式的過渡層,密度分布為

圖9 不同界面Atwood 數(shù)時,界面擾動幅值隨時間的變化Fig.9.Time evolution of interface perturbation amplitude at different Atwood numbers.

圖10 對比了不同Lm時,界面擾動幅值隨時間的變化.Lm越大,密度過渡層越寬.物理參數(shù)Ma=9.5,L=5.5,λy=20 μm,η=0.1,At=0.77.隨著密度過渡層變寬,界面擾動增長得到了抑制.e 指數(shù)形式的密度過渡層對不穩(wěn)定性的影響可通過在Atwood 數(shù)引入一個衰減因子來描述,可以將其整體看作等效Atwood 數(shù),燒蝕面的等效Atwood數(shù) 為=At/(1+kLm)=1/(1+kLm)[5],的降低可以減小密度擾動類RM 不穩(wěn)定性到界面的耦合.

圖10 不同界面過渡層寬度時,界面擾動幅值隨時間的變化Fig.10.Time evolution of interface perturbation amplitude with different density transition layer.

5 總結(jié)與討論

本文研究了沖擊波與密度擾動作用產(chǎn)生的類RM 不穩(wěn)定性增長規(guī)律以及不穩(wěn)定性到界面的耦合規(guī)律,獲得了一些新的物理認識.類RM 不穩(wěn)定性的產(chǎn)生機制是沖擊波與密度擾動作用時,由于斜壓效應(yīng)產(chǎn)生的渦量沉積.我們發(fā)現(xiàn),類RM 不穩(wěn)定性線性階段的增長速度δv∝kyΔuη.密度擾動的類RM 不穩(wěn)定性到界面的耦合有聲波耦合和渦合并兩種機制,當密度擾動距離界面較遠時只有聲波可耦合到界面,當密度擾動距離界面比較近時會發(fā)生渦合并.聲波耦合引起的界面擾動增長速度.界面上的Atwood 數(shù)為正時,界面上渦量和密度擾動的渦量方向相同,渦合并導(dǎo)致擾動速度增大.降低界面上的Atwood 數(shù)或增大界面上的密度標長可以減小密度擾動的類RM 不穩(wěn)定性到界面的耦合.

充分理解材料內(nèi)部密度擾動類RM 不穩(wěn)定性的增長規(guī)律及不穩(wěn)定性到界面的耦合機制,對理解ICF 內(nèi)爆中密度擾動對燒蝕面上不穩(wěn)定性及對內(nèi)爆性能的影響是重要的.關(guān)于界面Atwood 數(shù)對擾動到界面耦合的影響目前只是有一個定性的認識,仍需繼續(xù)開展定量的研究.靶內(nèi)擾動非線性的演化規(guī)律認識還不足.內(nèi)爆早期會產(chǎn)生沖擊波、稀疏波以及壓縮波,不同的波與密度擾動作用產(chǎn)生的不穩(wěn)定性及不穩(wěn)定性到界面的耦合也有待研究.