DOI: 10.12677/ijfd.2024.124009, PDF, HTML, XML,  被引量   
作者: 尤 一^*, 韩 豪, 李 晓, 吕世伟：海军潜艇学院基础部，山东 青岛
关键词: 光滑粒子流体动力学；简化有限差分法；顶盖驱动流；SPH； SFDI； Lid-Driven Flow

文章引用：尤一, 韩豪, 李晓, 吕世伟. 改进SPH方法在流动问题中的应用研究[J]. 流体动力学, 2024, 12(4): 88-96. https://doi.org/10.12677/ijfd.2024.124009

1. 引言

在船舶与海洋工程领域，经常可以观察到剧烈的自由表面流动现象，如液舱晃荡、甲板上浪和波浪砰击等。在这些过程中，会产生由流体与结构之间的相互作用而引起的整体或局部的冲击载荷，对这些冲击载荷的精确评估和计算对于船舶及海洋结构物的设计很重要。

目前，对于自由表面流动问题的数值研究主要是基于欧拉网格开展的，如FEM (有限元法) [1]和FVM (有限差分法) [2]。但是由于剧烈流动具有强非线性的复杂特性，使得传统的网格法在处理这方面问题时会存在很大的困难。这些困难主要在于如何通过网格方法精确有效地处理自由表面流动中运动界面的复杂性，当流动非常剧烈时，运动界面处会产生大变形甚至是非连续现象，使得界面的判别和捕捉会变得异常复杂，有时甚至会由于网格的大变形导致模拟的失败。因此近些年兴起的无网格计算方法逐渐引起了人们的关注。目前，光滑粒子流体动力学(Smoothed Particle Hydrodynamics, SPH)方法[3]作为一种无网格方法已被广泛应用于自由表面流动的模拟，如多相流、液舱晃荡、波浪冲击和流固耦合问题。SPH方法的拉格朗日特性使得其在自由液面流动的模拟过程中可以不受网格畸变的限制，即使是一些表面波的破碎、飞溅等现象也可以进行捕捉和预报。

尽管传统SPH方法已经在很多领域得到应用，但由于粒子的拉格朗日特性，容易在流场中产生非物理压力振荡[4]。在传统SPH的基础上，研究者们提出了很多改进的SPH方法以提高计算的准确性和稳定性，如黎曼SPH [5]和δ-SPH [6]。δ-SPH方法是目前应用比较广泛而且效果很好的方法，它是通过在连续性方程中添加适当的密度数值扩散项来减小压力场中的高频数值振荡。近年来，许多研究都证明了δ-SPH方案在模拟剧烈流动方面具有非常好的可靠性[7]。

尽管SPH方法近些年来取得了很大的进展，但仍存在很多挑战[8]，一个重要原因就是SPH方法本身所使用的离散格式的精度不足，尤其是当自由表面粒子的支持域被截断时，使得核近似和粒子近似的精度不能使守恒方程和动态边界条件精确满足，这些误差会使得计算的密度场产生波动，即便是很小的波动，在经过状态方程(EOS)的放大后也会形成比较大的压力场噪声，又反作用于粒子产生位移上的扰动，使得粒子分布变得杂乱，因此需要进一步改进插值方法来提高计算精度。近些年研究者们提出了很多高阶的插值方法应用于SPH方法的计算中，比较常见有最小二乘法和一阶相容核近似法，虽然这些方法都很好地起到了提升插值精度的作用，但是由于涉及到的矩阵运算量很大，使得计算效率并不理想，而且在自由液面的破碎点或飞溅处容易产生奇异矩阵，导致计算的失败，因此找到一种可以满足计算精度且对计算效率影响不大的插值方法应用于SPH方法中是目前SPH方法研究中非常重要的一项内容。本文将一种应用于MLPG_R法[9]中的简化有限差分法(SFDI)应用于δ-SPH方法中，来研究其对于流动问题计算精度的改进效果。

2. 计算原理

2.1. δ-SPH方法

传统的SPH方法主要是基于弱可压缩的假设，正压流体的控制方程主要是连续性方程和动量方程：

$\frac{D\rho }{Dt}+\rho \nabla \cdot u=0$ (1)

$\frac{Du}{Dt}=-\frac{1}{\rho }\nabla P+f$ (2)

其中： $\rho$ 是流体的密度，t是时间，u是速度矢量，P是压强，g是重力加速度，f是体积力。

在SPH方法中，一个任意函数 $f\left(x\right)$ 的积分可以表示为：

$f\left(x\right)={\displaystyle {\int }_{\Omega }f\left(x^{\prime }\right)\delta \left(x-x^{\prime }\right)\text{d}{x}^{\prime }}$ (3)

其中：Ω是积分域， $\delta \left(x-x^{\prime }\right)$ 是狄拉克函数，其形式如下：

$\delta \left(x-x^{\prime }\right)=\{\begin{array}{l}1,x=x^{\prime }\\ 0,x\ne x^{\prime }\end{array}$ (4)

在SPH方法中，狄拉克函数使用核函数 $W\left(x-x^{\prime },h\right)$ 来代替，其中h被称为光滑长度。于是函数的核近似公式可以写为：

$\langle f\left(x\right)\rangle ={\displaystyle {\int }_{\Omega }f\left(x^{\prime }\right)W\left(x-x^{\prime },h\right)\text{d}{x}^{\prime }}$ (5)

因此，目标i粒子的函数 $f\left(x\right)$ 基于离散粒子的近似可以写成：

$\langle f\left(x_i\right)\rangle ={\displaystyle {\sum }_{j=1}^N\frac{m_j}{{\rho }_j}f\left(x_j\right)W\left(x_i-x_j,h\right)}$ (6)

其中：目标j是i粒子支持域内的粒子，N是i粒子支持域内的粒子总数， $m_j$ 是j粒子的质量， ${\rho }_j$ 是j粒子的密度。

根据上述的核近似和粒子近似，i粒子的质量守恒方程和动量方程可以表示为：

$\frac{D{\rho }_i}{Dt}={\rho }_i{\displaystyle {\sum }_{j=1}^N\frac{m_j}{{\rho }_j}{u}_{ij}\cdot {\nabla }_i{W}_{ij}}+\delta h{C}_0{\displaystyle {\sum }_{j=1}^N\frac{m_j}{{\rho }_j}{\Psi }_{ij}\cdot {\nabla }_i{W}_{ij}}$ (7)

$\frac{D{u}_i}{Dt}=g-{\displaystyle {\sum }_{j=1}^N{m}_j\left(\frac{P_i}{{\rho }_i^2}+\frac{P_j}{{\rho }_j^2}+{\Pi }_{ij}\right){\nabla }_i{W}_{ij}}$ (8)

${\Psi }_{ij}=2\left({\rho }_j-{\rho }_i\right)\frac{r_{ji}}{{\left|r_{ij}\right|}^2}-\left(\nabla {\rho }_i^L+\nabla {\rho }_j^L\right)$ (9)

${\langle \nabla \rho \rangle }_i^L={\displaystyle {\sum }_{j=1}^N\frac{m_j}{{\rho }_j}\left({\rho }_j-{\rho }_i\right)L_i\cdot {\nabla }_i{W}_{ij}}$ (10)

$L_i={\left[{\displaystyle {\sum }_{j=1}^N\frac{m_j}{{\rho }_j}{r}_{ji}\otimes {\nabla }_i{W}_{ij}}\right]}^{-1}$ (11)

${\Pi }_{ij}=\{\begin{array}{ll}\frac{-{\alpha }_{\Pi }{c}_{ij}{\varphi }_{ij}+{\beta }_{\Pi }{\varphi }_{ij}^2}{{\rho }_{ij}},\hfill & u_{ij}\cdot r_{ij}<0\hfill \\ 0,\hfill & u_{ij}\cdot r_{ij}\ge 0\hfill \end{array}$ (12)

${\varphi }_{ij}=\frac{h_{ij}{u}_{ij}\cdot r_{ij}}{{\left|r_{ij}\right|}^2+0.01h_{ij}^2}$ (13)

其中系数δ控制着密度扩散的强度，一般取0.1， $C_0$ 是人工声速，一般取值为10倍的最大速度，r是位置矢量， ${\Pi }_{ij}$ 是人工粘性项， ${\alpha }_{\Pi }$ 和 ${\beta }_{\Pi }$ 是人工粘性系数。在公式(10)中， ${\langle \nabla \rho \rangle }^L$ 是被矫正矩阵L矫正后的重整化密度梯度。

粒子的压强是通过状态方程(EOS)来进行求解的：

$P_i=C_0^2\left({\rho }_i-{\rho }_0\right)$ (14)

其中 ${\rho }_0$ 是初始流体粒子的密度。

2.2. 简化有限差分法方案

在SPH方法中，求解控制方程需要用到散度项和梯度项，但是上文中提到的粒子近似的误差会对计算散度和梯度时的精确度造成比较大的影响，从而影响密度场和压力场的计算，因此改进粒子近似的插值方式对于计算精度的提升具有很大的作用。

Ma在MLPG_R法中应用的简化有限差分法(SFDI)是一种基于泰勒级数展开的二阶精度格式[9]，具体的形式为：

${\left(\frac{\partial f}{\partial x}\right)}_{r_i}={\displaystyle \underset{j=1,j\ne i}{\overset{N}{\sum }}\frac{n_{i,x_m}{B}_{ij,x_k}-n_{i,xy}{B}_{ij,y}}{n_{i,x_1}{n}_{i,x_2}-n_{i,xy}^2}\left(f_j-f_i\right)}$ (15)

${\left(\frac{\partial f}{\partial y}\right)}_{r_i}={\displaystyle \underset{j=1,j\ne i}{\overset{N}{\sum }}\frac{n_{i,x_1}{B}_{ij,y}-n_{i,xy}{B}_{ij,x}}{n_{i,x_1}{n}_{i,x_2}-n_{i,xy}^2}\left(f_j-f_i\right)}$ (16)

$n_{i,xy}={\displaystyle \underset{j=1,j\ne i}{\overset{N}{\sum }}\frac{\left(r_{j,x_m}-r_{i,x_m}\right)\left(r_{j,x_k}-r_{i,x_k}\right)}{{|r_j-r_i|}^2}W\left(|r_j-r_i|\right)}$ (17)

$B_{ij,x_m}=\frac{r_{j,x_m}-r_{i,x_m}}{{|r_j-r_i|}^2}W\left(|r_j-r_i|\right)$ (18)

3. 简化有限差分法的收敛性和精度分析

在本节中将会通过算例对SFDI方案在SPH方法中的收敛性和精度进行测试。

图1中给出了计算域的设置：边长为1 m的方形域，如图所示初始粒子均匀布置，测试函数的形式为 $f\left(x,y\right)={\text{e}}^{x+2y}$ ，分别使用普通插值和简化有限差分法来计算 $f\left(x,y\right)$ 对x的一阶偏导数，将计算结果与解析解进行对比，分析两者的收敛性和误差。计算中两种方案均分别使用0.02 m、0.01 m、0.005 m和0.0025 m的粒子尺寸，对应的计算域粒子总数为2500、10,000、40,000和160,000。

定义数值计算的平均误差为

${\text{Error}}_x=\frac{1}{N}{\displaystyle \underset{j=1}{\overset{N}{\sum }}\left|\frac{{\tilde{f}}_x\left(r_j\right)-f_x\left(r_j\right)}{{\tilde{f}}_x\left(r_j\right)}\right|}$ (19)

其中：N为计算域中的粒子数，i为目标粒子，j为i粒子的邻域粒子， ${\tilde{f}}_x\left(r_j\right)$ 为解析解。

表1中给出了使用两种插值方法计算的一阶偏导数在不同粒子尺寸下的计算误差对比，从对比中可以看出通过SFDI的应用，使得偏导数计算的误差大大减小。图2中给出了两种方法的计算误差随粒子尺寸的收敛性分析，从图中的对比可以看出，传统的偏导数模型收敛精度只能够达到一阶精度，而应用了SFDI方案后，精度可以达到二阶，误差大大降低，接下来将在流动问题的模拟中进一步对改进的方法进行测试和验证。

Figure 1. Particle distribution in the computational domain

图1. 计算域的粒子分布

Table 1. Errors of first-order partial derivatives calculated by two interpolation schemes using different particle sizes

表1. 在不同粒子尺寸下的两种插值方法计算一阶偏导数的计算误差

粒子尺寸(m)	传统方法的平均误差(100%)	SFDI的平均误差(100%)
0.02	0.05326	0.00530
0.01	0.02723	0.00138
0.005	0.01418	0.00035
0.0025	0.00757	0.00009

Figure 2. Comparisons for error convergence of first order partial derivative by using two schemes

图2. 两种方案计算的一阶偏导数误差的收敛性对比

4. 流动问题中的模型应用

在本节中将会通过一个二维的方腔顶盖驱动流算例来对改进的SPH方案进行测试。图3中给出了方腔的示意图，如图所示：方腔的边长为1 m，两个侧壁和底部都保持静止，方腔顶部以U = 1 m/s的速度驱动流体向右运动，粒子尺寸设置为0.0025 m，全域布置160,000粒子。

Figure 3. Schematic view of the square cavity

图3. 方腔的示意图

图4给出了时间t = 40s时计算域内的速度幅值的分布，此时流动状态已经处于稳定，从图中可以看出两种计算方案中速度的分布大体相同，但是如果对局部进行放大后可以看出传统的SPH方法经过40秒后所计算的粒子分布杂乱，出现了聚集现象，扰动明显，这是因为SPH方法本身的拉格朗日特性所引起的，由于粒子与粒子之间没有连接，是自由运动的，所以会有逐渐混乱的趋势，但是因为传统方法的计算精度不高，使得粒子的分布很快就出现了混乱的现象，而改进后的SPH方法由于使用了高精度的插值方案，提升了一阶偏导数的计算精度，使得连续性方程和动量方程的计算精度都得到提升，计算出的粒子图中粒子的分布均匀度有了明显的改善。

Figure 4. Comparisons of snapshots of the lid-driven flow. (a) Results of SPH; (b) Results of SPH with SFDI scheme

图4. 方腔顶盖驱动流的粒子对比图。(a) SPH结果；(b) SPH + SFDI结果

图5中给出了时间t = 40 s时，中纵线处x方向上的速度分量曲线对比图，图6中给出了时间t = 40s 时，中横线处y方向上的速度分量曲线对比图，分别为两种SPH方法的计算结果与实验数据[10]的对比。从对比结果中可以看出经过SFDI方案改进后的SPH结果与实验值有着更好的一致性，计算结果更加精确，这也证明了改进后的SPH方法在流动问题的模拟中具有更好的精度。

Figure 5. Comparisons of x-direction component of velocity at the midline (t = 40 s)

图5. 中纵线处x方向上的速度分量曲线对比(t = 40 s)

Figure 6. Comparisons of y-direction component of velocity at the middle horizontal line (t = 40 s)

图6. 中横线处y方向上的速度分量曲线对比(t = 40 s)

图7和图8分别给出了改进SPH方法在不同粒子尺寸(0.01 m、0.005 m和0.0025 m)下的计算结果(中纵线处x方向上的速度分量曲线和中横线处y方向上的速度分量曲线)，从对比中可以看出采用改进后的SPH方法计算的结果随着粒子尺寸的减小更加趋近于实验值，体现出了很好的收敛性。

Figure 7. Comparisons of x-direction component of velocity at the midline calculated by the improved SPH method using different particle sizes

图7. 不同粒子尺寸下改进SPH方法计算的中纵线处x方向上的速度分量曲线对比

Figure 8. Comparisons of y-direction component of velocity at the middle horizontal line calculated by the improved SPH method using different particle sizes

图8. 不同粒子尺寸下改进SPH方法计算的中横线处y方向上的速度分量曲线对比

5. 结论

本文采用一种改进的SPH方法应用于流动问题的模拟中，在对比了传统插值方法与新插值方法对于一阶偏导数的求解效果后发现改进后的插值方案能够取得误差更小的一阶偏导数结果，因此对于SPH方法中质量守恒方程和动量方程等控制方程中梯度项和散度项的求解也具有更大的优势。在此基础上，将改进后的SPH方法应用到流体问题的计算中，文中选取了经典案例顶盖驱动流作为验证案例，对比了传统SPH方法和改进后的SPH方法的计算结果，进一步验证了新算法的鲁棒性和精确性，证明了新算法在流动问题的模拟方面具有巨大的潜力。本文经过研究得出的结论如下：简化有限差分法可以应用于SPH方法中，并起到提升一阶偏导数计算精度的作用，通过得到更高精度的一阶偏导数可以提升连续性方程和动量方程的计算精度，从而提升SPH方法在流动问题模拟中的鲁棒性和精确性。

NOTES

^*通讯作者。

参考文献

[1]	Wu, G.X., Ma, Q.W. and Eatock Taylor, R. (1998) Numerical Simulation of Sloshing Waves in a 3D Tank Based on a Finite Element Method. Applied Ocean Research, 20, 337-355. https://doi.org/10.1016/s0141-1187(98)00030-3
[2]	Chen, B. and Nokes, R. (2005) Time-Independent Finite Difference Analysis of Fully Non-Linear and Viscous Fluid Sloshing in a Rectangular Tank. Journal of Computational Physics, 209, 47-81. https://doi.org/10.1016/j.jcp.2005.03.006
[3]	Monaghan, J.J. (1994) Simulating Free Surface Flows with SPH. Journal of Computational Physics, 110, 399-406. https://doi.org/10.1006/jcph.1994.1034
[4]	Wang, L., Khayyer, A., Gotoh, H., Jiang, Q. and Zhang, C. (2019) Enhancement of Pressure Calculation in Projection-Based Particle Methods by Incorporation of Background Mesh Scheme. Applied Ocean Research, 86, 320-339. https://doi.org/10.1016/j.apor.2019.01.017
[5]	Inutsuka, S. (2002) Reformulation of Smoothed Particle Hydrodynamics with Riemann Solver. Journal of Computational Physics, 179, 238-267. https://doi.org/10.1006/jcph.2002.7053
[6]	Antuono, M., Colagrossi, A. and Marrone, S. (2012) Numerical Diffusive Terms in Weakly-Compressible SPH Schemes. Computer Physics Communications, 183, 2570-2580. https://doi.org/10.1016/j.cpc.2012.07.006
[7]	Marrone, S., Antuono, M., Colagrossi, A., Colicchio, G., Le Touzé, D. and Graziani, G. (2011) δ-SPH Model for Simulating Violent Impact Flows. Computer Methods in Applied Mechanics and Engineering, 200, 1526-1542. https://doi.org/10.1016/j.cma.2010.12.016
[8]	Gotoh, H. and Khayyer, A. (2018) On the State-of-the-Art of Particle Methods for Coastal and Ocean Engineering. Coastal Engineering Journal, 60, 79-103. https://doi.org/10.1080/21664250.2018.1436243
[9]	Ma, Q.W. (2008) A New Meshless Interpolation Scheme for MLPG_R Method. CMES: Computer Modeling in Engineering & Sciences, 23, 75-89.
[10]	Ghia, U., Ghia, K.N. and Shin, C.T. (1982) High-Resolutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method. Journal of Computational Physics, 48, 387-411. https://doi.org/10.1016/0021-9991(82)90058-4