密度擾動的類Richtmyer-Meshkov 不穩定性增長及其與無擾動界面耦合的數值模擬*

孫貝貝 葉文華 張維巖

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

2) (北京應用物理與計算數學研究所,北京 100094)

1 引言

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

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

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

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

2 數值模擬初始設置

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

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

圖1 數值模擬初始設置示意圖Fig.1.Schematics of the initial configuration.

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

3 線性階段密度擾動的類RM 不穩定性

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

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

表1 不同算例的參數設置Table 1.Initial physical parameters in different cases.

圖2 算例SPI1 在1 ns (a)和3 ns (b)時的渦量場.藍色圈起部分為密度擾動區域,紅色虛線標注了沖擊波位置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 不穩定性并不存在明確的界面,無法像RM 不穩定性一樣找到確定的界面來統計擾動幅值.圖2 的渦量場分布表明可以使用流場中最大渦量的渦對寬度來表征密度擾動的寬度.選取ωc=ηωωmax為臨界值(ωc應大于由擾動沖擊波產生的渦量最大值),統計渦量大于臨界值的區域寬度作為渦對的寬度.沖擊波過后,密度擾動區域獲得擾動速度,沖擊波作用產生的渦對發生旋轉,密度擾動區域變寬.顯然,擾動速度越大密度擾動的寬度增大越快,渦對的寬度增加越快.因此,可以使用密度擾動區域渦對的寬度增長來表征類RM 不穩定性的增長.另外,沖擊波與密度擾動作用產生的渦量與y方向的切向速度密切相關[17?21,28],也可以使用y方向擾動速度的最大值來表征不穩定性的增長.圖3 給出了算例SPI1 的渦對寬度D以及y方向最大擾動速度隨時間的演化.本文主要研究線性階段的類RM 不穩定性,可以看到擾動寬度隨時間線性增大,擾動速度隨時間的演化與RM 不穩定性類似,在線性值附近振蕩.圖3(a)中實線由1.5 ns 后的數據進行線性擬合得到的.盡管ηω分別為0.2,0.4,0.6 和0.8 時擬合出的增長速度不同,但均可反應線性增長規律.本文后續的算例選取ηω=0.2 來統計密度擾動的寬度.

圖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.

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

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

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

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

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

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

圖6 y方向擾動速度云圖(a) L=5.53 μm;(b) L=3.53 μm;(c) L=1.53 μm.物理參數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.

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

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

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

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

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

5 總結與討論

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

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