殷建伟 潘昊 吴子辉 郝鹏程 胡晓棉

1)(北京理工大学机电学院,北京 100081)

2)(北京应用物理与计算数学研究所,计算物理重点实验室,北京 100094)

3)(中国工程物理研究院研究生部,北京 100088)

1 引 言

冲击波加载不同密度的流体界面时,界面扰动将增长演化为涡结构并最终形成湍流混合层,这种流体动力学现象称为Richtymer-Meshkov(RM)流动不稳定性[1−4],大量研究表明,流体的表面张力、黏性等耗散性质将削弱RM不稳定性的增长.对于固体介质,其RM流动过程受到固体物性的显著影响,呈现出不同的物理现象:在冲击波作用下,材料的强度致稳效应将导致扰动振幅增长一定程度后围绕某均值周期性振荡[5−9],若冲击波压力极高,如惯性约束聚变(ICF)内爆中加载靶丸的压力达到100 GPa量级[10],扰动将剧烈增长,伴随升温、破碎、喷射、多物质混合等复杂物理现象,是武器内爆动力学、物质界面微喷射与混合动力学过程重点关注的问题.

目前对固体中RM流动的理论研究基于弹塑性材料均匀性假设,未考虑材料的相变、断裂等复杂物理因素.Piriz[8]提出了刚塑性材料扰动面上尖钉振幅的极大值表达式,并在Buttler等[11,12]和Dimonte等[13]的实验中得到初步验证.Jensen等[14]在平面撞击实验中利用该表达式推算了铈在极端加载条件下的强度,结果表明不同方法得到的强度结果偏差约20%—50%.因此,如何改进极值振幅表达式,更准确地描述弹塑性固体RM流动的扰动增长规律是目前亟待解决的问题,也是理解强度致稳机制的核心问题,有助于提高对界面失稳早期阶段的物理认识.

本文采用AFE2D程序对炸药爆轰加载下金属材料的界面扰动增长过程进行数值模拟研究,定性分析了初始扰动的振幅与波长对扰动增长的影响以及材料强度的致稳效应,利用相关冲击波关系式及流体RM不稳定性的研究成果,进一步推导得到了极值振幅的理论关系式,该式能够描述扰动增长被抑制时振幅与初始扰动、材料性质之间的关系,通过分析关系式中各项的相关性,初步得到了线性关系式的适用范围.

2 界面扰动增长问题的数值模拟

考虑如图1所示的界面扰动问题,高能炸药爆轰驱动金属样品,该问题是Buttler等[11]开展的铜、锡等材料RM实验的简化模型,其中高能炸药厚度12 mm,金属样品厚6 mm,样品中的冲击波向介质/真空界面传播时加载扰动自由面,扰动面为正弦形,其初始波长为λ0,初始振幅为ξ0.

图1 冲击波加载界面扰动问题示意图Fig.1.Illustration of the surface perturbation problem loaded by the shock wave.

炸药爆轰采用时间起爆方式,状态方程为JWL类型状态方程[15]

其中R1,R2和ω是JWL模型参数;A和B是JWL模型系数,由高能炸药的CJ爆速DCJ、CJ爆压PCJ及JWL参数确定;V=ρ/ρ0是相对比容,ρ0是初始密度,ρ是材料密度;E是单位初始体积的内能.固体材料的状态方程为Mie-Grüneisen状态方程[16]

其中µ=ρ/ρ0−1;γ0为Grüneisen系数;b为Grüneisen系数的一阶体积修正量;c0,S1,S2和S3为材料冲击波-粒子速度关系式的系数.固体材料本构采用Prandtl-Reuss关联塑性流动律和von Mises屈服准则描述[7−9],即

下标采用Einstein求和约定,Dij为应变率张量;Sij为偏应力张量;G为材料剪切模量,Y为屈服强度,采用弹性理想塑性模型,保持G和Y不变.几种常见金属材料的状态方程参数与剪切模量见表1,炸药的爆轰参数与JWL状态方程参数见表2.

利用AFE2D程序数值模拟PBX9501炸药爆轰驱动铜飞片的实验(见文献[11]Cu-pRad0426实验),如图1所示,样品的自由面上引入kξ0(k=2π/λ0,初始扰动的波数)分别为0.12,0.35,0.75和1.5的初始扰动,初始波长均为550µm,铜样品中的冲击波压力约36 GPa.图2为实验的数值模拟结果,可以看到冲击波加载扰动面时,初始波谷位置首先启动,形成向外的突起尖钉(Spike)结构;初始波峰位置稍后启动,与尖钉等部位间的速度差导致其向材料内部凹陷,形成气泡(Bubble)结构;反相后尖钉与气泡结构进一步演化,初始扰动kξ0=0.12时尖钉基本不增长,kξ0=0.35时尖钉增长一段时间后被抑制,kξ0=0.75和1.5时增长不受抑制,尖钉演化为射流.

表1 金属材料的状态方程参数和剪切模量Table 1.Parameters of equation of state and shear modulus of several selected metals.

表2 PBX9501炸药的爆轰参数与JWL状态方程参数Table 2.Parameters of the detonation performance and JWL equation of state of the high explosives PBX9501.

定义尖钉的振幅ξSP为尖钉与无扰动区域(Bulk)之间在冲击波传播方向的相对位移,尖钉的增长速度vSP为尖钉与无扰动区域之间在冲击波传播方向的相对速度,即

气泡的振幅ξBU和增长速度vBU的定义类似,

模拟结果显示,kξ0=0.12时尖钉的最大振幅=max(ξSP)=10.55 µm,相比初始振幅=11µm几乎没有增长,但相位出现反转;kξ0=0.35时=155.3 µm. 在Buttler对CupRad0426实验数据的分析中[11],kξ0=0.12时=(7±4)µm,kξ0=0.35时=160 µm.表3对比了不同kξ0时尖钉增长速度的极大值=max(vSP),模拟结果与实验数据基本一致,表明本文采用的程序、模型和参数适用于界面扰动问题,可用于分析界面扰动的增长规律.

表3 PBX9501炸药驱动铜飞片实验的模拟结果Table 3.Simulation results of the Cu-pRad0426 experiment where the copper fl yer plate was driven by the high explosive PBX9501.

图2 Cu-pRad0426实验的模拟结果 (a)初始扰动界面;(b)冲击波首次加载界面后0.7µs时刻的模拟结果;(c)3.7µs时刻的模拟结果;(d)3.7µs时刻实验质子照相图像[11]Fig.2.The simulation results of the Cu-pRad0426 experiment:(a)Initial state;(b)plot of the simulation results at 0.7µs which relative to shockwave breakout at the free surface;(c)plot of the simulation results at 3.7µs;(d)proton radiographic picture of the experiment at 3.7µs[11].

3 界面扰动的增长规律

3.1 初始形态对扰动增长的影响

首先分析扰动初始形态对后期增长的影响,利用无量纲数kξ0表征初始扰动的振幅与波长之比,考虑图1所示的界面扰动问题,选取初始波长λ0=450—600 µm,初始振幅ξ0=15—38 µm,则kξ0的变化范围0.20—0.40,样品材料为OFHC铜,屈服强度Y=0.47 GPa,炸药为PBX9501.图3比较了冲击波首次加载界面后3.1µs时刻不同kξ0的扰动增长情况,尖钉的振幅随着kξ0增加逐步变大,尖钉形状在kξ0=0.20时仍接近正弦扰动形态,继续提高kξ0尖钉将逐渐拉伸、变窄,最终演化为射流.上述结果表明,对于相同的加载条件和样品材料,若初始扰动的振幅与波长比值越小,相应的初始界面越“平坦”,扰动增长越不明显,界面保持稳定;振幅与波长比值越大,相应的初始界面越“沟壑”,扰动越易于增长,甚至形成局部的微射流,界面失稳.

图3 不同kξ0值的界面扰动增长情况,冲击波首次加载界面后3.1µs时刻Fig.3.Growth results of the perturbated surfaces with different values of kξ0,3.1 µs relative to shockwave breakout at the free surface.

3.2 材料强度对扰动增长的致稳效应

相关研究表明,金属材料的动态屈服强度存在较强的应变硬化效应,并且在高应变率下表现出显著的应变率硬化效应[17−21].本文所考虑的冲击波加载扰动面问题中,扰动面上形成的尖钉和气泡等特征结构均为局部高应变区,局部应变率达到106—107s−1.然而目前高应变率范围内金属材料的强度数据极少,对应变率的硬化规律认识亦有待进一步完善,本文所采用的弹性理想塑性模型虽不能反映材料真实的应变与应变率硬化效应,但在均匀性假设下通过调节屈服强度模拟不同的硬化效果也不失为一种可行的近似方法,能够用于近似阐述材料强度的致稳效应.

选择扰动的kξ0=0.251,取OFHC铜的屈服强度为Y=0.17 GPa,令其逐步增加至0.57 GPa以模拟不同的硬化效应.图4为扰动界面的增长情况,强度较低时扰动面的尖钉演化为射流,提高材料强度后尖钉增长受到抑制,逐渐从尖锐的射流形态退化为类似正弦波形的扰动形态,定性地显示出材料强度对扰动增长的致稳效应,表明在相同的加载和初始扰动条件下,强度越高的材料扰动增长越不明显.

3.3 扰动增长规律的线性近似

由图3和图4可知,提高初始扰动的kξ0,扰动面的尖钉振幅将逐渐增大,在爆轰加载下扰动增长受到材料强度的显著抑制作用,因此当扰动增长被抑制时,尖钉增长幅度的控制因素应包含初始扰动kξ0与强度Y的某种组合形式.

图4 扰动增长情况,调节屈服强度以模拟不同程度的硬化效应,冲击波首次加载界面后3.1µs时刻Fig.4. Growth results of the surfaces with varied yield strength for the different equivalent effects of hardening,3.1µs relative to shockwave breakout at the free surface.

利用(4)式定义的尖钉振幅、增长速度和Dimonte等[7,13]的分析方法,类似于尖钉最大振幅和最大增长速度的定义,得到气泡最大振幅=max(ξBU)和气泡最大增长速度=max(vBU),利用初始波数无量纲化尖钉和气泡的最大振幅

利用材料初始密度ρ0、屈服强度Y无量纲化尖钉和气泡的最大增长速度,分别标记为无量纲数和

图5 (网刊彩色)尖钉和气泡的最大振幅与最大增长速度的关系(铜强度分别取0.27,0.32和0.37 GPa)Fig.5. (color online)Relations between the nondimensional maximum amplitude and growth velocity of the spikes and bubbles(yield strength of copper varied as 0.27,0.32 and 0.37 GPa).

图5和图6分别显示了铜、铝样品中尖钉和气泡的最大振幅与最大增长速度无量纲数之间的关系,选择初始扰动kξ0变化范围0.15—0.3以避免尖钉演化为射流,屈服强度分别取三组值.模拟结果显示,无量纲化的尖钉最大振幅与最大增长速度在一定范围内表现出近似线性关系,斜率C约为0.3,即

图6 (网刊彩色)尖钉和气泡的最大振幅与最大增长速度关系(铝强度分别取0.24,0.29和0.34 GPa)Fig.6. (color online)Relations between the nondimensional maximum amplitude and growth velocity of the spikes and bubbles(yield strength of aluminum varied as 0.24,0.29 and 0.34 GPa).

Piriz理论预测的斜率值为0.29[7],Dimonte利用PAGOSA程序拟合的数据点见图5[13],忽略截距后斜率值为0.24,本文结果中线性拟合通过原点,斜率约为0.3,更接近Piriz的理论预测值.

随着扰动kξ0增大,尖钉的增长速度逐步增加,最大振幅与增长速度之间表现出非线性关系,尖钉形态从扰动增长被抑制时的稳定态向持续增长的不稳定态过渡,直至尖钉头部收缩变窄形成射流,最终脱离样品主体形成喷射物.流体RM不稳定性分析研究表明[1,2,10,22],在线性近似假设下尖钉的最大增长速度可近似表示为初始扰动kξ0、冲击波经过后的界面起跳速度Δv和界面两侧Atwood数的乘积,即

其中Atwood数A=(ρ2−ρ1)/(ρ2+ρ1)表示界面两侧物质的密度差,对于固体/真空界面,A=−1.冲击波在弹塑性材料中传播时,界面起跳速度Δv可由冲击波强度和材料的Hugoniot冲击线确定[11,17,23−25],若介质的初始压力和速度为零,压力为P的冲击波传播时波后粒子速度与压力的关系为

其中η表示冲击波前后介质的压缩率,η=1−ρ0/ρ,与材料的声阻抗和冲击波速度有关;vp为波后粒子速度.对于介质/真空界面,自由面速度与波后粒子速度存在倍数关系,Δv≈2vp,将(9)和(10)式代入(8)式中得

因此尖钉振幅的增长因子为

(12)式表明,扰动增长被抑制时尖钉的最大振幅增长因子与参数组合kξ0/Y成线性关系,并且冲击波压力越高,尖钉增长越显著.

3.4 线性近似的适用范围

分析(12)式中制约尖钉增长因子的各因素:系数C是线性关系的斜率;压缩率η与冲击波压力P共同表征了冲击波的压缩效应,为扰动增长提供启动能量;kξ0描述了扰动界面的初始形态,为初始扰动的振幅与波长之比;屈服强度Y是材料的物性参数,扰动增长受其抑制.因此后三个因素之间是相互独立的,本文分析系数C与其他因素之间的相关性,讨论线性近似的适用范围.

考虑三种材料,分别是6061铝、OFHC铜和304不锈钢,状态方程为Mie-Grüneisen状态方程,本构模型为弹性理想塑性模型,参数见表1,炸药为PBX9501,相关模型参数见表2.每种材料的屈服强度取三组值以模拟不同程度的硬化效应,初始扰动kξ0的变化范围0.15—0.35.图7为不同材料、强度和kξ0时系数C随参数组合kξ0/Y的变化情况,图8为kξ0=0.223和0.279时铝材料扰动面上尖钉、气泡和无扰动区域的速度曲线,分别对应扰动增长被抑制和未被抑制的情况.

图7 (网刊彩色)无量纲的尖钉最大振幅与增长速度之比随参数组合kξ0/Y 的变化情况Fig.7.(color online)Plots of the ratios between the non-dimensional maximum amplitude and growth velocity of the spikes verse the combination of parameters kξ0/Y.

由图7可知,系数分布存在三个阶段,对应于界面扰动的稳定、稳定与不稳定之间的过渡和不稳定阶段.

1)参数组合kξ0/Y≤0.8 GPa−1,系数C基本为常数,扰动处于被抑制的稳定态,尖钉的最大振幅与增长速度之间为线性关系,增长因子由(12)式描述,并且图7中不同材料的系数一致,说明在稳定段系数C基本保持独立性,不受其它因素的影响.

2)参数组合0.8 GPa−1≤kξ0/Y≤1.0 GPa−1,系数C不再为常数,扰动增长情况由稳定态向不稳定态过渡.此阶段扰动可能会形成稳定的尖钉形态,也可能进一步发展成射流,取决于两个特征时间的竞争,即冲击波首次加载后在样品内部波系多次反射形成的二次冲击波加载界面所需时间twave,与冲击波首次加载后尖钉振幅从开始增长至被抑制所需的稳定时间tstable.若twave＞tstable,表明在二次冲击波到达前,强度的致稳效应已抑制了扰动增长,尖钉、气泡和未扰动区域的速度趋于一致,如图8所示,从被抑制到二次冲击波到达,这段时间内尖钉、气泡和无扰动区域之间的相对速度保持为零,意味着尖钉、气泡的振幅不再变化,界面形状保持不变;若twave≤tstable,表明在二次冲击波到达自由面时,强度的致稳效应未能将尖钉充分减速,尖钉与其他位置相比存在速度差,冲击波二次加载将破坏此时的尖钉和气泡形态,界面形状进一步演化.在过渡阶段,(8)式的线性近似关系与(12)式的尖钉增长规律依然可能成立,取决于扰动增长在二次冲击波达到前是否能被抑制.

图8 (网刊彩色)扰动增长稳定与不稳定时尖钉、气泡和未扰动区域的速度曲线Fig.8.(color online)Velocity pro files of the spikes,bubbles and the unperturbated areas in the stable and unstable growth of the perturbations.

3)参数组合kξ0/Y≥1.0 GPa−1,扰动增长显著,尖钉头部迅速演化为射流,界面失稳.

综上所述,在扰动增长的稳定态,系数C保持独立性,线性关系(8)式适用,描述尖钉增长规律的线性近似(12)式具备一定的材料普适性;在稳定态向不稳定态过渡阶段以及不稳定态,系数C受其他因素影响失去独立性,(8)式和(12)式不具备材料普适性;图7中系数C保持为常数的平台区域揭示了界面扰动从增长被抑制的线性关系区向形成射流的非线性区过渡时存在稳定性阈值,表明强度介质的RM扰动增长中也存在类似Rayleigh-Taylor扰动增长的稳定性阈值[26],初步研究结果表明阈值对应的物理状态kξ0/Y取值约0.8 GPa−1.

4 结 论

本文研究了爆轰加载下不同材料飞片RM流动的扰动增长规律,分析结果表明扰动增长被抑制时,尖钉振幅增长因子与初始扰动的kξ0成正比,与材料的强度Y成反比,理论关系式符合相关实验结果,与Piriz利用牛顿第二定律近似推导所得的结果基本一致,该理论关系式也为测量材料动态屈服强度提供了新的实验技术途径.对扰动增长规律影响因素的相关性分析表明,不同金属材料的扰动增长均存在稳定区域,稳定区内描述尖钉振幅增长因子与加载条件、初始扰动和材料性质之间关系的(12)式具备一定的材料普适性,稳定区范围kξ0/Y＜0.8 GPa−1,其中阈值0.8 GPa−1的物理意义有待进一步分析研究.

[1]Richtmyer R D 1960Commun.Pure Appl.Math.13297

[2]Meshkov E E 1969Sovit.Fluid Dyn.4151

[3]Mikaelian K O 2013Phys.Rev.E87031003

[4]Brouillette M 2002Annu.Rev.Fluid Mech.34445

[5]Piriz A R,Lopez Cela J J,Tahir N A,Hoff mann D H H 2006Phys.Rev.E74037301

[6]Piriz A R,Lopez Cela J J,Cortazar O D,Tahir N A,Hoff mann D H H 2005Phys.Rev.E72056313

[7]Piriz A R,Lopez Cela J J,Tahir N A,Hoff mann D H H 2008Phys.Rev.E78056401

[8]Piriz A R,Lopez Cela J J,Tahir N A 2009Phys.Rev.E80046305

[9]Lopez Ortega A,Lombardini M,Pullin D I,Meiron D I 2014Phys.Rev.E89033018

[10]Remington B A,Rudd E R,Wark J S 2015Phys.Plasmas22090501

[11]Buttler W T,Oro D M,Preston D L,Mikaelian K O,Cherne F J,Hixson R S,Mariam F G,Morris C,Stone J B,Terrones G,Tupa D 2012J.Fluid Mech.70360

[12]Buttler W T,Oro D M,Olsen R T,Cheren F J,Hammerberg J E,Hixson R S,Monfared S K,Pack C L,Rigg P A,Stone J B,Terrones G 2014J.Appl.Phys.116103519

[13]Dimonte G,Terrones G,Cheren F J,Germann T C,Dunpont V,Kadau K,Buttler W T,Oro D M,Morris C,Preston D L 2011Phys.Rev.Lett.107264502

[14]Jensen B J,Cheren F J,Prime M B,Fezzaa K,Iverson A J,Carlson C A,Yeager J D,Ramos K J,Hooks D E,Cooley J C,Dimonte G 2015J.Appl.Phys.118195903

[15]Sun Z F,Xu H,Li Q Z,Zhang C Y 2010Chin.J.High Pressure Phys.2455(in Chinese)[孙占峰,徐辉,李庆忠,张崇玉2010高压物理学报2455]

[16]Robinson A C,Swegle J W,1989J.Appl.Phys.662838

[17]Zhu J S,Hu X M,Wang P,Chen J,Xu A G 2010Adv.Mech.40400(in Chinese)[朱建士,胡晓棉,王裴,陈军,许爱国2010力学进展40400]

[18]Vogler T J,Chhabildas L C 2006Int.J.Impact Engng.33812

[19]Barton N R,Bernier J V,Becker R,Arsenlis A,Cavallo R,Marian J,Rhee M,Park H S,Remington B A,Olson R T 2011J.Appl.Phys.109073501

[20]Smith R F,Eggert J H,Rudd R E,Swift D C,Blome C A,Collins G W 2011J.Appl.Phys.110123515

[21]Park H S,Rudd R E,Cavallo R M,Barton N R,Arsenlis A,Belof J L,Blobaum K J M,El-dasher B S,Florando J N,Huntington C M,Maddox B R,May M J,Plechaty C,Prisbrey S T,Remington B A,Wallace R J,Wehrenberg C E,Wilson M J,Comley A J,Giraldex E,Nikroo A,Farrell M,Randall G,Gray III G T 2015Phys.Rev.Lett.114065502

[22]Wouchuk J G 2001Phys.Rev.E63056303

[23]Pan H,Wu Z H,Hu X M,Yang K 2013Chin.J.High Pressure Phys.27778(in Chinese)[潘昊,吴子辉,胡晓棉,杨堃2013高压物理学报27778]

[24]Pan H,Hu X M,Wu Z H,Dai C D,Wu Q 2012Acta Phys.Sin.61206401(in Chinese)[潘昊,胡晓棉,吴子辉,戴诚达,吴强2012物理学报61206401]

[25]Yu Y Y,Tan H,Hu J B,Dai C D,Chen D N,Wang H R 2008Acta Phys.Sin.572352(in Chinese)[俞宇颖,谭华,胡建波,戴诚达,陈大年,王焕然 2008物理学报 572352]

[26]Colvin J D,Legrand M,Remington B A,Schurtz G,Weber S V 2003J.Appl.Phys.935287