【Radioss每周干货】SPH 粒子法

SPH是Smooth Particle Hydrodynamics光滑粒子流体动力学的简称，它是基于插值理论的无网格数值方法。

它是用许多与自己周边粒子相互作用的粒子（质点）来描述物体（连续体），遵从连续介质力学的力学定律。

将物体的材料特性、速度和质量以一定的权重分配给每个粒子。在计算过程中，连续性的改变是由于物体的变形引起的。

SPH是将连续体动力学的守恒定律，以偏微分方程的形式，通过核函数W(r,h)将其转化为积分方程。

尽管粒子通常显示为离散点或球体，但在数学上，它们是连续的。

图 1 HV中粒子显示和数学上认为的粒子

核函数W(r,h)是用于描述粒子与其周围半径r范围内的粒子的影响度。

h是光滑长度，它有用户定义，那么粒子与其周围半径为2h内的粒子产生相互作用的影响。h定义的越大，那么粒子就与越多的周边粒子有相互作用，计算精度越高，但是计算资源就越多消耗。

所以需要优化粒子分布和定义适当的h，在精度和性能之间达到最佳折衷。

图 2 使用于SPH的核函数

SPH方法实际上也是一种Lagrangian方法，与传统的Lagrangian方法相比较，使用SPH在计算时不会出现由于单元畸变而引起的计算中没有沙漏能的问题、没有负体积问题、出现更少的时间步长问题。

图 3 球穿击板的SPH模拟方法和Lagrangian模拟方法

SPH特别适用于模拟具有非常大变形的现象，因为此种情况使用Lagrangian网格，单元的变形是非常大的，进而计算时间步长会变得非常小，计算变得非常没有效率，此时就需要换用SPH的方法。

鸟撞，液体晃动，波浪工程，弹道学，气体流动等都可以用SPH方法。

Radioss中的SPH方法与大多数函数兼容。例如，可以模拟两个物体相互作用，一个由Lagrangian有限元离散，另一个由SPH粒子离散。

使用Hypermesh可以轻松创建SPH粒子。可以在1D -> sph进入，也可以从Mesh-> Create-> SPH进入，基于surface和volume创建SPH粒子。

图 4 HyperMesh中的SPH网格划分功能块

Radioss中的粒子可以有三种分布方式：

 · 一种是六边形密网分布法 HCP，

 · 一种是面心立方堆积Face-Centered Cubic FCC

 · 以及还有一种是立方体顶点分布法。

由于HCP和FCC的分布方式的粒子最为紧凑，所以推荐使用这两种分布方式，而且这两种分布方式得到的结果也是非常接近的。

在HyperMesh中(如上图面板中)可以直接得到FCC的分布方式。而HCP的分布方式需要从Altair Connect上下载一个HyperMesh的脚本，然后可以生成相应的HCP粒子分布。 

图 5 SPH粒子的HCP/FCC分布方法

例如下面就是在HyperMesh中设置间距为1mm的FCC分布的粒子：

图 6 HyperMesh中设置SPH粒子的面板

建议一个模型中选用一种类型的粒子分布方式，并且所有粒子的间距（分布密度）要一样，不同尺寸的粒子之间的相互作用可能导致异常，有些情况下接触会发生问题。

粒子的单元属性使用/PROP/TYPE34(SPH)。

如果在HyperMesh中的粒子设置面板（如图 6 HyperMesh中设置SPH粒子的面板）设置了物体的密度”material density”，那么HyperMesh会根据产生的粒子数目自动计算每个粒子的质量mp。

它是这样计算的：

（这里n是粒子数目）

并且自动计算初始的光滑长度h0，也就粒子与其最近的粒子的距离。这个h0是这样计算的：

如果是FCC分布的粒子那么h取值稍微大一些的√2h0，最大可以取2h0。

所有下面的卡片在HyperMesh中定义好粒子间距和构件密度后自动产生的：

图 7 SPH粒子的单元属性卡片

在计算过程中实际的光滑长度h值会随着材料变形引起的密度变化而变化（遵从 p · h = 常量）。使用/H3D/SPH/ DIAMETER输出计算过程中每个时刻的光滑长度h值。

在粒子单元属性卡片中的alpha_cs这个参数是用于保证SPH粒子构件受拉时的稳定性。

alpha_cs的默认值是基于FCC粒子分布时自动设置的，建议受拉时设置alpha_cs=0.1 （也就是10%），但是这样可能会导致能量的吸收，建议在在HyperGraph中检查HE能量。

另外当流体受到负压时可以选用参数xi_stab来控制来稳定计算，这个参数使用时一般建议设置为0.3。

在粒子单元属性卡片中还有qa,qb(Artificial viscosity)人工粘度的参数，它是用于正确模拟冲击和高压缩区域粒子之间相互作用。一般默认使用：

qa = 2

qb = 1

这是用于使用于模拟冲击波/波传播情况，比如冲击速度接近声速在材料中传播的速度时。

在低马赫时，由于内部计算人工粘度时引入过度阻尼，所以系数必须降低，qa，qb应如下设置为不考虑阻尼的情况：

qa = 2.0e-30

qb = 1.0e-30

比如鸟撞中的SPH鸟模型需要设置这个参数。

这里要注意真实的物理粘度必须通过材料卡片设置，而不是在这里设置。

在Radioss中设置粒子的单元属性时还用一个用于定义全局的单元属性卡片/SPHGLO。

这里的Maxsph在v14.0.22版本后就不需要设置了，Radioss自动计算。

Nneigh是定义周边粒子的最大个数，当在定义的2h半径范围内搜索到超过Nneigh的粒子个数，那么就在这些粒子中选取最靠近的Nneigh个粒子当中的Lneigh个例子计算相互作用，所以Nneigh ≥ Lneigh。通常使用默认值就可以了。

Alpha_sort是用于计算搜索相邻粒子的安全半径，也就是每个粒子会搜索到大于安全半径内的相应粒子数，也就是扩大了（1+Alpha_sort）的搜索半径，默认的Alpha_sort=0.25。这个参数最好不要随意改动，增加这个参数，由于有更多的相邻粒子包含进去，会大大降低计算速度。

在Radioss中SPH粒子法模拟的物体可以使用绝大多的可以用于3D实体单元的材料卡片，具体可以参见工具书关于材料卡片兼容性的章节。

通常我们可以使用LAW4（HYD_JCOOK）来模拟金属材料（SPH粒子），常用于穿击金属板，切削金属等模拟；可以使用LAW6（HYD_VISC）来模拟流体，比如水，油等。

这里注意粘滞性是在材料卡片这里设置，比如在LAW6中设置，而不是在单元属性卡片中的qa,qb中设置。

那么各项异性的材料也可以用SPH粒子法，比如加劲混凝土使用SPH粒子法，由于各个方向上的配筋不同，加劲混凝土需要定义为正交各向异性（使用LAW24），此时需要在/PROP/SPH粒子单元属性的卡片中使用skew_ID定义相应的材料方向,skew中的局部x轴就是材料方向1，如果没有定义skew_ID那么就认为材料是各项同性。

另外材料卡片同样可以配合/EOS卡片来设置状态，使用/FAIL失效卡片来设置材料的失效。比如下面的使用SPH粒子法的金属铝材料卡片（单位制g-mm-mus）。

通常使用SPH粒子法在Radioss中失效并不是用单元删除，而是停止计算粒子与其周边粒子的相互作用。

SPH粒子的接触可以有粒子和粒子之间的接触，这个接触不需要特别的设置，它是一直有的。还有一种是粒子与Lagrangian构件之间的接触，那就可以用/INTER/TYPE7来设置，此时需要注意下面几点：

 · Lagrangian的壳体或者实体需要设置为主面master

 · SPH粒子的间距小于Lagrangian的壳体或者实体的网格，否则不能呈现光滑的接触面，而是会出现下面的凹凸不平的现象。

 · 设置接触时，Gap是粒子中心到壳的中性面的距离，或者粒子中心到实体单元的外表面的距离。

 · 只能设置Istf=0或者1000，这样接触刚度就使用Lagrangian的主面材料的刚度。

还有一种tied接触也可以设置，就是使用/INTER/TYPE2将贴近Lagrangian的主面粒子设置为tied接触中的从点（slave），而且必须设置Spotflag = 4，这样就不会传递旋转自由度。

另外还需要将向距离铁道接触表面2h范围（或更大）的粒子添加到/INTER/TYPE7接触中，以免这些粒子会与Lagrangian构件在变形过程中有穿透。

SPH粒子法的边界设置的确是与普通的Lagrangian中的边界设置有所不同。

比如使用/BCS时，Lagrangian中只要设置最外面的一层点的边界约束就可以，而SPH粒子法中这种约束还不够，最好约束2h范围内的粒子才可以。

通常SPH粒子法中使用/SPHBCS的边界约束。

使用/SPHBCS可以类似/BCS一样选取最外面一层的点。

Radioss计算时会产生出相应的一下虚拟点（ghost particles）， 比如下图中所示，黑色点选取作为边界，那么使用/BCS会使得内部的点（红色）在力的作用下，由于力没有平衡而越过边界。

但是如果使用/SPHBCS那么计算时，Radioss就会产生（蓝色的）虚拟点，在受力过程中有很好的力的平衡，而不会越过边界。

但是如果力太大还是越过了边界，那么参数Ilev可以用来定义这些越过边界的点是直接从计算中剔除，或者使用弹性冲击方程将这些点反弹回来。

有多少个虚拟点会产生那么由/SPHGLO中的Nghost 参数决定。另外还可在/SPHBCS中通过参数type定义粒子是在约束平面上固结（tied）还是滑移（slide）的。

除了上面的边界约束，还可以定义边界出入口，比如模拟注水和放水，只要在水管入口处使用/SPH/INOUT卡片定就可以。

以水管进出水为例，首先在水管的进出口用壳单元划分网格，这里要注意进口处壳单元的法向Normal顺着流体进去方向，而出口处壳单元的法向是逆着流体出去的方向。这些壳单元可以用void的材料和单元属性设置（/MAT/VOID+/PROP/VOID）。

然后在/SPH/INOUT卡片中入口设置Itye=1；而出口设置Itype=2，如果Itype=3那么就是设置没有反弹的出口Non-reflective frontiers (NRF)。如果这种Itype=4那么就是控制流量的出口。

在surf_ID中就使用出入口的壳网格形成的面。

这样就可以使得入口的位置当注入材料的质量等于粒子的质量时，将创建一个粒子。如何计算注入材料的质量？可以通过定义密度vs时间曲线，以及定义材料沿着入口法向的进入速度曲线来计算。所以需要定义I_Rho这个密度曲线和Iv_n这个速度曲线。

而粒子运动到出口位置，通过计算压力来决定多少粒子可以出去，所以卡片中需要定义I_ps这个压力曲线。

如果是选用NRF出口（Itype=3）那么类似Itype=2需要设置压力曲线，并且还需要设置lc这个长度，它是用于计算下面的截止频率fc的，一般取的大一点的数值，比如这个实例中就取整个SPH粒子域里面的水管模型的长度就可以了。

无论是任何一种Itype的设置都需要Disp_Gap来控制粒子间距，一般Disp_Gap可以设置2h就可以了。

如果需要这个水管进出水的SPH粒子模型可以联系我们。

SPH粒子法稳定计算的时间步长同样也有节点时间步长控制和单元时间步长控制。

单元时间步长计算如下：

通常粒子之间间距越大，稳定计算的时间步长就越大。使用单元时间步长控制（/DT/SPHCEL）时Δtsca一般建议使用0.3。

节点时间步长：

（mi是粒子的质量，

而Ki是粒子之间的刚度有核函数定义）

节点时间步长Δtsca控制时一般建议使用0.67。通常推荐使用节点时间步长控制，它比单元时间步长控制计算更加稳定。

在Radioss中还有一个功能是SOL2SPH，也就是在不失效时它就是普通的实体单元Solid，当实体单元失效时，该实体单元就会变成相应的SPH粒子。这样的方式比全部使用SPH粒子的方法计算要更加高效。

而且同时具有在大变形区域由于使用了SPH粒子所以单元畸变，没有负体积等等这些实体单元会碰到的问题，使用SOL2SPH和全部使用SPH的计算精度相当。

那么在定义了SOLSPH后，到底是如何触发实体单元到SPH粒子的呢？有下面几种情况：

1.    当定义了/DT/BRICK/DEL的控制，实体单元的变形导致时间步长得到满足所定义的最小时间步长那么这个实体单元不是删除而是触发成了SPH粒子。

   
   

2. 当实体单元满足所定义的材料失效准则，实体单元不是删除而是触发成了SPH粒子。

   
   

3. 当定义接触（/INTER/TYPE7）时，当任何一个粒子在它的搜索范围（详见前面SPH单元属性）内搜索到了不同于它的其他部件（part）的没有触发的SPH粒子，那么这个每个触发SPH粒子的实体单元就被触发了，里面的粒子就和这个粒子有了接触，因此当两个部件都使用SOL2SPH时，他们直接不需要再特别设置TYPE7的接触了。

当/SPHGLO中的Isol2sph设置为2时，那就不是不同部件（part）而是不同子集（subset）之间计算粒子接触而引起的SOL2SPH粒子触发。

那么SOL2SPH在模型中的设置也是非常简单。

以下面为例，如果定义part 2为SOL2SPH，那么在part 2使用的实体单元属性（/PROP/SOLID或者/PROP/SOL_ORTH）中定义Ndir和sphpart_ID这两个参数，在sphpart_ID中定义SPH的part信息（下面实例中是ID=33的SPH的part），这个SPH的part使用和实体的part一样的材料ID（下面实例中都使用ID=11的材料卡片）。

在SPH的part中需要定义SPH的单元属性（下面实例中即ID=15的SPH单元属性）。只有实体构件的实体单元属性和SPH单元属性就定义好了。

这里实体单元属性中的Ndir是指实体单元变成SPH粒子时，每个实体单元中每个方向的SPH粒子个数。比如下面的六面体单元中当Ndir=3就是六面体每个方向上（r,s,t）有3个SPH粒子产生。如果是四面体单元当Ndir=2那就是四面体每个方向（r,s,t）有2个SPH粒子产生。

当然目前SOL2SPH的使用时注意一些兼容性问题，比如只能用于六面体和四面实体单元，不能用于SPH的inlet和outlet，具体可以参见Radioss的工具书。