基于SPH方法的含非齐次纽曼边界污染物输运预测方法

本发明属于计算机流场仿真研究，具体涉及采用光滑粒子流体动力学方法在含非齐次纽曼边界条件的污染物输运问题上进行求解的一种方法。

背景技术：

1、在长距离输水线路中存在众多交叉建筑物和控制建筑物、水中和土壤中的各种可溶性物质、人类频繁活动形成的各种环境物质以及天然河道相互连通等诸多因素都会对水质产生影响。而且在实际问题中，要考虑到各类复杂边界条件对污染物输运过程的影响，这些问题导致工程面临着严峻的水环境治理考验。因此，需要更加准确的模型和方法进行求解，对污染物输运问题进行模拟研究。只有深入掌握污染物在水体中的输运扩散规律，才能确保在不同的外界环境下提出更加切实可行的控污处理方案，这不仅是维护水资源和生态环境的重要措施，也是保障工程项目水质安全和平稳运行的关键。

2、光滑粒子流体动力学(smooth particle hydrodynamics，sph)是一种数值方法，也是一种插值方法，它首先被引入天体物理学领域，现在广泛用于流体数值模拟。sph将流体看作由一群离散的粒子组成，在污染物输运过程中包括入流粒子、流体粒子以及出流粒子组成，每个粒子的物理属性通常有五个：位置、速度、密度、质量和浓度。众所周知，流体粒子的新位置可以通过流体粒子的原位置、流体粒子的速度求出，入流粒子、出流粒子均可以看做特殊的流体粒子，同样需要具有流体粒子所具有的物理属性。通过流体粒子的原浓度和扩散系数确定下一个时间步上该流体粒子的新浓度是保证模拟结果正确性的关键。在sph数值模拟过程中，所研究的流体流动的区域称为问题域。针对问题域规模较大的情况，可以对问题域中一个小区域内的关键性的流体流动问题进行研究，该小区域称为计算域。计算域的大小直接决定着计算量的大小，因此计算域不能无限地扩大，需要控制在一定的范围内。流体从入流边界流入计算域，从出流边界流出计算域，入流边界和出流边界也是计算域和问题域的相交部分。

3、之前的方法[1]成功的处理了污染物输运问题，但在处理过程中并没有考虑非齐次纽曼边界条件，而在实际问题中往往会遇到这类边界条件。但是边界条件的处理一直是sph方法需要攻克的难题。目前广泛采用的做法是在边界附近添加虚拟粒子以模拟边界条件的影响。基于该思想有两种边界处理方法，一种是使用边界外固定虚粒子的方法，其需要利用计算域内粒子对虚粒子进行插值，是一种外推的方法；另一种是使用边界上及边界外虚粒子方法，其可使用边界信息及计算域内粒子对外界虚粒子进行插值。然而，采用上述方法求解带有非齐次纽曼边界的污染物输运问题时，因为虚粒子的存在会导致计算负担增大，计算效率显著降低，另外计算结果的数值精度也不能达到理想的二阶精度，并不能对污染物输运过程进行有效的模拟。

4、参考文献：

5、[1]侯庆志,陈少康,苗春福,等.基于光滑粒子流体动力学的污染物输运预测方法:cn202210560376.2[p].cn202210560376.2[2024-04-12].

技术实现思路

1、本发明的目的是提供一种基于光滑粒子流体动力学的含非齐次纽曼边界的污染物输运预测方法，能够预测污染物在水中的流动过程中浓度发生的变化，得到符合现实的污染物传播过程，在实际工程领域具有指导作用，具有数值精度高、计算效率高等特点。本发明的技术方案如下：

2、一种基于sph方法的含非齐次纽曼边界污染物输运预测方法，包括以下步骤：

3、步骤一，系统参数初始化，参数信息包括计算域的顶点坐标，入流边界的顶点坐标，出流边界的顶点坐标，粒子间距，光滑长度，模拟总时间，时间步长和扩散系数；

4、步骤二，生成粒子

5、1)生成流体粒子，在流体域均匀分布流体粒子，并添加初始流体粒子信息；初始流体粒子信息包括流体粒子密度、流体粒子质量、流体粒子x方向速度、流体粒子y方向速度、流体粒子的类型标识号为1、流体粒子的浓度；

6、2)生成入流粒子，在流体域均匀分布入流粒子，并添加初始入流粒子信息；初始入流粒子信息包括入流粒子密度、入流粒子质量、入流粒子x方向速度、入流粒子y方向速度、入流粒子的类型标识号为2、入流粒子的浓度；

7、3)生成出流粒子，在流体域均匀分布出流粒子，并添加初始出流粒子信息；初始出流粒子信息包括出流粒子密度、出流粒子质量、出流粒子x方向速度、出流粒子y方向速度、出流粒子的类型标识号为3、出流粒子的浓度；

8、4)生成辅助边界粒子，在含有非齐次纽曼边界的边界处均匀分布辅助边界粒子，添加辅助边界粒子信息；辅助边界粒子信息包括辅助边界粒子密度、辅助边界粒子质量、辅助边界粒子x方向速度、辅助边界粒子y方向速度、辅助边界粒子的类型标识号为4、辅助边界粒子的浓度；

9、步骤三，列出描述污染物输运问题的控制方程的离散格式和出入口边界及非齐次纽曼边界的施加方式并迭代计算；

10、步骤四，输出结果：模拟总时间到，结束循环，并输出这一时刻所有流体粒子的位置、速度、浓度，所有入流粒子的位置、速度、浓度以及所有出流粒子的位置、速度和浓度，作为最终结果。

11、进一步的，步骤三的方法如下：

12、1)描述污染物输运问题的控制方程用sph数值方法离散后的格式为：

13、

14、

15、

16、

17、

18、其中，w(xi-xj，h)称为光滑核函数，是核函数的梯度，h称为光滑长度，大于粒子间距，核函数通常选取三次样条函数，用r表示xi和xj之间的相对距离，r＝|xi-xj|/h，则核函数w(xi-xj，h)表达式为：

19、

20、i表示当前流体粒子，其临近2h范围是其支持域，j是当前流体粒子支持域内的其中一个流体粒子，xi和xj分别是当前流体粒子i和其支持域内流体粒子j的位置，n是当前粒子i支持域内所有流体粒子的个数，n代表当前时刻，n+1代表下一时刻；mj是流体粒子j的质量，ρi和ρj分别是当前流体粒子i和其支持域内流体粒子j的密度，分别为计算过程中的中间变量，ci和cj分别是当前流体粒子i和其支持域内流体粒子j的浓度，d表示扩散系数，t是时间步长；

21、在公式(1)～(4)中，为示性函数，其表达式如下：

22、

23、

24、

25、

26、xi、yi分别表示当前流体粒子i在x方向和y方向上的位置分量，xj、yj分别表示支持域内流体粒子j在x方向和y方向上的位置分量；

27、在公式(1)～(4)中，

28、

29、

30、其中，和分别表示当前流体粒子i在x方向和y方向上核函数的梯度；

31、2)出入口边界及非齐次纽曼边界的施加方式；

32、出口边界施加自由出流边界条件；流体粒子流到出口域就改变其流体粒子类型标识号1为出流粒子类型标识号3，成为新的出流粒子，其他属性保持不变；

33、入口边界施加自由入流边界条件；入流粒子流到流体域就改变其入流粒子类型标识号2为流体粒子类型标识号1，成为新的流体粒子，其他属性保持不变；

34、非齐次纽曼边界的施加：通过辅助边界粒子来模拟非齐次纽曼边界条件，非齐次纽曼边界条件是浓度在边界上的导数值，被视为中间变量或辅助边界粒子的和根据已知的非齐次纽曼边界条件赋值计算得到，辅助边界粒子的浓度值与非齐次纽曼边界条件表示为如下单边导数的形式：

35、

36、其中，xa代表辅助边界粒子位置，f(xa)为辅助边界粒子的浓度值，f′(xa)是在边界处的非齐次纽曼边界条件，dx是粒子间距，f(xa+dx)和f(xa+2dx)分别是距离边界辅助点最近和第二近的流体粒子的浓度值信息，等式右侧三阶导数项是公式的截断误差；

37、3)迭代计算。

38、进一步的，所述迭代计算方法如下：当前时刻从0开始，计算过程如下：

39、①判断当前时刻是否达到模拟总时间；是，则结束计算；否，则重复计算过程②～④；

40、②对于计算域中的所有粒子均在初始化时已分配坐标，每个粒子的新位置等于原位置加上粒子本身速度乘以时间步长；对于更新位置后进入流体域的入流粒子，改变其粒子属性中的类型标识号，成为新的流体粒子，并在入口边界最外层的相应位置对入口粒子加以补充；对于更新位置后流出流体域的流体粒子，改变其粒子属性中的类型标识号，成为新的出流粒子；

41、③各流体粒子轮流作为当前流体粒子参与计算，用当前流体粒子在支持域内搜索到的邻近流体粒子和辅助边界粒子，根据方程(1)～(2)计算出当前流体粒子的和另外，辅助边界粒子的和根据已知的非齐次纽曼边界条件进行计算得到。

42、④各流体粒子轮流作为当前流体粒子参与计算，用当前流体粒子在支持域内搜索到的邻近流体粒子和辅助边界粒子，根据方程(3)～(5)计算出当前流体粒子在一个时间步后的新浓度，另外，辅助边界粒子的浓度根据方程(13)计算得到；

43、⑤当前时刻＝当前时刻+1个时间步长；重复计算过程①。

44、进一步的，步骤三中计算过程②所述的在入口边界最外层的相应位置对入口粒子加以补充的方法为：当入流粒子流入流体域时，在流体域内生成一个新的流体粒子，该新的流体粒子继承原来入流粒子的位置、速度、浓度属性；流入流体域的入流粒子其沿入口边界法向方向上的坐标值被改为沿入口边界法向方向上的坐标值减去入口边界的长度，重新生成新的入流粒子。

45、进一步的，步骤二中所述的生成粒子：入流粒子和出流粒子浓度为0。

46、进一步的，光滑长度h为粒子间距dx的1.33倍。

47、本发明的特点及有益效果是：本发明提供了一种计算效率高、数值精度高的处理含非齐次纽曼边界的污染物输运问题的方法，可以更为准确的预测污染物输运过程的浓度变化。