光滑粒子流体动力学方法SPH

SPH光滑粒子流体动力学方法光滑粒子流体动力学方法(smoothed particle hydrodynamics)无网格方法的主要思想无网格方法的主要思想：通过使用一系列任意分布的节点（或：通过使用一系列任意分布的节点（或粒子）来求解各种各样边界条件的积分方程或偏微分方程组粒子）来求解各种各样边界条件的积分方程或偏微分方程组(PDEs)从而得到精确稳定的数值解，这些节点或粒子之间不需要网格进行从而得到精确稳定的数值解，这些节点或粒子之间不需要网格进行连接。连接。 Lucy, Gingold(1977)分别提出了分别提出了SPH方法，最早用于天体物理方法，最早用于天体物理现象的模拟，随后别广泛地应用于连续固体力学和流体力学中。现象的模拟，随后别广泛地应用于连续固体力学和流体力学中。 11kg弹丸弹丸1418m/s撞靶速度下穿靶过程的数值模拟撞靶速度下穿靶过程的数值模拟1500m/s速度下弹体侵彻混凝土靶变形过程的数值模拟速度下弹体侵彻混凝土靶变形过程的数值模拟侵彻过程弹体温度分布云图侵彻过程弹体温度分布云图碎浪与弹性挡墙之间的相互作用碎浪与弹性挡墙之间的相互作用近似函数近似函数构造方法构造方法偏微分方程偏微分方程的离散形式的离散形式核估计方法核估计方法(Kernel Approximation, KA)移动最小二乘法移动最小二乘法(Moving Least Square, MLS)再生核估计方法再生核估计方法(Repuducing Kernel Method, RKM)径向基函数方法径向基函数方法(Radial Basic Function, RBF)单位分解方法单位分解方法(Patition of Unity, PU)强形式强形式以各种全局或局部加权余量法为统一框架的弱形式以各种全局或局部加权余量法为统一框架的弱形式两两条条主主线线无网格法无网格法(Meshfree Methods)强形式强形式是直接从微分方程及其定解条件出发，将近似函数及其导是直接从微分方程及其定解条件出发，将近似函数及其导数的估计形式带入基本方程、本构方程和初边值条件中去，联立数的估计形式带入基本方程、本构方程和初边值条件中去，联立方程进行求解。该方法思路简单，便于程序编制，应用范围广泛，方程进行求解。该方法思路简单，便于程序编制，应用范围广泛，在流体和固体的计算中都有所发展，适用于计算激波、高速冲击、在流体和固体的计算中都有所发展，适用于计算激波、高速冲击、爆轰、穿甲等冲击动力学问题。但此类算法的精度较低，稳定性爆轰、穿甲等冲击动力学问题。但此类算法的精度较低，稳定性较差，且边界条件的引入比较困难。较差，且边界条件的引入比较困难。弱形式弱形式就是从加权余量法或变分原理出发，把微分方程及其定解就是从加权余量法或变分原理出发，把微分方程及其定解条件转换成弱形式条件转换成弱形式(Weak Form)或或Galerkin形式，即用测试函数形式，即用测试函数(Test function)与控制方程相乘后在全局或部分区域内积分，并利与控制方程相乘后在全局或部分区域内积分，并利用高斯散度定理得到不同形式的弱形式，然后进行离散化求解。用高斯散度定理得到不同形式的弱形式，然后进行离散化求解。通过引进新的无网格近似函数构造方法或采用新的偏微分方程的通过引进新的无网格近似函数构造方法或采用新的偏微分方程的离散形式，就可以期待开发出更加高效和精确的无网格方法。离散形式，就可以期待开发出更加高效和精确的无网格方法。 SPH方程的构造常按两个关键步进行。方程的构造常按两个关键步进行。第一步为积分表示法，又称场函数近似法；第一步为积分表示法，又称场函数近似法；第二步为粒子近似法。第二步为粒子近似法。光滑粒子流体动力学光滑粒子流体动力学一种无网格粒子法，一种无网格粒子法，湖南大学出版社，湖南大学出版社，G. R. Liu，M. B. Liu著著场函数核近似法（积分表示法）场函数核近似法（积分表示法）函数核近似的标准表达式：函数核近似的标准表达式： xxxxxdhWff,h是定义光滑函数是定义光滑函数W的影响区域的光滑长度。的影响区域的光滑长度。W被称为光滑核函数被称为光滑核函数(smoothing kernel function)或光滑或光滑函数函数(smoothing function)，简称为核，简称为核(kernel)函数。函数。粒子近似法粒子近似法与与SPH核近似法相关的连续积分表示式，可转化为支核近似法相关的连续积分表示式，可转化为支持域内所有粒子叠加求和的离散化形式。持域内所有粒子叠加求和的离散化形式。 NjjjjjjjNjjjjjjNjjjjNjjjhWfmmhWfVhWfVhWfdhWff1111,1,1,xxxxxxxxxxxxxxxxx粒子近似法粒子近似法在粒子在粒子i处的函数的粒子近似式最终可写为：处的函数的粒子近似式最终可写为：ijjNjjjiWfmfxx1hWWjiij,xx 上式说明了粒子上式说明了粒子i处的任一函数值可通过应用光滑函数对处的任一函数值可通过应用光滑函数对其紧支域内所有粒子相对应的函数值进行加权平均近似。其紧支域内所有粒子相对应的函数值进行加权平均近似。ijijNjjjiWfmfxx1SPH计算公式计算公式光滑函数光滑函数最近相邻粒子搜索法最近相邻粒子搜索法(NNPS)人工粘度人工粘度边界处理边界处理交界面处理交界面处理光滑长度的更新光滑长度的更新SPH方程的求解方程的求解激波管问题激波管问题SPH程序结构程序结构目目 录录1、密度的粒子近似法、密度的粒子近似法 由于粒子的分配与光滑长度的变化主要依赖于密度，故在由于粒子的分配与光滑长度的变化主要依赖于密度，故在SPH法法中密度近似法非常重要。中密度近似法非常重要。 在在SPH法中有两种方法对密度进行展开，第一种方法是对密度直法中有两种方法对密度进行展开，第一种方法是对密度直接用接用SPH近似法，称为近似法，称为密度求和法密度求和法。第二种方法是。第二种方法是连续性密度法连续性密度法，通，通过应用过应用SPH近似法的概念对连续性方程进行转换而得到。近似法的概念对连续性方程进行转换而得到。NjijjjNjijjiWmWm11SPH计算公式计算公式ijjNjjjiWfmfxx1NjijjiWm1密度求和法密度求和法：改进方案（正则化）改进方案（正则化）此方法可提高自由边界处和相同材此方法可提高自由边界处和相同材料粒子密度不连续交界面处的精度料粒子密度不连续交界面处的精度iijNjjjjiixWvmdtd1011iijNjjjxWmNjiijjjiiiijNjijjixWmvxWvm11iijNjijjjiixWvmdtd1jiijvvv连续性密度法连续性密度法：xvdtdijijNjjjiWfmfxx1 对于广义流体问题的模拟，应用修正的密度求和法可对于广义流体问题的模拟，应用修正的密度求和法可得到较好的结果，对于具有强间断问题的模拟（如爆炸、得到较好的结果，对于具有强间断问题的模拟（如爆炸、高速冲击等），应优先选取连续性密度法。高速冲击等），应优先选取连续性密度法。2、动量方程的粒子近似法、动量方程的粒子近似法xdtdv1将以上两式相加可得：将以上两式相加可得：iijNjjjjiixWmdtdv11011iijNjjjiiiijNjjiijxWmxWmiijNjjijijixWmdtdv1将动量方程等号右端的梯度项直接应用将动量方程等号右端的梯度项直接应用SPH粒子近似法进行变换得：粒子近似法进行变换得：此外，有：此外，有：动量守恒方程：动量守恒方程：ijijNjjjiWfmfxx1此对称方程的优点为：可降低此对称方程的优点为：可降低粒子不一致问题产生的误差。粒子不一致问题产生的误差。3、能量方程的粒子近似法、能量方程的粒子近似法xvpdtdeiijNjijjjiiiixWvmpxvp1xpvpvxxvp1iijNjijjjjiiixWvpmxvp11iijNjijjijijiixWvppmxvp121iijNjijjijijixWvppmdtde121能量守恒方程：能量守恒方程：ijijNjjjiWfmfxx1iijNjijjjiixWvmdtd1dtdpxvpxvp22光滑函数光滑函数光滑函数的性质：光滑函数的性质：1,xxxdhW一、正则化条件一、正则化条件由于光滑函数的积分值等于由于光滑函数的积分值等于1，故此条件也称为归一化条件。，故此条件也称为归一化条件。二、当光滑长度趋向于零时具有狄拉克函数性质二、当光滑长度趋向于零时具有狄拉克函数性质xxxxhWh,lim0三、紧支性条件三、紧支性条件0,hWxxhxxMonaghan和和Lattanzio在三次样条函数的基础上提出了在三次样条函数的基础上提出了称为称为B样条函数的光滑函数样条函数的光滑函数 ：光滑函数光滑函数2, 021,26110,2132,332RRRRRRahRWd现有现有SPH文献中最为广泛应用的光滑函数文献中最为广泛应用的光滑函数在一维、二维和三维空间中分别有：在一维、二维和三维空间中分别有： ， 和和 。had12715h323hhRxx三次样条函数及其一阶导数三次样条函数及其一阶导数四次样条函数及其一阶导数四次样条函数及其一阶导数五次样条函数及其一阶导数五次样条函数及其一阶导数光滑函数一览表光滑函数一览表最近相邻粒子搜索法最近相邻粒子搜索法(NNPS)一般将包含在支持域中的粒子称为相关粒子的一般将包含在支持域中的粒子称为相关粒子的最近相邻粒子最近相邻粒子(NNP)。寻找最近相邻粒子的过程通常称为寻找最近相邻粒子的过程通常称为最近相邻粒子搜索最近相邻粒子搜索(NNPS)。在在SPH方法中常用的三种方法中常用的三种NNPS方法为：方法为：全配对搜索法全配对搜索法(all-pair search)链表搜索法链表搜索法(linked-list search algorithm)树形搜索法树形搜索法(tree search algorithm)全配对粒子搜索法全配对粒子搜索法对于给定的粒子对于给定的粒子i，应用全配对搜索法即是计算粒子，应用全配对搜索法即是计算粒子i到粒子到粒子j的距离。的距离。若该距离小于粒子若该距离小于粒子i的支持域的半径的支持域的半径 ，则粒子，则粒子j为粒子为粒子i的支持域内的支持域内的粒子。的粒子。h 在二维空间里应用链表在二维空间里应用链表搜索法搜索最近相邻粒子搜索法搜索最近相邻粒子在光滑长度为空间常量的情况下，应用链在光滑长度为空间常量的情况下，应用链表搜索法非常有效。表搜索法非常有效。在实现链表算法时，要在问题域上铺设一在实现链表算法时，要在问题域上铺设一临时网格，网格单元的空间大小应选取与临时网格，网格单元的空间大小应选取与支持域的空间大小一致。那么，对于给定支持域的空间大小一致。那么，对于给定的粒子的粒子i，其相邻粒子只能在同一网格单元，其相邻粒子只能在同一网格单元内，或者在紧密相邻的单元内。所以，当内，或者在紧密相邻的单元内。所以，当时，一维、二维和三维空间里的搜索范围时，一维、二维和三维空间里的搜索范围分别为分别为3，9，27个单元内。个单元内。链表搜索法存在的问题是，当光滑长度可链表搜索法存在的问题是，当光滑长度可变时，网格空间就不一定适应每一个粒子，变时，网格空间就不一定适应每一个粒子，此时若再应用链表搜索法，则搜索效率就此时若再应用链表搜索法，则搜索效率就会很低。会很低。2h 树形搜索法非常适宜求解具有可变光滑长度的问题。这种算法是通过粒子树形搜索法非常适宜求解具有可变光滑长度的问题。这种算法是通过粒子的位置来构造有序树。一旦树形结构构造起来后，便能高效地搜索最近相邻粒的位置来构造有序树。一旦树形结构构造起来后，便能高效地搜索最近相邻粒子。子。 树形搜索法将最大问题域递归分割成一个个卦限，直到每一个卦限内只包树形搜索法将最大问题域递归分割成一个个卦限，直到每一个卦限内只包含一个粒子为止。树形结构构造完成后，即可以开始进行最近相邻粒子搜索。含一个粒子为止。树形结构构造完成后，即可以开始进行最近相邻粒子搜索。 给定任一粒子给定任一粒子i, 并以并以i为中心，用边长为为中心，用边长为 的立方体将粒子包围起来，然的立方体将粒子包围起来，然后再检测粒子后再检测粒子i的搜索立方体空间是否与并列层次内的其他节点所占的空间有重的搜索立方体空间是否与并列层次内的其他节点所占的空间有重合的地方。若没有，则中止往下搜索；若有，则继续往下一层搜索。直到所搜合的地方。若没有，则中止往下搜索；若有，则继续往下一层搜索。直到所搜索到的当前节点处只有一个粒子为止。接着，检查此粒子是否在给定粒子的支索到的当前节点处只有一个粒子为止。接着，检查此粒子是否在给定粒子的支持域内，若是，则将其记为粒子的相邻粒子。持域内，若是，则将其记为粒子的相邻粒子。ih2 迄今为止，与迄今为止，与SPH相关的论著中，相关的论著中，Monaghan型的人工粘度是最型的人工粘度是最为广泛使用的人工粘度。它不仅将动能转换为热能，提供了冲击波阵为广泛使用的人工粘度。它不仅将动能转换为热能，提供了冲击波阵面必不可少的耗散，而且防止了粒子相互接近时的非物理穿透。具体面必不可少的耗散，而且防止了粒子相互接近时的非物理穿透。具体表达式如下：表达式如下： 0, 00,2ijijijijijijijijijcxvxv22ijijijijijhxxvjiijccc21jiij21jiijhhh21jiijvvvjiijxxx人工粘度人工粘度 在以上所有方程中在以上所有方程中 和和 为标准常数，一般取值在为标准常数，一般取值在1.0左右。因子左右。因子 用于防止粒子相互靠近时产生的数值发散。式中与用于防止粒子相互靠近时产生的数值发散。式中与 相关的项得到的是相关的项得到的是体积粘度，而与体积粘度，而与 相关的项是用于防止在高马赫数时粒子的相互穿透。相关的项是用于防止在高马赫数时粒子的相互穿透。式中所给出的人工粘度被引入到式中所给出的人工粘度被引入到SPH方程的压力项中。方程的压力项中。ijh1 . 0边界处理边界处理 由于在边界上或邻近边界处的粒子存在缺陷，即在积分的时候会由于在边界上或邻近边界处的粒子存在缺陷，即在积分的时候会被边界截断，故而被边界截断，故而SPH方法不能完全适用于整个区域。在边界上或邻方法不能完全适用于整个区域。在边界上或邻近边界处的粒子只受到边界内的粒子的影响作用，而边界外由于没有近边界处的粒子只受到边界内的粒子的影响作用，而边界外由于没有粒子，故而边界外不对粒子产生影响。这种单边影响作用会导致求解粒子，故而边界外不对粒子产生影响。这种单边影响作用会导致求解结果错误，因为在固定边界表面，虽然粒子速度为零，但是其他变量，结果错误，因为在固定边界表面，虽然粒子速度为零，但是其他变量，如密度，则不一定为零。如密度，则不一定为零。边界处理边界处理 在在Liu等的研究中，使用了虚粒子来处理固定边界条件，提出了等的研究中，使用了虚粒子来处理固定边界条件，提出了两种类型的虚粒子。第一种类型的粒子（型号两种类型的虚粒子。第一种类型的粒子（型号I）设置在固定边界上；）设置在固定边界上；第二种类型的粒子（型号第二种类型的粒子（型号II ）分布在边界的邻域内。）分布在边界的邻域内。型号型号II 的虚粒子具有与相应实粒子相同的密度和压力，但速度相反。的虚粒子具有与相应实粒子相同的密度和压力，但速度相反。型号型号I的虚粒子被引入到实粒子的核函数核粒子近似法中。的虚粒子被引入到实粒子的核函数核粒子近似法中。边界处理边界处理 当类型当类型I的虚粒子成为邻近边界处的实粒子的相邻粒子时，则会的虚粒子成为邻近边界处的实粒子的相邻粒子时，则会在沿着两粒子的中心线处对实粒子产生一个作用力。在沿着两粒子的中心线处对实粒子产生一个作用力。1, 01,0020021ijijijijnijnijijrrrrrxrrrrDPB式中：参数式中：参数n1和和n2一般取值分别为一般取值分别为12和和4。D是由具体问题而定的参数，是由具体问题而定的参数，一般取与速度最大值的平方相等的量级。截止半径一般取与速度最大值的平方相等的量级。截止半径r0在此问题的模拟在此问题的模拟分析中非常重要，在一般情况下，分析中非常重要，在一般情况下， r0的取值与粒子的初始间距的大小的取值与粒子的初始间距的大小相近。相近。 型号型号II的虚粒子没有固定的参数。它们是在每个计算步中由对应的的虚粒子没有固定的参数。它们是在每个计算步中由对应的实粒子对称产生的。可以应用型号实粒子对称产生的。可以应用型号II的虚粒子来处理固定边界和自由表的虚粒子来处理固定边界和自由表面。面。 通过数值算例的测试验证了应用虚粒子来处理边界的可靠性和有效通过数值算例的测试验证了应用虚粒子来处理边界的可靠性和有效性。其不仅仅提高了性。其不仅仅提高了SPH近似法在边界区域处的精度，而且防止了粒子近似法在边界区域处的精度，而且防止了粒子非物理穿透边界。非物理穿透边界。交界面处理交界面处理 由于由于SPH方法具有拉格朗日性质和粒子性质，在整个演变过程中，来自不方法具有拉格朗日性质和粒子性质，在整个演变过程中，来自不同介质的相互接触的粒子可能会随着运动而分离，甚至有可能不再相邻。同介质的相互接触的粒子可能会随着运动而分离，甚至有可能不再相邻。交界面处理交界面处理 如果求和仅仅是使用相同材料的粒子，则离得很近但是材料不同如果求和仅仅是使用相同材料的粒子，则离得很近但是材料不同的两粒子就不再是相邻粒子，因此，在交界面附近的区域，模拟就会的两粒子就不再是相邻粒子，因此，在交界面附近的区域，模拟就会遇到如一般边界的粒子缺陷问题。相反，如果求和使用不同材料的粒遇到如一般边界的粒子缺陷问题。相反，如果求和使用不同材料的粒子，那么离得很近但是材料不同的两粒子就可以认为是相邻粒子，因子，那么离得很近但是材料不同的两粒子就可以认为是相邻粒子，因此能够减少在此能够减少在SPH核近似和粒子近似中的边界缺陷问题。然而，这种核近似和粒子近似中的边界缺陷问题。然而，这种处理将会导致一些非物理穿透或掺杂的问题。在大部分的情况下非物处理将会导致一些非物理穿透或掺杂的问题。在大部分的情况下非物理穿透都不是很严重，仅仅有一或两层不同材料的粒子发生相互穿透。理穿透都不是很严重，仅仅有一或两层不同材料的粒子发生相互穿透。对于具有高强度载荷作用的问题，穿透对精度来说可能是致命的，有对于具有高强度载荷作用的问题，穿透对精度来说可能是致命的，有时甚至会造成时甚至会造成SPH程序运行的崩溃。程序运行的崩溃。 Liu等人提出了一种粒子对粒子的交界面算法，无论是对于高速等人提出了一种粒子对粒子的交界面算法，无论是对于高速冲击问题还是水下爆炸问题这种算法都是非常有效的。在此算法中，冲击问题还是水下爆炸问题这种算法都是非常有效的。在此算法中，物质交界面通过核近似求和来处理，即在求解守恒方程时考虑了不同物质交界面通过核近似求和来处理，即在求解守恒方程时考虑了不同材料粒子之间的相互作用。然而，这样处理还不是很完善，因为在水材料粒子之间的相互作用。然而，这样处理还不是很完善，因为在水下爆炸过程中，高压爆炸气体和周围的水介质之间的激烈相互作用通下爆炸过程中，高压爆炸气体和周围的水介质之间的激烈相互作用通常会导致交接面附近的非物理穿透和掺杂。因此在交界面附近的不同常会导致交接面附近的非物理穿透和掺杂。因此在交界面附近的不同材料的粒子，当它们彼此之间趋于穿透时就对其应用一个力的特殊处材料的粒子，当它们彼此之间趋于穿透时就对其应用一个力的特殊处理方法。如果以下条件满足就认为穿透发生：理方法。如果以下条件满足就认为穿透发生：12ijjirhhpe 所施加的惩罚力是成对地作用在两个相互接近的粒子上，力的所施加的惩罚力是成对地作用在两个相互接近的粒子上，力的方向是沿着粒子的中心线方向。方向是沿着粒子的中心线方向。式中的参数式中的参数 和和n1、n2 分别取分别取105，6和和4。 惩罚力与接触粒子求和的结合应用，虽然在交接面附近仍然存在惩罚力与接触粒子求和的结合应用，虽然在交接面附近仍然存在着数值振荡，但是却能够很好地预防在水下爆炸模拟中出现的非物理着数值振荡，但是却能够很好地预防在水下爆炸模拟中出现的非物理穿透问题。穿透问题。101221peperxpepepPBijijnnijp光滑长度的更新光滑长度的更新 在在SPH方法中，光滑长度非常重要，直接影响到求解的计算效率方法中，光滑长度非常重要，直接影响到求解的计算效率和精度。和精度。NjijijijiWvvmdtd1ijniNjnjnijnininiWvvmdhdtdh1dtddhdtdhnininini在下一个时间步，光滑长度变为：在下一个时间步，光滑长度变为：tdtdhhhninini1 Benz(1989)提出了一种对光滑长度提出了一种对光滑长度h进行动态变换的方法，即在进行动态变换的方法，即在连续性方程中将光滑函数对时间求导：连续性方程中将光滑函数对时间求导：tdtdhhhninini1 为了模拟具有极大密度不均匀的问题，为了模拟具有极大密度不均匀的问题，Liu等推导出了以下优化和等推导出了以下优化和松弛过程，此过程的目的是使每个粒子与数目近似为常量的最近相邻粒松弛过程，此过程的目的是使每个粒子与数目近似为常量的最近相邻粒子相互作用，由两步组成：预测步和修正步。子相互作用，由两步组成：预测步和修正步。预测步：预测步：修正步：修正步：一旦得到一旦得到 ，则可以确定当前的相邻粒子数，则可以确定当前的相邻粒子数 。若。若 与与 大大致相近，则其是可取的。若误差为百分之几，则要在致相近，则其是可取的。若误差为百分之几，则要在1.0附近对松弛因附近对松弛因子子 进行调整，以得到新的进行调整，以得到新的 和和 。重复此过程，直到两者大致相。重复此过程，直到两者大致相近为止。近为止。1nih1niN1niN0iN1nih1niN其中其中 为松弛因子，初始取值为为松弛因子，初始取值为1.0，然后在随后的优化和松弛过程步，然后在随后的优化和松弛过程步骤中在骤中在1.0附近进行调节。附近进行调节。优化和松弛过程优化和松弛过程 SPH方程的求解方程的求解 iiNjijijiijjjiijiNjijiijjjiijiNjijijijidtdWppmdtdeWppmdtdWmdtdvxvvvvv122122121 上述的上述的SPH方程是一系列与时间相关的联立的常微分方程，可以方程是一系列与时间相关的联立的常微分方程，可以使用标准的数值分析算法进行求解，如：使用标准的数值分析算法进行求解，如：leapfrog(LF)、Predictor-corrector以及以及Runge-Kutta(RK)等。在实际应用中，由于等。在实际应用中，由于leapfrog法对法对存储的需求量低，而且计算效率高，故很受欢迎。存储的需求量低，而且计算效率高，故很受欢迎。SPH程序结构程序结构 0 x10v5 . 2e1p001875. 0 x0 x25. 0795. 1e1795. 0p0075. 0 x模拟中共使用了模拟中共使用了400个粒子，所有粒子具有相同的质量。个粒子，所有粒子具有相同的质量。320个粒子个粒子均匀地分布在高密度区域，剩下的均匀地分布在高密度区域，剩下的80个粒子均匀地分布在低密度区个粒子均匀地分布在低密度区域。这样初始密度分布的目的是为了得到所要求的沿着管方向密度域。这样初始密度分布的目的是为了得到所要求的沿着管方向密度不连续分布形状。不连续分布形状。0v激波管问题激波管问题 t=0.20s时的密度分布图时的密度分布图t=0.20s时的速度分布图时的速度分布图t=0.20s时的压力分布图时的压力分布图