密度扰动的类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 数对扰动到界面耦合的影响目前只是有一个定性的认识,仍需继续开展定量的研究.靶内扰动非线性的演化规律认识还不足.内爆早期会产生冲击波、稀疏波以及压缩波,不同的波与密度扰动作用产生的不稳定性及不稳定性到界面的耦合也有待研究.