拉氏方法模拟二维多介质可压缩流体的运动

赵海波; 肖波; 柏劲松; 段书超; 王刚华; 阚明先; 陈芳

doi:10.11858/gywlxb.20170694

拉氏方法模拟二维多介质可压缩流体的运动

doi: 10.11858/gywlxb.20170694

中国工程物理研究院流体物理研究所, 四川绵阳 621999

基金项目:

科学挑战计划 TZ2016001

国家自然科学基金 11405167

国家自然科学基金 11571293

国家自然科学基金 11672276

国家自然科学基金 11532012

中国工程物理研究院科学技术发展基金 2015B0201023

详细信息

作者简介:
赵海波(1994-), 男, 硕士研究生, 主要从事计算流体力学研究.E-mail:18310737226@163.com

通讯作者:
柏劲松(1968-), 男, 研究员, 博士生导师, 主要从事计算力学研究

中图分类号: O354
计量
- 文章访问数: 7016
- HTML全文浏览量: 2867
- PDF下载量: 257
出版历程
- 收稿日期: 2017-12-21
- 修回日期: 2018-01-09

Simulation of Two-Dimensional Multi-Material Compressible Flows Using Lagrangian Methods

Institute of Fluid Physics, CAEP, Mianyang 621999, China

摘要

摘要: 在用拉格朗日方法模拟二维多介质可压缩流体的运动时，网格发生大变形往往是模拟不能正常进行下去的重要原因。网格动态局域重分可以有效改善网格的畸变程度，使计算得以持续。针对三角形计算网格提出了一种新的动态局域重分方法，包含"对角线交换""长边劈裂""短边融合"和"帽子戏法"4种基本操作，其中前3种操作不仅作用于同种介质内部，还可将其拓展到多介质界面处，与"帽子戏法"一起处理界面附近的大变形网格。在网格动态局域重分后，将旧网格上的物理量映射到新网格上，先计算出新三角形的质量和内能，再根据动量守恒和能量守恒对新三角形的格点速度及内能进行修正。利用该方法对冲击波与气泡相互作用和R-T不稳定性问题进行了数值模拟，取得了良好的效果。在R-T不稳定性算例中，采用同种介质和不同介质两种模型进行对比，模拟结果验证了该方法的有效性。
- 多介质流体 /
- 大变形运动 /
- 交错网格 /
- 网格动态局域重分 /
- 守恒重映
Abstract: When simulating the two-dimensional multi-material compressible flows using the Lagrangian method, the distortion of the mesh is often the reason why the process terminates.The dynamic local remeshing is an option to relax the distortion.In this study, we combined the "diagonal-swapping" "edge-splitting" "edge-merging" and "hat-trick" to deal with the distortive mesh.The 4 operations could act not only in single material, but also be extended to the multi-material interface to deal with the distortive mesh on the interface.After the remeshing, the quantities on the old mesh are projected onto the new one.During this process, the mass, momentum and energy are kept conservative.We verified the effectiveness and accuracy of this method by using it to simulate the shock/bubble interaction and the R-T instability.In the R-T instability simulation, the single material and the multi-material models are applied to prove that this method is able to deal with the distortive mesh both inside the material and on the interface.
- multi-material flow /
- large deformation motion /
- staggered grid /
- dynamic local remeshing /
- conservative remapping

HTML全文

采用拉氏方法模拟多介质可压缩流体运动，可以避免因为输运项带来的数值耗散，并且可以精确描述物质界面。但是当流体发生大变形时，拉氏方法将面临网格畸变的情况，导致模拟难以进行下去。当前模拟大变形多介质可压缩流体运动较常见的方法是在一般的任意拉格朗日-欧拉(简称ALE)方法^[1]基础上发展起来的多介质ALE方法^[2-5]。该方法通常包括3个步骤：首先, 利用拉氏方法计算下一时刻的物理量；然后, 调整网格格点坐标，改善网格单元品质；最后将旧网格上的物理量映射到新网格上去。在多介质ALE方法中，当物质界面发生极端变形(例如界面翻卷形成漩涡)时，为了保持界面处的网格单元品质，会允许界面处的网格发生跨物质运动，同时需要引入混合物质网格单元描述，此时该方法与欧拉方法已无太大区别。欧拉方法中常用的混合网格界面追踪方法如VOF(Volume of Fluid)^[6]或MOF(Moment-of-Fluid)^[7]等，也因此常常被用于多介质ALE方法以实现界面重构。总的来说，多介质ALE方法的算法复杂度较高，实现起来并不容易。本研究将讨论在纯拉氏方法中引入网格动态局域重分以消除畸变网格，并在重映过程中避免出现混合物质网格，实现对多介质大变形流体运动过程的纯拉氏模拟。

网格动态局域重分方法是近来人们研究较多的一种方法，国内外学者提出了各种网格动态局域重分算法^[8-10]来处理大变形网格或者实现网格自适应，这些算法一般会改变格点之间的连接方式，但不改变格点的位置，文献中将这类方法称为h-remeshing。在h-remeshing中，首先需要定义“质量函数”，即描述网格单元的几何品质，通常从网格的大小和形状两个角度去设置质量函数。然后，根据质量函数给出大变形网格或者需要处理网格的判断条件。最后根据判断条件给出具体的操作过程，如：针对三角形网格，通常包括“对角线交换”“长边劈裂”和“短边融合”^[11]；针对四面体网格，包括“空间翻转”“删除边长”和“插入格点”等^[9]。目前文献上关于网格动态局域重分基本都是针对同种介质内部的情形，将该方法拓展到多物质界面处会面临一些新的难点，包括如何尽可能地避免对原物质界面进行修改，在不得已修改物质界面的情况下如何处理横跨原物质界面的网格单元，等等。

本研究提出了一种应用于多介质流体模拟的动态局域重分方法，它包含4种基本操作：“对角线交换”“长边劈裂”“短边融合”和“帽子戏法”。该重分方法不仅作用于单介质内部而且可以作用于多物质界面处。网格重分完成后，则需要将旧网格上的物理量映射到新网格上。动态局域重分方法的物理量重映比一般的网格重分的物理量重映其实更简单，只需要对4种网格重分操作完成重映即可。对于在物质界面处产生的横跨多物质的新网格单元，本研究中采用一种单物质描述。在这种描述下，不可避免地需要在网格单元之间进行物质迁移。为了尽可能降低人为的物质迁移所带来的数值误差，对物质界面处的网格重分采用更严格的重分判据，降低操作次数，并且引入“帽子戏法”操作以减少物质迁移。基于以上方法，对冲击波与气泡相互作用和Rayleigh-Taylor(R-T)不稳定现象两个算例进行数值模拟，以此验证了该方法模拟多介质大变形运动的有效性。

1. 交错网格型拉氏方法的离散格式

计算网格为三角形网格，采用交错网格型拉格朗日方法，即：速度定义在网格格点上，质量和内能定义在网格中心，如图 1所示。离散格式与文献[10]中的一致，具体为

图 1 交错网格型拉格朗日方法^[10]

Figure 1. Staggered grid hydrodynamics Lagrangian discretization

下载: 全尺寸图片幻灯片

${m_c} = {\rm{con}}\;{\rm{st}}$

(1)

${m_p} = \frac{1}{3}\sum\limits_{c \in C\left( p \right)} {{m_c}}$

(2)

$\mathit{\boldsymbol{x}}_p^{n + 1} = \mathit{\boldsymbol{x}}_p^n + \Delta t\mathit{\boldsymbol{u}}_p^n + \frac{1}{2}{\left( {\Delta t} \right)^2}\frac{{\sum\limits_{c \in C\left( p \right)} {\mathit{\boldsymbol{f}}_{pc}^n} }}{{{m_p}}}$

(3)

$\mathit{\boldsymbol{u}}_p^{n + 1} = \mathit{\boldsymbol{u}}_p^n + \Delta t\frac{{\sum\limits_{c \in C\left( p \right)} {\mathit{\boldsymbol{f}}_{pc}^n} }}{{{m_p}}}$

(4)

$E_c^{n + 1} = E_c^n - \sum\limits_{c \in P\left( c \right)} {\mathit{\boldsymbol{f}}_{pc}^n \cdot \left( {\mathit{\boldsymbol{x}}_p^{n + 1} - \mathit{\boldsymbol{x}}_p^n} \right)}$

(5)

式中：m_c为三角形单元c的质量；m_p为格点p的质量；C(p)是包含该格点的所有三角形的集合；P(c)是属于三角形c的格点的集合；x、u、Δt分别为格点的位移、速度和时间步长；E_c为三角形单元c的内能；n和n+1分别表示第n和第n+1时刻；f_pc表示网格单元c对其格点p的作用力，具体表达式为

${\mathit{\boldsymbol{f}}_{pc}} = \left( {{p_c} + {q_c}} \right)\left( {\begin{array}{*{20}{c}} { - \frac{{{y_{{p^ - }\left( c \right)}} - {y_{{p^ + }\left( c \right)}}}}{2}}\\ {\frac{{{x_{{p^ - }\left( c \right)}} - {x_{{p^ + }\left( c \right)}}}}{2}} \end{array}} \right)$

(6)

式中：p_c表示网格单元的压强值；q_c表示作用在网格上的黏性力大小，具体可以写为

${q_c} = \rho {c_{\rm{v}}}\dot v$

(7)

式中：ρ为网格单元的密度，c_v为黏性系数， $\dot v$ 为体积相对变化率。

上述(1)式~(5)式是拉氏方法的离散格式：(1)式表示网格单元质量不随时间变化；(2)式是格点质量与包含该格点的网格单元质量的关系，本研究中假定每个网格单元将质量均分给3个格点；(3)式、(4)式决定了格点在Δt时间内的位移变化和速度变化；(5)式决定了网格单元在Δt时间内的内能变化。

每个单元对应的计算时间步长的表达式为

$\Delta t = \frac{{{C_{{\rm{safc}}}}{h_{\min }}}}{{{v_{\rm{s}}}}}$

(8)

式中：C_safe是时间因子，取为0.05；h_min是三角形单元中最短高的长度，v_s是网格单元的声速大小。最终选取所有网格单元时间步长中的最小值作为整体的Δt。

2. 网格动态局域重分算法

采用网格动态局域重分处理在拉氏模拟中遇到的大变形网格，具体包含4种基本操作：“对角线交换”“长边劈裂”“短边融合”“帽子戏法”。在执行过程中，会根据不同的判断标准选择相应的网格操作。

2.1 周期性网格

模拟过程中采用周期性网格。所谓周期性网格，是指物理边界处网格的空邻域被其相对物理边界处的网格所替换。如图 2所示，三角形△₁的原空邻域被其相对物理边界处的三角形△₂所替换。采用周期性网格将使每个单元不出现空邻域，使程序代码变得更简洁。

图 2 周期性网格

Figure 2. Periodical mesh

下载: 全尺寸图片幻灯片

2.2 网格重分基本操作

针对多介质流体运动的模拟，不涉及滑移现象和滑移算法，即：在界面处两种介质共用一层格点。在模拟的过程中，界面有可能发生大变形，导致模拟进行不下去，此时需要对界面处的大变形网格进行特殊处理。本研究中，对横跨原物质界面的网格采用单介质描述。此种处理有可能导致界面形状改变以及物质局域的人为迁移，我们将其视为一种让步，是保证不出现“混合网格”所需要付出的代价。为此，在算法中尽量限制了对物质界面处网格的处理次数，以降低方法所带来的误差。下面介绍详细的操作流程。

2.2.1 对角线交换

以三角形是否存在大于给定角度值的内角作为判断标准进行对角线交换操作。如图 3(a)所示，当△ABC的最大角∠ABC大于某个值(本研究取2π/3)时，则找到与△ABC的最长边AC相邻的△ACD，在这两个三角形组成的四边形ABCD中进行对角线交换操作。图 3(a)中，△ABC与△ACD包含相同的介质。如果两个三角形包含不同的介质，即AC是物质界面时，如果发生了对角线交换，则操作如图 3(b)所示。在这个过程中，有可能带来物质界面形状的变化，即界面形状由AC变为ADC，但是基于以下两点，可以将该影响降低。

图 3 对角线交换操作

Figure 3. Diagonal-swapping

下载: 全尺寸图片幻灯片

(1) 将该操作的要求提高，比如当∠ADC大于更大的值(本研究取3π/4)时才进行对角线交换，可以从两方面降低该操作对界面形状的影响：一是提高判断标准，则执行次数将减少；二是∠ADC越大，则ADC越接近AC，物质界面形状的改变越小。

(2) 增加网格数目，则单个网格描述的界面长度减少，从而降低界面操作对界面形状的影响。

2.2.2 长边劈裂

该操作发生的条件是三角形的最长边大于某个给定的长度。如图 4所示，当△ABC的最长边AC的长度大于给定值L₁(本研究取为2L, 其中L是初始时刻网格边长的平均长度)时，则在AC边的中点插入E点，将原来的2个三角形分裂为4个三角形。该操作不涉及界面形状的改变。当该操作发生在物质界面处时，则判断标准更加严格，即L₁值降低(本研究取为L)，在一定情况下使描述界面的格点数目增多。

图 4 长边劈裂操作

Figure 4. Edge-splitting

下载: 全尺寸图片幻灯片

2.2.3 短边融合

与“长边劈裂”相反，“短边融合”操作处理的是由于三角形边长过短而带来的畸变问题。在程序中，计算出三角形最短边的长度，再将最短边的长度与给定的长度值L₂(本研究取为0.5L)比较，如果最短边的长度小于L₂，则删除该边及相应的三角形并且改变格点之间的连接方式。如图 5(a)所示，当AF边小于L₂时，则删除AF边，相应的△AEF和△ABF也会被删除掉，原来与F点相连接的格点重新与A点连接以组成新的三角形。以上操作针对的是同种材料。当被删除格点周围的三角形材料不同时，操作如图 5(b)所示。该操作也可能带来物质界面形状的变化。与“对角线交换”操作类似，可以通过提高该操作的判断条件和加密网格将其影响降低。该过程可以看作对界面的一种“平滑”处理，网格越密，操作带来的误差越低。

图 5 短边融合操作

Figure 5. Edge-merging

下载: 全尺寸图片幻灯片

2.2.4 帽子戏法

该操作处理的是位于界面处的三角形网格。如图 6所示，AD是物质界面，△ABF、△BCF、△CDF和△ADF由同种材料构成，△ADE由另外一种材料构成。在该操作中，同样是以三角形是否存在大于某角度值的内角作为发生条件，但是与多物质界面处的对角线交换发生条件略有不同，比如当∠AFD大于某个角度(本研究取3π/4)且F点周围的三角形由同种介质构成时，将F点移动到AD边的中点G，同时连接EG。这样就在一定程度上改善了界面处三角形网格品质而不引起界面形状的改变。引入该操作的目的是为了在多物质界面处尽可能代替“对角线交换”，从而减少物质迁移带来的误差。

图 6 帽子戏法操作

Figure 6. Hat-trick on the interface

下载: 全尺寸图片幻灯片

2.3 网格重分操作流程

采用如图 7所示的流程将以上4种操作组织起来，数值模拟证明该流程可以有效地处理大变形网格。具体为：

图 7 网格动态局域重分流程图

Figure 7. Main flow of the remeshing algorithm

下载: 全尺寸图片幻灯片

(1) 利用三角形的角度和介质种类作为判断条件，执行一轮“对角线交换”和“帽子戏法”；

(2) 利用三角形的边长作为判断条件，执行一轮“长边劈裂”；

(3) 再次执行多次“对角线交换”和“帽子戏法”，以改善(2)中新生成网格的品质；

(4) 根据三角形的边长判断是否需要执行“短边融合”，如需要，则执行完后再次执行多次“对角线交换”和“帽子戏法”，然后返回步骤(2)；否则该流程结束。

在程序实现过程中，“对角线交换”和“帽子戏法”放在方法SwapWithHatTrickForOneTurn(伪代码见附加材料1)中，“长边劈裂”放在方法SplitForOneTurn(伪代码见附加材料2)中，“短边融合”放在方法MergeForOneTurn(伪代码见附加材料3)中。

3. 物理量重映

网格重分后，需要将旧网格上的物理量重新映射到新网格上去。在重映过程中，需要保证物理量守恒，即质量、动量、能量守恒。先求出新三角形的质量和内能，再根据动量和能量守恒对新三角形的格点速度及内能做一个修正。在重映的过程中，假设网格单元的质量密度和比内能为常数分布。在以上4种网格基本操作中：“长边劈裂”和“帽子戏法”的物理量重映在网格内统一处理，不需要区分单介质内或者物质界面处；而“对角线交换”和“短边融合”则需要做此区分，因为它们涉及到物质是否迁移。

3.1 “对角线交换”后的物理量守恒重映

如图 3所示，执行完“对角线交换”后，四边形ABCD的对角线由AC变为BD，旧三角形△ABC和△ACD变为新三角形△ABD和△BCD。接下来进行物理量重映，分两种情况处理。

(1) △ABC与△ACD包含相同的介质，如图 3(a)所示。

首先进行质量和内能重映，假设△ABC和△ACD的质量密度分别为为ρ_ABC和ρ_ACD，比内能为e_ABC和e_ACD，面积为s_ABC和s_ACD，则新△ABD的质量m_ABD和内能E_ABD分别为

${m_{ABD}} = {m_{ABG}} + {m_{ADG}} = {\rho _{ABC}}{s_{ABG}} + {\rho _{ACD}}{s_{ADG}}$

(9)

${E_{ABD}} = {E_{ABG}} + {E_{ADG}} = {e_{ABC}}{m_{ABG}} + {e_{ACD}}{m_{ADG}}$

(10)

新三角形△ABD的密度ρ_ABD为

${\rho _{ABD}} = {m_{ABD}}/{s_{ABD}}$

(11)

用同样的方法可以得到新三角形△BCD的密度和内能。得到了新三角形的质量，根据(2)式则可以求得格点的新质量，为了计算格点新速度，需要计算出在质量重映过程中动量的改变量。即先假设每个格点的速度不发生变化，则每个格点的质量变化所带来的动量改变量ΔM为

$\Delta \mathit{\boldsymbol{M}} = \sum\limits_{i = 2}^4 {{\mathit{\boldsymbol{v}}_i} \cdot \Delta {m_i}}$

(12)

式中:v_i是该格点的速度，Δm_i是该格点的质量变化值。然后利用动量改变量求出所有格点的速度修正值

${\mathit{\boldsymbol{v}}^{{\rm{Revised}}}} = \frac{{\Delta \mathit{\boldsymbol{M}}}}{{\sum\limits_{i = 1}^n {m_i^{{\rm{New}}}} }}$

(13)

则格点的新速度为

$\mathit{\boldsymbol{v}}_i^{{\rm{New}}} = \mathit{\boldsymbol{v}}_i^{{\rm{Old}}} + {\mathit{\boldsymbol{v}}^{{\rm{Revised}}}}\;\;\;\left( {i = 1,2,3,4} \right)$

(14)

同样地，根据能量守恒求三角形的内能修正量，首先计算出动能的改变量ΔE^k

$\Delta {E^{\rm{k}}} = \frac{1}{2}\sum\limits_{i = 1}^4 {m_i^{{\rm{New}}}{{\left( {\mathit{\boldsymbol{v}}_i^{{\rm{New}}}} \right)}^2}} - \frac{1}{2}\sum\limits_{i = 1}^4 {m_i^{{\rm{Old}}}{{\left( {\mathit{\boldsymbol{v}}_i^{{\rm{Old}}}} \right)}^2}}$

(15)

式中：m_i^New和m_i^Old分别为新、旧格点的质量，v_i^New与v_i^Old分别为新、旧格点的速度。根据能量守恒，动能改变量的大小等于内能改变量的大小，然后将内能改变量均分给所涉及的三角形

${E^{{\rm{Revised}}}} = - \Delta {E^{\rm{k}}}/2$

(16)

则新三角形修正后的内能为

${{E'}_i} = {E_i} + {E^{{\rm{Revised}}}}\;\;\left( {i = 1,2} \right)$

(17)

(2) △ABC与△ACD包含不同的介质，如图 3(b)所示。

在这种情况下，为了保证不出现“混合网格”，将△ACD中的物质转移到与之相邻且包含同种材料的三角形中去。假设△ADE和△CDF含有与△ACD相同的材料，则

${{m'}_{ADE}} = {m_{ADE}} + {m_{ACD}}/2,\;\;\;\;\;\;{{m'}_{CDF}} = {m_{CDF}} + {m_{ACD}}/2$

(18)

${{E'}_{ADE}} = {E_{ADE}} + {E_{ACD}}/2,\;\;\;\;\;\;{{E'}_{CDF}} = {E_{CDF}} + {E_{ACD}}/2$

(19)

式中：m^′和E^′表示新质量和新内能。

△ABD和△BCD则分别继承△ABG和△BCG的物质，即

${m_{ABD}} = {m_{ABG}},\;\;\;\;\;{m_{BCD}} = {m_{BCG}}$

(20)

${E_{ABD}} = {E_{ABG}},\;\;\;\;\;{E_{BCD}} = {E_{BCG}}$

(21)

得到新三角形的质量后，根据(2)式计算出格点的新质量。然后采用与上述类似的方法对速度和内能进行修正。特别提到，△ADE和△CDF的物质不总是与△ACD的物质相同，有可能只有一个三角形与△ACD的物质相同，在这种情况下，将△ACD所有的物质转移到该三角形中；而如果没有三角形与△ACD的物质相同，则△ACD中的物质将被丢弃。

3.2 “长边劈裂”后物理量守恒重映

如图 4所示，当AC边大于给定值的时候，则进行“长边劈裂”操作，由长边AC变为两个短边AE和CE。本研究中假设格点E位于AC边的中点处，且其速度为格点A和格点C速度之和的一半，即

$\left\{ \begin{array}{l} {\mathit{\boldsymbol{x}}_E} = \left( {{\mathit{\boldsymbol{x}}_A} + {\mathit{\boldsymbol{x}}_C}} \right)/2\\ {\mathit{\boldsymbol{v}}_E} = \left( {{\mathit{\boldsymbol{v}}_A} + {\mathit{\boldsymbol{v}}_C}} \right)/2 \end{array} \right.$

(22)

式中：x表示位置，v表示速度。三角形的数目由原来的2个变为4个。接下来进行物理量重映。

首先进行质量和内能重映。假设△ABC和△ACD的质量分别为m_ABC和m_ACD，内能为E_ABC和E_ACD，则新三角形△ABE、△BCE、△CDE、△ADE的质量和内能分别为

${m_{ABE}} = \frac{1}{2}{m_{ABC}},\;\;\;\;{m_{BCE}} = \frac{1}{2}{m_{ABC}}, \\ {m_{CDE}} = \frac{1}{2}{m_{ACD}},\;\;\;\;{m_{ADE}} = \frac{1}{2}{m_{ACD}}$

(23)

${E_{ABE}} = \frac{1}{2}{E_{ABC}},\;\;\;\;{E_{BCE}} = \frac{1}{2}{E_{ABC}}, \\ {E_{CDE}} = \frac{1}{2}{E_{ACD}},\;\;\;\;{E_{ADE}} = \frac{1}{2}{E_{ACD}}$

(24)

新三角形的密度为

${\rho _{ABE}} = {\rho _{BCE}} = {\rho _{ABC}},\;\;\;\;{\rho _{CDE}} = {\rho _{ADE}} = {\rho _{ACD}}$

(25)

求得新三角形的质量后，根据(2)式可以求得格点的新质量。同样地，根据动量守恒修正格点速度，根据能量守恒修正新三角形的内能。

3.3 “短边融合”后的物理量守恒重映

在图 5中，执行完“短边融合”操作后，原格点F、△ABF和△AEF将被删除。原来与F点连接的格点重新与格点A相连组成新的三角形。在执行“短边融合”的操作时，新旧三角形的数目会发生变化，因此在计算新三角形的物理量时需要提前计算重叠系数C_jⁱ。如图 8所示，C_jⁱ指新三角形△_i中包含旧三角形△_j的部分占旧三角形△_j的比例，此处比例指的是面积比例。例如：C₁¹表示新三角形△ABC中包含旧三角形△ABF的部分(红色区域)占△ABF面积的比例。由于本研究假设在网格单元中物理量(质量、内能)密度是常数，因此网格单元的面积之比等于网格内的物理量(质量、内能)大小之比。已知重叠系数可以很方便地计算新三角形的质量和内能。

图 8 计算重叠系数

Figure 8. Calculation of overlapping coefficients

下载: 全尺寸图片幻灯片

接下来进行物理量重映，分为两种情况。

(1) “短边融合”发生在同种介质中，如图 5(a)所示。

首先对质量和内能进行重映。假设旧三角形的质量和内能为m_j、E_j(j=1，2，…，k₁；k₁为旧三角形的个数)，新三角形的质量、面积和内能分别为mⁱ、sⁱ和Eⁱ(i=1, 2, …k₂；k₂为新三角形的个数)。根据重叠系数的定义，则新三角形的质量和内能分别为

${m^i} = \sum\limits_{j = 1}^{{k_1}} {C_j^i{m_j}}$

(26)

${E^i} = \sum\limits_{j = 1}^{{k_1}} {C_j^i{E_j}}$

(27)

新三角形的密度为

${\rho ^i} = {m^i}/{s^i}$

(28)

求得新三角形的质量后，根据(2)式可以求得格点的新质量。同样地，根据动量守恒修正格点速度，根据能量守恒修正新三角形的内能。

(2) “短边融合”发生在不同介质中，如图 9所示。

图 9 多物质界面处短边融合后物理量的重映过程

Figure 9. Remapping after edge-merging on the interface

下载: 全尺寸图片幻灯片

当“短边融合”发生在物质界面处时，为了保证不出现“混合网格”，采用如下策略进行物理量重映。

① 确定新三角形的物质类型, 本研究中规定新三角形的物质类型与其对应的旧三角形物质类型一致。例如，新三角形△ABC应与旧三角形△BCF的物质类型一致，新三角形△ACD应与旧三角形△CDF的物质类型一致。

② 对质量和内能进行初步的重映

${m^i} = \sum\limits_{\begin{array}{*{20}{c}} {j = 1}\\ {{T_j} = {T_i}} \end{array}}^{{k_1}} {C_j^i{m_j}}$

(29)

${E^i} = \sum\limits_{\begin{array}{*{20}{c}} {j = 1}\\ {{T_j} = {T_i}} \end{array}}^{{k_1}} {C_j^i{E_j}}$

(30)

式中：T表示物质类型。(29)式和(30)式对j的取值有限制，即：第j个旧三角形应与第i个新三角形的物质类型一致。例如在求新三角形△ADC的质量和内能时，只有旧三角形△ABF、△BCF、△CDF对△ADC有贡献，分别为图 9(b)中的△₃、△₄和△₆，而△₅和△₇则对新三角形△ADC没有贡献。

③ 为了保证物质守恒，在步骤②中，对(29)式和(30)式没有贡献的C_jⁱm_j和C_jⁱE_j将被转移到与第j个旧三角形对应的新三角形中。例如图 9(b)中的重叠部分△₅和△₇：△₇将被转移到与旧三角形△DEF对应的新三角形△ADE中；而△₅因其所在的旧三角形会被删除，所以并没有对应的新三角形让其转移，本研究采取的策略是将其均分到周边包含相同介质的三角形中，即新三角形△ADE中。如果该三角形周边没有包含相同介质的三角形，则其物质将被丢弃。

④ 上述步骤完成后，则完成了质量和内能在多物质界面处的重新分配，根据(2)式得到格点的新质量。类似地，根据动量守恒修正格点速度，根据能量守恒修正新三角形的内能。从图 9和上述处理过程可以看到，对多介质界面处的网格进行“短边融合”相当于对边界做了一个平滑处理，原边界从AFD变为AD，如果AF长度足够小，则这种处理带来的人工误差就越小。

3.4 “帽子戏法”操作后的物理量守恒重映

如图 6所示，“帽子戏法”实际上可以看作是“短边融合”和“长边劈裂”两种操作组合而成的，先对长边AD进行劈裂，再融合掉短边FG。本研究中假设G点是AD边的中点，且速度是A点与D点速度之和的一半，即

$\left\{ \begin{array}{l} {\mathit{\boldsymbol{x}}_G} = \left( {{\mathit{\boldsymbol{x}}_A} + {\mathit{\boldsymbol{x}}_D}} \right)/2\\ {\mathit{\boldsymbol{v}}_G} = \left( {{\mathit{\boldsymbol{v}}_A} + {\mathit{\boldsymbol{v}}_D}} \right)/2 \end{array} \right.$

(31)

因此，该方法的重映过程分为两个部分，即上半部分的“短边融合”后的重映和下半部分“长边劈裂”后的重映。因为“短边融合”发生在同种介质中，所以与3.3节中单一介质内的“短边融合”类似，需要先计算重叠系数，然后根据重叠系数计算新三角形的质量和内能。下半部分的重映过程较为简单，因为G点是AD边的中点，只需将△ADE的质量及内能均分到△AEG和△DEG中。得到三角形的质量和内能后，根据(2)式得到格点的新质量，然后利用动量守恒对格点的速度进行修正，利用能量守恒对三角形单元的内能进行修正。

4. 测试算例

4.1 冲击波与气泡相互作用

该算例来自文献[4, 11, 13]，几何尺寸如图 10所示，中间区域填充氦气，其余部分填充空气。区域右边界是可移动活塞，以速度u^*=-124.824 m/s向左运动。氦气和空气都采用多方气体状态方程描述，初始条件分别为：(ρ，p，e，λ)_Helium=(0.182, 10⁵, 8.479 2×10⁵, 1.648)，(ρ，p，e，λ)_Air=(1, 10⁵, 2.5×10⁵, 1.4)，其中ρ、p、e分别为密度、压力和比内能，λ为多方气体指数，数值均采用kg-m-s单位制。模拟中网格数量为350×89，黏性系数取为0.000 1。

图 10 激波与气泡相互作用

Figure 10. Shock/bubble interaction

下载: 全尺寸图片幻灯片

图 11是不同时刻t激波与气泡相互作用的模拟结果，其中左列是文献[11]中给出的压力云图与气泡边界线的模拟结果，右列是利用本研究方法模拟得到的结果。可见，二者符合较好，验证了本研究方法的准确性。

图 11 本研究模拟结果与文献[11]的对比

Figure 11. Comparison of simulation results between this study and Ref.[11]

下载: 全尺寸图片幻灯片

图 12给出了从激波开始作用于气泡到气泡即将离开原来位置的时间段内，本研究模拟结果与实验结果^[13]的对比，二者基本一致。图 13给出了本研究、文献[11]模拟结果和实验结果^[13]在激波作用后期(t=1 342.153×10^-6时刻)的对比。可以看到，本研究结果与文献[11]结果一致，但与实验有一些差异，其原因可能是实验条件与计算模拟中假设的条件有一定差异，对此将在后续进行更仔细的研究。

图 12 本研究模拟结果与实验^[13]的对比

Figure 12. Comparison of simulations results of this study with experiment results^[13]

下载: 全尺寸图片幻灯片

图 13 本研究模拟、文献[11]模拟以及实验^[13]结果在t=1 342.153×10^-6时刻的对比

Figure 13. Comparison of simulation result of this study, simulation result of Ref.[11], and experiment result^[13] at t=1 342.153×10^-6

下载: 全尺寸图片幻灯片

4.2 R-T不稳定性

该算例来自文献[4]。如图 14所示，模型由两种不同密度的流体组成：上层为重流体，密度为2；下层为轻流体，密度为1。两种流体都采用理想气体状态方程描述，多方指数λ=1.4。该模型区域大小为：Ω=[0, 1/3]×[0, 1], 界面处的曲线方程为：y(x)=0.5+0.01cos(6πx)。初始时各点都处于静止状态，所受到的压力为

图 14 单模R-T不稳定性

Figure 14. Single mode R-T instability

下载: 全尺寸图片幻灯片

${p_{\rm{d}}}\left( {x,y} \right) = 1 + \frac{1}{2}{\rho _{\rm{u}}}g + {\rho _{\rm{d}}}g\left( {\frac{1}{2} - y} \right)$

(32)

${p_{\rm{u}}}\left( {x,y} \right) = 1 + {\rho _{\rm{u}}}g\left( {1 - y} \right)$

(33)

式中：下标“u”“d”分别代表重流体、轻流体；g为重力加速度，此处取为0.1。

对于R-T不稳定性的模拟，可以采用单介质和多介质两种方式。在单介质模拟中，上、下两层液体被当作一种介质的两种状态；在多介质模拟中，上、下两层液体被当作不同介质。两者的区别会带来网格动态重分的区别，即：在多介质模拟中，由于物质界面的存在，界面处的网格重分将采用前述的专门处理界面处网格的方法。R-T不稳定发展过程中，随着“尖钉”和“气泡”的形成及卷曲，网格会出现极大的变形。该算例可以有效地测试本研究提出的网格局域动态重分方法改善单一介质内部及多介质界面处大变形网格的能力。

4.2.1 单介质的R-T不稳定性过程

首先针对单介质的R-T不稳定性进行模拟。采用3种不同的网格划分：24×72、48×144和100×300，边界条件为周期边界，黏性系数为0.000 1。图 15~图 17是不同网格对应的密度云图。

图 15 不同时刻24×72网格对应的密度云图(单介质)

Figure 15. Density maps at different time on 24×72 grid (Single material)

下载: 全尺寸图片幻灯片

图 16 不同时刻48×144网格对应的密度云图(单介质)

Figure 16. Density maps at different time on 48×144 grid (Single material)

下载: 全尺寸图片幻灯片

图 17 不同时刻100×300网格对应的密度云图(单介质)

Figure 17. Density maps at different time on 100×300 grid (Single material)

下载: 全尺寸图片幻灯片

根据3种网格模拟的单介质R-T不稳定性的演化过程，可以得出以下结论。

(1) 本研究提出的动态局域重分方法可以有效地处理单介质内的大变形网格，成功地模拟了R-T不稳定从线性阶段到非线性阶段再到湍流混合阶段的演化过程。这在没有采用网格动态局域重分的拉氏方法中很难模拟出来。

(2) 随着网格数目的增加，该方法可以更精细地刻画流体的运动过程，图 18是t=6时刻，网格分别为24×72、48×144和100×300的模拟结果。从图 18中可以看到，密网格可以较为精细地捕捉到由流体的切向速度差异所引起的K-H(Kelvin-Helmholtz)不稳定性。

图 18 t=6时刻不同网格对应的流体结构

Figure 18. Flow structures with 3 mesh grids at t=6

下载: 全尺寸图片幻灯片

图 19(a)~图 19(c)给出了单介质模拟在尖钉附近的网格形貌，此时的网格相比于初始网格，其拓扑结构和单元数目都已经发生了改变。图 20给出了初始网格为100×300时，不同时刻网格重分操作的发生次数。图 21是该模拟过程中相应的质量、动量及能量变化量(Δm/m、M_x/M_y及ΔE/E)，可以看到质量没有发生变化，x方向动量相对于y方向动量的变化极小，并且能量的最大损失量只占总能量的0.001 1%，可以忽略不计。

图 19 t=6时刻不同计算方式对应的网格

Figure 19. Mesh of different simulations at t=6

下载: 全尺寸图片幻灯片

图 20 单介质模拟每个操作发生的次数

Figure 20. Times of each operations in single material simulation

下载: 全尺寸图片幻灯片

图 21 单介质模拟质量、动量及能量的改变量

Figure 21. Variation of mass, momentum and energy in single material simulation

下载: 全尺寸图片幻灯片

4.2.2 多介质的R-T不稳定过程

从多介质的角度模拟R-T不稳定性的演化过程，以验证本研究提出的方法对界面处大变形网格处理的有效性。计算网格为100×300，边界条件为周期边界，黏性系数为0.000 1，模拟结果如图 22所示。

图 22 不同时刻100×300网格的密度云图(多介质)

Figure 22. Material maps at different time on 100×300 grid (Multi-material)

下载: 全尺寸图片幻灯片

根据多介质R-T不稳定性的模拟结果，可以得出如下结论。

(1) 网格动态局域重分方法较成功地应用到多介质界面处，改善了多介质界面处的大变形网格。这为处理多介质界面处的大变形网格提供了一个新思路。

(2) 图 23是t=6时多介质R-T不稳定性与单介质R-T不稳定性在相同的尺寸下模拟结果的对比，总体形态一致，在界面处略有不同。一方面是由于在多介质模拟中，界面采用了更严格的网格重分判据，使多介质界面处的网格重分与单介质内部相应位置的网格重分次数不同；另一方面是多介质处的重映方式不同。最终导致界面扰动的发展略有不同。

图 23 相同网格下多介质(左)与单介质(右)R-T不稳定性对比

Figure 23. Comparison between multi-material (Left) and single material (Right) R-T instability on the same mesh

下载: 全尺寸图片幻灯片

图 19(d)给出了多介质100×300网格模拟过程中尖钉附近的网格形貌。图 24是该模拟过程中，不同时刻网格重分操作的发生次数。图 25是相应的质量、动量及能量变化，可以看出：在t=8之前，能量误差与质量误差较小；在t=8之后，误差开始明显增加。对应图 22的发展过程可以看到，t=8之后，湍流开始充分发展，此时多物质界面处网格重分次数增多，物质丢失的可能性增大，因此质量和能量误差增加，这是由本研究采取的算法决定的。但是在t=8之前，湍流尚未形成或者尚未充分发展，质量和能量误差最高达到0.1%，该误差仍然很小。因此本研究算法对于湍流尚未充分发展阶段的模拟仍具有指导意义，但是当湍流充分发展后，数值误差较大，不再具有明显的指导意义。

图 24 多介质模拟每个操作发生的次数

Figure 24. Times of each operations in multi-material simulation

下载: 全尺寸图片幻灯片

图 25 多介质模拟质量、动量和能量的改变量

Figure 25. Variation of mass, momentum and energy in multi-material simulation

下载: 全尺寸图片幻灯片

5. 结论

本研究提出的网格动态局域重分方法不仅可以处理单一介质内部的大变形网格，同时也为处理多介质界面处的大变形网格提供了一个有效的策略。该方法由4种基本网格操作组成，即“对角线交换”“长边劈裂”“短边融合”和“帽子戏法”，每种方法对应不同的发生条件和物理量重映策略。利用该方法对冲击波与气泡相互作用和R-T不稳定性的演化过程进行了数值模拟，取得了良好的模拟结果，其中R-T不稳定性模拟分单一介质和多介质两种情况。通过模拟结果可以看出，该重分方法能够有效地处理网格畸变问题，完整模拟出流体运动的过程。

图 1 交错网格型拉格朗日方法^[10]

Figure 1. Staggered grid hydrodynamics Lagrangian discretization

下载: 全尺寸图片幻灯片

图 2 周期性网格

Figure 2. Periodical mesh

下载: 全尺寸图片幻灯片

图 3 对角线交换操作

Figure 3. Diagonal-swapping

下载: 全尺寸图片幻灯片

图 4 长边劈裂操作

Figure 4. Edge-splitting

下载: 全尺寸图片幻灯片

图 5 短边融合操作

Figure 5. Edge-merging

下载: 全尺寸图片幻灯片

图 6 帽子戏法操作

Figure 6. Hat-trick on the interface

下载: 全尺寸图片幻灯片

图 7 网格动态局域重分流程图

Figure 7. Main flow of the remeshing algorithm

下载: 全尺寸图片幻灯片

图 8 计算重叠系数

Figure 8. Calculation of overlapping coefficients

下载: 全尺寸图片幻灯片

图 9 多物质界面处短边融合后物理量的重映过程

Figure 9. Remapping after edge-merging on the interface

下载: 全尺寸图片幻灯片

图 10 激波与气泡相互作用

Figure 10. Shock/bubble interaction

下载: 全尺寸图片幻灯片

图 11 本研究模拟结果与文献[11]的对比

Figure 11. Comparison of simulation results between this study and Ref.[11]

下载: 全尺寸图片幻灯片

图 12 本研究模拟结果与实验^[13]的对比

Figure 12. Comparison of simulations results of this study with experiment results^[13]

下载: 全尺寸图片幻灯片

图 13 本研究模拟、文献[11]模拟以及实验^[13]结果在t=1 342.153×10^-6时刻的对比

Figure 13. Comparison of simulation result of this study, simulation result of Ref.[11], and experiment result^[13] at t=1 342.153×10^-6

下载: 全尺寸图片幻灯片

图 14 单模R-T不稳定性

Figure 14. Single mode R-T instability

下载: 全尺寸图片幻灯片

图 15 不同时刻24×72网格对应的密度云图(单介质)

Figure 15. Density maps at different time on 24×72 grid (Single material)

下载: 全尺寸图片幻灯片

图 16 不同时刻48×144网格对应的密度云图(单介质)

Figure 16. Density maps at different time on 48×144 grid (Single material)

下载: 全尺寸图片幻灯片

图 17 不同时刻100×300网格对应的密度云图(单介质)

Figure 17. Density maps at different time on 100×300 grid (Single material)

下载: 全尺寸图片幻灯片

图 18 t=6时刻不同网格对应的流体结构

Figure 18. Flow structures with 3 mesh grids at t=6

下载: 全尺寸图片幻灯片

图 19 t=6时刻不同计算方式对应的网格

Figure 19. Mesh of different simulations at t=6

下载: 全尺寸图片幻灯片

图 20 单介质模拟每个操作发生的次数

Figure 20. Times of each operations in single material simulation

下载: 全尺寸图片幻灯片

图 21 单介质模拟质量、动量及能量的改变量

Figure 21. Variation of mass, momentum and energy in single material simulation

下载: 全尺寸图片幻灯片

图 22 不同时刻100×300网格的密度云图(多介质)

Figure 22. Material maps at different time on 100×300 grid (Multi-material)

下载: 全尺寸图片幻灯片

图 23 相同网格下多介质(左)与单介质(右)R-T不稳定性对比

Figure 23. Comparison between multi-material (Left) and single material (Right) R-T instability on the same mesh

下载: 全尺寸图片幻灯片

图 24 多介质模拟每个操作发生的次数

Figure 24. Times of each operations in multi-material simulation

下载: 全尺寸图片幻灯片

图 25 多介质模拟质量、动量和能量的改变量

Figure 25. Variation of mass, momentum and energy in multi-material simulation

下载: 全尺寸图片幻灯片

参考文献(13)

[1]	HIRT C W, AMSDEN A A, COOK J L.An arbitrary Lagrangian-Eulerian computing method for all flow speeds[J].Journal of Computational Physics, 1974, 14(3):227-253. doi: 10.1016/0021-9991(74)90051-5
[2]	ANBARLOOEI H R, MAZAHERI K.Moment of fluid interface reconstruction method in multi-material arbitrary Lagrangian Eulerian (MMALE) algorithms[J].Computer Methods in Applied Mechanics and Engineering, 2009, 198(47):3782-3794. https://www.deepdyve.com/lp/wiley/moment-of-fluid-interface-reconstruction-method-in-axisymmetric-Aki0sZlem3
[3]	ZENG Q. MMALE numerical simulation for multi-material large deformation fluid flows[C]//Journal of Physics: Conference Series, 2014, 510: 012047.
[4]	LOUBRE R, MAIRE P H, SHASHKOV M J, et al.A reconnection-based arbitrary-Lagrangian-Eulerian method[J].Journal of Computational Physics, 2010, 229(12):4724-4761. doi: 10.1016/j.jcp.2010.03.011
[5]	田宝林, 刘妍, 申卫东, 等.可压缩多介质流动的整体ALE方法(GALE)[J].北京理工大学学报, 2010, 30(5):622-625. http://www.oalib.com/paper/4590079 TIAN B L, LIU Y, SHEN W D, et al.Global arbitrary Lagrangian-Eulerian method for compressible multimaterial flows[J].Transactions of Beijing Institute of Technology, 2010, 30(5):622-625. http://www.oalib.com/paper/4590079
[6]	BENSON D J.Volume of fluid interface reconstruction methods for multi-material problems[J].Applied Mechanics Reviews, 2002, 55(2):151-165. doi: 10.1115/1.1448524
[7]	DYADECHKO V, SHASHKOV M.Reconstruction of multi-material interfaces from moment data[J].Journal of Computational Physics, 2008, 227(11):5361-5384. doi: 10.1016/j.jcp.2007.12.029
[8]	LIU J.A second-order changing-connectivity ALE scheme and its application to FSI with large convection of fluids and near contact of structures[J].Journal of Computational Physics, 2016, 304:380-423. doi: 10.1016/j.jcp.2015.10.015
[9]	WICKE M, RITCHIE D, KLINGNER B M, et al.Dynamic local remeshing for elastoplastic simulation[J].ACM Transactions on Graphics, 2010, 29(4):1-11. http://graphics.berkeley.edu/papers/Wicke-DLR-2010-07/
[10]	王瑞利, 林忠, 魏兰.针对流体大变形问题网格邻域可变技术[J].计算物理, 2011, 28(4):501-506. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=jswl201104005 WANG R L, LIN Z, WEI L.Reconnection-based Lagrangian-local remeshing method for large deformations[J].Chinese Journal of Computational Physics, 2011, 28(4):501-506. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=jswl201104005
[11]	STÉPHANE D P.Metric-based mesh adaptation for 2D Lagrangian compressible flows[J].Journal of Computational Physics, 2011, 230(5):1793-1821. doi: 10.1016/j.jcp.2010.11.030
[12]	BARLOW A J, MAIRE P H, RIDER W J, et al.Arbitrary Lagrangian-Eulerian methods for modeling high-speed compressible multimaterial flows[J].Journal of Computational Physics, 2016, 322:603-665. doi: 10.1016/j.jcp.2016.07.001
[13]	HAAS J F, STURTEVANT B.Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[J].Journal of Fluid Mechanics, 1987, 181:41-76. doi: 10.1017/S0022112087002003