基于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 对液电效应冲击波峰值压力进行准确模拟，除了选取合适的等效方法，还应结合具体的放电条件，建立合适的数值模型，在满足计算要求的条件下实现冲击波峰值压力的快速计算。