基于网格的数值方法是研究冲击动力学问题的重要工具^[1]。其中，拉格朗日方法能够清晰描述物质界面的演化，但在处理大变形问题时网格会严重扭曲、变形，如不进行特殊处理，计算精度会受到严重影响，甚至得到非物理解或迫使计算终止；欧拉方法不会出现网格扭曲的问题，但在进行多介质问题计算时，如不进行特殊处理，物质界面会逐渐模糊，甚至偏离物理解^[2]。为了发挥拉氏和欧式方法的优点，1964年，W.F.Noh^[3]提出了耦合欧拉-拉格朗日(coupled Euler Lagrange, CEL)的概念。他将整个计算域划分成若干子域，一些子域划入拉格朗日计算域，其他子域划入欧拉计算域，通过协调两种计算域接触面上的力、速度等物理量实现耦合计算。这一思想奠定了CEL方法的基础。此后，CEL方法迅速发展，新的耦合算法层出不穷。其中，有代表性的算法包括：L.Olovsson^[4]结合LS-DANY的接触算法^[5]提出的一种罚函数算法(Penalty Method)，R.P.Fedkiw^[6]和M.Arienti等^[7]结合Level Set和Ghost Fluid方法提出的Ghost Fluid Eulerian Lagrangian (GEL)方法。目前CEL方法的思想已经在ABAQUS，ANSYS等大型商用软件中得到实现，并在侵彻^[8]、固体结构对爆炸冲击波响应^[9]等冲击动力学问题的研究中得到推广和应用。

但是，目前大部分CEL算法在耦合计算时，欧拉和拉格朗日计算域都是预先设定好的，随着流场的发展，拉格朗日计算域内仍可能出现网格扭曲甚至自相交等问题，拉氏方法原有的缺点依然存在，特别是类似材料发生拓扑变化，产生新的物质界面等拉氏方法处理困难的问题，此时预设区域的局限性就会凸现出来。例如，计算弹体侵彻问题时，原始CEL方法预先将弹体划入拉格朗日计算域，将靶板划入欧拉计算域。当入射速度较低时，弹体变形较小，预设区域的耦合算法能够得到较好效果^[8]。但当弹速较高，发生高速甚至超高速侵彻的情况下，弹头将快速钝化、损耗，磨损掉的部分与靶板粘结在一起导致弹体质量不断下降，直至侵彻结束(如图 1^[10])，此时预设区域的耦合算法将无法计算。为了模拟拉格朗日计算域内存在网格扭曲或自相交等问题的冲击动力学现象，需要开发一种新型的CEL耦合算法。

图  1  高速侵彻实验回收的弹体^[10]

Figure  1.  Recovered bullets in experiments of high velocity impact^[10]

下载: 全尺寸图片 幻灯片

本文中，提出一种将拉氏计算域映射到欧拉计算域的紧耦合CEL算法。通过这种映射，将欧拉-拉格朗日接触面上界面和物理量的协调问题转换为欧拉区域内的多介质计算问题。在映射过程结束后，从欧拉计算域积分得到欧拉-拉氏接触面上拉氏介质的边界力，完成耦合计算。采用新算法计算一个侵彻实验和一个结构对爆炸载荷的响应实验，通过与实验数据对比验证算法的有效性。

1.   CEL耦合方法

进行CEL耦合计算时，在欧拉-拉氏重叠区内，拉氏介质被插入到欧拉区域，为欧拉区域提供速度、位移等边界条件。在欧拉-拉氏接触面上，对欧拉介质进行应力积分，为拉氏计算域提供力边界条件。若拉氏域内的介质产生大变形，导致拉氏网格单元严重扭曲或自相交，则将该单元及其附近区域的拉氏介质转化为欧拉计算域内的介质。为了实现上述思想，编制支持将拉氏计算域映射到欧拉区域的CEL程序，程序的欧拉模块采用的是多介质弹塑性程序MEPH^[11]，拉氏模块采用的是简单、流行的有限元程序。CEL耦合算法的流程可简单描述如下。

(1) 清除前一步映射到欧拉计算域的拉格朗日介质。

(2) 将拉格朗日单元分解成四面体单元。

(3) 将拉格朗日介质插入欧拉网格：

      (a) 计算欧拉-拉氏重叠区域的体积；

      (b) 从拉格朗日计算模块获得拉氏介质物理量的分布情况；

      (c) 将获得的拉氏介质的物理量插入欧拉网格；

      (d) 将映射后的数据传给欧拉模块。

(4) 积分欧拉介质对拉氏介质的作用力：

      (a) 计算欧拉-拉氏接触面在每个欧拉单元内的面积；

      (b) 积分得到接触面上欧拉介质施加给拉氏单元的边界力；

      (c) 分解(b)步得到的力，获得拉氏单元所受的法向力；

      (d) 如果存在摩擦，采用Coulomb模型计算拉氏单元所受的摩擦力；

      (e) 将(a)~(d)算得的力重新分配到拉氏单元的结点，并传给拉格朗日模块。

(5) 统一欧拉和拉格朗日模块的时间步长，分别运行欧拉和拉格朗日模块。

(6) 进入下一时间步，从(1)重新开始。

从(1)~(6)不难发现，新耦合算法的关键是拉氏区向欧拉区的映射过程。在讨论映射算法前，首先要在欧拉程序中为插入的拉氏介质安排合理的材料编号。设在一个计算问题中拉氏区有N_L种介质，欧拉区有N_E种介质(不包含由拉氏区插入的介质)，在欧拉程序中这些介质不会被删除，故称这类介质为永久介质。由拉氏区插入到欧拉区的有N_LE种介质(N_LE≤N_L)，这些介质是被临时插入到欧拉区的，故称这类介质为临时介质。在欧拉程序中，若编号为M的介质为永久介质，则有1≤M≤N_E，若为临时介质，则有N_E<M≤N_E+N_LE。采用这样的编号系统可以方便的插入或删除临时介质，也容易区分永久介质和临时介质。

1.1   拉氏介质到欧拉区的映射

拉氏区到欧氏区的映射过程，实际上就是将拉氏介质插入到欧拉网格的过程。

(a) 设某拉氏单元与一欧拉单元相交，重叠区体积为V_o。设拉氏单元内密度均匀，为该单元的平均密度ρ_l，则该拉氏单元内介质应插入到这一欧拉单元的各物理量为：

式中：m、c、p、S、P、E、V分别表示质量、声速、压力、偏应力、动量、内能和体积，下标o、L和ref分别表示重叠区物理量、拉氏单元的物理量以及欧拉区的参考值，上标*表示附加了权值的物理量。

(b) 对所有欧拉-拉氏重叠区内的拉氏单元重复(a)。

(c) 计算永久介质M(1≤M≤N_E)的体积V_{E, M}、质量m_{E, M}、总内能E_{E, M}，空介质的体积V_{E, void}以及插入临时介质前欧拉单元中的加权声速c_E^*、加权偏应力p_E^*和永久介质M_{E, M}^*的加权压力：

式中：V_E为欧拉单元的体积；c_E为欧拉单元中的声速；e_{E, M}为欧拉单元中的质量内能；ϕ_{E, void}为空介质的体积；ϕ_{E, M}、ρ_{E, M}和e_{E, M}分别为欧拉单元上永久介质M(1≤M≤N_E)的体积分数、密度和比内能。

(d) 将(b)得到的数据插入到欧拉网格，设该拉氏介质在欧拉程序中对应编号为M(N_E < M≤M_E+N_LE)的临时介质，则临时介质M的体积V_{E, M}、质量m_{E, M}、加权压力p_{E, M}^*和总内能E_{E, M}，以及欧拉单元的加权声速c_E^*、加权偏应力S_E^*和动量P_E为：

式中：上角标“′”代表赋值迭代后的量。

(e) 在计算欧拉单元的声速、偏应力等物理量前，首先要考虑体积可忽略的介质和过盈单元。

若某介质体积可忽略，则用空介质替换该介质。若所有介质体积分数之和大于1(过盈单元，如图 2所示，此时V_L+V_E>V_cell，则V_{E, void}=V_cell-V_L-V_E < 0，不符合V_{E, void}≥0的物理规律，必须修正)，根据单元内介质的分布情况分类处理。若过盈单元仅含有永久介质或临时介质，在保障密度、比内能和压力不变的前提下，按介质的体积比修正各介质的体积、质量、内能和加权压力，使得V_{E, void}=0；若过盈单元同时含有永久介质和临时介质，根据下面的算法修正介质的体积、内能等物理量。

图  2  过盈单元

Figure  2.  Overfilled cell

下载: 全尺寸图片 幻灯片

计算介质M的声速c_{E, M}，得到介质M的体积模量K_{E, M}：

进一步计算介质M的体积修正因子:

式中：N_mat=N_E+N_LE，为单元内的介质总数，V_{E, M}为介质体积。介质M的体积修正值ΔV_{E, M}和压力修正值Δp_{E, M}分别为：

修正后的体积、压力分别为p_{E, M}+Δp_{E, M}和V_{E, M}+ΔV_{E, M}。修正介质M的内能:

式中：E′_{E, M}为E_{E, M}的修正值。重复(c)、(d)。

(f) 若空介质体积V_{E, void}>0，且V_{E, void}+V_L+V_E>V_cell，说明V_{E, void}偏大，则减小V_{E, void}，使得V_{E, void}+V_L+V_E=V_cell。

(g) 由式(2)~(3)中的加权声速、加权压力、加权偏应力、内能和体积计算声速、压力、偏应力、比内能和体积分数。将结果传回欧拉程序，结束映射过程。

CEL方法中，欧拉-拉氏接触面上，拉氏介质的物质界面很难与欧拉介质的物质界面重合。两种描述下的物质界面匹配问题是CEL方法的一大难题。而上述映射过程是基于介质体积的，相对独立于物质界面，降低了CEL方法对物质界面一致性的要求。

1.2   积分拉氏介质的边界力

欧氏区到拉氏区的映射，计算的是欧拉-拉氏接触面上，欧拉介质对拉氏介质施加的作用力。

(a) 若拉氏单元与一欧拉单元相交，则该欧拉单元被称为欧拉-拉氏混合单元(E/L单元)。首先要计算E/L单元内欧拉-拉氏接触面面积A_o和接触面上欧拉介质对拉氏介质施加的压力p和偏应力S。图 3阐述了p和S的计算过程，其中I为单位张量。图 3中拉氏介质覆盖了三个E/L单元e(i, j+1), e(i+1, j+1)和e(i+1, j+2)。n是E/L单元e(i+1, j+1)内拉氏面元的外法向。沿n方向外推一个欧拉单元的距离，可得仅有永久介质的欧拉单元-前(forward)单元(单元e(i+2, j))。本文中采用E/L单元压力p_mix和前单元压力p_for的加权平均值作为永久介质对拉氏面元的压力(权系数w_mix的计算公式如图 3所示，本文中w_E=1，φ_{E, mix}表示E/L单元中永久介质的体积分数之和)。E/L单元的偏应力是永久介质和临时介质加权平均的结果，不能表示永久介质对拉氏面元的作用。这里用前单元的偏应力S_for作为永久介质在拉氏面元上的偏应力。若计算气固耦合问题，永久介质通常为气体，此时拉氏面元所承受的压力等于p_for。计算结果表明，这种只有一阶精度的计算方法可以满足数值计算的需要(误差小于5%)。

图  3  拉氏介质边界上的压力、偏应力和应力

Figure  3.  Pressure, deviatoric stress and stress on the boundary of the Lagrangian material

下载: 全尺寸图片 幻灯片

(b) 计算应力σ的法向分量:

(c) 用法向应力分量计算拉氏面元受的法向力f_n:

式中：f_n为f_n的大小。

(d) 如果考虑摩擦，则计算摩擦力f_t:

式中：f、μ、m、Δt分别表示面元切向量、摩擦系数、面元平均质量和时间步长；f_t为f_t的大小；v_t为切向速度；v_r表示接触面上，拉氏与欧氏介质的相对速度。当拉氏介质侵入欧拉介质时，v_r=-v_L，v_L为拉氏面元的质心速度。这是由于此时拉氏介质比它周围的欧拉介质硬，而拉氏介质的表面更精确。当欧拉介质侵入拉氏介质时，v_r=v_E-v_L，v_E为E/L单元的形心速度。

(e) 至此可算得拉氏面元受永久介质的力f=f_n+f_t，此时f作用在面元的形心。将f重新分配到面元的结点，每个结点分得1/n的力(n为面元的结点数)。

(f) 对所有E/L单元内欧拉-拉氏接触面上的拉氏面元重复步骤(a)~(e)。

(g) 若考虑摩擦，还要计算E/L单元内永久介质的摩擦力f_{t, avg}并修正永久介质的应力状态。首先定义E/L单元内临时介质的界面法向n_avg：

式中：N_L为与一个E/L单元相交的拉氏单元的数量，j表示这些拉氏单元的编号，n_j表示编号为j的拉氏单元的外法向，A_{E, j}为编号为j的拉氏面元在该E/L单元内的面积。由此可知，永久介质的界面法向为-n_avg。进一步可由式(11)计算永久介质的界面切向量t_avg。定义摩擦力f_{t, avg}:

式中：f_{t, j}表示编号为j的拉氏单元所受的摩擦力。然后可得剪应力分量τ_avg：

最后根据n_avg将τ_avg从局部坐标系旋转到全局坐标系即可。

2.   算例

采用冲击动力学问题中常见的两个问题验证提出的CEL算法。算例中，金属材料采用Mie-Grüneisen状态方程和Johnson-Cook弹塑性模型，空气采用理想气体状态方程，炸药采用JWL状态方程。

2.1   侵彻实验

弹体侵彻是冲击动力学领域的基础问题之一。这里选用Piekutowski等^[12]的长杆弹侵彻实验。其中长杆弹为弹长88.9mm、直径12.9mm、弹头长21.4mm的钢制卵形弹体，靶板为半径304mm、厚26.3mm的铝板。本文选用了弹体与靶板正碰的五种工况，弹体的入射速度依次为341、396、508、730和863m/s。侵彻过程中，子弹形状几乎未发生变化，采用拉格朗日方法计算。靶板变形巨大并被穿透，采用欧拉方法计算。

初始时刻，拉氏与欧氏网格如图 4所示。随着子弹入射速度的不断增加，穿靶后子弹的剩余速度也在不断变化。表 1比较了实验数据与CEL方法的计算结果。参考表 1，入射速度相对较低时(341m/s)，计算结果的误差较大，超过了5%；而其他工况的计算结果均小于2%。这是计算低速侵彻问题时所用的欧拉网格分辨率不够所致，单独用不同的欧拉网格步长对工况1(入射速度为341m/s)进行计算，欧拉网格最小单元尺寸为2、1、0.5、0.25mm时，相对误差分别为40.74%、5.79%、2.42%、1.81%，这说明在实际计算时通过缩小欧拉单元的尺寸可以减小误差。

图  4  初始时刻网格示意图

Figure  4.  Initial grid

下载: 全尺寸图片 幻灯片

表  1  各工况弹体残余速度

Table  1.  Residual velocity of the bullet

工况	入射速度/(m·s-1)	剩余速度/(m·s-1)		相对误差/%
		实验	计算
1	341	164	173.5	5.79
2	396	266	270	1.50
3	508	415	422	1.69
4	730	665	674	1.35
5	863	802	815	1.62

下载: 导出CSV 

| 显示表格

穿靶后的靶板形状是实验提供的另一重要参考数据。图 5比较了Piekutowski等给出的工况2~5的照片(左侧黑白图片)与相近时刻CEL方法的计算结果(右侧彩色图片)。图中的计算结果与实验图片吻合较好。

图  5  弹体与靶板的变形情况对比

Figure  5.  Distortion of the bullet and the target

下载: 全尺寸图片 幻灯片

2.2   金属圆盘对爆炸冲击波的响应

薄壁结构对爆炸冲击波的响应也是冲击动力学领域的一个基础问题。这里选用A.Neuberger等^[13]的实验，实验设备如图 6所示。在质量为W的TNT炸药起爆后，冲击波在空气中传播并迅速到达距离炸药质心R处的钢制圆盘(盘厚d，直径D)。装置底端设置一梳子形状的测量装置，用来测量圆盘的最大位移δ。由于炸药产物为气体，变形大，采用欧拉方法描述。而薄壁结构厚度小，若采用欧拉方法需要数量庞大的精细单元，浪费大量计算资源，因此采用拉氏方法描述。

图  6  实验装置^[13]

Figure  6.  Experimental setup^[13]

下载: 全尺寸图片 幻灯片

初始时刻，拉氏与欧氏网格如图 7所示。表 2列出了实验中4种工况下各参数的取值，并比较了实验测量的圆盘最大位移与CEL方法、文献[13]的计算结果。所有计算结果与实验数据符合较好。位移测量装置是一个梳子状的金属结构，该装置测量的位移实际上是结构产生的塑性形变，未包含结构的弹性形变部分。因此，实验测量的最大位移应比真实结果小。表 2中CEL方法计算得到的最大位移均大于实验测量值，更好地反映了实际情况。而文献[13]的计算结果没有此规律。说明CEL方法在模拟气体-结构相互作用的问题时仍然能够保障足够的计算精度。

图  7  初始时刻圆盘附近的网格

Figure  7.  Initial grid around the plate

下载: 全尺寸图片 幻灯片

表  2  各工况实验与数值模拟得到的圆盘最大位移

Table  2.  Experimental and numerical maximal displacement of plates

工况	d/cm	D/cm	W/kg	R/cm	δ/d
					Exp.[13]	LS-DYNA[13]	CEL
1	1	50	1.094	10	4.85	4.98	4.93
2	2	100	8.752	20	5.35	5.15	5.49
3	1	50	1.094	6.5	7.45	7.72	7.67
4	2	100	8.752	13	8.25	8.15	8.29

下载: 导出CSV 

| 显示表格

3.   结论

提出一种可将拉氏计算域映射到欧拉计算域的CEL方法。并使用该算法模拟了两类典型的冲击动力学问题：弹体侵彻和结构对爆炸载荷的响应。通过与实验数据比较发现，本文提出的算法得到的计算结果与实验数据符合较好，能够反映固体结构的真实力学行为。

下一步，将对CEL方法进行更加系统的测试，并实现将拉氏介质的部分区域映射到欧拉网格。