无网格弱可压缩SPH数值算法及应用扩展/博士后文库

作者简介

内容简介

第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）问题中也有广泛的应用，尤其是在涉及自由表面流与结构变形耦合方