Numerical Simulations of Shock Problems With the Revised KDF-SPH Method
-
摘要: 基于光滑粒子流体动力学(smooth particle hydrodynamics, SPH)方法中的光滑核近似和Taylor级数展开原理，利用核函数矩对KDF-SPH(kernel derivative free SPH)方法进行了修正.为了验证修正方法的适用性和可行性，将该方法应用于不同情况下一维激波管问题的数值模拟，并对模拟结果进行分析.结果表明，修正方法能很好地捕捉到激波和接触不连续的位置和强度，修正方法不要求核函数可导性，不计算核函数矩，计算量更小，计算效率更高.Abstract: Based on the smooth kernel approximation and the Taylor series expansion of the smooth particle hydrodynamics (SPH) method, the kernel function moment was used to revise the KDF-SPH (kernel derivative free SPH) method. To prove the applicability and feasibility of the proposed revised scheme, the scheme was applied to the numerical simulations of 1D shock tube problems under different conditions, and the simulation results were analyzed. The results show that, the revised method can well capture the positions and strengths of shock waves and contact discontinuities. The revised method does not require the derivability of the kernel function, does not calculate the kernel function moment, and has a smaller computation cost with a higher calculation efficiency.
-
0. 引言
光滑粒子流体动力学(SPH)方法是一种无网格纯Lagrange型数值方法.它最初由Lucy[1]提出，用于解决三维开放空间的天体物理学问题，由于在这些问题中质点的运动与流体或气体很相似，SPH方法也被广泛应用于模拟流体结构相互作用、内部波、水力涡流、液面波、多相流[2-3]等情况.
在SPH方法中，场函数及其导数的核近似和粒子近似都需要计算核梯度，这导致了对核函数的一些特殊要求，如核函数应该是可微的并且足够平滑，这些要求通常限制了对核函数的选择.Feng等[4]在SPH方法的基础上提出了核函数矩的概念，通过计算核函数矩建立了一种避免求解核函数导数的SPH(KDF-SPH)方法，该方法由于不需要求解核函数导数因而使得核函数的选择更加灵活.
Huang等[5]提出了一种无核梯度SPH(KGF-SPH)方法，在整个计算过程中不必计算核梯度，当核函数不可微或由此产生的核梯度不够平滑时，增加了核函数的可选择性.Maatouk[6]通过Taylor级数展开，用类Newton迭代方法替换导数，构造了一种求解非线性系统的无导数SPH(DF-SPH)迭代方法.Imin等[7-8]针对传统SPH方法和DF-SPH方法，提出了一种计算函数一阶导数和二阶导数的高精度改进方法，并将改进方法应用于具体问题，取得了良好的计算结果，验证了改进措施的有效性.Huang等[9]提出了一种采用迭代粒子移位技术(PST)的KGF-SPH方法来模拟翼型周围的流动，该方法可以在不使用核梯度的情况下保持高精度.Garoosi等[10]建立了一种新的无核导数弱可压缩SPH(KDF-WCSPH)模型，用于模拟自由表面流动和对流换热.当用KGF-SPH方法模拟黏性不可压缩流时，通常使用速度的嵌套近似来近似Laplace项，这可能会引入嵌套近似的数值误差，也会导致处理边界条件的困难.Huang等[11]提出了一种改进的KGF-SPH方法，用Laplace算子的一种新的离散格式来模拟黏性、不可压缩流体流动，改进的KGF-SPH方法避免了一阶导数的嵌套逼近，保持了“无核梯度”的良好特性.
激波管问题是一类经典的数学物理问题，是描述激波传播过程中，非定常、非平衡的复杂物流运动行为的数学模型，要求考虑质量、动量和能量的守恒关系，以及流体热力学状态方程的影响.激波管问题本身有精确解，在数值模拟中，激波管问题是检验数值方法的重要手段.王建玲等[12]通过Euler方程组经典算例测试了一种三阶WENO差分格式的低耗散、高分辨特性.张成治等[13]通过MHD方程的激波算例测试了一种四阶熵稳定半离散有限体积格式，该格式具有低耗散、高精度、无振荡并且可以有效捕捉间断等特点.1983年，Monaghan和Gingold[14]首次尝试应用SPH方法模拟激波管问题，利用SPH近似方法描述Lagrange型Euler方程，并在动量方程中引入M&G83型人工黏性项处理激波中的非物理振荡现象.Li等[15]在多相冲击管状问题、冲击波对多流体界面的影响问题中验证和讨论了一种基于插值方法的可压缩SPH模型.Meng等[16]在涉及接触不连续性、极端冲击波和强稀疏波试验中验证了一种基于SPH方法的冲击捕获方案.Wang等[17]提出了一种SPH方法的WENO格式实现方法，该方法在具有不连续和小尺度结构的可压缩流动中能够准确捕捉冲击波.Sirotkin和Yoh[18]利用Riemann解算器对SPH方法进行改进，在粒子数为600的情况下对Sod问题、Sjögreen试验、冲击波试验和强冲击碰撞问题进行了模拟，均得到较好的模拟结果.徐建于和黄生洪[19]采用基于HLL Riemann求解器的SPH算法，对部分一维激波管问题在粒子总数为600的情况下进行模拟，证明了算法的可靠性，实现了对强激波和大密度比物质界面的有效分辨和追踪.
本文通过直接代入核函数矩的值对Feng等[4]提出的KDF-SPH方法进行了修正，并编写了相应程序模拟了不同情况下的一维激波管问题.修正方法在整个计算过程中不需要对核函数导数及核函数矩进行推导和求解.修正方法的有效性在数值模拟问题中得到了充分体现.本文其余部分结构如下：第1节中简要介绍了传统SPH方法的基本原理；第2节中介绍了修正KDF-SPH方法；第3节介绍了控制方程；第4节讨论了不同情况下一维激波管问题的数值模拟情况；第5节给出了全文的结论.
1. 传统SPH方法
SPH方程的构造通常分为两步：第一步是核近似，第二步是粒子近似.连续函数f(x)的核近似表达式为
$$ \langle f(x)\rangle=\frac{1}{h^{d}} \int_{\varOmega_{{W}}} f\left(x^{\prime}\right) W\left(\frac{x-x^{\prime}}{h}\right) \mathrm{d} x^{\prime}, $$ (1) 函数偏导数的核近似表达式为
$$ \left\langle\frac{\partial f(x)}{\partial x_{\alpha}}\right\rangle=\frac{1}{h^{d}} \int_{\varOmega_{W}} f\left(x^{\prime}\right) \frac{\partial}{\partial x_{\alpha}^{\prime}}\left[W\left(\frac{x-x^{\prime}}{h}\right)\right] \mathrm{d} x^{\prime}, $$ (2) 式中d为计算维数，$W\left(x-x^{\prime}, h\right)=\frac{1}{h^{d}} W\left(\frac{x-x^{\prime}}{h}\right)$为光滑核函数.光滑核函数取决于两个变量：两点的距离|x-x′|和光滑长度h.核函数影响区域ΩW在一维、二维和三维情形下分别为区间、圆和以x为中心、kh为半径的球，不同的核函数k不一样.
在SPH方法中把整个求解域离散成一系列任意分布的有质量和容量的有限多个粒子.如果核近似式中粒子j处微元体dx′的体积为ΔVj，质量为mj，密度为ρj，位置坐标为xj，记函数f(x)在粒子j处的值为f(xj)，假设以粒子i为中心的核函数影响范围内的粒子总数为N，由式(1)、式(2)可得函数f(x)在粒子i处的粒子近似式及一阶偏导数粒子近似式：
$$ \left\langle f\left(x_{i}\right)\right\rangle=\sum\limits_{j=1}^{N} \frac{m_{j}}{\rho_{j}} f\left(x_{j}\right) W_{i j}, $$ (3) $$ \left\langle\frac{\partial f(x)}{\partial x_{\alpha}}\right\rangle=\sum\limits_{j=1}^{N} \frac{m_{j}}{\rho_{j}} f\left(x_{j}\right) \frac{\partial}{\partial x_{\alpha}^{\prime}} W_{i j} , $$ (4) 式中Wij为关于粒子i和j的核函数.
2. 修正KDF-SPH方法
函数f(x′)在点x处进行Taylor级数展开，可得
$$ \begin{gather*} f\left(x^{\prime}\right)=f(x)+f^{\prime}(x)\left(x^{\prime}-x\right)+\frac{1}{2} f^{\prime \prime}(x)\left(x^{\prime}-x\right)^{2}+\frac{1}{6} f^{\prime \prime \prime}(x)\left(x^{\prime}-x\right)^{3}+\cdots= \\ f(x)+h f^{\prime}(x) \frac{x^{\prime}-x}{h}+\frac{h^{2}}{2} f^{\prime \prime}(x)\left(\frac{x^{\prime}-x}{h}\right)^{2}+\frac{h^{3}}{6} f^{\prime \prime \prime}(x)\left(\frac{x^{\prime}-x}{h}\right)^{3}+\cdots . \end{gather*} $$ (5) Taylor展开式(5)的两边同时乘$\frac{x^{\prime}-x}{h} W\left(x^{\prime}-x, h\right)$并在核函数影响区域Ω上取积分，可得
$$ \begin{gather*} \int_{\varOmega} \frac{x^{\prime}-x}{h} f\left(x^{\prime}\right) W \mathrm{~d} x^{\prime}=\int_{\varOmega} \frac{x^{\prime}-x}{h} f(x) W \mathrm{~d} x^{\prime}+\int_{\varOmega}\left(\frac{x^{\prime}-x}{h}\right)^{2} h f^{\prime}(x) W \mathrm{~d} x^{\prime}+ \\ \int_{\varOmega}\left(\frac{x^{\prime}-x}{h}\right)^{3} \frac{h^{2}}{2} f^{\prime \prime}(x) W \mathrm{~d} x^{\prime}+\int_{\varOmega}\left(\frac{x^{\prime}-x}{h}\right)^{4} \frac{h^{3}}{6} f^{\prime \prime \prime}(x) W \mathrm{~d} x^{\prime}+\cdots= \\ M_{1} f(x)+h M_{2} f^{\prime}(x)+\frac{h^{2}}{2} M_{3} f^{\prime \prime}(x)+\frac{h^{3}}{6} M_{4} f^{\prime \prime \prime}(x)+\cdots, \end{gather*} $$ (6) 式中
$$ \begin{align*} M_{i}= & \frac{1}{h} \int_{\varOmega}\left(\frac{x-x^{\prime}}{h}\right)^{i} W\left(\frac{x-x^{\prime}}{h}\right) \mathrm{d} x^{\prime}= \\ & \frac{1}{h} \int_{x-k h}^{x+k h}\left(\frac{x-x^{\prime}}{h}\right)^{i} W\left(\frac{x-x^{\prime}}{h}\right) \mathrm{d} x^{\prime} \stackrel{s=\frac{x-x^{\prime}}{h}}{\Rightarrow} \int_{-k}^{k}(s)^{i} W(s) \mathrm{d} s \end{align*} $$ (7) 是一维核函数矩，可以证明Mi只与核函数有关且并不依赖于h.
核函数的归一性可表示为
$$ \int_{\varOmega} W\left(x-x^{\prime}, h\right) \mathrm{d} x^{\prime}=1. $$ 该性质确保了在光滑函数支持域上的积分是归一的，此性质同样确保了连续函数积分表达式的零阶连续性(C0).根据式(7)可得
$$ M_{1}=\frac{1}{h} \int_{\varOmega}\left(\frac{x-x^{\prime}}{h}\right) W\left(\frac{x-x^{\prime}}{h}\right) \mathrm{d} x^{\prime}=\frac{1}{h} \int_{x-k h}^{x+k h}\left(\frac{x-x^{\prime}}{h}\right) W\left(\frac{x-x^{\prime}}{h}\right) \mathrm{d} x^{\prime} \stackrel{\frac{x-x^{\prime}}{h}}{\Rightarrow} \int_{-k}^{k} s W(s) \mathrm{d} s. $$ 由于核函数为偶函数(对称性)，那么sW(s)是奇函数，奇函数在对称区间上的积分值为0，则有M1=0，由式(6)可得
$$ f^{\prime}(x)=\frac{\left\langle\frac{x^{\prime}-x}{h} f(x)\right\rangle}{h M_{2}}+o\left(h^{2}\right) . $$ (8) 式(8)即为Feng等[4]提出的KDF-SPH方法计算一阶偏导数的表达式.KDF-SPH方法虽然相较于传统SPH方法不用计算核函数导数，但却需要计算核函数矩的值.
当粒子在核函数影响区域内均匀分布时，核函数矩M2是一个常数.M2的值可以通过计算的空间维数和选择的核函数来确定[9].根据Fulk和Quinn[20]对一维光滑粒子流体力学核的分析，如果在计算中使用三次B-样条(B-spline)核函数[21]，那么M2=1/3，式(8)可以表示为
$$ f^{\prime}(x)=\frac{3}{h}\left\langle\frac{x^{\prime}-x}{h} f(x)\right\rangle+o\left(h^{2}\right) . $$ (9) 式(9)是本文中计算一阶偏导数的修正KDF-SPH方法的表达式.
表 1给出了修正KDF-SPH方法、KDF-SPH方法、传统SPH方法计算量的对比.由表 1可以看出，修正KDF-SPH方法既不用求核导数，也不用计算核函数矩，计算量明显减少.
表 1 计算量对比Table 1. Comparison of calculationsrevised KDF-SPH KDF-SPH conventional SPH nuclear derivative needless needless needed kernel function moment needless needed needless 3. 控制方程
流体的状态变化满足三个物理守恒定律，即质量守恒、动量守恒和能量守恒.通过这三个守恒定律，可以推导出相应的连续性方程、动量方程以及能量方程.无黏可压理想流体一维Euler方程[22]为
$$ \left\{\begin{array}{l} \frac{\mathrm{d} \rho}{\mathrm{~d} t}=-\rho \nabla \cdot v, \\ \frac{\mathrm{~d} v}{\mathrm{~d} t}=-\frac{\nabla \cdot p}{\rho}, \\ \frac{\mathrm{~d} e}{\mathrm{~d} t}=-\frac{p}{\rho} \nabla \cdot v, \end{array}\right. $$ (10) 式中p，ρ，v和e分别为压强、密度、速度和内能.为了使方程闭合，需要通过添加状态方程来保证系统的完整性，其表达式为
$$ p=(\gamma-1) \rho e, $$ (11) 式中常数项系数γ=1.4.
在激波管问题中涉及到不同压力、密度的粒子之间的相互作用.涉及多相流、异重流等现象时，在不同类型的流体界面上容易出现数值不稳定现象，需要使用特殊的算法对控制方程进行处理.在本文中主要是通过在动量方程和能量方程的离散过程中添加人工黏性项，来减轻数值结果在接触不连续面处的振荡现象.广泛使用的M&G83型人工黏性Πij的表达式为
$$ \varPi_{i j}= \begin{cases}\frac{-\alpha \bar{c}_{i j} \mu_{i j}+\beta \mu_{i j}^{2}}{\bar{\rho}_{i j}}, & v_{i j} \cdot r_{i j} \leqslant 0, \\ 0, & v_{i j} \cdot r_{i j} \geqslant 0, \end{cases} $$ (12) 式中
$$ \begin{aligned} & \mu_{i j}=\frac{h_{i j} \nu_{i j} \cdot x_{i j}}{\left|x_{i j}\right|^{2}+\varphi^{2}}, \bar{c}_{i j}=\frac{1}{2}\left(c_{i}+c_{j}\right), \bar{\rho}_{i j}=\frac{1}{2}\left(\rho_{i}+\rho_{j}\right), h_{i j}=\frac{1}{2}\left(h_{i}+h_{j}\right), \\ & v_{i j}=v_{i}-v_{j}, x_{i j}=x_{i}-x_{j}, \bar{\rho}_{i j}=\frac{1}{2}\left(\rho_{i}+\rho_{j}\right), h_{i j}=\frac{1}{2}\left(h_{i}+h_{j}\right) ; \end{aligned} $$ α和β为常数项系数，一般取值在1.0左右；c为声速，由c2=eiγ(γ-1)计算；φ表示粒子矢量速度，用于防止粒子在相互靠近时产生的数值发散现象，φ=0.1hij.
由于质点的数量和质量保持常数，所以总质量守恒，利用式(3)有
$$ \rho_{i}=\sum\limits_{j=1}^{N} m_{j} W_{i j} . $$ (13) 在求解控制方程时，人工黏性常作为一个单独的因素加入动量方程和能量方程中.与使用传统SPH方法求解控制方程类似[14]，式(10)使用修正KDF-SPH方法求解，可表示为
$$ \left\{\begin{array}{l} \frac{\mathrm{d} \rho_{i}}{\mathrm{~d} t}=\frac{3}{h_{j}} \sum\limits_{j=1}^{N} \frac{x_{j}-x_{i}}{h_{i}}\left(v_{i}-v_{j}\right) m_{j} W_{i j}, \\ \frac{\mathrm{~d} v_{i}}{\mathrm{~d} t}=-\frac{3}{h_{j}} \sum\limits_{j=1}^{N} \frac{x_{j}-x_{i}}{h_{i}}\left(\frac{p_{i}}{\rho_{i}^{2}}+\frac{p_{j}}{\rho_{j}^{2}}+\varPi_{i j}\right) m_{j} W_{i j}, \\ \frac{\mathrm{~d} e_{i}}{\mathrm{~d} t}=\frac{3}{2 h_{j}} \sum\limits_{j=1}^{N} \frac{x_{j}-x_{i}}{h_{i}}\left(\frac{p_{i}}{\rho_{i}^{2}}+\frac{p_{j}}{\rho_{j}^{2}}+\varPi_{i j}\right)\left(v_{i}-v_{j}\right) m_{j} W_{i j} . \end{array}\right. $$ (14) 需要注意的是，在上式中，近似只使用核函数W本身，没有出现核导数Wx和Wy，无核导数使得光滑函数的选择更加灵活.
4. 不同情况的激波管问题数值模拟
如图 1所示，激波管问题[23]建立在一根封闭管道上，由隔膜分隔成两段，初始时刻t0，左侧气体的速度、密度和压力为vL，ρL，pL，右侧气体的速度、密度和压力为vR，ρR，pR，不考虑黏性，突然将隔膜抽出(或隔膜突然消失)，不同状态的气体分别会在管内产生激波、膨胀波和接触间断.那么t时刻后，管道中气体的速度、密度和压力将如何分布?
初始时刻，左右两侧气体的状态由质量相等、分布均匀的粒子表示，粒子间距受到密度之比影响，有
$$ \frac{\rho_{\mathrm{L}}}{\rho_{\mathrm{R}}}=\frac{\Delta x_{\mathrm{R}}}{\Delta x_{\mathrm{L}}}, $$ (15) 式中Δx为初始粒子间距.
本文中所有激波管问题，数值模拟中使用的初始条件和数值设置来自文献[18].设定粒子总数N=600，在x=0处设置隔膜.在模拟中对左右边界处均设置了部分虚粒子，在每个时间步下虚粒子的压力、密度、速度和能量与初始时刻保持不变.我们选取了4种不同条件下的数值模拟实验，它们包括Sod问题、Sjögreen测试、爆炸波问题和强冲击碰撞问题.分别使用传统SPH方法、KDF-SPH方法和修正方法进行模拟，由于模拟的粒子数量较多，为了便于区分清楚，选取了部分粒子的数据绘图，并与精确解进行对比，分析模拟结果.
4.1 Sod问题
Sod问题[24]是经典的测试问题之一，目的是验证不同数值方法在解决一维可压缩流的基本问题上的可靠性和精度.这个问题的初始条件涉及到左侧和右侧两个状态.在运行后，气体将从左侧初始状态重新分布到右侧，并生成一个左膨胀波、一个接触间断和一个右激波.在一维情景中，其初始条件为
$$ v_{\mathrm{L}}=0, \rho_{\mathrm{L}}=1.0, p_{\mathrm{L}}=1.0 ; v_{\mathrm{R}}=0, \rho_{\mathrm{R}}=0.125, p_{\mathrm{R}}=0.1. $$ 管道中设置左侧粒子数NL=530，右侧粒子数NR=70，左侧粒子间距ΔxL=0.001 5，右侧粒子间距ΔxR=0.012，光滑长度h=0.012，时间步长Δt=0.002 s，总模拟时间t=0.2 s.图 2是Sod问题的SPH粒子分布图，图 3是Sod问题在t=0.2 s时刻压力、密度、速度和能量模拟曲线与精确解的对比结果.
图 3 t=0.2 s时刻，Sod问题的SPH解、KDF-SPH解、修正解与精确解注 为了解释图中的颜色，读者可以参考本文的电子网页版本，后同.Figure 3. The SPH solution, the KDF-SPH solution, the revised solution and the exact solution of the Sod problem at t=0.2 s由图 3可见，在Sod问题的模拟中，修正方法的模拟结果比KDF-SPH方法的模拟结果更好.传统SPH方法和修正方法都能够准确捕获到左膨胀波和右激波，两者的模拟结果与精确解有很好的一致性.在接触间断附近，修正方法的模拟结果比传统SPH方法更好.这表明修正方法可以准确地模拟可压缩流动中的激波和接触不连续面等现象，具有不错的稳定性和可靠性.
4.2 Sjögreen测试
Sjögreen测试是另一个经典的可压缩流测试问题.它在一些方面与Sod问题类似，但是其精确解涉及到两个膨胀波而不是一个膨胀波和一个激波.当气体进入中间区域时，它会被压缩，从而形成左膨胀波和右膨胀波，这时中心位置附近形成接触不连续面.在一维情景中，其初始条件为
$$ v_{\mathrm{L}}=-2.0, \rho_{\mathrm{L}}=1.0, p_{\mathrm{L}}=0.4 ; v_{\mathrm{R}}=2.0, \rho_{\mathrm{R}}=1.0, p_{\mathrm{R}}=0.4. $$ 管道中设置左侧粒子数NL=300，右侧粒子数NR=300，左侧粒子间距ΔxL=0.002 2，右侧粒子间距ΔxR=0.002 2，光滑长度h=0.001 8，时间步长Δt=0.005 s，总模拟时间t=0.15 s.图 4给出了Sjögreen测试在t=0.15 s时刻压力、密度、速度和能量模拟曲线与精确解的对比结果.
图 4 t=0.15 s时刻，Sjögreen测试的SPH解、KDF-SPH解、修正解与精确解Figure 4. The SPH solution, the KDF-SPH solution, the revised solution and the exact solution of the Sjögreen test at t=0.15 s由图 4可见，除接触间断附近模拟结果与精确解存在一定偏差外，其余解均与精确解吻合良好.修正方法能够有效地捕捉到管内气体压力、密度、速度和能量的变化.其中，接触间断附近能量解与精确解的偏差在其他算法中也存在类似情况，Monaghan[25]认为：SPH方法也出现相似的问题是由于粒子质量对其自身密度的贡献总是限制了粒子的最小密度.
4.3 爆炸波问题
本小节模拟了Woodward和Colella[26]所描述的爆炸波问题.从气体冲击方面来看，该问题与Sod问题相似，不过初始条件更加严峻，产生的强烈冲击波对数值模拟方法的准确性和稳定性提出了更高的要求.精确解包括一个左膨胀波、一个接触不连续面和一个右侧的激波.在一维情景中，其初始条件为
$$ v_{\mathrm{L}}=0, \rho_{\mathrm{L}}=1.0, p_{\mathrm{L}}=1000 ; v_{\mathrm{R}}=0, \rho_{\mathrm{R}}=1.0, p_{\mathrm{R}}=0.01 . $$ 管道中设置左侧粒子数NL=404，右侧粒子数NR=196，左侧粒子间距ΔxL=0.002 6，右侧粒子间距ΔxR=0.002 6，光滑长度h=0.007，时间步长Δt=0.000 1 s，总模拟时间t=0.012 s.图 5给出了爆炸波问题在t=0.012 s时刻压力、密度、速度和能量模拟曲线与精确解的对比结果.
图 5 t=0.012 s时刻，爆炸波问题的SPH解、KDF-SPH解、修正解与精确解Figure 5. The SPH solution, the KDF-SPH solution, the revised solution and the exact solution of the blast wave problem at t=0.012 s由图 5可见，在初始1.0×105压强比条件下，修正方法能有效地捕捉密度的波形.同时，在压力与能量分布图中，我们可以清晰地看到在激波处，三种方法均出现了明显的数值振荡现象.在接触间断附近，传统SPH方法和KDF-SPH方法与精确解偏离较大，而修正方法与精确解吻合更为良好.综上所述，对于处理爆炸波问题，修正方法在膨胀波和接触间断处有着良好的抗扰动能力.
4.4 强冲击碰撞问题
这个问题描述了两个强激波的碰撞[27]，精确解包括一个左侧激波、一个接触不连续面和一个右侧激波.它对数值模拟方法的准确性和稳定性提出了非常高的要求，因为在激波相撞附近会产生非常强烈的物理现象，比如激波反弹和排斥效应，这些现象需要被准确地模拟才能得到可靠的模拟结果.在一维情景中，其初始条件为
$$ v_{\mathrm{L}}=19.598, \rho_{\mathrm{L}}=5.999, p_{\mathrm{L}}=460.894 ; v_{\mathrm{R}}=-6.196, \rho_{\mathrm{R}}=5.999, p_{\mathrm{R}}=46.095. $$ 管道中设置左侧粒子数NL=396，右侧粒子数NR=204，左侧粒子间距ΔxL=0.004 4，右侧粒子间距ΔxR=0.004 4，光滑长度h=0.006，时间步长Δt=0.000 035 s，总模拟时间t=0.035 s.图 6给出了强冲击碰撞问题在t=0.035 s时刻的压力、密度、速度和能量模拟曲线与精确解的对比结果.
图 6 t=0.035 s时刻，强冲击碰撞问题的SPH解、KDF-SPH解、修正解与精确解Figure 6. The SPH solution, the KDF-SPH solution, the revised solution and the exact solution of the strong shock problem at t=0.035 s由图 6可见，修正方法能有效地捕捉完整的管内气体密度和速度的变化.在压力与能量分布图中，除右激波处产生了明显的数值振荡现象，其余解均与精确解吻合良好.修正方法实现了对强激波和界面强间断的有效分辨和追踪，从而获得准确和可靠的模拟结果.
4.5 计算时长比较
为了比较修正方法、KDF-SPH方法和传统SPH方法的计算时间，对本节中的4种数值实验，保持其他条件不变，将粒子数量由N=600扩大为N=6 000，记录程序运行时长.
表 2给出了粒子数N=600以及N=6 000时修正方法、KDF-SPH方法和传统SPH方法程序运行时长比较.从表中可以看出，修正方法的程序运行时长比KDF-SPH方法和传统SPH方法都要小.当粒子数量增大时，差值也随之增大.修正方法比KDF-SPH方法和传统SPH方法的计算量更小，计算效率更高.
表 2 程序运行时长比较Table 2. Comparison of program running time costscase N t/s revised KDF-SPH KDF-SPH conventional SPH Sod problem 600 37.066 70 52.095 91 47.466 41 6 000 2 890.730 82 4 429.704 88 4 232.116 63 Sjögreen test 600 37.789 27 52.667 87 48.074 76 6 000 3 027.268 75 4 543.114 52 4 077.673 14 blast wave problem 600 46.298 74 63.743 67 57.654 41 6 000 2 974.784 74 4 603.566 58 4 017.884 46 strong shock problem 600 391.313 07 551.295 24 508.942 10 6 000 29 577.200 15 46 118.617 35 40 347.813 63 5. 结论
本文提出了一种修正KDF-SPH方法.传统SPH方法需要有核函数导数才能对场函数及其梯度进行近似，与传统SPH方法不同，修正方法不需要求解核函数导数.与KDF-SPH方法相比，修正方法不需要计算核函数矩.修正方法保留了SPH方法和KDF-SPH方法的优点，并且计算量更小，计算效率更高，特别是在处理大规模粒子系统时更加切实可行.
本文将修正方法应用于实际问题的数值实验.结果表明，修正方法在模拟这些问题时与精确解均较为吻合，具有很好的可靠性，修正方法同样能够很好地重现激波和接触不连续面的位置和强度.
在模拟实际结构时，往往需要大量的SPH粒子来表征系统.在传统SPH方法中，使用较多的粒子数量会导致计算成本大大增加.因此，研究修正KDF-SPH方法以提高计算效率，是非常必要的.修正KDF-SPH方法无疑是一个非常有前途的方向，未来还需要通过进一步的深入研究和改进，不断提高其计算效率和精度，将该方法应用到更加广泛的领域，并对不同领域的问题进行针对性的优化改进.
-
图 3 t=0.2 s时刻，Sod问题的SPH解、KDF-SPH解、修正解与精确解
注 为了解释图中的颜色，读者可以参考本文的电子网页版本，后同.
Figure 3. The SPH solution, the KDF-SPH solution, the revised solution and the exact solution of the Sod problem at t=0.2 s
图 4 t=0.15 s时刻，Sjögreen测试的SPH解、KDF-SPH解、修正解与精确解
Figure 4. The SPH solution, the KDF-SPH solution, the revised solution and the exact solution of the Sjögreen test at t=0.15 s
图 5 t=0.012 s时刻，爆炸波问题的SPH解、KDF-SPH解、修正解与精确解
Figure 5. The SPH solution, the KDF-SPH solution, the revised solution and the exact solution of the blast wave problem at t=0.012 s
图 6 t=0.035 s时刻，强冲击碰撞问题的SPH解、KDF-SPH解、修正解与精确解
Figure 6. The SPH solution, the KDF-SPH solution, the revised solution and the exact solution of the strong shock problem at t=0.035 s
表 1 计算量对比
Table 1. Comparison of calculations
revised KDF-SPH KDF-SPH conventional SPH nuclear derivative needless needless needed kernel function moment needless needed needless 
下载: 导出CSV
表 2 程序运行时长比较
Table 2. Comparison of program running time costs
case N t/s revised KDF-SPH KDF-SPH conventional SPH Sod problem 600 37.066 70 52.095 91 47.466 41 6 000 2 890.730 82 4 429.704 88 4 232.116 63 Sjögreen test 600 37.789 27 52.667 87 48.074 76 6 000 3 027.268 75 4 543.114 52 4 077.673 14 blast wave problem 600 46.298 74 63.743 67 57.654 41 6 000 2 974.784 74 4 603.566 58 4 017.884 46 strong shock problem 600 391.313 07 551.295 24 508.942 10 6 000 29 577.200 15 46 118.617 35 40 347.813 63
下载: 导出CSV
-
[1] LUCY L B. A numerical approach to the testing of the fission hypothesis[J]. The Astronomical Journal, 1977, 82: 1013. doi: 10.1086/112164 [2] HE F, ZHANG H S, HUANG C, et al. Numerical investigation of the solitary wave breaking over a slope by using the finite particle method[J]. Coastal Engineering, 2020, 156: 103617. doi: 10.1016/j.coastaleng.2019.103617 [3] HE F, ZHANG H S, HUANG C, et al. A stable SPH model with large CFL numbers for multi-phase flows with large density ratios[J]. Journal of Computational Physics, 2022, 453: 110944. doi: 10.1016/j.jcp.2022.110944 [4] FENG D Y, IMIN R. A kernel derivative free SPH method[J]. Results in Applied Mathematics, 2023, 17: 100355. doi: 10.1016/j.rinam.2023.100355 [5] HUANG C, LEI J M, LIU M B, et al. A kernel gradient free (KGF) SPH method[J]. International Journal for Numerical Methods in Fluids, 2015, 78(11): 691-707. doi: 10.1002/fld.4037 [6] MAATOUK K. Third order derivative free SPH iterative method for solving nonlinear systems[J]. Applied Mathematics and Computation, 2015, 270: 557-566. doi: 10.1016/j.amc.2015.08.083 [7] IMIN R, IMINJAN A, HALIK A. A new revised scheme for SPH[J]. International Journal of Computational Methods, 2018, 15(5): 1-17. [8] IMIN R, WEI Y, IMINJAN A. New corrective scheme for DF-SPH[J]. Computational Particle Mechanics, 2020, 7(3): 471-478. doi: 10.1007/s40571-019-00273-w [9] HUANG C, LONG T, LI S M, et al. A kernel gradient-free SPH method with iterative particle shifting technology for modeling low-Reynolds flows around airfoils[J]. Engineering Analysis With Boundary Elements, 2019, 106: 571-587. doi: 10.1016/j.enganabound.2019.06.010 [10] GAROOSI F, SHAKIBAEINIA A. Numerical simulation of free-surface flow and convection heat transfer using a modified weakly compressible smoothed particle hydrodynamics (WCSPH) method[J]. International Journal of Mechanical Sciences, 2020, 188: 105940. doi: 10.1016/j.ijmecsci.2020.105940 [11] HUANG C, LEI J M, LIU M B, et al. An improved KGF-SPH with a novel discrete scheme of Laplacian operator for viscous incompressible fluid flows[J]. International Journal for Numerical Methods in Fluids, 2016, 81(6): 377-396. doi: 10.1002/fld.4191 [12] 王建玲, 李小纲, 汪文帅. 一个改进的三阶WENO-Z型格式[J]. 应用数学和力学, 2021, 42(4): 394-404. doi: 10.21656/1000-0887.410203WANG Jianling, LI Xiaogang, WANG Wenshuai. An improved 3rd-order WENO-Z type scheme[J]. Applied Mathematics and Mechanics, 2021, 42(4): 394-404. (in Chinese) doi: 10.21656/1000-0887.410203 [13] 张成治, 郑素佩, 陈雪, 等. 求解理想磁流体方程的四阶WENO型熵稳定格式[J]. 应用数学和力学, 2023, 44(11): 1398-1412. doi: 10.21656/1000-0887.440178ZHANG Chengzhi, ZHENG Supei, CHEN Xue, et al. A 4th-order WENO-type entropy stable scheme for ideal magnetohydrodynamic equations[J]. Applied Mathematics and Mechanics, 2023, 44(11): 1398-1412. (in Chinese) doi: 10.21656/1000-0887.440178 [14] MONAGHAN J J, GINGOLD R A. Shock simulation by the particle method SPH[J]. Journal of Computational Physics, 1983, 52(2): 374-389. doi: 10.1016/0021-9991(83)90036-0 [15] LI M K, ZHANG A M, PENG Y X, et al. An improved model for compressible multiphase flows based on smoothed particle hydrodynamics with enhanced particle regeneration technique[J]. Journal of Computational Physics, 2022, 458: 111106. doi: 10.1016/j.jcp.2022.111106 [16] MENG Z F, ZHANG A M, WANG P P, et al. A shock-capturing scheme with a novel limiter for compressible flows solved by smoothed particle hydrodynamics[J]. Computer Methods in Applied Mechanics and Engineering, 2021, 386: 114082. doi: 10.1016/j.cma.2021.114082 [17] WANG P P, ZHANG A M, MENG Z F, et al. A new type of WENO scheme in SPH for compressible flows with discontinuities[J]. Computer Methods in Applied Mechanics and Engineering, 2021, 381: 113770. doi: 10.1016/j.cma.2021.113770 [18] SIROTKIN F V, YOH J J. A smoothed particle hydrodynamics method with approximate Riemann solvers for simulation of strong explosions[J]. Computers & Fluids, 2013, 88: 418-429. [19] 徐建于, 黄生洪. 圆柱形汇聚激波诱导Richtmyer-Meshkov不稳定的SPH模拟[J]. 力学学报, 2019, 51(4): 998-1011.XU Jianyu, HUANG Shenghong. Numerical simulation of cylindrical converging shock induced Richtmyer-Meshkov instability with SPH[J]. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(4): 998-1011. (in Chinese) [20] FULK D A, QUINN D W. An analysis of 1-D smoothed particle hydrodynamics kernels[J]. Journal of Computational Physics, 1996, 126(1): 165-180. doi: 10.1006/jcph.1996.0128 [21] LIU G R, LIU M B. Smoothed Particle Hydrodynamics: a Meshfree Particle Method[M]. Singapore: World Scientific Publishing, 2003. [22] SIGALOTTI L D G, LÓPEZ H, TRUJILLO L. An adaptive SPH method for strong shocks[J]. Journal of Computational Physics, 2009, 228(16): 5888-5907. doi: 10.1016/j.jcp.2009.04.041 [23] DANAILA I, JOLY P, KABER S M, POSTEL M. An Introduction to Scientific Computing[M]. New York: Springer-Verlag, 2007. [24] SOD G A. A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[J]. Journal of Computational Physics, 1978, 27(1): 1-31. doi: 10.1016/0021-9991(78)90023-2 [25] MONAGHAN J. SPH and Riemann solvers[J]. Journal of Computational Physics, 1997, 136(2): 298-307. doi: 10.1006/jcph.1997.5732 [26] WOODWARD P, COLELLA P. The numerical simulation of two-dimensional fluid flow with strong shocks[J]. Journal of Computational Physics, 1984, 54(1): 115-173. doi: 10.1016/0021-9991(84)90142-6 [27] TORO E F. Riemann Solvers and Numerical Methods for Fluid Dynamics[M]. Berlin: Springer-Verlag, 2009. -
下载:
计量
- 文章访问数: 211
- HTML全文浏览量: 88
- PDF下载量: 23
- 被引次数: 0
渝公网安备50010802005915号