一种考虑界面不连续的改进的有限粒子法

王 璐，杨 扬，徐 绯

(西北工业大学航空学院，陕西 西安 710072)

光滑粒子流体动力学(smoothed particle hydrodynamics，SPH)是一种典型的拉格朗日型无网格粒子法，具有概念简单、自适应性强等诸多优点，工程应用十分广泛[1-3]。但其发展至今仍然存在一定的数值缺陷，其中在粒子缺失较明显的模型边界或物理量变化较大的非连续界面附近表现出的核估计精度不足尤为显著[4-5]。研究发现，可以通过提高核函数的一致性提高边界区域的计算精度。为此，Liu等[6]在SPH的基础上通过重构核函数提出了再生核质点法(reproducing kernel particle method，RKPM)，而后Chen等[7]以泰勒级数展开为基础提出了一种对SPH估计式进行正则化的方法，即修正光滑粒子法(corrective smoothed particle method，CSPM)，并得到了成功应用[8]；随后， Liu等[9]基于泰勒级数展开提出了有限粒子法(finite particle method，FPM)，相比传统的SPH方法，FPM在核函数的选取上较宽松，边界区域计算精度高，高阶扩展性强。相比CSPM，FPM可以同时估计核函数及其导数值，一定程度上降低了CSPM中低阶导数用于高阶导数计算而产生的累积误差。然而FPM也存在一些固有缺陷：在处理不连续问题时，若与连续体内部处理方法一致，算法的准确性会大大降低；求解线性方程组一定程度上增加了计算复杂度，计算时间较长；系数矩阵是否奇异不可控，影响计算稳定性。因而使用FPM方法时需对界面进行更精细的处理[10-11]。基于以上缺陷，本文中针对界面不连续问题对FPM方法进行改进，提高界面处的精度和计算效率。

针对传统SPH方法，目前已发展了不同的界面处理方法[12-13]，其中较经典的是Liu等[14]基于CSPM提出的非连续SPH(discontinuous SPH，DSPH)方法，DSPH方法改善了传统SPH方法在非连续区域上的计算精度。本文中基于DSPH方法处理非连续区域的思想，从FPM出发，推导DSFPM(discontinuous special FPM)方程，在理论层面改进FPM算法在界面处的计算精度。

1 DSFPM

1.1 FPM基本理论

根据Liu等[9]基于泰勒展开和矩阵理论给出的FPM基本方程，经过具体推导，给出FPM中一维以及二维情形下函数及其一阶导数的粒子近似式：

(1)

(2)

式中:f(xi)为点xi处的函数值，fx(xi)为点xi处的函数导数值；类似地，f(xi,yi)、fx(xi,yi)和fy(xi,yi)分别为点(xi,yi)处的函数值及偏x和偏y导数；j为粒子i支持域内的粒子，Wij为基函数，Wij,x=∂Wij/∂x，Wij,y=∂Wij/∂y分别为基函数的偏导数，mj、ρj分别为粒子j的质量和密度，N为粒子i完整支持域内的粒子数。

从式(1)～(2)可以看出，式(1)～(2)等号右边第一项矩阵只和粒子i支持域内粒子的位置有关，与其函数值无关，粒子的分布决定了矩阵是否奇异和FPM计算的稳定性。

1.2 DSPH基本理论

不同于FPM在粒子i完整支持域内进行泰勒展开，DSPH基于界面处物理量的不连续性，在非连续区域两边分别进行泰勒展开。Xu等[15]和闫蕊等[16]在Liu等[14]提出的一维不连续SPH方法的基础上将DSPH扩展到多维情形，并且避免了多维情形下确定不连续位置，降低了计算的复杂性。式(3)～(4)为DSPH粒子估计式，如图1 所示，Ω为粒子的整个支持域，Ω1和Ω2为其子域，Ω1和Ω2之间不连续界面为Γ(Ω1∪Ω2∪Γ=Ω)。

图1 界面附近粒子支持域Fig.1 The supported domain of the particle near the interface

(3)

(4)

式中:fα(xi)=∂f(xi)/∂xα，α、β表示坐标方向(一维为x,二维为x、y)。方程(3)和(4)等号右边第一项与CSPM保持一致，第二项是考虑了不连续界面后产生的附加项，附加项中只提取了异侧粒子的位置、质量和密度信息，方程式很明显为隐式格式，保证了计算的稳定性但较耗时。

1.3 DSFPM推导

从式(1)可以看到，当遇到非连续物理时，FPM在完整支持域内求和的做法对于求解非连续问题的精度是远远不够的，本节将基于DSPH方法处理非连续问题的思想对FPM进行理论改进。

从方程式(3)出发，对其进行变形推导得到:

(5)

(6)

从式(5)～(6)可以看出DSPH公式经过变形后，本质是只在子域Ω1内对被估粒子进行一致化修正，这样被估粒子的物理信息就完全与异侧粒子分隔开。基于该想法，对FPM基本方程即式(1)～(2)进行下列修正[17]，得到DFPM(discontinuous FPM)的基本方程:

(7)

(8)

其中，DFPM只考虑了邻域内与粒子i同侧的粒子信息。该修正式既保证了高阶导数的计算精确性，也实现了不连续界面处粒子信息的准确估计。但是式(7)～(8)等号右边第一项矩阵仍不能准确保证计算的稳定性，并且矩阵的各项求解也使得计算过程变得繁琐费时，因此下面进一步对式(7)～(8)进行变形整理。

从式(7)出发，将式(7)整理为：

(9)

一维情形下，粒子所占体积为粒子间距dj，区域Ω1内的粒子数为N1，类似文献[18]的处理方法，在式(9)中，记：

(10)

(12)

式(9)变形为：

KDCf=KDF

(13)

图2 一维情形下DSFPM粒子选取方式Fig.2 Particle selection mode for 1-D DSFPM

图3 二维情形下DSFPM粒子选取方式Fig.3 Particle selection mode for 2-D DSFPM

可以看出，K只与核函数有关，D为粒子间距信息，C为区域Ω1内的粒子坐标信息，F反映了区域Ω1内的粒子函数值信息。此时，f为待求未知量。由于FPM中核函数选取自由，因此我们只考虑K为满秩的情形，如果在区域Ω1内只选取距离粒子i最近的2个粒子，如图2所示，红虚线为界面，此时，K为可逆方阵，D也为可逆方阵。在式(13)两边同时乘以D-1K-1，式(13)简化为：

Cf=F

(14)

(15)

类似地，二维情形下，如果在粒子(xi,yi)的支持域子域Ω1内选取距离粒子(xi,yi)最近的3个粒子，如图3所示，红虚线为界面，此时K和D亦均为可逆矩阵，同样可得到式(14)。具体表达式为：

(16)

式(15)～(16)即为DSFPM基本方程，相比于式(7)～(8)，DSFPM消除了核函数对于矩阵奇异性的影响，提高了计算的稳定性，计算过程大大简化。

2 DSFPM粒子近似精度分析

一般情况下，如果一种近似格式能够再生k阶多项式，那么则称该格式具有Ck阶连续性。类似地，我们定义如果近似格式能够再生k阶非连续多项式，那么该格式具有Ck阶非连续性。接下来，通过一系列数值算例对传统FPM、DSPH 和DSFPM等3种方法的粒子近似精度进行对比分析。在数值计算中，选用常用的三次B-样条函数作为光滑函数[19]。

2.1 C0阶非连续性分析

C0阶表征近似格式对于非连续常值函数的再生能力，考虑区间[0,1]上的非连续常值函数：

(17)

自x=0至x=1均匀分布21个粒子，粒子间距d=0.05，光滑长度h=1.2d。图4给出了FPM、DSPH和DSFPM等3种方法估计的函数值误差φ及其一阶导数误差ψ,红虚线为界面。从图4可以看出，DSPH和DSFPM均可以准确计算界面附近粒子函数及其导数值，FPM在界面处计算精度低。这说明FPM处理不连续问题有明显缺陷，而其余2种方法改善了界面处的精度缺陷。

图4 3种方法估计非连续常值函数误差及其导数误差Fig.4 Estimation errors for the discontinuous 0-order function and its derivative by three methods

2.2 C1阶非连续性分析

C1阶表征近似格式对于非连续线性函数的再生能力，考虑一维不连续函数：

(18)

自x=0至x=1均匀分布21个粒子，粒子间距及光滑长度不变。图5给出FPM、DSPH和DSFPM等3种方法近似的函数值误差φ及其一阶导数误差ψ。

从图5可以看出：(1)对线性函数值，FPM和DSPH均不能准确估计，FPM在界面(x=0.525)附近的误差值远大于DSPH相应的误差值。这说明DSPH在一定程度上改善了界面处的计算精度。而DSPH在边界处(x=0，x=1)不能准确估计函数值，这是由于边界粒子支持域内粒子信息缺失导致核函数估计精度不足。而DSFPM能够很好地估计界面和边界处线性函数值。(2)对线性函数一阶导数值，FPM在界面附近的计算误差较大，DSPH对边界和界面附近的函数值估计有误差导致DSPH估计一阶导数时仍不精确，这体现了DSPH利用低阶导数计算高阶导数带来的误差累积效应，DSFPM仍可以准确再生一阶导数值，大大改善了界面以及边界附近的计算精度。DSFPM不仅延续了FPM在边界处核估计的高精度性，也继承了DSPH对于非连续区域计算的修正。

图5 3种方法估计非连续一次函数误差及其导数误差Fig.5 Estimation errors for the discontinuous 1-order function and its derivative by three methods

2.3 C2阶非连续性分析

这一部分将验证DSFPM对二次函数的再生能力。考虑二维情形下，区域 [0,1]×[0,1]上的非连续二次函数：

(19)

在二维区域内均匀分布441个粒子，光滑长度h=1.2d。图6为估计二次函数及其一阶偏导(由于界面垂直于y轴,因此只关注偏y导数误差)的误差。从图6可以看出，DSFPM不能对二次函数进行准确再生。这说明DSFPM只能对一次非连续函数进行准确估计，具有C1阶非连续性。这是因为DSFPM推导过程中利用的FPM基本方程只进行了一阶泰勒展开，并未考虑更高阶项。如果将f(x,y)泰勒展开到更高阶项，那么可以得到更高阶精度的DSFPM，这也体现了DSFPM的高阶可扩展性。

图6 DSFPM计算二次非连续函数本身和其偏y导数误差分布Fig.6 Estimation error distribution for the discontinuous 2-order function and its partial y derivative by DSFPM

3 碰撞算例模拟

为了进一步验证DSFPM的工程实用我们对碰撞算例进行模拟。由于DSFPM所选取参与计算的粒子数有限，且参与计算粒子的构型不能退化为直线，因此对于大变形碰撞问题，DSFPM可能会产生精度不足或不稳定现象。为了避免DSFPM处理大变形问题的缺陷，采用如图7所示的算法流程图。一方面，在变形初始阶段仍使用DSFPM，若粒子应变达到一定阈值时，DSFPM将自动转化为DFPM，自适应地实现大变形碰撞问题；另一方面，当粒子的应变未达到设定的阈值，而此时DSFPM中所选取的有限数目粒子的构型不能满足稳定性，算法也转化为DFPM算法。这节将对小变形和大变形问题分别进行模拟。

图7 处理碰撞冲击问题的算法框图Fig.7 Algorithm diagram of the collision problem

3.1 离散公式

弹性体运动控制方程为：

(20)

式中：ρ、p、x和u分别为弹性体的密度、各向同性压力、位移和速度。σ、g、c和ρ0分别为应力、重力加速度、声速和参考密度。利用SPH方法离散后，用α、β表示坐标方向，控制方程(20)用指标法可写为：

(21)

σαβ=Sαβ-δαβp

(22)

(23)

(24)

(25)

3.2 小变形碰撞

这部分以铝块碰撞为例对小变形问题进行模拟，模型如图8所示，两铝块长20 mm,高10 mm,初始时左侧铝块以20 m/s的初始速度向右运动，右侧铝块静止，M点为右侧铝块碰撞界面处的中点。铝参考密度为2 785 kg/m3，声速取为5 328 m/s，杨氏模量为72 GPa，泊松比为0.3。粒子间距d=0.5 mm，光滑长度h=1.2d。

图8 铝块碰撞模型Fig.8 The model for the aluminum block collision

分别将式(4)、(6)、(16)中的f(xi,yi)替换为vα(xi,yi)，实现FPM、DSPH和DSFPM等3种方法对于速度梯度的估计。为了保证粒子相互作用的对称性，动量方程仍都采用式(21)进行计算。式(6)、(16)中只考虑了与粒子i同侧的粒子信息，但是由于2个铝块之间存在能量传递，因此在DSPH和DSFPM计算中加入接触力作为2个铝块相互作用的体现，当粒子i支持域内搜索到异侧粒子s时，粒子s对粒子i施加一定大小的作用力，接触力的具体表达式[20]为：

(26)

式中：ci为粒子i的声速，mi和ms分别为粒子i和粒子s的质量，W(ris)为核函数，ris=ri-rs。该接触力直接加到动量方程中，且满足mif(|ris|)=-msf(|rsi|)。

3.2.1速度分布

图9 3种方法模拟得到的铝块速度随时间的变化曲线Fig.9 Velocity-time cuvers of aluminum blocks by three methods

图10 有限元模拟得到的铝块速度随时间的变化曲线Fig.10 Velocity-time cuvers of aluminum blocks by finite element method

图11给出了FPM、DSPH和DSFPM模拟得到的在铝块碰撞过程中某一时刻的速度分布图。从图11可以看出，该时刻两侧铝块正在进行速度交换，DSFPM和DSPH模拟得到的速度连续且分布均匀，但FPM模拟得到的速度在2个铝块界面附近出现很明显的不连续，震荡较明显。这由于铝块速度相差较大，FPM的全支持域搜索会导致粒子速度梯度估计过大，造成FPM在界面附近速度梯度计算精度不足。

图11 3种方法模拟得到的铝块碰撞过程中的速度分布Fig.11 Velocity distribution of the aluminum blocks by three methods

3.2.2应力分布

2个铝块发生碰撞时，在铝块内部会产生冲击波，压缩(拉伸)波在弹性体内传递，图12给出碰撞过程中某一时刻3种方法模拟的x方向应力分布，该时刻2铝块应力为负，受到挤压。从图12可以看出:(1)FPM模拟的界面附近出现很严重的正负应力交错现象，说明FPM在界面处的精度最低。(2)DSPH虽然大大改善了界面处应力交错现象，但其在模型四周边界处出现了应力不均匀的缺陷，体现出DSPH在边界处核近似精度不足。(3)DSFPM模拟得到的应力不论是在碰撞界面还是模型边界都很均匀，说明了DSFPM处理不连续界面和边界的优势。

图12 3种方法模拟得到铝块碰撞过程中的x方向应力分布Fig.12 x-direction total stress distribution of the aluminum blocks simulated by three methods

图13 3种方法模拟得到点M处x方向应力随时间的变化曲线Fig.13 x-directional total stress changes at point M simulated by three methods

图13给出3种方法模拟得到右侧铝块点M的x方向应力随时间变化曲线。从图13可以看到，该点应力约在0.04 ms达到最低值，压缩应力达到最大，随后，应力迅速上升。最后DSFPM和DSPH模拟的应力便保持稳定，而FPM计算的应力迅速上升后达到正值，出现应力震荡现象，和图12的现象一致。

3.2.3计算效率

DSFPM计算耗时最短，为99.49 s；DSPH和FPM计算所用时间相当，分别为171.57、173.89 s。这是由于在粒子法计算过程中，粒子邻域搜索这一环节占计算总时间比重大，而根据框图7，DSFPM是基于初始构型对有限个粒子进行选取的，在计算过程中不需要进行邻域搜索，这大大降低了计算时间。DSPH和FPM都需进行邻域搜索，两者仅在速度梯度计算方法上有所差异，因此计算时间相差很小。对该算例，DSFPM可以节省约40%的计算时间。

3.3 大变形碰撞

基于框图7算法，对如下大变形问题进行模拟，如图14所示，左侧物块以50 m/s撞击右侧板，板两端固定，左侧物块长20 mm，高5 mm，密度为7 200 kg/m3，弹性模量E1=270 GPa，泊松比为0.3。右侧板长2.5 mm，高20 mm，密度为2 785 kg/m3，弹性模量E2为72 GPa，泊松比为0.3。粒子间距d=0.2 mm，光滑长度h=1.2d。

图14 物块冲击模型Fig.14 The model for block collosion

该算例中选取的应变阈值为0.01，当粒子应变超过该值时，DSFPM算法将被转化为DFPM算法。图15给出了不同时刻框图7算法模拟的速度分布。从图15可以看出，右侧板从中部开始变形，两端固支附近弯矩最大，最后板两端和中间附近出现裂口。该算法模拟的速度分布较均匀，体现了本文提出的算法在界面处计算精度的优势。

若其他条件不变，减小右侧板弹性模量E2，便可以得到更大变形的结果。图16为右侧板弹性模量为72 MPa时不同时刻的碰撞速度分布图。从图16可以看出，当变形增大时，框图7算法仍然能够模拟出连续均匀的速度。由此得知，DSFPM结合DFPM后，它不仅避免了DSFPM在处理大变形问题中参与粒子数较少的问题，也规避了DSFPM中有限粒子构型可能带来的稳定性问题，是处理大变形问题的有效方法。

图16 E2=72 MPa时，不同时刻框图7算法模拟的速度云图Fig.16 Velocity distribution of the blocks simulted by the algorithm in Fig.7 at different times when E2=72 MPa

4 结 论

针对FPM处理非连续物理场的缺陷，基于DSPH方法提出了DSFPM。通过对DSFPM进行近似精度分析发现DSFPM具有二阶精度，可以准确再生一次非连续函数。在铝块冲击碰撞小变形问题中，通过对比分析DSFPM、DSPH和FPM等3种方法模拟的铝块速度、应力以及计算时间，验证了DSFPM在不连续界面处的计算优势。DSFPM大大改善了界面区域的计算精度，提高了计算效率，是FPM方法的有效改进。由于DSFPM算法本身选取粒子的特殊性，DSFPM在大变形问题中具有一定局限性，为了使改进算法仍然适用于大变形，我们提出了DSFPM自适应转化为DFPM的方法，并通过冲击大变形算例验证了该算法的可行性。