Interface Compression Technique in PPM
-
摘要: 高精度多组分分段抛物线法（Piecewise Parabolic Method，PPM）在对可压缩多相流问题进行模拟计算时，在不同组分交界面上存在界面扩散。为此，通过引入包含界面压缩和密度修正的人工界面压缩方法，抑制界面扩散现象。采用一个界面函数表示运动的物质界面，在多组分质量守恒方程和输运方程中添加考虑人工压缩和人工黏性的压缩源项，并在伪时间内采用二阶中心差分法和两步Runge-Kutta方法进行离散求解，采用Strang型分裂格式实现了整体算法的时间二阶精度。一维与二维数值模拟试验表明，结合人工界面压缩之后的PPM能有效抑制界面上数值扩散问题，在长时间的数值模拟中，人工界面压缩能够将扩散界面厚度维持在一定网格之内且保持界面形状不改变，尤其对于涉及稀疏波的问题，如激波引起的水中气泡坍塌，界面压缩效果更为显著。Abstract: This paper describes an artificial interface compression technique for the multi-fluid piecewise parabolic method (PPM). The proposed approach enables the simulation of interfaces between compressible multi-fluid flows with high density ratios and strong shock waves. A compression source term incorporated both interface compression and density correction is added to the mass conservation equation. The compression source term is solved in pseudo-time steps using the interface compression technique and the advection part is solved by multi-fluid PPM. The Strang splitting algorithm achieves second-order accuracy by combining the solutions of the advection operator and the interface compression operator. Numerical tests on the interaction of shock waves with interfaces in compressible multi-fluid flows reveal that multi-fluid PPM combined with the artificial interface compression technique can effectively prevent the smearing phenomenon, which is often observed at the contact interface. For long-time simulations, artificial interface compression with interface sharpening can constrain the thickness of the diffused interface to a few cells and maintain the interface profile. This artificial interface compression technique works well with multi-fluid PPM and the effect is obvious. It is a significant step in the accurate simulation of the collapse of air cavities in water, which involves strong rarefaction waves.
-
近年来随着液化天然气(LNG，Liquefied Natural Gas)产业在全球迅速发展，天然气的液化技术和设备也在不断发展完善、日渐成熟。我国LNG领域内的相关研究起步较晚，许多技术远远落后发达国家水平，在天然气液化工艺及装置的生产等方面缺乏自主产权。因此，开展天然气液化工艺及装置的研究，对于实现液化装置的国产化、高效化有十分重要的意义[1-3]。
超声速旋流分离技术是一种新兴的天然气加工处理技术，被较为广泛地用于天然气脱水、脱重烃、脱酸等方面，近年来开始逐渐应用于天然气液化方面[4-7]。天然气超声速液化的原理是：高压天然气混合物在Laval喷管内达到一定的温度、压力条件，开始凝结成核，最终凝结成液滴，在后续工艺中进行进一步气液分离。与传统的天然气液化技术相比，具有结构工艺简单、支持无人操作(适用于海底天然气处理)、对水合物抑制剂依赖性小、投资和运行成本低等优势[8-9]。
为了探究入口复杂多变的压力条件对天然气超声速液化特性的影响，对甲烷-乙烷气体混合物的超声速凝结流动特性进行研究，在凝结成核与生长理论的基础上建立了适用于甲烷-乙烷双可凝气体混合物的凝结流动数学模型，重点研究了入口压力对天然气混合物在Laval喷管内主要流动与凝结参数的影响规律。
1. Laval喷管结构设计
Laval喷管结构主要包括入口段、渐缩段、喉部及扩张段4部分[10-11]。各部分参数如表 1所示，L0为入口长度，r1为渐缩段入口半径，L1为渐缩段长度，rcr为喉部截面半径，L2为渐扩段长度，r2为渐扩段出口半径。为尽量减小流场涡流的影响，渐缩段采用双三次曲线设计，喉部采用一段平缓光滑的圆弧作为过渡曲线，渐扩段采用等膨胀率设计，膨胀率取为10 000 s-1。考虑到实验加工方便，保证曲面的精度，且能够更加直观地观察Laval喷管内部的流场分布情况，所设计Laval喷管截面采用矩形截面，三维结构如图 1所示。
表 1 Laval喷管各部分参数Table 1. Parameters of Laval nozzleL0/mm r1/mm L1/mm rcr/mm L2/mm r2/mm 50.00 20.00 56.01 2.50 71.28 6.15 2. 超声速凝结流动数学模型及计算方法
2.1 数学模型
采用欧拉双流体模型开展数值计算，控制方程主要包括气相流动方程和液相流动方程。在无滑移假设及欧拉双流体模型的前提下分别建立气相及液相流动控制方程组，液滴数目守恒方程及液滴半径、数目、湿度关系式分别添加到对应源相方程中[12-14]。
气相流动控制方程组
∂ρv∂t+∂∂xj(ρvuj)=Sm (1) ∂ρc2∂t+∂∂xj(ρc2uj)=Sm,c2 (2) ∂∂t(ρvui)+∂∂xj(ρvujui)=−∂pv∂xi+∂∂xj[μ(∂uj∂xi+∂ui∂xj−23δij∂uj∂xj)]+∂∂xj(−ρv¯u′iu′j)+Su (3) ∂∂t(ρvE)+∂∂xj(ρvujE+ujpv)=∂∂xj(keff∂T∂xj+uiτeff)+Sh (4) 液相流动控制方程组
∂∂t(ρY)+∂∂xj(ρujY)=SY (5) ∂ρN∂t+∂∂xj(ρNuj)=J (6) rd=3√3Y/(4πρLN) (7) 式中:ui、uj为时均速度分量，m/s；ρv为气相密度，kg/m3；ρ为气液混合相密度，kg/m3；pv为时均压力，Pa；μ为黏度，kg/(m·s)；δij为Kronecker delta数；E为总能量，J/kg；keff为有效导热系数，W/(m·K)；τeff为有效应力张量，无量纲；Y为液相质量分数，无量纲；rd为液滴半径，m；drd/dt为液滴生长速率，m/s；N为液滴数目，kg-1。
成核模型采用文献[15-16]中提出的双组分气体自发凝结成核模型修正方法。液滴生长过程采用Gyarmathy液滴生长模型，模型中液滴与气体间的传热系数[17-18]为
kr=λvrd11+2√8π1.5Pvγ1+γKn (8) 依据传热、传质过程，可推导得到液滴生长速率计算模型
drddt=λvρLhLV(Ts−T)(1−rcrd)rd(1+2√8π1.5Pvγ1+γKn) (9) 式中:λv为气体导热系数，W/(m·K)；Pv为气体Prandtl数；γ为气体比热比；hLV为凝结潜热，J/kg；Ts为气体压力对应的饱和温度，K；Kn表示Kundsen数。由于双组分气体不存在压力对应的饱和温度Ts这一概念，将双组分相图中露点线类比于单组分中饱和曲线。
针对气体状态方程的选择，由于低温气体已偏离理想气体假设，本研究采用了NIST真实气体模型进行计算。
2.2 湍流模型
湍流发生时会导致流体之间相互交换动量、能量，也会造成浓度的改变。本研究建模时忽略相间速度的滑移，即液滴产生不影响湍流，因此只考虑气相的湍流方程。FLUENT中提供了以下几种湍流模型：S-A模型、标准k-ε模型、RNG k-ε模型、Realizable k-ε模型、k-ω模型以及雷诺应力模型。S-A模型主要应用流动分离区附近模拟，标准k-ε模型、RNG k-ε模型一般用于各向同性的均匀湍流，k-ω模型可用于带压梯度的流动模拟和跨声速激波模拟，雷诺应力模型主要用于龙卷风、燃烧室等强烈旋转流动的模拟。对Laval喷管内跨声速流动，采用k-ω模型可以获得较为理想的计算精度和计算速度，故本研究采用该模型进行数值计算。
2.3 计算方法
气体在Laval喷管中的流动属于高速可压缩流动，采用密度基进行求解，流动控制方程组、湍流动能方程、湍流耗散率方程均采用二阶迎风格式进行离散。
根据双组分气体在Laval喷管内的高速可压缩的流动特性，入口和出口边界设置为压力入口边界和压力出口边界条件，对于气体在Laval喷管内的超声速流动，由于所有的流动参数都可从Laval喷管内部外推得到，故在出口处不进行相应设置，固体壁面边界设置为无滑移、无渗流、绝热边界条件。
在数学模型中，由于气相方程添加了源相方程，液相方程定义了标量以及引入的真实气体方程，这些仅靠在FLUENT自带的模型和材料物性无法满足要求，需要编写相应的用户自定义函数(UDF)。本研究编写UDF时，分别定义DEFINE AJUST、DEFINE SOURCE和DEFINE PROPERTY 3个宏函数。DEFINE AJUST宏用来定义过饱和度、过冷度、成核速率、液滴生长率、液滴半径、液滴质量以及液滴表面张力等参数，DEFINE SOURCE宏用来定义控制方程中的质量、动量和能量源相，DEFINE PROPERTY用来定义数值计算中用到的真实气体的热力学参数如黏度系数、导热系数等。
3. 实验验证
为验证所建立的双组分气体凝结数学模型及数值计算方法的准确性，采用本研究所设计的Laval喷管结构，在中国石油大学(华东)超声速气体凝结流动实验系统开展了水-乙醇双可凝组分气体凝结相变实验研究。实验条件为：Laval喷管入口压力0.586 MPa，入口温度288.05 K，气体湿度98.1%，水与乙醇摩尔体积比84:16，气体流量为323.78标方每小时，实验测得的Laval喷管沿程压力分布如图 2所示，可以看出，压力分布实验结果与数值计算结果吻合较好，说明本研究所建立的双组分气体超声速凝结流动特性数学模型及数值计算方法具有一定的准确性和可靠性。
图 2 Laval喷管沿程压力分布数据对比Figure 2. Comparison of pressure distribution data in Laval nozzle4. 不同压力条件下天然气超声速液化特性
保持入口温度及组成(甲烷体积分数90%、乙烷体积分数10%)不变，研究不同的入口压力对Laval喷管内部甲烷-乙烷双组分气体凝结过程中压力、温度、成核率、液滴生长率、液滴半径、液相质量分数的影响。在数值计算中设定的入口温度为270 K，设定入口压力分别为5.5、6.0和6.5 MPa。Laval喷管内双组分气体凝结参数的变化趋势及对比如图 3~图 8所示。
从压力与温度分布可以看出，气体进入Laval喷管后压力、温度不断降低，当达到一定过冷度时，气体发生凝结并释放潜热，但凝结突跃现象对压力造成的影响并不显著，压力在Laval喷管渐扩段减小到了一个比较稳定状态，温度在减小到最小值后又略微上升，这主要是由于液滴凝结释放潜热引起的。随着入口压力的增大，出口压力略微升高，温升位置有所提前，出口温度也越高，这是因为，随着入口压力的增大，液滴成核与凝结量也随之增多，释放的潜热也就越多。
从成核率分布可以看出，保持其他条件一致，当压力发生变化时，成核速率的变化趋势几乎相同，在刚进入Laval喷管的一段距离内为零，在某一位置处开始，成核率从零开始突跃一直增大到峰值后迅速减小至零。随着入口压力从5.5 MPa增大到6.5 MPa，成核的发生位置(Wilson点)不断向前移动，逐渐向Laval喷管喉部靠拢，且成核率的最值逐渐增大。当压力为5.5 MPa时，成核发生位置为x=0.120 6 m，且在x=0.147 5 m处达到极限成核，为4.044×1020 m-3·s-1；当压力为6 MPa时，成核发生位置较5.5 MPa时向前移动，为x=0.119 9 m，极限成核位置也随压力的增大而前移，在x=0.139 3 m处达到8.062×1020 m-3·s-1；当压力继续增大为6.5 MPa时，成核发生位置较6 MPa时更加靠近喉部，为x=0.118 2 m，极限成核位置在x=0.132 8 m处，为9.015×1020 m-3·s-1。
从液滴半径分布可以看出，随着入口压力的增大，Laval喷管内平均液滴半径越大，出口液滴半径也随之增大，当压力为5.5、6.0、6.5 MPa时，对应最大液滴半径尺寸分别为415.86、447.88和477.44 nm。由此可知，压力的升高有利于液滴的生长。
从液滴生长率分布可以看出，液滴生长率在气体刚进入Laval喷管时一直为零，当液滴开始发生成核凝结时液滴生长率开始突增，变化到最大值后又迅速减小，最终减小为零。综合图 6和图 7还可以看出，随着入口压力的升高，在成核开始时液滴生长率较大，液滴半径增长速度较快，但一段距离后液滴生长率下降更大，液滴半径增长速度也明显放缓。
从液相质量分数分布可以看出，伴随着混合气体的凝结成核，液相质量分数也不断增大，且随着入口压力的升高，Laval喷管出口处的湿度值随之增大，当压力为5.5 MPa时，湿度的最大值为3.989 2%，当压力增大到6.5 MPa时湿度最终增大到7.382 0%。
5. 结论
(1) 建立了三维双组分天然气混合物超声速凝结流动数学模型，对Laval喷管内双组分混合物凝结流动进行了数值模拟，得出沿Laval喷管轴向的参数分布，通过开展双可凝组分气体凝结相变实验，对比发现数值模拟与实验结果基本一致，说明了所建立数学模型及计算方法的正确性。
(2) 利用数值模型研究了入口参数对天然气混合物超声速液化特性的影响，结果表明，保持Laval喷管入口温度及组成不变，增大入口压力，混合气体成核位置前移，成核率、平均液滴半径、液相质量分数均随之增大，即入口压力越大，混合气体在Laval喷管内越易发生凝结。在实际生产中，可以通过调节入口压力来促进天然气的凝结，提高Laval喷管的液化效率。
-
图 2 水中气泡坍塌问题示意图(初始状态)
Figure 2. Air cavity collapse in water test (Description of the initial conditions)
图 3 水中气泡坍塌问题的界面函数分布图
Figure 3. Air cavity collapse in water test （Mapping of the interface function）
图 4 水中气泡塌陷问题的密度（上）与界面函数（下）分布图(
Δx=Δy=0.00625 )Figure 4. Air cavity collapse in water test (Mapping of the density (top half) and the interface function (bottom half), computed with
Δx=Δy=0.00625 )
图 5 轴线上的界面函数分布(“●” 和“□”分别表示有、无AIC)
Figure 5. Profile of the interface function along the axis y = 0 with (●) and without (□) AIC
图 6 空气-R22气柱相互作用示意图(初始状态)
Figure 6. Air-R22 shock-cylinder interaction test（Description of the initial conditions）
图 8 空气-R22激波与气柱相互作用的数值纹影图
Figure 8. Numerical schlieren diagram for the air-R22 shock-cylinder interaction problem
表 1 水中气泡塌陷问题中状态方程参数及初始参数
Table 1. Equation of state parameters and the initial time data of an air cavity collapse in water
Material ρ/(kg·m–3) p/Pa u/(m·s–1) v/(m·s–1) γ Water (Post-shock) 1.325 1.915×104 68.52 0 4.4 Water (Pre-shock) 1 1 0 0 4.4 Air 0.001 1 0 0 1.4 
下载: 导出CSV
表 2 空气-R22气柱相互作用问题中的状态方程参数及初始参数
Table 2. Equation of state parameters and the initial time data of air-R22 shock-cylinder interaction
Material ρ/(kg·m–3) p/MPa u/(m·s–1) v/(m·s–1) γ Air (Post-shock) 1.686 0.159 –113.5 0 1.400 Air (Pre-shock) 1.225 0.101 0 0 1.400 R22 3.863 0.101 0 0 1.249
下载: 导出CSV
-
[1] MESHKOV E E. Instability of the interface of two gases accelerated by a shock wave [J]. Fluid Dynamics, 1969, 4(5): 101–104. [2] BROUILLETTE M. The Richtmyer-Meshkov Instability [J]. Annual Review of Fluid Mechanics, 2002, 34(34): 445–468. [3] LOMBARDINI M, PULLIN D I, MEIRON D I. Turbulent mixing driven by spherical implosions (Part 1): flow description and mixing-layer growth [J]. Journal of Fluid Mechanics, 2014, 748(2): 85–112. [4] LOMBARDINI M, PULLIN D I, MEIRON D I. Turbulent mixing driven by spherical implosions (Part 2): turbulence statistics [J]. Journal of Fluid Mechanics, 2014, 748(10): 113–142. [5] CLEMENS N T, MUNGAL M G. Large-scale structure and entrainment in the supersonic mixing layer [J]. Journal of Fluid Mechanics, 1995, 284(284): 171–216. [6] KAWAI S, LELE S K. Large-eddy simulation of jet mixing in supersonic crossflows [J]. American Institute of Aeronautics and Astronautics, 2010, 48(9): 2063–2083. doi: 10.2514/1.J050282 [7] JOHNSEN E, COLONIUS T. Shock-induced collapse of a gas bubble in shockwave lithotripsy [J]. The Journal of the Acoustical Society of America, 2008, 124(4): 2011–2020. doi: 10.1121/1.2973229 [8] KLASEBOER E, HUNG K C, WANG C, et al. Experimental and numerical investigation of the dynamics of an underwater explosion bubble near a resilient/rigid structure [J]. Journal of Fluid Mechanics, 2005, 537(537): 387–413. [9] RANJAN D, OAKLEY J, BONAZZA R. 3D shock-bubble interactions [J]. Physics of Fluids, 2013, 25(9): 117–140. [10] THEOFANOUS T G. Aerobreakup of newtonian and viscoelastic liquids [J]. Annual Review of Fluid Mechanics, 2011, 43(1): 661–690. doi: 10.1146/annurev-fluid-122109-160638 [11] COLELLA P, WOODWARD P R. The piecewise parabolic method (PPM) for gas-dynamical simulations [J]. Journal of Computational Physics, 1984, 54(1): 174–201. doi: 10.1016/0021-9991(84)90143-8 [12] 马东军, 孙德军, 尹协远. 高密度比多介质可压缩流动的PPM方法 [J]. 计算物理, 2001, 18(6): 517–522. doi: 10.3969/j.issn.1001-246X.2001.06.008MA D J, SUN D J, YIN X Y. Piecewise parabolic method for compressible flows of multifluids with high density ratios [J]. Chinese Journal of Computational Physics, 2001, 18(6): 517–522. doi: 10.3969/j.issn.1001-246X.2001.06.008 [13] BAI J S, WANG B, WANG T, et al. Numerical simulation of the Richtmyer-Meshkov instability in initially nonuniform flows and mixing with reshock [J]. Physical Review E, 2012, 86(6): 066319. doi: 10.1103/PhysRevE.86.066319 [14] BAI J S, ZOU L Y, WANG T, et al. Experimental and numerical study of shock-accelerated elliptic heavy gas cylinders [J]. Physical Review E, 2011, 82(2): 056318. [15] XIAO J X, BAI J S, WANG T. Numerical study of initial perturbation effects on Richtmyer-Meshkov instability in nonuniform flows [J]. Physical Review E, 2016, 94(1): 013112. doi: 10.1103/PhysRevE.94.013112 [16] SHYUE K M. A fluid-mixture type algorithm for compressible multicomponent flow with van der Waals equation of state [J]. Journal of Computational Physics, 1999, 156(1): 43–88. doi: 10.1006/jcph.1999.6349 [17] SHYUE K M. An efficient shock-capturing algorithm for compressible multicomponent problems [J]. Journal of Computational Physics, 1998, 142(1): 208–242. doi: 10.1006/jcph.1998.5930 [18] ALLAIRE G, CLERC S, KOKH S. A five-equation model for the simulation of interfaces between compressible fluids [J]. Journal of Computational Physics, 2002, 181(2): 577–616. doi: 10.1006/jcph.2002.7143 [19] JOHNSEN E, COLONIUS T. Implementation of WENO schemes in compressible multicomponent flow problems [J]. Journal of Computational Physics, 2006, 219(2): 715–732. doi: 10.1016/j.jcp.2006.04.018 [20] KOKH S, ALLAIRE G. Numerical simulation of 2-D two-phase flows with interface [C]//TORO E F. Godunov Methods.Boston, MA: Springer, 2001: 513–518. [21] II S, XIE B, XIAO F. An interface capturing method with a continuous function: the THINC method on unstructured triangular and tetrahedral meshes [J]. Journal of Computational Physics, 2014, 259: 260–269. [22] SHYUE K M, XIAO F. An Eulerian interface sharpening algorithm for compressible two-phase flow: the algebraic THINC approach [J]. Journal of Computational Physics, 2014, 268(2): 326–354. [23] XIAO F, HONMA Y, KONO T. A simple algebraic interface capturing scheme using hyperbolic tangent function [J]. International Journal for Numerical Methods in Fluids, 2005, 48(9): 1023–1040. doi: 10.1002/fld.975 [24] XIAO F, II S, CHEN C. Revisit to the THINC scheme: a simple algebraic VOF algorithm [J]. Journal of Computational Physics, 2011, 230(19): 7086–7092. doi: 10.1016/j.jcp.2011.06.012 [25] KOKH S, LAGOUTIÈRE F. An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model [J]. Journal of Computational Physics, 2010, 229(8): 2773–2809. doi: 10.1016/j.jcp.2009.12.003 [26] FRIESS M B, KOKH S. Simulation of sharp interface multi-material flows involving an arbitrary number of components through an extended five-equation model [J]. Journal of Computational Physics, 2014, 273(273): 488–519. [27] DELIGANT M, SPECKLIN M, KHELLADI S. A naturally anti-diffusive compressible two phases Kapila model with boundedness preservation coupled to a high order finite volume solver [J]. Computers and Fluids, 2015, 114(1): 265–273. [28] OLSSON E, KREISS G, ZAHEDI S. A conservative level set method for two phase flow II [J]. Journal of Computational Physics, 2007, 225(1): 785–807. doi: 10.1016/j.jcp.2006.12.027 [29] SHUKLA R K, PANTANO C, FREUND J B. An interface capturing method for the simulation of multi-phase compressible flows [J]. Journal of Computational Physics, 2010, 229(19): 7411–7439. doi: 10.1016/j.jcp.2010.06.025 [30] SHUKLA R K. Nonlinear preconditioning for efficient and accurate interface capturing in simulation of multicomponent compressible flows [J]. Journal of Computational Physics, 2014, 276(1): 508–540. [31] FREUND J B, SHUKLA R K, EVAN A P. Shock-induced bubble jetting into a viscous fluid with application to tissue injury in shock-wave lithotripsy [J]. The Journal of the Acoustical Society of America, 2009, 126(5): 2746–2756. doi: 10.1121/1.3224830 [32] SHU C W. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws [C]//QUARTERONI A. Advanced Numerical Approximation of Nonlinear Hyperbolic Equations: Lecture Notes in Mathematics （Vol.1697）. Berlin, Heidelberg: Springer, 1998: 325–432. [33] NOURGALIEV R R, DINH T N, THEOFANOUS T G. Adaptive characteristics-based matching for compressible multifluid dynamics [J]. Journal of Computational Physics, 2006, 213(2): 500–529. doi: 10.1016/j.jcp.2005.08.028 [34] HU X Y, KHOO B C. An interface interaction method for compressible multifluids [J]. Journal of Computational Physics, 2004, 198(1): 35–64. doi: 10.1016/j.jcp.2003.12.018 [35] SHYUE K M. A wave-propagation based volume tracking method for compressible multicomponent flow in two space dimensions [J]. Journal of Computational Physics, 2006, 215(1): 219–244. doi: 10.1016/j.jcp.2005.10.030 [36] HAAS J F, STURTEVANT B. Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities [J]. Journal of Fluid Mechanics, 1987, 181(1): 41–76. -
下载:
下载:
点击查看大图
计量
- 文章访问数: 10976
- HTML全文浏览量: 8060
- PDF下载量: 25
