Experiment and Simulation of Carbon Dioxide Hydrate Formation Mechanism under High Pressure

光滑粒子流体动力学(Smoothed Particle Hydrodynamics, SPH)方法中的搜索算法较耗时，即每个时间步都要对领域粒子进行搜索，粒子越多，耗时情况越突出，与有限元法相比，SPH方法的计算效率要低得多。为了解决SPH方法计算效率低的问题，Johnson等^[1-2]和Attaway等^[3]将有限元与SPH方法相结合，提出了SPH-FEM耦合的算法，即:在小变形区域使用有限元法，大变形区域仍使用SPH方法。该方法不仅提高了计算效率，而且适应性较强。

目前，采用SPH方法对爆炸焊接进行数值模拟的相关报道较少，而且多采用二维SPH方法。Tanaka^[4]采用SPH方法对爆炸焊接的斜碰撞过程进行了数值模拟，成功地模拟出射流、波形和涡旋，波长的模拟结果相对实验结果偏大。李晓杰等^[5]采用SPH方法及热塑性流体力学模型对爆炸复合板的斜碰撞过程中出现的界面波进行了数值模拟，模拟结果与张登霞等^[6-7]实验结果的一致性较好。刘江等^[8]利用AUTODYN软件中的SPH方法模拟了爆炸复合的斜碰撞，结合模拟中有效塑性变形、温度及剪切应力呈现的变化规律发现，爆炸复合的结合机理集塑性变形、熔化和扩散为一体。本研究将采用三维SPH方法对双面爆炸焊接过程进行模拟，将其结果与实验及理论结果进行对比，分析SPH-FEM耦合方法对爆炸焊接模拟的有效性。

1.   计算模型及参数选取

1.1   计算模型

以前期45钢/Q235钢双面爆炸焊接实验^[9]为基础，考虑到计算效率，利用LS-DYNA建立如图 1及图 2所示的两组双面爆炸焊接SPH-FEM耦合的三维真实计算模型，选用的炸药为乳化炸药(玻璃微球的质量分数为5%)，计算模型中基板和复板的材料、尺寸、间隙(δ)及药厚如表 1所示。起爆方式为点起爆。

图  1  计算模型Ⅰ(10 mm药厚)

Figure  1.  Calculation model Ⅰ with explosive thickness of 10 mm

下载: 全尺寸图片 幻灯片

图  2  计算模型Ⅱ(5 mm药厚)

Figure  2.  Calculation model Ⅱ with explosive thickness of 5 mm

下载: 全尺寸图片 幻灯片

表  1  计算模型中材料的相关参数

Table  1.  Related parameters of materials in calculation models

Calculationmodel	Flyer plate			Base plate		Gap δ/mm	Size of explosive/(mm×mm×mm)
	Material	Size/(mm×mm×mm)		Material	Size/(mm×mm×mm)
Ⅰ	45 steel	300×150×2		Q235	300×150×16	6	300×150×10
Ⅱ	45 steel	300×150×2		Q235	300×150×16	6	300×150×5

下载: 导出CSV 

| 显示表格

基、复板采用3D Solid 164实体单元，单元边长为0.1 cm；炸药划分为光滑粒子，粒子的大小Δr取为0.1 cm。考虑到模型的对称性，为了提高计算效率，采用1/2模型进行计算。单位制为cm-g-μs。

1.2   材料模型及参数设定

数值计算中乳化炸药采用高能燃烧模型^[10-11]及JWL状态方程^[12]。JWL状态方程表达式为

式中：A_JWL、B_JWL、R₁、R₂和ω为材料参数；p为爆轰产物压力，GPa；E₀为初始比内能，kJ/cm³；v为爆轰气体产物的相对比容，为无量纲量。炸药的相关参数见表 2，其中：ρ为密度，D为炸药爆速。

表  2  乳化炸药的JWL状态参数^[13]

Table  2.  JWL equation-of-state parameters of emulsion explosives^[13]

ρ/(g·cm-3)	D/(m·s-1)	AJWL/GPa	BJWL/GPa	R1	R2	ω	E0/(kJ·cm-3)
1.12	4 510	326.42	5.808 9	5.80	1.56	0.57	3.323

下载: 导出CSV 

| 显示表格

数值计算中，基、复板均采用Mie-Grüneisen状态方程^[14]和Johnson-Cook材料模型^[15]。Johnson-Cook材料模型的形式如下

式中：ε_p为有效塑性应变；$\dot{\varepsilon }_{\text{p}}^{*}={{\dot{\varepsilon }}_{\text{p}}}/\dot{\varepsilon }_{\text{p}}^{0}$为有效塑性应变率，其中$\dot{\varepsilon }_{\text{p}}^{0}$为参考应变率；A、B、C、m及n为与材料相关的常数；无量纲温度T^*表示为T^*=(T-T_r)/(T_m-T_r)，其中T_r为室温, T_m为熔点。45钢选用与Q235钢相同的Johnson-Cook材料模型参数，具体参数如表 3所示。

表  3  Q235钢的Johnson-Cook模型参数^[16]

Table  3.  Johnson-Cook parameters of Q235 steel^[16]

ρ/(g·cm-3)	G/GPa	A/GPa	B/GPa	C	n	m	Tm/K	Tr/K
7.83	77	0.792	0.51	0.014	0.26	1.03	1 793	294

下载: 导出CSV 

| 显示表格

2.   模拟结果与分析

2.1   10 mm药厚的模拟结果

2.1.1   碰撞点位移

图 3所示是爆炸焊接结束时复板的竖向位移云图。由图 3可看出，复板的位移大致均为6 mm，表明基、复板已完全复合。为了更直观地观察复板单元位移的变化情况，在复板上选择3个特征单元(431 806、437 359、444 788)，输出其位移-时间曲线，如图 4所示。由图 4可看出，特征单元的竖向位移均略大于间隙(6 mm)，这是由于在爆炸载荷作用下复板有一定程度的减薄率所致。

图  3  10 mm药厚下爆炸焊接结束时复板的z向位移云图

Figure  3.  z-direction displacement contour of flyer plate with explosive thickness of 10 mm at the end of explosive welding

下载: 全尺寸图片 幻灯片

图  4  10 mm药厚下复板上3个特征单元的z向位移-时间历程

Figure  4.  z-direction displacement histories of 3 characteristic elements with explosive thickness of 10 mm

下载: 全尺寸图片 幻灯片

2.1.2   复板碰撞速度

图 5所示是一对分别取自基板与复板结合界面处的特征单元(基板单元：798 751；复板单元：416 251)，特征单元的选取与前期实验^[9]中金相试样的取样位置一致。

图  5  10 mm药厚下的一对特征单元

Figure  5.  A pair of characteristic elements with explosive thickness of 10 mm

下载: 全尺寸图片 幻灯片

图 6所示是这对特征单元的速度-时间曲线。可以看出，基板在碰撞前有一个正的速度峰；该现象的产生如文献[17]所述，是由于爆轰产物不断堆积以及前碰撞点在待复合区产生的振动能所致。复板上所取单元的最大碰撞速度为897 m/s。

图  6  一对特征单元(见图 5)的速度-时间曲线

Figure  6.  Velocity-time curves of the pair of characteristic elements (see Fig. 5)

下载: 全尺寸图片 幻灯片

图 7所示是在复板结合界面处所选取的3个特征单元(410 476、416 251、420 976)。图 8所示是这3个特征单元的速度-时间曲线。

图  7  10 mm药厚下复板结合界面处的3个特征单元

Figure  7.  3 characteristic elements at the bonding interface of flyer plate with explosive thickness of 10 mm

下载: 全尺寸图片 幻灯片

图  8  10 mm药厚下3个特征单元的速度-时间曲线

Figure  8.  Velocity-time curves of 3 characteristic elements with explosive thickness of 10 mm

下载: 全尺寸图片 幻灯片

由图 8可以看出，随着距起爆端距离的增加，复板的碰撞速度增大。由文献[17]的结论可知，该现象是由于基板与复板的碰撞在金属板的待复合区产生了强烈振动引起的。

2.1.3   碰撞点压力分布

图 9所示是在结合界面处选取的3个特征单元(415 576、418 051、419 776)，单元415 576取在复板中心处，与前期实验^[7]中取样做金相观察的位置一致。图 10所示是3个特征单元的压力历程。

图  9  10 mm药厚下复板结合界面处的3个特征单元

Figure  9.  3 characteristic elements at the bonding interface of flyer plate with explosive thickness of 10 mm

下载: 全尺寸图片 幻灯片

图  10  10 mm药厚下3个特征单元的压力-时间曲线

Figure  10.  Pressure-time curves of 3 characteristic elements with explosive thickness of 10 mm

下载: 全尺寸图片 幻灯片

由图 10可以看出，随着距起爆端距离的增加，复板的碰撞压力增大。由文献[17]的结论可知，该现象是爆轰产物不断堆积以及前碰撞点在金属板待复合区振动能不断增加的共同作用结果。

2.2   5 mm药厚的模拟结果

2.2.1   碰撞点位移

图 11所示是爆炸焊接结束时复板的竖向位移云图。由图 11可看出，复板的位移大致均为6 mm，表明基、复板已完全复合。为了更加直观地观察复板单元位移的变化情况，在复板上选择3个特征单元(432 182、438 034、443 960)，输出其位移-时间曲线，如图 12所示。由图 12可看出，特征单元的竖向位移均略大于6 mm，但较10 mm药厚下的竖向位移小。这是由于5 mm药厚下的爆炸载荷作用比10 mm药厚下小，导致5 mm药厚下的复板减薄率比10 mm药厚下低。

图  11  5 mm药厚下爆炸复合结束时复板的z向位移云图

Figure  11.  z-displacement contour of flyer plate with explosive thickness of 5 mm at the end of explosive welding

下载: 全尺寸图片 幻灯片

图  12  5 mm药厚下复板上特征单元的z向位移-时间历程

Figure  12.  z-displacement histories of 3 characteristic elements with explosive thickness of 5 mm

下载: 全尺寸图片 幻灯片

2.2.2   复板碰撞速度

图 13所示是一对分别取自基板与复板结合界面处的特征单元(基板单元：799 201；复板单元：416 701)，特征单元的选取与前期实验^[9]中金相试样的取样位置一致。

图  13  5 mm药厚下的一对特征单元

Figure  13.  A pair of characteristic elements with explosive thickness of 5 mm

下载: 全尺寸图片 幻灯片

图 14所示是这对特征单元的速度-时间曲线, 可以看出，基板在碰撞前也有一个正的速度峰。复板上所取单元的最大碰撞速度为565 m/s。

图  14  一对特征单元(见图 13)的速度-时间曲线

Figure  14.  Velocity-time curves of the pair of characteristic elements (see Fig. 13)

下载: 全尺寸图片 幻灯片

图 15所示是在复板结合界面处所选取的3个特征单元(411 976、417 001、423 826)。图 16所示是这3个特征单元的速度-时间曲线。由图 16可以看出，随着距起爆端距离的增加，复板的碰撞速度增大。

图  15  5 mm药厚下复板结合界面处的3个特征单元

Figure  15.  3 characteristic elements at the bonding interface of flyer plate with explosive thickness of 5 mm

下载: 全尺寸图片 幻灯片

图  16  5 mm药厚下3个特征单元的速度-时间曲线

Figure  16.  Velocity-time curves of 3 characteristic elements with explosive thickness of 5 mm

下载: 全尺寸图片 幻灯片

2.2.3   碰撞点压力分布

图 17所示是在结合界面处选取的3个特征单元(416 326、418 801、422 776)，单元416 326取在复板中心处，与前期实验^[7]中取样做金相观察的位置一致。图 18所示是这3个特征单元的压力历程。由图 18可以看出，随着距起爆端距离的增加，复板的碰撞压力增大。

图  17  5 mm药厚下复板结合界面处的3个特征单元

Figure  17.  3 characteristic elements at the bonding interface of flyer plate with explosive thickness of 5 mm

下载: 全尺寸图片 幻灯片

图  18  5 mm药厚下3个特征单元的压力历程

Figure  18.  Pressure histories of 3 characteristic elements with explosive thickness of 5 mm

下载: 全尺寸图片 幻灯片

2.3   分析与讨论

由图 6可以看出，10 mm药厚下复板的最大碰撞速度为897 m/s。由图 14可以看出，5 mm药厚下复板的最大碰撞速度为565 m/s。利用前期工作^[18]中提到的3种理论公式(Gurney公式、Aziz公式、Deribas公式)计算了复板的碰撞速度，如表 4、表 5所示，并与数值模拟结果进行了比较。由表 4和表 5可以看出：Gurney公式和Aziz公式的计算结果均存在较大的偏差；而由Deribas公式计算的两组结果与数值模拟结果较接近，误差均未超过5%，且与前期实验结果较吻合，证明了SPH-FEM耦合算法的可靠性。

表  4  10 mm药厚下碰撞速度理论计算结果与数值模拟结果的比较

Table  4.  Comparison of collision velocity between theoretical calculation and numerical simulation with explosive thickness of 10 mm

Theoreticalformula	Massfraction	Collision velocity/(m·s-1)		Error/%
		Theoretical calculation[18]	Simulation
Gurney	0.75	1 089	897	-21.0
Aziz	0.75	711	897	20.0
Deribas	0.75	853	897	4.9

下载: 导出CSV 

| 显示表格

表  5  5 mm药厚下碰撞速度理论计算结果与数值模拟结果的比较

Table  5.  Comparison of collision velocity between theoretical calculation and numerical simulation with explosive thickness of 5 mm

Theoreticalformula	Massfraction	Collision velocity/(m·s-1)		Error/%
		Theoretical calculation[18]	Simulation
Gurney	0.45	863	565	-52.7
Aziz	0.45	480	565	15.0
Deribas	0.45	576	565	-1.9

下载: 导出CSV 

| 显示表格

由图 10可以看出，10 mm药厚下复板单元415 576处的碰撞压力为17.08 GPa。由图 18可以看出，5 mm药厚下复板单元416 326处的碰撞压力为11.25 GPa。

Ezra等提出的碰撞压力的计算公式为^[16]

式中：v_{s, 1}、v_{s, 2}分别表示复板、基板的声速，m·s^-1；ρ₁、ρ₂分别表示复板、基板的密度，g·cm^-3；v_p表示复板的碰撞速度，m·s^-1。

结合表 4和表 5中3种理论公式计算得到的碰撞速度，通过(3)式可得到复板的碰撞压力，表 6及表 7为其理论计算值与数值模拟结果的比较。可见：Gurney公式和Aziz公式的计算结果均存在较大的偏差；而由Deribas公式计算的两组结果与数值模拟结果较接近，误差均未超过5%，说明Deribas公式和SPH-FEM耦合方法对双面爆炸焊接具有较好的指导意义。

表  6  10 mm药厚下碰撞压力理论计算结果与数值模拟结果的比较

Table  6.  Comparison of collision pressure betweentheoretical calculation and numerical simulationwith explosive thickness of 10 mm

Theoreticalformula	Collision pressure/GPa		Error/%
	Calculation	Simulation
Gurney	22.08	17.08	-29.3
Aziz	14.42	17.08	15.6
Deribas	17.30	17.08	-1.3

下载: 导出CSV 

| 显示表格

表  7  5 mm药厚下碰撞压力理论计算结果与数值模拟结果的比较

Table  7.  Comparison of collision pressure betweentheoretical calculation and numerical simulationwith explosive thickness of 5 mm

Theoreticalformula	Collision pressure/GPa		Error/%
	Calculation	Simulation
Gurney	17.50	11.25	-55.6
Aziz	9.73	11.25	13.5
Deribas	11.68	11.25	-3.8

下载: 导出CSV 

| 显示表格

3.   结论

利用LS-DYNA软件和SPH-FEM耦合方法对前期双面爆炸焊接实验进行了三维数值模拟，并将模拟结果与实验及理论计算结果进行了对比，得到如下结论。

(1) 10 mm药厚和5 mm药厚下复板位移均略大于间隙值6 mm，这是由于爆轰载荷作用下复板有一定的减薄率所致。

(2) 10 mm药厚下，复板中部的最大碰撞速度为897 m/s，碰撞压力为17.08 GPa；5 mm药厚下，复板中部的最大碰撞速度为565 m/s，碰撞压力为11.25 GPa。通过与3种理论公式(Gurney公式、Aziz公式、Deribas公式)计算得到的碰撞速度进行比较发现，数值模拟结果与Deribas公式的计算结果较接近，误差较小，且与实验结果较吻合，证明了SPH-FEM耦合方法用于双面爆炸复合模拟的有效性，同时Deribas公式和SPH-FEM耦合方法对双面爆炸复合具有较好的指导意义。

(3) 10 mm药厚和5 mm药厚下复板的碰撞速度及碰撞压力均随着距起爆端距离的增加而增大，该现象是由于爆轰产物的不断堆积和前碰撞点在金属板待复合区振动能的不断增加共同作用的结果。