概述

光滑粒子流体动力学（SPH）方法是一种无网格数值方法。通常的有限元分析中需要定义节点和单元，而该方法用点的集合来描述给定的部件，无须定义单元。在SPH法中这些点通常被称为粒子或拟颗粒。

图1中对比了两种方法。两个离散模型描述的都是瓶子里装的液体。左边的模型是由流体占据的传统四面体网格；在右边，同样的流体体积是由离散点的集合表示的。注意，后者情况下没有网格连接这些点（粒子），它们无需像左边传统的有限元定义多节点单元从而保持连通性。在ABAQUS中除了直接定义SPH粒子外，还支持先定义传统的连续单元，然后在分析开始时或在分析过程中将单元网格自动转换成粒子。

图1 有限单元和SPH颗粒的分布

光滑质点流体动力学（SPH）是一种纯拉格朗日方法，它允许通过插值性质直接离散化一个给定的连续性方程组而无需定义空间网格。SPH的主要优势是无固定网格，对于流体流动、结构大变形和自由表面等难题，该方法处理得相对自然恰当。

SPH的核心并非基于在压缩中彼此碰撞或在张力作用下表现出粘性行为的离散颗粒(球)。相反，它是将连续偏微分方程组巧妙离散化的一种方法，这一点与有限元法非常相似。SPH利用插值来近似域中任意点的场变量值，粒子某个变量值通过对相邻粒子对应的值叠加求和来近似，这些粒子以下角标j来区分，其核函数为W（非零）如下。

SPH的核心是核函数，它可以被理解为一种在一定光滑长度h范围内其他临近粒子对研究粒子影响程度的权函数，如图2所示。其中，光滑长度h决定了对某个点的插值产生影响的粒子数目。

图2核函数示意图

SPH方法自提出（Gingold和Monaghan，1977）以来得到了大量的理论支持，与该方法相关的出版物非常之多。大量的参考文献见下文。

该方法可以使用ABAQUS / Explicit中的所有材料（包括用户自定义材料），可以定义其他拉格朗日模型中的任何初始条件和边界条件，也支持与其他拉格朗日部件的相互接触，从而扩大了该方法应用范围。

对于变形不明显的问题，该方法不如拉格朗日有限元分析准确，对于变形较大的问题，该方法不如耦合欧拉-拉格朗日（CEL）方法准确。当模型的大部分节点都是SPH粒子时，不宜使用多CPU进行分析。

应用

SPH分析对涉及极端变形的问题非常有效。例如：液体晃动，波浪工程，弹道学，喷涂，气体流动，二次撞击造成的闭塞和破碎等。在许多情况下，CEL法和SPH法都可以使用。在许多CEL分析中，材料所占体积明显小于空穴，因而计算成本可能很高。此时，SPH分析是首选方法。例如，在较大空间内追踪初次撞击产生的碎片直到二次撞击发生时，若采用CEL法分析则计算成本很大，而采用SPH分析不存在任何额外成本。

人工粘度

SPH中的人工粘度与有限元方法中的体积粘度（bulk viscosity）意义相同。和其他拉格朗日有限单元一样，粒子单元使用线性和二次粘度来抑制计算响应中的高频噪声。在少数情况下，ABAQUS的默认值不合适，你可以通过关键字更改人工粘度。

Input 输入文件使用：          

使用下列语句来指定线性和二次人工粘度的比例因子：

*SECTION CONTROLS

, , , 线性粘度比例因子, 二次粘度比例因子

初始条件

“Initial conditions in Abaqus/Standard and Abaqus/Explicit,” Section 34.2.1中讲述了所有可用于Explicit分析的初始条件。与力学分析相关的初始条件均可用于SPH分析。

边界条件

边界条件的定义请参照 “Boundary conditions in Abaqus/Standard and Abaqus/Explicit,” Section 34.3.1.

载荷

“Applying loads: overview,” Section 34.4.1描述了所有可用的显式动力分析的载荷类型。节点集中载荷可用，重力载荷是SPH中唯一可使用的分布式载荷。

材料选择

任何在Abaqus/Explicit中可用的材料模型均支持SPH分析。

单元

SPH分析通过单节点PC3D单元实现。这些单元通过定义粒子在空间中的位置来构造研究的部件。

定义PC3D单元的方式类似于定义质点群。分布在模型表面或内部的粒子的坐标与有限单元节点的坐标含义相似。为使计算更加精确，应尽可能均匀地将这些粒子布置在空间中。

除了直接定义PC3D单元外，还可以先定义传统连续介质有限单元（C3D8R, C3D6, C3D4），然后在分析开始时或在分析过程中将单元网格自动转换成粒子。具体讨论见“Finite element conversion to SPH particles,” Section 15.2.2。

ABAQUS / Explicit中的SPH方法是通过三次样条插值多项式实现的，其建立在经典的光滑粒子流体动力学理论之上，参考见下文。

SPH方法不支持二维分析。轴对称可以利用楔形扇区和对称边界条件模拟。PC3D单元不支持沙漏或畸变控制，单元间没有面或边连接。

SPH核插值器

默认情况下，ABAQUS / Explicit中的SPH分析是通过一个三次样条插值多项式实现的。此外，还可以选择二次（Johnson et al，1996）或五次（Wendland,1995）插值。

ABAQUS中采用的是经典的光滑粒子流体动力学理论，当然还可以采用平均流修正的XSPH方法（参见Monaghan,1992），也可用 Randles 和Libersky修正的归一化SPH（NSPH）方法。

以上这些控制选项将在“Using section controls for smoothed particle hydrodynamics (SPH)” in “Section controls,” Section 27.1.4中讨论。

计算粒子体积

目前还无法自动计算粒子占据的体积。你需要指定一个特征长度用于计算粒子的体积，继而计算粒子的质量。假设节点均匀地分布在空间中，每个粒子都与一个以粒子为中心的小立方体相关联。当堆叠在一起，这些立方体将充满整个部件。特征长度是立方体边长的一半。从实际角度来看，一旦创建了节点，就可以将两个节点间距离的一半作为特征长度。此外，如果部件的质量和密度已知，可以利用部件的体积除以粒子总数得到每个粒子所关联的小立方体的体积。当特征长度是小立方体体积立方根的一半时比较合理。如果你将模型定义数据输出到数据（.dat）文件中，那么你可以查阅到模型中相关集合的质量。(参见“Model and history definition summaries” in “Output,” Section 4.1.1).

Input 输入文件使用：    

使用下列语句来定义SPH部件:

*ELEMENT, TYPE=PC3D, ELSET=particle_body

element number, node number

Repeat the data line as often as necessary.

*SOLID SECTION, ELSET=particle_body, MATERIAL=material_name

characteristic length associated with particle volume

光滑长度计算

虽然每个粒子单元在模型中只使用一个节点定义，SPH方法在球形影响区域内计算相邻粒子的贡献。这个球形影响区域的半径即为光滑长度。光滑长度与上文所说的特征长度相互独立，它控制着SPH法的插值性质。默认情况下，光滑长度自动计算。变形过程中，粒子产生相对运动，因此与某个粒子相邻的粒子会（通常）发生变化。在每个Explicit增量步中，基于粒子云重新计算内部相连关系和运动学参数（如正应变和切应变，变形梯度等），然后应力以与减缩积分单元类似的方式计算，继而基于SPH准则反算出粒子云中单元的节点力。

默认情况下，ABAQUS / Explicit在分析开始计算出一个光滑长度，使得在每个影响区域内的平均粒子数是30-50个。光滑长度在分析过程中保持不变。因此每个单元内粒子的平均数可能在分析过程中增加或减少，其取决于模型主要是压缩还是膨胀，如果分析以压缩为主，那么与某个单元关联的粒子总数可能超过所允许的最大值，导致分析终止。默认情况下，每个单元关联的粒子数最大为140。

控制选项在“Using section controls for smoothed particle hydrodynamics (SPH)” in “Section controls,” Section 27.1.4.中讨论。

SPH区域

在分析的开始，ABAQUS会计算出一个矩形区域，在边界内的粒子将被追踪。这个固定的矩形区域比模型总尺寸大10%，其中心位于模型的几何中心。随着分析的进行，如果一个粒子飞出区域，它将成为一个自由飞行质点，不再参与SPH计算。如果之后粒子又重新飞入区域内，它将再次参与计算。

你可以修改区域边界的尺寸，具体讨论在 “Using section controls for smoothed particle hydrodynamics (SPH)” in “Section controls,” Section 27.1.4.

约束

由于PC3D是拉格朗日单元，其节点也可以定义与其他单元、连接体等的约束。由于这些单元没有面或边，基于面定义的单元不能使用PC3D单元。因此，基于单元面定义的约束（如fasteners）不能用于粒子单元。

相互作用

使用粒子定义的部件可以与其他有限单元定义的部件建立接触关系。接触关系类似于基于节点的面（与粒子相关联）和基于单元的面（或解析刚体）之间的接触关系。SPH支持通用接触和接触对。任何基于节点的面定义的接触类型和准则都可用，包括cohesive behavior。支持定义不同的接触特性。默认情况下，粒子不是通用接触域中的一部分，这与一维单元（1-node elements）类似（例如point masses）。默认的粒子接触厚度与在截面定义中指定的特征长度值相同；因此，为了实现接触，粒子的接触行为等效成一个球体，其半径等于小立方体（与粒子相关联的）的内接球的半径。

对于PC3D单元，接触厚度不能为0，否则计算难以收敛。推荐使用默认的或指定合理的接触厚度。

允许在使用PC3D单元定义的不同部件间建立接触关系。但是，只有当相互碰撞的SPH部件是由相同的流体材料组成时，接触才有意义。例如水滴落入装水的桶中。在与固体相关的应用中，比如子弹穿透装甲板，建模时一个部件必须使用常规有限元定义。

粒子和欧拉体之间不能定义接触。

Input 输入文件使用：

使用下列语句来定义网格单元或解析刚体面与粒子定义的面之间的接触:

*CONTACT

*CONTACT INCLUSIONS

node-based particle surface, element-based/analytical_surface

输出

PC3D单元能够输出连续介质单元所有的与力学相关的参数：应力，应变，能量，状态变量，场，用户自定义变量。Abaqus/Explicit分析中的所有节点变量均可输出。

限制

SPH法有以下限制:

· 该方法在变形不太严重,单元不扭曲时不如一般的拉格朗日有限元分析准确，在较大的变形时不如耦合欧拉-拉格朗日分析准确。SPH方法主要用于传统有限元方法或CEL法使用受限时，或者计算成本很高时。

· 当材料处于拉应力状态，粒子的运动将变得不稳定，导致所谓的拉伸失稳。这种不稳定性，与标准SPH法的插值技术密切相关，在模拟固体拉伸时尤其明显。这使得粒子趋向于聚集在一起而表现出破碎行为（fracture-like behavior）。

· 使用粒子单元定义的部件的质量分布和使用连续单元（如C3D8R）定义的部件的质量分布有所不同。当使用粒子单元时，部件中所有粒子的体积都是相同的。因此，部件中所有粒子的节点质量是一样的。如果节点的布置不是一个常规的立方排列，质量分布将会不准确，在部件的自由表面上尤其明显。

· 面载荷不能用于PC3D单元。但是，分布载荷，例如压力，可以先施加在其他有限元表面然后通过接触关系作用到粒子单元上。

· 使用粒子单元定义的部件，如果它们没有使用相同的截面属性，那么它们不能定义接触关系。因此，你不能使用SPH法构建出不同材料的混合体。

· 目前Abaqus/CAE不支持SPH，可以使用Abaqus/CAE先生成质量单元（mass elements），写入输入文件（.inp），然后手动编辑输入文件，将质量单元转换为粒子。也可以利用ABAQUS将单元转化成粒子，参照 “Finite element conversion to SPH particles,” Section 15.2.2。此外，也可以先创建C3D8R单元，然后利用脚本将单元转换为粒子，参照www.3ds.com/support/knowledge-base。

· 对于通过一个实体截面定义的部件，重力荷载和质量缩放必须施加在该实体截面的所有单元上，不能仅定义到某个单元子集上。

对于具有下列特征之一的模型，ABAQUS将采用单线程进行计算：

· 分析过程中，单元根据相关准则逐渐转化为粒子时

· PC3D单元对应着多个实体截面

· 采用NSPH时

·  预定义场变量（包括温度）随材料属性变化

如果使用多线程（CPU）进行计算，SPH分析将受到以下限制：

· 不支持SPH从面节点的接触输出（Contact output）。

· 不支持单元的历史输出。

· 不支持能量历史输出，模型总体能量除外。

· 多线程运算时，不能启用动态负载平衡（Dynamic load balancing）功能

· 建议每个并行域至少有10000个粒子。

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删