基于无网格人工粘性项处理浅水流动干湿动边界的方法

1.本发明涉及洪水预测领域，尤其涉及一种基于无网格人工粘性项处理浅水流动干湿动边界的方法。


背景技术：

2.浅水波方程(shallow water equations，swes)是研究近海潮汐、河口涌潮及溃坝洪流等现象的重要数值模型，然而二维浅水波方程是由一阶非线性偏微分方程组成的双曲守恒律系统，方程本身具有很强的梯度特性，在求解浅水波方程中存在有许多挑战，如处理不连续水流、处理干湿界面、处理河床高程项和混合流态的模拟等。近四十年来，关于浅水波方程数值模拟的研究开发了各种数值算法，大多数研究者采用的最基本的是基于网格的方法，如有限差分法(finite difference method，fdm)、有限元法(finite element method，fem)和有限体积法(finite volume method，fvm)等。近年来出现的无网格方法，它基于点的近似，部分或彻底地摆脱了网格间的拓扑连接关系。因此它能够避免求解流动问题时时间演化过程中的网格重构，无网格法这一特性使它在处理上述大梯度问题时有着广阔的应用前景。
3.间断性问题是浅水波方程数值模拟中的难点，如果没有一个稳定的方案，当遇到不连续水流时，在间断附近会产生非物理性振荡，导致结果不稳定。
4.随着计算技术的进步，高分辨率激波捕捉方法已经成功地应用于求解浅水波方程的齐次形式。然而，这些数值模型以牺牲算法复杂度为代价来抑制解决方案中的振荡行为。为了解决强间断问题和强激波问题，需要建立一个稳定有效的数值格式，能够减小间断处的非物理性振荡，并在光滑区域获得高阶形式的精确解，同时又能减少计算成本。无网格人工粘性项(meshless artificial viscosity，mav)是一种有效的方法，其通过在间断区域给方程增加适量粘性来最小化间断附近的非物理振荡。


技术实现要素：

5.针对现有数值模拟存在的不足和空白，本发明提出一种基于无网格人工粘性项处理浅水流动干湿动边界的方法。旨在建立一个能够精确有效地求解浅水波方程的数值模型，使得该数值模型有较强的激波捕捉能力，并应用于模拟浅水流动中的干湿变化过程。
6.一种基于无网格人工粘性项处理浅水流动干湿动边界的方法，包括以下步骤：
7.s1、选择合适的粘度参数，对保守形式的二维浅水波方程的间断区间采用无网格法人工粘性项增加粘性，得到具有无网格粘性项的控制方程；
8.s2、采用maccormack法对无网格粘性项的控制方程的时间项进行离散；
9.s3、采用广义有限差分法对无网格粘性项的控制方程不同方向支持域下的空间项进行离散；
10.s4、将s3中离散后的空间项代入s2中对时间项离散后的控制方程，并结合无网格粘性项的控制方程边界条件，求解得到计算域中各计算点的物理量；
11.s5、采用时间推进法处理湿节点，完成干湿边界的模拟运动；其中当计算域中的计算点水深等于或大于预设阈值深度h
thre
时，该节点是湿节点。否则，该节点为干节点。
12.本发明提供的有益效果是：本发明采用广义有限差分法并改变以往的圆形支持域形状，根据前后向差分选取相应的支持域方向对空间项进行离散；采用maccormack方法离散控制方程的时间项；在上述基础上结合所提出的基于无网格法发展的人工粘性方法。本发明的有益效果包括：提高了浅水波方程求解的精度和效率，从而获得较强的激波捕捉能力，并应用于模拟浅水流动中的干湿变化过程。
附图说明
13.图1是本发明方法流程示意图；
14.图2为本发明实施例一中溃坝水流通过下游的底部凸起渠道图；
15.图3为本发明实施例一中的水深数值结果与实验数据的比较图；
16.图4为本发明实施例二中溃坝水流通过底部三个圆丘渠道示意图；
17.图5为本发明实施例二中不同时刻和位置处过渠道下游三个圆丘的溃坝水流结果比较图。
具体实施方式
18.为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地描述。
19.请参考图1，图1是本发明方法流程示意图；
20.本发明提供一种基于无网格人工粘性项处理浅水流动干湿动边界的方法，包括以下步骤：
21.s1、选择合适的粘度参数，对保守形式的二维浅水波方程的间断区间采用无网格法人工粘性项增加粘性，得到具有无网格粘性项的控制方程；
22.需要说明的是，步骤s1中，保守形式的二维浅水波方程具体如下：
[0023][0024]
其中，
[0025][0026]
式中，流动变量为u＝(h,uh,vh)
t
，h为水深，u和v分别为x和y方向上的垂向平均速度分量，g是重力加速度，s
0x
、s
0y
和s
fx
、s
fy
分别为x和y方向上的底坡和摩阻源项。
[0027]
需要说明的是，步骤s1中，无网格法的人工粘性项的形式，表示如下：
[0028]
du＝d
x
u dyu＝μ
x
μyꢀꢀꢀꢀꢀꢀ
(3)
[0029]
式中，d为耗散算子，du代表人工粘性项，d
x
u dyu为对应于两个坐标方向的值；
[0030]
无网格粘性项的控制方程如下：
[0031][0032]
其中
[0033]
其中μ
x
,μy的计算如下所示：
[0034][0035][0036]
ε
xi
＝κiν
xi
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0037]
ε
yi
＝κiν
yi
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0038][0039][0040]
式中，κ是一个正的粘度参数，ε
x
和εy为粘性系数。上标l、r、d和u分别代表所选的局部支持域的方向，下标i表示计算域中的第i个点，ns为第i点邻近点的个数，下标j表示第i点的ns个临近点中的第j点，hi为第i点处的水深，h
i,j
为第i点的ns个临近点中的第j点处的水深。
[0041]
s2、采用maccormack法对无网格粘性项的控制方程的时间项进行离散；
[0042]
需要说明的是，步骤s2采用maccormack法对无网格粘性项的控制方程的时间项进行离散，具体为：maccormack法的预测步骤如下式：
[0043][0044]
maccormack法的校正步骤如下式：
[0045][0046]
s3、采用广义有限差分法对无网格粘性项的控制方程不同方向支持域下的空间项进行离散；
[0047]
需要说明的是，步骤s3具体为，采用广义有限差分法对式(11)、(12)的空间项进行离散，如下：
[0048][0049][0050]
s4、根据无网格粘性项的控制方程边界条件，求解得到各节点的物理量；
[0051]
需要说明的是，步骤s4具体为：
[0052]
无网格粘性项的控制方程边界条件的边界条件为：
[0053][0054]
式中，v＝(u,v)t为速度矢量，n＝(n
x
,ny)t是向外的法向量，τ＝(τ
x
,τy)t是由右手定则定义的切向向量；
[0055]
让计算域的内部点满足控制方程式，边界点满足边界条件，形成稀疏代数方程组，如下：
[0056]
[c]
3n
×
3n
{u}
3n
×1＝{f}
3n
×1ꢀꢀꢀꢀꢀ
(16)
[0057]
求解式(16)得到预测步和校正步各点的物理量hi，uhi，vhi(i＝1,2,
···
,n)。
[0058]
s5、采用时间推进法处理湿节点，完成干湿边界的模拟运动。
[0059]
最后，利用人工粘性方法来精确地捕捉激波，保持数值稳定性，同时能简便高效地解决复杂地形上的湿/干边界变化，精确地模拟浅水流动中干湿边界的移动。
[0060]
本技术再次对模拟浅水流动中干湿边界的移动流程总结如下：
[0061]
step1：在计算域中布置内部点及边界点并设置初始水位、水工建筑物位置及阈值深度h
thre
等条件；
[0062]
step2：设置maccormack预测步循环和校正时间步循环(maccormack法是一种两步预测-校正方法，预测步采用向后差分离散，校正步采用向前差分离散)，当前时间步计算完毕后即可自动进入下一时间步；
[0063]
step3：从第一个点开始判断是否为湿节点，如是湿节点则执行step4，否则执行step10；
[0064]
step4：采用maccormack法对添加人工粘性项的二维浅水波方程时间项进行离散，再对离散后的方程进行空间项离散具体步骤step5～step7；
[0065]
采用maccormack法对添加人工粘性项的二维浅水波方程进行离散，其中添加人工粘性项的二维浅水波方程格式如下：
[0066][0067]
预测步骤：
[0068][0069]
校正步骤：
[0070][0071]
式中：b和f表示后向差分和前向差分，
[0072][0073]
式中：h为水深，u和v分别为x和y方向上的垂向平均速度分量，g是重力加速度，s
0x
、s
0y
和s
fx
、s
fy
分别为x和y方向上的底坡和摩阻源项，并表示为：
[0074][0075][0076]
式中：zb为河床高程，n为曼宁粗糙系数。
[0077]
step5：将第i节点子区域ns(ns为第i节点邻近点的个数)内任意一点的物理量用泰勒级数展开，表达式为：
[0078][0079]
式中：xj表示子区域ns内第j点的位置坐标(xj,yj)；xi表示第i节点的位置坐标；hj表示两个计算点之间的横坐标距离；kj表示两个计算点之间的纵坐标距离；e为物理量；e(xi)为xi节点的物理量；e(xj)为xj节点的物理量。
[0080]
step6：定义残差函数b(e)，表达式为：
[0081][0082]
式中：ω(hj,kj)为计算点物理量对计算点间距离的权重函数。
[0083]
step7：根据移动最小二乘法的思想，使残差函数b(e)对偏导数向量de取最小值，求得节点物理量的一至二阶偏导数，即：
[0084]
[0085][0086][0087][0088]
式中：到分别表示各点上的物理量权重值；e为物理量；j表示子区域ns内第j点；ns为第i节点邻近点的个数。
[0089]
step8：将节点物理量的一至二阶偏导数，代入预测-校正方程组，并结合控制方程及边界条件形成代数方程组[c]
3n
×
3n
{u}
3n
×1＝{f}
3n
×1，式中，[c]为关于变量u的稀疏系数矩阵，{f}为方程的非齐次项。
[0090]
step9：求解方程组得到流动变量u＝(h
n 1
,uh
n 1
,vh
n 1
)
t
，式中，h为水深，u和v分别为x和y方向上的垂向平均速度分量；
[0091]
step10：判断h
n 1
是否大于大于阈值深度，如不是则执行步骤骤s11，如是则可判定为湿节点；
[0092]
step11：湿节点变为干节点；
[0093]
step12：取h
n 1
＝hn，uh
n 1
＝uhn，vh
n 1
＝vhn，并判断是否有相邻的湿节点水位高于该干节点，如是则执行步骤s11，如不是则可判定为干节点；
[0094]
step13：干节点变为湿节点；
[0095]
step14：判断是否已处理所有节点，如是则执行step13，如不是则从下一个点执行step3；
[0096]
step15：判断是预测步还是校正步，如是预测步则执行步骤s16，如是校正步则执行step17；
[0097]
step16：执行校正步后执行step3；
[0098]
step17：判断是否到达终点时间，如是则结束，如不是则进入下一时间步执行step2。
[0099]
以下结合具体实施例对本发明的实现方式进行详细的描述。
[0100]
实施例一
[0101]
如图2所示，本实施例为过渠道下游三角堰的溃坝水流模拟，在闸址下游4m、10m、13m和20m处放置了四个水深测量仪。在干河床条件下，整个通道在闸门开启时是干的，渠道末端是开放边界。
[0102]
本案例计算总点数n＝4575，时间步长
△
t＝0.005s，c＝0.1，干湿转换的阈值水深为hthre＝0.001m。总模拟时间为溃坝后的40s内。图3分别给出了四个测量仪(g4、g10、g13和g20)位置处水深时程图的数值结果。将模拟结果与实验数据进行比较，可以看出模拟结果与实验数据很接近，误差基本较小，并与加权最小二乘局部多项式逼近法(wls local polynomial approximation)和sph的研究成果吻合较好。
[0103]
实施例二
[0104]
如图4所示，本实施例为过渠道下游三个圆丘的溃坝水流模拟。上游区域的初始水
深为1.875m，下游区域最初是干燥的，渠道的曼宁系数设为0.018，整个渠道是封闭的，四边都是无滑移边界条件。本实施例计算总点数n＝25351，时间步长
△
t＝0.001s，c＝0.2，干湿转换的阈值水深为h
thre
＝0.001m。
[0105]
为验证模拟结果的准确性，图5分别给出了t＝8s和t＝28s时沿圆丘水平中轴线y＝7.5m和y＝15m的自由面剖面，并与加权最小二乘局部多项式逼近法(wls local polynomial approximation)的研究成果进行比较，吻合较好，该案例表明该数值模型较好地模拟了地形中陡坡的干湿变化过程，提供稳定可靠的溃坝洪水模拟能力，能够准确有效地模拟复杂地形上的流动现象。
[0106]
本发明的有益效果是：。
[0107]
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。