基于LS-DYNA 的液電效應沖擊波數值模擬*

余 慶，張 輝，楊睿智

（中國石油大學（北京）石油工程學院，北京 102249）

自蘇聯Yutkin 把液電效應作為沖擊動力源并廣泛應用于工業加工以來，液電效應研究進入迅速發展階段，學者們開展了大量的研究工作，強大的沖擊波已經被用于體外碎石、液電成形、油田解堵、壓裂等領域。截至目前，液電效應的產生機理仍是人們研究的重點，液相放電擊穿理論主要分為場致電離理論和氣泡理論，其中“氣泡擊穿理論”更為普遍接受。基于該理論，可以將液相放電等離子體的形成過程概括為：電極加熱→氣層（泡）形成→氣層擊穿→形成等離子體通道，推動周圍液體介質形成沖擊波。實際應用中，通常把氣層擊穿之前的過程稱為預擊穿過程，之后的過程稱為主放電過程。由于在液體擊穿過程中伴隨著強烈的熱、電、機械等物理化學變化，因而對于液電效應擊穿機理的研究并不充分，至今仍缺乏一個被研究者普遍接受的擊穿機理解釋理論。而從實驗角度出發研究液相等離子體的特性和機理，通常受到測量儀器的限制，不能準確地獲取等離子體通道內部的特性，并且受環境因素影響嚴重，這使得測量數據的準確性進一步下降。

綜上所述，由于實驗條件的限制以及擊穿理論的不成熟，至今鮮有成熟的商業數值模擬軟件能夠有效模擬液電效應的整個過程，一般模擬方式是通過假設等離子體通道阻抗模型，并結合等離子體通道能量守恒方程、LCR 等效電路方程、沖擊波近似方程進行一維的數值模擬計算。然而，由于通道電阻隨通道內溫度、數密度等因素的變化而變化，準確估計等離子體通道的電阻變化十分困難。工程上通常利用數值模擬軟件對一些具體工況進行模擬，分析液電效應產生的沖擊波對結構物的作用效果，而僅僅是一維的模擬計算數據尚不足以滿足工程需要，因此，開展將液電效應產生的沖擊波運用于已有數值模擬軟件中的研究具有重要意義。

本文中，介紹兩種基于商用顯示動力學軟件LS-DYNA 間接模擬液電效應產生沖擊波的方法：水下爆炸等效方法和理想氣體等效方法，并對計算結果進行比較，給出改進建議，以期為液電效應產生沖擊波的工程應用提供參考依據。

1 液電效應沖擊波特性分析

在介紹液電效應沖擊波等效方法之前，有必要對液電效應沖擊波的形成與能量特性進行分析。圖1為液電效應產生沖擊波的示意圖。可以看出，液電效應產生沖擊波的過程是一個瞬時、連續且復雜的能量轉換過程。首先向電容器充電，當電容充電電壓達到要求后，觸發控制開關，這樣存儲在電容器內的電能迅速注入到位于液電介質中的正負極兩端，位于正負電極之間的水介質逐漸被加熱、氣化，形成等離子體通道并在通道與液體介質界面處向外輻射沖擊波，并且過程中還伴隨著光、熱輻射和熱傳導。在間隙擊穿之后，放電回路可以視為RLC 等效放電回路，圖2 為液相放電等效電路，圖中為儲能電容，和分別為回路電阻和等離子體通道電阻，和分別為回路電感和等離子體通道電感，=+為電路總電感。這里將看作常數，并忽略的影響，基于基爾霍夫定律就可列出微分方程：

圖1 沖擊波形成示意圖Fig. 1 Schematic diagram of shock wave generation

圖2 液相放電等效電路Fig. 2 Equivalent circuit of high-voltage discharge in a liquid

在形成沖擊波的整個過程中主要涉及4 個能量概念：（1）儲存在電容器內的初始能量；（2）液體介質擊穿時儲存在電容器內的能量；（3）注入到等離子體通道內的沉積能量；（4）由沉積能量轉換為等離子體通道內的機械能，即沖擊波能量。這些能量之間的關系14]可以表示為：

式中：η為充電電容儲能轉換為擊穿時電容電能的效率，η為放電能量轉換為等離子體通道沉積能量的效率，η為等離子通道沉積能量轉換為沖擊波能量的效率，η 為總能量轉換效率，為初始電壓，為擊穿電壓，()為等離子體通道電功率，τ 為主放電過程持續時間，()為放電間隙兩端測量的電流，為壓力傳感器到放電通道的水平距離，ρ為液體密度，為液體介質中沖擊波的傳播速度，為測得的沖擊壓力。

為了能夠對上述概念有清楚的認識，這里使用文獻[15-16]中的典型實驗波形，如圖3 所示，圖3(a)描述了高壓探頭距離放電間隙中心17 cm 處等離子體通道的電流、電壓、沖擊波變化，圖3(b)描述了以擊穿時刻作為初始時刻，由式(3) 計算得到的沉積功率與沉積能量變化。通過圖像數據，利用式(1)～(6)分別計算，可得到η=73.38%，η=47.6%，η=28.2%，η=9.85%。由于不同的實驗條件下得到的效率存在明顯差異，因此應根據具體的實驗數據來確定。

圖3 典型的實驗結果Fig. 3 Typical experimental results

而關于液電效應沖擊波峰值壓力的大小一般有近似擬合公式：

式中：為液電效應沖擊波的峰值壓力；主要由傳感器距離決定；α 為常數，根據實際的實驗結果擬合得到。Touya 等根據實驗數據給出了參數與的關系，得到：

式中：的單位為MPa；的單位為mm；的單位為kJ。

2 水下爆炸等效

水下爆炸等效主要是基于能量等效的原則，不同的是一種是爆炸能量等效，一種是沖擊波能量等效。

2.1 爆炸能量等效

爆炸能量等效主要是基于注入到等離子體通道內的能量與炸藥爆炸能量相等的原則來建立相關關系，具體關系式為：

式中：為炸藥爆炸的能量；為TNT 炸藥的質量，kg；為TNT 炸藥的爆炸熱值，一般取值為4.52 MJ/kg。

根據該方法，只要測得放電間隙兩端的電流與電壓值，就可以根據波形得到液電效應注入到間隙內的能量，然后求得炸藥質量。質量給定后，根據等離子體通道的具體形態，將炸藥設置成圓柱形或球形，即可通過LS-DYNA 軟件模擬不同爆心距下水下爆炸沖擊波的峰值壓力。關于水下爆炸沖擊波的研究較多，其中Cole 等、Zamyshlyaev 等關于無限水中爆炸激波的峰值壓力經驗公式應用較廣：

式中：為炸藥中心距離測點的距離，m；為炸藥質量，kg；為沖擊波峰值壓力，MPa；為炸藥包半徑。

假設為0.472 kJ，通過式(9)計算得到TNT 炸藥質量約為1.044 2×10kg。假定炸藥形態為球形，則可以建立如圖3 所示的有限元模型。為了減小計算量，采用1/8 模型，水域的四周設置無反射邊界，水域的對稱面設置法向約束，水域為長方體， 尺寸為12 cm×12 cm×24 cm，單元使用多物質ALE 算法；球形炸藥通過關鍵字*INITIAL-VOLUME-FRACTION-GEOMETRY 添加到水域中；水域和炸藥采用ALE 單元構建。

圖4 水下爆炸數值模型Fig. 4 Numerical models of underwater explosion

TNT 炸藥采用Mat_High_Explosive_Burn 材料模型和JWL 狀態方程，該方程可以描述為：

式中：為爆轟產物的壓力；為單位體積的內能；為比容，即單位體積裝藥產生的爆轟產物的體積；、、、、ω 為JWL 狀態方程參數。具體參數為：密度ρ=1 630 kg/cm，C-J 爆轟壓力=27 GPa，=374 GPa，=7.33 GPa，=4.15，=0.95，ω=0.3，=7 GJ/m。

水采用空物質材料本構模型和Grüneisen 狀態方程，表達式為：

式中：、、為-曲線斜率的系數；為對一階體積的修正；ρ為初始密度；γ為Grüneisen 狀態方程參數；=ρ/ρ–1。具體參數為：ρ=1 020 kg/m，=1 484 m/s，、、、γ、分別為1.979、0、0、0.11、3。

2.2 沖擊波能量等效

不同于爆炸能量等效，沖擊波能量等效主要是基于爆炸產生的沖擊波能量與液電效應產生的沖擊波能量相等來建立相關關系的，具體表達式如下：

式中：為液電效應產生的沖擊波能量，右端項為爆炸產生的沖擊波能量，單位為kJ。

根據該方法，只要利用實驗測得液中沖擊波的壓力時程曲線，或者通過式=η=ηη計算得到沖擊波的能量，即可得到炸藥質量。

同樣假設為0.472 kJ，η為28.2%，通過式(13)計算得到TNT 炸藥質量為3.194 5×10kg。就可以利用有限元模型模擬水下爆炸產生沖擊波的過程。

3 理想氣體等效

不同于將等離子體通道等效于炸藥，等離子體通道還可以看作絕熱膨脹的理想氣體。則等離子體通道內的壓力可以由下面的狀態方程描述：

式中：為等離子體通道內的壓力；ρ為當前等離子體通道內的密度；ρ為初始等離子體通道內的密度；γ 為絕熱指數，這里取1.25。

為了在LS-DYNA 中描述等離子體通道的材料與狀態方程，這里分別使用了關鍵字*MAT_NULL 和關鍵字*EOS_LINEAR_POLYNOMIALWITH_ENERGY_LEAK，這樣就可以表征等離子體通道內的壓力變化，具體形式為：

式中：、、、、均為0；==γ－1，使其形式滿足式(14)；為初始單位體積內能，可以根據沉積功率時程曲線來自定義；=ρ/ρ－1；初始等離子體通道密度ρ=1 000 kg/m。

水的材料與狀態方程的選取以及有限元模型均與水下爆炸等效部分中的模型相同。為了保證計算精度的同時節省計算成本，等離子體通道的半徑選擇為1 mm。同樣假設為0.472 kJ，為了便于輸入到LS-DYNA 中，沉積功率曲線常常簡化為三角波，如圖5 所示，這樣就可以模擬不同爆心距下的沖擊波壓力。

圖5 等離子體通道內的沉積功率Fig. 5 Electrodeposition power-time curves in the plasma part

4 分析與討論

4.1 網格尺寸的選取

一般而言，基于LS-DYNA 的水下爆炸沖擊波峰值壓力的準確模擬與所建立的數值模型密切相關，網格密度、人工黏性以及狀態方程參數等因素均會對峰值壓力產生影響，其中網格尺寸對峰壓的影響尤為顯著，網格太大會影響計算精度，網格太小則會產生不必要的計算成本，因此選擇合適的網格尺寸進行數值模擬至關重要。本文中主要根據水下爆炸經驗公式（式(10)）與數值模擬結果的擬合程度來確定網格尺寸，值得注意的是，這種確定網格尺寸的方式主要基于爆炸等效方法，而未考慮理想氣體等效方法，但這里仍然給出了網格尺寸對基于理想氣體等效方法模擬結果的影響規律，計算結果如圖6 所示。從圖6(a)～(b)中可以看到，隨著網格尺寸減小，不同爆心距下，水下爆炸產生的沖擊波峰壓逐漸增大，與理論公式越來越吻合，說明網格尺寸減小使得計算精度逐漸提高。當網格尺寸從0.10 cm減小到0.08 cm 后，兩種情況下的計算結果差別不大，但在實際的計算過程中，計算網格數目和計算時間卻分別增長了1.95 倍和2.50 倍，因此選擇0.10 cm 的網格尺寸進行數值模擬。不同于水下爆炸等效，基于理想氣體等效方法計算得到的沖擊波峰壓隨著網格尺寸的減小呈先增大后減小的趨勢，峰壓計算結果也比采用水下爆炸等效方法小2 個數量級，該特性多半與沉積能量的輸入方式和傳遞效率有關。

圖6 不同網格大小下峰值壓力的比較Fig. 6 Comparison of the peak pressures under different mesh sizes

4.2 不同等效方法的比較

從圖7 中可以看出，沖擊波以等離子體通道中心為原點，以球面波的形式迅速向外傳播，初始壓力較大（呈現紅色），隨著沖擊波向外擴散，沖擊波壓力逐漸衰減（逐漸變為綠色）。對比圖7(a) 和(c)、圖7(b)和(c)可以看出，采用爆炸等效方法產生的沖擊波壓力顯著高于理想氣體等效方法產生的沖擊波壓力，約高1～2 個量級。通過圖7(a)～(b)可以發現，采用爆炸能量等效方法比采用沖擊波能量等效方法產生的沖擊波壓力值大。而從沖擊波的傳播速度來看，爆炸等效方法產生的沖擊波速度大于理想氣體等效方法，這點可通過=155 μs 時沖擊波的移動距離觀察到。

為了定量分析采用3 種等效方法模擬沖擊波的效果，給出了距放電間隙中心17 cm 處的沖擊波壓力時程曲線和距放電中心不同位置處的沖擊波峰值壓力變化曲線，分別如圖8 和圖9 所示。從圖8 中可以看出，采用爆炸等效方法模擬的沖擊波波速較高，峰值更大，衰減更快，而采用理想氣體等效方法模擬的沖擊波則相反，波速較低，峰值小，衰減慢；采用爆炸能量等效方法比采用沖擊波能量等效方法模擬得到的沖擊波峰值壓力高，接近于其2 倍，這主要是由于前者計算得到的炸藥量更大，根據式(10)可知，峰值壓力的大小與炸藥質量成正相關，炸藥質量越大，峰值壓力更高，這些特點也與圖7 觀察到的現象一致；而且可以發現，3 種等效方法模擬得到的波形在衰減階段都存在明顯的振蕩，說明數值模擬夸大了實際壓力波動，這可能是由于數值模擬中對水的黏滯性和導熱性造成的能量耗散考慮不足導致的，可引入人工黏性來改進，但是引入人工黏性又會引起峰壓下降，本文中主要以沖擊波峰壓作為主要分析因素，因此不再調整人工黏性系數。從圖9 中可以看出，隨著與放電中心距離的增大，不管采用哪種方法模擬得到的沖擊波峰值壓力都是逐漸減小的。同一位置處，爆炸能量等效的峰值壓力最高，沖擊波能量等效的峰值壓力次之，理想氣體等效的峰值壓力最低，而且對比Touya 經驗公式（式(8)），發現采用等效爆炸方法模擬出的沖擊波峰值壓力普遍要比Touya 經驗公式要高，這與文獻[20]中的模擬結果是一致的。只有當爆心距較大時，采用沖擊波能量等效方法模擬得到的峰值壓力才越來越趨近于Touya 經驗公式計算得到的峰值壓力。

圖7 不同等效方法下的沖擊波傳播過程Fig. 7 Propagation process of the shock wave using different equivalent methods

圖8 不同等效方法下沖擊波壓力時程曲線Fig. 8 Shock wave pressure-time history curves using different equivalent methods

圖9 距離放電中心不同位置處的沖擊波峰值壓力Fig. 9 Shock wave peaks at different points from the discharge center

4.3 方法改進與修正

通過對比，發現上述兩類等效方法各有優缺點，采用爆炸等效，關鍵在于準確估計出炸藥質量，而炸藥質量的準確估計取決于能量轉換效率以及等離子體通道內沉積能量的準確計算。在實際工程作業中，由于工況的復雜性、測量手段的局限性以及放電擊穿的隨機性，會導致沉積能量和效率的計算誤差較大；采用理想氣體等效方法模擬雖然能夠詳細地描述等離子體通道內的注入能量，在液電成形方向模擬效果良好，但是模擬得到的峰值壓力卻不理想。為此，以在相同爆心距下的峰值壓力相等為基準，對上述等效方法進行改進。

對于爆炸等效方法，可按照以下步驟進行改進。

（1）通過調整數值模型的網格大小、人工黏性系數等參數，使得數值模擬得到的水下爆炸沖擊波峰值壓力與式(10)較吻合；

（2）聯立Touya 經驗公式（式(8)）（如果條件允許也可通過實驗得到水中不同位置處的峰壓與放電能量之間的關系式來代替）和爆炸經驗公式（式(10)），得到如下關系式：

（3）根據沖擊波與結構物的作用區域來確定式(16)中的或。例如：當爆炸激波與結構物的作用區域在5～15 cm 之間時，和的取值范圍分別為0.37～0.57 和0.70～0.82，可分別取0.47 和0.76。

（4）計算某一擊穿能量下的炸藥質量。求得和，給定擊穿能量即可根據式(16)計算得到炸藥質量和球形炸藥半徑，例如：=472 J，則=/η≈992 J，計算得到炸藥質量為1.565 8×10kg，炸藥半徑為2.483 2 mm。

（5）將炸藥半徑輸入到LS-DYNA 中，即可計算不同爆心距下的沖擊波峰壓。如圖9 所示，可以看到，修改后的等效爆炸模型在10～15 cm 處的峰值壓力略底于Touya 經驗公式，但是整體十分吻合。

而對于理想氣體等效方法，可采取增加輸入沉積功率的方式來進行改進，如圖9 所示。通過大量的試算，發現當輸入到LS-DYNA 中的沉積功率為測量計算得到的沉積功率的30 倍，甚至100 倍時，計算得到的沖擊波峰值壓力有明顯提高，但也依然小于Touya 經驗公式計算得到的沖擊波峰壓，而此時輸入的峰值沉積功率達到4 056 MW，與實際差距較大，說明采用理想氣體等效方法模擬沖擊波峰值壓力效果不太理想。

4.4 沉積能量對峰壓的影響

上述模型的建立過程中，僅僅給出了等離子體通道沉積能量為471 J 時采用不同等效方法模擬沖擊波峰壓的結果對比，而沉積能量是決定沖擊波峰壓大小的關鍵因素之一，因此有必要探究不同沉積能量下各等效方法的結果差異。依然采用文獻[15-16]中的實驗數據，利用LS-DYNA 計算了沉積能量分別為176 和301 J 時的沖擊波峰壓，計算結果如圖10 所示。結合圖9 和圖10 可以看到，不同沉積能量下，采用爆炸能量等效方法計算得到的峰值壓力相較于其他等效方法和Touya 經驗公式總是最高的；隨著沉積能量減小，采用沖擊波能量等效方法計算得到的峰值壓力經歷了先高于、然后趨近、再低于Touya 經驗公式計算峰壓的變化過程，這主要是由液電效應放電特性決定的，電極間距離減小導致的沉積能量的減小會導致η減小，但η基本保持不變，而Touya 經驗公式主要是由放電能量（和η）確定的，沖擊波能量等效方法計算得到的炸藥質量主要由和η確定，在放電能量變化不大時，Touya 經驗公式計算得到的峰壓變化也不大，但是由于電極間距離減小導致的沉積能量的減小會使得沖擊波能量減小，進而計算得到的炸藥質量也會減小，因而會出現沖擊波能量等效方法計算得到的峰壓逐漸小于Touya 經驗公式的現象；采用理想氣體等效方法相比其他等效方法模擬計算得到的峰壓總是最低的，小1～2 個數量級；采用改進的等效爆炸方法在不同沉積能量下均能與經驗公式擬合較好，說明改進方法是可行的。

圖10 不同沉積能量下沖擊波峰值壓力曲線Fig. 10 Shock wave peaks under different deposited energies

5 結 論

為了有效模擬液電效應產生的沖擊波以滿足工程需要，介紹了兩種基于LS-DYNA 的數值模擬方法，水下爆炸等效方法（分為爆炸能量等效與沖擊波能量等效）和理想氣體等效方法，并在等離子體通道沉積能量為472 J 的情況下對模擬方法進行比較與改進，分析了不同沉積能量下各等效方法峰壓計算結果的差異，得到以下主要結論。

（1）模擬得到沖擊波峰值壓力排序為，爆炸能量等效的沖擊波峰值壓力最高，沖擊波能量等效的沖擊波峰值壓力次之，理想氣體等效沖擊波峰值壓力最低，理想氣體等效模擬的峰壓相較于前2 種方法小1～2 個數量級；沖擊波波速大小排序為，爆炸能量等效與沖擊波能量等效的沖擊波波速相等，且高于理想氣體等效的沖擊波波速。

（2）放電能量變化不大時，電極間距離減小引起的沉積能量減小會使得3 種等效方法（爆炸能量等效、沖擊波能量等效、理想氣體等效）模擬得到的峰壓均有不同程度的減小；3 種等效方法得到的峰壓的大小順序不發生變化；改進后的等效爆炸方法能夠適應沉積能量的變化，與Touya 經驗公式擬合較好。

（3）基于LS-DYNA 對液電效應沖擊波峰值壓力進行準確模擬，除了選取合適的等效方法，還應結合具體的放電條件，建立合適的數值模型，在滿足計算要求的條件下實現沖擊波峰值壓力的快速計算。