删除或更新信息，请邮件至freekaoyan#163.com(#换成@)

基于气液相界面捕捉的统一气体动理学格式

本站小编 Free考研考试/2022-01-01

王昭¹, 严红¹^,²^,³^,^*^,

1西北工业大学动力与能源学院, 西安 710072

2西北工业大学陕西省航空发动机内流动力学重点实验室, 西安 710072

3先进航空发动机协同创新中心, 北京 100191;

UNIFIED GAS-KINETIC SCHEME FOR TWO PHASE INTERFACE CAPTURING

WangZhao¹, YanHong¹^,²^,³^,^*^,

1School of Power and Energy, Northwestern Polytechnical University, Xi’an 710072, China;

2Shaanxi Key Laboratory of Internal Aerodynamics in Aero-Engine, Northwestern Polytechnical University, Xi’an 710072, China;

3Collaborative Innovation Center for Advanced Aero-Engine, Beijing 100191, China;

中图分类号:O359
文献标识码:A

通讯作者:^*严红, 教授, 主要研究方向: 超声速流动. E-mail:yanhong@nwpu.edu.cn^*严红, 教授, 主要研究方向: 超声速流动. E-mail:yanhong@nwpu.edu.cn
版权声明:2018力学学报期刊社力学学报期刊社所有
基金资助:

国防基础科研科学挑战计划资助项目(TZ2016001).

展开

摘要
气液相界面运动的研究无论是在科学还是工程领域都是非常重要的. 其中, 非平衡流动的计算尤其受到关注. 基于此, 我们构造了捕捉气液相界面的统一气体动理学格式. 由于统一气体动理学格式将自由输运和粒子碰撞耦合起来更新宏观物理量和微观分布函数, 故而可以求解非平衡流动. 具体思路是, 通过将范德瓦尔斯状态方程所表达的非理想气体效应引入统一气体动理学格式之中来捕捉气液相界面, 两相的分离与共存通过范德瓦尔斯状态方程描述. 由于流体在椭圆区域是不稳定的, 因此气液相界面可以通过蒸发和凝结过程自动捕捉. 如此, 一个锋锐的相界面便可以通过数值耗散和相变而得到. 利用该方法得到麦克斯韦等面积律(Maxwell construction)对应的数值解, 并与其相应的理论解相比较, 二者符合良好. 而后, 通过对范德瓦尔斯状态方程所描述的液滴表面张力进行数值计算, 验证了Laplace定理. 此外, 通过模拟两个液滴的碰撞融合过程, 进一步证明了该格式的有效性. 但是, 由于范德瓦尔斯状态方程的特性, 其所构造的格式仅适用于液/气两相密度比小于5的情况.

关键词：统一气体动理学格式;范德瓦尔斯状态方程;相界面捕捉;麦克斯韦等面积律;液滴表面张力
Abstract
The study of interfacial motion of gas-liquid phase is very important in science and engineering. Considering the non-equilibrium flow calculation, a unified gas-kinetic scheme for gas-liquid two phase interface capturing is presented in this paper. Since the free transport and particle collision are coupled to update the macroscopic variables and microscopic distribution functions, the unified gas-kinetic scheme can solve the non-equilibrium flow. The van der Waals (vdW) equation of state (EOS) is included to describe the coexistence of gas and liquid and the phase transition between them. Because of the characteristics of vdW EOS, the interface between gas and liquid can be captured naturally through condensation and evaporation processes. As a result, the new scheme can solve the gas-liquid two phase problems. Finally, the proposed method is used to obtain the numerical solution of Maxwell construction, which agrees well with the corresponding theoretical solution. Then, the Laplace’s theorem is verified by numerical calculation of the surface tension of the droplet corresponding to the van der Waals state equation. In addition, the collision of the two droplets is simulated, which proves the validity of the scheme further. However, due to the characteristics of the van der Waals equation of state, the constructed scheme is only applicable to the case where the liquid/gas two-phase density ratio is less than 5.

Keywords：unified gas-kinetic scheme;van der Waals equation of state;interface capturing;Maxwell construction;surface tension

-->0
PDF (9455KB)元数据多维度评价相关文章收藏文章
本文引用格式导出EndNote Ris Bibtex 收藏本文-->
王昭, 严红. 基于气液相界面捕捉的统一气体动理学格式[J]. 力学学报, 2018, 50(4): 711-721 https://doi.org/10.6052/0459-1879-17-364
Wang Zhao, Yan Hong. UNIFIED GAS-KINETIC SCHEME FOR TWO PHASE INTERFACE CAPTURING[J]. Acta Mechanica Sinica, 2018, 50(4): 711-721 https://doi.org/10.6052/0459-1879-17-364
气液相界面运动的研究无论是在科学还是工程领域都是非常重要的[1-2]. 为了精确捕捉气液相界面, 一般有两种思路: (1)界面追踪, 该方法需要求解界面对流输运方程, 最具代表性的有流体体积法(VOF)[3-4]和水平集方法[5]. 以上两种方法均是以求解 Navier-Stokes (N-S)方程为基础;(2)自动捕捉, 亦即通过求解气液相变引起的分离, 来捕捉气液相界面. 该方法相对于界面追踪更有优势. 首先不需要求解界面的对流输运方程, 加之它是以波尔兹曼方程为基础, 对物理现象的描述更接近本质, 所以计算结果更加精确.
由波尔兹曼方程离散所得的格子波尔兹曼方程, 在特定条件下(低马赫数), 通过 Chapman-Enskog 展开可以推导出不可压 N-S 方程[6]. 以此为基础的格子波尔兹曼方法由于不再受理想气体的限制, 在近几十年里得到充分的发展和广泛的应用[7-23]. 为了能够自动捕捉相界面, 需要对非理想状态下, 分子之间的作用力加以建模. 基于此, ****们已经建立了多种模型 [7-10]. 其中, 单组分伪势模型由于形式简洁受到广泛应用. Zhang 和Chen [11]在 2003 年对单组分伪势模型引入范德瓦尔斯(van der Waals, vdW)状态方程, 并提出直接采用势函数的方式计算粒子间相互作用力来提高伪势模型的稳定性. 为提高模型模拟的密度比, 2006 年, Yuan等 [12] 引入不同非理想气体状态方程, 发现模型可以在很高密度下(>1000)仍然保持稳定. 2007 年, Kupeshtokh 等 [13]发现, 当采用修正平衡态分布函数作用力方法时, 模拟所得的结果受松弛时间的影响, 这与物理实际不吻合. 为解决这个问题, 他们提出了一种新的作用力模型 [13], 即有限差分模型(EDM). 此模型基本解决了单组份伪势模型在大密度比条件下不稳定的问题, 但是其采用对比态方程中涉及物理单位与实际单位的转换关系不满足格子波尔兹曼模拟的相似性原理, 会使其在应用中出现很多问题. 另外, 模型的提出大多源于经验, 缺乏理论基础. 因此, 针对以上问题, 胡安杰 [14]对这个模型进行理论分析, 并且基于其模型做出了改进. 此外, 2012 年Gong 和Cheng [24]通过将密度和温度分布函数引入格子波尔兹曼方法, 并嵌入了Peng-Robinson 状态方程, 研究了气泡的生长和脱离水平壁面问题. 最近, 该团队又通过修正伪势模型研究了水平面的薄膜沸腾传热问题 [25], 为格子波尔兹曼方法在传热方面的应用做出很大贡献.
如前所述, 格子波尔兹曼方程适用于低速不可压连续流体的模拟, 虽然通过高阶矩的方法试图将格子波尔兹曼方法推广至稀薄流动 [26-29], 但是限于弱非平衡态. 近年来, 基于气体动理学BGK 方程, 离散速度空间的方法发展起来了一套适用于全流域的统一气体动理学格式(unified gas-kinetic scheme, UGKS) [30-32]. 该格式的精要在于将自由输运和粒子碰撞耦合起来更新宏观物理量和微观气体分布函数. 不同于其他诸多方法的关键是, 在计算网格界面通量时采用动理学方程的积分解. 该积分解描述的是两尺度的物理过程, 其一为粒子碰撞所带来的平衡态的积分所体现的流体力学尺度, 其二则是粒子自由输运所体现的动理学尺度. 二者之间的权重取决于网格分辨率和时间步长. 因此, 该格式适用于全流域的计算.
本文通过将范德瓦尔斯状态方程引入到统一气体动理学格式之中, 建立了基于气液相界面捕捉的UGKS. 利用该方法, 得到麦克斯韦等面积律(Maxwell construction)的数值解, 观察其与理论解的符合程度. 而后, 通过对范德瓦尔斯状态方程所描述的液滴表面张力进行数值计算, 验证拉普拉斯定理, 此外, 通过模拟两个液滴的碰撞过程, 进一步证明该格式的有效性.

1 问题描述

对于实际气体而言, 由于考虑了分子间的相互作用力, 故而其状态方程与理想气体状态方程是有本质区别的, 其中, 最早提出关于实际气体状态方程的便是著名的范德瓦尔斯状态方程. 由于在临界温度下, 该方程其所描述的压力?密度曲线是非单调的, 故而可以实现气液两相共存.
通过临界参数加以无量纲化的范德瓦尔斯方程为

p = \frac{8 ρ T}{3 - ρ} - 3 ρ^{2}

此时, 恒温下, 范德瓦尔斯状态方程所描述的曲线如图1所示. 通过麦克斯韦等面积律则可求得两相共存的相参数, 亦即

ρ_{l}

ρ_{g}

p

. 此外,

ρ_{α}

ρ_{β}

表征的是两个极值点处的密度.

显示原图|下载原图ZIP|生成PPT
图1范德瓦尔斯状态方程的曲线示意图.
-->Fig.1The schematic diagram for van der Waals equation of state.
-->

流体密度可按照如下分类[33]: (1)

\frac{1}{ρ} < \frac{1}{ρ_{l}}

时为液相; (2)

\frac{1}{ρ_{l}} < \frac{1}{ρ} < \frac{1}{ρ_{α}}

时为亚稳态域;(3)

\frac{1}{ρ_{α}} < \frac{1}{ρ} < \frac{1}{ρ_{β}}

时为不稳定椭圆域(混合物); (4)

\frac{1}{ρ_{β}} < \frac{1}{ρ} < \frac{1}{ρ_{g}}

时为亚稳态域; (5)

\frac{1}{ρ_{g}} < \frac{1}{ρ}

时为气相.
当流体密度处于椭圆域时, 由于

? p / ? ρ

的斜率为负, 任何小的扰动都会被放大. 假设某个网格位于椭圆域, 其密度有一个正的扰动, 便会导致压力下降, 此时邻近网格的流体便会向该网格汇聚, 如此该网格的密度便会增加, 换言之正的扰动被放大. 故而该区域的流体混合物或蒸发为气相或凝结为液相. 因此, 范德瓦尔斯状态方程有内在的物理机制使得两相实现分离, 并且能够捕捉尖锐的相界面. 为了使用该机制, 发展的界面捕捉格式必须精确预测气液两相的密度. 换言之, 即在没有显含麦克斯韦等面积律的情况下, 格式能够有固有的耗散机制捕捉正确密度跳跃所对应的物理解.
综上所述, 将范德瓦尔斯状态方程所表达的非理想气体效应引入统一气体动理学格式之中, 便可以得到尖锐的相界面以及正确的物理解.

2 基于气液相界面捕捉的统一气体动理学格式

2.1 动理学模型

对于单组分二维气液两相流而言, 其流场控制方程为

f_{t} + u f_{x} + v f_{y} = \frac{g - f}{τ}

这里, 平衡态

g = ρ (\frac{λ}{<})

. 其中,

τ

为松驰时间,

ρ

为密度,

(U, V)

为宏观速度,

K

为内自由度, 此外,

ξ^{2} = {ξ_{1}}^{2} + {ξ_{2}}^{2} + ? + {ξ_{K}}^{2}

. 其中,

ξ_{1}, ξ_{2}, ?, ξ_{K}

表征的是在每个内自由度上的微观速度. 需要强调的是, 此处的

λ

须由非理想气体状态方程, 即式(1) 的范德瓦尔斯状态方程来确定[33]. 其形式如下

λ = \frac{ρ}{2 p} = \frac{3 - ρ}{6 ρ^{2} - 18 ρ + 16 T} = Λ (ρ)

2.2 简化动理学模型

对于二维模型而言, 定义如下的简化气体分布函数

H = \int ? fd ξ

宏观量可如下表示

W = (\begin{array}{l} ρ \\ ρ U \\ ρ V \end{array}) = \int ? (\begin{array}{l} H \\ uH \\ vH \end{array}) dudv

此时的动理学方程变为

H_{t} + u H_{x} + v H_{y} = \frac{H_{eq} - H}{τ}

其中,

H_{eq} = ρ (\frac{λ}{<})

2.3 基于气液相界面捕捉的统一气体动理学格式

统一气体动理学格式是一种有限体积方法, 二维情况下的物理域可划分为控制体

Ω_{i, j}

. 时空离散的分布函数为

H_{i, j}^{n} = H_{i, j, α, β}^{n} = H (t_{n}, x_{i}, y_{j}, u_{α}, v_{β})

此时, 宏观量的表达为

W_{i, j}^{n} = \sum ? (\begin{array}{l} κ_{α, β} H_{i, j}^{n} \\ κ_{α, β} H_{i, j}^{n} u_{α} \\ κ_{α, β} H_{i, j}^{n} v_{β} \end{array})

其中,

κ_{α, β}

为粒子速度在

(u_{α}, v_{β})

处的数值积分权重.
通过守恒性原理, 易得分布函数随时间的演化为

\begin{array}{l} H_{i, j, α, β}^{n + 1} & = H_{i, j, α, β}^{n} + \frac{1}{Δ x} \int_{t^{n}}^{t^{n + 1}} ? (f_{i - 1 / 2, α, β}^{H} - f_{i + 1 / 2, α, β}^{H}) dt + \\ \frac{1}{Δ y} \int_{t^{n}}^{t^{n + 1}} ? (g_{j - 1 / 2, α, β}^{H} - g_{j + 1 / 2, α, β}^{H}) dt + \frac{Δ t}{2} ? \\ (\frac{H_{eq, i, j, α, β}^{n + 1} - H_{i, j, α, β}^{n + 1}}{τ^{n + 1}} + \frac{H_{eq, i, j, α, β}^{n} - H_{i, j, α, β}^{n}}{τ^{n}}) \end{array}

宏观量随时间的演化为

\begin{array}{l} W_{i, j}^{n + 1} = & W_{i, j}^{n} + \frac{1}{Δ x} (F_{i - 1 / 2, j} - F_{i + 1 / 2, j}) + \\ \frac{1}{Δ y} (G_{i, j - 1 / 2} - G_{i, j + 1 / 2}) \end{array}

其中

F_{i - 1 / 2, j} = \int_{t_{n}}^{t_{n + 1}} ? \int ? u ψ H_{i - 1 / 2, j} dudvdt

G_{i, j - 1 / 2} = \int_{t_{n}}^{t_{n + 1}} ? \int ? v ψ H_{i, j - 1 / 2} dudvdt

ψ = (1, u, v)^{T}

f_{i - 1 / 2, α, β}^{H} = \int ? u H_{i - 1 / 2, j} dudv, g_{j - 1 / 2, α, β}^{H} = \int ? v H_{i, j - 1 / 2} dudv

首先求解式(8), 得到n+1步的宏观量, 从而得到对应的平衡态, 以此来更新式(7), 得到n+1 步的分布函数. 故而, 发展UGKS 的核心便在于如何构造界面的气体分布函数.
下面对构造界面气体分布函数加以详细说明. 为了简单起见, 本文采用方向分裂技术来构造格式, 即对于二维流动而言, 认为在y方向上分布函数分布是均匀的, 亦即

H_{y} = 0

. 于是, 在x 方向上, 动理学方程变为

H_{t} + u H_{x} = \frac{H_{eq} - H}{τ}

上式存在局部解析解, 以时间步

t^{n} = 0

的界面

x_{i + 1 / 2} = 0

为例. 式(9)的解随着时间的演化如下

\begin{array}{l} H (0, t, u) = & \frac{1}{τ} \int_{0}^{t} ? H_{eq} (x', t', u') e^{- (t - t') / τ} dt' + \\ e^{- t / τ} H_{0} (x^{n}, 0, u^{n}) \end{array}

其中,

x' = - u (t - t')

.
初始时刻界面附近的分布函数

H_{0}

为

H_{0} (x, 0, u_{α}) = \{\begin{array}{l} H_{i + 1 / 2, α}^{L} + σ_{i, α}^{H} x, x ? 0 \\ H_{i + 1 / 2, α}^{R} + σ_{i + 1, α}^{H} x, x > 0 \end{array}

亦即

H_{i + 1 / 2, α}^{L} = H_{i, α} + (x_{i + 1 / 2} - x_{i}) σ_{i, α}^{H}

其中

H_{i + 1 / 2, α}^{R} = H_{i + 1, α} - (x_{i + 1} - x_{i + 1 / 2}) σ_{i + 1, α}^{H}

\begin{array}{l} H_{eq} (x, t, u_{α}) = & H_{eq}^{0} [1 + (1 - H [x]) a_{H}^{L} x + \\ H [x] a_{H}^{R} x + A_{H} t] \end{array}

而界面附近的平衡态分布函数近似由泰勒展开为

H_{eq}^{0}

其中,

x = 0

为

t = 0

\begin{array}{l} H (0, t, & u_{α}) = M t_{1} H_{eq}^{0} + M t_{2} (α_{H}^{L} H [u_{α}] + α_{H}^{R} (1 - H [u_{α}])) ? \\ u_{α} H_{eq}^{0} + M t_{3} A_{H} H_{eq}^{0} + \\ e^{- t / τ} [(H_{i + 1 / 2, α}^{L} + (- u_{α} t) σ_{i, α}^{H}) H [u_{α}] + \\ (H_{i + 1 / 2, α}^{R} + (- u_{α} t) σ_{i + 1, α}^{H}) (1 - H [u_{α}])] \end{array}

处的Maxwell分布. 通过积分可得

M t_{1} = 1 - e^{- t / τ}

其中,

M t_{2} = τ (- 1 + e^{- t / τ}) + t e^{- t / τ}

M t_{3} = τ (t / τ - 1 + e^{- t / τ})

α_{H}^{L}, α_{H}^{R}, A_{H}

.
上式中各系数

σ_{i}^{H}, σ_{i + 1, α}

和网格分布函数梯度

Δ t

的确定方法具体见文献[34].
综上所述, 对基于气液相界面捕捉的统一气体动理学格式的构造方法作如下总结:
(1) 将原始的动理学方程(2)变为简化分布函数所对应的方程(9).
(2) 网格界面的简化气体分布函数通过动理学方程的积分解(10)而得到.
(3) 通过式(7)和式(8)分别更新分布函数和宏观量.
(4) 对以上三步进行循环, 便可求得流场的时间演化解.

3 算例验证

本节通过求出麦克斯韦等面积律数值解, 并与其对应的理论解相比较, 从而验证格式的有效性. 而后, 通过对范德瓦尔斯状态方程所描述的液滴表面张力的求取, 验证拉普拉斯定理. 最后模拟两个液滴的碰撞过程. 其中, 插值所使用的限制器为van Leer限制器. 时间步长

τ

通过CFL条件来确定. 松弛时间

τ = χ_{1} Δ t + χ_{2} \frac{| p^{l} - p^{r} |}{p^{l} + p^{r}}

定义如下

p^{l} = Λ (ρ^{l}), p^{r} = Λ (ρ^{r})

其中,

χ_{1} = 0.05, χ_{2} = 200

, 上标l和r表示界面的左右侧. 此外, 在所有计算中,

T = 0.9

3.1 麦克斯韦等面积律

气液两相均有稳态、亚稳态与完全非稳态三种状态. 下面采用本文提出的方法首先计算

(X, Y) \in [0,0.05] \times [0,1]

时气液两相密度处于不同初值下的数值解, 观察其在定常情况下是否能够回到麦克斯韦等面积律所描述的理论解. 而后得出不同温度下的气液共存密度, 绘制成气液共存曲线, 并与理论解进行对比. 对于UGKS而言, 计算域为

\times

, 划分为 1

\begin{array}{l} ρ (y) = & ρ_{ini, g} + \frac{ρ_{ini, l} - ρ_{ini, g}}{2} [\tanh \frac{2 (y - 15 L)}{Υ} - \\ \tanh \frac{2 (y - 45 L)}{Υ}] \end{array}

60 个均匀网格. X 与Y 方向均采用周期性边界条件. 初始速度为零, 密度分布为[36]

ρ_{ini, g}, ρ_{ini, l}

其中, 初始的气液密度分布

L

赋为麦克斯韦理论解,

Υ

为Y方向网格边长,

Υ = 2 L

为初始的相界面厚度, 取为

(ρ_{ini, g} = 0.4257, ρ_{ini, l} = 1.6573)

. 为了避免表面张力对密度分布的影响, 初始的相界面是水平的.
例1. 初始值为麦克斯韦等面积律的理论解所对应的气液相密度, 即

ρ_{g} = 0.4194

.
如图2所示, 虚线为初始密度和压力分布, 实线便是数值计算所得到的结果, 其具体数值为

ρ_{l} = 1.6550

(ρ_{ini, g} = 0.3, ρ_{ini, l} = 1.85)

, 与理论解的相对误差分别为 1.48%和0.14%. 其中, 液相密度与理论解相符较好, 而气相密度则偏差较大, 与文献[12]的结果是一致的. 这是由于气相密度较小, 故而液相发生很小的数值振荡便会引起气相的较大偏差. 此外, 在相界面附近, 压力发生了急剧变化, 这与物理实际是不相符的, 究其原因, 这是由于范德瓦尔斯方程并非真实的状态方程所导致的. 范德瓦尔斯方程在两相共存区与实际情况偏差较大, 其所描述的相变线并非实际的相变过程线. 相界面处的密度分布在两相之间, 由图1 可知, 其对应的压力分布或大于气液两相的共存压力, 或小于气液两相的共存压力, 故而会在相界面处出现急剧变化, 这与文献[36]中图1(b) 所示结果是一致的.

显示原图|下载原图ZIP|生成PPT
图2例1: 虚线为初始分布, 实线为麦克斯韦等面积律的数值解.
-->Fig.2Case1: The dotted line is the initial distribution. The solid line is the numerical result of Maxwell Construction.
-->

例2. 气液相初始密度均处于稳态, 具体值为

\frac{1}{ρ_{α}} < \frac{1}{ρ_{ini, g}} < \frac{1}{ρ_{β}}

. 其数值结果如图3 所示. 不同于例1的是Y= 0.5处由液相蒸发为气相, 这应该是由于初始液相的压力远大于平衡压力, 如图3(b) 所示, 因此随着时间的演化, 液相通过相界面将压力扰动传到气相, 故而远离相界面的液相(Y= 0.5)压力便会急剧下降, 从而变为气相. 其压力与两侧的液体平衡, 如图3(b)中实线所示. 如果降低初始液相的压力, 该现象便会消失. 总之, 无论最终存在几个相界面, 气液两相的压力是相等的, 且收敛在理论解附近.

显示原图|下载原图ZIP|生成PPT
图3例2: 虚线为初始分布, 实线为麦克斯韦等面积律的数值解.
-->Fig.3Case2: The dotted line is the initial distribution. The solid line is the numerical result of Maxwell Construction.
-->

例3. 气液相初始密度均处于完全不稳定椭圆域, 即

\frac{1}{ρ_{α}} < \frac{1}{ρ_{ini, l}} < \frac{1}{ρ_{β}}

, 并且

(ρ_{ini, g} = 0.8, ρ_{ini, l} = 1.2)

. 具体值为

? p / ? ρ

.
其数值结果如图4所示. 可以看到, 当气液两相初始密度均处于完全不稳定椭圆域时, 会通过蒸发和液化实现相分离, 并最终稳定在理论解附近.

显示原图|下载原图ZIP|生成PPT
图4例3: 虚线为初始分布, 实线为麦克斯韦等面积律的数值解.
-->Fig.4Case3: The dotted line is the initial distribution. The solid line is the numerical result of Maxwell Construction.
-->

以上3种初值所对应的数值解对比如表1所示, 结果表明, 数值解最终都可以回到适定的气液相密度, 即使麦克斯韦等面积律没有在格式中显式体现. 数值与理论的相对误差小于5%, 且网格界面可以用2到3个网格来捕捉, 从而证明该格式是稳定且有效的. 这也是由于范德瓦尔斯状态方程的物理特性所导致的. 对例2而言, 气液相初始密度均处于稳态. 由于两相压力不等, 液相压力大于理论解, 而气相压力小于理论解. 为了达到平衡, 如图1所示, 液相参数便沿着范德瓦尔斯曲线向下游变化, 与此同时, 气相参数沿着曲线向上游运动, 于是在理论解附近达到平衡. 再如例3所示, 当流体密度处于椭圆域时, 由于

p / ρ < 0

的斜率为负, 任何小的扰动都会被放大. 假设某个网格位于椭圆域, 其密度有一个正的扰动, 便会导致压力下降, 此时邻近网格的流体便会向该网格汇聚, 如此该网格的密度便会增加, 换言之, 正的扰动被放大. 故而该区域的流体混合物或蒸发为气相或凝结为液相. 而气相和液相的压力会在理论解附近达到平衡. 总之, 即使初始值不在理论解附近, 随着物理过程的演化, 也会收敛到理论解附近.
Table 1
表1
表1T = 0.9时不同初值的数值解对比
Table 1Numerical coexistence solution of different initial values at T = 0.9

Initial	Coexistence	Relative
densities	densities	error/%
(gas/liquid)	(gas/liquid)
0.425 7/1.657 3	0.419 4/1.655 0	1.48/0.14
0.3/1.85	0.404 4/1.653 7	5.00/0.22
0.8/1.2	0.421 3/1.655 9	1.03/0.08

新窗口打开
为了进一步确认该格式有效性, 计算不同温度下, 以理论解为初值的数值解, 绘制出气液两相的共存曲线, 并与理论解相比较.
例4. 如上所述, 分别计算T = 0.847 5, 0.86, 0.88, 0.90, 0.92, 0.94, 0.96, 0.98八个温度下的数值解, 并与理论解对比, 其结果如图5 所示. 温度1为临界情况, 对应的便是图5 中的顶点, 此时流体由单一的气相变为气液共存相. 温度越低, 共存气液两相的密度比越大, 故而共存曲线为抛物线型. 显然, 数值解与理论解相符良好.

显示原图|下载原图ZIP|生成PPT
图5范德瓦尔斯状态方程所描述的共存曲线.
-->Fig.5Coexistence curve of vdW equation of state.
-->

需要强调的是, 由于范德瓦尔斯方程的刚性, 由式(1)可知, 当T < 0.844时, 图1所示的范德瓦尔斯曲线上, 相界面处的压力会出现负值, 从而导致$p/\rho < 0$, 亦即

λ < 0

. 而格式的构造过程中, 计算通量时会出现

\sqrt[]{λ}

, 因此要求

λ > 0

. 因此利用该方程所构造的格式仅适用于T > 0.844的气液两相流计算. 通过理论解可知, T = 0.844 时的气液共存密度比为5. 随着温度的降低, 共存密度比会越大. 故而, 利用范德瓦尔斯方程构造的格式所适用的密度比小于5.

3.2 液滴表面张力

例5. 利用该格式计算范德瓦尔斯状态方程所描述的液滴表面张力, 验证其是否满足拉普拉斯定理. 对于UGKS而言, 计算域为1

\times

1, 温度T= 0.9, 速度为零. 初始的密度分布为[35]

\begin{array}{l} ρ (x, y) = & \frac{ρ_{ini, l} + ρ_{ini, g}}{2} - \frac{ρ_{ini, l} + ρ_{ini, g}}{2} . \\ \tanh \frac{2 (\sqrt[]{(x - x_{1})^{2} + (y - y_{1})^{2}} - R_{0})}{Υ} \end{array}

其中,

ρ_{ini, l} = 1.6573

ρ_{ini, g} = 0.4257

(x_{1}, y_{1})

,为计算域中心, 取为(0.5, 0.5),

R_{0}

为初始液滴半径.
根据拉普拉斯定理, 液滴内外的压力差与液滴半径成反比, 亦即

Δ P = P_{in} - P_{out} = σ / R

其中,

σ

为表面张力,

R

为液滴半径, 亦即液滴中心到半密度点的距离, 顾名思义, 半密度点的密度为

ρ = (ρ_{in} + ρ_{out}) / 2

Δ P

R_{0} = 0.1,0.15,0.2,0.25,0.3

为内外压差.
通过对

\times

五种情形在网格划分为180

1.75 \times 1 0^{- 3}, 1.2 \times 1 0^{- 3}, 8.5 \times 1 0^{- 4}, 6.7 \times 1 0^{- 4}, 6 \times 1 0^{- 4}

180时进行数值模拟可知, 液滴的内外压差分别为

σ = 1.75 \times 1 0^{- 4}

. 以上结果与文献[36] 中图3 所示的结果在量级上是一致的. 易得其值与半径成反比, 如图6 所示, 计算可得表面张力

\times

, 满足拉普拉斯定理, 从而进一步证明该格式的有效性.

显示原图|下载原图ZIP|生成PPT
图6液滴内外压差与半径之间的关系.
-->Fig.6The relationship between the pressure difference and the radius of droplet .
-->

从物理上来讲, 由于范德瓦尔斯状态方程所描述的分子之间存在范德瓦尔斯力, 亦即相互作用势, 导致在相界面处势函数的法向二阶导数不等于切向二阶导数, 从而根据表面张力的定义[37-38], 可知由于范德瓦尔斯力的存在, 在物理上液滴表面会出现相应的表面张力, 进一步证明了上述数值结果的可靠性. 事实上, 之前已有诸多文献采用范德瓦尔斯模型中的内聚力来处理表面张力[39-41]. 因为从微观角度看, 表面张力就是由于液体表面的分子所受上层空间气相分子引力小于内部液相分子对它的吸引力造成, 合力方向指向液体内部.

3.3 液滴碰撞

例6. 利用该格式计算两个液滴的碰撞, 计算域为1

(ρ_{ini, g} = 0.4257, ρ_{ini, l} = 1.6573)

1, 温度T = 0.9. 初始的气液密度分布为置为麦克斯韦等面积律所对应的理论解, 亦即

\times

. 计算域划分为 100{Invalid MML}100. 两个液滴的中心坐标分别为(0.38,0.5)和(0.62,0.5), 其半径均为0.1, 如图7(a) 所示. 两个液滴各自以U = 0.3 的初始速度对碰. 图7 为液滴形状随时间之演化, 可以看到液滴的碰撞和融合, 其过程与文献[33]相类似. 当两个液滴相向运动靠近时, 介于二者之间的气体就会被压缩形成一层气膜, 同时压强升高, 液滴发生微小变形. 因为液滴运动缓慢, 气膜内的气体在液滴接触之前有足够的时间流出, 所以液滴间就会发生微小变形后聚合 [42]. 由于耗散的作用, 液滴最终形态与标准的圆形还有些差异 [39]. 此外, 如图8 所示, 时间t = 1.6时所得X和Y方向的中心线的密度分布中, 两个网格点便可以捕捉一个尖锐的相界面.

显示原图|下载原图ZIP|生成PPT
图7两个液滴碰撞的演化过程.
-->Fig.7Time evolution of the collision of two droplets.
-->

显示原图|下载原图ZIP|生成PPT
图8密度分布: (a)图7(f)中沿X方向中心线, (b)图7(f)中沿Y方向中心线.
-->Fig.8The density distribution (a) along the central line of Fig.7 (f) in the X direction and (b) along the central line of Fig.7 (f) in the Y direction.
-->

需要强调的是, 该格式没有计算能量的演化, 亦即仅关注了相变所带来的传质问题, 而没有考虑相应的传热问题, 这也是下一步需要努力的方向.

4 结论

本文将范德瓦尔斯状态方程所表达的非理想气体效应引入统一气体动理学格式之中来捕捉气液相界面. 利用该方法得到麦克斯韦等面积律的数值解, 并与其对应的理论解相比较, 二者符合较好. 而后, 通过对范德瓦尔斯状态方程所描述的液滴表面张力进行数值计算, 验证了拉普拉斯定理. 此外, 通过模拟两个液滴的碰撞过程, 进一步证明了该格式的有效性. 但是, 由于范德瓦尔斯状态方程的特性, 其所构造的格式仅适用于液/气两相密度比小于5的情况. 下一步将致力于构造能够适用于大密度比(> 1000)计算的数值格式.
The authors have declared that no competing interests exist.

参考文献原文顺序
文献年度倒序
文中引用次数倒序
被引期刊影响因子

[1]	刘文超, 刘曰武. 低渗透煤层气藏中气?水两相不稳定渗流动态分析 . 力学学报, 2017, 49(4): 828-835 [本文引用: 1] (Liu Wenchao, Liu Yuewu.Dynamic analysis on gas-water two-phase unsteady seepage flow in low-permeable coalbed gas reservoirs .Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(4): 828-835 (in Chinese)) [本文引用: 1]
[2]	李康, 刘娜, 何志伟等. 一种基于双界面函数的界面捕捉方法 . 力学学报, 2017, 49(6): 1290-1300 (Li Kang, Liu Na, He Zhiwei, et al.A new interface capturing method based on double interface functions .Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(6): 1290-1300 (in Chinese))
[3]	Hirt CW, Nichols BD.Volume of fluid (VOF) method for the dynamics of free boundaries .J Comput Phys, 1981, 39(1): 201-225 [本文引用: 1]
[4]	张洋, 陈科, 尤云祥等. 壁面约束对裙带气泡动力学的影响 . 力学学报, 2017, 49(5): 1050-1058 (Zhang Yang, Chen Ke, You Yunxiang, et al.Confinement effect on the rising dynamics of a skirted bubble .Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(5): 1050-1058 (in Chinese))
[5]	Sussman M, Smereka P, Osher S.A level set approach for computing solutions to incompressible two-phase flow .J Comput Phys, 1994, 114(1): 146-159 [本文引用: 1]
[6]	郭照立, 郑楚光. 格子Boltzmann 方法的原理及应用. 北京: 科学出版社, 2009 [本文引用: 1] (Guo Zhaoli, Zheng Chuguang.Principle and Application of Lattice Boltzmann Method. Beijing: Science Press, 2009 (in Chinese)) [本文引用: 1]
[7]	Shan X, Chen H.Lattice Boltzmann model for simulating flows with multiple phases and components .Phys Rev E, 1993, 47(3): 1815-1819 [本文引用: 2]
[8]	Shan X, Chen H.Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation .Phys Rev E, 1994, 49(4): 2941-2948
[9]	Gunstensen AK, Rothman DH, aleski SZ, et al. Lattice Boltzmann model of immiscible fluids .Phys Rev A, 1991, 43(8): 4320-4327
[10]	Swift MR, Osborn WR, Yeomans JM.Lattice Boltzmann simulation of nonideal fluids .Phys Rev Lett, 1995, 75(5): 830-833
[11]	Zhang R, Chen H.Lattice Boltzmann method for simulations of liquid-vapor thermal flows .Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics, 2003, 67(6): 66711 [本文引用: 1]
[12]	Yuan, P, Schaefer L. Equations of state in a lattice Boltzmann model .Phys Fluids, 2006, 18(4): 42101 [本文引用: 2]
[13]	Kupershtokh AL, Karpov DI, Medvedev DA, et al.Stochastic models of partial discharge activity in solid and liquid dielectrics .IET Science, Measurement & Technology, 2007, 1(6): 303-311 [本文引用: 2]
[14]	胡安杰. 多相流动格子Boltzmann 方法研究. [博士论文] . 重庆: 重庆大学, 2015 [本文引用: 1] (Hu Anjie.Study on lattice Boltzmann method for multiphase flow. [PhD Thesis] . Chongqing: Chongqing University, 2015 (in Chinese)) [本文引用: 1]
[15]	Li Q, Luo KH, Kang QJ, et al.Lattice Boltzmann methods for multiphase flow and phase-change heat transfer .Progress in Energy and Combustion Science, 2016, 52: 62-105
[16]	Chikatamarla SS, Karlin IV.Entropic lattice Boltzmann method for multiphase flows .Physical Review Letters, 2015, 114(17): 174502
[17]	Wang Y, Shu C, Huang HB, et al.Multiphase lattice Boltzmann flux solver for incompressible multiphase flows with large density ratio .Journal of Computational Physics, 2015, 280: 404-423
[18]	Liu H, Kang Q, Leonardi CR, et al.Multiphase lattice Boltzmann simulations for porous media applications .Computational Geosciences, 2016, 20(4): 777-805
[19]	Lycett-Brown D, Luo KH.Improved forcing scheme in pseudopotential lattice Boltzmann methods for multiphase flow at arbitrarily high density ratios .Physical Review E, 2015, 91(2): 023305
[20]	Shao JY, Shu C.A hybrid phase field multiple relaxation time lattice Boltzmann method for the incompressible multiphase flow with large density contrast .International Journal for Numerical Methods in Fluids, 2015, 77(9): 526-543
[21]	郭宇隆. 基于格子Boltzmann 方法的气液混合流体模拟. [硕士论文] . 广州: 华南理工大学, 2015 (Guo Yulong.Study on gas-liquid mixed fluid simulation based on lattice Boltzmann method [Master Thesis] . Guangzhou: South China University of Technology, 2015 (in Chinese))
[22]	史冬岩, 王志凯, 张阿漫. 一种模拟气液两相流的格子波尔兹曼改进模型 . 力学学报, 2013, 46(2): 224-233 (Shi Dongyan, Wang Zhikai, Zhang Aman.Improved model of simulating gas-liquid two-phase flow by lattice Boltzmann .Chinese Journal of Theoretical and Applied Mechanics, 2013, 46(2): 224-233 (in Chinese))
[23]	曾建邦, 李隆键, 廖全等. 格子Boltzmann 方法在相变过程中的应用 . 物理学报, 2010, 59(1): 178-185 (Zeng Jianbang, Li Longjian, Liao Quan, et al.Application of lattice Boltzmann method in phase transition process .Acta Physica Sinica, 2010, 59(1): 178-185 (in Chinese))
[24]	Gong S, Cheng P.A lattice Boltzmann method for simulation of liquid-vapor phase-change heat transfer .International Journal of Heat and Mass Transfer, 2012, 55(17-18): 4923-4927 [本文引用: 1]
[25]	Dong L, Gong S, Cheng P.Direct numerical simulations of film boiling heat transfer by a phase-change lattice Boltzmann method .International Communications in Heat and Mass Transfer, 2018, 91: 109-116 [本文引用: 2]
[26]	Guo Z, Zhao T S, Shi Y.Physical symmetry, spatial accuracy, and relaxation time of the lattice Boltzmann equation for micro gas flows .Journal of Applied Physics, 2006, 99(7): 074903 [本文引用: 1]
[27]	Zhang Y, Qin R, Emerson DR.Lattice Boltzmann simulation of rarefied gas flows in microchannels .Physical Review E, 2005, 71(4): 047702
[28]	Guo Z, Zheng C, Shi B.Lattice Boltzmann equation with multiple effective relaxation times for gaseous microscale flow .Physical Review E, 2008, 77(3): 036707
[29]	Shan X, Yuan XF, Chen H.Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation .Journal of Fluid Mechanics, 2006, 550(1): 413-441
[30]	Xu K, Huang JC.A unified gas-kinetic scheme for continuum and rarefied flows .Journal of Computational Physics, 2010, 229(20): 7747-7764
[31]	毛枚良, 江定武, 李锦等. 气体动理学统一算法的隐式方法研究 . 力学学报, 2015, 47(5): 822-829 (Mao Meiliang, Jiang Dingwu, Li Jin, et al.Study on implicit implementation of the unified gas kinetic scheme .Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(5): 822-829 (in Chinese))
[32]	江定武, 毛枚良, 李锦等. 气体动理学统一算法中相容性条件不满足引起的数值误差及其影响研究 . 力学学报, 2015, 47(1):163-168 (Jiang Dingwu, Mao Meiliang, Li Jin, et al.Study on the numerical error introduced by dissatisfying the conservation constraint in UGKS and its effects .Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(1): 163-168 (in Chinese))
[33]	Xu K.A kinetic method for hyperbolic-elliptic equations and its application in two-phase flow .Journal of Computational Physics, 2001, 166(2): 383-399 [本文引用: 2]
[34]	Wang Z, Yan H, Li Q, et al.Unified gas-kinetic scheme for diatomic molecular flow with translational, rotational, and vibrational modes .Journal of Computational Physics, 2017, 350: 237-259 [本文引用: 1]
[35]	Huang H, Sukop M, Lu X.Multiphase Lattice Boltzmann Methods: Theory and Application . Beijing: John Wiley & Sons, 2015 [本文引用: 1]
[36]	Hu A, Li L, Chen S, et al.On equations of state in pseudo-potential multiphase lattice Boltzmann model with large density ratio .International Journal of Heat and Mass Transfer, 2013, 67: 159-163 [本文引用: 3]
[37]	Hu A, Li L, Uddin R.Force method in a pseudo-potential lattice Boltzmann model .Journal of Computational Physics, 2015, 294: 78-89 [本文引用: 1]
[38]	Shan X.Pressure tensor calculation in a class of nonideal gas lattice Boltzmann models .Physical Review E, 2008, 77(6): 066702
[39]	强洪夫, 陈福振, 高巍然. 基于SPH 方法的低韦伯数下三维液滴碰撞聚合与反弹数值模拟研究 . 工程力学, 2012, 29(2): 21-28 [本文引用: 1] (Qiang Hongfu, Chen Fuzhen, Gao Weiran.Simulation of coalescence and bouncing of three-dimensional droplet collisions with low weber numbers based on SPH method .Journal of Engineering Mechanics, 2012, 29(2): 21-28 (in Chinese)) [本文引用: 1]
[40]	Nugent S, Posch HA.Liquid drops and surface tension with smoothed particle applied mechanics .Physical Review E, 2000, 62: 4968-4975
[41]	Meleán Y, Sigalotti LDG.Coalescence of colliding van der Waals liquid drops .International Journal of Heat and Mass Transfer, 2005, 48: 4041-4061
[42]	夏盛勇. 三氧化二铝液滴碰撞机理及模型研究. [博士论文] . 西安: 西北工业大学, 2015 [本文引用: 1] (Xia Shengyong.Physics and model of alumina droplet collisions. [PhD Thesis] . Xi’an: Northwestern Polytechnical University, 2015 (in Chinese)) [本文引用: 1]

基于气液相界面捕捉的统一气体动理学格式

本站小编 Free考研考试/2022-01-01

UNIFIED GAS-KINETIC SCHEME FOR TWO PHASE INTERFACE CAPTURING

1 问题描述