[发明专利]含截面变化的明渠非恒定流无网格粒子模拟方法

权利要求书

1.一种含截面变化的明渠非恒定流无网格粒子模拟方法，其特征是，步骤如下：

步骤一，初始化系统的相关变量和运行参数；

步骤二，生成粒子信息；

步骤三，列出求解方程并迭代计算：

求解的一维浅水波方程为：

其中，t为时间，x为位置，A为截面面积，u为水流速度，H为水深，g为重力加速度，S₀为河床坡度，S_f为摩擦项；根据运动流体的物质导数可以将方程(1)和(2)转化为拉格朗日形式：

对于任意的一维光滑函数f(x)，其积分表达形式为：

其中δ(x-ξ)为狄拉克函数，Ω为包含x的区间长度，在光滑粒子流体动力学方法中，若用紧支光滑函数W(x-ξ,h)代替狄拉克函数，则f(x)的近似积分表达式为：

函数导数的积分表达式为：

又因为

[f′(ξ)]W(x-ξ,h)＝[f(ξ)W(x-ξ,h)]′-f(ξ)W′(x-ξ,h) (8)

所以，对于紧支函数W(x-ξ,h)有

为了将积分形式转化为求和形式，进行粒子近似得：

其中N为光滑函数紧支域中的粒子个数，x_j为j粒子的位置坐标，△x_j为粒子对应的长度，二维情况下为粒子对应的面积，三维情况下为粒子对应的体积，导数W′(x-x_j,h)与粒子j相关；

因为

所以，在粒子i处函数梯度的粒子近似式写为：

进一步，在粒子i处，明渠横截面面积函数A的粒子近似式写为：

其中V_j为粒子对应的体积，函数AH导数的粒子近似式写为：

利用式(14)将方程(4)离散，得

其中Π_ij为Monaghan型人工粘度；

通过方程(13)和(15)分别计算出每个时间步下每个流体粒子的面积和速度，进而求解出粒子的高度和位置；

步骤四，输出结果：

①每完成一个时间步的计算就更新其结果，即模拟每一个时刻明渠水流的水深、横截面积和流速；

②完成时间步的循环，输出最终结果。

步骤一所述初始化系统的相关变量和步骤二所述的生成粒子信息，具体包括：

①初始化与问题相关的变量信息和运行参数；

②生成流体粒子信息，在流体域初始化粒子分布，并添加初始信息；

③生成虚粒子信息，根据距离流体域边界2h的流体粒子设置虚粒子的初始信息。

2.如权利要求1所述的含截面变化的明渠非恒定流无网格粒子模拟方法，其特征是，进一步地，所述初始化相关变量信息和运行参数具体设置如下：

计算域是长度为50m的一维空间，大坝位于计算域的中心25m处，0～25m为有水区域，高度为1m，25～50m为干床区域，高度为0，不考虑河底摩擦和河床坡度，即s₀＝s_f＝0；

基于浅水波方程的SPH离散格式(13)和(15)，模拟一维溃坝的拉格朗日粒子方法为：

步骤一，初始化，初始化系统的相关变量和运行参数，具体包括：

计算域长度为50m的一维空间，大坝位于计算域的中心25m处，0～25m为有水区域，粒子间距△x₀＝0.01m，初始粒子高度为1m，25～50m为干床区域,粒子的初始速度都为0m/s，质量为10kg，密度为1000kg/m³，计算时间步长为0.001s，计算时间为5s，采用可变光滑长度进行计算，初始光滑长度h₀＝1.0△x，取三次样条函数作为核函数；

步骤二，生成粒子信息，具体包括：

初始化粒子步骤中，总共生成流体粒子2500个，不包括虚粒子，在计算域的左侧布置虚粒子，保证距左端2h范围内的流体粒子能够进行正常运算，虚粒子的质量、密度、高度和速度与流体粒子相同；

通过方程(13)和(15)分别计算出每个时间步下每个流体粒子的面积和速度，进而求解出粒子的高度和位置，具体计算过程为：

①循环每个时间步；

②初始化计算域中的粒子后，搜索目标粒子的邻域粒子，首先通过方程(13)计算该粒子的面积，其次根据该粒子的面积和已知的明渠形状求解其高度，最后更新该粒子的面积和高度；

③通过方程(15)计算出目标粒子的速度，再通过速度计算粒子位置的变化，最后更新该粒子的速度和位置；

④通过上一步计算出的粒子位置来确定是否更改粒子的类型；

⑤循环计算域内的每一个流体粒子，重复做②③④中的操作。