论文部分内容阅读
摘 要：光滑粒子动力学（SPH）方法是一种纯无网格粒子方法，目前已广泛应用于非牛顿流体流动问题。该文在探讨了光滑粒子动力学方法应用于非牛顿流体流动问题的研究现状的基础上，分析了该方法数值模拟中存在精度低和张力不稳定的问题及产生原因，并提出了相应的解决方案。
关键词：光滑粒子动力学 非牛顿 数值模拟 张力不稳定
中图分类号：O35 文献标识码：A 文章编号：1674-098X（2014）07（a）-0250-02
1 光滑粒子动力学方法应用简介
光滑粒子动力学（Smoothed Particle Hydrodynamics，SPH）方法是一种基于Lagrange描述的纯无网格粒子方法。它不需要使用任何网格，而是将连续体离散成有限个粒子，这些粒子携带所有物理量（密度、压力、速度、温度等）的信息。其次通过核函数将连续介质的控制方程转化为积分形式。这样整个流场变成了一系列粒子的表达，每个粒子的物理量都由周围粒子的物理量插值得到。
SPH方法最初由Lucy[1]、Gingold和Monaghan[2]提出并用于求解三维开放空间天体物理问题。1994年，Monaghan[3]首次将SPH方法用于自由表面流的数值模拟。随后，许多学者运用SPH方法对牛顿流体的自由表面流动问题[4-5]进行了模拟研究，大大推动了SPH方法在流体力学领域的发展。目前，SPH方法在聚合物自由面方面的研究也取得了一定进展。2005年，Ellero和Tanner[6]运用SPH方法模拟了低Reynolds数下粘弹性熔体剪切流动问题。2007年、2008年Rafiee[6]运用SPH方法模拟了非牛顿熔体的自由表面流问题。上述研究均是对简单流动进行研究。2006年，Fang等人[6]首次采用SPH方法研究了粘弹性液滴落到刚性固壁上的复杂流动问题，讨论了SPH方法中的人工应力参数的选取。2012年，本文作者等[6]对传统SPH方法进行改进的基础上研究了聚合物自由面问题。2010年，Fan等人[6]首次运用SPH方法研究了高粘性幂律流体的充模问题。
总之，随着计算机技术的快速发展，SPH方法在聚合物流动问题中的应用越来越广泛。
2 光滑粒子动力学方法中存在的缺点
与其他方法相比，SPH方法也是计算流体力学领域中一种较新的数值方法。随着应用范围越来越广，SPH方法的一些缺点也逐渐显现出来。SPH方法的缺点主要体现以下几方面[4，7-8]。
2.1 数值精度低、稳定性差
SPH方法在粒子分布均匀的区域内具有二阶精度，但在粒子分布不均匀区域或边界处至多具有一阶精度。经研究发现，SPH方法精度低的主要原因是：在粒子分布不均匀区域，某流体粒子影响域内的粒子数比较少；而对于边界附近，流体粒子的影响域被边界截断，边界外不存在粒子，从而出现比较严重的粒子缺失现象。而流动过程中必然出现粒子分布不均匀的情况。提高SPH方法的精度成为促进SPH方法进一步发展的重要因素。
2.2 压力振荡严重
SPH方法虽然在模拟自由表面流动问题中具有优势，尤其是对于大变形问题。然而，由SPH方法得到的压力却存在严重振荡。这是因为：除了上述SPH方法数值精度比较低的原因，另一个原因就是弱可压SPH方法中始终采用压力状态方程来求解压力。虽然有学者针对此问题提出了改进措施，如Colagrossi[6]采用密度初始化方法、耗散项方法等在一定程度上改善了压力振荡的情况，但是上述方法始终采用压力状态方程来求解压力，无法从根本上解决压力振荡问题。
2.3 存在张力不稳定
张力不稳定性或称“拉伸不稳定性”，该术语来源于传统SPH方法求解固体力学中材料拉伸问题时出现的一种非物理粒子簇集现象[7]。许多学者[6]在固体变形问题的数值模拟中对其出现的原因进行了研究和分析，但迄今为止关于张力不稳定性产生的本质原因没有一个统一的理论分析，也未能够提出普遍应用的克服方法。目前，应用比较广泛的克服张力不稳定性的方法是Monaghan和Gray等人[6]提出的人工应力法。
3 解决方案
上述SPH方法的缺点严重限制了SPH方法的推广与发展。因此，对SPH方法的改进有着重要的理论意义和应用价值。
针对传统SPH方法数值精度低的问题，众所周知，核函数在SPH方法中扮演着举足轻重的角色。因此，本文借助Taylor展开方法[4，8]，在核函数导数的左侧乘以一个修正矩阵，且此修正矩阵中不需要求解核函数导数，以此来提高SPH方法的数值精度，即
针对压力振荡严重的问题，本文考虑结合不可压缩SPH方法和弱可压SPH方法的优点，借助于有限差分方法的思想，对压力进行迭代求解，直至求得的压力计算得到的密度与初始密度的误差在允许范围之内时，迭代终止。
针对张力不稳定性问题，本文根据Swegle等人[7]提出的张力不稳定性结论，给出了一种容易施加的人工应力方法，其具体离散形式为
其中，为一个正常数（），为流体的应力张量分量，核函数的二阶导数分量为。上述人工应力方法离散的形式比已有的人工应力形式简单，且同时考虑了拉伸与压缩两种情况引起的不稳定。
4 结语
本文首先对SPH方法应用于聚合物流动问题的研究现状进行了阐述，然后根据相关文献及作者的数值模拟经验分析了运用SPH方法模拟聚合物流动问题时存在精度低和张力不稳定性问题及产生的原因，最后针对存在的问题，提出了相应的解决方案。
参考文献
[1] L.B.Lucy.A numerical approach to the testing of the ssion hypothesis[J].Astronomical Journal，1977（83）：1013-1024.
[2] R.A.Gingold，J.J.Monaghan. Smoothed particle hydrodynamics theory and application to non-spherical stars[J].Monthly Notices of the Royal Astronomical Society，1977（181）：375-389.
[3] J.J. Monaghan.Simulating free surface ows with SPH[J].Journal of Computational Physics，1994（110）：399-406.
[4] G.R.Liu，M.B.Liu.Smoothed Particle Hydrodynamics：A Mesh-free Particle Method[M].Singapore：World Scientific，2003.
[5] 刘谋斌，常建忠.光滑粒子动力学方法中粒子分布与数值稳定性分析[J].物理学报，2010（59）：3654-3662.
[6] 蒋涛.聚合物自由面问题的修正对称SPH方法模拟研究[D].博士研究生学位论文，西安：西北工业大学，2012.
[7] J.Swegle，W.D.L.Hicks， S.W.Attaway.Smooth particle hydrodynamics stability analysis[J].Journal of Computational Physics，1995， （116）：123-134.
[8] Jiang T，Ouyang J，Ren J，Yang B，Xu X.A mixed corrected symmetric SPH（MC-SSPH）method for computational dynamic problems[J].Computer Physics Communications，2012（183）：50-62.
关键词：光滑粒子动力学 非牛顿 数值模拟 张力不稳定
中图分类号：O35 文献标识码：A 文章编号：1674-098X（2014）07（a）-0250-02
1 光滑粒子动力学方法应用简介
光滑粒子动力学（Smoothed Particle Hydrodynamics，SPH）方法是一种基于Lagrange描述的纯无网格粒子方法。它不需要使用任何网格，而是将连续体离散成有限个粒子，这些粒子携带所有物理量（密度、压力、速度、温度等）的信息。其次通过核函数将连续介质的控制方程转化为积分形式。这样整个流场变成了一系列粒子的表达，每个粒子的物理量都由周围粒子的物理量插值得到。
SPH方法最初由Lucy[1]、Gingold和Monaghan[2]提出并用于求解三维开放空间天体物理问题。1994年，Monaghan[3]首次将SPH方法用于自由表面流的数值模拟。随后，许多学者运用SPH方法对牛顿流体的自由表面流动问题[4-5]进行了模拟研究，大大推动了SPH方法在流体力学领域的发展。目前，SPH方法在聚合物自由面方面的研究也取得了一定进展。2005年，Ellero和Tanner[6]运用SPH方法模拟了低Reynolds数下粘弹性熔体剪切流动问题。2007年、2008年Rafiee[6]运用SPH方法模拟了非牛顿熔体的自由表面流问题。上述研究均是对简单流动进行研究。2006年，Fang等人[6]首次采用SPH方法研究了粘弹性液滴落到刚性固壁上的复杂流动问题，讨论了SPH方法中的人工应力参数的选取。2012年，本文作者等[6]对传统SPH方法进行改进的基础上研究了聚合物自由面问题。2010年，Fan等人[6]首次运用SPH方法研究了高粘性幂律流体的充模问题。
总之，随着计算机技术的快速发展，SPH方法在聚合物流动问题中的应用越来越广泛。
2 光滑粒子动力学方法中存在的缺点
与其他方法相比，SPH方法也是计算流体力学领域中一种较新的数值方法。随着应用范围越来越广，SPH方法的一些缺点也逐渐显现出来。SPH方法的缺点主要体现以下几方面[4，7-8]。
2.1 数值精度低、稳定性差
SPH方法在粒子分布均匀的区域内具有二阶精度，但在粒子分布不均匀区域或边界处至多具有一阶精度。经研究发现，SPH方法精度低的主要原因是：在粒子分布不均匀区域，某流体粒子影响域内的粒子数比较少；而对于边界附近，流体粒子的影响域被边界截断，边界外不存在粒子，从而出现比较严重的粒子缺失现象。而流动过程中必然出现粒子分布不均匀的情况。提高SPH方法的精度成为促进SPH方法进一步发展的重要因素。
2.2 压力振荡严重
SPH方法虽然在模拟自由表面流动问题中具有优势，尤其是对于大变形问题。然而，由SPH方法得到的压力却存在严重振荡。这是因为：除了上述SPH方法数值精度比较低的原因，另一个原因就是弱可压SPH方法中始终采用压力状态方程来求解压力。虽然有学者针对此问题提出了改进措施，如Colagrossi[6]采用密度初始化方法、耗散项方法等在一定程度上改善了压力振荡的情况，但是上述方法始终采用压力状态方程来求解压力，无法从根本上解决压力振荡问题。
2.3 存在张力不稳定
张力不稳定性或称“拉伸不稳定性”，该术语来源于传统SPH方法求解固体力学中材料拉伸问题时出现的一种非物理粒子簇集现象[7]。许多学者[6]在固体变形问题的数值模拟中对其出现的原因进行了研究和分析，但迄今为止关于张力不稳定性产生的本质原因没有一个统一的理论分析，也未能够提出普遍应用的克服方法。目前，应用比较广泛的克服张力不稳定性的方法是Monaghan和Gray等人[6]提出的人工应力法。
3 解决方案
上述SPH方法的缺点严重限制了SPH方法的推广与发展。因此，对SPH方法的改进有着重要的理论意义和应用价值。
针对传统SPH方法数值精度低的问题，众所周知，核函数在SPH方法中扮演着举足轻重的角色。因此，本文借助Taylor展开方法[4，8]，在核函数导数的左侧乘以一个修正矩阵，且此修正矩阵中不需要求解核函数导数，以此来提高SPH方法的数值精度，即
针对压力振荡严重的问题，本文考虑结合不可压缩SPH方法和弱可压SPH方法的优点，借助于有限差分方法的思想，对压力进行迭代求解，直至求得的压力计算得到的密度与初始密度的误差在允许范围之内时，迭代终止。
针对张力不稳定性问题，本文根据Swegle等人[7]提出的张力不稳定性结论，给出了一种容易施加的人工应力方法，其具体离散形式为
其中，为一个正常数（），为流体的应力张量分量，核函数的二阶导数分量为。上述人工应力方法离散的形式比已有的人工应力形式简单，且同时考虑了拉伸与压缩两种情况引起的不稳定。
4 结语
本文首先对SPH方法应用于聚合物流动问题的研究现状进行了阐述，然后根据相关文献及作者的数值模拟经验分析了运用SPH方法模拟聚合物流动问题时存在精度低和张力不稳定性问题及产生的原因，最后针对存在的问题，提出了相应的解决方案。
参考文献
[1] L.B.Lucy.A numerical approach to the testing of the ssion hypothesis[J].Astronomical Journal，1977（83）：1013-1024.
[2] R.A.Gingold，J.J.Monaghan. Smoothed particle hydrodynamics theory and application to non-spherical stars[J].Monthly Notices of the Royal Astronomical Society，1977（181）：375-389.
[3] J.J. Monaghan.Simulating free surface ows with SPH[J].Journal of Computational Physics，1994（110）：399-406.
[4] G.R.Liu，M.B.Liu.Smoothed Particle Hydrodynamics：A Mesh-free Particle Method[M].Singapore：World Scientific，2003.
[5] 刘谋斌，常建忠.光滑粒子动力学方法中粒子分布与数值稳定性分析[J].物理学报，2010（59）：3654-3662.
[6] 蒋涛.聚合物自由面问题的修正对称SPH方法模拟研究[D].博士研究生学位论文，西安：西北工业大学，2012.
[7] J.Swegle，W.D.L.Hicks， S.W.Attaway.Smooth particle hydrodynamics stability analysis[J].Journal of Computational Physics，1995， （116）：123-134.
[8] Jiang T，Ouyang J，Ren J，Yang B，Xu X.A mixed corrected symmetric SPH（MC-SSPH）method for computational dynamic problems[J].Computer Physics Communications，2012（183）：50-62.