【流体】基于平滑粒子流体动力学和空间哈希模拟相互作用流体颗粒的网格附Matlab代码-优快云博客

本文链接：https://blog.youkuaiyun.com/matlab_daizuo/article/details/145695701

✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。

🍎 往期回顾关注个人主页：Matlab科研工作室

🍊个人信条：格物致知,完整Matlab代码及仿真咨询内容私信。

🔥 内容介绍

本项目旨在利用光滑粒子流体动力学 (Smoothed Particle Hydrodynamics, SPH) 方法结合空间哈希 (Spatial Hashing) 技术，模拟一个由相互作用的流体粒子构成的网格。通过调整流体静止密度 (p_rho)、平滑半径 (h) 和粘度系数 (p_mu) 等参数，能够模拟不同类型的流体。然而，目前该程序存在一些亟待解决的问题，主要表现为粒子受到的力和相应的速度过大，导致所有粒子在运行一段时间后，迅速飞出模拟边界。本文将深入探讨这些问题，分析其潜在原因，并提出可能的调试方向和解决方案。

SPH方法是一种无网格拉格朗日方法，它将连续的流体视为一系列离散的粒子，并通过这些粒子之间的相互作用来模拟流体的宏观行为。每个粒子都携带自身的物理属性，如质量、密度和速度。粒子的相互作用力通过平滑核函数进行加权平均，从而近似计算流体的各种物理量。空间哈希技术则被用于优化邻域搜索，加速寻找每个粒子周围的其他粒子，从而提高模拟效率。

理论上，通过合理设置参数和算法设计，SPH可以有效地模拟各种流体现象，例如水、烟雾、气体等。然而，在实际应用中，数值不稳定性和参数敏感性往往会导致模拟结果出现偏差甚至错误。本项目中出现的粒子飞出边界现象，正是由于模拟过程中产生的力过大所致。

造成力过大的原因可能是多方面的，需要从代码实现的各个环节进行排查：

1. 参数设置不合理:

静止密度 (p_rho):
静止密度是SPH模拟中的关键参数，它直接影响压力梯度的计算。如果p_rho设置过小，会导致压力梯度被高估，从而产生过大的排斥力。反之，如果p_rho设置过大，则会导致压力梯度被低估，可能导致粒子聚集或者出现其他不稳定现象。因此，需要仔细调整p_rho的值，使其与实际模拟的流体密度相符。
平滑半径 (h):
平滑半径决定了粒子之间的相互作用范围。如果h设置过小，会导致邻域粒子数量不足，从而影响压力和粘度力的计算精度。如果h设置过大，则会导致计算复杂度增加，并且可能导致粒子之间的相互作用过于强烈。因此，需要根据粒子间距和模拟精度要求，选择合适的h值。
粘度系数 (p_mu):
粘度系数控制了流体内部的粘性力。如果p_mu设置过大，会导致粘性力过强，阻碍粒子的自由运动。如果p_mu设置过小，则会导致流体流动过于湍急，容易出现不稳定现象。因此，需要根据实际模拟的流体粘度，选择合适的p_mu值。

2. 力计算公式错误:

压力梯度计算:
压力梯度是驱动流体运动的主要力量。本项目可能使用的压力梯度计算公式包括经典的SPH公式、修正的SPH公式 (例如： WCSPH, PCISPH) 等。需要仔细检查压力梯度计算公式的推导过程和代码实现，确保其正确无误。特别需要关注压力梯度中包含的密度梯度项，该项的计算精度直接影响压力的准确性。此外，边界粒子的压力计算也需要特别处理，因为边界附近的粒子邻域信息不完整，容易导致压力梯度计算错误。
粘性力计算:
粘性力主要用于模拟流体的粘性特性。常用的粘性力计算公式基于拉普拉斯算子。同样需要仔细检查粘性力计算公式的推导过程和代码实现，确保其正确无误。尤其需要关注粘性力计算公式中包含的粒子速度差项，该项的计算精度直接影响粘性力的准确性。
其他力 (例如：表面张力):
如果项目中还包含其他力 (例如：表面张力)，也需要检查其计算公式和代码实现，确保其正确无误。

3. 时间步长选择不合理:

时间步长是决定模拟稳定性的关键参数。如果时间步长过大，会导致数值积分误差增大，从而导致能量溢出和不稳定现象。反之，如果时间步长过小，则会导致模拟效率降低。因此，需要根据Courant-Friedrichs-Lewy (CFL) 条件和其他稳定性条件，选择合适的时间步长。常用的时间步长控制方法包括固定时间步长和自适应时间步长。自适应时间步长可以根据模拟过程中的状态动态调整时间步长，从而提高模拟效率和稳定性。

4. 边界条件处理不当:

边界条件的处理对SPH模拟的稳定性至关重要。本项目可能使用的边界条件包括周期性边界条件、固定边界条件、自由表面边界条件等。需要仔细检查边界条件的实现，确保其能够正确地处理边界附近的粒子。例如，对于固定边界条件，需要确保边界粒子能够正确地阻碍流体粒子的运动。对于自由表面边界条件，需要确保自由表面能够保持平滑和稳定。常用的边界处理方法包括鬼粒子法、惩罚力法等。

5. 代码实现错误:

代码实现错误是导致模拟结果出现偏差的常见原因。需要仔细检查代码的各个环节，例如变量类型定义、数组下标访问、循环控制等，确保没有出现错误。同时，还需要使用调试工具 (例如：GDB) 对代码进行单步调试，查找潜在的错误。

调试方向和解决方案：

参数调优:
首先，需要对 p_rho, h, p_mu 等参数进行仔细调优，找到一组能够稳定运行的参数组合。可以尝试使用参数扫描方法， systematically 调整每个参数的值，并观察模拟结果的变化。
力计算公式验证:
重新推导并验证压力梯度和粘性力的计算公式，确保其正确无误。可以使用简单的测试用例，例如单个粒子受到压力梯度的作用，或者两个粒子之间的粘性力作用，来验证力计算公式的正确性。
时间步长调整:
根据CFL条件和其他稳定性条件，选择合适的时间步长。可以尝试使用自适应时间步长，根据模拟过程中的状态动态调整时间步长。
边界条件测试:
对边界条件进行单独测试，验证其能够正确地处理边界附近的粒子。可以使用简单的测试用例，例如让粒子与边界发生碰撞，或者让粒子在边界附近自由运动，来验证边界条件的正确性。
代码审查和调试:
对代码进行仔细审查，查找潜在的错误。可以使用调试工具 (例如：GDB) 对代码进行单步调试，查找程序中的错误，例如内存泄漏、数组越界等。
可视化分析:
将模拟结果进行可视化分析，可以更直观地观察粒子的运动轨迹和速度变化。常用的可视化工具包括 Paraview, Visit 等。通过可视化分析，可以帮助我们更好地理解模拟过程中的问题，并找到解决问题的方法。

总结:

本项目中出现的粒子飞出边界问题，很可能由参数设置不合理、力计算公式错误、时间步长选择不合理、边界条件处理不当以及代码实现错误等多种因素共同导致。解决该问题需要从以上几个方面进行全面排查和调试，并结合可视化分析工具，逐步缩小问题范围，最终找到并解决导致模拟结果出现偏差的根本原因。通过对参数进行调优、验证力计算公式、调整时间步长、测试边界条件、进行代码审查和调试，并进行可视化分析，相信能够最终解决本项目中存在的粒子飞出边界问题，实现稳定可靠的流体模拟。

⛳️ 运行结果

🔗 参考文献

📣 部分代码

Nx = floor(xMax/h);     % This will create the amount of bins in the x direction.Ny = floor(yMax/h);     % Same as above but for the y coordinate.dx = xMax/Nx;           % Actual bin dimensions.dy = yMax/Ny;           % Same as above but for the y coordinate.numBins = Nx * Ny;      % The bins that contain the particles in the system.IDs(1:N) = struct('pos', posVal, 'vel', velVal,'force',forceVal,'den',0,'neigh',[]);BinIDs(1:numBins) = struct('part',particles,'adj',adjBins);