格子Boltzmann模型结合MPR方程模拟流体饱和气液密度

引用本文

闵琪, 段远源, 王晓东. 格子Boltzmann模型结合MPR方程模拟流体饱和气液密度. 清华大学学报: 自然科学版, 2014, 54(5): 619-623[Qi MIN, Yuanyuan DUAN, Xiaodong WANG. Lattice Boltzmann method for the fluid saturation density based on the volume translated Peng-Robinson equation of state[J]. 2014, 54(5): 619-623] 复制到剪切板

Permissions

格子Boltzmann模型结合MPR方程模拟流体饱和气液密度

闵琪¹, 段远源², 王晓东³

1. 清华大学核能与新能源技术研究院,北京 100084

2.清华大学热科学与动力工程教育部重点实验室,北京 100084

3.华北电力大学新能源电力系统国家重点实验室,北京 102206

作者简介: 闵琪(1985-), 女(汉), 江西, 助理研究员。E-mail:minq86@mail.tsinghua.edu.cn

基金:国家自然科学基金资助项目 (51206093,21176133);

摘要

格子Boltzmann伪势能两相模型由于缺少与真实流体相关的物理参数,因此无法实现对某种特定流体的模拟。本文将比容平移后的Peng-Robinson状态方程(MPR)引入模型中,通过状态方程中的无量纲参数偏心因子和临界压缩因子,建立模拟流体与真实流体的关联,使得模型可以区别各种真实流体,并模拟了氩、氧气、烷烃等8种流体在不同温度下的饱和气液相密度。结果表明,MPR和Peng-Robinson方程(PR)的格子Boltzmann两相模型均能很好的描述8种流体的饱和气相密度; 而MPR方程能够很好的再现氩、氮气、氧气等非极性流体的饱和液相密度,对于其它流体,MPR方程较PR方程对液相密度的描述有所改进,但仍与实验数据有一定差异。总体上MPR方程能够更好地模拟不同流体的饱和气液相密度。

关键词: 格子Boltzmann; 气液两相流; 状态方程; 比容平移; 饱和液相密度

中图分类号:TK121 文献标志码:A 文章编号:1000-0054(2014)05-0619-05

Lattice Boltzmann method for the fluid saturation density based on the volume translated Peng-Robinson equation of state

Qi MIN¹, Yuanyuan DUAN², Xiaodong WANG³

1. Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China

2. Key Laboratory for Thermal Science and Power Engineering of the Ministry of Education, Tsinghua University, Beijing 100084, China

3. State Key Laboratory of Alternate Electrical Power System with Renewable Energy Sources, North China Electric Power University, Beijing 102206, China

Fund:

Abstract

The lattice Boltzmann pseudo-potential model has few parameters for real fluids. In this study, the volume translated Peng-Robinson equation of state (MPR) was incorporated into the lattice Boltzmann pseudo-potential model with two parameters, the acentric factor and the critical compressibility, used to determine the real fluid parameters. The MPR was used to calculate the dimensionless liquid and vapor saturated densities of eight kinds of fluids with the results compared with the Peng-Robinson(PR) equation of state and experimental data. Both the MPR and the PR equation of state can accurately describe the dimensionless saturated vapor density. The simulation results for the saturated liquid density given by the MPR agree well with experimental data for simple non-polar fluids such as Ar, N₂ and O₂. For more complex fluids, the saturated liquid density simulated by both the MPR and the PR equation of state were not good, but the MPR result was better than the PR equation of state result. Generally speaking,the MPR gives better predictions than the PR equation of state for simulations of fluid saturated densities.

Keyword: lattice Boltzmann; liquid vapor flow; equation of state; volume translation; saturation density

Show Figures

气液相变和两相流在自然界以及工程实际中极为常见,其物理过程复杂,数值模拟方法是揭示其物理规律的有效方法。近年来得到广泛发展的格子Boltzmann方法以颗粒微团为模拟单元,利用非平衡态统计物理方法,计算效率高,模拟尺度接近宏观实际,同时可以反映流体运动的微观机制和运动规律。格子Boltzmann方法中用于处理两相流动问题的模型包括^[1]: 伪势能模型法^[2]、自由能法^[3,4]和颜色模型法^[5,6]。其中由Shan和Chen^[7]提出的伪势能模型法,参照了分子间势能函数的思想,直接对颗粒间相互作用力提出势能模型,该方法能够反映多相流体动力学的物理本质,模型物理意义明确,得到了广泛的应用。

格子Boltzmann伪势能模型中,颗粒间相互作用势能模型 ψ( x)决定了流体所满足状态方程的形式,因此 ψ( x)是控制气液两相性质的关键物理关系。学者曾提出多种不同形式的势能函数 ψ( x), 如 ψ( x) = ρ₀[1-exp( -ρ/ρ₀)]; ψ( x) =ρ; ψ( x) =g $\begin{matrix} {ρ_{0}}^{2} \end{matrix}$ ρ² /[2( ρ₀ ₊ρ)²]^{[7,8,9,10,11,12]}等形式,其中 ρ为流体密度, ρ₀和 g是模型参数。但这些势能函数所衍生出的状态方程中没有与真实流体相关的物理量,因此模型无法模拟某种特定流体。也有学者提出将立方型状态方程、维里方程等引入格子Boltzmann模型中。Yuan等^[13]对van der Waals (VDW)、 Redlich-Kwong (RK)、 Redlich-Kwong Soave (RKS)、 Peng-Robinson (PR)和Carnahan-Starling (CS)状态方程在格子Boltzmann模型中的应用进行了大量对比研究。其中PR和RKS状态方程含有偏心因子ω, 由于该参数为无量纲参数,可以取流体的实际值,因此ω将模型与真实流体关联起来,通过改变ω, 可以模拟不同流体的气液相平衡状态。

然而引入PR或RKS状态方程后的格子Boltzmann模型仅能识别流体的偏心因子,对于其它的流体性质仍无法区别,因此在格子Boltzmann两相模型中引入更多的与真实流体相关的物理量对模型的进一步发展具有重要意义。Lin等^[14]的研究结果表明,通过对PR方程引入比容平移量,能够较好地修正PR状态方程对流体饱和液相密度描述不精确的缺陷,该比容平移量取决于流体无量纲参数临界压缩因子。

本文将结合比容平移后的PR状态方程(MPR), 建立与真实流体偏心因子和临界压缩因子同时相关的格子Boltzmann两相模型,并用该模型模拟真实流体的饱和气液相密度,与PR方程的模拟结果进行比较,同时与实验数据进行对比。

1 格子Boltzmann方法理论基础

格子Boltzmann方法由分布函数演化方程、平衡态分布函数和离散速度模型三部分组成。在无外力情况下,BGK (Bhantagar-Gross-Krood)^[15]近似的格子Boltzmann方程为:

$\begin{matrix} \begin{matrix} f_{i} (x + e_{i} Δt, t + Δt) \\ = f_{i} (x, t) - \frac{Δt}{τ_{v}} [f_{i} (x, t) - f_{i}^{eq} (x, t)] . \end{matrix} \end{matrix}$ (1)

其中: f_i( x, t)是空间中 x点在 t时刻的密度分布函数; 格子的长度 Δx定义为1 lu, 离散时间 Δt定义为1 ts; e_i为粒子离散速度在 i方向上的分量; $\begin{matrix} \frac{Δt}{τ_{v}} \end{matrix} \begin{matrix} [f_{i}^{} (x, t) - f_{i}^{eq} (x, t)] \end{matrix}$ 是BGK碰撞项; τ_v是松弛时间。

演化方程中 e_i由速度离散模型确定,二维九速(D2Q9)模型离散速度如下:

在D2Q9模型中平衡态密度分布函数应用了Maxwell分布的展开形式^[16]

$\begin{matrix} \begin{matrix} f_{i}^{eq} (x, t) = \\ w_{i} ρ [1 + \frac{e_{i} \cdot u}{c_{s}^{2}} + \frac{{(e_{i} \cdot u)}^{2}}{2 c_{s}^{4}} - \frac{{(u)}^{2}}{2 c_{s}^{2}}] . (3) \end{matrix} \end{matrix}$

其中: u为速度, c_s为格子声速, w_i为权重系数。当 i=0时, w_i=4 /9; 当 i=1,2,3,4时, w_i =1 /9; 当 i=5,6,7,8时, w_i =1 /36。

Shan和Chen^[7]提出的伪势能模型中微团间相互作用力 F( x, t)提出势能模型:

其中: ψ( x)为 x处的势能函数, G为模型参数。势能函数控制的状态方程满足:

$\begin{matrix} P = ρ / 3 + G / 6 {[ψ (ρ)]}^{2} . (5) \end{matrix}$

2 MPR状态方程

PR方程^[17]的原始形式如下:

$\begin{matrix} p = \frac{RT}{v - b} - \frac{aα (T)}{v^{2} + 2 bv - b^{2}} . (6) \end{matrix}$

其中: p为压力; R通用气体常数; T为热力学温度; 为比容; 和b为状态方程系数

$\begin{matrix} a = 0.45724 \frac{R^{2} T_{c}^{2}}{p_{c}}; b = 0.0778 \frac{R T_{c}}{p_{c}} . (7) \end{matrix}$

其中: T _c和p _c分别为临界温度和压力。α(T)函数的定义如下:

$\begin{matrix} α (T) = {[1 + m (1 - T_{r}^{0.5})]}^{2} . (8) \end{matrix}$

其中: T _r=T/ _c为对比温度。

当偏心因子ω≤0.49时, m为:

$\begin{matrix} m = 0.3765 + 1.5423 ω - 0.2699 ω^{2} . (9) \end{matrix}$

当偏心因子ω>0.49时, m为:

$\begin{matrix} m = 0.3796 + 1.4850 ω - 0.1644 ω^{2} + 0.0167 ω^{3} . (10) \end{matrix}$

应用比容平移理论,引入比容平移量c, MPR状态方程如下:

$\begin{matrix} p = \frac{RT}{v + c - b} - \frac{a}{{(v + c)}^{2} + 2 b (v + c) - b^{2}} . (11) \end{matrix}$

其中c为在临界点的修正c _c与温度相关的函数f(t)的乘积形式^[14]:

$\begin{matrix} \begin{matrix} c = c_{c} f (t), \\ c_{c} = {(v_{PR})}_{c} - {(v_{\exp})}_{c} = (0.3074 - Z_{c}) \frac{R T_{c}}{p_{c}}, \\ f (t) = β + (1 - β) \exp (0.5 χ), β = - 2.8431 \\ \exp [- 64.2184 (0.3074 - Z_{c})] + 0.1735, \\ χ = - 99.2558 + 301.6201 Z_{c} . (12) \end{matrix} \end{matrix}$

3 气液相变模拟结果

本文采用周期边界条件模拟了50×50格子区域内的流体相变过程,分别结合 PR方程和 MPR方程模拟流体的相变过程,格子单位下 PR和 MPR方程的共同参数为R=1, a=2/49, b=2/21^[13]。可以计算得到临界参数T _c=0.07292, P _c=0.059 57, ρ _c=2.657 3。通过设定流体的初始密度,使得流体处于亚稳态,受微小扰动后发生两相分离,最终形成稳定的气液相平衡状态,本文中初始密度ρ₀=2。模拟中认为t时刻与t-1000时刻的流场最大速度之差在10^-6量级时,系统达到稳定状态。

本模拟中选取了8种流体: 氩( Ar)、氮气( N₂)、氧气( O₂)、二氧化碳( CO₂)、甲烷( CH₄)、乙烷( C₂ H₆)、丙烷( C₃ H₈)、正丁烷( C₄ H₁₀)。图1 a-1 h为8种流体的饱和气液相密度随温度变化关系的模拟结果和实验结果的对比,横纵坐标采用的是无量纲化后的对比密度和对比温度。对于 PR方程, ρ _r=ρ/ρ _c; 对于 MPR方程,需考虑比容平移量,ρ _r=ρ/(ρ _c+1/ c)。从图中可以看出,引入 PR或 MPR方程的格子 Boltzmann模型都能够较好地模拟流体的气相密度。但 PR方程对液相密度的模拟结果与实验结果有较大的差异,而 MPR方程对液相密度的计算明显较 PR方程更接近实验结果。表1列出了 MPR方程和 PR方程的液相密度模拟结果与实验值之间的绝对偏差,其中 $\begin{matrix} {ρ_{r}}^{sim} \end{matrix}$ 为模拟结果, $\begin{matrix} {ρ_{r}}^{\exp} \end{matrix}$ 为实验值。可以看出对于氩、氮气、氧气类分子结构较简单的非极性流体, MPR方程能够很好地再现流体的饱和液相密度,而对于其它流体, MPR方程虽然较 PR方程对液相密度的描述有所改进,但仍与实验数据有一定差异。

	Figure Option View Download New Window Download As Powerpoint Slide
	图1 PR和MPR状态方程控制下流体饱和气液密度的模拟结果与实验值比较

表1 MPR和PR方程饱和液相密度模拟绝对偏差比较

δ_a	Ar	N₂	O₂	CO₂
MPR	0.046	0.056	0.082	0.155
PR	0.363	0.397	0.438	0.298
δ_a	CH₄	C₂H₆	C₃H₈	C₄H₁₀
MPR	0.082	0.186	0.168	0.264
PR	0.394	0.484	0.474	0.469

注: δ_a= $\begin{matrix} {ρ_{r}}^{sim} \end{matrix}$ - $\begin{matrix} {ρ_{r}}^{\exp} \end{matrix}$

表1 MPR和PR方程饱和液相密度模拟绝对偏差比较

对于 PR方程,偏心因子数ω是模型中表征不同流体的唯一无量纲参数; 而 MPR方程有2个与实际流体关联的无量纲参数偏心因子ω和临界压缩因子Z _c,提高了格子 Boltzmann两相模型对流体的识别能力,并且通过引入比容平移量,提高了格子 Boltzmann两相模型对流体饱和液相密度的模拟能力。

4 结论

本文通过在格子Boltzmann伪势能两相模型基础上引入MPR状态方程,建立了新的格子Boltzmann两相模型,该两相模型中含有两个与真实流体关联的无量纲参数偏心因子和临界压缩因子,提高了格子Boltzmann模型与真实流体的相关性。模拟结果表明,结合MPR状态方程的格子Boltzmann模型比结合PR状态方程的格子Boltzmann模型能够更好地描述流体饱和液相密度。并且MPR方程对于氩、氮气、氧气等分子结构较简单的非极性流体的饱和液相密度的模拟结果较其他分子结构较复杂流体更接近实验值。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	郭照立, 郑楚光, 李青. 流体动力学的格子Boltzmann方法 [M]. 北京: 科学出版社, 2009. GUO Zhaoli, ZHENG Chuguang, LI Qing. Theory and Applications of Lattice Boltzmann Method [M]. Beijing: Science Press, 2009. (in Chinese) [本文引用:1]
[2]	李维仲, 李爽. 用格子Boltzmann方法模拟液滴合并过程[J]. 热科学与技术, 2007, 6(3): 379-393. LI Weizhong, LI Shuang. Simulation of droplets coalescence process by lattice Boltzmann method[J]. Journal of Thermal Science and Technology, 2007, 6(3): 379-393. (in Chinese) [本文引用:1] [CJCR: 0.511]
[3]	Swift M R, Orlandini E, Osborn W R, et al. Lattice Boltzmann simulations of liquid-gas and binary fluid systems[J]. Physical Review. E, 1996, 54(5): 5041-5052. [本文引用:1] [JCR: 2.326]
[4]	Inamuro T, Ogata T, Tajima S, et al. A lattice Boltzmann method for incompressible two-phase flows with large density differences[J]. Journal of Computational Physics, 2004, 198(2): 628-644. [本文引用:1] [JCR: 2.485]
[5]	Grunau D, Chen S Y, Eggert K. A lattice Boltzmann model for multiphase fluid-flows[J]. Physics of Fluids. A, Fluid Dynamics, 1993, 5(10): 2557-2562. [本文引用:1]
[6]	Gunstensen A K, Rothman D H, Zaleski S, et al. Lattice Boltzmann Model of Immiscible Fluids[J]. Physical Review. A, 1991, 43(8): 4320-4327. [本文引用:1] [JCR: 2.991]
[7]	Shan X W, Chen H D. Lattice boltzmann model for simulating flows with multiple phases and components[J]. Physical Review. E, 1993, 47(3): 1815-1819. [本文引用:3] [JCR: 2.326]
[8]	Raiskinmaki P, Koponen A, Merikoski J, et al. Spreading dynamics of three-dimensional droplets by the lattice- Boltzmann method[J]. Computational Materials Science, 2000, 18(1): 7-12. [本文引用:1] [JCR: 1.879]
[9]	Hyvaluoma J, Raiskinmaki P, Jasberg A, et al. Evaluation of a lattice-Boltzmann method for mercury intrusion porosimetry simulations[J]. Future Generation Computer Systems, 2004, 20(6): 1003-1011. [本文引用:1] [JCR: 2.639]
[10]	Martys N S, Chen H D. Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method[J]. Physical Review. E, 1996, 53(1): 743-750. [本文引用:1] [JCR: 2.326]
[11]	Pan C, Hilpert M, Miller C T. Lattice-Boltzmann simulation of two-phase flow in porous media[J]. Water Resources Research, 2004, 40(1). [本文引用:1] [JCR: 3.709]
[12]	Qian Y H, Orszag S A. Scalings in diffusion-driven reaction a+B-]C-numerical simulations by lattice BGK models[J]. Journal of Statistical Physics, 1995, 81(1-2): 237-253. [本文引用:1] [JCR: 1.284]
[13]	Yuan P, Schaeter L. Equations of state in a lattice Boltzmann model[J]. Physics of Fluids, 2006, 18(4): 042101. [本文引用:2] [JCR: 2.04]
[14]	林鸿. 气液界面张力的密度梯度理论模型和实验研究 [D]. 北京: 清华大学, 2006. LIN Hong. Gradient Theory Modeling and Experimental Investigation of the Surface Tension [D]. Beijing: Tsinghua University, 2006. (in Chinese) [本文引用:2]
[15]	Bhatnagar P L, Gross E P, Krook M. A model for collision processes in gases . 1. Small amplitude processes in charged and neutral one-component systems[J]. Physical Review, 1954, 94(3): 511-525. [本文引用:1]
[16]	Sukop M C, Thorne D T. Lattice Boltzmann Modeling[M]. New York, USA: Springer, 2006. [本文引用:1]
[17]	Peng D Y, Robinson D B. A new two-constant equation of state[J]. Industrial and Engineering Chemistry, Fundamentals, 1976, 15(1): 59-64. [本文引用:1]
[18]	Tegeler C, Span R, Wagner W. A new equation of state for argon covering thefluid region for temperatures from the melting line to 700 K at pressures up to 1000 MPa[J]. Journal of Physical and Chemical Reference Data, 1999, 28(3): 779-850. [本文引用:1] [JCR: 3.108]
[19]	Dewan R K, Mehta S K. Correlation between topological features and surface tension of binary liquid mixtures[J]. Monatshefte für Chemie, 1990, 121(8-9): 593-600. [本文引用:1] [JCR: 1.347]
[20]	Setzmann U, Wagner W. A new equation of state and tables of thermodynamic properties for methane covering the range from the melting line to 625 K at pressures up to 1000 MPa[J]. Journal of Physical and Chemical Reference Data, 1991, 20(6): 1061-1151. [本文引用:1] [JCR: 3.108]
[21]	Friend D G, Ingham H, Ely J F. Thermophysical properties of ethane[J]. Journal of Physical and Chemical Reference Data, 1991, 20(2): 275-347. [本文引用:1] [JCR: 3.108]
[22]	Miyamoto H, Watanabe K. A thermodynamic property model for fluid-phase propane[J]. International Journal of Thermophysics, 2000, 21(5): 1045-1072. [本文引用:1] [JCR: 0.623]
[23]	Miyamoto H, Watanabe K. Thermodynamic property model for fluid-phase n-Butane[J]. International Journal of Thermophysics, 2001, 22(2): 459-475. [本文引用:1] [JCR: 0.623]

2009

0.0