非凸区域上SPH计算结果后处理方法研究

郑 俊，于开平，张嘉钟，魏英杰

(哈尔滨工业大学航天学院，150001哈尔滨，zhengjun324203@126.com)

SPH(smoothed particle hydrodynamics)是一种基于Lagrangian描述的无网格方法［1］，其计算结果的后处理一般是直接将存储于粒子上的数据与颜色矩阵关联来获得物质域上物理量分布［2］，无法得到基于网格的算法后处理的连续的等值线(面)图、云图、流线图，且不易进行微积分运算以及切片上的数据可视化等.Massidda［3］指出无网格计算可视化是一项新的技术.Daniel［4］开发的面向SPH后处理的开源软件Splash，需在Linux下编译，主要面向天文领域，且该软件包对于水动力学等其它领域的适用性还有待拓展.Biddiscombe等［5］基于Paraview开发了可处理SPH数据的程序，但其处理时还是直接在粒子上，并不方便，功能也不全面.无网格法计算可视化的理论［6-9］，主要是采用Delaunay三角化方法，将无网格点集进行三角化，获得非结构化网格，然后基于非结构化网格进行后处理，可得到传统计算可视化所得的各种图形.但这些研究只针对于凸区域上的无网格点集，而非凸区域上的无网格计算可视化，鲜有报导.但在SPH所计算的液体飞溅、高速撞击中的材料碎片等问题中，计算域都是非凸的，对非凸区域上的无网格点集进行Delaunay三角化将获得一些并不属于物质域的单元，称这些单元为空白单元;如何方便的过滤掉这些单元，国内还没有较系统的工作.文献［6］采用在程序中事先指定边界的方法来去除这些空白单元，但在SPH中，边界处于复杂变化中，不可能通过事先指定.因此如何在复杂的非凸区域上进行Delaunay三角化并去除空白单元，从而为SPH算法提供后处理所需的网格，是非凸区域上SPH算法计算可视化急需解决的问题.

本文以凸区域上离散点集Delaunay三角化为基础，提出一种“单元称重”算法去除非凸区域上离散点集最小凸包三角化后形成的一些空白单元，准确得到了非凸区域SPH后处理所需要的网格，结合后处理软件Tecplot实现了非凸区域上SPH计算结果的可视化.

SPH作为一种基于物质描述的无网格算法，在自由表面流动等领域得到了非常好的应用［2，10］，以至在流体的计算机实时动画领域也开始采用该算法［11-12］，计算结果用很逼真的形象显示出来.动画领域对计算结果的处理，主要关注物质域表面的提取［10，13-14］及色泽处理，这些技术也是值得SPH算法后处理借鉴的.

1 SPH基本原理

SPH利用函数f与光滑函数w的加权积分来近似f与Dirac函数δ的卷积，从而可以近似得到f在计算域某点函数及空间导数的近似值［1］，离散后可得到

其中τj为粒子体积.利用它们，便可以将流体的动力学方程进行离散.由于本文侧重的是计算可视化的探讨，后文算例计算所需要的方程和参数会给出相应文献出处，这里不再赘述.本文中所有算例的光滑函数为分段五次样条函数［1］，初始光

2 物质域离散点集的Delaunay三角化及“单元称重”法

对于任意一个离散点集，都存在一个最小凸包包含这些点.对该最小凸包进行Delaunay三角化，可形成由该离散点集的点作为节点的三角网格(单元)集，而且每1个三角单元的外接圆都不包含其它的节点;而由这些外接圆的圆心形成的多边形称为Voronoi图［15］.离散点集的Delaunay三角化和Voronoi图是计算几何中最重要的算法，这里不再赘引.

2.1 凸区域离散点集上的Delaunay三角化

采用SPH粒子离散物质域之后，可以将SPH粒子集合的最小凸包进行Delaunay三角化，从而得到非结构化网格(或单元集)，并将SPH粒子作为单元节点，然后利用有限元插值方法，得到单元内任何一点处的物理量的值.对于凸区域上的离散点集(粒子集合)，也可借助Matlab中Delaunay算法来对其最小凸包进行三角化［7-8］.凸区域上离散点集的最小凸包的Delaunay三角化较为简单，且是非凸区域上离散点集三角化的基础.

2.2 非凸区域离散点集上的Delaunay三角化及“单元称重”算法

对非凸区域点集直接进行Delaunay三角化会得到一些空白单元，这些单元包含任何物质.如图1(a)所示，正方形区域的中心空洞中并没有任何物质点，但是该点集最小凸包是包含所有粒子的正方形区域，因此该凸包也包含该正方形点集中心的空洞区域;当对该凸包进行三角化时，该空洞也同时被三角化了，如图1(b)所示.

图1 非凸区域的点集及其三角化后形成的空白单元

去除这些空白区域，文献［6］通过事先指定边界—比如指定空洞的周界，然后除掉处于该周界内的单元，剩下的单元就可以作为后处理的网格了.该方法对于静变形，或材料(或计算域)不破碎、不卷曲的情况是可考虑应用的.但是在SPH领域，特别是波浪的飞溅、卷曲问题中，该方法无法应用.因此，本文提出一种“单元称重”算法，可以去除各种复杂非凸区域上的离散点集三角化后形成的空白单元，并且其可行性在具有材料飞溅和粒子不均匀性的高速撞击问题中得到了验证.

2.2.1 空白单元过滤的“单元称重”法

“单元称重”法的基本思想是，不属于物质域的空白单元，其质量比其它单元的质量或其节点(粒子)的质量要小很多，这类单元可以从单元集中去除.因此该算法的核心问题是，如何寻找1个方法来计算单元的质量，并提出1个比较的标准将它与周围的粒子的质量进行比较，然后将其去除.解决这个问题，可以考虑将SPH求和近似算法式(1)结合三角单元的外心(也就是上文提到的相应三角网格的Voronoi多边形的顶点)，具体算法如下.

每1个单元的位置用该单元的外心的位置re来表示，其中re为外心的矢径，下标“e”表示单元的序号.而用外心处的加权质量来表征该单元的质量.因为该单元的外心正好是某Voronoi多变形的顶点，因此该单元的外接圆中将不含任何其它单元的其它节点，这确保了对该单元质量影响权重最大的是属于该单元本身的节点;而该单元的所有节点到单元外心的距离都相等(都等于外径)，确保了它的节点对其外心的影响权重相等.这些性质非常重要，这也是Delaunay三角化的优点.

要利用SPH求和近似得到外心处的加权质量，需要知道该单元外心处的光滑长度值.对任1个单元，可以取光滑长度为该单元所有节点的光滑长度的代数平均，即.其中j为该单元上节点的序号，h为光滑长度;利用式(1)可得单元加权质量me为

其中k为该单元外心的紧支域内支持粒子的序号.这些粒子可能不仅仅包含该单元本身上的节点，还可能包括周围其它单元的节点;也有可能该外心的紧支域内不含任何支持粒子，甚至不含单元本身的节点—此时其加权质量为0.

通过下式—即“单元称重”算法的定义式，来过滤掉不属于物质域内的空白单元

因此式(4)可以表达为

式(4)可理解为，若计算所得单元的加权质量me小于该单元节点质量的平均值¯m的一定量，那么该单元可以被过滤掉;式(6)可以理解为，若一个单元的加权质量因子βe小于事先设定的标准值βs，则这个单元应该被过滤掉.对于图1(b)中非凸区域的离散点集三角化后形成的空白单元，利用“单元称重”法就可以去除，如图2所示.

由图2可发现，当βs不同时，可以过滤掉一部分共同的空白单元.βs较小时，过滤后周界显得光滑，如图2(a);βs变大时，周界显得较不光滑，如图2(c).这是因为βs较小时，只有那些加权质量me特别小的单元被过滤掉，被过滤掉的单元较少，所以边界处显得光滑;而βs变大时，有些加权质量不是特别小的单元也被过滤掉;βs越大，被过滤掉的单元越多.注意到，式(3)中粒子对单元外心的权重wek与单元的外径成反比关系，所以外径越大的单元(也就是单元空间尺寸越大)，其加权质量自然越小.

特别是当临界加权质量因子βs=0时，只有那些加权质量为0的单元被过滤掉，此时表现为因为单元外径太大，该单元上的3个节点都对单元的外心的权重wek为0.考虑到SPH中粒子间存在相互作用的概念是基于wek不为0，所以可以将这些节点对其单元的加权求和称为对外心的作用.注意到SPH中粒子间的互相作用在数学上表现为wek不为0，但在物理意义上表现为粒子间的作用力不为0.而在连续介质力学中，力的产生是因为相邻的物质发生相对运动或变形，即力就是物质的相互作用.因此如果两点之间不互相作用，那么可以认为这两个点上的质量粒子没有相互接触.所以如果三角单元上的节点对其外心没有作用，即wek为0，那么可认为该单元上的粒子与其外心间没有质量接触;而没有质量接触，则表明该单元内部有质量间断，该单元就自然可被过滤掉.可见单元称重算法是一种符合物理意义的算法.而βs＞0的含义是把单元内粒子间互相作用极其微弱的单元也过滤掉.

图2 βs不同时，“单元称重”法去除空白单元的情况

下面对单元称重法为什么能过滤掉空白单元进行更深层次考察.主要出发点是考察空白单元的尺寸和加权质量因子βe间的关系.

2.2.2 “单元称重”法过滤空白单元的可行性分析

首先给出单元称重算法可行的直观解释.由图1可以发现空白单元的3个节点中只要有一对节点的连线跨越了空白区域，则这个连线的距离就会比较大.而该空白单元的3个节点中，只要有2个节点之间的连线距离比较大，那么外接圆半径就会很大，而外径很大时，根据光滑函数的衰减特性可知，这3个节点对该单元外心处的权重wek就会迅速减小并接近于0.而图2(a)中处于边界附近的小单元依然可以保存下来，原因也正在于此，即它们之间的连线并未完全跨越空白区域.

定义单元的特征长度为外接圆半径，并定义单元的特征尺度为特征长度与其光滑长度之比，记为γe.显然，只要某单元的特征长度接近于光滑函数wek的紧支域半径，那么单元的加权质量就为0，那些特别大的空白单元就立即被过滤掉;而那些单元特征长度大于光滑长度而接近于紧支域半径的，因为光滑函数的衰减性质，由式(3)所得的加权质量会非常小，于是由式(4)中βs的大小决定它们是否被过滤掉.

为了给出式(4)的定性分析，在一个等边三角形的单元上，根据式(3)和(4)，绘制出该单元本身加权质量因子βe随单元特征尺度变化的大致关系，如图3所示.由图3可知，随γe增长，曲线、横坐标线及γe竖线围成的曲边三角区域越来越小.特别是当尺度因子γe接近2.5时，也就是单元尺度接近于2.5倍的光滑长度时，此时βs≈0，即只有那些加权质量me接近于0的单元被过滤掉，剩下的大多数单元被保存下来.随着尺度因子γe的缩小，并接近1.5时，被过滤的单元数量几乎不变，如图中的最大曲边三角形，此时βs≈1，也就是那些单元加权质量me小于粒子平均质量¯m的单元基本上都被过滤了(即那些粒子对其外心作用比较弱的甚至稍强的也都被过略掉了)，保留下的单元其质量都大于其节点的平均质量.这可作为图2中随着临界βs增大，被过滤的单元几乎不变的原因，因为剩下的可被过滤的单元越来越少.

图3 单元加权质量随其特征尺度的变化关系

“单元称重”法中反映一个单元是否可被过滤的参数，即加权质量因子βe，它反映了一个单元内部所含质量的多少，是单元内部粒子之间作用强弱的表现.从以上分析中还可得知，它也是单元尺度γe小或大的间接体现.也就是说，本质上讲，式(6)是根据被过滤单元尺度γe来判定哪些大尺度的单元是可以被过滤掉的—这与SPH中粒子间是否是粒子作用对是由它们之间距离来判断的思想是一致的.

临界加权质量因子βe的合理范围，可以根据图3中的曲线来估计.一般地，如果尺度因子为2以上，即特征长度为2倍光滑长度以上的单元被过滤掉，那么根据图3，就可以得到此时βs≈0.15.

3 SPH后处理的有限元插值方法及数值算例

形成网格后，就可由网格节点上所存储的物理量的数值进行有限元插值，从而得到网格单元内各点处的相应数值.完成这个过程，可选择自主开发程序，也可利用现有软件，如Tecplot［16］.而要利用成熟的后处理软件Tecplot精确处理SPH这些数据，就需要将SPH数据转换到离散该物质域的有限元格式，从而实现有限元插值.图4为液滴问题［17］的几个时刻的处理结果.三维时的情是类似的.

非凸区域上的SPH算法后处理同样如此，但是需要将离散点集三角化之后，利用“单元称重”法将空白单元过滤掉，再将SPH数据结合过滤后的单元集用有限元格式导入到Tecplot中处理.若图2区域上各粒子上具有F=sin(x2+y2)cos(x2-y2)的物理量，则其利用以上SPH算法的后处理方法得到的曲面与等值线云图如图5所示.Tecplot功能远不限于此;它还具有微积分、切片、光照等功能.

对于大多数基于SPH的研究，本文方法已足够使用.特别是去除非凸区域上离散点集三角化的空白单元，“单元称重”法简洁且适应性强.值得注意的是，“单元称重”法不需要事先判断该点集是否属于非凸区域，而是直接进行单元过滤.因此只要程序中加入该算法，不但不对凸区域上点集三角化后形成的的单元有影响，而且可直接对非凸集三角化后的单元进行过滤.

图4 凸区域上利用有限元插值方法得到的SPH后处理结果

图5 非凸区域上利用有限元插值方法得到的SPH后处理结果

3.1 “单元称重”算法在粒子非均匀分布与材料飞溅情况的适应性

以上所计算的算例粒子分布都比较均匀，如图4～5.为了检验“单元称重”法对实际复杂情况的适用性，本文计算了一个高速撞击问题，在该算例中具有较严重的粒子分布不均匀性与材料飞溅问题.控制方程同文献［17］，计算域如图6(a)～(c)所示，粒子的间距为5.0×10-2m，而时间步长为5.0 μs，取人工粘性常数αII=1，βII=2［17］.球形液滴以1 000 m/s初始速度向左运动.2.5 ms和12.5 ms时粒子分布、三角化后的单元集以及经过加权质量因子βs=1.0e-3过滤后的单元集，如图6所示.

图6 不同时刻粒子分布及在其最小凸包上三角化形成的单元集和过滤后的情形，βs=1.0e－3

4 结论

本文提出一种“单元称重”算法，用于去除非凸区域上离散点(粒子)集的三角化后形成的空白单元.“单元称重”法只过滤非凸区域离散点集的三角化后形成的空白单元，而对凸区域上离散点集三角化形成的单元没有影响，该方法简洁且对复杂的非凸区域适应性较强.利用“单元称重”法过滤后的单元集，结合SPH计算所得数据，并将数据文件写成有限元格式，导入到Tecplot中进行直接后处理，获得传统科学计算可视化所需的计算图形，为SPH算法的后处理，提供了一种实用简洁的途径.该方法的适用性和可行性在具有材料飞溅和粒子不均性的数值算例中得到了验证.

［1］LIU G R，LIU M B.光滑粒子流体动力学—一种无网格粒子法［M］.韩旭，杨刚，强洪夫，译.长沙:湖南大学出版社，2005:1-200.

［2］OGER G，DORING M，ALESSANDRINI B，et al.Twodimensional SPH simulations of wedge water entries［J］.Journal of Computational Physics，2006，213(2): 803-822.

［3］MASSIDDA L.ARMANDO，a SPH code for CERN-some theory，a short tutorial，the code description and some examples［R］.［S.l.］:European Organization for Nuclear Research，2008.

［4］DANIEL P.Visualisation of SPH data using SPLASH-v1.9.1［EB/OL］.2007-9-12.http://www.astro. ex.ac.uk/people/dprice/splash/index.html.

［5］BIDDISCOMBE J，GRAHAM D，MARUZEWSKI P.Interactive visualization and exploration ofSPH data［C］//Second International Workshop of Smoothed Particle Hydrodynamics European Research Interest Community.Madrid:［s.n.］，2007.

［6］文建波，周进雄，张红艳，等.基于Delaunay三角化的无网格法计算结果后处理［J］.应用力学学报，2003，4(20):600-601.

［7］田仲可.基于面片句柄对象的无网格法可视化研究［J］.计算机工程与设计，2008，29(13):3513-3515.

［8］史宝军，袁明武，陈永强.无网格方法数值结果的可视化方法与实现［J］.工程力学，2004，21(6):51-55.

［9］ZHOU J X，WEN J B，ZHANG H Y，et al.A nodal integration and post-processing technique based on Voronoi diagram for Galerkin meshless methods［J］.Comput Methods Appl Mech Engrg，2003，192:3831-3843.

［10］BECKER M，TESCHNER M.Weakly compressible SPH for free surface flows［C］//Eurographics ACM SIGGRAPH Symposium on Computer Animation.San Diego:［s.n.］，2007:63-72.

［11］柳有权，刘学慧，朱红斌，等.基于物理的流体模拟动画综述［J］.计算机辅助设计与图形学学报，2005，17(12):2581-2589.

［12］KOUMOUTSAKOS P，COTTET G H，ROSSINELLI D. Flow simulations using particales bridging computer graphics and CFD［R］.http://www.cse-lab.ethz.ch/ teaching/classes/mulsup.html.

［13］ROSENTHAL P，ROSSWOG S，LINSEN L.Direct surface extraction from smoothed particle hydrodynamics simulation data［C］//Proceedings of the.4th High-End VisualizationWorkshop. Berlin:LehmannsMedia，2007:50-61.

［14］KIM J，CHA D，CHANG B，et al.Practical animation of turbulent splashing water［C］//In Proceedings of the 2006 ACM SIGGRAPH/Eurographics Symposium on Computer Animation.Eurographics:［s.n.］，2006: 335-344.

［15］周培德.计算几何—算法分析与设计［M］.北京:清华大学出版社，2000:1-100.

［16］TECPLOT.Tecplot 10 Release 6 on line help［EB/ OL］.www.tecplot.Com/support.

［17］FANG J，PARRIAUX A，RENTSCHLER M，et al. Improved SPH methods for simulating free surface flows of viscous fluids［J］.Applied Numerical Mathematics，2009，59(2):251-271.