本发明公开了一种基于改进的有限差分格式的水流数值模拟方法，可用于计算数学、模拟潮汐流、溃坝、洪水等水利工程等技术领域。
背景技术：
浅水波方程(又称为Saint-Venant方程)被广泛应用于海洋、河流、湖泊等水利工程领域，也可用来数值模拟溃坝、洪水等自然灾害问题。近二十年，浅水波方程的数值求解得到广泛关注并得到很好的发展。数值求解浅水波方程的主要困难在于：一是其属于非线性双曲守恒律方程的一种，其解包含激波、强间断等复杂结构，所以对数值格式的要求很高；二是方程右边的源项很复杂，跟河床的地势有关，所以在静水状态下数值格式需要满足一定的守恒性。
许多学者对此进行研究，并成功的应用于相关的实际问题。1993年，Alcrudo和Carcia-Navarro将有限体积Godunov格式成功应用于求解二维浅水波方程。Bermudez和Vazquez构造了一种迎风格式求解浅水方程，并首次定义了数值格式对于浅水方程必须满足的“C-property”性质。此外，还有Leveque构造的波传播方法，Capilla和Balaguer-Beser设计的有限体积中心格式，Castro等人的Roe格式等。
上述用于求解浅水波方程的基本都是一阶或二阶格式，而近几十年针对双曲守恒律方程的高阶数值求解格式一直是重点的研究内容并得到了很好的发展。1987年，Harten和Osher首次提出了ENO(Essentially Non-Oscillatory基本无振荡)格式并成功应用于计算流体力学等相关领域。ENO格式的主要思想是在所有候选模板中选出最光滑的重构目标单元的值，从而实现高精度和基本无振荡的性质。而因为其存在模板选择过程从而使得计算过程中会造成计算结果的浪费，所以在此基础上构造出了WENO(Weighted ENO)格式。1994年，WENO格式首次被提出并构造出一维情形的三阶有限体积形式。1996年，Jiang和Shu构造了3阶和5阶有限差分WENO格式，并给出了WENO格式的一般计算框架，包括重构多项式的构造，线性权的计算，光滑指示器的计算以及非线性权和最终重构值的计算。此后诸多学者投入到ENO和WENO格式的理论研究和工程应用中。但直到2002年，Vokovic和Sopta才将上述ENO和WENO格式应用于求解一维浅水波方程。并继续于2004年，与Crnjaric-Zic合作将有限体积WENO和中心WENO格式应用于求解浅水波方程和Open-Channel流动方程。Xing和Shu随后构造了满足“C-property”性质的有限差分WENO格式，并通过数值算例验证了其有效性。
有限差分WENO格式构造过程中的最优线性权是唯一确定的，否则会失去数值格式的最优设计设计精度，其具体的值计算相对复杂，尤其在不均匀网格中。因此，Zhu和Qiu于2016年设计了一种线性权可以任取的有限差分WENO格式用于求解计算流体力学中的Euler方程，并且可以保证数值格式在任意线性权下都能保持光滑区域的最优精度。基于此思想，本发明将其应用于求解浅水波方程，并结合“well-balanced”技术保证数值格式的“C-property”性质，最后利用3阶满足TVD性质的Runge-Kutta时间离散方法数值模拟某一区域内的水流状态，并通过圆柱溃坝等数值实验来验证方法的可行性。
技术实现要素：
本发明针对现有技术中的不足，提供一种基于改进的有限差分格式的水流数值模拟方法。
为实现上述目的，本发明采用以下技术方案：
一种基于改进的有限差分格式的水流数值模拟方法，其特征在于，建立空间直角坐标系，用均匀网格剖分计算区域，在形成的计算网格中，构造5阶精度的基本加权无振荡格式进行计算区域内的数值计算模拟，具体包括以下步骤：
步骤一、在水域数值计算区域内，生成直角坐标系并进行网格剖分；
步骤二、在形成的计算网格上构造具有5阶精度的有限差分WENO格式离散浅水波方程组中的空间导数项，其构造过程中的线性权可以任意取值；
步骤三、采用well-balanced方法平衡流通量和河床梯度的变化；
步骤四、采用经典的Runge-Kutta时间离散方法对浅水波方程组进行时间方向上的推进，从而得到任意时刻的水流深度、各个方向的速度以及水流量，从而数值模拟出实际问题中各个时刻的水流状态。
为优化上述技术方案，采取的具体措施还包括：
所述步骤一中，进行网格剖分后的相应网格点坐标记为(xi，yj)，其中x，y表示空间变量，下标i，j表示网格序号。
所述步骤二中，描述水流运动的浅水波方程组为：
将上述方程组统一写成：Ut+F(U)x+G(U)y＝S，其中，t表示时间变量，x，y表示空间变量，U＝(h，hu，hv)T表示守恒变量向量，F(U)，G(U)表示通量向量，F(U)x表示F(U)对x求导，G(U)y表示G(U)对y求导，S＝(0，-ghbx，-ghby)T表示河床效应项，h，u，v，b都是关于时间变量t和空间变量x，y的函数，分别表示水流深度、水平方向速度、竖直方向速度以及河床高度，bx表示b对x求导，by表示b对y求导，T表示转置，g表示重力加速度。
针对该浅水波方程组，线性权可以任取的有限差分WENO格式的构造过程如下：
第1步、对方程组中的空间导数项F(U)x和G(U)y进行差分近似有：
其中，(xi，yj)表示序号为i，j的点的坐标，分别表示x方向和y方向点(xi±1/2，yj)和(xi，yj±1/2)处的数值通量；
第2步、对每个坐标点(xi，yj)处的通量进行Lax-Friedrichs通量分裂：
F(Ui，j)＝F(Ui，j)++F(Ui，j)-，
G(Ui，j)＝G(Ui，j)++G(Ui，j)-，
其中，Ui，j表示U在点(xi，yj)处的函数值，αl表示x方向通量矩阵的最大特征值，αk表示y方向通量矩阵的最大特征值；
第3步、通过F(Ui，j)±的值构造线性权可以任取的新有限差分WENO格式，得到F(Ui+1/2，j)±处的高阶近似值，同理，通过G(Ui，j)±的值得到G(Ui+1/2，j)±的值。
所述第3步具体包括如下步骤：
步骤3.1、先固定y＝yj，取5个单元点的模板T0＝{Ii-2，j，Ii-1，j，Ii，j，Ii+1，j，Ii+2，j}，构造一个四次多项式另取2个包含2个单元点的模板T1＝{Ii-1，j，Ii，j}和T2＝{Ii，j，Ii+1，j}，构造两个线性多项式和分别满足如下约束条件：
将F±(Ui，j)记为则的具体表达式为：
步骤3.2、任取一组和为1的正数γ0，γ1，γ2作为线性权；
步骤3.3、根据三个模板构造出的多项式计算光滑指示器
其中，n表示对应多项式的序号，α表示求和指标，r表示对应多项式的次数，表示多项式对自变量x求α次偏导数；
则的具体表达式为：
步骤3.4、根据任意取值的线性权和光滑指示器，计算非线性权计算公式如下：
其中，表示计算过渡值；
步骤3.5、根据步骤3.1得到的线性项式和步骤3.4得到的非线性权，得到最终F(Ui+1/2，j)±的高阶近似值：
步骤3.6、再固定x＝xi，得到G(Ui，j+1/2)±的高阶近似值。
所述步骤三中，采用well-balanced方法的具体步骤如下：
第1步，分裂浅水波方程组中的河床效应项S如下：
其中，S1，S2表示分裂后的两个源项；
第2步，运用有限差分WENO求解分裂后的两个源项S1，S2；
第3步，得到水流控制方程的半离散有限差分格式：
其中L(Ui，j)＝-F(U)x-G(U)y+S1+S2，表示通过有限差分WENO格式求解后得到的高阶近似值。
所述步骤四中，时间步长为Δt，tn表示第n时间层，表示U(xi，yj，tn)的值，在tn时间层进行计算得到然后利用满足TVD性质的Runge-Kutta时间离散方法离散半离散有限差分格式从而得到时空全离散有限差分格式，如下：
其中，和为计算中间的过渡值，空间方向和时间方向的离散完成后得到的时空全离散有限差分格式为关于时间层的迭代公式，已知初始水流状态的水深h，x方向的水流速度u，y方向的水流速度v，然后根据迭代公式不断循环，即可得到水流在某一时刻或稳定状态时的物理量的高精度数值近似值。
本发明的有益效果是：在均匀的坐标网格中构造线性权可以任意取值的有限差分WENO格式，并结合“well-balanced”技术用于求解具有光滑或间断河床的浅水流问题。与传统的有限差分WENO数值方法相比，改进的五阶有限差分WENO格式具有如下优点：1)线性权可以任意取值，只需满足加起来等于1的正数即可，并且在光滑区域都能达到最优设计精度；2)线性权跟网格无关，任意的网格下都能任意取值而不掉精度；3)收敛性更好，对于定常问题能够快速收敛到机器精度。另外，本发明还结合了“well-balanced”技术，从而能够更好的处理光滑或间断河床的浅水流问题。
附图说明
图1是河床与河流的剖面图。
图2a到图2e是光滑河床水面扰动问题不同时刻的水深等值线图。
图3是间断河床的三维示意图。
图4a-4b是间断河床水面扰动问题不同时刻的水深等值线图。
图5a-5f是平河床圆柱溃坝问题不同时刻的三维示意图。
图6是驼峰河床圆柱溃坝问题不同时刻的三维示意图。
具体实施方式
现在结合附图对本发明作进一步详细的说明。
考虑二维浅水流控制方程，即浅水波方程：
为方便表述，将其简写为：
Ut+F(U)x+G(U)y＝S (2)
其中，t表示时间变量，x，y表示空间变量，U＝(h，hu，hv)T表示守恒变量向量，F(U)，G(U)表示通量向量，F(U)x表示F(U)对x求导，G(U)y表示G(U)对y求导，S＝(0，-ghbx，-ghby)T表示河床效应项，bx表示b对x求导，by表示b对y求导，h，u，v，b分别表示水流深度、水平方向速度、竖直方向速度以及河床高度，且都是关于时间变量t和空间变量x，y的函数，T表示转置，g表示重力加速度，一般取为9.812m/s2。
值得注意的是，除了b(x，y)只跟空间变量有关以外，上述所有的物理量都是关于时间变量和空间变量的函数，即随时间和空间位置的变化而变化。
令xi+1/2-xi1/2＝Δx，yj+1/2-yj-1/2＝Δy，网格中心为网格单元为Ii，j＝[xi-1/2，xi+1/2]×[yj-1/2，yj+1/2]，其中下标i，j为坐标序号。
定义Ui，j表示U在网格单元Ii，j内的值，即Ui，j＝U(xi，yj)。tn表示第n时间层，Δt表示时间步长，则tn+1＝tn+Δt。表示Ui，j在第n时间层的值，下面的过程全在第n层时间进行，所以暂时不考虑时间变量t，即将简写为Ui，j。
首先，针对上述浅水波方程组，构造线性权可以任取的有限差分WENO格式离散通量关于空间变量的导数项，具体过程如下：
第1步，对方程组中的空间导数项F(U)x和G(U)y进行差分近似有：
其中，(xi，yj)表示序号为i，j的点的坐标，分别表示x方向和y方向点(xi±1/2，yj)和(xi，yj±1/2)处的数值通量。
第2步，对每个坐标点(xi，yj)处的通量进行Lax-Friedrichs通量分裂：
F(Ui，j)＝F(Ui，j)++F(Ui，j)-， (5)
G(Ui，j)＝G(Ui，j)++G(Ui，j)-， (6)
其中，Ui，j表示U在点(xi，yj)处的函数值，αl表示x方向通量矩阵的最大特征值，αk表示y方向通量矩阵的最大特征值。
第3步，通过上述F(Ui，j)±的值构造线性权可以任取的新有限差分WENO格式，得到F(Ui+1/2，j)±处的高阶近似值，类似可以通过G(Ui，j)±的值得到G(Ui+1/2，j)±的值，具体步骤如下：
步骤3.1先固定y＝yj，取5个单元点的模板T0＝{Ii-2，j，Ii-1，j，Ii，j，Ii+1，j，Ii+2，j}，构造一个四次多项式另取2个包含2个单元点的模板T1＝{Ii-1，j，Ii，j}和T2＝{Ii，j，Ii+1，j}，构造两个线性多项式和分别满足如下约束条件：
为方便表述，将F±(Ul，j)记为则的具体表达式为：
步骤3.2任取一组和为1的正数γ0，γ1，γ2作为线性权。相比于原来的有限差分WENO格式，这里的线性权可以任意取值而不影响格式在光滑区域的5阶精度，而不是唯一确定并且跟网格有关的正数。
步骤3.3根据三个模板构造出的多项式计算光滑指示器
其中，n表示对应多项式的序号，α表示求和指标，r表示对应多项式的次数，表示多项式对自变量x求α次偏导数。
进一步，的具体表达式为：
步骤3.4根据任意取值的线性权和光滑指示器计算非线性权计算公式如下：
步骤3.5根据步骤3.1得到的多项式和步骤3.4得到的非线性权，得到最终F(Ui+1/2，j)±的高阶近似值：
步骤3.6再固定x＝xi，重复步骤3.1-3.5，类似可得到G(Ui，j+1/2)±的高阶近似值。
其次，为减小河床地势对数值格式的影响，采用的well-balanced技术，具体操作如下：
第1步，分裂浅水波方程组中的河床效应项S如下：
第2步，运用上述的改进的有限差分WENO求解分裂后的两个源项S1，S2；
第3步，进而得到水流控制方程的半离散有限差分格式：
其中，L(Ui，j)＝-F(U)x-G(U)y+S1+S2，表示通过改进的有限差分WENO格式求解后得到的值。
最后，利用满足TVD性质的Runge-Kutta时间离散方法离散半离散有限差分格式，从而得到时空全离散有限差分格式，如下：
其中，和为计算中间的过渡值。空间方向和时间方向的离散完成后得到的时空全离散有限差分格式为关于时间层的迭代公式，已知初始水流状态的水深h，x方向的水流速度u，y方向的水流速度v，然后根据迭代公式不断循环，即可得到水流在某一时刻或稳定状态时的物理量的高精度数值近似值。
下面给出四个算例作为本发明所公开方法的具体实施例。
实施例一、光滑河床的水面扰动问题。
在区域[0，2]×[0，1]内求解二维浅水波方程，河床是一个椭圆状的驼峰，如下：
初始状态的水深和流量分别为：
hu(x，y，0)＝hv(x，y，0)＝0。
在区域x∈[0.05，0.15]内加点小扰动，平静的水面会形成波动往右传播。采用200×100的计算网格，计算时间分别为t＝0.12s，0.24s，0.36s，0.48s，0.60s。图1是河床与河流剖面图，图2a到图2e画出了各个时刻运用本发明方法计算出的水面俯视图，从图中明显可以看出扰动往右传播并经过驼峰河床的效果图。
实施例二、间断河床的水面扰动问题。
在区域[0，2]×[0，1]内求解二维浅水波方程，河床是一个包含间断的驼峰，如下：
初始状态的水深和流量分别为：
hu(x，y，0)＝hv(x，y，0)＝0。
在区域x∈[0.05，0.15]内加点小扰动，平静的水面会形成波动往右传播并经过间断的驼峰河床。我们采用200×100的计算网格，计算时间分别为t＝0.30s，0.45s。图3画出了河床的地形图，图4a到图4b画出了各个时刻运用本发明方法计算出的水面俯视图，从图中明显可以看出扰动往右传播并经过驼峰河床的效果图。
实施例三、平河床圆柱溃坝问题。
计算区域[0，50m]×[0，50m]内有一半径为11m的圆形水坝，水坝内的水深为10m，水坝外水深1m。由于人为或自然因素导致坝体瞬间崩塌，导致水坝内的水迅速往外扩散。运用本发明方法模拟这一过程。计算网格为200×200，计算时间分别为t＝0.0s，0.2s，0.4s，0.6s，0.8s，1.0s。图5a到图5f刻画了这一圆形溃坝问题的过程。
实施例四、驼峰河床的圆柱溃坝问题。
本问题描述的是一个底部具有一个驼峰状河床的圆形溃坝问题，计算区域为[0，2]×[0，2]，计算网格为200×200，计算时间为t＝0.15s。河床的表达式如下：
初始状态的水深及流量分别为：
hu(x，y，0)＝hv(x，y，0)＝0.0.
图6画出了终止时刻的水深以及河床的结果图。
以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。