对圣维南方程组进行离散的方法很多。例如特征线法、有限差分法、有限单元法等。由于近年来电子计算机的广泛应用和其性能的不断提高,数值模拟和计算技术飞速进展,计算方法多种多样,读者可参阅有关文献,这里将简要介绍有限差分解法。
图13.2.4
差分法的基本思想是在空间s和时间t两个方面将本来属于连续性的问题离散化,将连续的解域离散为在网格点上求解,用数值分析中差商逼近微商的方法将连续的微分方程组离散为网格点上的差分方程组,然后求解差分方程组,以得到结点上的未知量。用有限差分法求偏微分方程的数值解时,首先应将偏导项化为相应的差商形式。设人为划分的矩形网格如图13.2.4所示,图中i表示空间位置,n表示时间。设u(s,t)是地点s和时间t的函数,是待求的未知量,则U ni表示定义在网格点A 处地点为s=iΔs、时间为t=nΔt的U 值。Δs和Δt分别为距离和时间的间隔,称为空间步长和时间步长。
使用有限差分法时,首先应将函数U (s,t)在s=iΔs、t=nΔt点的空间偏导数写成如下差商形式:
式中O(Δs)表示与Δs同阶的截断误差。此空间偏导数也可以写成:
可见差商形式是多种多样的。式 (13.2.13)、式 (13.2.14)和式(13.2.15)分别表示前差、后差和中差。式 (13.2.15)中的O(Δs2)表示与Δs2同阶的截断误差。
对函数U (s,t)在s=iΔs、t=nΔt点的时间偏导数的处理方法是类同的,例如采用前差形式时为:
将偏微分方程中每一项偏导数均用相应的差商代替,就能写出该微分方程在每一网格点(i,n)i=1,2,…;n=1,2,…处相应的差分方程式。对一组网格点写出各点相应的差分方程时就形成了差分方程组。求解差分方程组,可解得U 在不同位置和不同时刻的值,这就是偏微分方程的数值解。但是,由于写出差分方程时忽略了差商表达式中的截断误差项,所以解出的U 是差分解,并非偏微分方程的微分解,而是它的近似值。
用不同形式的差商逼近圣维南方程组中的偏导数时,会得到不同的差分方程(差分格式)。被采用的差分格式应满足数值解的相容性、收敛性和稳定性。所谓相容性是指一个微分方程采用某种差分格式化为相应差分方程后,当步长Δs和Δt趋向零时,这个差分方程在任一网格点都应收敛于该微分方程。收敛性则指差分方程的解应收敛于微分方程的解。稳定性是指在逐时段进行求解时,计算中的舍入误差和初始误差(包括初始值和边界条件等的误差)被控制在一个有限的范围内,而不是无限增长。详细地讨论这一问题已超出本书范围,必要时读者可参阅有关参考书。
主要的差分格式有两大类:一种是显格式,另一种是隐格式。显格式指求时段末瞬时网格点上的未知量时可直接从时段初瞬时各网格点上的已知量推得,例如图13.2.5上的未知量可根据已知量,,,…直接求得。隐格式则指未知量不能由已知量,,,…直接求出而必须由含有未知量,,,…组成的方程组联立求解。隐格式虽然求解过程复杂,但稳定性比显格式好,显格式受稳定性要求的约束,时间步长Δt不能取大。隐格式的时间步长不受这方面限制。当然Δt取得太大,精度会受到影响。总之,差分格式种类较多。下面介绍求解圣维南方程组比较成功的几种差分格式。
图13.2.5
1.显式差分格式
(1)扩散格式。在图13.2.5所示的i、n结点上,对s的偏导数用中心差商:
对t的偏导数:
式中a为权重系数,0≤a≤1。当a=1时,相当于对t取向前差商;当a=0时,为纯扩散格式,此时式(13.2.18)变为:
经验表明,取较小的a值,例如a=0.1,常能得到较好的成果。
将式(13.2.17)、式 (13.2.18)分别代入连续方程和运动方程可得相应的差分方程。
代入连续方程式(13.2.8)得:
式中B为河段的平均水面宽度。整理得:
式中 n——粗糙系数。
代入运动方程式(13.2.10),得:
整理得:
利用差分方程式(13.2.19)及式 (13.2.20)求解n+1时刻的水力要素时,n时刻的水力要素已经求出,相当于式右端为已知,求解时自左至右取i=2,3,…可显式地求出n+1时刻各结点上的水力要素值。
图13.2.6
扩散格式的收敛稳定条件是(证明从略):
(2)蛙跳格式。这种格式对时间t和距离s都用中心差商(图13.2.6):
将以上格式代入连续方程(13.2.8)和运动方程(13.2.10)可得相应的差分方程:
利用上面两式求n+1时刻的水力要素zn+1i、vn+1i时,n、n-1时刻的水力要素均已求出,因此,也可自左至右取i=2,3,…逐点显式地求出n+1时刻的水力要素值。
蛙跳格式的收敛稳定条件也为:
以上介绍的是扩散格式和蛙跳格式用于计算内点的情况。边界点的情况与内点不一样,多了一个表达边界条件的方程式。如果在边界上也像内点一样,仍由圣维南方程组建立两个相应的差分方程,则将造成方程式多于未知数的矛盾。因此,边界点的计算需另作考虑。有一种考虑是将连续方程和动量方程综合为一个偏微分方程,然后化为差分方程,这样就减少了一个差分方程,矛盾得以解决。下面以蛙跳显格式为例说明这一方法。
设上游边界条件为水位过程线:z=z(s0,t),则此时边界点只有一个未知数v,而描述水流规律的方程式有两个,多于未知数个数,发生矛盾。斯托克(Stoker J.J.)提出把连续方程和动量方程通过统一量纲合并为一个方程式的办法解决上述矛盾。实践证明这是一个可取的处理方法。(www.daowen.com)
如图13.2.7所示,边界点计算相当于内点计算中所取蛙跳网格点的一半,故作差商逼近时,距离偏导数只能用前差,不能用中差,其他与内点同。考虑到上述特点,经过差分逼近,上式可变换为:
图13.2.7
求得vK后,如尚需求流量QK,则QK=AK(zK)vK(过水断面面积用相应于zK的AK值)。
将上式转换为差分方程,可得:
式中zK及vK都是未知数。可结合边界条件QK=AK(zK)vK用逐渐逼近法求解。
下游边界点的计算可仿此进行,不再重复。
2.四点隐格式
在将微分方程的偏导项化为差商时,形式是多样的。采用不同形式的差商形成的差分方程也不同,判别差商形式好坏的标准应是解的精度、收敛性、稳定性以及计算工作量等。用四点隐格式离散圣维南方程组,上述几方面都是比较理想的。
(1)差分格式。四点隐格式将计算点取在时—空网格的中心点 (见图13.2.8中M点),即差分方程是针对图中的点1,2,…,M,…,N-1来建立的,但变量仍定义在网格点上,设M 点是t=nΔt到t= (n+1)Δt、s=iΔs到s= (i+1)Δs的网格中心点。t=nΔt时刻所有网格点处的参数值均为已知值,t= (n+1)Δt时刻所有网格点处的参数值均为待求量。若对M 点建立差分格式,按四点隐格式的法则,应采用如下的差分格式:
图13.2.8
式中f代表任一参数诸如断面平均流速v、水深h、断面面积A、水面宽度B、渠底坡i及摩阻坡降J等。f的上标n表示时间,f的下标i表示距离。∂f/∂s和∂f/∂t分别是f值对距离和时间的偏导数。
(2)差分方程组。利用式(13.2.25)、式 (13.2.26)及式 (13.2.27),对t=nΔt到t= (n+1)Δt时段所有网格中心点1,2,…,M,…,N-1写出圣维南方程组(13.2.11)相应的差分方程式,对每一个网格中心点,都可得到两个差分方程组,例如对M 点有:
式中的摩阻坡降J,根据谢才公式有:
式中 C——谢才系数;
n——粗糙系数;
R——水力半径。
故式(13.2.30)可写成:
设t=nΔt到t= (n+1)Δt时段从上游边界到下游边界共具有(N-1)个网格,则t= (n+1)Δt横线上共有N 个网格点,每个网格点有两个未知量,即水深h和断面平均流速v,故共有待求的未知量2N 个。由 (N-1)个网格可得出2 (N-1)个差分方程式,尚缺两个方程式,故必须由边界条件提供两个补充方程,例如:上游边界给出流量过程线,即
A——断面面积,是水深h的函数。
下游边界给出水深过程线,即
这样,共有2N 个差分方程式,组成差分方程组,其形式如下:
上边界条件:G0(h1,v1)=0
下边界条件:
为书写简便,以上各式中变量均略去了上标n+1。
根据式(13.2.32)及式(13.2.33):
根据式(13.2.28)和式(13.2.29),式 (13.2.34)中的函数F和函数G 分别为:
式中a、b、c、d、e、p、a′、b′、c′、d′、e′、h′、K′均为与t=n·Δt时刻的水力参数有关的已知常数项。其中:
为书写方便,以上常数项中的水力参数均略去了上标n。
免责声明:以上内容源自网络,版权归原作者所有,如有侵犯您的原创版权请告知,我们将尽快删除相关内容。