计算域可变的CEL方法*

付 峥，刘 军，冯其京，王 政，张树道

(北京应用物理与计算数学研究所，北京 100094)

计算域可变的CEL方法*

付 峥，刘 军，冯其京，王 政，张树道

(北京应用物理与计算数学研究所，北京 100094)

为了更好地兼顾拉氏方法和欧拉方法各自的特长，提出一种可将拉氏介质映射到欧拉计算域的耦合欧拉-拉格朗日(CEL)方法。通过这种映射，将欧拉-拉格朗日重叠区域的接触面协调问题转换为欧拉区域内的多介质计算问题，简化了CEL方法的构造过程。通过与侵彻实验和结构对爆炸冲击波响应实验的比较，验证了新算法，计算结果与实验数据符合较好。

CEL方法；映射算法；计算域；冲击动力学

基于网格的数值方法是研究冲击动力学问题的重要工具[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]Fig.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)不难发现，新耦合算法的关键是拉氏区向欧拉区的映射过程。在讨论映射算法前，首先要在欧拉程序中为插入的拉氏介质安排合理的材料编号。设在一个计算问题中拉氏区有NL种介质，欧拉区有NE种介质(不包含由拉氏区插入的介质)，在欧拉程序中这些介质不会被删除，故称这类介质为永久介质。由拉氏区插入到欧拉区的有NLE种介质(NLE≤NL)，这些介质是被临时插入到欧拉区的，故称这类介质为临时介质。在欧拉程序中，若编号为M的介质为永久介质，则有1≤M≤NE，若为临时介质，则有NE

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

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

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

(1)

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

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

(2)

式中：VE为欧拉单元的体积；cE为欧拉单元中的声速；eE,M为欧拉单元中的质量内能；φE,void为空介质的体积；φE,M、ρE,M和eE,M分别为欧拉单元上永久介质M(1≤M≤NE)的体积分数、密度和比内能。

(3)

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

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

图2 过盈单元Fig.2 Overfilled cell

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

计算介质M的声速cE,M，得到介质M的体积模量KE,M：

(4)

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

(5)

式中：Nmat=NE+NLE，为单元内的介质总数，VE,M为介质体积。介质M的体积修正值ΔVE,M和压力修正值ΔpE,M分别为：

(6)

(7)

修正后的体积、压力分别为pE,M+ΔpE,M和VE,M+ΔVE,M。修正介质M的内能:

(8)

(f) 若空介质体积VE,void>0，且VE,void+VL+VE>Vcell，说明VE,void偏大，则减小VE,void，使得VE,void+VL+VE=Vcell。

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

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

1.2 积分拉氏介质的边界力

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

(a) 若拉氏单元与一欧拉单元相交，则该欧拉单元被称为欧拉-拉氏混合单元(E/L单元)。首先要计算E/L单元内欧拉-拉氏接触面面积Ao和接触面上欧拉介质对拉氏介质施加的压力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单元压力pmix和前单元压力pfor的加权平均值作为永久介质对拉氏面元的压力(权系数wmix的计算公式如图3所示，本文中wE=1，φE,mix表示E/L单元中永久介质的体积分数之和)。E/L单元的偏应力是永久介质和临时介质加权平均的结果，不能表示永久介质对拉氏面元的作用。这里用前单元的偏应力Sfor作为永久介质在拉氏面元上的偏应力。若计算气固耦合问题，永久介质通常为气体，此时拉氏面元所承受的压力等于pfor。计算结果表明，这种只有一阶精度的计算方法可以满足数值计算的需要(误差小于5%)。

图3 拉氏介质边界上的压力、偏应力和应力Fig.3 Pressure, deviatoric stress and stress on the boundary of the Lagrangian material

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

(9)

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

(10)

式中：fn为fn的大小。

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

(11)

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

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

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

(12)

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

(13)

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

(14)

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

2 算例

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

图4 初始时刻网格示意图Fig.4 Initial grid

2.1 侵彻实验

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

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

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

表1 各工况弹体残余速度

图5 弹体与靶板的变形情况对比Fig.5 Distortion of the bullet and the target

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

图6 实验装置[13]Fig.6 Experimental setup[13]

图7 初始时刻圆盘附近的网格Fig.7 Initial grid around the plate

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

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

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

3 结 论

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

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

[1] 高凌天,刘凯欣,刘颖.无网格方法在冲击动力学中的应用[C]∥第三届全国计算爆炸力学会议.青岛,2006.

[2] 王瑞利.一种物理量重映方法的研究[J].数值计算与计算机应用,2002,23(4):296-302. Wang Ruili. An research in remapping techniques[J]. Numerical Computation and Computer Application, 2002,23(4):296-302.

[3] Noh W F. CEL: A time-dependent two-space-dimensional coupled eulerian-lagrangian code[J]. Methods in Computational Physics, 1964,3:117-179.

[4] Olovsson L. On the arbitrary Lagrangian-Eulerian finite element method[D]. Linköping, Sweden: Linköpings University, 2000.

[5] Hallquist J O. LS-DYNA theoretical manual[M]. Livermore, CA: Livermore Software Technology Corporation, 1998:1-5.

[6] Fedkiw R P. Coupling an Eulerian fluid calculation to a Lagrangian solid calculation with the ghost fluid method[J]. Journal of Computational Physics, 2002,175(1):200-224.

[7] Arienti M, Hung P. A level set approach to Eulerian-Lagrangian coupling[J]. Journal of Computational Physics, 2003,185(1):213-251.

[8] Brown K H, Burns S P, Christon M A. Coupled Eulerian-Lagrangian methods for earth penetrating weapon applications: SAND2002-1014[R]. Office of Scientific and Technical Information Technical Reports, 2002.

[9] 宁建国,王猛.关于计算爆炸力学的进展与现状[J].力学与实践,2012,34(1):10-19. Ning Jianguo, Wang Meng. Review on computational explosion mechanics[J]. Mechanics in Engineering, 2012,34(1):10-19.

[10] 何翔,徐翔云,孙桂娟,等.弹体高速侵彻混凝土的效应实验[J].爆炸与冲击,2010,30(1):1-6. He Xiang, Xu Xiangyun, Sun Guijuan, et al. Experimental investigation on projectiles’ high-velocity penetration into concrete targets[J]. Explosion and Shock Waves, 2010,30(1):1-6.

[11] 刘军,何长江,梁仙红.三维弹塑性流体力学自适应欧拉方法研究[J].高压物理学报,2008,22(1):72-78. Liu Jun, He Changjiang, Liang Xianhong. An Eulerian adaptive mesh refinement method for three dimensional elastic-plastic hydrodynamic simulations[J]. Chinese Journal of High Pressure Physics, 2008,22(1):72-78.

[12] Piekutowski A J, Forrestal M J, Poormon K L, et al. Perforation of aluminum plates with ogive-nose steel rods at normal and oblique impacts[J]. International Journal of Impact Engineering, 1996,18(7/8):877-887.

[13] Neubergera A, Peles S, Rittel D. Scaling the response of circular plates subjected to large and close-range spherical explosions.Part I: Air-blast loading[J]. International Journal of Impact Engineering, 2007,34(5):859-873.

(责任编辑 王小飞)

A CEL method with changeable computational domain

Fu Zheng, Liu Jun, Feng Qijing, Wang Zheng, Zhang Shudao

(InstituteofAppliedPhysicsandComputationalMathematics,Beijing100094,China)

In order to improve the large deformation problem within the Lagrangian domain, a CEL method with the capability of mapping Lagrangian materials to the Eulerian domain was presented. The contact problem in the Eulerian-Lagrangian overlapping region was converted to the multi-materials problem in the Eulerian region. The construction of the CEL method was also simplified. The method was verified by the calculation of two experiments (a steel bullet impacting an aluminous plate experiment and a structure response to the blast wave experiment). It is found that the numerical results agree well with the experimental data.

CEL method; mapping algorithm; computational domain; impact dynamics

10.11883/1001-1455(2017)03-0528-08

2015-07-03；

2015-10-23

国家自然科学基金项目(11472060，11371069)； 中国工程物理研究院科学技术发展基金项目(2014B0201030)

付 峥(1984- )，男，助理研究员，fu_zheng@iapcm.ac.cn。

O384 国标学科代码： 13035

A