粒子法SPH防止粒子穿透的方法

1月前浏览223

粒子法SPH（Smoothed Particle Hydrodynamics）在模拟流体、固体和多物质系统时，可能会遇到粒子穿透边界的问题，尤其是在处理自由表面流动、刚性壁面以及不同介质交界处。防止粒子穿透的方法是SPH技术发展中的关键改进之一，以下是一些主要策略及其优缺点：

1、固定边界条件：

方法：在与固体边界的接触面上，为粒子定义合适的力或压力，确保它们受到足够的排斥力以保持在边界内部。这通常通过在边界附近施加一个虚拟粒子层或者直接计算边界对实际粒子的反作用力来实现。

优点：简单直观，可以有效避免粒子穿过物理边界。

缺点：需要精确地确定和计算边界上的力分布，过度约束可能导致模拟失真，特别是在高速冲击或复杂几何形状的情况下。

2、再生核函数（RKF）：

方法：设计特殊的再生核函数，在边界区域增强粒子间的相互作用，从而自然地阻止粒子穿透边界而不需额外的约束条件。

优点：自然地融入到SPH框架中，无需复杂的边界条件处理，能较好地模拟边界附近的流场特性。

缺点：再生核函数的设计较为复杂，可能增加计算成本，并且在某些情况下可能难以达到理想的效果。

3、镜像粒子法：

方法：在实体边界外部创建“镜像”粒子，它们的位置和速度相对于边界是对称的，这样可以保证内部粒子与外部镜像粒子之间的相互作用相当于与实际固体边界的相互作用。

优点：物理意义明确，能够准确反映边界条件，适用于复杂几何结构。

缺点：增加了粒子的数量和计算量，对于动态变化的边界或自适应网格可能会带来挑战。

4、人工粘性或阻尼项：

方法：在粒子接近边界时，添加额外的人工粘性或阻尼项，使得粒子速度快速减小并转向，避免穿透。

优点：实施相对简单，可减少穿透现象。

缺点：引入非物理参数，可能影响流体动力学行为的真实性，尤其是在低雷诺数流动模拟中。

5、粒子与边界交互模型：

方法：采用更加细致的边界处理模型，如基于接触力模型，允许粒子和边界进行局部的碰撞响应计算，模拟真实的力学接触过程。

优点：可以模拟更真实复杂的边界效应，包括弹性碰撞、摩擦等。

缺点：计算成本较高，且算法实现复杂度增加，对实时性和稳定性要求较高。

6、耦合有限元/SPH方法：

方法：将SPH与有限元或其他网格基数值方法结合，利用有限元法处理边界条件严格而准确，而在流体核心区域使用SPH方法。

优点：结合了两种方法的优点，既保证了边界条件的正确处理，又保留了SPH在流体模拟方面的灵活性。

缺点：增加了耦合算法的复杂性，需要解决两者之间的数据传递和同步问题。

综上所述，不同的防止粒子穿透边界的方法各有优劣，实际应用中往往需要根据具体模拟场景、精度需求及计算资源等因素综合选择和优化。随着研究的深入和技术的发展，新的防止穿透策略不断涌现，力求在保持SPH方法优势的同时，提升其在处理边界条件时的表现。

来源：CFD饭圈

Fluent仿真实例-大涡模拟大风吹过圆柱体的噪声

流体流过圆柱体产生的噪声案例描述：空气以69.2m/s的速度吹向直径为1.9cm的圆柱体，用Fluent仿真此时产生的噪声。基于圆柱体直径的Reynolds数大概是90000。其他尺寸参数见下图。对于声学仿真，推荐使用LES湍流模型，因为LES模型求解所有涡旋尺度比网格尺度大的涡旋，能较好预测到噪声。1、启动软件并导入网格1.1启动Fluent软件，选择2D双精度版本，单核求解。1.2导入网格文件“cylinder2d.msh.gz”，网格下载在文章底部。为了改善求解速度，将网格重新读取编录，操作：Mesh->Reorder->Domain在文本窗口中显示Fluent采用了ReverseCuthill-McKee方法进行。2、求解器设置3、模型设置3.1湍流模型-大涡LES模型在2D求解器中，LES模型是隐藏的，就是你打开湍流模型面板是找不到的。在文本窗口中输入下面命名“(rpsetvar'les-2d?'#t)”，键盘回车键。命令输入要英文状态，括号也要输入，还有一点就是不能复制黏贴输入，只能手动敲键盘输入才有效，本人亲测过了，Fluent版本是15.0。再次打开湍流模型，就发现LES已经出现可选了。此时会弹出一个warning提示框，点击OK即可。4、边界条件4.1inlet边界，边界类型为velocity-inlet。4.2outlet边界，边界类型为pressure-outlet。保留默认设置。5、求解设置5.1离散方案设置。5.2松弛因子设置，将pressure松弛因子调到0.7。5.3残差设置。5.4初始化5.5此时，保持case和data文件，命名为“cylinder2dt0.00”。设置计算步数，先计算20步。5.6迭代完20步后，设置面监控升阻力。先设置阻力监控具体参数设置如下。设置升力监控。设置计算参考值。保持case和data文件，文件名“cylinder2dt0.00”，覆盖之前保持的。5.7输入迭代步数4000.点击calculate计算。计算完后流动的时间是t=0.02s。此时，升力系数如下。此时，阻力系数如下。保持case和data文件，命名为“cylinder2dt0.02”。6、设置声学模型操作：Models->Acoustics->Edit...将压力松弛因子设置为1.点击计算，将会继续迭代4000步，计算完后，保持case和data文件，命名为“cylinder2dt0.04”。修改声学模型的参数，操作：Models->Acoustics->Edit...将SourceCorrelationLength设置为0.095m，此值等于圆柱体直径0.019m的5倍。如果是3D模型计算，此值不需要设置。设置声学信号。操作：Solution->RunCalculations->AcousticSignals...设置两个信号接收器receivers。其中一个在Y方向的35D（D为圆柱体直径）就是-0.665m处，另外一个在128D就是-2.432m处。计算和读取声学信号进两个接收器中。7、声学后处理7.1显示声压随时间变化图。选择两个接收器。点击plot按钮，显示声压随时间变化曲线。可以看到，红色线（接收器2）的时间相对黑色线（接收器1）平移到大约5e-3s，而且也弱很多。这是因为接收器2比接收器1离圆柱体更加远。7.2傅里叶转换。操作：Results–>Plots–>FFT。参数设置如下接收器1的声压谱显示如下。来源：CFD饭圈