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

王璐; 杨扬; 徐绯

doi:10.11883/bzycj-2017-0390

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

doi: 10.11883/bzycj-2017-0390

西北工业大学航空学院, 陕西西安 710072

基金项目:

国家自然科学基金项目 11272266

国家自然科学基金项目 11702220

航空科学基金项目 2016ZD53038

中央高校基本科研业务费专项资金项目 3102017zy066

详细信息

作者简介:
王璐(1995-), 女, 博士研究生, wanglulu@mail.nwpu.edu.cn

通讯作者:
徐绯(1970-), 女, 博士, 教授, 博士生导师, xufei@nwpu.edu.cn

中图分类号: O382
计量
- 文章访问数: 6290
- HTML全文浏览量: 2192
- PDF下载量: 56
- 被引次数: 2
出版历程
- 收稿日期: 2017-10-30
- 修回日期: 2018-01-10
- 刊出日期: 2019-02-05

An improved finite particle method for discontinuous interface problems

School of Aeronautics, Northwestern Polytechnical University, Xi'an 710072, Shaanxi, China

摘要

摘要: 有限粒子法（finite particle method，FPM）作为SPH（smoothed particle hydrodynamics）方法的重要改进，有效提高了边界区域粒子的近似精度，但是当FPM处理多物理场时，在不连续界面附近的计算精度会大大降低，并且FPM必须满足的矩阵非奇异性也提高了对界面处理的要求。本文中基于DSPH（discontinuous SPH）方法，提出了一种考虑界面不连续的改进FPM—DSFPM（discontinuous special FPM）法，旨在改善FPM在界面不连续处的计算精度，从而进一步提高其计算效率和稳定性。首先，分析了DSFPM的核近似精度。其次，根据不同的工程问题，给出DSFPM处理小变形和大变形问题的算法流程。利用DSFPM、DSPH和FPM等3种方法对弹性铝块小变形碰撞冲击算例进行了模拟，通过对比分析铝块的速度和应力以及计算时间验证了DSFPM算法在非连续界面处计算精度和计算效率的优势。最后，通过结合DSFPM和DFPM（discontinuous FPM）实现了对于大变形问题的模拟。
- FPM /
- 界面 /
- DSFPM /
- 计算精度 /
- 计算效率
Abstract: The finite particle method(FPM) is an important improvement for the smoothed particle hydrodynamics(SPH) method, which effectively improves the calculation accuracy of boundary particles. However, when the discontinuous physical field is solved by the FPM, the accuracy in the vicinity of the discontinuous interface is greatly reduced, and the non-singularity of the matrix must be satisfied in the FPM, which requires an elaborate handling of the interface. Based on the discontinuous SPH(DSPH) method, this paper proposed an improved FPM-discontinuous special FPM(DSFPM), which considers the discontinuous interface, aiming to improve the computational accuracy at the interface and further improve the efficiency and stability of the FPM. In this paper, the estimation accuracy of the DSFPM was analyzed firstly, and then the algorithm flow diagram of the DSFPM to deal with the small deformation and large deformation problems was demonstrated. Next, the DSFPM, DSPH and FPM were used to simulate the small deformation problem-elastic aluminum blocks impact. By comparing the velocity and stress of the aluminum blocks and computational time, we verified the accuracy and computational efficiency of the DSFPM. Finally, the simulation of the large deformation problem was realized by a combining method with the DSFPM and DFPM.
- FPM /
- interface /
- DSFPM /
- accuracy /
- computational cost

HTML全文

光滑粒子流体动力学(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中一维以及二维情形下函数及其一阶导数的粒子近似式：

$\left[ {\begin{array}{*{20}{c}} {f\left( {{x_i}} \right)}\\ {{f_x}\left( {{x_i}} \right)} \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} {\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}} }&{\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}\left( {{x_j} - {x_i}} \right)} }\\ {\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,x}}} }&{\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,x}}\left( {{x_j} - {x_i}} \right)} } \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} {\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{x_j}} \right){W_{ij}}} }\\ {\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{x_j}} \right){W_{ij,x}}} } \end{array}} \right]$

(1)

$\begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {f\left( {{x_i},{y_i}} \right)}\\ {{f_x}\left( {{x_i},{y_i}} \right)}\\ {{f_y}\left( {{x_i},{y_i}} \right)} \end{array}} \right] = {{\left[ {\begin{array}{*{20}{c}} {\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}} }&{\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}\left( {{x_j} - {x_i}} \right)} }&{\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}\left( {{y_j} - {y_i}} \right)} }\\ {\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,x}}} }&{\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,x}}\left( {{x_j} - {x_i}} \right)} }&{\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,x}}\left( {{y_j} - {y_i}} \right)} }\\ {\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,y}}} }&{\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,y}}\left( {{x_j} - {x_i}} \right)} }&{\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,y}}\left( {{y_j} - {y_i}} \right)} } \end{array}} \right]}^{ - 1}} \times }\\ {\left[ {\begin{array}{*{20}{c}} {\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{x_j},{y_j}} \right){W_{ij}}} }\\ {\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{x_j},{y_j}} \right){W_{ij,x}}} }\\ {\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{x_j},{y_j}} \right){W_{ij,y}}} } \end{array}} \right]} \end{array}$

(2)

式中:f(x_i)为点x_i处的函数值，f_x(x_i)为点x_i处的函数导数值；类似地，f(x_i, y_i)、f_x(x_i, y_i)和f_y(x_i, y_i)分别为点(x_i, y_i)处的函数值及偏x和偏y导数；j为粒子i支持域内的粒子，W_ij为基函数，W_{ij, x}=∂W_ij/∂x，W_{ij, y}=∂W_ij/∂y分别为基函数的偏导数，m_j、ρ_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所示，Ω为粒子的整个支持域，Ω₁和Ω₂为其子域，Ω₁和Ω₂之间不连续界面为Γ(Ω₁∪Ω₂∪Γ=Ω)。

$f\left( {{\mathit{\boldsymbol{x}}_i}} \right) = \frac{{\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{\mathit{\boldsymbol{x}}_j}} \right){W_{ij}}} }}{{\sum\limits_{j \in \mathit{\Omega }} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}} }} + \frac{{f\left( {{\mathit{\boldsymbol{x}}_i}} \right)\sum\limits_{j \in {\mathit{\Omega }_2}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}} }}{{\sum\limits_{j \in \mathit{\Omega }} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}} }}$

(3)

${f_\alpha }\left( {{\mathit{\boldsymbol{x}}_i}} \right) = \frac{{\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}\left( {f\left( {{\mathit{\boldsymbol{x}}_j}} \right) - f\left( {{\mathit{\boldsymbol{x}}_i}} \right)} \right){W_{ij,\beta }}} }}{{\sum\limits_{j \in \mathit{\Omega }} {\frac{{{m_j}}}{{{\rho _j}}}\left( {x_j^\alpha - x_i^\alpha } \right){W_{ij,\beta }}} }} + \frac{{{f_\alpha }\left( {{\mathit{\boldsymbol{x}}_i}} \right)\sum\limits_{j \in {\mathit{\Omega }_2}} {\frac{{{m_j}}}{{{\rho _j}}}\left( {x_j^\alpha - x_i^\alpha } \right){W_{ij,\beta }}} }}{{\sum\limits_{j \in \mathit{\Omega }} {\frac{{{m_j}}}{{{\rho _j}}}\left( {x_j^\alpha - x_i^\alpha } \right){W_{ij,\beta }}} }}$

(4)

图 1 界面附近粒子支持域

Figure 1. The supported domain of the particle near the interface

下载: 全尺寸图片幻灯片

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

1.3 DSFPM推导

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

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

(5)

(6)

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

$\left[ {\begin{array}{*{20}{c}} {f\left( {{x_i}} \right)}\\ {{f_x}\left( {{x_i}} \right)} \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}} }&{\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}\left( {{x_j} - {x_i}} \right){W_{ij}}} }\\ {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,x}}} }&{\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}\left( {{x_j} - {x_i}} \right){W_{ij,x}}} } \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{x_j}} \right){W_{ij}}} }\\ {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{x_j}} \right){W_{ij,x}}} } \end{array}} \right]$

(7)

$\begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {f\left( {{x_i},{y_i}} \right)}\\ {{f_x}\left( {{x_i},{y_i}} \right)}\\ {{f_y}\left( {{x_i},{y_i}} \right)} \end{array}} \right] = {{\left[ {\begin{array}{*{20}{c}} {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}} }&{\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}\left( {{x_j} - {x_i}} \right)} }&{\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}\left( {{y_j} - {y_i}} \right)} }\\ {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,x}}} }&{\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,x}}\left( {{x_j} - {x_i}} \right)} }&{\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,x}}\left( {{y_j} - {y_i}} \right)} }\\ {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,y}}} }&{\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,y}}\left( {{x_j} - {x_i}} \right)} }&{\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,y}}\left( {{y_j} - {y_i}} \right)} } \end{array}} \right]}^{ - 1}} \times }\\ {\left[ {\begin{array}{*{20}{c}} {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{x_j},{y_j}} \right){W_{ij}}} }\\ {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{x_j},{y_j}} \right){W_{ij,x}}} }\\ {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{x_j},{y_j}} \right){W_{ij,y}}} } \end{array}} \right]} \end{array}$

(8)

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

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

$\left[ {\begin{array}{*{20}{c}} {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}} }&{\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}\left( {{x_j} - {x_i}} \right){W_{ij}}} }\\ {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,x}}} }&{\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}\left( {{x_j} - {x_i}} \right){W_{ij,x}}} } \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {f\left( {{x_i}} \right)}\\ {{f_x}\left( {{x_i}} \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{x_j}} \right){W_{ij}}} }\\ {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{x_j}} \right){W_{ij,x}}} } \end{array}} \right]$

(9)

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

$\begin{array}{*{20}{c}} {\mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}} }&{\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}\left( {{x_j} - {x_i}} \right){W_{ij}}} }\\ {\sum\limits_{j \in {\Omega _1}} {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij,x}}} }&{\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}\left( {{x_j} - {x_i}} \right){W_{ij,x}}} } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{W_{i1}}}&{{W_{i2}}}& \cdots &{{W_{i{N_1}}}}\\ {{W_{i1,x}}}&{{W_{i2,x}}}& \cdots &{{W_{i{N_1},x}}} \end{array}} \right] \times }\\ {\left[ {\begin{array}{*{20}{c}} {{d_1}}&{}&{}&{}\\ {}&{{d_2}}&{}&{}\\ {}&{}& \ddots &{}\\ {}&{}&{}&{{d_{{N_1}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 1&{{x_1} - {x_i}}\\ 1&{{x_2} - {x_i}}\\ \vdots&\vdots \\ 1&{{x_{{N_1}}} - {x_i}} \end{array}} \right] = \mathit{\boldsymbol{KDC}}} \end{array}$

(10)

$\mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{x_j}} \right){W_{ij}}} }\\ {\sum\limits_{j \in {\mathit{\Omega }_1}} {\frac{{{m_j}}}{{{\rho _j}}}f\left( {{x_j}} \right){W_{ij,x}}} } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{W_{i1}}}&{{W_{i2}}}& \cdots &{{W_{i{N_1}}}}\\ {{W_{i1,x}}}&{{W_{i2,x}}}& \cdots &{{W_{i{N_1},x}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{d_1}}&{}&{}&{}\\ {}&{{d_2}}&{}&{}\\ {}&{}& \ddots &{}\\ {}&{}&{}&{{d_{{N_1}}}} \end{array}} \right]\\ \left[ {\begin{array}{*{20}{c}} {f\left( {{x_1}} \right)}\\ {f\left( {{x_2}} \right)}\\ \vdots \\ {f\left( {{x_{{N_1}}}} \right)} \end{array}} \right] = \mathit{\boldsymbol{KDF}}$

(11)

$\mathit{\boldsymbol{f}} = \left[ {\begin{array}{*{20}{c}} {f\left( {{x_i}} \right)}\\ {{f_x}\left( {{x_i}} \right)} \end{array}} \right]$

(12)

式(9)变形为：

$\mathit{\boldsymbol{KDCf}} = \mathit{\boldsymbol{KDF}}$

(13)

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

$\mathit{\boldsymbol{Cf}} = \mathit{\boldsymbol{F}}$

(14)

$\left[ {\begin{array}{*{20}{c}} 1&{{x_1} - {x_i}}\\ 1&{{x_2} - {x_i}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {f\left( {{x_i}} \right)}\\ {{f_x}\left( {{x_i}} \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {f\left( {{x_1}} \right)}\\ {f\left( {{x_2}} \right)} \end{array}} \right]$

(15)

图 2 一维情形下DSFPM粒子选取方式

Figure 2. Particle selection mode for 1-D DSFPM

下载: 全尺寸图片幻灯片

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

$\begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} 1&{{x_1} - {x_i}}&{{y_1} - {y_i}}\\ 1&{{x_2} - {x_i}}&{{y_2} - {y_i}}\\ 1&{{x_3} - {x_i}}&{{y_3} - {y_i}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {f\left( {{x_i},{y_i}} \right)}\\ {{f_x}\left( {{x_i},{y_i}} \right)}\\ {{f_y}\left( {{x_i},{y_i}} \right)} \end{array}} \right]}\\ { = \left[ {\begin{array}{*{20}{c}} {f\left( {{x_1},{y_1}} \right)}\\ {f\left( {{x_2},{y_2}} \right)}\\ {f\left( {{x_3},{y_3}} \right)} \end{array}} \right]} \end{array}$

(16)

图 3 二维情形下DSFPM粒子选取方式

Figure 3. Particle selection mode for 2-D DSFPM

下载: 全尺寸图片幻灯片

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

2. DSFPM粒子近似精度分析

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

2.1 C⁰阶非连续性分析

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

$f\left( x \right) = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} 2\\ 5 \end{array}&\begin{array}{l} x \le 0.525\\ x < 0.525 \end{array} \end{array}} \right.$

(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种方法估计非连续常值函数误差及其导数误差

Figure 4. Estimation errors for the discontinuous 0-order function and its derivative by three methods

下载: 全尺寸图片幻灯片

2.2 C¹阶非连续性分析

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

$f\left( x \right) = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} 2x\\ 5x + 5 \end{array}&\begin{array}{l} x \le 0.525\\ x > 0.525 \end{array} \end{array}} \right.$

(18)

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

图 5 3种方法估计非连续一次函数误差及其导数误差

Figure 5. Estimation errors for the discontinuous 1-order function and its derivative by three methods

下载: 全尺寸图片幻灯片

从图 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对于非连续区域计算的修正。

2.3 C²阶非连续性分析

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

$f\left( {x,y} \right) = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} 2xy + {x^2}\\ xy + {y^2} + 5 \end{array}&\begin{array}{l} y \le 0.525\\ y > 0.525 \end{array} \end{array}} \right.$

(19)

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

图 6 DSFPM计算二次非连续函数本身和其偏y导数误差分布

Figure 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 处理碰撞冲击问题的算法框图

Figure 7. Algorithm diagram of the collision problem

下载: 全尺寸图片幻灯片

3.1 离散公式

弹性体运动控制方程为：

$\left\{ \begin{array}{l} \frac{{{\rm{d}}\rho }}{{{\rm{d}}t}} = - \rho {\rm{div}}\left( \mathit{\boldsymbol{u}} \right)\\ \rho \frac{{{\rm{d}}\mathit{\boldsymbol{u}}}}{{{\rm{d}}t}} = \frac{{\partial \mathit{\boldsymbol{\sigma }}}}{{\partial \mathit{\boldsymbol{x}}}} + \mathit{\boldsymbol{g}}\\ p = {c^2}\left( {\rho - {\rho _0}} \right)\\ \frac{{{\rm{d}}\mathit{\boldsymbol{x}}}}{{{\rm{d}}t}} = \mathit{\boldsymbol{u}} \end{array} \right.$

(20)

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

$\left\{ \begin{array}{l} \frac{{{\rm{d}}{\rho _i}}}{{{\rm{d}}t}} = {\rho _i}\sum\limits_j {\frac{{{m_j}}}{{{\rho _j}}}\left( {u_i^\alpha - u_j^\alpha } \right){W_{ij,\alpha }}} \\ \frac{{{\rm{d}}u_i^\alpha }}{{{\rm{d}}t}} = \sum\limits_j {{m_j}\left( {\frac{{\sigma _i^{\alpha \beta }}}{{\rho _i^2}} + \frac{{\sigma _j^{\alpha \beta }}}{{\rho _j^2}} + {\mathit{\Pi }_{ij}}{\delta ^{\alpha \beta }}} \right){W_{ij,\beta }}} + {F^\alpha }\\ {p_i} = {c^2}\left( {{\rho _i} - {\rho _{0i}}} \right)\\ \frac{{{\rm{d}}x_i^\alpha }}{{{\rm{d}}t}} = u_i^\alpha \end{array} \right.$

(21)

式中：j为粒子i支持域内的粒子，W_{ij, α}和W_{ij, β}为核函数在α、β方向上的偏导数，δ^αβ为狄拉克δ函数，F^α为粒子所受的外力，Π_ij为人工黏性；σ^αβ为粒子总应力，包含各项同性压力和偏应力(见式(22)); 偏应力率由式(23)给出，其中，μ为剪切模量， ${\dot \varepsilon ^{\alpha \beta }}$ 为应变率，R^αβ为旋转率张量，根据连续介质力学， ${\dot \varepsilon ^{\alpha \beta }}$ 、R^αβ可由式(24)~(25)得到，式中均包含速度梯度。由此可以看出，速度梯度对连续方程和加速度方程的计算起着决定性作用，准确地得到界面附近的速度梯度是模拟弹性体的关键。

${\sigma ^{\alpha \beta }} = {S^{\alpha \beta }} - {\delta ^{\alpha \beta }}p$

(22)

$\frac{{{\rm{d}}{S^{\alpha \beta }}}}{{{\rm{d}}t}} = 2\mu \left( {{{\dot \varepsilon }^{\alpha \beta }} - \frac{1}{3}{\delta ^{\alpha \beta }}{{\dot \varepsilon }^{\alpha \beta }}} \right) + {S^{\alpha \gamma }}{R^{\beta \gamma }} + {R^{\alpha \gamma }}{S^{\gamma \beta }}$

(23)

${{\dot \varepsilon }^{\alpha \beta }} = \frac{1}{2}\left( {\frac{{\partial {v^\alpha }}}{{\partial {x^\beta }}} + \frac{{\partial {v^\beta }}}{{\partial {x^\alpha }}}} \right)$

(24)

${R^{\alpha \beta }} = \frac{1}{2}\left( {\frac{{\partial {v^\alpha }}}{{\partial {x^\beta }}} - \frac{{\partial {v^\beta }}}{{\partial {x^\alpha }}}} \right)$

(25)

3.2 小变形碰撞

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

图 8 铝块碰撞模型

Figure 8. The model for the aluminum block collision

下载: 全尺寸图片幻灯片

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

$f\left( {\left| {{\mathit{\boldsymbol{r}}_{is}}} \right|} \right) = 0.01c_i^2\frac{{{\mathit{\boldsymbol{r}}_{is}}}}{{{{\left| {{\mathit{\boldsymbol{r}}_{is}}} \right|}^2}}}W\left( {{\mathit{\boldsymbol{r}}_{is}}} \right)\frac{{2{m_s}}}{{{m_i} + {m_s}}}$

(26)

式中：c_i为粒子i的声速，m_i和m_s分别为粒子i和粒子s的质量，W(r_is)为核函数，r_is=r_i-r_s。该接触力直接加到动量方程中，且满足m_if(|r_is|)=-m_sf(|r_si|)。

3.2.1 速度分布

理论上讲，2个铝块发生碰撞之后，两者的动能相互传递，左侧铝块静止，而右侧铝块以20 m/s向右运动。给出了采用3种方法计算得到的左右铝块平均速度曲线。采用DSFPM、DSPH、FPM模拟得到的右侧铝块在完成碰撞后速度分别为19.989 77、19.918 44、14.01 m/s，可以看出FPM误差较大，且碰撞后较长时间才达到稳定状态。根据动量定理：P= $\int {\mathit{\boldsymbol{F}}{\rm{d}}\mathit{t}}$ =mΔv，提取DSFPM和DSPH右侧铝块每一时刻所受的冲击力，经过积分运算，得到P_DSFPM=11.13 kg·(m/s)，P_DSPH=10.95 kg·(m/s)，而理想状态下右侧铝块动量增量mΔv=11.14 kg·(m/s)，这说明DSFPM更好地实现了速度交换过程。我们也对比了同等模型下，采用有限元方法所得到的速度变化曲线，如图 10所示，有限元模拟结果在碰撞后产生的残余应力使得物块的速度出现震荡现象，而且可以看出碰撞结束后有限元模拟得到的速度并没有达到理论上的20 m/s和0 m/s，一直处于震荡状态，相对而言，DSPH和DSFPM模拟结果比有限元模拟结果更稳定。

图 9 3种方法模拟得到的铝块速度随时间的变化曲线

Figure 9. Velocity-time cuvers of aluminum blocks by three methods

下载: 全尺寸图片幻灯片

图 10 有限元模拟得到的铝块速度随时间的变化曲线

Figure 10. Velocity-time cuvers of aluminum blocks by finite element method

下载: 全尺寸图片幻灯片

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

图 11 3种方法模拟得到的铝块碰撞过程中的速度分布

Figure 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方向应力分布

Figure 12. x-direction total stress distribution of the aluminum blocks simulated by three methods

下载: 全尺寸图片幻灯片

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

图 13 3种方法模拟得到点M处x方向应力随时间的变化曲线

Figure 13. x-directional total stress changes at point M simulated by three methods

下载: 全尺寸图片幻灯片

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/m³，弹性模量E₁=270 GPa，泊松比为0.3。右侧板长2.5 mm，高20 mm，密度为2 785 kg/m³，弹性模量E₂为72 GPa，泊松比为0.3。粒子间距d=0.2 mm，光滑长度h=1.2d。

图 14 物块冲击模型

Figure 14. The model for block collosion

下载: 全尺寸图片幻灯片

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

图 15 E₂=72 GPa时，不同时刻框图 7算法模拟的速度分布图

Figure 15. Velocity distribution of the blocks simulted by the algorithm in Fig. 7 at different times when E₂=72 GPa

下载: 全尺寸图片幻灯片

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

图 16 E₂=72 MPa时，不同时刻框图 7算法模拟的速度云图

Figure 16. Velocity distribution of the blocks simulted by the algorithm in Fig. 7 at different times when E₂=72 MPa

下载: 全尺寸图片幻灯片

4. 结论

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

图 1 界面附近粒子支持域

Figure 1. The supported domain of the particle near the interface

下载: 全尺寸图片幻灯片

图 2 一维情形下DSFPM粒子选取方式

Figure 2. Particle selection mode for 1-D DSFPM

下载: 全尺寸图片幻灯片

图 3 二维情形下DSFPM粒子选取方式

Figure 3. Particle selection mode for 2-D DSFPM

下载: 全尺寸图片幻灯片

图 4 3种方法估计非连续常值函数误差及其导数误差

Figure 4. Estimation errors for the discontinuous 0-order function and its derivative by three methods

下载: 全尺寸图片幻灯片

图 5 3种方法估计非连续一次函数误差及其导数误差

Figure 5. Estimation errors for the discontinuous 1-order function and its derivative by three methods

下载: 全尺寸图片幻灯片

图 6 DSFPM计算二次非连续函数本身和其偏y导数误差分布

Figure 6. Estimation error distribution for the discontinuous 2-order function and its partial y derivative by DSFPM

下载: 全尺寸图片幻灯片

图 7 处理碰撞冲击问题的算法框图

Figure 7. Algorithm diagram of the collision problem

下载: 全尺寸图片幻灯片

图 8 铝块碰撞模型

Figure 8. The model for the aluminum block collision

下载: 全尺寸图片幻灯片

图 9 3种方法模拟得到的铝块速度随时间的变化曲线

Figure 9. Velocity-time cuvers of aluminum blocks by three methods

下载: 全尺寸图片幻灯片

图 10 有限元模拟得到的铝块速度随时间的变化曲线

Figure 10. Velocity-time cuvers of aluminum blocks by finite element method

下载: 全尺寸图片幻灯片

图 11 3种方法模拟得到的铝块碰撞过程中的速度分布

Figure 11. Velocity distribution of the aluminum blocks by three methods

下载: 全尺寸图片幻灯片

图 12 3种方法模拟得到铝块碰撞过程中的x方向应力分布

Figure 12. x-direction total stress distribution of the aluminum blocks simulated by three methods

下载: 全尺寸图片幻灯片

图 13 3种方法模拟得到点M处x方向应力随时间的变化曲线

Figure 13. x-directional total stress changes at point M simulated by three methods

下载: 全尺寸图片幻灯片

图 14 物块冲击模型

Figure 14. The model for block collosion

下载: 全尺寸图片幻灯片

图 15 E₂=72 GPa时，不同时刻框图 7算法模拟的速度分布图

Figure 15. Velocity distribution of the blocks simulted by the algorithm in Fig. 7 at different times when E₂=72 GPa

下载: 全尺寸图片幻灯片

图 16 E₂=72 MPa时，不同时刻框图 7算法模拟的速度云图

Figure 16. Velocity distribution of the blocks simulted by the algorithm in Fig. 7 at different times when E₂=72 MPa

下载: 全尺寸图片幻灯片

参考文献(20)

[1]	MONAGHAN J J. Simulating free surface flows with SPH[J]. Journal of Computational Physics, 1994, 110:399-406. DOI: 10.1006/jcph.1994.1034.
[2]	杨秀峰, 刘谋斌.SPH方法Delaunay三角刨分与自由液面重构[J].计算力学学报, 2016, 33(4):594-598. DOI: 10.7511/jslx201604027. YANG Xiufeng, LIU Moubin. Delaunay triangulation and free surface extraction for SPH method[J]. Chinese Journal of Computational Mechanics, 2016, 33(4):594-598. DOI: 10.7511/jslx201604027.
[3]	龙厅, 胡德安, 韩旭.FE-ISPH与FE-WCSPH模拟流固耦合问题的比较研究[C]//中国计算力学大会.贵阳, 2014: 547. http://cpfd.cnki.com.cn/Article/CPFDTOTAL-AGLU201408004093.htm
[4]	刘谋斌, 宗智, 常建忠.光滑粒子动力学方法的发展与应用[J].力学进展, 2011, 41(2):219-236. DOI: 10.6052/1000-0992-2011-2-lxjzJ2010-078. LIU Moubin, ZONG Zhi, CHANG Jianzhong. Developments and applications of smoothed particle hydrodynamics[J]. Advances in Mechanics, 2011, 41(2):219-236. DOI: 10.6052/1000-0992-2011-2-lxjzJ2010-078.
[5]	傅学金, 强洪夫, 杨月诚.固体介质中SPH方法的拉伸不稳定性问题研究进展[J].力学进展, 2007, 37(3):375-388. DOI: 10.3321/j.issn:1000-0992.2007.03.005. FU Xuejin, QIANG Hongfu, YANG Yuecheng. Advances in the tensile instability of smoothed particle hydrodynamics applied to solid dynamics[J]. Advances in Mechanics, 2007, 37(3):375-388. DOI: 10.3321/j.issn:1000-0992.2007.03.005.
[6]	LIU W K, JUN S, LI S, et al. Reproducing kernel particle methods for structure dynamics[J]. International Journal for Numerical Methods in Engineering, 1995, 38(10):1655-1679. DOI: 10.1002/nme.1620381005.
[7]	CHEN J K, BERAUN J E. A generalized smoothed particle hydrodynamics method for nonlinear dynamic problem[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 190:225-239. DOI: 10.1016/S0045-7825(99)00422-3.
[8]	章杰, 苏少卿, 郑宇, 等.改进SPH方法在陶瓷材料层裂数值模拟中的应用[J].爆炸与冲击, 2013, 33(4):401-407. DOI: 10.3969/j.issn.1001-1455.2013.04.011. ZHANG Jie, SU Shaoqing, ZHENG Yu, et al. Application of modified SPH method to numerical simulation of ceramic spallation[J]. Explosion and Shock Waves, 2013, 33(4):401-407. DOI: 10.3969/j.issn.1001-1455.2013.04.011.
[9]	LIU M B, LIU G R. Restoring particle consistency in smoothed particle hydrodynamics[J]. Applied Numerical Mathematics, 2006, 56(1):19-36. DOI: 10.1016/j.apnum.2005.02.012.
[10]	郑兴, 段文洋.K2_SPH方法及其对二维非线性水波的模拟[J].计算物理, 2011, 28(5):659-666. DOI: 10.3969/j.issn.1001-246X.2011.05.004. ZHENG Xing, DUAN Wenyang. K2_SPH Method and application for 2D nonlinear water wave simulation[J]. Chinese Journal of Computational Physics, 2011, 28(5):659-666. DOI: 10.3969/j.issn.1001-246X.2011.05.004.
[11]	刘谋斌, 杨秀峰, 邵家儒.高精度SPH方法及其在海洋工程中的应用[C]//颗粒材料计算力学会议论文集.兰州, 2014: 39-41. http://cpfd.cnki.com.cn/Article/CPFDTOTAL-AGLU201408001007.htm
[12]	ADAMI S, HU X Y, ADAMS N A. A generalized wall boundary condition for smoothed particle hydrodynamics[J]. Journal of Computational Physics, 2012, 231(21):7057-7075. DOI: 10.1016/j.jcp.2012.05.005.
[13]	LIU M B, SHAO J R, CHANG J Z. On the treatment of solid boundary in smoothed particle hydrodynamics[J]. Science China:Technological Sciences, 2012, 55(1):244-254. DOI: 10.1007/s11431-011-4663-y.
[14]	LIU M B, LIU G R, LAM K Y. A one-dimensional meshfree particle formulation for simulating shock waves[J]. Shock Wave, 2003, 13:201-211. DOI: 10.1007/s00193-003-0207-0.
[15]	XU F, ZHAO Y, YAN R, et al. Multi-dimensional discontinuous SPH method and its application to metal penetration analysis[J]. International Journal for Numerical Methods in Engineering, 2013, 93:1125-1146. DOI: 10.1002/nme.4414.
[16]	闫蕊, 徐绯, 张岳青.DSPH方法的有效性验证及应用[J].爆炸与冲击, 2013, 33(2):133-139. DOI: 10.3969/j.issn.1001-1455.2013.02.004. YAN Rui, XU Fei, ZHANG Yueqing. Validation of DSPH method and its application to physical problems[J]. Explosion and Shock Waves, 2013, 33(2):133-139. DOI: 10.3969/j.issn.1001-1455.2013.02.004.
[17]	宋俊豪, 张超英, 梁朝湘, 等.RDSPH:一种适用于一维非连续条件的新SPH方法[J].广西师范大学学报(自然科学版), 2009, 27(3):9-13. DOI: 10.3969/j.issn.1001-6600.2009.03.003. SONG Junhao, ZHANG Chaoying, LIANG Chaoxiang, et al. A new one-dimensional smoothed particle hydrodynamics method in simulating discontinuous problem[J], Journal of Guangxi Normal University (Natural Science Edition), 2009, 27(3):9-13. DOI: 10.3969/j.issn.1001-6600.2009.03.003.
[18]	YANG Yang, XU Fei, ZHANG Meng, et al. An effective improved algorithm for finite particle method[J]. International Journal of Computational Methods, 2016, 13(4):1641009. DOI: 10.1142/S0219876216410097.
[19]	MONAGHAN J J. Smoothed particle hydrodynamic[J]. Annual Review of Astronomy and Astrophysics, 1992, 30(1):543-574. DOI: 10.1146/annurev.aa.30.090192.002551.
[20]	MONAGHAN J J, KAJTAR J. SPH particle boundary forces for arbitrary boundaries[J]. Computer Physics Communications, 2009, 180(10):1811-1820. DOI: 10.1016/j.cpc.2009.05.008.

施引文献

期刊类型引用(2)

1.	王敏，马小晶，陈尉，张博文. 基于FPM方法的入水问题数值模拟研究. 应用力学学报. 2023(01): 199-210 . 百度学术
2.	王璐，徐绯，杨扬. 完全拉格朗日SPH在冲击问题中的改进和应用. 力学学报. 2022(12): 3297-3309 . 百度学术

其他类型引用(0)

资源附件(0)

访问统计

图(16)

计量

文章访问数: 6290
HTML全文浏览量: 2192
PDF下载量: 56
被引次数: 2

1. DSFPM
1.1 FPM基本理论
1.2 DSPH基本理论
1.3 DSFPM推导
2. DSFPM粒子近似精度分析
2.1 C0阶非连续性分析
2.2 C1阶非连续性分析
2.3 C2阶非连续性分析
3. 碰撞算例模拟
3.1 离散公式
3.2 小变形碰撞
3.3 大变形碰撞
4. 结论

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

doi: 10.11883/bzycj-2017-0390

作者简介: 王璐(1995-), 女, 博士研究生, wanglulu@mail.nwpu.edu.cn

通讯作者: 徐绯(1970-), 女, 博士, 教授, 博士生导师, xufei@nwpu.edu.cn

计量

出版历程

An improved finite particle method for discontinuous interface problems

1. DSFPM

1.1 FPM基本理论

1.2 DSPH基本理论

1.3 DSFPM推导

2. DSFPM粒子近似精度分析

2.1 C0阶非连续性分析

2.2 C1阶非连续性分析

2.3 C2阶非连续性分析

3. 碰撞算例模拟

3.1 离散公式

3.2 小变形碰撞

3.2.1 速度分布

3.2.2 应力分布

3.2.3 计算效率

3.3 大变形碰撞

4. 结论

期刊类型引用(2)

其他类型引用(0)

计量

出版历程

目录

1. DSFPM

1.1 FPM基本理论

1.2 DSPH基本理论

1.3 DSFPM推导

2. DSFPM粒子近似精度分析

2.1 C0阶非连续性分析

2.2 C1阶非连续性分析

2.3 C2阶非连续性分析

3. 碰撞算例模拟

3.1 离散公式

3.2 小变形碰撞

3.3 大变形碰撞

4. 结论

作者简介:
王璐(1995-), 女, 博士研究生, wanglulu@mail.nwpu.edu.cn

通讯作者:
徐绯(1970-), 女, 博士, 教授, 博士生导师, xufei@nwpu.edu.cn

2.1 C⁰阶非连续性分析

2.2 C¹阶非连续性分析

2.3 C²阶非连续性分析

2.1 C⁰阶非连续性分析

2.2 C¹阶非连续性分析

2.3 C²阶非连续性分析