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