SPH法理论基础

1、摘要1第一章 绪论21.1 基于网格的方法21.1.1 拉格朗日网格31.1.2 欧拉网格31.1.3 基于网格的数值方法的局限性31.2 无网格方法41.3光滑粒子动力学方法4第二章 SPH方法在流体动力学问题中的应用52.1 光滑粒子动力学原理52.2 光滑粒子动力学基本方程52.2.1 函数的积分表达52.2.2 函数的粒子表达72.2.3光滑函数72.3 拉格朗日型的Navier-Stokes 方程72.4 Navier-Stokes方程的SPH表达式8第三章 SPH的计算实施方法93.1 粒子的密度近似93.2 核函数103.3 状态方程113.4 人工粘度123.5 边界处理133

2、.6 时间积分13摘要计算机数值仿真逐渐成为解决现代工程和科学问题的一条重要途径。数值仿真能为理论提供测试和检验，有助于对复杂的物理问题加深认识，甚至还能帮助解释和发现新现象。基于网格的数值方法虽然已经有广泛的应用，但是在很多方面仍存在不足之处，比如在计算流体动力学的大变形、运动物质交界面、自由表面等问题时，由于网格产生畸变导致计算误差过大或无法进行，从而使其在许多问题的应用上受到限制。近年来，无网格法倍受关注，这种方法在许多应用中都优于传统的基于网格的有限元法、有限差分法以及有限体积法等数值方法。本文主要研究新一代无网格方法-光滑粒子流体动力学方法（SPH）。第一章 绪论1.1 基于网格的方

3、法数值计算方法通常可以分为两种：基于网格的方法和无网格方法。通常对于物理控制方程的描述有两种基本方法：欧拉描述法和拉格朗日描述法。欧拉描述法是对空间的描述方法，其典型代表是有限差分法；拉格朗日描述法是对物质点的描述方法，其典型代表是有限元法。欧拉描述和拉格朗日描述对应着两种不同的区域离散化网格：欧拉网格和拉格朗日网格。针对不同类型的问题，这两种网格在数值方法中都得到广泛应用。1.1.1 拉格朗日网格拉格朗日网格，基于拉格朗日网格的数值方法在整个计算过程中网格是固定附着于物质上的，网格会随着物质的运动而运动，所以在物质点上的所有场变量的整个时间历程都可以很容易地追踪，常见的方法如有限元法等。基于

4、拉格朗日网格方法的优点是：由于在相关的偏微分方程里不存在迁移项，所以程序在方案设计上会变得相对简单而且运行较快；由于只需在问题域内布置网格，问题域外不需要布置，所以计算效率很高；不规则或者复杂的几何形状可以用不规则的网格来处理。由于具有以上这些优点，拉格朗日方法得到广泛的应用，并且能成功地求解计算固体力学问题。然后，基于拉格朗日网格的方法难以应用于具有极大网格变形的情况，因为其公式的形式是以网格为基础的，当网格变形太大的时候，公式的精度和求解都会受到很大的影响。另外，由于时间步长是由最小单元尺寸所控制，若网格太小就会影响计算的效率，甚至会导致计算失败。1.1.2 欧拉网格欧拉网格，相对于拉格朗

5、日网格，欧拉网格刚好相反，它是固定在模拟对象所处的空间上，模拟对象在固定网格单元上运动。因此，在物质流过网格时，所有网格节点和网格单元依然固定在空间上而且不会随时间的改变而改变。通过模拟质量、动量和能量经过网格单元边界的通量，可以计算质量、速度和能量等等的分布。在整个计算过程中网格单元的形状和体积都保持不变。由于欧拉网格在时间和空间上都是固定的，物体的大变形不会引起网格本身的任何变化，所以再以物质流动为主体的计算流体力学领域里，较为广泛的应用欧拉法。但是，欧拉法依然有很多缺点。由于欧拉法不能用固定网格来追踪物质的运动，所以很难分析物质上固定点的场变量的变化情况而只能得到空间上固定的欧拉网格的场

6、变量的变化情况；由于欧拉法追踪的是流过网格单元边界的质量、动量和能量的通量，所以自由表面、变形边界和运动物质交界面的位置就很难精确确定；由于欧拉法需要在整个计算区域上都覆盖网格，所以有时为了提高计算效率而使用较粗糙的网格，这样就会降低离散化区域的求解精度。1.1.3 基于网格的数值方法的局限性传统的基于网格的数值方法如有限差分法和有限元法已经广泛地应用在计算流体力学和计算固体力学的各个领域中，是现在进行区域离散化和数值离散化模拟的主要方法。但是基于网格的数值方法在很多方面仍存在不足之处，比如在计算流体动力学中的大变形、运动物质交界面、自由表面等问题时，由于网格产生畸变导致计算误差过大或无法进行

7、，使其在许多问题的应用上受到限制。在基于网格的数值方法中，数值模拟的先决条件就是在问题域生成网格。对于基于欧拉网格的方法，在固定的欧拉网格上要精确确定自由表面、变形边界、运动交界面和不均匀物质之间的位置是一项非常困难的工作。并且欧拉方法也不适用于研究如粒子流动这类问题。即需要在固定的物质体内监控材料特性的问题。对于拉格朗日网格法如有限元法，进行模拟前必须要在研究对象上建立网格，这项操作常常占用很大的计算工作量。当所研究对象是一系列离散的物质点时，同样不适合使用基于网格的方法，用连续的基于网格的方法对这些离散系统进行模拟通常不是好的选择。在计算大变形问题时，比如高速碰撞、金属加工成型、动态裂纹扩

8、展、流固耦合和应变局部化等，基于网格的方法遇到的困难更大。原因是发生大变形时，网格畸变过大，使得计算中断，有限元方法通常采用单元“销蚀”法或重分网格技术来克服这种困难。但是单元“销蚀”法本身缺乏物理依据，纯粹是为了使计算进行下去的一种数值手段。而网格重分技术在节点重新分配物理量时，很难保证系统动量、能量守恒，因而导致计算的精度下降。此外，网格重分技术不是很容易实现，为了更好地解决大变形问题，必须对网格有新的处理方法或去除网格，所以各种无网格方法相继被提出来。1.2 无网格方法无网格方法产生于三十多年以前，当时发展很慢，直到近几年，由于无网格方法的近似函数不依赖于网格，在分析涉及大变形的问题中被

9、认为优于传统的基于网格的有限差分法和有限元法，才受到了众多研究者的高度重视，成为国际计算仿真界的研究热点之一，并得到了迅速的发展【1】。无网格方法的主要思想是：通过使用一系列任意分布的节点（或粒子）来求解具有各种边界条件的积分方程或偏微分方程组，从而得到精确稳定的数值解，这些节点或粒子之间不需要网格进行连接。由于无网格法采用基于点的近似，可以彻底或部分地消除网格，不需要网格的初始划分和重构，因此不仅可以保证计算的精度，而且可以减小计算的难度。经过几十年的研究，目前已发展了十余种无网格方法。最早提出的无网格方法是光滑粒子动力学方法（Smoothed Particle Hydrodynamics，

10、简称SPH），它是由Lucy、Gingold和Monaghan在1977年分别提出的，并且在天体领域得到成功地应用【2,3】。但是由于精度和稳定性问题，该方法最初并未得到广泛的应用。80年代，Monaghan等人在该方法的研究与应用中作出突出贡献【4,5】，将光滑粒子动力学方法应用到连续固体力学和流体力学中，他们将SPH方法解释为核函数法，模拟了流场中的激波强间断现象。90年代，Swegle、Dyka和Chen等提出了SPH方法不稳定的原因及稳定化方法【6-8】。Johnson和Beissel等人也提出了一些用来改善应变计算的方法【9】。随着研究的深入，SPH方法被应用于水下爆炸数值仿真【10

11、】、高速碰撞中材料动态响应数值仿真等领域【11-13】。近几年来，我国的学者也开始关注SPH计算方法，中科院的张锁春对SPH方法进行了综述【14】，国防科学技术大学的邓方刚研究了SPH在柱坐标体系下的应用【15】。国防科大的贝新源、岳宗五等将SPH方法用于高速碰撞问题研究【16】。1.3光滑粒子动力学方法光滑粒子流体动力学（SPH-Smoothed particle hydrodynamics）方法是近30年发展起来的一种新的纯拉格朗日无网格粒子法，它最初是用于解决三维开放空间的天体物理学问题【2,3】，现已被广泛地研究和扩展，并被应用于具有材料强度的动态响应问题和具有大变形的流体动力学问题。

12、SPH方法和它的派生方法是现在粒子法的主要类型。SPH方法通过对邻近粒子进行加权平均而得到稳定的光滑近似性质，该方法应用在流体动力学问题的范围内。由于SPH方法中的自适应性、粒子性质和拉格朗日性质的和谐结合，使其在工程和科学领域都得到了实际应用。国外对SPH方法研究的比较早，已经被应用到很多领域中。SPH方法最初在天体研究中应用较多，如双子星和行星碰撞、超新星、星系的形成和瓦解、黑洞和中子星合并、白矮星的单个或多重爆炸甚至是宇宙的进化等问题的研究。第二章 SPH方法在流体动力学问题中的应用2.1 光滑粒子动力学原理流体动力学问题的求解主要是解基于密度、速度、能量等变量场的偏微分方程组。除了一些

13、简单问题可以得到偏微分方程组的解析解以外，大部分的流体动力学问题只能寻求其数值解。光滑粒子动力学的原理如下：1） 使用SPH方法求解流体动力学问题时，问题域用一系列任意分布的粒子来表示，粒子之间没有任何连接，因此SPH方法是无网格性质的，通常为了求解方便，初始粒子都均匀排列；2） 场函数用积分表示法来近似，在SPH方法中称为核近似法；3） 然后应用支持域内的相邻粒子对应的值叠加求和取代场函数的积分表达式来对场函数进行粒子近似，由于在每一个时间步内都要进行粒子近似，支持域内的有效粒子为当前时刻支持域内的粒子，因此SPH方法具有自适应性；4） 将粒子近似法应用于所有偏微分方程组的场函数相关项中，将

14、偏微分方程组进行离散，可知SPH是一种纯拉格朗日方法；5） 粒子被附上质量后，则意味着这些粒子是真实的具有材料特性的粒子；最后应用显式积分法得到所有粒子的场变量随时间的变化值。从以上分析可以看到，SPH方法是一种纯拉格朗日的具有无网格、自适应属性的流体动力学求解方法。2.2 光滑粒子动力学基本方程使用SPH方法求解流体动力学问题一般分两步进行。第一步是积分表示，即使用积分表达式对函数进行近似，函数的积分表达式可以通过对核函数的影响区域积分进行近似。第二步是粒子表示，即通过对最近相邻粒子的值进行累加求和来近似离散点的函数值。2.2.1 函数的积分表达标准SPH是基于核估计发展起来的，在SPH方法

15、中，对任意一个函数f，其积分近似表达式为其中，<f(x)>为函数f(x)的近似值，x为位置矢量，为包含x的几分体，W(x-x,h)为光滑函数，其依赖于两点之间的距离|x-x和光滑长度h。光滑函数在SPH近似法中其着重要的作用，它决定了函数表达式的精度和计算效率。光滑函数W常常选用偶函数，其应满足3个条件。第一个条件为正则化条件，如下所示由于光滑函数的积分值等于1，故此条件也称为归一化条件。第二个条件是当光滑长度趋向于零时具有狄拉克函数性质这里，第三个条件是紧支性条件式中，K是与点x处光滑函数相关的常数，并确定光滑函数的有效范围（非零）。此有效范围称作点x处光滑函数的支持域内的积分。

16、因此，一般来说积分域就是支持域。根据分布积分和散度定理，将第一个式子经过变换，对f（x）的导数的估计为由以上方程得知，SPH估计将函数的空间导数转化为光滑函数的空间导数，因此，可以根据核函数求取任意场函数的空间导数。2.2.2 函数的粒子表达SPH方法最终要将函数积分近似表达式转化为支持域内所有粒子叠加求和的离散化形式，流体的问题域被离散化为有限数量的粒子，其中每个粒子具有独立的质量、密度及其它物理属性。经过离散化，函数的积分表达式可以写成如下的粒子近似式其中，j和i表示粒子编号，mj和j分别j粒子的质量和密度，N是粒子的总数，Wij=W(xi-xj，h)=W(|xi-xj |，h)。上式说明

17、，在SPH方法中，对任意一个函数f，它在某一位置处的函数值可以通过应用光滑函数W对其光滑长度为h的紧支域内所有粒子插值求和的形式来表达【17】。光滑长度h在SPH方法中非常重要，它在计算过程中决定了计算的实用性、有效性和可靠的适应性。光滑长度h直接影响到计算结果的效率和结果的精度。若h太小，则在以kh为半径的支持域内没有足够的粒子对所求粒子施加影响，这样会导致计算结果精度较低。如果h太大，太多的粒子对所求粒子产生影响，则可能将粒子所有细节信息和局部特性忽略，导致粒子趋于雷同，同样会影响到结果的精度。具体光滑长度应该选择多大与所研究的问题相关，没有一个广义的光滑长度。2.2.3光滑函数SPH通过

18、光滑函数来进行积分表示。光滑函数非常重要，因为它不仅决定函数表达式的形式，定义粒子支持域的大小，而且还决定核近似和粒子近似的一致性和计算效率。常用的光滑函数有三种，分别是高斯型光滑函数，三次样条光滑函数和高次样条函数。可参见参考文献【3,18,19,】2.3 拉格朗日型的Navier-Stokes 方程对控制方程的描述有两种形式：欧拉描述法和拉格朗日描述法。欧拉描述法是对空间点的描述，而拉格朗日描述法是对物质点的描述。流体的流动过程可以用欧拉方程以及适当的状态方程来表示，这些方程表示了质量、动量和能量守恒。如果用上标和来表示坐标方向和求和重复指标，采用时间全导数形式，N-S方程组可以表示为：连

19、续性方程：动量方程：能量方程：在以上方程中，为总应力张量，有两部分组成：一部分是压力p，另外一部分是粘性应力。对于牛顿流体，粘性剪应力与剪应变率成比例，且比例系数为粘性系数。其中如果将压力与粘性应力分开写，我们可以得到如下的能量方程在上方式子中，各字母代表的意义如下：V -速度向量E-内能-密度p-压力t-时间2.4 Navier-Stokes方程的SPH表达式根据光滑函数的概念和SPH方法的粒子近似原理，流体流动的运动方程可以写成如下的离散SPH形式：其中 i-参考粒子。 J-粒子i的相邻粒子。第三章 SPH的计算实施方法3.1 粒子的密度近似由于在SPH方法中粒子的运动是由相邻粒子之间的密

20、度差以及压力差产生的，所以密度近似在SPH方法中非常重要。在一般的SPH方法中主要有两种方式对密度进行展开：一种是密度求和法，另一种是连续密度法。在SPH方法中粒子的运动是由相邻粒子之间的密度差以及压力差产生的，所以密度近似在SPH方法中非常重要。在一般的SPH方法中主要由两种方式对密度进行展开：一种是密度求和法，另一种是连续密度法。在SPH方法的实际应用中，密度求和法似乎使用得更为广泛，主要是因为密度求和法体现了SPH近似法的本质。现在已经有不少改进此种算法精度的修正方案。其中一种有效地改进是将邻近粒子的SPH光滑函数累加式，即上式的等号右端正则化。上式可提高自有边界处和密度不连续材料交界面

21、处的精度，此式也非常适合模拟非连续广义流体的流动问题。密度求和法的缺陷是所需要的计算量大，因为必须在计算其他参数前计算密度，而且还要计算自身的光滑函数。近似密度的另一种粒子近似法为连续性密度法，这种近似法通过应用SPH近似法的概念对连续性方程进行转换而得。常用的连续性密度方程为：连续性密度近似法不需要先计算密度，因此可以提高计算效率，并适于并行处理算法的编程。具体应用哪一种密度算法要根据研究的问题而定，对于广义流体问题的模拟，应用修正的密度求和法可得到较好的结果。对于具有强间断问题的模拟，如爆炸、高速冲击等，应优先选取连续性密度法。3.2 核函数核函数在SPH近似法中起着重要作用，因为它控制了

22、插值的方式和粒子影响域的效能，从而决定了函数表达式的精度和计算效率。由于SPH插值基于积分插值理论，用核函数来近似狄拉克函数，因此核函数应该具有某些特性，如非负性，紧支性，归一性，衰减性和狄拉克函数条件。在实际应用中，通常使用两种形式的核函数，一种是如下所示的高斯核函数。公式中的d在一维空间中为，在二维空间中为，在三维空间中为。另一种是三次样条核函数，到目前为止，三次样条函数是最广泛应用的核函数。在一维空间中，在二维空间中，在本文的三维空间中，。3.3 状态方程在求解可压缩流体问题的标准SPH方法中，粒子的运动是由于压力梯度的作用而产生的，而粒子的压力是通过状态方程由粒子自身的密度和内能来计算

23、的。然而，在不可压缩流体问题中，流体实际的状态方程限制了时间步长的大小，即时间步长不能太小。有效地计算动量方程中的压力项是仿真计算不可压缩流体的一个主要任务。事实上，理论不可压缩流体实际上是可压缩的，所以可以引进人工压缩率。人工压缩率的引进就是把所有理论不可压缩流都考虑为实际上是可压缩的。因此，可用仿真可压缩状态方程去仿真不可压缩流。不同的物质应该具有不同的状态方程。在本文中应用较多的是如下所示的状态方程:式中：是常数，在许多情况下取=7；0是参照密度；B是由具体问题而定的参数，用于限制密度的最大改变量。在许多情况下，一般用B作为初始压力。上式中减去1能消除自由表面流动的边缘效应。由上式可以看

24、出，密度的微小变化都会导致压力值产生较大的波动。另外一个可选择的人工状态方程为式中：c为声速。1997年，Morris使用此人工状态方程通过SPH方法数值仿真了低雷诺数不可压缩流动【20】。Zhu等使用此人工状态方程通过SPH方法用于微孔尺度的数值模型，用于数值仿真穿越多孔介质的流动。在人工压缩率技术中，声速是一个需要慎重考虑的因素。如果使用实际的声速，则意味着用理想不可压缩的人工流体来近似真实的流体。相对密度差为：式中：Vb和M分别是流体的整体速度和马赫数。由于真实的声速相当大，所以相应得到的马赫数非常小，相对密度差几乎可以忽略。因此为了用人工压缩流体近似真实流体，应该使用比真实值要小很多的

25、声速。所以对声速大小的要求为：一方面，必须足够大，以至于人工可压缩流体的特性与真实流体充分接近；另一方面，应足够小，可使时间步的增量在容许范围内。3.4 人工粘度对于某些数值问题，误差跳跃很大，会产生数值振荡，为了提高算法的稳定性，并且防止粒子间相互接近时的非物理穿透，Monaghan和Gingold在SPH方法的动量方程中引入了人工粘度【21】，它不是物理的粘度，只是数值的修正。式中，在式中，为常数，取值一般为1.0左右【22】。因子=0.1hij用于防止粒子相互靠近时产生的数值发散。c和v分别表示声速和粒子的速度矢量。将人工粘度引入到动量方程中，最终动量方程的SPH粒子近似式可写为如下形式

26、：3.5 边界处理在SPH方法的应用中，边界条件的处理既是该方法的优点，也是目前的薄弱环节，因为SPH方法最初是应用在没有边界条件限制的领域内，后来逐渐扩展到一些具有边界的工业应用中。我们使用虚粒子来处理固定边界条件，通过设置边界虚粒子可以方便地模拟固体边界，同时也要防止粒子的非物理穿透。根据不同的数值仿真条件，我们在SPH数值模拟中一共使用了两种类型的虚粒子第一种类型的虚粒子设置在固定边界上，与Monaghan所使用的相似【4】。二种类型的虚粒子分布在边界的外部，与Libersky所使用的相似【23】。第二种类型的虚粒子按以下的方式构造，即给定一个实粒子i，则在边界外与实粒子对称处分布一个虚

27、粒子，这些虚粒子具有与相对应实粒子相同的压力和密度，但速度方向相反。第二种类型的虚粒子通常在边界条件不断变化的场合下使用【24】。3.6 时间积分可以使用一些标准的计算方法如跳蛙法(LF)、预估校正法、龙格一库塔法(RK)对普通微分方程进行数值积分来求解每一个粒子的物理变化。我们使用跳蛙法进行时间积分，跳蛙法的优点是计算时所需要的存储量低，而且在每一次计算中只需要进行一次优化估值。粒子的密度、速度、内能和位移按照下面的公式递推: 应该注意的是，跳蛙法是条件稳定的，稳定的条件为CFL(Courant一Friedrichs一Levy)条件，此条件一般会导致时间步长与光滑长度互成比例。在本文中，时间

28、步长按如下方式取值:式中，为Courant系数，一般取值在0.3左右【4】。【1】 宋康祖 紧支函数无网格方法研究【D】。清华大学博士学位论文。2000.【2】 Lucy L B.Numerical approach to testing the fission hypothesisj,Astronomical Journal.1977(82):1013-1024【3】 Gingold R A,Monaghan J J.Smoothed Particle Hydrodynamics:Theory and Application to Non-spherical starsC.Monthly N

29、otices of the Royal Astronomical Society.1977(181)375-389.【4】 Monaghan J J.An introduction to SPHj.Comput.Phys.Commun.1988,48:89-96.【5】 Monaghan J J.Why particle methods workj. SIAM j.Sci.Stat.Comput.1982,3(4);422-433【6】 Johnson G R,Beissel S R.Normalized smoothing functions for SPH impact computati

30、onsJ.Int J Num Meth Eng,1996,39:2725-2741.【7】 Swegle J W,Hicks D L, Attaways S W. Smoothed particle hydrodynamics stability analysisJ.J Comput Phys,1995,116:123-134.【8】 Dyka C T. Addressing tension instability in SPH methodsR.technical Report Nrl/MR/6384,NRL,1994.【9】 Johnson G R,stryk R A. Beissel S

31、 R.SPH for high velocity impact computationsJ. Comput.Methods Appl.Mech.Engrg.1996,139:347-373.【10】 Swegle J W, Attaway S W. On the feasibilility of using smoothed particle hydrodynamics for underwater explosion calculationsj.Comput Mtch,1995,17:151-168.【11】 Johnson G R,stryk R A,Beissel S R.SPH for

32、 high velocity impact computationsJ.Computer Methods in Applied Mechanics and Engineering,1996,139:347-373【12】 Libersky L D,petschek A G,et al. High strain lagrangian hydrodynamics:A three-dimensional SPH code for dynamic material responseJ.J Comput Phys,1993.109:67-75【13】 Johnson G R,Beissel S R.Normalized smoothing functions for SPH impact computationsJ.International Journal for Numerical Methods in Engineering,1996,39:2725-2741.【14】 张锁春.光滑质点流体动力学（SPH）方法（综述）J.计算