圆柱形汇聚激波诱导 Richtmyer-Meshkov不稳定的 SPH 模拟1)

doi:10.6052/0459-1879-19-041

圆柱形汇聚激波诱导 Richtmyer-Meshkov不稳定的 SPH 模拟 ¹⁾

徐建于, 黄生洪^,²⁾

中国科学技术大学工程科学学院,合肥 230026

NUMERICAL SIMULATION OF CYLINDRICAL CONVERGING SHOCK INDUCED RICHTMYER-MESHKOV INSTABILITY WITH SPH ¹⁾

Xu Jianyu, Huang Shenghong^,²⁾

School of Engineering Science, University of Science and Technology of China, Hefei 230026, China

通讯作者: 2) 黄生洪,副教授,主要研究方向：先进SPH算法研究.E-mail:hshnpu@ustc.edu.cn

收稿日期: 2019-02-1 接受日期: 2019-05-21 网络出版日期: 2019-07-18

基金资助:

1) 挑战专题TZ2016001

国家自然科学基金. U1530125

Received: 2019-02-1 Accepted: 2019-05-21 Online: 2019-07-18

作者简介 About authors

摘要

汇聚激波诱导不同物质界面的Richtmyer-Meshkov(RM)不稳定现象在惯性约束核聚变领域有重要的学术意义和工程背景.基于网格离散的宏观流体力学方法由于数值扩散问题往往需要高阶精度算法才能准确追踪界面演化,且对大变形和破碎合并等复杂界面追踪也极为困难.光滑粒子流体动力学(smoothed particlehydrodynamics,SPH)方法采用纯拉格朗日算法,可以有效克服上述难点.但经典SPH算法需采用人工黏性处理强间断,在激波间断处往往会出现严重的非物理振荡,对于涉及强冲击不稳定性问题,很难达到理想的模拟效果.本文采用基于HLL黎曼求解器的SPH算法,实现了对强激波和大密度比物质界面的有效分辨和追踪.一维数值校核证明了代码的可靠性、健壮性,并进一步模拟了二维圆柱形汇聚冲击波冲击四边形轻/重气界面诱导的RM不稳定性问题,与已有实验结果进行了对比,发现模拟结果与实验结果吻合.通过分析界面演化过程中的密度及压力变化,发现本文所采用的方法可准确地追踪激波与界面作用的复杂界面和波系演化规律.研究结果为进一步理解和解释汇聚冲击条件下的RM不稳定性机理奠定了基础.

关键词： SPH ; 黎曼求解器 ; Richtmyer-Meshkov不稳定性 ; 激波 ; 圆柱形汇聚 ; 数值模拟

Abstract

The Richtmyer-Meshkov (RM) instability induced by converging shock waves at interfaces of different substances has an important academic significance and engineering background in the field of inertial confinement fusion. The macroscopic fluid dynamics method based on grid discretization requires high order precision algorithm to track the interface evolution accurately because of numerical diffusion problem, and it is extremely difficult to track the complex interface evolution such as large deformation and fragmentation merging, etc. Smoothed particle hydrodynamics (SPH) method is a pure Lagrangian algorithm, which can effectively overcome the addressed difficulties. However, the classical SPH algorithm requires artificial viscosity to smooth the strong discontinuities, otherwise large non-physical oscillations may occur. For the problem involving strong shock instability, it is difficult to achieve ideal results. In this paper, the SPH algorithm based on HLL Riemann solver is adopted to effectively distinguish and track the strong shock wave and the material interface with a large density ratio. The reliability and robustness of the code were validated by four classical 1D shock tube tests, and it is found that the smoothing effect of the density algorithm used by SPH on the contact discontinuity can be improved by reducing the initial particle spacing. The smaller the initial particle spacing, the higher the numerical simulation accuracy, but it cannot be completely eliminated. However, the position of the interface is actually marked by the media properties of the particles, and does not affect the discrimination of the interface position under the SPH Lagrangian algorithm. Then the 2D cases of RM instability induced by cylindrical converging shock wave impacting at the quadrilateral light/heavy gas interface were simulated. It is found that the simulation results are quantitatively in good agreement with the existing experimental results. By analyzing the density and pressure changes in the process of interface evolution, it is also found that the models and methods adopted can accurately track the complex interfaces and shock waves evolution patterns during the RM instability process. The relevant results lay a foundation for further understanding and explanation of RM instability mechanism under extreme converging shock conditions.

Keywords： smoothed particle hydrodynamics ; Riemann solver ; Richtmyr-Meshkov instability ; shock ; cylindrical converging ; numerical simulation

PDF (29026KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

徐建于, 黄生洪. 圆柱形汇聚激波诱导 Richtmyer-Meshkov不稳定的 SPH 模拟 ¹⁾. 力学学报[J], 2019, 51(4): 998-1011 DOI:10.6052/0459-1879-19-041

Xu Jianyu, Huang Shenghong. NUMERICAL SIMULATION OF CYLINDRICAL CONVERGING SHOCK INDUCED RICHTMYER-MESHKOV INSTABILITY WITH SPH ¹⁾. Chinese Journal of Theoretical and Applied Mechanics[J], 2019, 51(4): 998-1011 DOI:10.6052/0459-1879-19-041

引言

两种不同密度的轻重流体交界面,在受到激波冲击作用时,会发生扰动,导致界面失稳,这种现象被称为Richtmyer-Meshkov(RM)不稳定性.RM不稳定性在各种流体运动中普遍存在^[1-3],尤其是惯性约束核聚变(inertialconfinement fusion, ICF)问题中.惯性约束核聚变是由激光辐照含有聚变原料的靶丸,使表面烧蚀层物质快速变为向外喷射的等离子体,同时形成向靶丸中心汇聚的圆柱形或球形激波.当汇聚激波通过靶丸的不同物质界面时,由于加工精度等原因,界面存在一些初始微小不光滑结构(扰动),因此发生RM不稳定性现象.这可能会使反应温度降低,甚至中止反应.所以,研究ICF极端条件下多介质界面不稳定性是实现聚变点火的关键基础课题之一,具有重要的工程意义和学术价值.

RM不稳定性最初由Richtmyer^[4]于1960年提出,并进行了严格的理论推导,得到首个RM不稳定性的线性冲击模型,后来Meshkov^[5]利用实验对其进行了验证.ICF问题中流体动力学不稳定性发生在曲面几何中,研究汇聚冲击和可控界面间的RM不稳定性对实现ICF点火非常关键,人们由此开展了一系列实验,其中激波管成为重要实验工具.

Takayama等^[6-8]设计了一种垂直环形同轴激波管,可产生柱状曲面激波,研究曲面激波在气柱界面上传播而引起的RM不稳定性.Benjamin和Fritz^[9]使用凝聚相爆轰的办法生成了一种可用来研究RM不稳定性的强激波,消除了激波管由于机械强度的限制而难以产生高马赫激波的缺点.文献^[10,11,12]设计了一种激波管试验段的曲面壁轮廓,此种装置可将平面激波转变为柱形激波,研究了入射激波马赫数、汇聚角等对激波汇聚过程的影响.Si等^[13]利用实验研究了圆柱形汇聚冲击波与多边形重气柱的相互作用,界面由空气与六氟化硫重气形成,界面形状分别为正八边形、正方形以及等边三角形.实验中在高低压腔室形成激波,激波在试验段环形管中向内汇聚,从而对中心处的物质界面产生扰动.实验结果获得了界面的变形过程及振幅变化曲线.文献^[14,15,16]设计了一种半环形激波管,研究了受圆柱汇聚冲击波影响的单模界面的RM不稳定性.

除了实验研究,人们还对汇聚RM不稳定性进行数值模拟.相比于实验研究,数值模拟可获得更多精确的结果及流场信息.Tian等^[17]在柱坐标系下采用二维可压缩欧拉方程模拟了圆柱形汇聚冲击波扰动正弦曲线界面的RM不稳定性,结果表明,Atwood数绝对值越大,则界面振动幅度越大.Zhang等^[18]利用五阶和七阶WENO格式开发了一套程序,计算了不同强度的平面激波冲击气柱界面的RM不稳定性.Lombardini等^[19-20]在扰动的球形界面处进行湍流混合的大涡模拟,观察了界面受到球形爆炸冲击波影响的演化现象.如前所述,这些数值模拟方法大多都是基于网格离散算法,有欧拉网格和欧拉、拉格朗网格两种方法对物理方程进行描述.网格计算方法在解决流体力学问题方面已较为成熟^[21-22],但对于一些特殊问题,往往具有较大的局限性.欧拉网格方法中,在固定的网格上确定运动交界面、自由表面等具有较大困难.欧拉、拉格朗日网格则需要网格变形或运动,对于大变形问题,可能会因网格畸变而出现错误.上述问题需要采用一些特殊方法进行处理^[23-26].不管哪种网格离散,由于对流项离散导致的数值扩散问题往往需要较高阶精度算法才能获得较准确的界面演化结果.即使如此,界面也需要若干网格来离散,以抑制非物理振荡.显然,对于RM不稳定性问题研究,采用能精确追踪运动界面的数值模拟方法具有重要意义.

本文试图利用光滑粒子流体动力学(smoothed particle hydrodynamics,SPH)方法来研究汇聚RM不稳定性问题.目前使用此方法来进行汇聚RM不稳定性问题的研究还很少.SPH方法是一种无网格方法^[27-30],用于模拟流体流动.该方法于1977年由Lucy^[31]和Gingold等^[32]提出,一开始是服务于天体物理学领域,用来求解相关问题.在此方法中,材料由一系列可自由运动的粒子组成,这些粒子具有材料的质量、密度、速度、压力等物理性质,并按照守恒方程控制的规律运动.该方法可以模拟任意变形问题,由于采用拉格朗日方法描述流体,无对流项离散,理论上无数值耗散,且无须在网格上追踪界面,适合复杂界面跟踪、大变形等应用.国内也开展了不少关于SPH的研究工作.1995年,贝新源等^[33]发表一篇介绍SPH的论文,这是国内最早的关于SPH方法的论述.1996年,张锁春^[34]也对SPH方法进行了综述.国内主要将SPH方法应用于冲击动力学及计算流体动力学方面的研究.贝新源等^[35]将三维SPH算法应用于高速斜碰撞问题.蔡清裕等^[36]将有限元和SPH方法相结合,模拟了刚性弹丸侵彻混凝土的问题.李金柱等^[37]利用SPH方法模拟了超高速碰撞现象.孙鹏楠等^[38]利用SPH方法模拟了自由气泡上浮过程,研究其运动特性.闫蕊等^[39]利用SPH方法研究了空气对平板入水的影响.饶登宇等^[40]采用SPH方法对温度场、水分场的耦合模型进行数值求解.任选其等^[41]基于SPH/FEM耦合方法对浅水冲击的液体喷溅特性开展研究,取得了较好的研究效果.孙红^[42]建立了刀具土壤切削有限元模型,利用SPH方法对其开沟过程进行数值模拟.何鲜峰等^[43]利用采用SPH方法,模拟了溢流坝泄流过程.不过,经典的SPH算法在激波间断上容易产生非物理压力振荡,往往需要添加人工黏性来抑制振荡.物质界面两侧容易发生粒子穿透等问题,需要特殊处理.2000年,Parshikov^[44]依据黎曼解的方法,提出粒子接触算法,利用黎曼解来求解粒子互相作用过程中的物理量,可有效抑制各物理量在界面上的非物理振荡,并且此算法避免了显式加入人工黏性项、人工热量项等问题,程序简洁,且有效提升了强激波间断的模拟精度,并且基于黎曼解的接触间断算法也有效避免了粒子穿透等问题,为复杂冲击界面演化模拟提供了新的解决思路.近年来已有较多的应用开展.2006年,徐志宏等^[45]利用基于黎曼解的粒子接触算法模拟了激波管与飞片碰撞过程中波的传播问题.2009年,Hu等^[46]利用HLLC黎曼求解器处理了可压缩多相流流动的界面问题.2010年,Rogers等^[47]利用基于黎曼解的SPH-ALE的方法模拟了防波堤和波浪之间的相互作用,取得良好效果.不过,SPH粒子接触算法在汇聚RM不稳定性研究领域还鲜有应用.

2013年,Sirotkin 等^[48]推导出基于Local Lax-Friedrichs (LLF),Harten Lax van Leer(HLL)和精确黎曼求解器的SPH方程,并比较了几个求解器的性能,认为HLL和LLF近似黎曼求解器在强间断分辨和捕捉方面表现最优.本文采用基于HLL近似黎曼求解器的SPH算法,编制了计算程序并进行校核,在获得稳定、准确计算结果的前提下,针对Si等^[13]的实验,模拟了二维圆柱形汇聚冲击波冲击多边形轻/重气界面诱导RM不稳定现象的演化过程.所得结果与实验结果进行了对比,获得了定性及定量的校核结果.之后进一步探究了汇聚冲击波冲击多边形轻/重气界面诱导RM不稳定性的物理机制.

1 SPH数值方法

1.1 基于HLL黎曼求解器的动量方程与能量方程

根据文献^[48]提出的基于HLL黎曼求解器的SPH方法,得到用于数值模拟的动量方程和能量方程

(1)$$ \frac{{d}{{u}}_i }{{d}t} =-\frac{2}{\rho _i }\sum\limits_{j = 1}^N {v_j [w_i p_i + w_j p_j-w_{ij} (u_j-u_i )]\nabla _i W_{ij} } $$

(2)$$ \frac{{d}E_i }{{d}t} = \frac{2}{\rho _i }\sum\limits_{j = 1}^N {v_j [w_i p_i u_i + w_j p_j u_j-w_{ij} (E_j-E_i )]W_q } $$

其中

(3)$$ W_q = {{e}}_{ij} \cdot \nabla _i W_{ij}, {{e}}_{ij} = {{r}}_{ij} / r_{ij} , {{r}}_{ij} = {{r}}_i-{ {r}}_j $$

(4)$$w_i = \frac{b_j }{b_i + b_j },w_j = \frac{b_i }{b_i + b_j },w_{ij} = \frac{b_i b_j }{b_i + b_j } $$

(5)$$ b_K = \max (C_k ,\overline C ),K = i,j $$

(6)$$ \overline C = \frac{\sqrt {\rho _i } C_i + \sqrt {\rho _j } C_j }{\sqrt {\rho _i } + \sqrt {\rho _j } } ~ $$

$C$为拉格朗日声速, $C$由下式计算得出

(7)$$ C = \sqrt {\gamma p\rho } $$

其中,${{u}}$为速度矢量,$E$为比总能,$E = e + \dfrac{1}{2}u^2$,$e$为比内能,$p$,$u$,$\rho $分别为粒子压力、速度和密度. $v_j = {m_j }/ {\rho _j }$,表示粒子$j$的体积.

$\nabla _i W_{ij} $为光滑函数 $W_{ij} $的梯度, $\nabla _i W_{ij}$表示如下

(8)$$ \nabla _i W_{ij} = \frac{x_i-x_j }{r_{ij} }\frac{\partial W_{ij} }{\partial r_{ij} } = \frac{x_{ij} }{r_{ij} }\frac{\partial W_{ij} }{\partial r_{ij} } $$

式中,$r_{ij} $为粒子$i$和粒子$j$之间的距离. 光滑函数$W_{ij} = W(r_{ij} ,h_{ij} )$,$h_{ij} = (h_i + h_j ) / 2$为光滑长度.

1.2 连续性方程

由于密度在SPH方法中决定着光滑长度的变化以及粒子分布,因此密度计算在SPH方法中极其重要. 本文中密度计算采用密度求和法

(9)$$ \rho _i = \sum\limits_{j = 1}^N {m_j W_{ij} } $$

密度求和法体现SPH近似的本质,其可保证积分在整个问题域内严格遵守质量守恒定律,保证总质量守恒.

但是,需要指出的是,文献^[48]给出的上述密度计算公式,在界面附近并不区分流体介质类型. 这就使得流体密度在界面附近会出现一定的平滑. 这是文献^[48]为保证界面附近不出现非物理压强振荡而设计的一种简单化处理方法. 这看起来似乎平滑了界面,模糊了接触间断,但由于是拉格朗日算法,界面的位置实际可由粒子的介质类型来标记,因此,实际上并不会模糊对接触间断位置即界面的判别. 当然,这里也存在离散误差,这是由SPH离散算法决定的. 可以通过加密粒子的初始间距来降低其影响. 显然,粒子初始间距越小,精度越高. 本文将在第一个校核算例中明确这一特征.

1.3 光滑函数

$W_{ij}$ 为光滑函数,其满足单调性、光滑性、规范性、对称性、收敛性、正定性几个条件.本文使用分段五次光滑函数

(10)$$ W(q) = \alpha _{d} \left\{\!\!\begin{array}{lll} (3 \!- \!R)^5 - 6(2 - R)^5 + 15(1-R)^5,&\mbox{ 0} \le R{ < 1} \\ (3-R)^5-6(2-R)^5,&\mbox{ 1} \le R{ < 2} \\ (3-R)^5,&\mbox{ 2} \le R{ < 3} \\ 0,&R \ge 3\\ \end{array}\right. $$

$\alpha_{d}$在一维/二维/三维空间中数值分别为 $1 / 120 / h$, $7 / 478 / \pi h^2$,$3 / 359 / \pi h^3$. $R = r_{ij} / h$, $r_{ij}$ 是粒子$i$和粒子$j$之间的距离.

1.4 光滑长度

光滑长度决定了在计算某个粒子的物理量时,能够对该粒子起作用的邻近粒子数目. 粒子$i$的光滑长度为$h_i $,本文采用下式计算

(11)$$ h_i = \mu \left(\frac{m_i }{\rho _i }\right)^{1 / n} = \mu v_i ^{1 / n} $$

式中,$n$代表维度($n = 1,2,3)$,$\mu $为参数,本文中$\mu = 1.4$.

1.5 状态方程

本文所涉及的流体介质均采用理想气体状态方程及声速计算如下所示

(12)$$ p = (\gamma-1)\rho e $$

(13)$$ c = \sqrt {\gamma (\gamma-1)e} $$

式中,$\gamma $为气体常数.

1.6 固壁边界处理模型

在实际模拟中,有时需要给流体施加固壁条件. 在SPH数值模拟中,将固壁离散为边界粒子,这些边界粒子与其接触的流体粒子的初始粒子间距相同,光滑长度也相同. 这些边界粒子不参与实际计算,只是当靠近它的流体粒子运动到其影响范围时,对这些流体粒子施加一个沿两粒子中心线的强排斥力,防止流体粒子穿透固壁边界. 本文固壁边界使用由Monaghan^[49]于1994年提出的边界排斥力法、分子间相互作用的Lennard-Jones力形式

(14)$$ f(r) = \left\{ \begin{array}{lll} D\left[\left(\dfrac{r_0 }{r}\right)^{p_1 }-\left(\dfrac{r_0 }{r}\right)^{p_2 }\right]\dfrac{{ {r}}}{r^2},&r \le r_0\\ 0, &r > r_0 \\ \end{array}\right. $$

式中,$D$,$p_1 $,$p_2 $为可调参数. $p_1 > p_2 $,本文中$p_1 = 12$,$p_2 = 6$. $r$为两粒子之间的距离,$r_0 $为初始粒子间距.

1.7 时间积分

SPH方程是可以进行数值积分的运动方程. 常用的数值分析算法有蛙跳法,龙格、库塔法和预测、校正法. 由于蛙跳法计算效率高,存储需求量低,故本文采用蛙跳法进行数值积分.

蛙跳法在第一个时间步长内,积分形式如下

(15)$$ \left. \begin{array}{llll} t = t_0 + \Delta t \\ \rho _i (t_0 + \Delta t / 2) = \rho _i (t_0 ) + \dfrac{\Delta t}{2} \cdot \dfrac{{d}\rho _i (t_0 )}{{d}t} \\ e_i (t_0 + \Delta t / 2) = e_i (t_0 ) + \dfrac{\Delta t}{2} \cdot \dfrac{{d}e_i (t_0 )}{{d}t} \\ v_i (t_0 + \Delta t / 2) = v_i (t_0 ) + \dfrac{\Delta t}{2} \cdot \dfrac{{d}v_i (t_0 )}{{d}t} \\ x_i (t_0 + \Delta t) = x_i (t_0 ) + \Delta t \cdot v_i (t_0 + \Delta t / 2) \end{array} \right\} $$

在以后每一时间步的开端,粒子的密度、能量、速度按下列积分推进

(16)$$ \left. \begin{array}{lll} \rho _i (t) = \rho _i (t-\Delta t / 2) + \dfrac{\Delta t}{2} \cdot \dfrac{{d}\rho _i \left( {t-\Delta t} \right)}{{d}t} \\ e_i (t) = e_i (t-\Delta t / 2) + \dfrac{\Delta t}{2} \cdot \dfrac{{d}e_i \left( {t-\Delta t} \right)}{{d}t} \\ v_i (t) = v_i (t-\Delta t / 2) + \dfrac{\Delta t}{2} \cdot \dfrac{{d}v_i \left( {t-\Delta t} \right)}{{d}t} \end{array}\right\} $$

在后面时间步的最后,粒子密度、能量、速度、位置按下列积分推进

(17)$$ \left. \begin{array}{lll} t = t + \Delta t\\ \rho _i (t + \Delta t / 2) = \rho _i (t-\Delta t / 2) + \Delta t \cdot \dfrac{{d}\rho _i \left( t \right)}{{d}t}\\ e_i (t + \Delta t / 2) = e_i (t-\Delta t / 2) + \Delta t \cdot \dfrac{{d}e_i \left( t \right)}{{d}t}\\ v_i (t + \Delta t / 2) = v_i (t-\Delta t / 2) + \Delta t \cdot \dfrac{{d}v_i \left( t \right)}{{d}t}\\ x_i (t + \Delta t) = x_i (t) + \Delta t \cdot v_i (t + \Delta t / 2) \end{array} \right\} $$

1.8 时间步长

蛙跳法需要满足一定的稳定条件. 稳定条件是CFL(Courant-Friedrichs-Levy)条件. 本文中,时间步长由下式计算

(18)$$ \Delta t = \min \left( {\frac{\xi h_i }{h_i \nabla \cdot v_i + c_i + 1.2(\alpha _\Pi c_i + \beta _\Pi h_i \left| {\nabla \cdot v_i } \right|)}} \right) $$

式中,$\xi$为Courant系数,其值一般在0.3左右. 系数$\alpha_{\prod} = 1$,$\beta _{\prod}= 2$.

2 一维激波管的SPH模拟

为验证本文自编SPH程序的可靠性,首先模拟了四类经典激波管问题,分别为Sod问题、强激波间断试验、Sjögreen测试和强冲击碰撞问题. 激波管里充满气体,气体被隔膜分成左右两个部分,每个部分有各自的气体状态. 当突然把隔膜去掉,左右不同状态的气体将会产生作用,管内会产生冲击波、膨胀波和不连续接触间断. 激波管问题有精确解,便于与数值解进行对比. 本模拟中,所有粒子质量都相同. 有

(19)$$ \frac{\rho _{L} }{\rho _{R} } = \frac{\Delta x_{R} }{\Delta x_{L} } $$

式中,$\rho $为粒子密度,$\Delta x$为初始粒子间距.

对本文中所有激波管问题,分别使用3个不同的初始粒子间距进行模拟,并与精确解进行对比,探究程序的可靠性以及初始粒子间距对模拟结果的影响.

2.1 Sod问题

$t=0$时刻初始条件为： $ \rho _{L} \!=\! 1.0$ , $ u_{L} \!=\! 0.0$ , $ p_{L} \!=\! 1.0$ , $ \Delta x_{L} \!=\! 0.0015$ , $ N_{L} = 360$ ; $ \rho _{R} \!=\! 0.125$ , $ u_{R} \!=\! 0.0$ , $ p_{ R} \!=\! 0.1$ , $ \Delta x_{R} \!=\! 0.0117, N_{R} \!=\! 840$. 其中,$N_{L} $,$N_{R} $为左右状态气体的粒子数.

在此基础上,分别将初始粒子间距缩小为原先的1/5和1/10,图1给出了Sod问题在$t=0.2$ s时刻激波的密度、压力、能量、速度模拟曲线与精确解的对比结果. 其中,绿色、蓝色、红色虚线分别表示初始粒子间距为$\Delta x$,$\dfrac1 5\Delta x$,$\dfrac{1}{10}\Delta x$的激波物理量模拟曲线,黑色实线为精确解. 由图对比可见：

(1)各工况均与精确解对比良好,且随着初始粒子间距减小,激波间断的精度逐渐提高.

(2)从密度图来看,密度显示的接触间断确实存在一定的密度平滑现象,且粒子初始间距越大越明显. 随着粒子初始间距减小,平滑程度减弱. 但粒子初始间距减小到一定程度后,也不能完全消除其平滑趋势.

(3)但由物质ID标示的接触间断与理论解位置是一致的,且粒子初始间距大小对界面位置标示影响很小,三个间距结果几乎重合,只有细小的差别. 可见,由物质ID标示的接触间断位置受到间断密度平滑算法的影响很小,仍可以较精确地追踪流体界面位置.

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 $t=0.2$时刻Sod激波密度、压力、速度、能量的SPH解(虚线,颜色为不同粒子间距工况,双点化线为粒子物质属性标示的\\ 接触间断位置)与精确解(黑色实线)

Fig. 1 1D SPH solution (dashed line, color represents different$\Delta x$cases; double dot dash line is material ID of particles, which marks the contact discontinuity) and exact solution (black solid line) of density, pressure, velocity and energy of Sod's shock at $t=0.2$}

2.2 强激波间断试验

为验证本文采取的SPH算法对强激波间断问题计算的准确性,开展了$1.0\times$10$^{5}$压比1D强激波测试. 其在$t=0$时刻初始条件为： $ \rho _{L} = 1.0, u_{L} = 0.0, p_{L} = 1000, \Delta x_{L} = 0.0026, N_{L} = 404; $ $ \rho _{R} = 1.0, u_{R} = 0.0, p_{R} = 0.01, \Delta x_{R} = 0.0026, N_{R} = 196; $ 其中,$N_{L} $,$N_{R} $为左右状态气体的粒子数.

在此基础上,分别将初始粒子间距缩小为原先的1/5和1/10,图2给出了爆炸波测试问题在$t=0.012$ s时刻激波的密度、压力、能量、速度模拟曲线与精确解的对比结果. 其中,绿色、蓝色、红色虚线分别表示初始粒子间距为$\Delta x$,$\dfrac1 5\Delta x$,$\dfrac{1}{10}\Delta x$的激波物理量模拟曲线,黑色实线为精确解.

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 $t=0.012$时刻强激波的密度,压力,速度,能量的SPH解(虚线)与精确解(黑色实线)

Fig. 2 1D SPH solution (dash line) and exact solution (black solid line) of density, pressure, velocity and energy of strong shock at $t=0.012$

可见,强激波条件下,各物理量与精确解对比良好,且随着初始粒子间距减小,强激波间断的精度逐渐提高. 在初始$1.0\times$10$^{5}$压比条件下,产生的强激波间断被1D SPH清晰地捕捉到,且间断附近未产生非物理振荡,说明本文SPH算法可有效模拟强激波间断问题.

2.3 Sjögreen测试

Sjögreen测试问题在$t=0$时刻初始条件为： $ \rho _{L} = 1.0$, $u_{L} =-2.0, p_{L} = 0.4$, $\Delta x_{L} = 0.0011$, $N_{L} = 600$; $ \rho _{R} = 1.0$, $u_{R} = 2.0$, $p_{ R} = 0.4$, $\Delta x_{R} = 0.0011$, $N_{R} = 600$; 其中,$N_{L} $,$N_{R} $为左右状态气体的粒子数.

在此基础上,分别将初始粒子间距缩小为原先的1/5和1/10,图3给出了Sjögreen测试问题在$t=0.15$ s时刻激波的密度、压力、能量、速度模拟曲线与精确解的对比结果. 其中,绿色、蓝色、红色虚线分别表示初始粒子间距为$\Delta x$,$\dfrac1 5\Delta x$,$\dfrac{1}{10}\Delta x$的激波物理量模拟曲线,黑色实线为精确解.

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 $t=0.15$时刻激波的密度,压力,速度,能量的SPH解(虚线)与精确解(黑色实线)

Fig. 3 1D SPH solution (dotted line) and exact solution (black solid line) of density, pressure, velocity and energy of shock at $t=0.15$

由图可见,除接触间断附近能量和压力与理论解存在一定偏差外,其余解均与理论解吻合良好. 随着初始粒子间距减小,激波间断的精度逐渐提高.

需要指出的是,接触间断附近能量解与理论解的差异在其他算法的求解中也存在类似的现象,特别是网格类算法. 对于SPH方法在这种情况下也产生相似问题,Monaghon^[50]认为这是由于粒子质量对其自身密度的贡献总是限制了粒子最小的密度. Cha和Whitworth^[51]猜测是黎曼求解器造成的,可以用不同的黎曼求解器来解决这个问题. 但是,本文中的黎曼求解器和他们使用的求解器不同,所以这个问题应该也不是求解器造成的. 这有待于未来进一步解决.

2.4 强冲击碰撞

强冲击碰撞问题在$t=0$时刻初始条件为 $ \rho _{L} = 5.999$, $u_{L} = 19.598$, $p_{L} = 460.894$, $\Delta x_{L} = 0.0044$, $N_{L} = 396$; $ \rho _{R} = 5.999$, $u_{R} =-6.196$, $p_{R} = 46.095$, $\Delta x_{R} = 0.0044$, $N_{R} = 204$; 其中,$N_{L} $,$N_{R} $为左右状态气体的粒子数.

在此基础上,分别将初始粒子间距缩小为原先的1/5和1/10,图4给出了强冲击碰撞问题在$t=0.035$ s时刻激波的密度、压力、能量,速度模拟曲线与精确解的对比结果. 其中,绿色、蓝色、红色虚线分别表示初始粒子间距为$\Delta x$,$\dfrac1 5\Delta x$,$\dfrac{1}{10}\Delta x$的激波物理量模拟曲线,黑色实线为精确解.

由图可见,除接触间断附近密度有所平滑外,其余均与理论解吻合良好. 并且随着初始粒子间距减小,强激波间断的精度逐渐提高.

综上可见,本文采用的SPH算法,实现了对强激波和界面强间断的有效分辨和追踪. 当需要界面的实际位置时,拉格朗日粒子方法优于欧拉网格方法. 上述一维数值校核充分证明了代码的可靠性、健壮性和准确性. 由于该方法所有SPH公式以矢量形式给出,且在任何地方都能产生平滑解,适合高维扩展.

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 $t=0.035$时刻强激波的密度,压力,速度,能量的SPH解(虚线)与精确解(黑色实线)

Fig. 4 1D SPH solution (dotted line) and exact solution (black solid line) of density, pressure, velocity and energy of strong shock at $t=0.035$

3 柱形汇聚激波冲击多边形界面的SPH模拟

3.1 物理模型

将前述发展的SPH程序,由一维推广到二维. 对Si等^[13]实验获得的柱形汇聚激波冲击多边形动界面过程进行数值模拟,界面使用正方形,界面尺寸以及激波强度与实验保持一致. 计算模型见图5.

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 计算模型示意图

Fig. 5 Schematic diagram of the computational model

图5中,$r_0 = 470~\mbox{mm}$,$r_1 = 600~\mbox{mm}$,$r_{s} = 603~\mbox{mm}$. 正方形区域为${SF}_6 $气体,点$a$为正方形顶点,$oa = 20~\mbox{ mm}$,为正方形顶点到其中心的距离,点$b$为正方形边界中心点. ${ I}$区为低压空气,$\Pi $区为高压空气,由高低压产生激波并向中心汇聚,对${SF}_6 $物质边界产生扰动. 最外围为固壁边界.

实验中,激波在距中心70 mm的位置开始向内汇聚,初始马赫数$M_0 = 1.17$. 为防止计算域过小而导致高低压边界对实际模拟产生影响,将计算域扩大,高低压边界置于距中心470 mm处. 由距中心70 mm处的马赫数1.17,根据CCW关系式^[52]可计算出距中心470 mm处的马赫数$M_1 =1.072$. 由入射激波马赫数$M_1 $与初始压比$P_{41} $的关系式计算出初始高低压之比. 计算公式如下

(20)$$ P_{41} \!=\! \left[ {1 + \frac{2\gamma }{\gamma + 1}(M_1 ^2-1)} \right]\left[ {1 - \frac{\gamma-1}{\gamma + 1}a_{14} (M_1-{M_1}^{-1})} \right]^{-\frac{2\gamma }{\gamma-1}} $$

式中,取常数空气$\gamma = 1.4$,初始压室中,初始温度比$a_{14} = 1$. {SF}$_{6}$气体与空气的参数如表1所示. 实际模拟中,所有粒子质量都相同,由于是二维模型,故

(21)$$ \frac{\Delta x_{Air} }{\Delta x_{{SF}_6 } } = \sqrt {\frac{\rho _{{SF}_6 } }{\rho _{Air} }} $$

表1 气体参数

Table 1 Gas parameters

新窗口打开| 下载CSV

3.2 模拟结果与讨论

利用SPH程序模拟二维圆柱形汇聚冲击波冲击多边形轻/重气界面诱导RM不稳定现象的演化过程,模拟中分别使用了两种粒子分布情况,如表2所示.

表2 粒子分布情况

Table 2 Particle distribution

新窗口打开| 下载CSV

下面对模拟中界面演化过程进行定量分析.图6给出了从激波扰动界面开始,界面侧面中心到焦点$o$的归一化位移随无量纲时间的变化情况,并与参考文献^[13]的实验曲线进行对比.

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 界面侧面中心到焦点$o$的归一化位移随无量纲时间的变化情况

Fig.6 Normalized displacement of the interface side center to the focal point $o$ as a function of dimensionless time

由图6可以看出,模拟所得曲线与实验所得曲线较为符合,且Case 2的模拟结果更为精确. 可见,要精确追踪界面位置,SPH需要一定的粒子密度,这类似于欧拉网格类算法中网格加密对结果精度的影响,即需要保证一定的离散精度.

图中$Ls$和$Lc$分别为界面侧面中心和初始界面顶点到焦点$o$的距离,$Lc=20~\mbox{ mm}$. 无量纲时间$tc = Lc / W_{i0} $,$W_{i0} = 444.8~\mbox{ m} /\mbox{s}$,为入射冲击波到达初始界面顶点的速度. 当汇聚冲击波到达界面顶点,然后在界面传播(透射和反射),界面受激波冲击影响向内收缩,$Ls$开始减小. 当汇聚激波到达界面中心,$Ls$达到最小值,之后激波开始向外反射,界面侧面中心向外运动,$Ls$逐渐增大.

图7为模拟与实验所得界面演化对比图. 其中,图像左边为Case 2数值模拟得到的界面,右边为参考文献^[13]的实验利用偏光技术拍摄所获得的界面,其中散光部分为界面肥皂膜破碎后形成的散射微小液滴,标示了界面的演化.

由上述模拟与实验所获得的界面对比图可以看出,模拟结果与实验结果具有较好的一致性. 界面在激波作用下出现尖钉状、气泡状相位反转现象. SPH模拟获得的界面清晰、准确,无数值扩散带来的模糊边界图像.

图8和图9分别给出了Case 2界面演化过程中密度及压力变化. 可以看出,汇聚激波具有指向界面中心的速度,到达四边形{ SF}$_{6}$界面前形状为圆形,到达四边形${SF}_{6}$界面四个顶点时,产生反射和透射激波. 透射激波继续向界面中心汇聚,不过由于界面的作用,已不是圆形,而是近似四边形. 相邻界面的透射激波相交形成透射反射激波. 当汇聚激波完全透过界面,继续向中心传递的透射激波形状变为弯边方形,反射激波则完全离开界面,并向外围传递. 随着透射激波逐渐向界面中心汇聚,界面也被迅速压缩向中心汇聚,界面侧面中心移动速度最快,界面顶点移动速度较慢,从而在界面侧面中心形成气泡状结构,在界面顶点形成尖钉状结构.

当透射激波在界面中心$o$点汇聚,相邻的透射反射激波在界面中心线处相遇,形成十字形非均匀密度区. 透射激波在中心汇聚后形成二次反射激波,以弯边方形波阵面向外传播. 随着二次反射激波向外运动,首先遇到汇聚过程中被压缩不规则的{ SF}$_{6}$界面边界,产生二次透射激波及三次反射激波,继续与界面作用,流场中的涡量分布被改变, 进而出现尖钉状、气泡状相位反转现象. 可见,本文采用的SPH算法准确地追踪了激波与界面相互作用的复杂界面和波系演化规律.

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 模拟界面与参考文献^[13]的实验获得界面对比图

Fig. 7 Comparison of simulated interface and experimental interface evolution of Ref.^[13]

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 密度云图

Fig. 8 Density contour evolution

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 密度云图(续)

Fig. 8 Density contour evolution (continued)

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 压力云图

Fig. 9 Pressure cloud map

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 压力云图(续)

Fig. 9 Pressure cloud map (continued)

4 结论

采用基于HLL黎曼求解器的SPH算法,实现了对强激波和界面强间断的有效分辨和追踪,主要获得如下结论：

(1)基于HLL黎曼求解器的SPH程序可有效模拟激波间断及接触界面复杂演化问题,一维数值验证所得结果与精确解均较为吻合,证明了程序的可靠性,合理性以及健壮性.

(2)本文SPH采用的密度算法会一定程度平滑接触间断处的密度差异,可通过减小初始粒子间距来改善这一影响. 初始粒子间距越小,则数值模拟精度越高,但不能完全消除. 不过,界面的位置实际是由粒子的介质属性来标记,在SPH拉格朗日算法下不会影响界面位置的判别.

(3) 对于汇聚激波冲击多边形界面问题,SPH模拟获得了清晰界面演化,没有数值扩散,界面演化过程与实验结果十分吻合,初步证明了SPH算法可有效模拟此类问题. 随着初始粒子间距减小,界面的追踪精度也随之提高.