SPH法理论基础.docx

SPH法理论基础.docx

《SPH法理论基础.docx》由会员分享，可在线阅读，更多相关《SPH法理论基础.docx（14页珍藏版）》请在冰豆网上搜索。

SPH法理论基础

摘要

计算机数值仿真逐渐成为解决现代工程和科学问题的一条重要途径。

数值仿真能为理论提供测试和检验，有助于对复杂的物理问题加深认识，甚至还能帮助解释和发现新现象。

基于网格的数值方法虽然已经有广泛的应用，但是在很多方面仍存在缺乏之处，比方在计算流体动力学的大变形、运动物质交界面、自由外表等问题时，由于网格产生畸变导致计算误差过大或无法进展，从而使其在许多问题的应用上受到限制。

近年来，无网格法倍受关注，这种方法在许多应用中都优于传统的基于网格的有限元法、有限差分法以及有限体积法等数值方法。

本文主要研究新一代无网格方法---光滑粒子流体动力学方法〔SPH〕。

第一章绪论

1.1基于网格的方法

数值计算方法通常可以分为两种：

基于网格的方法和无网格方法。

通常对于物理控制方程的描述有两种根本方法：

欧拉描述法和拉格朗日描述法。

欧拉描述法是对空间的描述方法，其典型代表是有限差分法；拉格朗日描述法是对物质点的描述方法，其典型代表是有限元法。

欧拉描述和拉格朗日描述对应着两种不同的区域离散化网格：

欧拉网格和拉格朗日网格。

针对不同类型的问题，这两种网格在数值方法中都得到广泛应用。

1.1.1拉格朗日网格

拉格朗日网格，基于拉格朗日网格的数值方法在整个计算过程中网格是固定附着于物质上的，网格会随着物质的运动而运动，所以在物质点上的所有场变量的整个时间历程都可以很容易地追踪，常见的方法如有限元法等。

基于拉格朗日网格方法的优点是：

由于在相关的偏微分方程里不存在迁移项，所以程序在方案设计上会变得相对简单而且运行较快；由于只需在问题域布置网格，问题域外不需要布置，所以计算效率很高；不规那么或者复杂的几何形状可以用不规那么的网格来处理。

由于具有以上这些优点，拉格朗日方法得到广泛的应用，并且能成功地求解计算固体力学问题。

然后，基于拉格朗日网格的方法难以应用于具有极大网格变形的情况，因为其公式的形式是以网格为根底的，当网格变形太大的时候，公式的精度和求解都会受到很大的影响。

另外，由于时间步长是由最小单元尺寸所控制，假设网格太小就会影响计算的效率，甚至会导致计算失败。

1.1.2欧拉网格

欧拉网格，相对于拉格朗日网格，欧拉网格刚好相反，它是固定在模拟对象所处的空间上，模拟对象在固定网格单元上运动。

因此，在物质流过网格时，所有网格节点和网格单元依然固定在空间上而且不会随时间的改变而改变。

通过模拟质量、动量和能量经过网格单元边界的通量，可以计算质量、速度和能量等等的分布。

在整个计算过程中网格单元的形状和体积都保持不变。

由于欧拉网格在时间和空间上都是固定的，物体的大变形不会引起网格本身的任何变化，所以再以物质流动为主体的计算流体力学领域里，较为广泛的应用欧拉法。

但是，欧拉法依然有很多缺点。

由于欧拉法不能用固定网格来追踪物质的运动，所以很难分析物质上固定点的场变量的变化情况而只能得到空间上固定的欧拉网格的场变量的变化情况；由于欧拉法追踪的是流过网格单元边界的质量、动量和能量的通量，所以自由外表、变形边界和运动物质交界面的位置就很难准确确定；由于欧拉法需要在整个计算区域上都覆盖网格，所以有时为了提高计算效率而使用较粗糙的网格，这样就会降低离散化区域的求解精度。

1.1.3基于网格的数值方法的局限性

传统的基于网格的数值方法如有限差分法和有限元法已经广泛地应用在计算流体力学和计算固体力学的各个领域中，是现在进展区域离散化和数值离散化模拟的主要方法。

但是基于网格的数值方法在很多方面仍存在缺乏之处，比方在计算流体动力学中的大变形、运动物质交界面、自由外表等问题时，由于网格产生畸变导致计算误差过大或无法进展，使其在许多问题的应用上受到限制。

在基于网格的数值方法中，数值模拟的先决条件就是在问题域生成网格。

对于基于欧拉网格的方法，在固定的欧拉网格上要准确确定自由外表、变形边界、运动交界面和不均匀物质之间的位置是一项非常困难的工作。

并且欧拉方法也不适用于研究如粒子流动这类问题。

即需要在固定的物质体监控材料特性的问题。

对于拉格朗日网格法如有限元法，进展模拟前必须要在研究对象上建立网格，这项操作常常占用很大的计算工作量。

当所研究对象是一系列离散的物质点时，同样不适合使用基于网格的方法，用连续的基于网格的方法对这些离散系统进展模拟通常不是好的选择。

在计算大变形问题时，比方高速碰撞、金属加工成型、动态裂纹扩展、流固耦合和应变局部化等，基于网格的方法遇到的困难更大。

原因是发生大变形时，网格畸变过大，使得计算中断，有限元方法通常采用单元“销蚀〞法或重分网格技术来克制这种困难。

但是单元“销蚀〞法本身缺乏物理依据，纯粹是为了使计算进展下去的一种数值手段。

而网格重分技术在节点重新分配物理量时，很难保证系统动量、能量守恒，因而导致计算的精度下降。

此外，网格重分技术不是很容易实现，为了更好地解决大变形问题，必须对网格有新的处理方法或去除网格，所以各种无网格方法相继被提出来。

1.2无网格方法

无网格方法产生于三十多年以前，当时开展很慢，直到近几年，由于无网格方法的近似函数不依赖于网格，在分析涉及大变形的问题中被认为优于传统的基于网格的有限差分法和有限元法，才受到了众多研究者的高度重视，成为国际计算仿真界的研究热点之一，并得到了迅速的开展【1】。

无网格方法的主要思想是：

通过使用一系列任意分布的节点〔或粒子〕来求解具有各种边界条件的积分方程或偏微分方程组，从而得到准确稳定的数值解，这些节点或粒子之间不需要网格进展连接。

由于无网格法采用基于点的近似，可以彻底或局部地消除网格，不需要网格的初始划分和重构，因此不仅可以保证计算的精度，而且可以减小计算的难度。

经过几十年的研究，目前已开展了十余种无网格方法。

最早提出的无网格方法是光滑粒子动力学方法〔SmoothedParticleHydrodynamics，简称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-13】。

近几年来，我国的学者也开场关注SPH计算方法，中科院的锁春对SPH方法进展了综述【14】，国防科学技术大学的邓方刚研究了SPH在柱坐标体系下的应用【15】。

国防科大的贝新源、岳宗五等将SPH方法用于高速碰撞问题研究【16】。

1.3光滑粒子动力学方法

光滑粒子流体动力学〔SPH-Smoothedparticlehydrodynamics〕方法是近30年开展起来的一种新的纯拉格朗日无网格粒子法，它最初是用于解决三维开放空间的天体物理学问题【2,3】，现已被广泛地研究和扩展，并被应用于具有材料强度的动态响应问题和具有大变形的流体动力学问题。

SPH方法和它的派生方法是现在粒子法的主要类型。

SPH方法通过对邻近粒子进展加权平均而得到稳定的光滑近似性质，该方法应用在流体动力学问题的围。

由于SPH方法中的自适应性、粒子性质和拉格朗日性质的和谐结合，使其在工程和科学领域都得到了实际应用。

国外对SPH方法研究的比拟早，已经被应用到很多领域中。

SPH方法最初在天体研究中应用较多，如双子星和行星碰撞、超新星、星系的形成和瓦解、黑洞和中子星合并、白矮星的单个或多重爆炸甚至是宇宙的进化等问题的研究。

第二章SPH方法在流体动力学问题中的应用

2.1光滑粒子动力学原理

流体动力学问题的求解主要是解基于密度、速度、能量等变量场的偏微分方程组。

除了一些简单问题可以得到偏微分方程组的解析解以外，大局部的流体动力学问题只能寻求其数值解。

光滑粒子动力学的原理如下：

1）使用SPH方法求解流体动力学问题时，问题域用一系列任意分布的粒子来表示，粒子之间没有任何连接，因此SPH方法是无网格性质的，通常为了求解方便，初始粒子匀排列；

2）场函数用积分表示法来近似，在SPH方法中称为核近似法；

3）然后应用支持域的相邻粒子对应的值叠加求和取代场函数的积分表达式来对场函数进展粒子近似，由于在每一个时间步都要进展粒子近似，支持域的有效粒子为当前时刻支持域的粒子，因此SPH方法具有自适应性；

4）将粒子近似法应用于所有偏微分方程组的场函数相关项中，将偏微分方程组进展离散，可知SPH是一种纯拉格朗日方法；

5）粒子被附上质量后，那么意味着这些粒子是真实的具有材料特性的粒子；最后应用显式积分法得到所有粒子的场变量随时间的变化值。

从以上分析可以看到，SPH方法是一种纯拉格朗日的具有无网格、自适应属性的流体动力学求解方法。

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

使用SPH方法求解流体动力学问题一般分两步进展。

第一步是积分表示，即使用积分表达式对函数进展近似，函数的积分表达式可以通过对核函数的影响区域积分进展近似。

第二步是粒子表示，即通过对最近相邻粒子的值进展累加求和来近似离散点的函数值。

2.2.1函数的积分表达

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

其中，为函数f（x）的近似值，x为位置矢量，Ω为包含x的几分体，W（x-x’,h）为光滑函数，其依赖于两点之间的距离|x-x’︳和光滑长度h。

光滑函数在SPH近似法中其着重要的作用，它决定了函数表达式的精度和计算效率。

光滑函数W常常选用偶函数，其应满足3个条件。

第一个条件为正那么化条件，如下所示

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

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

这里，

第三个条件是紧支性条件

式中，K是与点x处光滑函数相关的常数，并确定光滑函数的有效围〔非零〕。

此有效围称作点x处光滑函数的支持域的积分。

因此，一般来说积分域Ω就是支持域。

根据分布积分和散度定理，将第一个式子经过变换，对f〔x〕的导数的估计为

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

2.2.2函数的粒子表达

SPH方法最终要将函数积分近似表达式转化为支持域所有粒子叠加求和的离散化形式，流体的问题域被离散化为有限数量的粒子，其中每个粒子具有独立的质量、密度及其它物理属性。

经过离散化，函数的积分表达式可以写成如下的粒子近似式

其中，j和i表示粒子编号，mj和ρj分别j粒子的质量和密度，N是粒子的总数，Wij=W（xi-xj，h）=W（|xi-xj|，h）。

上式说明，在SPH方法中，对任意一个函数f，它在某一位置处的函数值可以通过应用光滑函数W对其光滑长度为h的紧支域所有粒子插值求和的形式来表达【17】。

光滑长度h在SPH方法中非常重要，它在计算过程中决定了计算的实用性、有效性和可靠的适应性。

光滑长度h直接影响到计算结果的效率和结果的精度。

假设h太小，那么在以kh为半径的支持域没有足够的粒子对所求粒子施加影响，这样会导致计算结果精度较低。

如果h太大，太多的粒子对所求粒子产生影响，那么可能将粒子所有细节信息和局部特性忽略，导致粒子趋于雷同，同样会影响到结果的精度。

具体光滑长度应该选择多大与所研究的问题相关，没有一个广义的光滑长度。

2.2.3光滑函数

SPH通过光滑函数来进展积分表示。

光滑函数非常重要，因为它不仅决定函数表达式的形式，定义粒子支持域的大小，而且还决定核近似和粒子近似的一致性和计算效率。

常用的光