剖面二维非恒定悬移质泥沙扩散方程的数值方法

1、剖面二维非恒定悬移质泥沙扩散方程的数值方法1 引言 数学模拟方法正在成为 研究 河流泥沙 问题 的重要手段。 目前 ，一维数学模型 发展 较成熟，已广泛 应用 于模拟长河段的长期变形，但它只能给出河段平均冲淤深度的沿程变化，如需了解短河段的河床变形细节，则要采用二维以至三维数学模型。不论是一维数学模型还是平面二维维数学模型，都不能反映含沙量沿垂线的分布状况，并忽略了含沙量沿垂线分布对垂线平均含沙量变化过程的 影响 。要解决这类问题，必须建立剖面二维数学模型。这种模型主要通过解剖面二维泥沙扩散方程来研究悬移质泥沙沿水深的分布及含沙量的变化过程，对水电站进口和其它引水工程的引水口高程的确定都能提供

2、较好的数值模拟。 等都做了有益的尝试；求扩散方程的数值解曾经因为缺乏高效率的计算工具而难以实现，直到 60 年代后，随着 计算机 的广泛应用，在各种复杂边界条件下求扩散方程的数值解不但成为可能，而且得到迅速的发展，在这方面，曹志先、崔侠 等做了大量 工作 ，取得了很多成果。 数值方法相对于解析方法在求解偏微分方程上有着明显的优势，即简单灵活、计算方便快捷，但要寻找一种精度高、稳定性好、计算方便的差分格式也并非易事。本文拟在前人研究的基础上着重讨论剖面二维泥沙扩散方程的数值解问题，希望能提供一种精度高、稳定性好、计算方便的数值解。 2 基本方程 剖面二维泥沙扩散方程的形式为 (1) 式中 x,y

3、 为水流方向和铅直方向的维轴； u,v 分别为沿水流方向和铅直方向的时均流速； sx , sy 分别为水流方向和铅直方向的泥沙扩散系数； ,S 分别为泥沙静水沉速和含沙量。 对于式 (1) 的求解，研究者一般会对它进行不同程度的简化，为此引入以下假定中的一种或几种： A. 非恒定流可以概化为梯级式恒定流，即 ; B. 在一个时段内，认为泥沙运动可以概化为处于恒定状态，即 ; C. 在二维流动中，纵向扩散系数与方程其他项相比，可以忽略不计，即认为方程右端第一项可以忽略； D. 认为悬移质泥沙粒径均一，即 =const; E. 认为水流为二维均匀流，即 v=0 。 为简单起见，我们讨论的范围限于水

4、流条件为二维非恒定均匀流，悬移质泥沙粒径均匀，为此引入假定 C 、 D 、 E 。这时，泥沙扩散方程为 (2) 目前，对 s 的变化 规律 研究得不很充分，一般假定 s = m (3) 其中 m 为动量传递系数，为修正值。由勃兰特尔掺长 理论 可得 m = u*y/(1-y/h) (4) 式中为卡门常数， u* 为摩阻流速。 对于 u ，我们取卡曼 - 勃兰特尔对数流速分布公式 (umax-u)/u*=(1/ )ln(h/y) (5) 令 W= + ( u*/h)(1-2y/h), 则式 (2) 可变形为 (6) 3 差分方程 3.1 网格的剖分 为建立差分方程，首先必须剖分网格。我们取时间步

5、长 t= ， X 方向的空间步长 x=h1 ， Y 方向的空间步长为 y=h2, 这样形成如下网格 Dh=(xj,yl,tn) xj=jh1,yl=lh2,tn=n 其中 j=0,.,N; l=0,.,M;n 0; 3.2 构造差分格式 通过对流方程和扩散方程的差分格式的构造，我们可以得到对流扩散方程的差分格式。由于隐式格式稳定性好，考虑 Crank-Nicholson 型隐式格式。为此，引入差分算子记号 2 y Snj,l=Snj,l+1-2Snj,l+Snj,l-1 LxSnj,l=Snj+1,l-Snj-1,l LySnj,l=Snj,l+1-Snj,l-1 为了看得更清楚，暂且取 h1

6、=h2=h. 对式 (6) 离散，则 C-N 格式为 Sn+1j,l-Snj,l/ +u/4hLx(Snj,l+Sn+1j,l)= s /2h2 2 y (Sn+1j,l+Snj,l)+W/4hLy(Sn+1j,l+Snj,l) (7) C-N 格式的精度是二阶的，绝对稳定。但对于二维 问题 ，由 (7) 导出的方程组，其系数矩阵不是三对角矩阵，不能用追赶法求解。因此，考虑构造交替方向的隐式格式 ( 命名为 Z-C 格式 ) (Sn+1/2j,l-Snj,l)/( /2)+u/2hLxSnj,l= s /h2 2 y S n+1/2 j,l +W/2hLySn+1/2 j,l (8) (Sn+

7、1j,l-Sn+1/2j,l)/ /2+u/2hLxSn+1j,l= s /h2 2 y Sn+1/2j,l+W/2hLySn+1/2j,l (9) 可以看出， 计算 Sn+1j,l 是由两步组成的，每一步仅是一个方向的隐式，故用两次追赶法即可。 3.3 精度 分析 现在，我们考虑 Z-C 格式的精度。先设法消去过渡值 Sn+1/2j,l ，为此，将 (8) 和 (9) 两式相加，可得 Sn+1j,l-Snj,l+ u/2hLx(Sn+1j,l+Snj,l)= s /2h2 2 y Sn+1/2j,l+ W/4hLxSn+1/2j,l (10) 将 (8) 和 (9) 两式相减，可得 Sn+1

8、/2j,l=(Sn+1j,l+Snj,l)/2+ s /4h2 2 y (Snj,l-Sn+1j,l)+ W/8hLy(Snj,l-Sn+1j,l) (11) 把式 (11) 代入 (10) ，变形整理，可得 (Sn+1j,l-Snj,l)(1- u s /8h3Lx 2 y - 2 uW/16h2L x Ly)=(Snj,l+Sn+1j,l)( u s /2h2 2 y + W/4hLy- u/4hLx) (12) 设 S(x,y,t) 是 (12) 的精确解，并假定 S(x,y,t) 关于 t 三次连续可微，关于 x,y 四次连续可微，那么利用 Taylor 级数展开可得 S(xj,yl,

9、tn+1)-S(xj,yl,tn)/ (1- 2 s /8h3Lx 2 y - 2 uW/16h2LxLy)-S(xj,yl,tn+1)-S(xj,yl,tn)/ ( s /2h2) 2 y + W/4hLy- u/4hLx)=O( 2 +h2) (13) 由此可见， Z-C 格式具有二阶精度。 3.4 稳定性分析 现在，我们来讨论 Z-C 格式的稳定性。为此，把式 (12) 变形整理得 (1- s /2h 2 2 y - W/4hLy)(1+ u/4hLx)Sn+1j,l=(1+ s /2h 2 y 2 + W/4hLy)(1- u/4hLx)S n j,l (14) 由式 (14) 可得出

10、过渡因子为 G( ,k)=(1-2 s /h 2 sin2k2h/2+i W/2hsink2h)(1-i u/2hsink1h)/(1+2 s /h2sin2k2h/2-i W/2hsink2h)(1+i u/2hsink1h) (15) 令 a=2 s /h2sin2k2h/2,b= W/2hsink2h,c= u/2hsink1h, 则 |G|2=|1-a2-b2+2bi/(1+a)2+b2|2 |1-c2-2ci/1+c2|2=1 (1-a2-b2)2+4b2/(1+a)2+b22=1-4a2+4a+4b2+4a3/(1+a)2+b22 (16) 显然，对于任意的 ,h,|G( ,k)|

11、21, 所以 Z-C 格式是绝对稳定的。 我们考虑初边值 问题 。 (1) 初始条件 用 Rouse 公式给出含沙量沿垂线分布 S(x,y,0)=Sa 5(a/h-a)z(h-y/y)z (17) 式中 z= /ku* 为悬浮指标， Sa 为近底含沙量， h 为水深，一般取 a=0.01 0.05h 。 (2) 水面条件 (y=h) (18) (3) 底部边界条件 (y=a) (19) 式中 Sa* 为近底挟沙力，即输沙平衡时的近底售沙量 Sa. (4) 近底含沙量计算近底含沙量在求解泥沙扩散方程时具有边界条件性质，它选取的正确与否，意味着所给边界条件是否正确。实际工程中一般缺乏实测资料，近底

12、含沙量不易测定。这里，我们利用水流挟沙力和含沙量沿垂线分布公式来反求近底含沙量 。 已知断面平均挟沙力为 S*=k(u3/gh )m (20) 假定 Sa*= S* (21) 输沙平衡时，含沙量沿垂线分布用 Rouse 公式 (17) 表示，用 (17) 表达挟沙力的垂线分布 ，然后沿垂线积分得断面平均挟沙力为 (22) 将 (22) 与 (21) 比较，可得 (23) 4.2 计算步骤 为方便计算，将式 (8) 和 (9) 式变形整理，并对 X ， Y 方向取不同的空间步长 ( s /2h22- W/4h2)Sn+1/2j,l+1-( s /h22+1)Sn+1/2j,l+( s /2h22

13、- W/4h2)Sn+1/2j,l+1= u/4h1LxSnj,l-Snj,l (24) - u/4h1Snj-1,l+Snj,l+ u/4h1Snj+1,l= s /2h22 2 y Sn+1/2j,l+ W/4h2LySn+1/2j,l+Sn+1/2j,l (25) 在一个时间层 ( 第 n 层 ) 内，计算分两步进行： 第一步，对式 (24) 用追赶法求第 n+1/2 层的过渡值。 令 C1=- u/4h1,C2=1,C3= u/4h1,E1= s /2h22,E2= W/4h2 D1= s /2h22- W/4h2,D2=- s /h22-1,D3= s /2h22+ W/4h2 l=

14、1 时， D1Sn+1/2j,0+D2Sn+1/2j,1+D3Sn+1/2j,2=C3LxSnj,1-Snj,1 l=2 时， D1Sn+1/2j,0+D2Sn+1/2j,1+D3Sn+1/2j,2=C3LxSnj,1-Snj,1 l=M 时 ,D1Sn+1/2j,M-1+D2Sn+1/2j,M+D3Sn+1/2j,M+1=C3L x Snj,M-Snj,M 令 Hl=-Snj,l+C3LxSnj,l (1lM) H1=-Snj,1+C3LxSnj,1-D1Sn+1/2j,0 HM=-Snj,M+C3LxSnj,M-D3Sn+1/2j,M+1 其中 Sn+1/2j,0 和 Sn+1/2j,M+

15、1 由边界条件给出，则用矩阵形式表示为 （ 26 ） 第二步，再对 (25) 式用追赶法求第 n+1 层的值 令 Fj=E1 2 y Sn+1/2 j,l +E2LySn+1/2j,l+Sn+1/2j,l (1jN) F1=E1 2 y Sn+1/21,l+E2LySn+1/21,l+Sn+1/21,l-C 1 S n+1 0,l FN=E1 2 y Sn+1/2N,l+E2L y Sn+1/2N,l+Sn+1/2N,l-C3Sn+1N+1,l 其中 Sn+10,l 和 Sn+1N+1 ， l 由边界条件给出，则同理可得矩阵方程 (27) 这样，按此步骤一层层地计算。 4.3 数值模拟合理性

16、分析 受所掌握的实测资料的限制， 目前 尚无法对本文提出的算法与含沙量沿垂线分布的实测值进行对比。我们用库里阿雷克沉沙池 的实测资料作了垂线平均值沿程变化的比较。该沉沙池的主要数据为：池深 h=1.53m; 平均流速 u=0.12m/s; 泥沙沉速 =0.0176cm/s; 悬浮指标 Z=0.01 。计算时取卡门常数 =0.4 ， a=0.05h 。表 1 给出了计算值和实测值，结果表明，计算值和实测值比较符合。 表 1 断面平均含沙量验证 ( 单位： kg/m3) Verifications of cross-sectional average sediment concentrations

17、 距离 (m) 0 200 400 600 800 1000 1200 计算值 3.012 2.573 2.071 1.639 1.209 0.849 0.539 实测值 3.012 2.650 1.810 1.520 1.141 0.822 0.561 为了进一步分析含沙量垂线分布计算结果的合理性，我们对另一组较粗的泥沙 ( =0.616cm/s,Z=0.4) 进行了对比计算。图 1, 图 2 分别为两组沙的计算结果。从图中可以看出，计算结果符合含沙量沿垂线分布的一般 规律 ，粗沙分布不均匀，细沙分布较均匀；近底浓度相对较大，水面浓度相对较小，不存在 Rouse 公式中水面含沙量为 0 的缺

18、陷；含沙量沿程衰减的特性较为明显。图 3 为较粗一组泥沙的相对含沙量沿垂线分布的沿程变化情况。图 3 表明，尽管进口断面按 Rouse 公式给出了含沙量沿垂线的分布，但由于该断面实际处于不平衡输沙状态，这种分布并不是稳定的。在距进口 200m 处，泥沙的分布调整到一种不平衡输沙状态，随着泥沙的沿程淤积，水流输沙向平衡方向 发展 ，垂线平均含沙量趋向于水流挟沙力，而含沙量沿垂线分布向平衡时的分布状态 (Rouse 公式 ) 发展。由于这种发展是趋向于稳定状态，因此愈接近下游，分布愈靠近 Rouse 公式。计算结果表明，本文提出的计算 方法 是合理可行的。 图 1 含沙量垂线分布的沿程变化 (Z=0.01) 图 2 含沙量垂线分布的沿程变化 (Z=0.4) Changes of vertical distributions of sediment concentrations(Z=0.01) Chang