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

1、剖面二维非恒定悬移质泥沙扩散方程的数值方法摘要：通过讨论剖面二维非恒定泥沙扩散方程的数值方法，建立了一种用于求解含沙量分布沿程变化的差分格式(Z-C格式)并通过一个具体的数值例子说明了计算的方法步骤。关键词：扩散方程差分格式精度稳定性1引言数学模拟方法正在成为研究河流泥沙问题的重要手段。目前，一维数学模型发展较成熟，已广泛应用于模拟长河段的长期变形，但它只能给出河段平均冲淤深度的沿程变化，如需了解短河段的河床变形细节，则要采用二维以至三维数学模型。不论是一维数学模型还是平面二维维数学模型，都不能反映含沙量沿垂线的分布状况，并忽略了含沙量沿垂线分布对垂线平均含沙量变化过程的影响。要解决这类问题，

2、必须建立剖面二维数学模型。这种模型主要通过解剖面二维泥沙扩散方程来研究悬移质泥沙沿水深的分布及含沙量的变化过程，对水电站进口和其它引水工程的引水口高程的确定都能提供较好的数值模拟。数值方法相对于解析方法在求解偏微分方程上有着明显的优势，即简单灵活、计算方便快捷，但要寻找一种精度高、稳定性好、计算方便的差分格式也并非易事。本文拟在前人研究的基础上着重讨论剖面二维泥沙扩散方程的数值解问题，希望能提供一种精度高、稳定性好、计算方便的数值解。2基本方程剖面二维泥沙扩散方程的形式为(1)式中x,y为水流方向和铅直方向的维轴；u,v分别为沿水流方向和铅直方向的时均流速；sx,sy分别为水流方向和铅直方向的

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

4、传递系数，为修正值。由勃兰特尔掺长理论可得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网格的剖分为建立差分方程，首先必须剖分网格。我们取时间步长t=，X方向的空间步长x=h1，Y方向的空间步长为y=h2,这样形成如下网格Dh=(xj,yl,tn)xj=jh1,yl=lh2,tn=n其中j=0,.,N;l=0,.,M;n0;3.2构造差分格式通过对流方程和扩散方程的差分格式的构造，我们可以得到对流扩散方程

5、的差分格式。由于隐式格式稳定性好，考虑Crank-Nicholson型隐式格式。为此，引入差分算子记号2ySnj,l=Snj,l+1-2Snj,l+Snj,l-1LxSnj,l=Snj+1,l-Snj-1,lLySnj,l=Snj,l+1-Snj,l-1为了看得更清楚，暂且取h1=h2=h.对式(6)离散，则C-N格式为Sn+1j,l-Snj,l/+u/4hLx(Snj,l+Sn+1j,l)=s/2h22y(Sn+1j,l+Snj,l)+W/4hLy(Sn+1j,l+Snj,l)(7)C-N格式的精度是二阶的，绝对稳定。但对于二维问题，由(7)导出的方程组，其系数矩阵不是三对角矩阵，不能用追赶

6、法求解。因此，考虑构造交替方向的隐式格式(命名为Z-C格式)(Sn+1/2j,l-Snj,l)/(/2)+u/2hLxSnj,l=s/h22ySn+1/2j,l+W/2hLySn+1/2j,l(8)(Sn+1j,l-Sn+1/2j,l)/2+u/2hLxSn+1j,l=s/h22ySn+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,

7、l+Snj,l)=s/2h22ySn+1/2j,l+W/4hLxSn+1/2j,l(10)将(8)和(9)两式相减，可得Sn+1/2j,l=(Sn+1j,l+Snj,l)/2+s/4h22y(Snj,l-Sn+1j,l)+W/8hLy(Snj,l-Sn+1j,l)(11)把式(11)代入(10)，变形整理，可得(Sn+1j,l-Snj,l)(1-us/8h3Lx2y-2uW/16h2LxLy)=(Snj,l+Sn+1j,l)(us/2h22y+W/4hLy-u/4hLx)(12)设S(x,y,t)是(12)的精确解，并假定S(x,y,t)关于t三次连续可微，关于x,y四次连续可微，那么利用Ta

8、ylor级数展开可得S(xj,yl,tn+1)-S(xj,yl,tn)/(1-2s/8h3Lx2y-2uW/16h2LxLy)-S(xj,yl,tn+1)-S(xj,yl,tn)/(s/2h2)2y+W/4hLy-u/4hLx)=O(2+h2)(13)由此可见，Z-C格式具有二阶精度。3.4稳定性分析现在，我们来讨论Z-C格式的稳定性。为此，把式(12)变形整理得(1-s/2h22y-W/4hLy)(1+u/4hLx)Sn+1j,l=(1+s/2h2y2+W/4hLy)(1-u/4hLx)Snj,l(14)由式(14)可得出过渡因子为G(,k)=(1-2s/h2sin2k2h/2+iW/2hs

9、ink2h)(1-iu/2hsink1h)/(1+2s/h2sin2k2h/2-iW/2hsink2h)(1+iu/2hsink1h)(15)令a=2s/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)|21,所以Z-C格式是绝对稳定的。4数值计算4.1边界条件我们考虑初边值问题。(1)初始条件用Rouse公式给出含沙量沿垂线

10、分布S(x,y,0)=Sa5(a/h-a)z(h-y/y)z(17)式中z=/ku*为悬浮指标，Sa为近底含沙量，h为水深，一般取a=0.010.05h。(2)水面条件(y=h)(18)(3)底部边界条件(y=a)(19)式中Sa*为近底挟沙力，即输沙平衡时的近底售沙量Sa.(4)近底含沙量计算近底含沙量在求解泥沙扩散方程时具有边界条件性质，它选取的正确与否，意味着所给边界条件是否正确。实际工程中一般缺乏实测资料，近底含沙量不易测定。这里，我们利用水流挟沙力和含沙量沿垂线分布公式来反求近底含沙量。已知断面平均挟沙力为S*=k(u3/gh)m(20)假定Sa*=S*(21)输沙平衡时，含沙量沿垂

11、线分布用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-W/4h2)Sn+1/2j,l+1=u/4h1LxSnj,l-Snj,l(24)-u/4h1Snj-1,l+Snj,l+u/4h1Snj+1,l=s/2h222ySn+1/2j,l+W/4h2LySn+1/2j,l+Sn+1/2j,l(25)在一个时间层(第n层

12、)内，计算分两步进行：第一步，对式(24)用追赶法求第n+1/2层的过渡值。令C1=-u/4h1,C2=1,C3=u/4h1,E1=s/2h22,E2=W/4h2D1=s/2h22-W/4h2,D2=-s/h22-1,D3=s/2h22+W/4h2l=1时，D1Sn+1/2j,0+D2Sn+1/2j,1+D3Sn+1/2j,2=C3LxSnj,1-Snj,1l=2时，D1Sn+1/2j,0+D2Sn+1/2j,1+D3Sn+1/2j,2=C3LxSnj,1-Snj,1l=M时,D1Sn+1/2j,M-1+D2Sn+1/2j,M+D3Sn+1/2j,M+1=C3LxSnj,M-Snj,M令Hl=

13、-Snj,l+C3LxSnj,l(1lM)H1=-Snj,1+C3LxSnj,1-D1Sn+1/2j,0HM=-Snj,M+C3LxSnj,M-D3Sn+1/2j,M+1其中Sn+1/2j,0和Sn+1/2j,M+1由边界条件给出，则用矩阵形式表示为（26）第二步，再对(25)式用追赶法求第n+1层的值令Fj=E12ySn+1/2j,l+E2LySn+1/2j,l+Sn+1/2j,l(1jN)F1=E12ySn+1/21,l+E2LySn+1/21,l+Sn+1/21,l-C1Sn+10,lFN=E12ySn+1/2N,l+E2LySn+1/2N,l+Sn+1/2N,l-C3Sn+1N+1,l

14、其中Sn+10,l和Sn+1N+1，l由边界条件给出，则同理可得矩阵方程(27)这样，按此步骤一层层地计算。4.3数值模拟合理性分析受所掌握的实测资料的*，目前尚无法对本文提出的算法与含沙量沿垂线分布的实测值进行对比。我们用库里阿雷克沉沙池的实测资料作了垂线平均值沿程变化的比较。该沉沙池的主要数据为：池深h=1.53m;平均流速u=0.12m/s;泥沙沉速=0.0176cm/s;悬浮指标Z=0.01。计算时取卡门常数=0.4，a=0.05h。表1给出了计算值和实测值，结果表明，计算值和实测值比较符合。表1断面平均含沙量验证(单位：kg/m3)Verificationsofcross-secti

15、onalaveragesedimentconcentrations距离(m)020040060080010001200计算值3.0122.5732.0711.6391.2090.8490.539实测值3.0122.6501.8101.5201.1410.8220.561为了进一步分析含沙量垂线分布计算结果的合理性，我们对另一组较粗的泥沙(=0.616cm/s,Z=0.4)进行了对比计算。图1,图2分别为两组沙的计算结果。从图中可以看出，计算结果符合含沙量沿垂线分布的一般规律，粗沙分布不均匀，细沙分布较均匀；近底浓度相对较大，水面浓度相对较小，不存在Rouse公式中水面含沙量为0的缺陷；含沙量沿

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

17、angesofverticaldistributionsofsedimentconcentrations(Z=0.01)Changesofverticaldistributionsofsedimentconcentrations(Z=0.4)5结语本文建立了求解二维非恒定泥沙扩散方程的一种差分格式(Z-C格式)。这种格式具有如下特点：1.精度较高(具有二阶精度)。2.稳定性好(无条件稳定)。3.计算较方便(每一时段利用两次追赶法即可)。图3相对含沙量的垂线分布变化(Z=0.4)Changesofverticaldistributionsofrelativesedimentconcentrations