ABSTRACT

In this work, the semi-Lagrangian method is extended to the simulation of incompressible multi-medium flow on MAC staggered grid method. We track the fluid interface by solving the level set equation, and the NGFM is used to define the fluid interface condition such that the calculation can be carried out the same as that in the single medium flow. Several classic examples are simulated to test the algorithm given in this work. The numerical results show that the given algorithm can track the interface accurately and capture the physical details of the multi-medium flow.

Keywords:Incompressible Flow, Semi-Lagrangian Method, Marker and Cell Method, Multi-Medium Flow, NGFM Method

不可压缩多介质流动问题的数值模拟

丁灏^*，王春武

南京航空航天大学理学院，江苏 南京

收稿日期：2018年2月7日；录用日期：2018年2月20日；发布日期：2018年2月27日

摘 要

本文将基于MAC (Marker and cell)交错网格上的半拉格朗日方法推广应用于求解不可压缩多介质流动问题。通过求解Level set方程捕捉界面的位置，利用NGFM (New Ghost Fluid Method)方法定义界面边界条件，将多介质问题转化为单介质问题进行计算。给出了一种求解不可压缩多介质流动问题的简单且高效的数值算法。本文针对多个经典算例进行数值试验，数值结果表明本文所给出的方法可以有效的模拟流体的运动细节。

关键词 :不可压缩流体，半拉格朗日方法，MAC方法，多介质流动，NGFM方法

Copyright © 2018 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

1965年，Harlow和Welch提出了基于交错网格的MAC (Marker and Cell)方法 [1] ，将压力定义在网格中心，速度定义在网格边界中心，该方法计算简单，用于N-S方程的数值模拟，是最早用来计算N-S方程的方法之一。半拉格朗日方法是一种将拉格朗日方法和欧拉方法相结合的计算方法，系统的半拉格朗日方法最早是由A. Staniforth和J. Cote于1990年提出的 [2] ，该方法最大的优点是简单且无条件稳定，即数值计算时，可以在较大的范围内选择时间步长。对于天气预报等需要长时间数值模拟的问题，可以大大提高计算效率，所以该方法一出现其应用就在大气科学领域获得了巨大的成功。

对于多介质流动来说，由于流场中存在着两种不同的介质，所以在界面两端密度和粘性等流体属性会发生大梯度的变化，如果采用单介质流动的数值模拟方法，在界面处会造成数值不稳定，甚至使得计算无法进行，因此合理有效的界面处理方法是解决多介质问题的关键。1999年Fedkiw等人提出了虚拟流体方法(Ghost Fluid Method, GFM) [3] ，通过在界面处定义界面边界条件，将多介质问题转化为单介质问题求解，该方法简单且易于推广到高维问题。为了准确模拟气液、气固界面问题，Fedkiw等人又提出了新的虚拟流体方法(new Ghost Fluid Method, GFM) [4] ，充分考虑了刚性较强流体的特性，有效地减少了因界面定义引入的误差，提高了数值模拟的准确性。

本文基于MAC交错网格利用半拉格朗日离散对流方程，对不可压缩水气流动问题进行数值模拟，考虑到水气界面处初始密度、压力或速度不连续，状态方程不同，不能直接进行计算，本文利用NGFM方法来定义界面边界，通过求解Level set方程 [5] [6] 跟踪界面的位置，对溃坝、水滴下落等问题进行了模拟。

2. 流体控制方程

考虑不可压缩Navier-Stokes (N-S)方程：

$u_t+u\cdot \nabla u+\frac{\nabla p}{\rho }=g+\gamma \nabla \cdot \nabla u$ (1)

$\nabla \cdot u=0$ (2)

其中 $u$ 表示速度， $p$ 为压力， $\rho$ 为密度， $g$ 为外力， $\gamma$ 为粘性系数。

为了便于求解N-S方程，本文采用Chorin’s投影法将复杂的N-S方程分裂成4个方程，分别进行离散计算：

$\frac{\partial u}{\partial t}+u\cdot \nabla u=0$ (3)

$\frac{\partial u}{\partial t}=\gamma \nabla \cdot \nabla u$ (4)

$\frac{\partial u}{\partial t}=g$ (5)

$\frac{\partial u}{\partial t}+\frac{1}{\rho }\nabla p=0$ (6)

3. 数值计算方法

本文基于MAC网格进行离散(图1)，即速度 $u$ 定义在在网格边界上，压力 $p$ 定义在网格中心。利用Chorin’s投影方法将方程(1)分裂成对流项、粘性项、外力项、压力项，然后分别进行离散计算。

3.1. 外力项和粘性项的计算

对于外力项方程(5)，使用欧拉方法进行时间离散，计算可以得到中间速度 ${{\displaystyle u}}^A$ 。

对于粘性项，采用中心差分对 $\gamma \nabla u$ 进行离散，使用欧拉方法对时间项进行离散，可以得到:

$u_{i+1/2,j}^v=u_{i+1/2,j}^A+\Delta t\gamma \left(\frac{u_{i+3/2,j}^A+u_{i+1/2,j+1}^A+u_{i-1/2,j}^A+u_{i+1/2,j-1}^A-4u_{i+1/2,j}^A}{\Delta x^2}\right)$ (7)

$v_{i,j+1/2}^v=v_{i,j+1/2}^A+\Delta t\gamma \left(\frac{v_{i,j+3/2}^A+v_{i+1,j+1/2}^A+v_{i,j-1/2}^A+v_{i-1,j+1/2}^A-4v_{i,j+1/2}^A}{\Delta y^2}\right)$ (8)

通过式(7)和(8)可以得到中间流体速度 ${{\displaystyle u}}^v$ 。

图1. 交错网格

3.2. 对流项计算

对流方程方程组(3)的标量一般形式为

$\frac{\partial q}{\partial t}+u\cdot \nabla q=0$ (9)

其中状态量q可以用来表示速度、密度或温度。利用根据拉格朗日方法的思想，把每一个网格点看作一个流体粒子，流体粒子会随着流体运动，跟踪并计算每一时刻每个流体粒子的状态，得到方程的解。下面首先来确定上一时刻流体粒子所在的位置。

流体粒子运动满足方程

$\frac{\text{d}x}{\text{d}t}=u\left(x\right)$ (10)

方程(10)是对流方程(3)的特征方程，利用欧拉公式离散可以得到

$x^{n+1}=x^n+u\left(x^{n+1}\right)\Delta t$ (11)

$x^n=x^{n+1}-u\left(x^{n+1}\right)\Delta t$ (12)

其中 $u\left(x^{n+1}\right)$ 为 $t^n$ 时刻处的速度。

利用流体的速度，通过求解方程(12)，可以求出流体粒子在时刻的位置 (如图2所示)，格式的CFL条件保证了流体粒子移动的距离不会超过一个网格间距。根据时刻的位置处流体粒子以及周围粒子的状态(对于速度，这里使用求解粘性项得到的中间速度)，可以插值得到流体粒子在处的状态，这样可以得到中间速度。

3.3. 压力项计算

对压力项进行离散时，先对在时间方向进行离散：

(13)

等式两边同时求散度，

(14)

结合，可以得到：

图2. 粒子轨迹图

(15)

对于梯度，采用欧拉差分进行离散，而对于，则用中心差分进行离散得到：

(16)

式(16)是一个关于压力的线性方程组，是一个对称半正定矩阵，用预处理共轭梯度法 [7] 求解线性方程得到压力，然后再通过(17)-(18)可以计算得到下一时间层的速度。

(17)

(18)

3.4. 时间步长

在时间步长上，由于将N-S方程分裂成几项进行离散，时间步长需要分别去考虑，对于对流项，时间步长需满足

(19)

对于粘性项，则有

(20)

取最小的时间步长

(21)

4. 界面处理方法

由于流场中存在着两种不同的介质，所以在界面两边流体参数和状态方程不同，如果采用单介质流动的数值模拟方法，在界面处会造成数值不稳定，甚至使得计算无法进行，本文利用NGFM方法为每一种流体在界面处定义边界条件，将多介质问题转化为单介质问题再利用格式(7)-(18)进行计算，可以得到流体的流动状态。

本文通过求解Level set方程

(22)

来跟踪界面的位置，其中为定义了符号的距离函数，为流体速度。方程(22)的时间和空间离散分别用三阶TVD Runge-Kutta方法和五阶WENO格式 [8] 。

5. 数值实验

本文利用所给出的算法对溃坝和水滴自由下落问题进行了数值模拟，计算区域的边界为固壁边界，采用无滑移边界条件。

算例1 计算区域为一个4 m × 2.2 m的封闭矩形区域，网格数400 220，在计算区域的最左端有一个长度为1 m，高度为2 m的水坝，水坝中的水受到重力，状态为，气体的状态。水随着水坝溃解由静止开始向下自由滩落。图3是t = 0 s，0.9 s，1.4 s，2.0 s时刻液面的位置图。

图4是VOF方法和本文方法计算得到的波前位置，通过与文献 [9] 计算结果比较，可见本文给出的方法可以有效的对溃坝问题进行模拟。

算例2 计算区域为一个3 m × 2 m的封闭矩形区域，网格数300 × 200，在距离右边界1.5 m和底部1 m处有一个水滴，水滴半径为1/3 m，密度，粘性，气体的密度，粘性。水滴受到重力的作用自由下落。

图5为t = 0 s，0.15 s，0.3 s，0.4 s，0.45 s，0.5 s，0.55 s，0.6 s时刻的水滴运动界面图。可见，水滴在受到重力自由下落的过程中，受到空气阻力的作用产生了向两边运动的速度，开始发生微小形变。落地后，水滴与地面发生碰撞，慢慢平摊下来。

图6是水滴平均高度变化图，从图中可以看出水滴的平均高度随着时间抛物下降，到0.38 s时水滴接触地面，水滴下降的趋势开始变缓。图7是水滴垂直方向的平均速度变化图。图8是水滴的体积变化图，在计算过程中有微小的质量丢失情况。通过与文献 [10] 中的计算结果进行对比，模拟效果较好。

图3. 溃坝过程图

图4. VOF和本文计算的波前位置对比图

图5. 水滴下落与地面交互过程图

图6. 水滴的平均高度随时间变化图

图7. 水滴垂直方向速度随时间变化图

图8. 水滴体积随时间变化图

6. 结论

本文利用NGFM方法定义界面边界，充分考虑界面两侧流体的性质，利用基于MAC交错网格上的半拉格朗日方法对不可压缩多介质流动问题进行了数值模拟，数值结果表明算法可以有效地反映流体运动的细节。

文章引用

丁 灏,王春武. 不可压缩多介质流动问题的数值模拟
Numerical Simulation of Incompressible Multi-Medium Flow Problem[J]. 现代物理, 2018, 08(01): 18-26. http://dx.doi.org/10.12677/MP.2018.81003

参考文献 (References)