• ISSN 1001-1455  CN 51-1148/O3
  • EI Compendex、CA收录
  • 力学类中文核心期刊
  • 中国科技核心期刊、CSCD统计源期刊
高级检索 E-mail Alert

球坐标系下多介质混合物模型的数值模拟

吴宗铎 赵勇 严谨 宗智 高云

引用本文:
Citation:

球坐标系下多介质混合物模型的数值模拟

    作者简介: 吴宗铎(1984- ),男,博士,讲师,wuzongduo0@aliyun.com;
    通讯作者: 赵勇, fluid@126.com
  • 中图分类号: O.382; O359+.1

Numerical simulation about the multi-component mixture model under spherical coordinate system

    Corresponding author: ZHAO Yong, fluid@126.com ;
  • CLC number: O.382; O359+.1

  • 摘要: 利用多介质混合模型在求解球坐标系下的Riemann问题时,需要考虑界面处压力平衡性弱、奇点处理、状态方程复杂等多个难点。本文将原始基于体积分数的Mie-Grüneisen多介质混合模型扩展到球坐标系下,并对多个细节进行了修正和改进,包括:在界面处对热力学参数进行修正、采用质量分数导出新输送方程、利用质量分数加权计算偏导数、采用相邻网格点的物理量定义奇点等。经过改进后的计算模型,可以得到无振荡的数值解,而且可以准确捕捉到冲击波和界面的位置。另外,使用改进后的质量分数模型比原始的体积分数模型得到的计算结果更准确。
  • 图 1  平面冲击波界面处p的数值特征

    Figure 1.  The numerical property of p at the interfacefor plane shock wave

    图 2  冲击波与界面位置随时间的变化

    Figure 2.  Time evolutions of position of shock wave and interface

    图 3  冲击波到达3倍R0时,压力与速度的分布

    Figure 3.  The distributions of pressure and velocity when shock wave reaches 3R0

    图 4  冲击波到达不同位置处的压力曲线

    Figure 4.  Pressure curves at the time instants when the shock wave reaches different positions

    图 5  冲击波在不同位置处的压力峰值变化

    Figure 5.  Peak value variation of pressure of the shock wave at different positions

    图 6  不同计算模型下压力随时间变化曲线

    Figure 6.  Time evolutions of pressure for different numerical models

    表 1  状态方程参数

    Table 1.  Coefficients of EOS

    材料 ρ0/(kg·m−3) A/GPa B/GPa R1 R2 ω
    TNT炸药 1 630 371.2 3.21 4.15 0.95 0.35
    下载: 导出CSV
    材料 ρ0/(kg·m−3) a1/GPa a2/GPa a3/GPa b0 b1 b2 b3
    1 000 2.19 9.224 8.767 0.394 1.393 7 0 0
    下载: 导出CSV
  • [1] LIUT G, KHOOB C, YEOK S. Ghost fluid method for strong shock impacting on material interface [J]. Journal of Computational Physics, 2003, 190(2): 651–681. DOI: 10.1016/S0021-9991(03)00301-2.
    [2] SHYUE K M. A fluid-mixture type algorithm for compressible multicomponent flow with Mie-Grüneisen equation of state [J]. Journal of Computational Physics, 2001, 171(2): 678–707. DOI: 10.1006/jcph.2001.6801.
    [3] 张宝銔, 张庆民, 黄风雷. 爆轰物理学 [M]. 北京: 北京理工大学出版社, 2001: 160, 377-383.
    [4] ABGRALL R. How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach [J]. Journal of Computational Physics, 1996, 125(1): 150–160. DOI: 10.1006/jcph.1996.0085.
    [5] SAUREL R, ABGRALL R. A simple method for compressible multifuid flows [J]. SIAM Journal on Scientific Computing, 1999, 21(3): 1115–1145. DOI: 10.1137/S1064827597323749.
    [6] ABGRALL R, KARNI S. Computations of compressible multifluids [J]. Journal of Computational Physics, 2001, 169(2): 594–623. DOI: 10.1006/jcph.2000.6685.
    [7] SAUREL R, ABGRALL R. A multiphase Godunov method for compressible multifluid and multiphase flows [J]. Journal of Computational Physics, 1999, 150(2): 425–467. DOI: 10.1006/jcph.1999.6187.
    [8] 柏劲松, 陈森华, 李平. 多介质流体非守恒律欧拉方程组的数值计算方法 [J]. 爆炸与冲击, 2001, 21(4): 265–271.
    BO Jinsong, CHEN Senhua, LI Ping. Numerical methods of multicomponent flows of non-conservative Euler equations [J]. Explosion and Shock Waves, 2001, 21(4): 265–271.
    [9] 柏劲松, 陈森华, 李平, 等. 多介质可压缩流体动力学界面捕捉方法 [J]. 爆炸与冲击, 2004, 24(1): 37–43.
    BAI Jingsong, CHEN Senhua, LI Ping. Interface capturing method for compressible multi-fluid dynamics [J]. Explosion and Shock Waves, 2004, 24(1): 37–43.
    [10] 梁姗, 刘伟, 袁礼. 七方程可压缩多相流模型的HLLC格式及应用 [J]. 力学学报, 2012, 44(5): 884–895. DOI: 10.6052/0459-1879-12-022.
    LIANG Shan, LIU Wei, YUAN Li. An HLLC scheme for the seven-equation multiphase model and its application to compressible multicomponent flow [J]. Chinese Journal of Theoretical and Applied Mechanics, 2012, 44(5): 884–895. DOI: 10.6052/0459-1879-12-022.
    [11] LIANG S, LIU W, YUAN L. Solving seven-equation model for compressible two-phase flow using multiple GPUs [J]. Computers and Fluids, 2014, 99(7): 156–171. DOI: 10.1016/j.compfluid.2014.04.021.
    [12] 刘娜, 陈艺冰. 多介质流体力学计算的谱体积方法 [J]. 爆炸与冲击, 2017, 37(1): 114–119. DOI: 10.11883/1001-1455(2017)01-0114-06.
    LIU Na, CHEN Yibing. High order spectral volume method for multi-component flows [J]. Explosion and Shock Waves, 2017, 37(1): 114–119. DOI: 10.11883/1001-1455(2017)01-0114-06.
    [13] ELEUTERIO F T. Riemann solvers and numerical methods for fluid dynamics [M]. 3rd ed. Berlin Heidelberg: Springer, 2009: 1-40. DOI: 10.1007/978-3-540-49834-6.
    [14] 师华强, 宗智, 贾敬蓓. 水下爆炸冲击波的近场特性 [J]. 爆炸与冲击, 2009, 29(2): 125–130. DOI: 10.11883/1001-1455(2009)02-0125-06.
    SHI Huaqiang, ZONG Zhi, JIA Jingbei. Short-range characters of underwater blast waves [J]. Explosion and Shock Waves, 2009, 29(2): 125–130. DOI: 10.11883/1001-1455(2009)02-0125-06.
    [15] LIU T G, KHOO B C, YEO K S. The numerical simulations of explosion and implosion in air: use of a modified Harten's TVD scheme [J]. International Journal for Numerical Methods in Fluids, 1999, 31(4): 661−680. DOI: 10.1002/(SICI)1097-0363(19991030)31:4<661::AID-FLD866>3.0.CO;2-G.
    [16] LIUT G, KHOOB C, YEOK S. The simulation of compressible multi-medium flow: I: A new methodology with test application to 1D gas-gas and gas-water cases [J]. Computers and Fluids, 2001, 30(3): 291–314. DOI: 10.1016/S0045-7930(00)00022-0.
    [17] FLORES J, HOLT M. Glimm’s method applied to underwater explosion [J]. Journal of Computational Physics, 1981, 44(2): 377–387. DOI: 10.1016/0021-9991(81)90058-9.
    [18] CHIESUM J E, SHIN Y S. Explosion gas bubbles near simple boundaries [J]. Shock and Vibration, 1997, 4(1): 11–25. doi: 10.1155/1997/615415
    [19] LIU G R, LIU M B. 光滑粒子动力学:一种无网格粒子法[M]. 韩旭, 等译. 长沙: 湖南大学出版社, 2005.
    [20] ZAMYSHLYAYEV B V, YAKOVLEV Y S. Dynamic loads in underwater explosion: AD-757183[R]. USA: Naval intelligence Support Center, 1973.
  • [1] 刘娜陈艺冰 . 多介质流体力学计算的谱体积方法. 爆炸与冲击, 2017, 37(1): 114-119. doi: 10.11883/1001-1455(2017)01-0114-06
    [2] 张军赵宁任登凤谭俊杰 . Level set方法在自适应Cartesian网格上的应用. 爆炸与冲击, 2008, 28(5): 438-442. doi: 10.11883/1001-1455(2008)05-0438-05
    [3] 徐爽赵宁王春武王东红 . 水/气多介质问题的界面处理方法. 爆炸与冲击, 2015, 35(3): 326-334. doi: 10.11883/1001-1455-(2015)03-0326-09
    [4] 柏劲松陈森华李平 . 多介质流体非守恒律欧拉方程组的数值计算方法. 爆炸与冲击, 2001, 21(4): 265-271.
    [5] 柏劲松陈森华李平张展冀 . 多介质可压缩流体动力学界面捕捉方法. 爆炸与冲击, 2004, 24(1): 37-43.
    [6] 王泽平郑坚安钢恽寿榕 . 延性介质动态损伤模型. 爆炸与冲击, 1992, 12(2): 136-144.
    [7] 陈宏朱卫兵张小彬孙润鹏郭金鑫 . 真实虚拟流方法在多介质可压缩流动模拟中的应用. 爆炸与冲击, 2013, 33(1): 29-37. doi: 10.11883/1001-1455(2013)01-0029-09
    [8] 柏劲松王涛邹立勇李平 . 可压缩多介质粘性流体和湍流的大涡模拟. 爆炸与冲击, 2010, 30(3): 262-268. doi: 10.11883/1001-1455(2010)03-0262-07
    [9] 陈二云赵改平戴韧马大为 . 可压缩多介质流动的间断有限元方法. 爆炸与冲击, 2010, 30(4): 419-423. doi: 10.11883/1001-1455(2010)04-0419-05
    [10] 柏劲松李平王涛谢彬钟敏陈森华 . 可压缩多介质粘性流体的数值计算. 爆炸与冲击, 2007, 27(6): 515-521. doi: 10.11883/1001-1455(2007)06-0515-07
    [11] 张忠珍王继海 . 界面不稳定性引起湍流混合的二阶封闭模型. 爆炸与冲击, 1996, 16(3): 232-242.
    [12] 胡湘渝张德良 . 氢氧混合气体爆轰波的真实化学反应模型数值模拟. 爆炸与冲击, 2002, 22(1): 1-7.
    [13] 杨玟王丽丽周海兵张树道 . 用浮阻力模型研究Richtmyer-Meshkov不稳定性诱导混合. 爆炸与冲击, 2015, 35(3): 423-427. doi: 10.11883/1001-1455(2015)03-0423-05
    [14] 郑建国马东军孙德军尹协远 . Mie-Grneisen状态方程可压缩多流体流动的PPM方法. 爆炸与冲击, 2006, 26(2): 156-162. doi: 10.11883/1001-1455(2006)02-0156-07
    [15] 高占鹏蔡宗良张杰鲍忠兴 . 玻璃钢及松软介质的本构模型. 爆炸与冲击, 1997, 17(3): 237-248.
    [16] 袁凤英闻利群张景林 . 炸药粒度、粒度级配与质量密度的关系及球列模型研究. 爆炸与冲击, 1999, 19(4): 304-307.
    [17] 安波蒋建伟蒋浩征 . 高速钢球对水介质侵彻时瞬时空腔形成的数值模拟. 爆炸与冲击, 1998, 18(3): 245-250.
    [18] 程俊霞胡晓棉 . 轴对称坐标系下含曲率的水平集方程在非结构网格上的数值方法. 爆炸与冲击, 2012, 32(2): 150-156. doi: 10.11883/1001-1455(2012)02-0150-07
    [19] 王明洋赵跃堂邓国强钱七虎 . 爆炸波在递增硬化介质中的传播近似模型及分析解. 爆炸与冲击, 1997, 17(4): 353-358.
    [20] 欧阳春赵国志杜中华李文彬 . 弹丸垂直侵彻钢筋混凝土介质的工程解析模型. 爆炸与冲击, 2004, 24(3): 273-277.
  • 加载中
图(6)表(2)
计量
  • 文章访问数:  61
  • HTML全文浏览量:  217
  • PDF下载量:  9
  • 被引次数: 0
出版历程
  • 收稿日期:  2018-03-12
  • 录用日期:  2018-06-27
  • 刊出日期:  2019-05-01

球坐标系下多介质混合物模型的数值模拟

    作者简介:吴宗铎(1984- ),男,博士,讲师,wuzongduo0@aliyun.com
    通讯作者: 赵勇, fluid@126.com
  • 1. 广东海洋大学海洋工程学院,广东 湛江 524088
  • 2. 大连海事大学交通运输工程博士后流动站,辽宁 大连 116026
  • 3. 大连理工大学船舶工程学院,辽宁 大连 116024
  • 4. 西南石油大学油气藏地质及开发工程国际重点实验室,四川 成都 610500

摘要: 利用多介质混合模型在求解球坐标系下的Riemann问题时,需要考虑界面处压力平衡性弱、奇点处理、状态方程复杂等多个难点。本文将原始基于体积分数的Mie-Grüneisen多介质混合模型扩展到球坐标系下,并对多个细节进行了修正和改进,包括:在界面处对热力学参数进行修正、采用质量分数导出新输送方程、利用质量分数加权计算偏导数、采用相邻网格点的物理量定义奇点等。经过改进后的计算模型,可以得到无振荡的数值解,而且可以准确捕捉到冲击波和界面的位置。另外,使用改进后的质量分数模型比原始的体积分数模型得到的计算结果更准确。

English Abstract

  • 多介质的激波间断问题一直以来都是计算流体力学方面的一个备受关注的问题,目前广泛应用的一个处理方法是利用Level Set函数追踪气水界面,同时结合近似Riemann解下的虚拟流体质点完成对冲击波参数的计算[1]。虽然得到的结果较理想,但是计算过程复杂,且计算依赖近似Riemann解,无法在JWL (Jones-Wilkins-Lee)及Mie-Grüneisen等复杂形式状态方程下展开更深入的研究[2-3]

    相比之下,利用多介质混合模型,可以将流场视为一个整体并用体积分数或质量分数来区分流体,从而在不直接求近似Riemann解的前提下完成计算。多介质混合物模型早期由Abgrall提出[4-5],通过非守恒方程来计算由气体或低压缩度的固液体组成的流场[5-6]。随后,Abgrall将该混合物模型扩展成一个七方程模型,能适应一般形式状态方程的流场[7]。柏劲松等[8]也利用此类多介质混合物方程来模拟两种以上介质的冲击过程,且配合等效方程可以适应一些复杂形式的状态方程[9]。梁珊等[10]和Liang等[11]将Abgrall的七方程模型与多GPU的并行计算相结合来提高计算精确度。为了让算法更好的适应并行计算,刘娜等[12]将谱体积方法运用到多介质混合物模型中,提升了计算效率。

    虽然该多介质混合物模型在采用直角坐标系处理平面波的问题时,取得不错的进展,但是对于常用于处理球面波的球坐标系下的多介质问题,仍然具备一定的局限性。第一,多介质混合物模型建立在速度和压力在界面处维持平衡的基础,仅对平面波适用,但是球面波的压力和速度都会随着传播半径的增加而衰减;第二,球坐标系下,原点位置处为冲击波运动的奇点位置,必须作出合适的处理;第三,为了增加工程实用性,计算模型需要能适应复杂形式的状态方程,而不仅仅是理想气体方程等形式简单的状态方程。

    考虑到以上因素,本文将直角坐标系下的基于Mie-Grüneisen状态方程下的多介质混合模型延伸到球坐标系下。为了压力和速度不平衡问题以及坐标原点的奇点问题,在原计算模型下作了参数修正并用质量分数替换体积分数,并对奇点处数值差分作了特殊处理—利用相邻网格点参数赋值给奇点参数但奇点处速度设定为0。最后,将通过数值算例验证其数值稳定性和计算准确性。

    • 考虑多介质的可压缩流体组成的混合流场,其运动形式可以用欧拉方程来描述:

      $ \left\{ \begin{aligned} &\frac{{\partial \rho }}{{\partial t}} + \frac{\partial }{{\partial x}}\left( {\rho u} \right) = 0\\ &\frac{{\partial \left( {\rho u} \right)}}{{\partial t}} + \frac{\partial }{{\partial x}}\left( {\rho {u^2} + p} \right) = 0\\ &\frac{{\partial \left( {\rho E} \right)}}{{\partial t}} + \frac{\partial }{{\partial x}}\left( {\rho Eu + pu} \right) = 0 \end{aligned} \right. $

      式中:ρpuE分别为密度、压力、速度和单位质量能量[13]。单位质量的总能量为内能和动能的总和,且E=e+u2/2。

      流场中介质的物理特性,可以用Mie-Grüneisen状态方程来表述:

      $ p - {p_{{\rm{ref}}}}(\rho ) = \rho \varGamma (\rho ,e)[e - {e_{{\rm{ref}}}}(\rho )] $

      式中:Γ为Mie-Grüneisen系数,preferef为参考点的压力和单位质量内能,与密度ρ有关。

    • 界面处,流体的压力p和速度u均保持平衡[4],此时pu在界面处不发生变化,而密度ρ与内能e则出现数值间断。因此∂p/∂x和∂u/∂x约为0,式(1)中的质量守恒与能量守恒可以分别表示为:

      $ \frac{{\partial \rho }}{{\partial t}} + u\frac{{\partial \rho }}{{\partial x}} = 0 $

      $ \frac{{\partial (\rho e)}}{{\partial t}} + u\frac{{\partial (\rho e)}}{{\partial x}} = 0 $

      将式(2)代入式(4),有:

      $ \frac{\partial }{{\partial t}}\left( {\frac{{p - {p_{\rm {ref}}}}}{\varGamma } + \rho {e_{\rm {ref}}}} \right) + u\frac{\partial }{{\partial x}}\left( {\frac{{p - {p_{\rm {ref}}}}}{\varGamma } + \rho {e_{\rm {ref}}}} \right) = 0 $

      由于Γpreferef都只和密度ρ有关的函数,且式(5)对于任意的表达式pref(ρ)和eref(ρ)均成立。这里首先假设preferef都为0,那么有[2,9]

      $ \frac{\partial }{{\partial t}}\left( {\frac{1}{\varGamma }} \right) + u\frac{\partial }{{\partial x}}\left( {\frac{1}{\varGamma }} \right) + \rho \left[ {\frac{\partial }{{\partial \rho }}\left( {\frac{1}{\varGamma }} \right)} \right]\frac{{\partial u}}{{\partial x}} = 0 $

      同理,可依次类推得到关于preferef的非守恒形式的方程:

      $ \frac{\partial }{{\partial t}}\left( {\frac{{{p_{{\rm{ref}}}}}}{\varGamma }} \right) + u\frac{\partial }{{\partial x}}\left( {\frac{{{p_{{\rm{ref}}}}}}{\varGamma }} \right) + \rho \left[ {\frac{\partial }{{\partial \rho }}\left( {\frac{{{p_{\rm {ref}}}}}{\varGamma }} \right)} \right]\frac{{\partial u}}{{\partial x}} = 0 $

      $ \frac{\partial }{{\partial t}}\left( {\rho {e_{{\rm{ref}}}}} \right) + u\frac{\partial }{{\partial x}}\left( {\rho {e_{{\rm{ref}}}}} \right) + \rho \left[ {\frac{\partial }{{\partial \rho }}\left( {\rho {e_{{\rm{ref}}}}} \right)} \right]\frac{{\partial u}}{{\partial x}} = 0 $

    • 由于式(6)~(8)的推导过程需要基于pu在界面处连续这一前提条件,但是在左右波面处,pu均出现间断面,式(3)和(4)不再成立。此时,对(6)作如下考虑:

      $ \begin{split} \frac{\partial }{{\partial t}}\left( {\frac{1}{\varGamma }} \right) + u\frac{\partial }{{\partial x}}\left( {\frac{1}{\varGamma }} \right) + \rho \left[ {\frac{\partial }{{\partial \rho }}\left( {\frac{1}{\varGamma }} \right)} \right]\frac{{\partial u}}{{\partial x}} =& \left[ {\frac{\partial }{{\partial \rho }}\left( {\frac{1}{\varGamma }} \right)} \right]\frac{{\partial \rho }}{{\partial t}} + u\left[ {\frac{\partial }{{\partial \rho }}\left( {\frac{1}{\varGamma }} \right)} \right]\frac{{\partial \rho }}{{\partial x}} + \rho \left[ {\frac{\partial }{{\partial \rho }}\left( {\frac{1}{\varGamma }} \right)} \right]\frac{{\partial u}}{{\partial x}} \\ {\rm{ = }}&\left[ {\frac{\partial }{{\partial \rho }}\left( {\frac{1}{\varGamma }} \right)} \right]\left(\frac{{\partial \rho }}{{\partial t}} + \frac{{\partial \rho u}}{{\partial x}}\right) {\rm{ = }}0 \end{split} $

      因此,在左右波面处虽然pu均间断,但由于在同类介质中1/Γ等参数的导数保持连续,因此,式(6)依然成立。同理,式(7)和(8)在左右波面处也成立。

    • 对于多介质流场,必须利用输送方程来捕捉界面的运动。流场中,某种介质的界面随时间变化的运动方程,称为输送方程。一般利用体积分数zi来建立介质的输送方程:

      $ \frac{{\partial {{\textit z}_i}}}{{\partial t}} + u\frac{{\partial {{\textit z}_i}}}{{\partial x}} = 0\;\;\;\;\;\;\;i = 1,2, \cdots $

      式中:i为介质的种类序号。非守恒变量1/Γpref/Γρeref均为独立变量,但它们对的ρ偏导数可表达成关于体积分数zi的加权叠加。这样,综合式(1)、(6)~(8)和(10),可完成直角坐标系下的多介质流场的求解。

    • 球坐标系下,则需要考虑更多的问题。

      第一,在采用直角坐标系下计算平面波问题时,pu在界面处维持不变且偏导数∂p/∂x与∂u/∂x基本为0。该数值特性一直从界面处向左右延伸,直到波面处才形成间断,如图1所示。但是采用球坐标系处理球面波问题时,当冲击波沿着半径r向外扩散时,传播的形式则不再是平面冲击波,而是球面冲击波。此时,虽然pu依然在界面处连续,但是球面波的扩散过程中pu会随着半径增加而发生衰减,并非保持不变。因此,球面冲击波在界面处的偏导数∂p/∂r与∂u/∂r不再是0,而是一个稳定的负数。这种情况下,式(3)~(4)非守恒性方程的建立的理论基础并不十分准确。

      图  1  平面冲击波界面处p的数值特征

      Figure 1.  The numerical property of p at the interfacefor plane shock wave

      第二,球坐标系下的内外介质产生相互作用时,内部的稀疏波(或激波)汇聚到一起后又反向朝着外部扩散,原点处便形成了一个奇点。奇点处的稀疏波形成了独特的反射机制。因此奇点的物理参数变化需仔细考虑。

    • 这里将球坐标系下的控制方程表示成守恒变量U与通量F的关系为[13-14]:

      $ \frac{{\partial {U}}}{{\partial t}} + \frac{{\partial {F}}}{{\partial r}} = {S}({U}) $

      其中: $ {U} = \left[ {array}{l}\rho \\\rho u\\\rho E{array} \right],\;{F} = \left[ {array}{l}\rho u\\\rho {u^2} + p\\\rho Eu + pu{array} \right],\;{S} = - \displaystyle\frac{\alpha }{r}\left[ {array}{l}\rho u\\\rho {u^2} + p\\\rho Eu + pu{array} \right]$

      式中:系数α=2。式(11)可以改写成:

      $ \frac{{\partial \tilde{ U}}}{{\partial t}} + \frac{{\partial \tilde{ F}}}{{\partial r}} = {S}(\tilde{ U}) $

      其中: $ \tilde{ U} = {r^2}{U},\;\tilde{ F} = {r^2}{F},\;{S}(\tilde{ U}) = {[0,\;2p/r, \;0]^{\rm T}}$

      此时,将非守恒形式的方程组添加到式(11)中,即可得到球坐标系下的控制方程。

    • 需要注意的是,非守恒形式的方程式(6)~(8)的推导,是以界面处于平衡状态为前提的。在球坐标系下,由于冲击波是沿着球面向四周扩散。随着扩散范围的增加,几乎所有的冲击波参数都会随着半径增加而降低。此时,pu在界面处将呈现出持续的减小趋势。这样,式(6)~(8)无法直接使用到式(12)中。在直角坐标系和球坐标系下,界面间断的关系分别为

      直角坐标系 $\displaystyle\frac{{\partial p}}{{\partial x}},\;\displaystyle\frac{{\partial u}}{{\partial x}} \to {\rm{0,}}\quad\;\;\;\displaystyle\frac{{\partial \rho }}{{\partial x}}, \;\displaystyle\frac{{\partial (\rho e)}}{{\partial x}},\;\displaystyle\frac{\partial }{{\partial x}}\left( {\displaystyle\frac{1}{\varGamma }} \right),\;\displaystyle\frac{\partial }{{\partial x}}\left( {\displaystyle\frac{{{p_{{\rm{ref}}}}}}{\varGamma }} \right),\;\displaystyle\frac{\partial }{{\partial x}}\left( {\rho {e_{{\rm{ref}}}}} \right) \to \infty $

      球坐标系 $\displaystyle\frac{{\partial p}}{{\partial r}}{\rm{\text{<} 0}},\; \displaystyle\frac{{\partial u}}{{\partial r}}{\rm{ \text{<} 0,}} \quad\;\;\;\displaystyle\frac{{\partial \rho }}{{\partial r}},\displaystyle\frac{{\partial (\rho e)}}{{\partial r}},\displaystyle\frac{\partial }{{\partial r}}\left( {\displaystyle\frac{1}{\varGamma }} \right),\displaystyle\frac{\partial }{{\partial r}}\left( {\displaystyle\frac{{{p_{{\rm{ref}}}}}}{\varGamma }} \right),\displaystyle\frac{\partial }{{\partial r}}\left( {\rho {e_{{\rm{ref}}}}} \right) \to \infty $

      考虑到三个非守恒变量1/Γpref/Γρeref均为密度ρ的函数,这里每隔一个时间步对它们做修正:

      $ \frac{1}{\varGamma }{\rm{ = }}\frac{1}{{\varGamma ({\rho _i})}},\;\;\; \frac{{{p_{{\rm{ref}}}}}}{\varGamma } = \frac{{{p_{{\rm{ref}}}}({\rho _i})}}{{\varGamma ({\rho _i})}},\;\;\;\rho {e_{{\rm{ref}}}} = {\rho _i}{e_{{\rm{ref}}}}({\rho _i})\;\;\;\;\;\;\frac{{{\rho _i}}}{\rho } \text{>} 0.99\;\;\;\;\;i = 1,2, \cdots $

      式中:ρi为介质i的密度。而的其他参数如ρpuE的计算,仍然借助当前时间步下的Γpreferef。而式(6)~(8)则保留在系统方程中,并且每隔一个时间步,非守恒变量1/Γpref/Γρeref都会得到修正。当时间步足够密时,非守恒形式在球坐标系下带来的误差可以控制在一定的范围内。

    • 由于每隔一个时间步,都会需要利用ρi进行修正。这里用质量分数代替体积分数来构造输送方程。利用式(12)中的质量守恒关系可得:

      $ \frac{{\partial ({r^2}\rho {y_i})}}{{\partial t}} + \frac{{\partial ({r^2}\rho {y_i}u)}}{{\partial r}} = 0\;\;\;\;\;\;i = 1,2, \cdots $

      根据式(14)得到的质量分数yi可以很容易得到各介质的密度ρi (i=1, 2, ···)。而非守恒方程(6)~(8)中的偏导数,则转换成与质量分数相关的加权叠加:

      $ \begin{aligned} &\phi {\rm{ = }}\frac{\partial }{{\partial \rho }}\left( {\frac{1}{\varGamma }} \right){\rm{ = }}{y_1}{\left. {\left( { - \frac{{{{\varGamma }_1'}}}{{\varGamma _1^2}}} \right)} \right|_{\rho {y_1}}} + {y_2}{\left. {\left( { - \frac{{{{\varGamma }_2'}}}{{\varGamma _2^2}}} \right)} \right|_{\rho {y_2}}} + \cdots \\ &\varphi {\rm{ = }}\frac{\partial }{{\partial \rho }}\left( {\frac{{{p_{{\rm{ref}}}}}}{\varGamma }} \right){\rm{ = }}{y_1}{\left. {\left( {\frac{{{\varGamma _1}{{p}_{{\rm{ref,}}1}'} - {p_{{\rm{ref,1}}}}{{\varGamma }_1'}}}{{\varGamma _1^2}}} \right)} \right|_{\rho {y_1}}} + {y_2}{\left. {\left( {\frac{{{\varGamma _2}{{p}_{{\rm{ref}},2}'} - {p_{{\rm{ref}},2}}{{\varGamma }_2}'}}{{\varGamma _2^2}}} \right)} \right|_{\rho {y_2}}} + \cdots \\ &\psi {\rm{ = }}\frac{\partial }{{\partial \rho }}\left( {\rho {e_{{\rm{ref}}}}} \right){\rm{ = }}{y_1}{\left. {\left( {{\rho _1}{{e}_{{\rm{ref,}}1}'} + {e_{{\rm{ref,}}1}}} \right)} \right|_{\rho {y_1}}} + {y_2}{\left. {\left( {{\rho _2}{{e}_{{\rm{ref,}}2}'} + {e_{{\rm{ref,}}2}}} \right)} \right|_{\rho {y_2}}} + \cdots \end{aligned} $

      式中:ϕφψ为非守恒变量1/Γpref/Γρeref偏导数,下角标1、2代表流场中介质的序号。

      相比式(10)给出的基于体积分数的输送方程,基于质量分数的输送方程能考虑到每一种介质的密度ρi,从而得到更为合理的偏导数ϕφψ

    • 对于球坐标系来说,奇点是一个无法避免的问题。参考文献[15]中对奇点的处理方式,奇点处质点的运动速度为0,即u=0,但是奇点处的质量和能量守恒关系依然成立:

      $ \begin{aligned} \frac{{\partial \rho }}{{\partial t}} + 3\frac{{\partial (\rho u)}}{{\partial r}} = 0, \quad\;\quad\; \frac{{\partial \rho E}}{{\partial t}} + 3\frac{{\partial (\rho Eu + pu)}}{{\partial r}} = 0 \end{aligned} $

      展开后,将u=0代入,有:

      $ \begin{aligned} \frac{{\partial \rho }}{{\partial t}} + 3\rho \frac{{\partial u}}{{\partial r}} = 0, \quad\;\quad\; \frac{{\partial \rho E}}{{\partial t}} + 3(\rho E + p)\frac{{\partial u}}{{\partial r}} = 0 \end{aligned} $

      由于在奇点处主要涉及冲击波的汇聚和发散,不涉及界面的运动,因此,可不考虑输送方程。参数Γpreferef由密度ρ直接计算得到。

      文献[15]中,对偏导数∂u/∂r采用了如下的计算方式:

      $ \frac{{\partial u}}{{\partial r}} = 2\frac{{{u_1}}}{{\Delta r}} - \frac{{{u_2}}}{{2\Delta r}} $

      式中:下角标1和2代表网格点序号。但考虑到本文的方程体系非守恒性较强,将u1u2都替换为奇点处的速度0,这样偏导数∂u/∂r=0。再根据式(17)计算得到奇点处的其他参数。

    • 在有限体积差分下,第j个网格点的变量Uj,可按如下格式计算:

      $ \frac{{\tilde { U}_j^{(n + 1)} - \tilde { U}_j^{(n)}}}{{\Delta t}} + \frac{{{{\tilde { F}}_{j + 1/2}} - {{\tilde { F}}_{j - 1/2}}}}{{\Delta r}} = {\tilde { S}_j} $

      式中:Δt和Δr分别为时间和空间步长,上标(n)代表第n步时间序号。对于物理量${\tilde { U}_j}$,第jj+1个网格点之间的通量${\tilde { F}_{j + 1/2}}$和源项${\tilde { S}_j}$分别为:

      $ {{\tilde{ U}}_j} = r_j^2\left[ \begin{array}{l} {\rho _j}\\ {\rho _j}{u_j}\\ {\rho _j}{E_j}\\ 1/{\varGamma _j}\\ {p_{ref\;j}}/{\varGamma _j}\\ {\rho _j}{e_{{\rm{ref}}\;j}}\\ \rho {y_i}_j \end{array} \right]\;\;\;\; {{\tilde{ F}}_{j + 1/2}} = r_j^2\left[ \begin{array}{l} {\rho _{j + 1/2}}{u_{j + 1/2}}\\ {\rho _{{\rm{j + }}1/2}}u_{{\rm{j + }}1/2}^2 + {p_{{\rm{j + }}1/2}}\\ ({\rho _{j + 1/2}}{E_{j + 1/2}} + {p_{j + 1/2}}){u_{j + 1/2}}\\ {u_j} \cdot {(1/\varGamma )_{j + 1/2}} + {\rho _j}{\phi _j} \cdot {u_{j + 1/2}}\\ {u_j} \cdot {({p_{{\rm{ref}}}}/\varGamma )_{j + 1/2}} + {\rho _j}{\varphi _j} \cdot {u_{j + 1/2}}\\ {u_j} \cdot {({\rho}{e_{{\rm{ref}}}})_{j + 1/2}} + {\rho _j}{\psi _j} \cdot {u_{j{\rm{ + }}1/2}}\\ {(\rho {y_i})_{j + 1/2}} \cdot {u_{{\rm{j + }}1/2}} \end{array} \right]\;\;\;\; {{\tilde{ S}}_j} = 2{r_j}\left[ \begin{array}{l} 0\\ {p_j}\\ 0\\ 0\\ 0\\ 0\\ 0 \end{array} \right] $

      式中:角标j+1/2表示通量中的参数,角标j表示网格点的参数;在球面冲击波中,通量${{\tilde{ F}}_{j + 1/2}}$表示沿半径r方向通过j点的物理变量。

    • 时间步长Δt,满足收敛的条件为:

      $ \Delta t = {C_{{\rm{CFL}}}}\frac{{\min (\Delta r)}}{{\max (\left| {{u_j}} \right| + {c_j})}} $

      式中:CCFL是为了保证迎风型计算格式收敛而设置的一个0~1之间的系数,满足收敛的条件为CCFL<1, 并且当CCFL取得越小时间步长越小,计算越精细, 这里取CCFL=0.4; ujcj表示第j个网格点的质点速度和声速。声速cj的表达式为:

      $ {{c}_{j}}=\sqrt{{\left.\left( {{H}_{j}}-\frac{1}{2}{{u}_{j}}{}^{2}-{{p}_{j}}{{\phi }_{j}}+{{\varphi }_{j}}-{{\psi }_{j}} \right)\right/}{\left( \frac{1}{{{\varGamma }_{j}}} \right)}\;} $

      式中:H=E+p/ρ

    • 这里考虑一个无因次的气水作用问题[16]。初始时刻,中心气团的物理状态为:p=8.3×103ρ=1.27,u=0,外部水的物理状态为:p=1.0,ρ=1.27,u=0。气体和水的状态方程分别为

      $ \begin{array}{c} p = (\gamma - 1)\rho e,\quad\;\;\gamma = 1.4;\\ p = (\gamma - 1)\rho e - \gamma B,\quad\;\;\gamma = 7.0,\quad\;\;B = 3\;000 \end{array} $

      式中:γB为描述水和气体在低压缩度时的热力学参数。气团初始半径为R0。由于中心气团的压力和密度高于周边水的。因此,高压的气团会在水中形成冲击波并推动界面向外部扩散。对冲击波波面和介质界面的运动位置的变化情况进行记录(波面为水中压力的最大值位置,界面位置为从外向内气体的质量分数刚超过0.5的位置),并绘制了冲击波和界面随时间变化的曲线,如图2所示。从图2可以看出,本文所得到的冲击波和界面的运动位置与Liu等[16]和Flores等[17]的结果都比较接近。三个数据的计算结果,在气水作用的初期差别不大。但是随着冲击波和界面向外扩散,本文的计算结果与Flores等[17]的计算结果相比,误差略偏大,但仍然在一个合理范围内。

      图  2  冲击波与界面位置随时间的变化

      Figure 2.  Time evolutions of position of shock wave and interface

      图3为冲击波到达3倍气团初始半径时候的压力和速度曲线。此时,向内收缩的稀疏波已达到中心原点处。由于气团已扩散开,中心处的压力开始逐渐下降,但中心处的速度仍然为0。速度曲线中,由中心向外的第一个拐点为稀疏波的端点,第二个拐点则为界面所在位置。另外,从图3中可以看出,处于数值奇点的原点处,数值稳定性较好。证明本文所用的奇点处理方式能效果良好。

      图  3  冲击波到达3倍R0时,压力与速度的分布

      Figure 3.  The distributions of pressure and velocity when shock wave reaches 3R0

    • 假设一TNT的实心药球质量为50 kg,将空间步长取为炸药球半径的1/100。利用改进后的多介质混合模型计算球坐标系下的冲击波运动情况。

      炸药的状态方程为JWL状态方程[3]

      $ p = A\left( {1 - \frac{{\omega \theta }}{{{R_1}}}} \right){{\rm e}^{ - {R_1}/\theta }} + B\left( {1 - \frac{{\omega \theta }}{{{R_2}}}} \right){{\rm e}^{ - {R_2}/\theta }} + \omega \theta {\rho _0}e $

      式中:ABR1R2ω均为常数,由圆筒试验标定得到[3]θ=ρ/ρ0ρ0为物质的初始密度。水的多项式状态方程为[18]

      $ p = \left\{ {\begin{array}{*{20}{l}} {{a_1}\mu + {a_2}{\mu ^2} + {a_3}{\mu ^3} + \left( {{b_0} + {b_1}\mu + {b_2}\mu + {b_3}\mu } \right){\rho _0}E}&\quad{\mu \text{>} 0}\\ {{a_1}\mu + \left( {{b_0} + {b_1}\mu } \right){\rho _0}E}&\quad{\mu \text{<} 0} \end{array}} \right. $

      式中:μ=θ−1;a1a3, b0b3为展开后的多项式系数。状态方程(23)~(24)的参数取值见1[19]

      材料 ρ0/(kg·m−3) A/GPa B/GPa R1 R2 ω
      TNT炸药 1 630 371.2 3.21 4.15 0.95 0.35

      表 1  状态方程参数

      Table 1.  Coefficients of EOS

      材料 ρ0/(kg·m−3) a1/GPa a2/GPa a3/GPa b0 b1 b2 b3
      1 000 2.19 9.224 8.767 0.394 1.393 7 0 0

      本算例中,将计算范围扩大至20倍以上的药球半径。这样,可以得到冲击波运动到不同位置处的压力,如图4所示。从图4可以看出,由于在偏导数ϕφψ的处理方式差异,两种计算模型得到的压力分布有差别。但都具备较好的数值稳定性。当冲击波运动到较远位置时,会出现多个压力峰值。

      图  4  冲击波到达不同位置处的压力曲线

      Figure 4.  Pressure curves at the time instants when the shock wave reaches different positions

      在实际工程中,最受重视的是第一次的峰值。图5为冲击波峰值压力随运动距离变化的曲线。作为对比,将给初始的体积分数模型[2]和改进后的质量分数模型同Zamyshlyayev的经验公式[20]曲线作了对比。对比发现,改进后的质量分数模型可以得到更为准确的计算结果。在压力随运动距离发生急剧下降的过程中,利用体积分数得到的计算结果,得到的结果比经验公式的结果略微偏大。

      图  5  冲击波在不同位置处的压力峰值变化

      Figure 5.  Peak value variation of pressure of the shock wave at different positions

      图6为利用体积分数与质量分数计算模型分别得到的压力时程曲线。从不同的计算结果中可以看出,利用改进后的质量分数模型,得到的计算结果好于体积分数模型。各个半径处的压力时程,质量分数模型均吻合得比较好。由于Zamyshlyayev的经验公式仅考虑冲击波压力最大峰值未考虑压力的二次峰值,本文也仅模拟了冲击波压力的最大峰值。

      图  6  不同计算模型下压力随时间变化曲线

      Figure 6.  Time evolutions of pressure for different numerical models

    • 针对球坐标下多介质界面捕捉困难的问题,本文对原始的Mie-Grüneisen多介质混合物模型作了热力学参数修正、输送方程变换、特殊奇点处理等多项数值改进。改进后的基于质量分数的模型,可以在常规状态方程下准确捕捉界面和冲击波的运动,且可以在适应JWL等形式复杂的状态方程。而适当的奇点处理也保证了数值计算在奇点处的稳定性。另外,数值算例证明,相比原始的体积分数模型,改进后的质量分数模型可以得到较为准确的计算结果。

参考文献 (20)

目录

    /

    返回文章
    返回