[发明专利]一种基于优化求解平面二维浅水方程的恒定流模拟方法

权利要求书

1.一种基于优化求解平面二维浅水方程的恒定流模拟方法，其特征在于，包括以下步骤：

步骤1，获取研究区域的范围；

步骤2，采用均匀矩形网格对所述研究区域进行网格剖分，共划分得到M×N个网格；其中，x方向网格M个，y方向网格N个，网格编号记为(i，j)；其中，1≤i≤M，1≤j≤N，每个网格尺寸为Δx×Δy，Δx为网格在x方向长度，Δy为网格在y方向长度，每个网格的网格中心存储待求变量；所述待求变量包括流速和水位；

步骤3，设置边界条件、初始条件和计算收敛条件：

所述研究区域在进口边界设置进口流量，在出口边界设置出口边界水位；

所述初始条件包括：每个网格的初始水位设置为所述出口边界水位，每个网格的初始水深设置为初始水位减去河底高程，各个网格的初始流速统一设定为该网格所在断面的断面平均流速；

所述计算收敛条件设定为：连续方程余量二范数与进口流量的比值小于允许值；

步骤4，采用有限体积法离散平面二维浅水方程中x方向动量方程和y方向动量方程，并进行亚松弛处理，得到离散的动量方程；求解所述离散的动量方程，得到各个网格的临时流速；

步骤4中，各个网格的临时流速，具体通过以下方式获得：

步骤41，笛卡尔坐标系下恒定平面二维浅水方程、x方向动量方程和y方向动量方程分别为：

式中：

u为x方向流速，v为y方向流速，h为水深，z为水位，g为重力加速度，Г为紊动扩散系数，为x方向摩阻项，为y方向摩阻项，n为糙率系数；

步骤42，采用有限体积法离散x方向动量方程(2)和y方向动量方程(3)，对流扩散项采用乘方格式离散，水位偏导数项采用中心差分格式离散，摩阻项采用半隐式格式处理，并进行亚松弛处理，得到离散动量方程为：

式中

其中：

P、E、W、N、S为不同相对位置的网格，网格P为中心位置网格，网格P正北方向的相邻网格为网格N，网格P正南方向的相邻网格为网格S，网格P正西方向的相邻网格为网格W，网格P正东方向的相邻网格为网格E；

e、w、n、s分别为网格P与网格E、网格P与网格W、网格P与网格N、网格P与网格S的分界面，记网格P编号为(i，j)，则网格E、网格W、网格N、网格S编号分别为(i+1，j)、(i-1，j)、(i，j+1)、(i，j-1)，分界面e、分界面w、分界面n、分界面s编号分别为(i+1/2，j)、(i-1/2，j)、(i，j+1/2)、(i，j-1/2)；

上标“0”代表上一轮迭代的计算值，为上一轮迭代中网格P处水深、x方向流速和y方向流速，分别为上一轮迭代中网格E、网格W、网格N、网格S处水位；

上标“*”代表临时值，分别为本轮迭代中网格P的x方向临时流速和y方向临时流速，分别为本轮迭代中网格E的x方向临时流速和y方向临时流速，分别为本轮迭代中网格W的x方向临时流速和y方向临时流速，分别为本轮迭代中网格N的x方向临时流速和y方向临时流速，分别为本轮迭代中网格S的x方向临时流速和y方向临时流速；

Δx、Δy分别为网格在x、y方向长度；ω_u为速度松弛因子；D_e、D_w、D_n、D_s分别为分界面e、分界面w、分界面n、分界面s处扩散通量，F_e、F_w、F_n、F_s分别为分界面e、分界面w、分界面n、分界面s处对流通量；为动量方程主对角系数，分别为动量方程中网格E、网格W、网格N、网格S对网格P的影响系数；

步骤43，采用交替方向迭代法分别求解离散动量方程(4)和离散动量方程(5)，得到各个网格的临时流速；

步骤5，构建优化的界面流速插值公式，采用所述优化的界面流速插值公式，计算界面流速；

步骤5，具体通过以下方式得到界面流速：

步骤51，优化的流速修正公式为：

式中：

分别为网格P修正后x方向流速和y方向流速；分别为网格P在x方向假拟流速和y方向假拟流速，式中z'_E、z'_W、z'_N、z'_S分别为网格E、网格W、网格N、网格S水位修正值；B_P、D_P分别为网格P在x方向的流速系数和流速修正系数；C_P、E_P分别为网格P在y方向的流速系数和流速修正系数；

步骤52，根据所述优化的流速修正公式，构建优化的界面流速插值公式为：

式中：

分别为分界面e、分界面w、分界面n、分界面s处的界面流速；分别为网格E、网格W、网格N、网格S处假拟流速；B_E、B_W、C_N、C_S分别为网格E、网格W、网格N、网格S流速系数；

采用所述优化的界面流速插值公式，得到界面流速；

步骤6，采用有限体积法离散平面二维浅水方程中连续方程，计算连续方程余量二范数，利用界面流速与水位之间的优化修正关系，获得水位修正方程，求解水位修正方程获得各个网格的水位修正值；

步骤6具体为：

步骤61，采用的有限体积法离散连续方程为：

式中：

分别为分界面e、分界面w、分界面n、分界面s处修正后界面流速，其与水位之间的优化修正关系为：

式中：

z'_P为网格P水位修正值，D_E、D_W、E_N、E_S分别为网格E、网格W、网格N、网格S流速修正系数；

步骤62，综合公式(10)～(12)，水位修正方程为：

式中

式中：

为水位修正方程主对角系数，分别为水位修正方程中网格E、网格W、网格N、网格S对网格P的影响系数；

步骤62，记网格P编号为(i，j)，连续方程余量二范数||A||计算公式为：

步骤63，采用交替方向迭代法求解水位修正方程(13)，获得各个网格的水位修正值；

步骤7，采用步骤6得到的水位修正值，更新步骤4得到的临时流速和步骤5得到的界面流速；采用步骤6得到的水位修正值，更新对应网格水位；

步骤8，判断连续方程余量二范数与进口流量的比值是否满足步骤3设定的计算收敛条件，如果满足，则迭代收敛，输出恒定条件下流速和水位分布；

如果不满足，则以更新后网格水位和流速作为初值，重复步骤4～步骤7，进入下一轮迭代，直到满足计算收敛条件。