2.2　光滑粒子动力学基本方程

使用SPH方法求解流体动力学问题一般分两步进行。第一步是积分表示,即使用积分表达式对函数进行近似,函数的积分表达式可以通过对核函数的影响区域积分进行近似。第二步是粒子表示,即通过对最近相邻粒子的值进行累加求和来近似离散点的函数值^[97]。

2.2.1　函数的积分表达

标准SPH是基于核估计发展起来的,在SPH方法中,对任意一个函数f,其积分近似表达式为

(2⁃1)

式中,<f(x)>为函数f(x)的近似值,x为位置矢量,Ω为包含x的积分体,W(x-x',h)为光滑函数,其依赖于两点之间的距离|x-x'|和光滑长度h。

光滑函数在SPH近似法中起着重要的作用,它决定了函数表达式的精度和计算效率。光滑函数W常常选用偶函数,其应满足3个条件。

第一个条件为正则化条件

(2⁃2)

由于光滑函数的积分值等于1,故此条件也称为归一化条件。

第二个条件是当光滑长度趋向于零时具有狄拉克函数性质

W(x-x',h)=δ(x-x')(2⁃3)

这里

δ(x-x')=(2⁃4)

第三个条件是紧支性条件

W(x-x',h)=0,当|x-x'|>κh 时(2⁃5)

式中,κ是与点x处光滑函数相关的常数,并确定光滑函数的有效范围(非零)。此有效范围称作点x处光滑函数的支持域。应用紧支性条件,在整个问题域内的积分即转换为在光滑函数的支持域内的积分。因此,一般来说积分域Ω就是支持域。

根据分部积分和散度定理,式(2⁃1)经过变换,对f(x)的导数的估计为

(2⁃6)

由上述方程得知,SPH估计将函数的空间导数转化为光滑函数的空间导数,因此,可以根据核函数求取任意场函数的空间导数。

2.2.2　函数的粒子表达

SPH方法最终要将函数积分近似表达式转化为支持域内所有粒子叠加求和的离散化形式,流体的问题域被离散为有限数量的粒子,其中每个粒子具有独立的质量、密度及其他物理属性。经过离散化,函数的积分表达式(2⁃1)可以写成如下的粒子近似式

(2⁃7)

式中,j和i表示粒子编号,m_j和ρ_j分别表示j粒子的质量和密度,N是粒子的总数,W_ij=W(x_i-x_j,h)=W(|x_i-x_j|,h)。

式(2⁃7)说明,在SPH方法中,对任意一个函数f,它在某一位置处的函数值可以通过应用光滑函数W对其光滑长度为h的紧支域内所有粒子插值求和的形式来表达^[98],如图2⁃1所示。

光滑长度h在SPH方法中非常重要,它在计算过程中决定了计算的实用性、有效性和可靠的适应性。光滑长度h直接影响到计算结果的效率和结果的精度。若h太小,则在以κh为半径的支持域内没有足够的粒子对所求粒子施加影响,这样会导致计算结果精度较低。如果h太大,太多的粒子对所求粒子产生影响,则可能将粒子所有细节信息和局部特性忽略,导致粒子趋于雷同,同样会影响到结果的精度。具体光滑长度应该选择多大与所研究的问题相关,没有一个广义的光滑长度。

式(2⁃6)经过粒子离散后,转化为

(2⁃8)

式中,Δ_iW_ij= =  ,r_ij=|x_i-x_j|。由于W存在紧支域,因此方程中的求和实际上只对i粒子影响区域内的粒子进行,而不是在整个计算域上,如图2⁃1所示。

2.2.3　光滑函数

无网格方法中的一个关键问题就是,如何在一系列不需要预先定义网格来提供节点连接的任意分布的离散点上有效地进行函数近似。函数近似的方法主要有三类:积分表示法、级数表示法和微分表示法。SPH方法通过光滑函数来进行积分表示。光滑函数非常重要,因为它不仅决定函数表达式的形式、定义粒子支持域的大小,而且还决定核近似和粒子近似的一致性、精度和计算效率。许多不同的光滑函数已被应用于SPH方法,许多不同的文献讨论了光滑函数的不同条件和性质。常用的光滑函数有三种,分别是高斯型光滑函数、三次样条光滑函数和高次样条函数。

2.2.3.1　高斯型光滑函数

Gingold和Monaghan最早使用了以下高斯型光滑函数来仿真非球形星体^[4],如图2⁃2所示。

W(R,h)=α_d(2⁃9)

式中,α_d在一维空间中为,在二维空间中为,在三维空间中为。

高斯型光滑函数是充分光滑的,即使对于高阶导数,也是一种较好的选择,这是因为它很稳定并且精度很高,特别是对于不规则粒子分布的情况。但是,它并不是真正严密的,因为理论上它是不可能为零的,除非R趋向于无穷大。但是由于它数值上趋于零的速度是很快的,所以实际上是紧密的。高斯型光滑函数计算量要求很大,因为光滑函数趋于零时要经历的距离比较大,特别是对于光滑函数的高阶导数。这样就会产生一个很大的支持区域,区域中用于粒子近似的粒子也会很多,导致运算时间较长。

Morris对高斯型光滑函数进行了修正^[99],修正的高斯光滑函数比原始的高斯型光滑函数更稳定、具有更高的效率,以下是二维的修正高斯光滑函数的表达式。

(2⁃10)

式中,δ为截断界限,在二维空间里面一般取δ=3h。

2.2.3.2　三次样条光滑函数

三次样条光滑函数又称B⁃样条光滑函数,由Monaghan和Lattanzio在三次样条函数的基础上首次使用^[100],如图2⁃3所示。

W(R,h)=α_d×(2⁃11)

在一维空间中α_d=,在二维空间中α_d=,在三维空间中α_d=。到目前为止,三次样条光滑函数是应用最广泛的光滑函数,这是因为在狭窄的紧支域中,三次样条光滑函数与高斯型光滑函数类似。但是,三次样条光滑函数的二阶导数是分段线性函数,相应地,其稳定性就会比那些光滑的核函数要差一些。另外,由于其光滑函数是分段的,所以其使用就会比只有一段的光滑函数要困难一些。

2.2.3.3　高次样条函数

Morris提出了更加近似于高斯型光滑函数并且更加稳定的高次样条函数(四次和五次)^[101],其中五次样条光滑函数如下所示

W(R,h)=α_d×(2⁃12)

式中,α_d的值在一维空间中为,在二维空间中为,在三维空间中为。