《无网格弱可压缩SPH数值算法及应用扩展》介绍了一种无网格数值模拟方法，即“弱可压缩光滑粒子动力学方法”。该方法是光滑粒子动力学（smoothed particle hydrodynamics，SPH）方法的一个分支，其原理简单、数值算法简洁，在处理自由表面流动、界面流动、流固耦合等问题时具有独特的优势。《无网格弱可压缩SPH数值算法及应用扩展》基于作者在弱可压缩SPH数值算法和应用方面的研究基础，对该方法的数学原理、数值算法和计算机编程方法进行了归纳和整理。《无网格弱可压缩SPH数值算法及应用扩展》针对具体的科学和工程问题，包括多相界面流动、表面张力、液滴冲击、毛细现象、颗粒冲蚀、磨料射流等，运用弱可压缩SPH方法建立相应的数值计算模型，形成了一些适用于特定问题的新型算法。《无网格弱可压缩SPH数值算法及应用扩展》在附录中提供了《无网格弱可压缩SPH数值算法及应用扩展》部分算例的计算机代码，它们采用Fortran 90语言编写。

第1章 绪论
　　1.1 数值模拟概述
　　计算机数值模拟已经成为解决复杂工程和科学问题的重要手段，它一般包含以下四个步骤：
　　（1）建立数学模型：将所关心的物理问题采用数学方程进行描述，例如，描述黏性流体运动的纳维-斯托克斯（Navier-Stokes）方程；
　　（2）离散化：将计算域分解为离散的计算单元，利用数值近似方法，将微分形式的数学方程转化为可求解的代数方程；
　　（3）编程求解：建立数值求解算法及流程，按照算法流程编制计算机程序，通过计算机运算完成离散化代数方程的求解；
　　（4）结果处理及分析：根据分析的需求提取求解后的计算数据，对数据进行可视化处理，如绘制流场的压力云图、速度云图等。
　　图1.1给出了数值模拟的一般流程。
　　图1.1 数值模拟的一般流程
　　数学模型是基于对物理现象的观察，并引入一些简化或假设条件之后建立的。在连续介质力学领域，数学模型通常表达为偏微分形式的控制方程（governing equations），并且需要指定边界条件和初始条件。一般来讲，控制方程、边界条件和初始条件共同决定了计算问题的解。为了对控制方程进行求解，需要将计算域离散化，即将计算域拆解为有限数量的离散单元。一种常见的离散模式是将计算域分解为网格和节点，网格和节点之间通过特定的拓扑关系进行约束，工程中常用的有限单元法（finite element method，FEM）和有限体积法（finite volume method，FVM）就是采用这种模式，它们被归类为“网格法”。
　　数学模型的求解过程需要采用计算机运算来实现，采用计算机语言（如：Fortran，C++等），按照算法流程编制计算机程序，完成对数学模型的计算机求解。计算精度、效率、鲁棒性以及程序的可读性是计算机编程时要考虑的主要指标。计算完成后，可以将计算结果与理论解和实验数据进行对比，对计算结果的准确性进行验证，该过程称为模型验证（model validation）。
　　目前，常用的数值模拟方法，大多是在连续介质力学（continuum mechanics）理论框架内提出的，基于连续介质力学的数值方法是计算力学的一个重要分支。按照研究对象分类，可以将数值模拟方法分为计算固体力学（computational solid mechanics，CSM）和计算流体力学（computational fluid dynamics，CFD）。以CFD为例，流体的控制方程表达为质量、动量和能量守恒方程，即流体介质的质量、动量和能量在流动过程中保持守恒。在守恒方程的基础上，结合材料本构关系、边界条件及初始条件，完成对特定流体问题的求解。
　　1.2 网格方法及其局限性
　　按照计算域离散点、线、面的拓扑关系，数值模拟方法可以分为网格法和无网格法。传统的网格法包括有限差分法（finite difference method，FDM）、有限单元法（FEM）、有限体积法（FVM）等。根据控制方程的描述方法，数值模拟方法又可以分为欧拉方法和拉格朗日方法。在这两类方法中，网格节点和场变量信息的关联方式不同：欧拉方法采用空间描述方式，计算节点一般是固定的，不随时间变化；拉格朗日方法是一种随体描述方法，材料属性是固着在计算节点上的，节点的位置会随着材料的变形或流体的流动而改变。
　　在网格法中，欧拉模式和拉格朗日模式分别对应了两种不同的网格类型，也称为欧拉网格和拉格朗日网格。欧拉网格是固定的，即计算过程中网格和节点是固定不动的；拉格朗日网格是随体网格，网格和节点与材料是固结在一起的。欧拉网格通常用于计算流体力学分析，拉格朗日网格通常用于计算固体力学分析。FEM是应用*广泛的拉格朗日方法之一，它在计算固体受力产生的变形时，相邻节点的移动会引起网格单元的膨胀或压缩，由于单元内材料质量是固定不变的，质量不会穿过网格边界，因此网格变形的程度与材料变形是一致的。由以上可知，拉格朗日方法具有以下特点：
　　（1）网格依附在材料上，便于在时间历程上追踪材料上任意点的场变量变化情况；
　　（2）网格节点可以位于边界或界面上，对自由表面、移动边界问题具有较好的适应性；
　　（3）通过配置不规则或适应性好的四面体网格，可以适应具有复杂几何形状的计算域。
　　拉格朗日方法的以上特点，使得拉格朗日方法通常用于计算固体力学分析。但是，拉格朗日网格方法不适用于材料大变形问题，材料大变形会导致严重的网格变形或扭曲，从而影响计算的稳定性和精度。
　　与拉格朗日方法不同，欧拉方法采用固定的网格系统，因此也就不会出现拉格朗日网格所具有的网格变形问题，这也是欧拉法被广泛应用于流体流动模拟的原因之一。但是，欧拉方法也有以下缺点：
　　（1）材料与网格不是固结在一起的，存在相对运动，不利于获得场变量信息随时间的变化规律；
　　（2）欧拉方法可以追踪穿过网格边界的质量、动量和能量通量，但很难精确地确定自由表面、可动边界和多相界面的位置信息；
　　（3）欧拉网格应该完全覆盖流体的计算域，而且需要足够大的计算域来覆盖流体可能流过的区域，如果为了降低计算量而使用较粗糙的网格，将不利于对局部流动细节的捕捉。
　　传统的网格方法，例如FDM，FEM和FVM等，已经广泛应用于工业领域和自然科学中，解决各类流体、固体以及流固耦合问题。尽管网格法已经取得了巨大成功，它仍旧存在一些固有的缺陷，限制了它在某些领域或问题中的应用。众所周知，采用网格法进行数值模拟的第一步是网格划分。以FVM为例，针对不规则或形状复杂的计算域划分网格从来都不是轻松的工作，通常需要复杂的数学变换来处理网格适应问题。因此，为了应对这些问题，学者开始提出无网格数值模拟理念和方法，解决一些传统网格方法难于或无法解决的问题，对于完善数值模拟理论、扩展计算机数值模拟应用具有重要的意义。
　　1.3 弱可压缩SPH方法
　　光滑粒子动力学（smoothed particle hydrodynamics，SPH）方法是应用*广泛的无网格方法之一，它采用粒子来离散计算域，对具有复杂几何形状的计算域以及材料大变形问题有很好的适用性。SPH方法*早由Lucy[1]，Monaghan[2]等学者提出并应用在天体物理学中，后来经过了Benz[3]等的进一步完善。在20世纪90年代，SPH方法得到了迅速发展。Monaghan[4]等将SPH方法应用于自由表面流动；Morris等[5]将SPH方法应用于低雷诺数黏性流动。目前，SPH在流体力学分析方面的应用非常广泛，涉及多相流[6]、不可压缩流[7]、多孔介质流动[8]、自由表面流动[9]等。在固体力学领域，SPH方法的应用包括高速冲击[10]、结构断裂[11]、刀具切割等。
　　在SPH方法发展历程中，衍生出了不可压缩SPH方法和弱可压缩SPH方法。弱可压缩SPH方法采用状态方程来计算压力。在弱可压缩SPH方法中，流体密度是可以发生微弱变化的，但是，密度的改变量被限制在一个较小的范围内，此时流体表现出的流动行为与不可压缩流体是相似的，从而实现以“弱可压缩”假设来处理“不可压缩”流动问题。由于弱可压缩SPH方法不需要求解泊松方程，因此它的数值算法更简洁。
　　此外，弱可压缩SPH方法有以下显著的优势[12]：
　　（1）SPH是一种拉格朗日方法，通过追踪粒子的运动，可以轻易获得流动系统随时间的变化规律，在自由表面流、多相流和可动边界问题的数值模拟方面具有优势；
　　（2）粒子的运动不受点、线、面拓扑关系的约束，避免了网格变形、扭曲带来的数值不稳定性问题，在处理固体大变形以及对流占优流体运动方面，具有独特优势；
　　（3）SPH采用“粒子”离散计算域，对具有复杂几何形状的计算域有更好的适应性；
　　（4）SPH粒子代表一定质量的介质，如果在计算中保持粒子的质量不变，则质量守恒定律自动满足；
　　（5）由于SPH方法的拉格朗日性质，在自由表面或多相流模拟时，不需进行界面追踪或重构，对相界面变化剧烈的情况有较好的适应性。
　　但是，弱可压缩SPH方法也有以下不足：
　　（1）SPH方法存在张力不稳定（tensile instability）问题，在求解流体负压或者固体拉伸问题时，会出现粒子非物理集聚或断裂现象。在SPH模拟中，可以通过施加足够大的背景压力来避免形成负压区域，或在动量方程中加入阻碍粒子集聚的人工应力项，或使用高阶核函数和核梯度校正可以在一定程度上缓解张力不稳定性问题，但是会影响SPH算法的精度和计算效率。
　　（2）SPH数值计算涉及大量的参数调整和选择，这些参数对计算精度和稳定性影响较大，如何量化和评估这些参数对计算结果的影响，是在使用SPH方法建模和模拟时需要重点解决的问题。
　　（3）SPH方法采用核积分和粒子积分完成数值近似过程，其数值近似精度受粒子分布和边界的影响很大，引入修正算法可以提升数值近似精度，但会增加算法的实施难度和计算效率。
　　（4）SPH模型一般采用显式时间积分算法求解，时间步长受Courant-Friedrichs-Lewy（CFL）条件约束。在计算过程中需要实时地进行邻域粒子搜索（nearest neighboring particle searching，NNPS），计算效率一般低于网格类方法。
　　1.4 SPH方法的数值模拟应用
　　1.4.1 自由表面流动
　　20世纪70年代，SPH*初被提出，应用于天体物理学领域[1，2]。由于星系运动通常没有传统意义上的“边界”，“星球流体”的“密度”会随时间和空间产生跨数量级的变化。正是如此，学者发现SPH方法在处理密度变化和无边界流动问题方面具有优势。这种优势源于SPH方法的拉格朗日特性：粒子代表介质，介质随粒子运动，不受流动边界的约束。在20世纪80年代，SPH除了继续应用在天体物理学领域外，也开始逐步应用于连续介质力学领域，如固体力学和结构断裂问题。直到20世纪90年代中期，Monaghan首次尝试将SPH用于“普通”流体，如自由表面流动（free surface flow）问题[4]。Monaghan采用压力状态方程来模拟不可压缩流动，流体的不可压缩性通过使用足够大的人工声速来实现。Monaghan的尝试表明，SPH在模拟自由表面流动问题时不需要对自由液面进行特殊的处理，这与其他发展更早的方法相比，具有明显的优势。因此，自20世纪90年代中期以来，SPH方法迅速地应用于各种涉及自由表面流动的问题，其他学者也在Monaghan模型的基础上提出了许多改进算法。
　　目前，SPH已经被广泛应用于工程和科学中的各类流体力学问题，下面介绍一些代表性工作：
　　（1）波浪与近海和海洋工程结构的相互作用[13-15]；
　　（2）自由表面流体绕流或淹没水工结构的流动过程模拟[16，17]；
　　（3）储罐或容器等罐体中的液体晃动[18]；
　　（4）水射流冲击问题[19，20]；
　　（5）浮船的水动力学问题[21-23]；
　　（6）洪水、海啸、滑坡等自然灾害模拟[24，25]；
　　（7）气液多相流及泥沙冲刷问题[26]。
　　SPH方法得到广泛应用的原因在于它能够轻松处理涉及复杂流动、材料大变形、多相界面等方面的问题。在实际应用中，编程人员对基本算法进行适当调整或改进，就可以处理所面对的流体或固体力学问题，得出与实验和参考数据吻合的结果。而对于传统的数值方法，在不借助商业软件的前提下，使用者需要首先对相对复杂的数值算法有深入的了解，在处理复杂问题时首先解决算法层面的问题。
　　1.4.2 流固耦合
　　SPH方法在流固耦合（fluid-structure interaction，FSI）问题中也有广泛的应用，尤其是在涉及自由表面流与结构变形耦合方

目录
《博士后文库》序言
序
前言
第1章 绪论 1
1.1 数值模拟概述 1
1.2 网格方法及其局限性 2
1.3 弱可压缩SPH方法 4
1.4 SPH方法的数值模拟应用 5
1.4.1 自由表面流动 5
1.4.2 流固耦合 6
1.4.3 多相流 7
1.4.4 表面张力和毛细现象 8
1.5 本书内容梳理 9
1.5.1 物理现象和工程问题 10
1.5.2 本书主要脉络 12
参考文献 13
第2章 SPH方法基础理论 18
2.1 引言 18
2.2 连续介质控制方程 19
2.3 数值近似原理 20
2.3.1 核近似和粒子近似 20
2.3.2 矢量函数散度的粒子近似 22
2.3.3 标量函数梯度的粒子近似 23
2.4 光滑函数和光滑长度 24
2.4.1 常用的光滑函数 24
2.4.2 光滑长度 25
2.5 控制方程的离散化 27
2.5.1 质量守恒方程 27
2.5.2 动量守恒方程 27
2.6 状态方程和人工声速 29
2.6.1 状态方程 29
2.6.2 人工声速值的预估 31
2.7 数值近似精度分析 33
2.7.1 核近似的一致性 34
2.7.2 粒子近似的一致性 34
2.8 数值精度的改进方法 36
2.8.1 修正的光滑粒子法（CSPH格式） 36
2.8.2 有限粒子法 37
2.8.3 修正的光滑粒子法和有限粒子法的对比分析 38
2.8.4 密度场重新初始化 45
2.8.5 核梯度修正 46
2.9 边界处理 48
2.9.1 固壁边界 48
2.9.2 周期性边界 52
2.10 时间积分格式 54
2.10.1 蛙跳法 54
2.10.2 预测校正法 54
本章小结 55
参考文献 55
第3章 多相界面流动的数值算法 57
3.1 引言 57
3.2 两相流模型及算法 58
3.2.1 控制方程及状态方程 58
3.2.2 离散方程及数值技术 61
3.2.3 时间积分格式及算法流程 64
3.2.4 模型验证与结果分析 65
3.3 三相流方程的调整 71
3.3.1 状态方程和界面张力 71
3.3.2 模型验证与结果分析 73
3.4 气-液-固接触角算法 76
3.4.1 模型描述 77
3.4.2 气液固三相接触点 78
3.4.3 虚拟界面构建 80
3.4.4 接触点的确定 82
3.4.5 模型验证与结果分析 83
本章小结 98
参考文献 98
第4章 自由液面的表面张力算法 99
4.1 引言 99
4.2 基本方程 99
4.3 表面张力算法 101
4.3.1 相间作用力模型介绍 101
4.3.2 连续表面力模型 105
4.3.3 表面粒子的鉴定 106
4.3.4 表面曲率和法向量的计算 108
4.3.5 表面张力模型验证 109
4.4 液滴冲击固壁和液面 111
4.4.1 液滴冲击固壁 111
4.4.2 液滴冲击液面 114
4.5 液滴冲击弹性基底建模 116
4.5.1 控制方程 117
4.5.2 离散方程 118
4.5.3 接触算法 120
4.5.4 时间积分格式及算法流程 120
4.5.5 计算结果分析 122
4.6 液滴冲击粗糙基底 129
4.6.1 问题描述 129
4.6.2 模型描述 130
4.6.3 液滴的铺展及回弹过程 132
4.6.4 液滴的刺穿行为分析 135
本章小结 138
参考文献 138
第5章 颗粒冲击损伤问题的数值算法 140
5.1 引言 140
5.2 影响因素及参数定义 140
5.2.1 影响因素 140
5.2.2 颗粒碰撞假设 141
5.2.3 颗粒运动参数及方位角 142
5.3 颗粒运动方程 143
5.3.1 运动方程 143
5.3.2 颗粒建模 144
5.3.3 颗粒形状描述 145
5.4 材料本构模型 148
5.4.1 延性材料 148
5.4.2 脆性材料 150
5.5 弹塑性固体SPH方程 152
5.5.1 控制方程 152
5.5.2 本构模型的实施 153
5.5.3 数值求解流程 154
5.5.4 结果与讨论 156
5.6 粒子分裂及动态细化算法 159
5.6.1 概述 159
5.6.2 粒子分裂准则 159
5.6.3 粒子分裂模式 160
5.6.4 误差分析 161
5.6.5 分裂粒子的变量赋值 162
5.6.6 结果与讨论 162
5.7 多颗粒重复冲击 167
5.7.1 多颗粒冲击算法流程 167
5.7.2 不同形状颗粒的冲蚀特性对比 170
5.8 颗粒冲击实验及模型验证 175
5.8.1 实验装置及参数 175
5.8.2 结果对比与分析 177
本章小结 180
第6章 颗粒-射流耦合冲蚀模型和算法 182
6.1 引言 182
6.2 磨料颗粒的建模 183
6.2.1 颗粒的碰撞过程 183
6.2.2 磨料颗粒模型 185
6.2.3 特征参数 186
6.2.4 颗粒形状 186
6.3 耦合模型建立及求解 189
6.3.1 水射流流体控制方程 190
6.3.2 靶体的控制方程 191
6.3.3 颗粒的运动方程 191
6.3.4 颗粒分布随机算法 193
6.3.5 射流与靶体碰撞的接触算法 196
6.3.6 耦合模型的求解 198
6.4 计算结果与讨论 201
6.4.1 纯颗粒冲击金属表面 201
6.4.2 磨料水射流切割脆性材料 204
6.4.3 磨料水射流切割塑性材料 207
6.4.4 磨料水射流切割塑性和脆性材料的对比 210
本章小结 213
第7章 数值算法的应用扩展 214
7.1 油气水混合液的重力分离过程模拟 214
7.1.1 两气泡上浮融合 214
7.1.2 分散油滴沉降分离过程 217
7.1.3 油水重力分离器模拟 220
7.1.4 油气水三相流动 223
7.2 液滴冲击弹性基底力学行为研究 226
7.2.1 液滴冲击两端固定梁 227
7.2.2 液滴冲击微悬臂梁 232
7.3 角型颗粒冲击延性金属材料的冲蚀机理研究 236
7.3.1 冲击角和方位角的影响 237
7.3.2 冲击速度的影响 242
7.3.3 颗粒的动能损失 246
7.3.4 复杂形状角型颗粒冲击 249
7.3.5 小结 255
7.4 角型颗粒冲击脆性材料的冲蚀机理研究 256
7.4.1 实验与模拟结果对比 256
7.4.2 表面裂纹和损伤凹坑 257
7.4.3 影响因素分析 261
7.5 磨料水射流切割特性研究 266
7.5.1 磨料水射流切割延性材料 266
7.5.2 磨料水射流切割脆性材料 269
7.5.3 与实验结果对比 270
7.6 面向水中爆炸问题的粒子搜索算法 274
7.6.1 数值模型 274
7.6.2 状态方程 277
7.6.3 气液界面处理 278
7.6.4 可变光滑长度 278
7.6.5 计算结果及分析 279
7.6.6 搜索算法优化 285
7.7 喷丸颗粒冲击模拟及表面覆盖率预测 292
7.7.1 概述 292
7.7.2 数值模型 292
7.7.3 单个弹丸冲击 294
7.7.4 大量弹丸冲击 298
7.8 浮体在波浪水槽中的运动行为分析 307
7.8.1 概述 307
7.8.2 基本方程 307
7.8.3 数值水槽模型 309
7.8.4 浮体在水槽中的运动分析 316
本章小结 323
第8章 计算机程序开发 325
8.1 数值计算流程 325
8.2 程序结构及变量定义 326
8.3 粒子的搜索及配对 328
8.3.1 粒子配对、存储及求和 328
8.3.2 邻域粒子搜索 330
8.4 计算实例及代码 334
8.4.1 溃坝流动 334
8.4.2 水弹性耦合 340
本章小结 344
与本书相关的学术论文 346
附录 349
附录1 MLS密度重整化源代码 349
附录2 LU分解的源代码 352
附录3 核梯度修正源代码 356
附录4 链表式搜索源代码 358
编后记 361