光滑粒子流体动力学方法SPH资料讲解
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
后再检测粒子i的搜索立方体空间是否与并列层次内的其他节点所占的空间有重 合的地方。若没有，则中止往下搜索；若有，则继续往下一层搜索。直到所搜 索到的当前节点处只有一个粒子为止。接着，检查此粒子是否在给定粒子的支 持域内，若是，则将其记为粒子的相邻粒子。
2020/6/15
人工粘度
迄今为止，与SPH相关的论著中，Monaghan型的人工粘度是最
2020/6/15
型号II的虚粒子没有固定的参数。它们是在每个计算步中由对应的 实粒子对称产生的。可以应用型号II的虚粒子来处理固定边界和自由表 面。
通过数值算例的测试验证了应用虚粒子来处理边界的可靠性和有效 性。其不仅仅提高了SPH近似法在边界区域处的精度，而且防止了粒子 非物理穿透边界。
2020/6/15
全配对粒子搜索法 对于给定的粒子i，应用全配对搜索法即是计算粒子i到粒子j的距离。
若该距离小于粒子i的支持域的半径 h，则粒子j为粒子i的支持域内
的粒子。
2020/6/15
2h 在二维空间里应用链表 搜索法搜索最近相邻粒子
2020/6/15
在光滑长度为空间常量的情况下，应用链表 搜索法非常有效。
光滑函数
光滑函数的性质： 一、正则化条件
Wxx,hdx1
由于光滑函数的积分值等于1，故此条件也称为归一化条件。 二、当光滑长度趋向于零时具有狄拉克函数性质
liW m x x ,h x x
h 0
三、紧支性条件
W xx,h0 xx h
2020/6/15
光滑函数
Monaghan和Lattanzio在三次样条函数的基础上提出了
2020/6/15
边界处理
当类型I的虚粒子成为邻近边界处的实粒子的相邻粒子时，则会 在沿着两粒子的中心线处对实粒子产生一个作用力。
PBij
D
r0 rij
n1
r0 rij
n2
xij ri2j
,
0,
r0 rij
1
r0 rij
1
式中：参数n1和n2一般取值分别为12和4。D是由具体问题而定的参数， 一般取与速度最大值的平方相等的量级。截止半径r0在此问题的模拟 分析中非常重要，在一般情况下， r0的取值与粒子的初始间距的大小 相近。
弱形式就是从加权余量法或变分原理出发，把微分方程及其定解 条件转换成弱形式(Weak Form)或Galerkin形式，即用测试函数 (Test function)与控制方程相乘后在全局或部分区域内积分，并利 用高斯散度定理得到不同形式的弱形式，然后进行离散化求解。
通过引进新的无网格近似函数构造方法或采用新的偏微分方程的 离散形式，就可以期待开发出更加高效和精确的无网格方法。
Lucy, Gingold(1977)分别提出了SPH方法，最早用于天体物理 现象的模拟，随后别广泛地应用于连续固体力学和流体力学中。
2020/6/15
2020/6/15
2020/6/15
2020/6/15
11kg弹丸1418m/s撞靶速度下穿靶过程的数值模拟
2020/6/151500m/s速度下弹体侵彻混凝土靶变形过程的数值模拟
2020/6/15
五次样条函数及其一阶导数
光滑函数一览表
2020/6/15
最近相邻粒子搜索法(NNPS)
一般将包含在支持域中的粒子称为相关粒子的最近相邻粒子(NNP)。 寻找最近相邻粒子的过程通常称为最近相邻粒子搜索(NNPS)。 在SPH方法中常用的三种NNPS方法为： 全配对搜索法(all-pair search) 链表搜索法(linked-list search algorithm) 树形搜索法(tree search algorithm)
2020/6/15
边界处理
在Liu等的研究中，使用了虚粒子来处理固定边界条件，提出了 两种类型的虚粒子。第一种类型的粒子（型号I）设置在固定边界上； 第二种类型的粒子（型号II ）分布在边界的邻域内。
型号II 的虚粒子具有与相应实粒子相同的密度和压力，但速度相反。 型号I的虚粒子被引入到实粒子的核函数核粒子近似法中。
2020/6/15
边界处理
由于在边界上或邻近边界处的粒子存在缺陷，即在积分的时候会 被边界截断，故而SPH方法不能完全适用于整个区域。在边界上或邻 近边界处的粒子只受到边界内的粒子的影响作用，而边界外由于没有 粒子，故而边界外不对粒子产生影响。这种单边影响作用会导致求解 结果错误，因为在固定边界表面，虽然粒子速度为零，但是其他变量， 如密度，则不一定为零。
在实现链表算法时，要在问题域上铺设一临 时网格，网格单元的空间大小应选取与支持 域的空间大小一致。那么，对于给定的粒子 i，其相邻粒子只能在同一网格单元内，或 者在紧密相邻的单元内。所以，当时，一维、 二维和三维空间里的搜索范围分别为3，9， 27个单元内。
链表搜索法存在的问题是，当光滑长度可变 时，网格空间就不一定适应每一个粒子，此 时若再应用链表搜索法，则搜索效率就会很 低。
p xvii pii jN 1mjj vijW xiij p x vii 1i jN 1mj pjj vij W xiij
ddiet1 2jN 1mj pi ipj j vijW xiij
2020/6/15
p x vii 1 2jN 1mj pi ip j j vij W xiij
粒子近似法
在粒子i处的函数的粒子近似式最终可写为：
f xi
N mj f xj j1 j
Wij
f xi
N
mj f
xj
iWij
j1 j
W ij W x i xj,h
上式说明了粒子i处的任一函数值可通过应用光滑函数对 其紧支域内所有粒子相对应的函数值进行加权平均近似。
2020/6/15
目录
➢SPH计算公式 ➢光滑函数
➢最近相邻粒子搜索法(NNPS) ➢人工粘度 ➢边界处理
➢交界面处理 ➢光滑长度的更新
2020/6/15
➢SPH方程的求解 ➢SPH程序结构 ➢激波管问题
SPH计算公式
1、密度的粒子近似法
由于粒子的分配与光滑长度的变化主要依赖于密度，故在SPH法 中密度近似法非常重要。
为广泛使用的人工粘度。它不仅将动能转换为热能，提供了冲击波阵
面必不可少的耗散，而且防止了粒子相互接近时的非物理穿透。具体
表达式如下：
ij cijiijj i2j ,
0,
vij xij 0 vij xij 0
ij
hij vij xij
xij 2 2
vij vi vj
cij
1 2
ci
cj
xij xi xj
W被称为光滑核函数(smoothing kernel function)或光滑 函数(smoothing function)，简称为核(kernel)函数。 h是定义光滑函数W的影响区域的光滑长度。
2020/6/15
粒子近似法
与SPH核近似法相关的连续积分表示式，可转化为支 持域内所有粒子叠加求和的离散化形式。
2020/6/15
此对称方程的优点为：可降低 粒子不一致问题产生的误差。
3、能量方程的粒子近似法
能量守恒方程：
de dt
p
v x
f xi
N
mj f
xj
iWij
j1 j
p x v p2 x vp2ddt
di
dt
i
N j1
mj
j
vij
W xiij
p x v 1 x pvv xp
在SPH法中有两种方法对密度进行展开，第一种方法是对密度直 接用SPH近似法，称为密度求和法。第二种方法是连续性密度法，通 过应用SPH近似法的概念对连续性方程进行转换而得到。
密度求和法：
f xi
N mj f
j1 j
xj
Wij
N
i mjWij j1
2020/6/15
改进方案（正则化）
2020/6/15
f x f x W x x , h d x
N
f x j W x x j , h V j
j1
N
f
j1
xjW
x x j,h
1 j
j V j
N
f
j1
xjW
x x j,h
1 j
mj
N mj j 1 j
f
xjW
x x j,h
N
m jW ij
i
j1
N j1
m
j j
W
ij
此方法可提高自由边界处和相同材 料粒子密度不连续交界面处的精度
连续性密度法：
d
dt
v x
dditi
N j1
m jj vj
W xiij
f xi
N
mj f
xj
iWij
j1 j
1
N mj
j1 j
Wij xi
0
i jN 1m jj vi W xiijivi jN 1m jj W xiij
强形式 以各种全局或局部加权余量法为统一框架的弱形式
2020/6/15
强形式是直接从微分方程及其定解条件出发，将近似函数及其导 数的估计形式带入基本方程、本构方程和初边值条件中去，联立 方程进行求解。该方法思路简单，便于程序编制，应用范围广泛， 在流体和固体的计算中都有所发展，适用于计算激波、高速冲击、 爆轰、穿甲等冲击动力学问题。但此类算法的精度较低，稳定性 较差，且边界条件的引入比较困难。
ij
1 2
i
j
hij
1 2
hi
hj
2020/6/15
在以上所有方程中 和 为标准常数，一般取值在1.0左右。因子 0用.1hij 于防止粒子相互靠近时产生的数值发散。式中与 相关的项得到的是体
积粘度，而与 相关的项是用于防止在高马赫数时粒子的相互穿透。式
中所给出的人工粘度被引入到SPH方程的压力项中。
2020/6/15
SPH方程的构造常按两个关键步进行。 第一步为积分表示法，又称场函数近似法； 第二步为粒子近似法。
光滑粒子流体动力学——一种无网格粒子法， 湖南大学出版社，G. R. Liu，M. B. Liu[著]
2020/6/15
场函数核近似法（积分表示法）
函数核近似的标准表达式：
fxfxW xx,hdx
N
mj f
xj
iWij
j1 j
将动量方程等号右端的梯度项直接应用SPH粒子近似法进行变换得：
dvi
dt
1
i
N
mj j1
j
Wij
j xi
此外，有：
N
m j
j1
iij W xiij
ii jN 1m jj W xiij 0
将以上两式相加可得：
div
dt
jN 1mj iij j W xiij
称为B样条函数的光滑函数 ：
WR,h ad 01632,2R2R312, R3,
0 R1
1 R2 R2
x x R
h
在一维、二维和三维空间中分别有：a d
1 h
， 15 7h 2
和
3。 2h 3
现有SPH文献中最为广泛应用的光滑函数
2020/6/15
三次样条函数及其一阶导数
四次样条函数及其一阶导数
SPH光滑粒子流体动力学方法 (smoothed particle hydrodynamics)
2020/6/15
报告人：马天宝
2013.Baidu Nhomakorabea.25
无网格方法的主要思想：通过使用一系列任意分布的节点（或
粒子）来求解各种各样边界条件的积分方程或偏微分方程组(PDEs) 从而得到精确稳定的数值解，这些节点或粒子之间不需要网格进行 连接。
2020/6/15
侵彻过程弹体温度分布云图
2020/6/15
碎浪与弹性挡墙之间的相互作用
无网格法(Meshfree Methods)
近似函数
两 构造方法 条
主 线
偏微分方程 的离散形式
核估计方法(Kernel Approximation, KA) 移动最小二乘法(Moving Least Square, MLS) 再生核估计方法(Repuducing Kernel Method, RKM) 径向基函数方法(Radial Basic Function, RBF) 单位分解方法(Patition of Unity, PU)
di
dt
i
N mj
j1 j
vij
W xiij
vij vi vj
2020/6/15
对于广义流体问题的模拟，应用修正的密度求和法可 得到较好的结果，对于具有强间断问题的模拟（如爆炸、 高速冲击等），应优先选取连续性密度法。
2020/6/15
2、动量方程的粒子近似法
动量守恒方程：
dv dt
1
x
f xi
树形搜索法非常适宜求解具有可变光滑长度的问题。这种算法是通过粒子 的位置来构造有序树。一旦树形结构构造起来后，便能高效地搜索最近相邻粒 子。
树形搜索法将最大问题域递归分割成一个个卦限，直到每一个卦限内只包 含一个粒子为止。树形结构构造完成后，即可以开始进行最近相邻粒子搜索。
给定任一粒子i, 并以i为中心，用边长为 2hi的立方体将粒子包围起来，然
2020/6/15
人工粘度
迄今为止，与SPH相关的论著中，Monaghan型的人工粘度是最
2020/6/15
型号II的虚粒子没有固定的参数。它们是在每个计算步中由对应的 实粒子对称产生的。可以应用型号II的虚粒子来处理固定边界和自由表 面。
通过数值算例的测试验证了应用虚粒子来处理边界的可靠性和有效 性。其不仅仅提高了SPH近似法在边界区域处的精度，而且防止了粒子 非物理穿透边界。
2020/6/15
全配对粒子搜索法 对于给定的粒子i，应用全配对搜索法即是计算粒子i到粒子j的距离。
若该距离小于粒子i的支持域的半径 h，则粒子j为粒子i的支持域内
的粒子。
2020/6/15
2h 在二维空间里应用链表 搜索法搜索最近相邻粒子
2020/6/15
在光滑长度为空间常量的情况下，应用链表 搜索法非常有效。
光滑函数
光滑函数的性质： 一、正则化条件
Wxx,hdx1
由于光滑函数的积分值等于1，故此条件也称为归一化条件。 二、当光滑长度趋向于零时具有狄拉克函数性质
liW m x x ,h x x
h 0
三、紧支性条件
W xx,h0 xx h
2020/6/15
光滑函数
Monaghan和Lattanzio在三次样条函数的基础上提出了
2020/6/15
边界处理
当类型I的虚粒子成为邻近边界处的实粒子的相邻粒子时，则会 在沿着两粒子的中心线处对实粒子产生一个作用力。
PBij
D
r0 rij
n1
r0 rij
n2
xij ri2j
,
0,
r0 rij
1
r0 rij
1
式中：参数n1和n2一般取值分别为12和4。D是由具体问题而定的参数， 一般取与速度最大值的平方相等的量级。截止半径r0在此问题的模拟 分析中非常重要，在一般情况下， r0的取值与粒子的初始间距的大小 相近。
弱形式就是从加权余量法或变分原理出发，把微分方程及其定解 条件转换成弱形式(Weak Form)或Galerkin形式，即用测试函数 (Test function)与控制方程相乘后在全局或部分区域内积分，并利 用高斯散度定理得到不同形式的弱形式，然后进行离散化求解。
通过引进新的无网格近似函数构造方法或采用新的偏微分方程的 离散形式，就可以期待开发出更加高效和精确的无网格方法。
Lucy, Gingold(1977)分别提出了SPH方法，最早用于天体物理 现象的模拟，随后别广泛地应用于连续固体力学和流体力学中。
2020/6/15
2020/6/15
2020/6/15
2020/6/15
11kg弹丸1418m/s撞靶速度下穿靶过程的数值模拟
2020/6/151500m/s速度下弹体侵彻混凝土靶变形过程的数值模拟
2020/6/15
五次样条函数及其一阶导数
光滑函数一览表
2020/6/15
最近相邻粒子搜索法(NNPS)
一般将包含在支持域中的粒子称为相关粒子的最近相邻粒子(NNP)。 寻找最近相邻粒子的过程通常称为最近相邻粒子搜索(NNPS)。 在SPH方法中常用的三种NNPS方法为： 全配对搜索法(all-pair search) 链表搜索法(linked-list search algorithm) 树形搜索法(tree search algorithm)
2020/6/15
边界处理
在Liu等的研究中，使用了虚粒子来处理固定边界条件，提出了 两种类型的虚粒子。第一种类型的粒子（型号I）设置在固定边界上； 第二种类型的粒子（型号II ）分布在边界的邻域内。
型号II 的虚粒子具有与相应实粒子相同的密度和压力，但速度相反。 型号I的虚粒子被引入到实粒子的核函数核粒子近似法中。
2020/6/15
边界处理
由于在边界上或邻近边界处的粒子存在缺陷，即在积分的时候会 被边界截断，故而SPH方法不能完全适用于整个区域。在边界上或邻 近边界处的粒子只受到边界内的粒子的影响作用，而边界外由于没有 粒子，故而边界外不对粒子产生影响。这种单边影响作用会导致求解 结果错误，因为在固定边界表面，虽然粒子速度为零，但是其他变量， 如密度，则不一定为零。
在实现链表算法时，要在问题域上铺设一临 时网格，网格单元的空间大小应选取与支持 域的空间大小一致。那么，对于给定的粒子 i，其相邻粒子只能在同一网格单元内，或 者在紧密相邻的单元内。所以，当时，一维、 二维和三维空间里的搜索范围分别为3，9， 27个单元内。
链表搜索法存在的问题是，当光滑长度可变 时，网格空间就不一定适应每一个粒子，此 时若再应用链表搜索法，则搜索效率就会很 低。
p xvii pii jN 1mjj vijW xiij p x vii 1i jN 1mj pjj vij W xiij
ddiet1 2jN 1mj pi ipj j vijW xiij
2020/6/15
p x vii 1 2jN 1mj pi ip j j vij W xiij
粒子近似法
在粒子i处的函数的粒子近似式最终可写为：
f xi
N mj f xj j1 j
Wij
f xi
N
mj f
xj
iWij
j1 j
W ij W x i xj,h
上式说明了粒子i处的任一函数值可通过应用光滑函数对 其紧支域内所有粒子相对应的函数值进行加权平均近似。
2020/6/15
目录
➢SPH计算公式 ➢光滑函数
➢最近相邻粒子搜索法(NNPS) ➢人工粘度 ➢边界处理
➢交界面处理 ➢光滑长度的更新
2020/6/15
➢SPH方程的求解 ➢SPH程序结构 ➢激波管问题
SPH计算公式
1、密度的粒子近似法
由于粒子的分配与光滑长度的变化主要依赖于密度，故在SPH法 中密度近似法非常重要。
为广泛使用的人工粘度。它不仅将动能转换为热能，提供了冲击波阵
面必不可少的耗散，而且防止了粒子相互接近时的非物理穿透。具体
表达式如下：
ij cijiijj i2j ,
0,
vij xij 0 vij xij 0
ij
hij vij xij
xij 2 2
vij vi vj
cij
1 2
ci
cj
xij xi xj
W被称为光滑核函数(smoothing kernel function)或光滑 函数(smoothing function)，简称为核(kernel)函数。 h是定义光滑函数W的影响区域的光滑长度。
2020/6/15
粒子近似法
与SPH核近似法相关的连续积分表示式，可转化为支 持域内所有粒子叠加求和的离散化形式。
2020/6/15
此对称方程的优点为：可降低 粒子不一致问题产生的误差。
3、能量方程的粒子近似法
能量守恒方程：
de dt
p
v x
f xi
N
mj f
xj
iWij
j1 j
p x v p2 x vp2ddt
di
dt
i
N j1
mj
j
vij
W xiij
p x v 1 x pvv xp
在SPH法中有两种方法对密度进行展开，第一种方法是对密度直 接用SPH近似法，称为密度求和法。第二种方法是连续性密度法，通 过应用SPH近似法的概念对连续性方程进行转换而得到。
密度求和法：
f xi
N mj f
j1 j
xj
Wij
N
i mjWij j1
2020/6/15
改进方案（正则化）
2020/6/15
f x f x W x x , h d x
N
f x j W x x j , h V j
j1
N
f
j1
xjW
x x j,h
1 j
j V j
N
f
j1
xjW
x x j,h
1 j
mj
N mj j 1 j
f
xjW
x x j,h
N
m jW ij
i
j1
N j1
m
j j
W
ij
此方法可提高自由边界处和相同材 料粒子密度不连续交界面处的精度
连续性密度法：
d
dt
v x
dditi
N j1
m jj vj
W xiij
f xi
N
mj f
xj
iWij
j1 j
1
N mj
j1 j
Wij xi
0
i jN 1m jj vi W xiijivi jN 1m jj W xiij
强形式 以各种全局或局部加权余量法为统一框架的弱形式
2020/6/15
强形式是直接从微分方程及其定解条件出发，将近似函数及其导 数的估计形式带入基本方程、本构方程和初边值条件中去，联立 方程进行求解。该方法思路简单，便于程序编制，应用范围广泛， 在流体和固体的计算中都有所发展，适用于计算激波、高速冲击、 爆轰、穿甲等冲击动力学问题。但此类算法的精度较低，稳定性 较差，且边界条件的引入比较困难。
ij
1 2
i
j
hij
1 2
hi
hj
2020/6/15
在以上所有方程中 和 为标准常数，一般取值在1.0左右。因子 0用.1hij 于防止粒子相互靠近时产生的数值发散。式中与 相关的项得到的是体
积粘度，而与 相关的项是用于防止在高马赫数时粒子的相互穿透。式
中所给出的人工粘度被引入到SPH方程的压力项中。
2020/6/15
SPH方程的构造常按两个关键步进行。 第一步为积分表示法，又称场函数近似法； 第二步为粒子近似法。
光滑粒子流体动力学——一种无网格粒子法， 湖南大学出版社，G. R. Liu，M. B. Liu[著]
2020/6/15
场函数核近似法（积分表示法）
函数核近似的标准表达式：
fxfxW xx,hdx
N
mj f
xj
iWij
j1 j
将动量方程等号右端的梯度项直接应用SPH粒子近似法进行变换得：
dvi
dt
1
i
N
mj j1
j
Wij
j xi
此外，有：
N
m j
j1
iij W xiij
ii jN 1m jj W xiij 0
将以上两式相加可得：
div
dt
jN 1mj iij j W xiij
称为B样条函数的光滑函数 ：
WR,h ad 01632,2R2R312, R3,
0 R1
1 R2 R2
x x R
h
在一维、二维和三维空间中分别有：a d
1 h
， 15 7h 2
和
3。 2h 3
现有SPH文献中最为广泛应用的光滑函数
2020/6/15
三次样条函数及其一阶导数
四次样条函数及其一阶导数
SPH光滑粒子流体动力学方法 (smoothed particle hydrodynamics)
2020/6/15
报告人：马天宝
2013.Baidu Nhomakorabea.25
无网格方法的主要思想：通过使用一系列任意分布的节点（或
粒子）来求解各种各样边界条件的积分方程或偏微分方程组(PDEs) 从而得到精确稳定的数值解，这些节点或粒子之间不需要网格进行 连接。
2020/6/15
侵彻过程弹体温度分布云图
2020/6/15
碎浪与弹性挡墙之间的相互作用
无网格法(Meshfree Methods)
近似函数
两 构造方法 条
主 线
偏微分方程 的离散形式
核估计方法(Kernel Approximation, KA) 移动最小二乘法(Moving Least Square, MLS) 再生核估计方法(Repuducing Kernel Method, RKM) 径向基函数方法(Radial Basic Function, RBF) 单位分解方法(Patition of Unity, PU)
di
dt
i
N mj
j1 j
vij
W xiij
vij vi vj
2020/6/15
对于广义流体问题的模拟，应用修正的密度求和法可 得到较好的结果，对于具有强间断问题的模拟（如爆炸、 高速冲击等），应优先选取连续性密度法。
2020/6/15
2、动量方程的粒子近似法
动量守恒方程：
dv dt
1
x
f xi
树形搜索法非常适宜求解具有可变光滑长度的问题。这种算法是通过粒子 的位置来构造有序树。一旦树形结构构造起来后，便能高效地搜索最近相邻粒 子。
树形搜索法将最大问题域递归分割成一个个卦限，直到每一个卦限内只包 含一个粒子为止。树形结构构造完成后，即可以开始进行最近相邻粒子搜索。
给定任一粒子i, 并以i为中心，用边长为 2hi的立方体将粒子包围起来，然