王 濤，汪 兵，林健宇，鐘 敏，柏勁松，李 平，陶 鋼

（1. 中國工程物理研究院流體物理研究所，四川 綿陽 621999；2. 南京理工大學能源與動力工程學院，江蘇 南京 210094）

當沖擊波或慣性力加載不同物質(zhì)間的擾動界面時，界面會失穩(wěn)并發(fā)生擾動增長的現(xiàn)象，前者稱為Richtmyer-Meshkov（RM）不穩(wěn)定性[1-2]，后者稱為Rayleigh-Taylor（RT）不穩(wěn)定性[3-4]。其中RT 不穩(wěn)定性與加載方向有關，只有當輕介質(zhì)加速重介質(zhì)時才會發(fā)生。界面不穩(wěn)定性在諸多領域有重要的應用背景，如慣性約束聚變[5-6]、超新星爆炸[7]、小行星撞擊[8-9]、地球內(nèi)核運動及板塊構(gòu)造[10-11]等，因此具有重要的研究意義。

Wang 等[12-18]、Bai 等[19-21]、Xiao 等[22]長期研究流體介質(zhì)界面不穩(wěn)定性及湍流混合演化規(guī)律和統(tǒng)計特性。Liu 等[23]通過直接數(shù)值模擬，研究了RM 不穩(wěn)定性和湍流混合中的能量傳遞過程和機制。李俊濤等[24-25]也通過數(shù)值模擬，研究了弱沖擊波加載V 形界面導致的RM 不穩(wěn)定性中的渦動力學行為特性及作用機制。Luo 等[26-27]、Si 等[28]、Ding 等[29]、Lei 等[30]在匯聚沖擊波誘導下的氣體界面不穩(wěn)定性實驗研究方面成績卓越，特別是首次測得了柱形單模態(tài)air/SF6界面的擾動振幅。在我國，有關界面不穩(wěn)定性的研究主要集中在流體介質(zhì)方面，很少涉及金屬材料界面不穩(wěn)定性，而界面不穩(wěn)定性也會在金屬材料中發(fā)生。而且，金屬材料的復雜物性，如強度、相變、損傷、微結(jié)構(gòu)等，使金屬材料界面不穩(wěn)定性比流體界面不穩(wěn)定性更復雜，因此，金屬材料界面不穩(wěn)定性是一種跨尺度多物理耦合的復雜流動問題。

金屬材料界面不穩(wěn)定性的理論研究，主要基于能量平衡[31-33]或力平衡[34-36]，推導擾動增長所滿足的色散關系，但是，這些線性理論分析均采用了理想塑性本構(gòu)關系，預測擾動增長受到較大限制。我們也基于能量平衡，但采用了Steinberg-Guinan（SG）和Johnson-Cook（JC）本構(gòu)模型及變壓力加載歷程，針對有限厚度金屬平板問題，推導了相應的擾動增長方程，將線性分析方法進行了拓展，可用于平面幾何金屬不穩(wěn)定性線性和非線性段增長的預測[37]。并且，對爆轟驅(qū)動和激光等離子體驅(qū)動的金屬界面不穩(wěn)定性實驗進行了分析，結(jié)果和實驗吻合較好，但是，這種線性分析方法依然有較大的局限性，如它還無法考慮加載過程等。

金屬材料界面不穩(wěn)定性實驗研究始于20 世紀70 年代。首先，由Barnes 等[38-39]通過爆轟驅(qū)動鋁平板實現(xiàn)，并采用高能X 射線裝置觀察界面擾動增長。后來，平面爆轟驅(qū)動金屬材料界面不穩(wěn)定性實驗[40-42]均采用相似的裝置。在美國、俄羅斯，開展了大量有關金屬材料界面不穩(wěn)定性的研究工作，包括界面擾動增長規(guī)律及影響因素的研究等。何長江等[43]通過數(shù)值模擬，研究了初始擾動波長對平面爆轟驅(qū)動金屬鋁RT 不穩(wěn)定性發(fā)展的影響。郝鵬程等[44]、劉軍等[45]通過對內(nèi)爆加載金屬鋼殼的數(shù)值模擬，研究了材料強度和初始波長對界面擾動增長的影響。Olson 等[46]在爆轟實驗中，研究了銅材料的微結(jié)構(gòu)和加工對RT 擾動增長的影響，發(fā)現(xiàn)單晶晶向和材料加工導致的應變硬化會影響擾動增長，而多晶材料的晶粒尺寸和晶界強化在所研究的加載條件下對擾動增長無影響。Jensen 等[47]采用高分辨率的原位診斷技術(shù)，研究了沖擊條件下鈰的RM 不穩(wěn)定性問題。Cherne 等[48]采用分子動力學模擬方法，研究了銅樣品自由面上不同形狀的二維溝槽擾動對銅/真空界面RM 不穩(wěn)定性發(fā)展的影響。

金屬材料界面不穩(wěn)定性增長是一種高壓高應變率條件下的大變形行為，與材料強度密切相關，宏觀表現(xiàn)就是材料強度可以抑制擾動增長。近幾年，基于這個思想，利用金屬材料界面不穩(wěn)定性增長，研究高壓高應變率條件下的材料強度[42,49-53]，進而發(fā)展了適用性更強的材料強度模型[54]。

由以上分析，目前對于金屬材料界面不穩(wěn)定性問題，雖然已經(jīng)開展了較多的研究，但是還不充分不深入，仍然有很多問題有待解決。比如，理論研究很不完備，實驗研究缺乏對物理細節(jié)的直接認識，而且高分辨率診斷設備和極端加載裝置限制了實驗工作的全面開展。因此，數(shù)值模擬成為一種有效便捷的研究手段。我們長期致力于流體界面不穩(wěn)定性及湍流混合問題的數(shù)值模擬研究[12-22]。在此基礎上，我們發(fā)展了適用于金屬材料大變形行為數(shù)值模擬的高保真度爆轟與沖擊動力學歐拉程序。本文中，利用自研的數(shù)值模擬歐拉程序，研究柱形匯聚幾何中內(nèi)爆驅(qū)動下金屬材料界面不穩(wěn)定性的動力學行為。

1 計算方法

針對考慮爆轟和材料彈塑性行為的多物質(zhì)、大變形及強沖擊波物理問題，發(fā)展了高保真度的爆轟與沖擊動力學歐拉有限體積計算程序。其守恒型控制方程組為：

式中：下標i、j 分別代表x、y、z 三個方向，遵循張量運算法則，V 為控制體體積，S 為控制體表面積，ni為其外法向單位矢量分量，ρ、uk(k=i, j)、p、E 分別為密度、速度、壓強和質(zhì)量總能量，sij為偏應力張量分量。

采用維數(shù)分裂技術(shù)，將方程組（1）描述的物理問題分解為多個一維問題，對于每個一維問題，采用PPM 方法，對單元內(nèi)物理量的分布進行插值和重構(gòu)；然后，通過歐拉型兩步算法進行計算，即物理量的Lagrange 推進求解，再將拉氏網(wǎng)格上的物理量映射到靜止的歐拉網(wǎng)格上。材料強度效應、炸藥爆轟過程和人工黏性在Lagrange 步中實現(xiàn)。多物質(zhì)界面采用體積分數(shù)方法進行捕捉。

1.1 狀態(tài)方程

數(shù)值模擬中，炸藥用JWL 狀態(tài)方程描述：

E=ρ0E′E′

式中：相對比容V=ρ0/ρ，A、B、R1、R2、ω 為常數(shù)， 為體積內(nèi)能， 為比內(nèi)能，cV為比熱。

金屬材料用Mie-Grüneisen 狀態(tài)方程描述：

式中：相對壓縮度 μ =ρ/ρ0?1 ，ρ0為初始密度，c 為零壓聲速，γ0為Grüneisen 系數(shù)，a、S1、S2、S3為常數(shù)。

硅橡膠采用凝聚介質(zhì)實用狀態(tài)方程：

式中：c0為正常態(tài)的聲速，γ 為材料參數(shù)。

1.2 本構(gòu)方程

金屬材料采用SG 本構(gòu)模型，它可以很好地描述高壓高應變率下的材料強度特性。SG 本構(gòu)模型在彈塑性本構(gòu)方程中引入壓力、溫度和應變率項，同時壓力與應變率對流動應力的耦合效應具有可分離變量特性。因為SG 本構(gòu)模型中流動應力依賴于壓力，所以材料的本構(gòu)方程與狀態(tài)方程之間存在某種耦合關系。這個耦合關系反映了高壓下金屬材料的壓力硬化特性。SG 本構(gòu)模型的動態(tài)屈服強度和剪切模量分別為：

式中：Y0和G0分別為初始屈服強度和剪切模量，β 和n 分別為材料應變硬化系數(shù)和硬化指數(shù)，A 為壓力硬化系數(shù)，η=ρ/ρ0為材料壓縮比，B 為溫度軟化系數(shù)，為零溫下的比內(nèi)能，C 為材料常數(shù)。

2 結(jié)果與分析

在進行柱面內(nèi)爆驅(qū)動不銹鋼金屬材料界面不穩(wěn)定性計算和分析前，先通過對文獻[40]中的爆轟驅(qū)動T-6061 鋁的RT 不穩(wěn)定性實驗的模擬，確認本計算程序的可靠性。該實驗采用HMX 炸藥爆轟加載1.5 mm 厚的鋁樣品，加載壓力峰值為30 GPa。初始擾動波長λ0=2 mm，初始擾動振幅a0=0.15，0.11 mm。圖1 為用自研的歐拉程序計算得到的擾動振幅和文獻[40]中采用ARES 程序計算的結(jié)果，及實驗結(jié)果（橫軸是自由面運動位移）?？梢钥闯?，數(shù)據(jù)吻合較好，說明自研的歐拉程序模擬這類爆轟驅(qū)動金屬材料界面不穩(wěn)定性是準確的。

圖 1 爆轟驅(qū)動鋁實驗的擾動振幅Fig. 1 Perturbation amplitudes of experiments driven by explosion

2.1 柱面內(nèi)爆計算模型和參數(shù)

以Frachet 等[55]的柱面內(nèi)爆驅(qū)動不銹鋼金屬材料的雙模態(tài)擾動界面不穩(wěn)定性實驗為基礎，開展柱面內(nèi)爆驅(qū)動下金屬材料界面不穩(wěn)定性發(fā)展演化的動力學行為特性數(shù)值模擬。

計算模型如圖2 所示：外層為炸藥，其外直徑和厚度分別為200 和50 mm，由于文獻[55]中沒有說明何種炸藥，本文中采用TNT；中心物質(zhì)為直徑92 mm 的硅橡膠；中間不銹鋼殼外直徑為100 mm，厚度為4 mm，且內(nèi)界面上預置有振幅0.5 mm、模數(shù)為13 和29 mm?1的雙模態(tài)初始擾動。

T N T 炸藥J W L 狀態(tài)方程參數(shù)分別為：ρ0=1.63 g/cm3，pcj=21.0 GPa，Dcj=6.93 km/s，A=371.2 GPa，B=3.231 GPa，R1=4.15，R2=0.95，ω=0.3。不銹鋼Mie-Grüneisen 狀態(tài)方程參數(shù)分別為：ρ0=7.896 g/cm3，c=4.569 km/s，γ0=2.17，a=0.47，S1=1.49，S2=S3=0；不銹鋼SG 本構(gòu)模型參數(shù)分別為：Y0=0.34 GPa，Ymax=2.5 GPa，G0=77.0 GPa，β=43.0，n=0.35，A=0.022 6 GPa?1，B=0.455 kK?1。硅橡膠凝聚介質(zhì)實用狀態(tài)方程參數(shù)分別為：ρ0=1.15 g/cm3，c0=1.343 km/s，γ=1.6。

圖 2 柱面內(nèi)爆驅(qū)動金屬材料界面不穩(wěn)定性計算模型Fig. 2 Computational model of metal interface instability driven by cylindrical implosion

2.2 加載和界面運動特性

沖擊波(SW)由低阻抗介質(zhì)向高阻抗介質(zhì)方向加載物質(zhì)界面時，會產(chǎn)生一個透射沖擊波和一個反射沖擊波；沖擊波由高阻抗介質(zhì)向低阻抗介質(zhì)方向加載物質(zhì)界面時，會產(chǎn)生一個透射沖擊波和一個反射稀疏波(RW)。炸藥起爆后，產(chǎn)生的柱形爆轟波向內(nèi)聚心運動，到達鋼殼外界面并作用后，向內(nèi)產(chǎn)生透射沖擊波，鋼殼也開始聚心運動，由于質(zhì)量守恒，鋼殼徑向厚度逐漸增大（見圖3）；該透射沖擊波約在7.8 μs時到達鋼殼內(nèi)界面，形成由高阻抗介質(zhì)向低阻抗介質(zhì)方向的沖擊波加載，內(nèi)界面加載壓力瞬間達到10 GPa（見圖4），同時向硅橡膠內(nèi)產(chǎn)生透射沖擊波和鋼殼內(nèi)的反射稀疏波；鋼殼內(nèi)的反射稀疏波向外運動并加載外界面后又向內(nèi)反射壓縮波(CW)，壓縮波和稀疏波會在鋼殼內(nèi)多次反射，并在內(nèi)外界面上分解而交替產(chǎn)生稀疏波和壓縮波，進而交替加載鋼殼界面，可見于鋼殼內(nèi)三角形的胞格結(jié)構(gòu)（見圖3）、加載壓力和內(nèi)界面速度振蕩區(qū)域（沖擊波加載和稀疏波卸載效應）（見圖4、5，界面運動速度以徑向向外為正）；硅橡膠內(nèi)的透射沖擊波聚心反彈后，在約20.5 μs 時到達并沖擊(RS)鋼殼內(nèi)界面，內(nèi)界面加載壓力瞬間達到約50 GPa，鋼殼聚心運動停止，開始向外反彈運動，其徑向厚度之后又逐漸減小，并向硅橡膠內(nèi)反射沖擊波和向鋼殼內(nèi)透射沖擊波，該反射沖擊波會再次聚心反彈加載，但已不在本算例時間內(nèi)，而鋼殼內(nèi)的透射沖擊波也會在內(nèi)外界面上分解而交替產(chǎn)生壓縮波和稀疏波。

由從鋼殼內(nèi)界面運動速度（見圖5）和加速度剖面（見圖6）看出，內(nèi)界面首次受到?jīng)_擊加載后擾動開始反相（峰谷轉(zhuǎn)換），發(fā)生RM 不穩(wěn)定性，接著又受到三次有效的向心反射壓縮波加載，在約12 μs 前，雖然有輕介質(zhì)加速度重介質(zhì)的RT 不穩(wěn)定性，但以RM 不穩(wěn)定性為主導。在12 μs 到?jīng)_擊波聚心反彈加載前，由于中心區(qū)域壓力的升高，界面聚心運動加速減速，擾動發(fā)展由輕介質(zhì)加速重介質(zhì)的RT 不穩(wěn)定性主導。沖擊波聚心反彈加載后及之后內(nèi)界面受到四次有效的向心反射稀疏波加載，擾動發(fā)展由RM 不穩(wěn)定性主導。

圖 3 密度顯示的一維波譜圖Fig. 3 One dimensional wave diagram displayed by density

圖 4 不銹鋼殼內(nèi)外界面加載壓力Fig. 4 Loading pressures on inner and outer interface

圖 5 不銹鋼殼內(nèi)界面運動速度Fig. 5 Velocity of inner interface

圖 6 不銹鋼殼內(nèi)界面加速度Fig. 6 Acceleration of inner interface

2.3 擾動界面演化

圖7 為不同網(wǎng)格分辨率（146、73、36 μm）下鋼殼內(nèi)界面雙模態(tài)擾動各模態(tài)（模數(shù)13 和29 mm?1）振幅增長過程，粗網(wǎng)格與中等網(wǎng)格、中等網(wǎng)格與細網(wǎng)格結(jié)果的最大相對誤差分別為3%、1.5%，滿足收斂要求，所以計算中采用36 μm 的網(wǎng)格。
圖8 為5、10、15、20、25 μs 時密度場圖像，顯示了鋼殼在內(nèi)爆沖擊波和聚心反彈等加載下的運動過程和內(nèi)界面擾動演化過程。圖9 為17.8、22.8 和25.3 μs 時鋼殼內(nèi)界面擾動演化的實驗圖像和數(shù)值模擬圖像。圖10 為對雙模態(tài)擾動演化做譜分析得到的0、10、15、20、25 和30 μs時各模態(tài)的擾動振幅，在多次沖擊加載下，所有模態(tài)都被激發(fā)，但主導界面演化的還是模數(shù)13 和29 mm?1兩個模態(tài)。圖11 為雙模態(tài)擾動振幅增長歷史，數(shù)值模擬結(jié)果和實驗差別較大，原因是文獻[55]中并沒有說明內(nèi)爆采用的是何種炸藥，而本模擬中采用TNT，TNT 爆炸后加載在鋼殼上的壓力和文獻[55]中的壓力有較大差別。由圖11 還可看出：低模數(shù)擾動增長比高模數(shù)擾動增長速度快，這與純流體情況剛好相反，是由金屬材料彈塑性機制導致；而且，對于低模數(shù)擾動，在擾動反相后至沖擊波聚心反彈加載前，擾動振幅以近似線性規(guī)律增長，而高模數(shù)擾動由于較強的非線性作用，非線性增長規(guī)律明顯。

圖 7 網(wǎng)格收斂性分析Fig. 7 Grid convergence

圖 8 密度場Fig. 8 Images of density fields

圖 9 鋼殼內(nèi)界面的擾動Fig. 9 Inner perturbed interface of steel shell

圖 10 雙模態(tài)擾動演化的譜分析Fig. 10 Spectral analysis of dual mode perturbation evolution

圖 11 雙模態(tài)擾動的振幅增長曲線Fig. 11 Amplitude growth curves of dual mode perturbation

2.4 初始擾動振幅、波長、鋼殼厚度和幾何特性的影響

本節(jié)主要研究初始擾動振幅、波長、鋼殼厚度和幾何特性對柱形內(nèi)爆驅(qū)動下鋼殼內(nèi)界面擾動增長的影響。

首先，研究初始擾動振幅的影響。保持初始擾動模數(shù)為13 mm?1，初始振幅分別為0.05、0.10、0.25、0.50、0.75 和1.00 mm。圖12 為不同初始擾動振幅時的振幅增長，較大的初始振幅意味著較大的振幅增長，而且擾動反相的時間與初始振幅無關。

其次，研究初始擾動波長的影響。保持初始振幅為0.5 mm，初始擾動的模數(shù)(波長的倒數(shù))分別為13、16、20、24、29、35、40、45、50、56、60 mm?1。圖13 為不同初始擾動模數(shù)時的振幅增長，高模數(shù)（小波長）擾動反相時間明顯比低模數(shù)（大波長）擾動反相時間短，即擾動反相時間與擾動模數(shù)相關。反相后高模數(shù)擾動振幅增長速度比低模數(shù)擾動振幅增長速度快，但是低模數(shù)擾動振幅會以該速度近似線性增長，直至聚心反彈沖擊波加載時截止；而高模數(shù)擾動振幅增長速度會逐漸減小，其振幅及增長速度會小于低模數(shù)擾動；而且，隨著初始擾動模數(shù)的增大，振幅增長在后期會停止，說明柱形內(nèi)爆驅(qū)動金屬材料界面不穩(wěn)定性增長存在截止波長。

圖 12 不同初始擾動振幅時的振幅增長曲線Fig. 12 Amplitude growth curves for different initial perturbation amplitude

圖 13 不同初始擾動模數(shù)時的振幅增長曲線Fig. 13 Amplitude growth curves for different initial perturbation mode number

接著，研究鋼殼初始厚度的影響。初始擾動模數(shù)分13 和60 mm?1兩種，初始振幅均為0.5 mm，初始厚度分別為4、7 和10 mm。圖14 為不同鋼殼初始厚度時振幅增長情況：隨著鋼殼初始厚度的增大，擾動振幅增長速度逐漸減??；而且，對于模數(shù)為60 mm?1的初始擾動，不同厚度時的擾動振幅幾乎相同，后期增長也趨于停止，即鋼殼厚度會抑制擾動增長，而且存在一個截止厚度。

最后，研究幾何特性對沖擊驅(qū)動鋼殼界面不穩(wěn)定性增長的影響。保持初始時刻的炸藥厚度、鋼殼厚度、硅橡膠厚度、擾動振幅和模數(shù)均相同的情況下，分為柱形匯聚幾何和平面幾何構(gòu)型。初始振幅為0.5 mm，模數(shù)為13 mm?1。圖15 為不同幾何構(gòu)型下的振幅增長情況，匯聚幾何構(gòu)型下的擾動振幅比平面情況下大，這可以由兩種情況下界面運動速度(見圖16)和界面加速度(見圖17、6)解釋。首先，幾何匯聚帶來幾何收縮效應（BP 效應），由于能量匯聚使加載壓力較大，界面加速運動較快，因此擾動增長速度也較快；其次，當界面運動減速時，幾何收縮會造成界面加速減速，而平面情況下界面是減速減速，因此減速運動階段的RT 不穩(wěn)定性發(fā)展在匯聚情況下更劇烈；再者，匯聚情況下沖擊波聚心反彈加載時間早于平面情況，而且反彈加載強度也較大，因此之后的RM 不穩(wěn)定性發(fā)展在匯聚情況下也更劇烈。所以，整體上匯聚幾何中界面不穩(wěn)定性發(fā)展速度快。

圖 14 不同鋼殼初始厚度時的振幅增長曲線Fig. 14 Amplitude growth curves for different initial thickness of steel shell

圖 15 不同幾何構(gòu)型下的振幅增長曲線Fig. 15 Amplitude growth curves for different geometrical configuration

圖 16 不同幾何構(gòu)型下的界面運動速度Fig. 16 Interface velocities for different geometrical configuration

圖 17 平面幾何中的界面加速度Fig. 17 Accelerations of interface in planar geometry

3 結(jié) 論

基于歐拉有限體積法，發(fā)展了適用于多物質(zhì)、大變形及強沖擊波條件下流體彈塑性問題數(shù)值模擬的高保真度爆轟與沖擊動力學計算程序。利用該數(shù)值模擬程序，研究了柱形匯聚幾何中內(nèi)爆驅(qū)動下金屬材料界面不穩(wěn)定性發(fā)展的動力學行為特征，得出以下結(jié)論。

（1）首次沖擊后，鋼殼開始聚心運動，在約12 μs 前，還受到了三次有效的向心反射壓縮波加載，在此階段，界面發(fā)展以RM 不穩(wěn)定性為主；在12 μs 后到?jīng)_擊波聚心反彈加載之前，界面聚心運動處于一種加速減速狀態(tài)，界面發(fā)展由RT 不穩(wěn)定性主導；在沖擊波聚心反彈加載后，界面又受到四次有效的向心反射稀疏波加載，此時界面發(fā)展又由RM 不穩(wěn)定性主導。

（2）初始擾動振幅較大，振幅增長較大；初始擾動模數(shù)較大，振幅增長較小，且存在一個截止模數(shù)（波長）；鋼殼厚度會抑制擾動增長，也存在一個截止厚度；在相同初始條件下，相較于平面幾何構(gòu)型，匯聚幾何結(jié)構(gòu)中由于幾何收縮效應和界面減速階段較強的RT 不穩(wěn)定性，使擾動增長速度更快。

金屬材料和流體最大的區(qū)別在于金屬材料的強度效應，而強度對擾動發(fā)展具有抑制作用。另外，由此帶來的材料斷裂破壞對擾動增長甚至混合都會產(chǎn)生影響，而且材料的斷裂破壞機制也復雜多變，其他如金屬材料的相變、微觀結(jié)構(gòu)特征等，都會影響金屬材料界面不穩(wěn)定性的發(fā)展，這都使金屬材料界面不穩(wěn)定性問題更加復雜。下一步，我們擬研究材料強度對金屬材料界面不穩(wěn)定性發(fā)展演化的影響，包括硬化、軟化等機制。