基于PPM的界面压缩方法研究

陈芳 李平 刘坤 柏劲松 林健宇 季路成

引用本文:
Citation:

基于PPM的界面压缩方法研究

    作者简介: 陈 芳(1990-),女,博士研究生,主要从事可压缩多相流模拟研究.E-mail:cocochen1929@163.com;
    通讯作者: 李平, lp0703@263.net
  • 中图分类号: O359.1

Interface Compression Technique in PPM

    Corresponding author: LI Ping, lp0703@263.net ;
  • CLC number: O359.1

  • 摘要: 高精度多组分分段抛物线法(Piecewise Parabolic Method,PPM)在对可压缩多相流问题进行模拟计算时,在不同组分交界面上存在界面扩散。为此,通过引入包含界面压缩和密度修正的人工界面压缩方法,抑制界面扩散现象。采用一个界面函数表示运动的物质界面,在多组分质量守恒方程和输运方程中添加考虑人工压缩和人工黏性的压缩源项,并在伪时间内采用二阶中心差分法和两步Runge-Kutta方法进行离散求解,采用Strang型分裂格式实现了整体算法的时间二阶精度。一维与二维数值模拟试验表明,结合人工界面压缩之后的PPM能有效抑制界面上数值扩散问题,在长时间的数值模拟中,人工界面压缩能够将扩散界面厚度维持在一定网格之内且保持界面形状不改变,尤其对于涉及稀疏波的问题,如激波引起的水中气泡坍塌,界面压缩效果更为显著。
  • 图 1  一维纯对流问题

    Figure 1.  Solution of the one-dimensional advection test

    图 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  水中气泡塌陷问题的密度(上)与界面函数(下)分布图(${\Delta x =\Delta y = 0.006\;25}$)

    Figure 4.  Air cavity collapse in water test (Mapping of the density (top half) and the interface function (bottom half), computed with ${\Delta x =\Delta y = 0.006\;25}$)

    图 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)

    图 7  不同时刻的界面分布

    Figure 7.  Mapping of the interface function at different time

    图 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/Pau/(m·s–1)v/(m·s–1)γ
    Water (Post-shock)1.3251.915×10468.5204.4
    Water (Pre-shock)11004.4
    Air0.0011001.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/MPau/(m·s–1)v/(m·s–1)γ
    Air (Post-shock)1.6860.159–113.501.400
    Air (Pre-shock)1.2250.101001.400
    R223.8630.101001.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.008
    MA 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.
  • [1] 柏劲松陈森华 . 一种捕捉可压缩流体多重交汇界面的改进型LS方法. 高压物理学报, 2003, 17(4): 268-274 . doi: 10.11858/gywlxb.2003.04.005
    [2] 柏劲松陈森华 . 重新初始化的LS方法跟踪二维可压缩多介质流界面运动. 高压物理学报, 2003, 17(1): 22-28 . doi: 10.11858/gywlxb.2003.01.004
    [3] 赵海涛王成宁建国 . 可压缩多介质问题的高精度数值模拟. 高压物理学报, 2013, 27(2): 261-267. doi: 10.11858/gywlxb.2013.02.014
    [4] 柏劲松陈森华钟敏 . 可压缩密实介质多流体高精度欧拉算法. 高压物理学报, 2002, 16(3): 204-212 . doi: 10.11858/gywlxb.2002.03.008
    [5] 赵海波肖波柏劲松段书超王刚华阚明先陈芳 . 拉氏方法模拟二维多介质可压缩流体的运动. 高压物理学报, 2018, 32(4): 042303-1-042303-13. doi: 10.11858/gywlxb.20170694
    [6] 周南 . 脉冲辐照热激波与压缩应力波. 高压物理学报, 1994, 8(3): 190-199 . doi: 10.11858/gywlxb.1994.03.006
    [7] 汤文辉张若棋经福谦胡金彪 . 冲击压缩下金属/窗口界面热弛豫过程的热阻模型研究. 高压物理学报, 1993, 7(4): 247-253 . doi: 10.11858/gywlxb.1993.04.002
    [8] 蔡灵仓刘福生经福谦 . 固态氦的理论等温压缩线. 高压物理学报, 1998, 12(2): 145-149 . doi: 10.11858/gywlxb.1998.02.013
    [9] 莫建军孙承纬 . 200 GPa压力范围内铝和铜的等熵压缩线计算. 高压物理学报, 2006, 20(4): 386-390 . doi: 10.11858/gywlxb.2006.04.008
    [10] 王桂吉谭福利孙承纬赵剑衡王刚华莫建军张宁汪小松吴刚韩梅 . 40GPa压力范围内铜和铝的准等熵压缩线. 高压物理学报, 2009, 23(4): 266-270 . doi: 10.11858/gywlxb.2009.04.005
    [11] 邹立勇刘仓理庞勇罗喜胜柏劲松杨基明 . 激波作用下SF6气泡界面演化和射流发展的数值模拟. 高压物理学报, 2013, 27(1): 90-98. doi: 10.11858/gywlxb.2013.01.013
    [12] 翟志刚董平罗喜胜 . 激波冲击V形界面的RM不稳定性实验研究. 高压物理学报, 2017, 31(6): 718-726. doi: 10.11858/gywlxb.2017.06.006
    [13] 邹立勇刘金宏谭多望黄文斌柏劲松刘仓理 . 弱激波冲击无膜重气柱和气帘界面的实验研究. 高压物理学报, 2010, 24(4): 241-247 . doi: 10.11858/gywlxb.2010.04.001
    [14] 张嘉炜黄生洪 . 极端冲击下激波诱导附加电场加速金属/气体界面的经验模型. 高压物理学报, 2019, 33(1): 012301-1-012301-9. doi: 10.11858/gywlxb.20180607
    [15] 王刚华柏劲松孙承纬莫建军王桂吉赵剑衡谭福利胡熙静 . 准等熵压缩流场反演技术研究. 高压物理学报, 2008, 22(2): 149-152 . doi: 10.11858/gywlxb.2008.02.007
    [16] 林华令 . 有限元法模拟混合物的冲击压缩特性. 高压物理学报, 1998, 12(1): 40-46 . doi: 10.11858/gywlxb.1998.01.007
    [17] 邹立勇刘金宏王建黄文斌谭多望刘仓理 . 一个用于界面不稳定性实验研究的竖式激波管. 高压物理学报, 2008, 22(2): 197-202 . doi: 10.11858/gywlxb.2008.02.014
    [18] 李巧燕施尚春杨金文孙悦 . 钕铁硼的冲击压缩特性. 高压物理学报, 2007, 21(2): 210-214 . doi: 10.11858/gywlxb.2007.02.016
    [19] 倪小军马宏昊沈兆武蒋耀港李磊 . 流-固界面冲击波入射压力的测试及数值模拟. 高压物理学报, 2013, 27(4): 528-534. doi: 10.11858/gywlxb.2013.04.010
    [20] 肖佳欣柏劲松王涛 . 密度非均匀流场中冲击加载双模态界面失稳现象的数值模拟. 高压物理学报, 2018, 32(1): 012301-1-012301-9. doi: 10.11858/gywlxb.20170501
  • 加载中
图(8)表(2)
计量
  • 文章访问数:  346
  • 阅读全文浏览量:  318
  • PDF下载量:  6
出版历程
  • 收稿日期:  2018-10-18
  • 录用日期:  2018-11-16
  • 网络出版日期:  2019-09-18
  • 刊出日期:  2019-10-01

基于PPM的界面压缩方法研究

    作者简介:陈 芳(1990-),女,博士研究生,主要从事可压缩多相流模拟研究.E-mail:cocochen1929@163.com
    通讯作者: 李平, lp0703@263.net
  • 1. 中国工程物理研究院流体物理研究所,四川 绵阳 621999
  • 2. 北京理工大学宇航学院,北京 100081

摘要: 高精度多组分分段抛物线法(Piecewise Parabolic Method,PPM)在对可压缩多相流问题进行模拟计算时,在不同组分交界面上存在界面扩散。为此,通过引入包含界面压缩和密度修正的人工界面压缩方法,抑制界面扩散现象。采用一个界面函数表示运动的物质界面,在多组分质量守恒方程和输运方程中添加考虑人工压缩和人工黏性的压缩源项,并在伪时间内采用二阶中心差分法和两步Runge-Kutta方法进行离散求解,采用Strang型分裂格式实现了整体算法的时间二阶精度。一维与二维数值模拟试验表明,结合人工界面压缩之后的PPM能有效抑制界面上数值扩散问题,在长时间的数值模拟中,人工界面压缩能够将扩散界面厚度维持在一定网格之内且保持界面形状不改变,尤其对于涉及稀疏波的问题,如激波引起的水中气泡坍塌,界面压缩效果更为显著。

English Abstract

  • 可压缩多相(多介质)流动中激波与界面的相互作用通常会产生许多复杂且重要的现象,例如:Richtmyer-Meshkov不稳定性现象[1-2]、流动混合现象[3-4]等。这些现象具有广泛的应用,如:超声速混合[5]、超声速燃烧[6]、激波碎石技术[7]、水下爆炸[8]、激波与气泡相互作用[9]以及激波与液滴相互作用[10]等。然而,人们对于涉及高密度比和强冲击波的激波与界面相互作用的现象并未完全理解。现有的针对可压缩多相流动的数值模拟方法也无法完全准确地解释以上现象,因此亟待提高数值方法的计算精度以研究激波与界面的相互作用。

    现有的模拟可压缩多相流动的高精度数值方法有很多,其中,由Colella和Woodward提出的分段抛物线方法(Piecewise Parabolic Method,PPM)是一种鲁棒性很好的间断捕捉方法。该方法采用分段抛物插值,在区域边界上采用双波近似黎曼求解器,能够有效地处理可压缩流动中的间断问题,如激波与接触间断[11]。在光滑区域,PPM具有空间三阶精度和时间二阶精度。近年来,人们针对PPM进行了一系列的研究,主要包括数值算法[12]和实际应用[13-15]两方面。算法方面,在PPM中引入了5方程模型[16],使其能够模拟具有刚性气体状态方程的多组分流动(Multi-Fluid Flows),随后发展的两步Lagrange-Remap型多组分PPM(Multi-Fluid PPM,MFPPM)[11]可以计算具有一般形式状态方程的可压缩多流体问题[12]。柏劲松等[13-15]基于MFPPM,采用大涡模拟研究了界面不稳定性问题和湍流混合等问题。但在这些研究中可以发现,由于数值扩散的存在,不同流体介质间的界面会随着计算时间的增加而变宽,降低了MFPPM的计算精度。此外,采用MFPPM计算一些实际不可混溶问题时数值扩散非常大,以致丧失其准确性。因此,在准确模拟可压缩流动界面问题中,如何抑制界面的数值扩散是研究的一个重点,同时这也是扩散界面方法研究的热点之一。

    扩散界面方法是现有的研究可压缩多相流动问题的一类重要研究方法。该方法将流体间的界面视为具有一定厚度的数值混合区域。如何利用合适的热力学状态方程和界面间断条件正确表达混合区域的特征是研究的一个难点。扩散界面方法常采用质量分数[17]、体积分数[18]或比热[19]捕捉界面。在此类方法中界面的扩散问题可以通过界面的锐化处理进行控制[20]。目前主要有:THINC方法[21-24]、反扩散方法[25-27]和界面压缩方法[28-30]。其中,Shukla等[29-30]在考虑人工压缩和人工黏性后,提出了包含界面压缩和密度修正的人工界面压缩(Artificial Interface Compression,AIC)方法,能够在长时间的数值计算中保持界面宽度在一定计算网格厚度内,同时维持界面形状不变。该方法用界面函数捕捉运动的物质界面,并限制了界面压缩的范围,能够保证即使发生界面的合并或破碎等拓扑变化时也不会出现局部极值,能够准确模拟具有大变形的可压缩多相流动问题。

    本研究在原MFPPM程序中引入AIC以实现界面压缩,并由此模拟计算一维和二维的可压缩两流体流动问题,验证方法的有效性和准确性。

    • 对于两相流动问题,定义单调且光滑的界面函数$\phi $来表示流体间的界面,$0 \leqslant \phi \leqslant 1$。在界面两边的流体中该函数值分别为0和1,也可以将$\phi $视为流体1的体积分数,则对于流体2有${\phi _2}{\rm{ = }}1 - \phi $。界面的运动可以用界面函数的对流方程来表示。

      控制方程采用5方程守恒模型[12]

      $ \frac{{\partial {{U}}}}{{\partial t}} + {{L}}\left( {{U}} \right) = {{S}}\left( {{U}} \right)\\,\;\; {{U}} = \left( \begin{gathered} {\rho _1}\phi \\ {\rho _2}\left( {1 - \phi } \right) \\ \rho {{u}} \\ \rho e \\ \phi \end{gathered} \right),\;\;{{L}}\left( {{U}} \right) = \left( \begin{gathered} \nabla \cdot \left( {{\rho _1}\phi {{u}}} \right) \\ \nabla \cdot \left[ {{\rho _2}\left( {1 - \phi } \right){{u}}} \right] \\ \nabla \cdot \left( {\rho {{uu}}} \right) + \nabla p \\ \nabla \cdot \left[ {\left( {\rho e +p} \right){{u}}} \right] \\ {{u}} \cdot \nabla \phi \end{gathered} \right),\;\;{{S}}\left( {{U}} \right) = \left( \begin{gathered} {\rho _1}{u_0}{{n}} \cdot \nabla \left[ {{\varepsilon _h}\left| {\nabla \phi } \right| - \phi \left( {1 - \phi } \right)} \right] \\ - {\rho _2}{u_0}{{n}} \cdot \nabla \left[ {{\varepsilon _h}\left| {\nabla \phi } \right| - \phi \left( {1 - \phi } \right)} \right] \\ 0 \\ 0 \\ {u_0} {{n}}\cdot\nabla \left[{{\varepsilon _h}\left| {\nabla \phi } \right| - \phi \left( {1 - \phi } \right)} \right] \end{gathered} \right) $

      式中:${\rho _1}$${\rho _2}$${p_1}$${p_2}$分别为两流体的密度和压力,${{u}} = {\left( {u,v,w} \right)^{\rm T}}$$x$$y$$z$ 3个方向的速度,$e$为内能,${u_0}$为常数。方程组(1)描述的是不可混溶、组分间无相互反应的无粘流动,且忽略了表面张力和其他外力的影响。混合区的密度、压力和总内能分别为

      $ \rho = \sum\limits_i {{\phi _i}{\rho _i}} ,\,\,\,\;\;p=\sum\limits_i {{p_i}{\phi _i}} ,\,\,\,\;\;\rho e =\sum\limits_i {{\phi _i}{\rho _i}{e_i}} ,\,\,\,\,\,\,\,\,\,i = 1,2 $

      式中:$\,{\rho _i}$${p_i}$${e_i}$分别表示第i个组分的密度、压力和内能。状态方程为

      $ p = {p_i}\left( {{\rho _i},{\rho _i}{e_i}} \right),\,\,\,\,\,\,\,\,\,{\rm{ }}i = 1,2 $

      当两流体采用理想气体状态方程或刚性气体状态方程时,5方程模型可简化为连续方程为${\rho _t} + \nabla \cdot \left( {\rho {{u}}} \right) = 0$的4方程模型。该4方程模型已被成功应用于分析激波与界面相互作用的研究中[192931]。值得注意的是,5方程模型不仅仅可以与本研究中提到的状态方程相结合,还能与其他更复杂甚至列表式的状态方程联立求解,如文献[29]。完整的控制方程模型由(1)式和(3)式给出。

      选择如下的界面函数形式

      $ \phi = \frac{1}{2}\left[ {1 \pm \tanh \left(\frac{d}{{2{\varepsilon _h}}}\right)} \right] $

      式中:${\varepsilon _h}$是与网格尺寸相关的一个小数,${\varepsilon _h}=\dfrac{{\Delta x}}{2} = \dfrac{{\Delta y}}{2} = \dfrac{{\Delta z}}{2}$$d$表示到初始界面中心位置的符号距离($d = x - {x_0}$$d = r - {r_0}$${x_0}$是初始界面的中心,${r_0}$为气泡的初始半径)。

    • 采用Strang型分裂格式实现整体算法的时间二阶精度,具体做法是:先将含有源项的控制方程分裂为不含源项的对流方程和含源项的压缩方程,采用界面压缩方法对含源项的压缩方程进行计算得到密度和界面函数,然后采用PPM计算不含源项的对流方程,之后交换压缩方程和对流方程的计算顺序,再进行计算,经过后一个计算得到的变量值就是所需的初始变量在下一时刻的值。

      $ {{{U}}^{n + 2}}={{Cp}}_{xy}^{\Delta t}{{p}}_{yx}^{\Delta t}{{C}}{{{U}}^n} $

      式中:${{{U}}^n}$表示在$n$时刻的原始变量值;${{C}}$表示在伪时间内的压缩处理;${{p}}_{xy}^{\Delta t}$${{p}}_{yx}^{\Delta t}$表示两个时间步$2\Delta t$内的PPM计算,只计算$\phi $${\rho _i}$

    • PPM是一种高精度计算格式,在光滑区具有空间三阶精度和时间二阶精度,主要计算步骤如下:(1)利用分段抛物插值定义密度、速度和压力等物理量;(2)利用双波近似黎曼求解器,在网格边界上采用牛顿法获得时间平均的压力和速度;(3)进行拉氏求解;(4)将计算结果映射回静止的欧拉网格中。变量的抛物插值详见文献[11]。在光滑区域内,PPM能给出准确解,但在间断上解的精度将下降。

    • 不同的压缩源项对应不同的界面函数方程,本研究采用

      $ {\phi _{{\tau _1}}} = {{n}} \cdot \nabla \left[ {{\varepsilon _h}\left| {\nabla \phi } \right| - \phi \left( {1 - \phi } \right)} \right] $

      式中:${{n}}$表示界面法向,一维时,${{n}}{\rm{ = }} \pm 1$${\tau _1} = {u_0}t$ 表示伪时间。该压缩源项并不会影响质量守恒性,且其压缩的有效性也已被证实,详见文献[29-30]。对应的密度修正方程为

      $ \frac{{\partial {\rho _i}}}{{\partial {\tau _2}}} = H\left( \phi \right){{n}} \cdot \left[ {\nabla \left( {{\varepsilon _h}{{n}} \cdot \nabla {\rho _i}} \right) - \left( {1 - 2\phi } \right)\nabla {\rho _i}} \right],\,\,\,\,\,\,\,\,\,i = 1,2 $

      采用$H\left[ \phi \right){\rm{ = }}{\left( {\phi \left( {1 - \phi } \right)} \right]^2}$函数来限制界面压缩的作用区域。由于界面函数值在计算中的大部分网格区域内为常数,只在界面附近的小部分区域发生变化,因此限定界面压缩作用区域能够有效地减少计算量。

      在界面法线方向,界面函数的梯度$\nabla \phi $变化很快,容易产生大的计算误差,因此,可以采用变化平缓的替代函数$\psi $来计算,使得在近界面附近,界面法向的计算不受界面函数梯度的影响。本研究采用$\psi = \dfrac{{{\phi ^\alpha }}}{{{\phi ^\alpha } + (1 - {\phi ^\alpha })}}$$0 < \alpha < 1$。随着指数$\alpha $值的增加,$\psi $变陡,取 $\alpha=0.1$,则 $0.1 \leqslant \psi \leqslant 0.9$

      在伪时间内采用二阶中心差分方法对界面压缩和密度修正项进行空间离散处理[29],时间推进采用两步Runge-Kutta方法[32]。迭代停止的条件是变量$\phi $$\rho $的无穷范数分别小于一个小数10–5和10–4。由于AIC采用的是二阶中心差分方法,压缩只作用在界面附近区域内,在远界面处界面函数值保持为0或1,因此AIC不会影响远界面处的变量值,从而不会影响PPM的精度。

      值得注意的是,在高密度比情况下,可能出现局部密度为负的现象,对于压缩项,当$\varepsilon _h^{} \geqslant 0.5h$h为网格尺寸)时,能够满足密度为非负。本研究在一维、二维模拟中,分别采用$\varepsilon _h^{1{\rm D}}{\rm{ = }}0.5h$$\varepsilon _h^{2{\rm D}}{\rm{ = }}0.72h$和均匀网格$\Delta x = \Delta y = \Delta z = h$。另外,$\Delta {\tau _1} \leqslant {h / 2}$$\Delta {\tau _2} \leqslant 2 \cdot {h^2}$,满足稳定性条件。

    • 本节给出了一些分别采用MFPPM以及MFPPM结合AIC的一维、二维数值试验结果。为了方便,所有算例都进行无量纲化处理,且都采用均匀结构网格。

    • 首先考虑一个一维纯对流问题,以证明AIC在MFPPM中的可行性。两流体均采用刚性气体状态方程${p_i} = \left( {{\gamma _i} - 1} \right){\rho _i}{e_i} + p_i^\infty $,其中:${\gamma _1} = 4.4$${\gamma _2} = 1.4$$p_1^\infty = 6 \times {10^8}$$p_2^\infty = 0$。描述的是一维的水泡(流体1)被空气(流体2)包围,压力均匀分布,速度连续。计算区域长度为1,计算网格数为100,采用周期性边界条件,初始条件为

      $\left( {\rho,u,p} \right) = \left\{ \begin{aligned} & (1,100,1)\;\;\;\quad 0.3 \leqslant x \leqslant 0.7 \\ & (0.1,100,1)\;\;\quad {\rm otherwise} \end{aligned} \right.$

      $ \phi = \frac{1}{2}\left[ {1 + \tanh \left( {\frac{{x - 0.3}}{{2\varepsilon _h^{1{\rm D}}}}} \right)} \right],\,\,\,\,\,\,\,\,\,{\rm{ }}0 \leqslant x \leqslant 0.5 $

      $ \phi = \frac{1}{2}\left[ {1 - \tanh \left( {\frac{{x - 0.7}}{{2\varepsilon _h^{1{\rm D}}}}} \right)} \right],\,\,\,\,\,\,\,\,\,{\rm{ }}0.5 < x \leqslant 1 $

      图1给出了$t = {\rm{0}}{\rm{.12}}$时刻的计算结果,其中,“□”和“●”分别代表原MFPPM的计算结果(Non-AIC)和MFPPM结合AIC后的计算结果(AIC)。从图1中可以看出,两者都与精确解(黑色直线)符合得较好,但AIC的结果更接近精确解。该算例证明了控制方程以及数值方法的正确性和有效性。

      图  1  一维纯对流问题

      Figure 1.  Solution of the one-dimensional advection test

    • 现考虑马赫数为1.72的激波诱导水中气泡塌陷问题,该问题具有广泛的物理研究和实际应用价值,如激波碎石、爆炸中热点的形成研究等[33-34]。模型示意图如图2所示,初始界面分布为

      图  2  水中气泡坍塌问题示意图(初始状态)

      Figure 2.  Air cavity collapse in water test (Description of the initial conditions)

      $ \phi = \frac{1}{2}\left[ {1 + \tanh \left( {\frac{{r - 1}}{{2\varepsilon _h^{2{\rm D}}}}} \right)} \right] $

      在气泡中,$\phi = 0$;在水中,$\phi = 1$。状态方程同上,其参数见表1,气泡初始半径为${r_0} = 1$

      Materialρ/(kg·m–3)p/Pau/(m·s–1)v/(m·s–1)γ
      Water (Post-shock)1.3251.915×10468.5204.4
      Water (Pre-shock)11004.4
      Air0.0011001.4

      表 1  水中气泡塌陷问题中状态方程参数及初始参数

      Table 1.  Equation of state parameters and the initial time data of an air cavity collapse in water

      图3给出了不同时刻的界面函数分布,从左至右分别对应0.010、0.015、0.018和0.020时刻。图3中,上、下图分别是有、无经过界面压缩的结果。每幅图片的上、下半部分分别表示$800 \times 400$$1600 \times 800$ 均匀计算网格的结果。从图3上、下两幅图的对比可以看出,界面的压缩效果很明显,尤其在上游界面有稀疏波的地方,在 $y{\rm{ = }}0$ 的轴线上压缩最为显著。MFPPM结合AIC计算得到的计算网格数为 $1600 \times 800$ 的结果如图4所示。计算结果显示无非物理振荡,且结果关于$y = 0$ 对称,与文献[3035]符合较好,验证了AIC方法的正确性和有效性。AIC的界面压缩特征使得界面形状更尖锐。图3图4中的黑实线表示气泡的初始位置。

      图  3  水中气泡坍塌问题的界面函数分布图

      Figure 3.  Air cavity collapse in water test (Mapping of the interface function)

      图  4  水中气泡塌陷问题的密度(上)与界面函数(下)分布图(${\Delta x =\Delta y = 0.006\;25}$)

      Figure 4.  Air cavity collapse in water test (Mapping of the density (top half) and the interface function (bottom half), computed with ${\Delta x =\Delta y = 0.006\;25}$)

      $t = 0.015$时刻,界面函数的等值线($\phi {\rm{ = }}0.05,0.5,0.95$)以及沿$y = 0$的剖面如图5所示。其中,左图的计算网格数是$200 \times 100$,右图是$1600 \times 800$。可以看出,经过AIC后的界面(蓝色)基本捕获在5个网格内。在远界面处MFPPM的计算结果与MFPPM结合AIC的计算结果一致。这说明界面压缩只在接触间断上有效,在激波和稀疏波附近,退化为MFPPM。界面函数的分布形状取决于混合区的密度分布,界面的数值宽度由定义的界面函数确定,这也解释了为何会出现图5中的分布形态和界面宽度。此外还可以看出,粗网格有AIC的界面比细网格无AIC的界面还要略窄。显而易见,在MFPPM中引入AIC可在整个计算过程中有效地抑制界面的扩散,且保持界面形状不变。

      图  5  轴线上的界面函数分布(“●” 和“□”分别表示有、无AIC)

      Figure 5.  Profile of the interface function along the axis y = 0 with (●) and without (□) AIC

    • 该算例是为了模拟计算Haas等[36]的一个实验,关于本算例的研究也有很多,详见文献[35]。模型示意图如图6所示,入射激波马赫数为1.22,初始界面分布为

      图  6  空气-R22气柱相互作用示意图(初始状态)

      Figure 6.  Air-R22 shock-cylinder interaction test(Description of the initial conditions)

      $ \phi = \frac{1}{2}\left[ {1 + \tanh \left( {\frac{{r - 25}}{{2\varepsilon _h^{2{\rm D}}}}} \right)} \right] $

      在R22气柱中,$\phi = 0$;在空气中,$\phi = 1$。采用理想气体状态方程,初始时刻的参数见表2,R22气柱初始半径为${r_0} = 0.025$。流场关于流动方向是对称的,因此只模拟了上半部分流动区域。

      Materialρ/(kg·m–3)p/MPau/(m·s–1)v/(m·s–1)γ
      Air (Post-shock)1.6860.159–113.501.400
      Air (Pre-shock)1.2250.101001.400
      R223.8630.101001.249

      表 2  空气-R22气柱相互作用问题中的状态方程参数及初始参数

      Table 2.  Equation of state parameters and the initial time data of air-R22 shock-cylinder interaction

      上下边界采用反射性边界条件。图7显示了不同时刻的R22气柱,其中实线为初始时刻气柱的位置,网格尺寸为${{\Delta x}/ D} = {{\Delta y} / D} = 0.002\;97$$\Delta x$$\Delta y$$D$分别表示xy方向的网格大小及圆柱初始直径。每幅图的上半部分是有AIC的结果,下半部分是MFPPM的计算结果。为了看得更清晰,这里只展示了气柱附近的局部情况。

      图  7  不同时刻的界面分布

      Figure 7.  Mapping of the interface function at different time

      图8给出了不同时刻的数值纹影图。每幅图的上半部分为有AIC的结果,下半部分为MFPPM的计算结果,红色实线表示初始时刻气柱的位置。从图8中可以看到入射激波、透射激波和反射激波。波的相互作用问题已在前人的研究中彻底进行了分析[33, 36]。采用MFPPM计算得到的结果与文献[35]中结果符合较好,包括大尺寸结构和界面的演化。另外,两个模拟结果的气柱位置与结构基本相同,在原MFPPM程序的模拟结果中出现了涡对结构扩散,而采用AIC后,在整个计算过程中,即使在高速射流形成后,界面依然保持其厚度。

      图  8  空气-R22激波与气柱相互作用的数值纹影图

      Figure 8.  Numerical schlieren diagram for the air-R22 shock-cylinder interaction problem

      经AIC处理后得到的不同涡对形状以及更好的不规则结构与实验[36]以及文献模拟结果[35]吻合较好,此处不再列出。当激波经过气柱上游界面后,会首先在气柱迎风方向的上游界面附近发展界面不稳定性,当激波经过气柱之后,界面上会出现卷起。激波经过之后界面的演化和发展情况表明,Kelvin-Helmholtz不稳定性引起了界面的变形、翻转和破碎,最后形成湍流混合。

    • 为控制计算可压缩两相流时界面上的数值扩散,在PPM中引入了人工界面压缩方法。采用5方程准守恒模型和有压缩源项的人工界面压缩(AIC)实现了对扩散界面的有效处理,采用Strang型算子分裂技术实现了整体格式的时间二阶精度。远界面上,界面函数值为常数(0或1),因此空间精度可以保持为原PPM的三阶精度。在PPM中引入AIC适用于高密度比、强冲击波情形。数值模拟试验表明该方法能够有效地处理间断,尤其是接触间断。PPM结合AIC能够顺利捕捉不同介质间的界面,在长时间的数值模拟中,能够控制界面厚度在一定数值网格范围之内,并保持界面形状不发生改变。总的来说,在PPM中引入AIC处理,有效地解决了计算中的界面扩散问题,为准确模拟可压缩两相流动提供了一种思路。其他应用方面以及涉及三维情形的研究还需投入更多精力。

参考文献 (36)

目录

    /

    返回文章
    返回