内爆加载金属界面不稳定性的数值分析*

郝鹏程,冯其京,胡晓棉,2

(1.北京应用物理与计算数学研究所，北京 100094;2.北京应用物理与计算数学研究所计算物理国家重点实验室，北京 100088;3.中国工程物理研究院研究生部，北京 100088)

内爆加载金属界面不稳定性的数值分析*

郝鹏程1,2,3,冯其京1,胡晓棉1,2

(1.北京应用物理与计算数学研究所，北京 100094;2.北京应用物理与计算数学研究所计算物理国家重点实验室，北京 100088;3.中国工程物理研究院研究生部，北京 100088)

采用自行研制的多介质弹塑性流体力学欧拉程序，对柱形内爆加载金属界面不稳定性进行了数值研究，数值模拟结果与文献实验数据吻合较好。数值结果表明：材料强度对界面不稳定性发展有不可忽略的抑制作用；材料屈服强度对较高模数不稳定性增长的抑制较强，而剪切模量对不稳定性发展的影响相似但敏感性相对较弱；金属界面不稳定性增长存在最不稳定模数，最不稳定模数随屈服强度增加而减小，并近似与屈服强度的对数呈线性关系；随着壳的厚度减小，扰动增长加快。

流体力学；界面不稳定性；弹塑性；内爆；材料参数

轻介质支撑或者加速重介质，介质界面会发生Rayleigh-Taylor不稳定性(RTI)，不稳定性诱发界面扰动增长，重介质形成尖钉侵入轻介质，轻介质气泡进入重介质；而当激波掠过介质界面，会发生Richtmyer-Meshkov不稳定性(RMI)，RMI的发展结果往往成为RTI发展的初始条件[1-2]。不稳定性现象广泛存在于天体物理、地球物理、惯性约束聚变(ICF)等领域。无黏流体的不稳定性增长最快，存在多种制稳机制，比如黏性、表面张力、密度梯度、材料强度、烧蚀等。对于弹塑性介质，J.W.Miles[3]给出了无限厚弹性介质RT不稳定性增长的截止波长；D.C.Drucker[4]给出了弹性介质RT不稳定性增长的截止振幅；A.I.Lebedev等[5]给出了有限厚弹性介质RT不稳定性发展截止扰动振幅与波长、厚度之间的关系。J.F.Barnes等[6]通过爆轰加载铝片实验，确认了金属界面RT不稳定性的发展存在截止波长和振幅。J.W.Swegle等[7]对弹塑性介质RTI进行了充分的理论分析以及数值模拟工作，认为弹塑性介质界面不稳定性发展不仅依赖于扰动波长、振幅、介质厚度，还与压力加载方式等因素相关。K.Mikaelian[8]提出黏性流体RMI的指数增长模型，J.D.Colvin等[9]将塑性介质应力与黏性流体类比，分析弹塑性介质的界面不稳定性，A.R.Piriz等[1-2]对弹塑性介质RTI以及RMI进行了有意义的理论解析探索。界面不稳定性是实验研究材料弹塑性参数及其动态响应特性的重要手段之一[10-11]。内爆过程中界面不稳定性是重要的研究方向，是影响ICF靶丸压缩、点火状态的重要因素之一。V.Frachet等[12]通过套筒实验，研究了柱形内爆加载下流体以及金属介质界面不稳定性，H.S.Park等[13]在更高压力、应变率下研究金了属界面不稳定性。刘军等[14]对文献[12]中的实验进行了数值模拟。

本文中，使用自行研制的多介质弹塑性流体力学欧拉程序(MEPH2Y)[15]，对柱形内爆加载金属界面不稳定性进行数值模拟，分析影响金属界面不稳定性发展的因素，探索金属界面不稳定的发展规律。

1 计算方法

1.1 计算模型

图1 柱形内爆计算模型Fig.1 Computational model in cylindrical implosion

采用与文献[12]实验相似的几何构型(如图1所示)，外层为钢壳，向中心分别为高能炸药、薄钢壳、硅橡胶。钢外壳，厚度2 cm；炸药，外壁半径10 cm，厚5 cm，从外壁面起爆爆炸药；内层钢壳，厚4 mm，在内层钢壳的内壁刻振幅0.5 mm的正弦曲线。

1.2 数值方法

在界面不稳定性问题中存在界面大变形现象，因此，在数值模拟中采用自行研制的欧拉程序。所使用的多介质弹塑性流体力学程序(MEPH2Y)采用交替方向方法进行时间推进，在每个方向上进行拉氏计算，而后重映到原始网格，界面采用VOF方法重构界面。该流体力学程序能够模拟流体和弹塑性介质，能够用于炸药爆轰、射流射弹形成及其侵彻过程等多介质复杂流场模拟，具有较强的解决实际工程计算问题的应用能力。

1.3 状态方程及本构模型

金属采用凝聚介质实用状态方程：

式中：p为压力，ρ为密度，c0为初始声速，ρ0为初始密度，γ为比热比。材料参数分别为：钢材料，ρ0=7.8 g/cm3，c0=3.8 km/s，γ=4.075；硅橡胶，ρ0=1.1 g/cm3，c0=0.5 km/s，γ=1.6。爆轰计算中，炸药使用CJ爆轰结合Wilkins反应率的时间起爆方式，爆轰产物采用理想气体状态方程。炸药参数为：初始密度1.8 g/cm3，爆速7.5 km/s，比热比3.0。

材料采用理想弹塑性本构模型，使用von Mises准则确定介质的弹塑性状态，材料剪切模量G为常数，材料屈服强度Y满足：Y=max(0,Y0(1-e/em))(Y0为初始屈服强度，em为融化能)。该本构模型中，剪切模量与状态无关，屈服强度仅与内能相关，即仅有软化的表现。材料本构参数分别为：钢壳，G=83 GPa，Y0=2 GPa，em=1.3 kJ/g；硅橡胶基本为流体模型，Y0=0.01 GPa，G=0.01 GPa，em=0.1 kJ/g；炸药为流体。

2 结果及分析

2.1 一维数值模拟结果

图2 金属壳外壁压力Fig.2 Pressure on outer interface of metal shell

图3 金属壳内壁速度Fig.3 Velocity on inner interface of metal shell

通过对柱坐标一维情况的数值分析，能够认识无扰动情况下界面运动及受力情况，为进一步分析界面不稳定性发展提供参考。数值模拟中，采用0.1 mm均匀网格。图2为钢壳外壁的压力，图3为金属壳内的界面速度，时间原点为爆轰波到达钢壳外壁的时刻。爆轰波在金属外壁的压力峰值约为45 GPa，波在金属壳内往返，能够清晰地看到6 μs内有3次冲击波作用，壳层内界面得到加速，界面扰动将发生RM不稳定性。而后直到中心反射的冲击波到达前，界面几乎为匀减速运动，加速度指向较重的钢壳，界面扰动的不稳定性增长将以RT不稳定性为主。

2.2 二维数值模拟结果

在二维模拟中，金属壳的初始扰动模数n取13和29。图4为模拟得到的扰动模数为29的壳在某时刻的界面，从中可以定性地看到，材料强度对界面不稳定性发展有较大影响。

图5为不同初始扰动的扰动增长结果，包括纯流体以及弹塑性计算。在纯流体计算中，29模数的比13模数的扰动发展快，这与流体的扰动发展理论相符。对于相同模数扰动的界面，考虑弹塑性计算中的扰动比纯流体计算中的扰动发展得慢，印证了材料强度是不稳定性发展的致稳因素。在弹塑性计算中，29模数的扰动未得到发展，而13模数的扰动发展虽然比纯流体计算慢，但界面扰动仍然得到一定发展。这个弹塑性计算结果与实验数据定性相符，即流体计算中发展较快的高模数在弹塑性计算中反而发展得较慢，表明金属材料的不稳定性发展与扰动模数及材料强度均相关。另外，弹塑性计算的结果与实验值有一定的差异，应该与计算中所使用的状态方程、本构关系等材料特性的差异相关。

图4 密度等值云图(n=29)Fig.4 Density contour plot

图5 金属界面扰动增长Fig.5 Growth of metal interfacial perturbation amplitude

2.3 材料强度与不稳定性的发展

图6 材料屈服强度对界面扰动增长的影响Fig.6 Effect of yield strength on perturbation growth

图7 材料剪切模量对界面扰动增长的影响Fig.7 Effect of shear strength on perturbation growth

固定材料剪切模量和屈服强度中的一个，单独调节另一个。为了排除其他因素的干扰，材料本构中屈服强度保持不变，即Y=Y0。图6给出模数13、29在不同屈服强度下的扰动发展情况，可以看到，材料强度抑制不稳定性发展，屈服强度越大扰动发展越慢。对不同模数的扰动，屈服强度的影响不同：对于13模扰动，材料强度对扰动发展影响相对较小；而对29模扰动发展影响显著，当Y0=2.0 GPa时，扰动甚至没有发展。而图7给出模数13当Y0=1.5 GPa时在不同剪切模量下的扰动发展情况，可以看出，与屈服强度相比较，剪切模量在10～1 000 GPa这样大范围内对不稳定性发展影响甚微。实际上，从材料屈服强度以及剪切模量在控制方程中的角色，可以理解他们对于材料所表现出的强度的贡献。屈服强度控制材料等效偏应力张量，而剪切模量则仅控制材料偏应力张量的增长速率，因此，在较高的应变率加载下，材料的屈服强度决定了材料所体现出来的强度。另外，在扰动发展初期的RM不稳定性发展阶段，内部压强较高，材料强度相对较小，扰动增长受到材料强度的影响较弱，材料强度的作用主要体现在后期的RT不稳定性主导阶段。

从图6可知，屈服强度对不同模数的扰动影响不同。这里，以扰动模数为变量分析不稳定性发展的机制，并在给定的屈服强度下，将后期发展最快的扰动模数称为最不稳定模数。图8展示了不同屈服强度下不同模数的扰动随时间增长情况，所考察的材料屈服强度Y0=0～2.0 GPa、扰动模数n=8～40。计算结果显示：纯流体计算(见图8(a))中，最不稳性模数为最大扰动模数；中等程度的屈服强度(见图8(b)～(c))下，随着扰动模数增高，扰动发展速率先加快后减慢；最大的屈服强度(见图8(d))下，最不稳定模数为最小的模数。图9为最不稳定模数随材料屈服强度的变化曲线，随着屈服强度的增大，最不稳定模数单调减小，如图9中小图所示，最不稳定模数与屈服强度的对数近似呈线性关系：nmax≈12.2-17.0 lg(Y0/GPa)。实际上，如图8所示，在扰动发展初期的RM不稳定性发展阶段，受材料强度影响较弱，随扰动模数的增加扰动增长速率均加快。在后期RT不稳定性发展阶段，界面减速是失稳因素，在该因素作用下高模数扰动增长较快；而材料强度是致稳因素，且强度愈大对高模数扰动增长的抑制作用愈强；两个相反因素作用的结果，使得在给定强度下扰动发展存在最不稳定模数。这个最不稳定模数下扰动增长最快，也可以理解为最危险的模数，或者说是该条件下的特征模数。

图8 不同屈服强度下扰动模数与扰动增长的关系Fig.8 Effect of perturbation mode on perturbation growth under different yield strength

图9 最不稳定模数与屈服强度的关系Fig.9 Relationship between yield strength and most instable mode

2.4 金属壳厚度对不稳定性发展的影响

将金属壳厚度4与2.5 mm的扰动发展情况进行了对比，如图10所示。相同扰动模数的薄壳扰动发展均快于厚壳， 在这里，以Y0=0.5 GPa为例，其他强度下扰动发展情况类似。另外，如图11所示，2.5 mm厚的最不稳定模数亦随屈服强度的对数线性减小，nmax≈12.1-11.8 lg(Y0/GPa)，其斜率绝对值比4.0 mm厚的小，即较薄的壳其最不稳定模数随强度的增长变化较缓。

图10 不同金属壳厚度下的扰动增长Fig.10 Effect of shell thickness on perturbation growth

图11 不同金属壳厚度下最不稳定模数与屈服强度的关系Fig.11 Relationship between yield strength and most instable mode with different shell thicknesses

3 结 论

利用自行研制的弹塑性流体力学的欧拉程序，对柱形内爆加载下金属界面不稳定开展了数值研究，结果表明本文的计算结果与文献实验结果吻合较好，进而研究了材料参数对于内爆加载金属界面不稳定发展的影响规律。

对于内爆加载下金属界面不稳定问题，壳体扰动模数、材料强度对界面不稳定性发展均有影响。对于纯流体，扰动模数越高不稳定发展越快；材料强度对界面不稳定性发展有不可忽略的抑制作用，材料屈服强度对高模数不稳定增长的抑制作用较大，而材料剪切模量对不稳定性发展的影响相对较弱；对于相同扰动的金属壳，存在最不稳定模数，且该模数随材料屈服强度的增加而减小，近似与屈服强度的对数呈线性关系；随着壳的厚度减小，扰动增长加快。

感谢秦承森研究员的建设性讨论。

[1] Piriz A R, Lopez Cela J J, Tahir N A, et al. Richtmyer-Meshkov instability in elastic-plastic media[J]. Physical Review E, 2008,78(5):056401.

[2] Piriz A R, Lopez Cela J J, Tahir N A. Linear analysis of incompressible Rayleigh-Taylor instability in solids[J]. Physical Review E, 2009,80(4):046305.

[3] Miles J W. Taylor instability of a flat plate[R]. San Diego: General Dynamics, GAMD-7335, 1966.

[4] Drucker D C. “Taylor instability” of the surface of an elastic-plastic plate[C]∥Nemat-Nasser C. Mechanics Today, Vol. 5. New York: Pergamon Press, 1980:37-47.

[5] Lebedev A I, Nizovtsev P N, Rayevsky V A. Rayleigh-Taylor instability in solids[C]∥4th International Workshop on the Physics of Compressible Turbulent Mixing (IWPCTM). Cambridge, England, 1993.

[6] Barnes J F, Blewett P J, McQueen R G, et al. Taylor instability in solids[J]. Journal of Applied Physics, 1974,45(2):727-732.

[7] Swegle J W, Robinson A C. Acceleration instability in elastic-plastic solids. I. Numerical simulations of plate acceleration[J]. Journal of Applied Physics, 1989,66(7):2838-2858.

[8] Mikaelian K. Effect of viscosity on Rayleigh-Taylor and Richtmyer-Meshkov instabilities[J]. Physical Review E, 1993,47(1):375-383.

[9] Colvin J D, Legrand M, Remington B A, et al. A model for instability growth in accelerated solid metals[J]. Journal of Applied Physics, 2003,93(9):5287-5301.

[10] Atchison W L, Zocher M A, Kaul A M. Studies of material constitutive behavior using perturbation growth in explosive and magnetically driven liner systems[J]. Russian Journal of Physical Chemistry B, 2008,2(3):387-401.

[11] Dimonte G, Terrones G, Cheme F J, et al. Use of the Richtmyer-Meshkov instability to infer yield stress at high-energy density[J]. Physical Review Letters, 2011,107(26):264502.

[12] Frachet V, Geleznikoff F, Guix R, et al. Rayleigh Taylor instability in cylindrical configuration[C]∥2nd International Workshop on the Physics of Compressible Turbulent Mixing (IWPCTM). Pleasanton, CA, 1989.

[13] Park H S, Lorenz K T, Cavallo R M, et al. Viscous Rayleigh-Taylor instability experiments at high pressure and strain rate[J]. Physical Review Letters, 2010,104(13):135504.

[14] 刘军，冯其京，周海兵.柱面内爆驱动金属界面不稳定性的数值模拟研究[J].物理学报,2014,63(15):155201. Liu Jun, Feng Qijing, Zhou Haibing. Simulation study of interface instability in metals driven by cylindrical implosion[J]. Acta Physica Sinica, 2014,63(15):155201.

[15] 冯其京，郝鹏程，杭义洪，等.聚能装药的欧拉数值模拟[J].爆炸与冲击,2008,28(2):138-143. Feng Qijing, Hao Pengcheng, Hang Yihong, et al. Eulerian numerical simulation of a shaped charge[J]. Explosion and Shock Waves, 2008,28(2):138-143.

(责任编辑 丁 峰)

A numerical study of the instability of the metal shell in the implosion

Hao Pengcheng1,2,3, Feng Qijing1, Hu Xiaomian1,2

(1.InstituteofAppliedPhysicsandComputationalMathematics,Beijing100094,China;2.NationalKeyLaboratoryofComputationalPhysics,InstituteofAppliedPhysicsandComputationalMathematics,Beijing100088,China;3.GraduateSchool,ChinaAcademyofEngineeringPhysics,Beijing100088,China)

In this work the instability of the metal shell in the cylindrical implosion was studied numerically using a multi-component elastic-plastic hydrodynamic Eulerian code. Agreeing with those of the experiments, the numerical results show that the material strength restrains the growth of the interfacial perturbation with an effect not to be overlooked. The material yield strength has an obvious restraining effect on the higher mode of perturbation, while that of the shear module is similar but less sensitive. There exists a most instable mode number, which decreases as the yield strength increases and is approximately linear with the logarithm of the yield strength, and the perturbation grows faster as the shell grows thinner.

fluid mechanics; interfacial instability; elastic-plastic simulation; implosion; material parameters

10.11883/1001-1455(2016)06-0739-06

2015-04-30； < class="emphasis_bold">修回日期：2015-08-17

2015-08-17

国家自然科学基金项目(11402029,11372052,11371069,U1430235)； 中国工程物理研究院科学技术发展基金项目(2015B0101021)

郝鹏程(1980— )，男，博士研究生，副研究员,hao_pengcheng@iapcm.ac.cn。

O357.41 <国标学科代码：1302531 class="emphasis_bold"> 国标学科代码：1302531 文献标志码：A国标学科代码：1302531

A