剖面二维非恒定悬移质泥沙扩散方程的数值方法(张耀新 吴卫民)

热度109票 浏览105次 【共0条评论】【我要评论 时间:1999年4月01日 13:46

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

-Su3@+Y8di+U0

张耀新 吴卫民水利论文lI7ZVOU

(武汉水利电力大学泥沙室)水利论文w)W%I9qq:y#~

摘 要通过讨论剖面二维非恒定泥沙扩散方程的数值方法,建立了一种用于求解含沙量分布沿程变化的差分格式(Z-C格式)并通过一个具体的数值例子说明了计算的方法步骤。水利论文Q3mf jSY Lg

关键词扩散方程 差分格式 精度 稳定性

"yj2hj ?3eZ0

1 引言

W+L1cb'a^ LpOmQ0

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

  泥沙扩散方程实际上是一个变系数的二阶线性偏微分方程,这样的方程在各种复杂边界条件下求解是极为困难的。求扩散方程的解析解在数学上存在着难以克服的困难,往往只能通过对方程的简化,才能得到一些简单边界条件下的解析解,在这方面,A.A.Kalinske、野满隆治、W.E.Dobbins、俞维强、张启舜、韦直林[2]等都做了有益的尝试;求扩散方程的数值解曾经因为缺乏高效率的计算工具而难以实现,直到60年代后,随着计算机的广泛应用,在各种复杂边界条件下求扩散方程的数值解不但成为可能,而且得到迅速的发展,在这方面,曹志先、崔侠[4]等做了大量工作,取得了很多成果。

\7Y n2Ti0

  数值方法相对于解析方法在求解偏微分方程上有着明显的优势,即简单灵活、计算方便快捷,但要寻找一种精度高、稳定性好、计算方便的差分格式也并非易事。本文拟在前人研究的基础上着重讨论剖面二维泥沙扩散方程的数值解问题,希望能提供一种精度高、稳定性好、计算方便的数值解。

"?4iJ;A$o4@3R4t~n(u0

2 基本方程

:zz ]}+M"r0
    剖面二维泥沙扩散方程的形式为

990207e1.gif (1861 字节)

Ni3S,p-H'~#|E7@@c0

(1)

8Ts#DI6bt l`5Q,Z0

式中 x,y为水流方向和铅直方向的维轴;u,v分别为沿水流方向和铅直方向的时均流速;εsxsy分别为水流方向和铅直方向的泥沙扩散系数;ω,S分别为泥沙静水沉速和含沙量。水利论文$c`Z_?&j6^/@

  对于式(1)的求解,研究者一般会对它进行不同程度的简化,为此引入以下假定中的一种或几种:水利论文A'X)Y'W3b1V:KU

  A.非恒定流可以概化为梯级式恒定流,即

990207e2.gif (1061 字节);水利论文6T3AyH8V&q Ew(A.n

  B.在一个时段内,认为泥沙运动可以概化为处于恒定状态,即水利论文H;~1u {~

990207e3.gif (978 字节);

m4AOPh`W*C j%n$K4o0

  C.在二维流动中,纵向扩散系数与方程其他项相比,可以忽略不计,即认为方程右端第一项可以忽略;水利论文hAp2cYK

  D.认为悬移质泥沙粒径均一,即ω=const;

0xr$v7[jd t[0

  E.认为水流为二维均匀流,即v=0。

y&G*f {b^-vs0i4E0

  为简单起见,我们讨论的范围限于水流条件为二维非恒定均匀流,悬移质泥沙粒径均匀,为此引入假定C、D、E。这时,泥沙扩散方程为水利论文v*cizfmA/h

990207e4.gif (1488 字节)

AGn7jx%Um&G'Y"`0

(2)水利论文tO WPb}z

目前,对εs的变化规律研究得不很充分,一般假定

εs=βεm水利论文+Q?*zuHgG4i

(3)

ByQwM0

其中εm为动量传递系数,β为修正值。由勃兰特尔掺长理论可得水利论文 l4zGA9wL

εm=κu*y/(1-y/h)

HCk#@9qr}/P0

(4)

nn1I;n5OA~ bmb0
式中κ为卡门常数,u*为摩阻流速。

对于u,我们取卡曼-勃兰特尔对数流速分布公式

;eH&G,?o K&^S/|0

(umax-u)/u*=(1/κ)ln(h/y)水利论文$rX@^bG0~4@?

(5)

q6z'[&Xin \$[g9Wc0

令W=ω+β(κu*/h)(1-2y/h),则式(2)可变形为

ipp5Az4`(b$JH0

990207e5.gif (1365 字节)水利论文v(X,V\$\nfH

(6)

@ frrD b$@\Q0
3 差分方程

3.1 网格的剖分

U}h%tgY0

 为建立差分方程,首先必须剖分网格。我们取时间步长Δt=τ,X方向的空间步长Δx=h1,Y方向的空间步长为Δy=h2,这样形成如下网格

:[ ^ kIDM0

Dh={(xj,yl,tn)|xj=jh1,yl=lh2,tn=nτ|}水利论文s5s k%x@G9e_ w

其中j=0,...,N; l=0,...,M;n≥0;水利论文0UZ ` s0I

3.2 构造差分格式

stU1e$G/D` B0

  通过对流方程和扩散方程的差分格式的构造,我们可以得到对流扩散方程的差分格式。由于隐式格式稳定性好,考虑Crank-Nicholson型隐式格式。为此,引入差分算子记号水利论文z'];@l5Q UcA

δ2ySnj,l=Snj,l+1-2Snj,l+Snj,l-1

3Ya*Gng4P'Av0

LxSnj,l=Snj+1,l-Snj-1,l

:FD@G1j'x6q:G0

LySnj,l=Snj,l+1-Snj,l-1水利论文NJ8n(aa$k

  为了看得更清楚,暂且取h1=h2=h.对式(6)离散,则C-N格式为

C9YGk7KIHE0

Sn+1j,l-Snj,l/τ+u/4hLx(Snj,l+Sn+1j,l)=εs/2h2δ2y(Sn+1j,l+Snj,l)+W/4hLy(Sn+1j,l+Snj,l)水利论文 P0e vP B{%p/z"WD

(7)水利论文rk"Q]5n$Rk`

C-N格式的精度是二阶的,绝对稳定。但对于二维问题,由(7)导出的方程组,其系数矩阵不是三对角矩阵,不能用追赶法求解。因此,考虑构造交替方向的隐式格式(命名为Z-C格式)水利论文t!WxX)J

(Sn+1/2j,l-Snj,l)/(τ/2)+u/2hLxSnj,ls/h2δ2ySn+1/2j,l+W/2hLySn+1/2j,l

Im9T1v~,\&L7d2OA]0

(8)水利论文`"O5M2D@_| Wj5~;`

(Sn+1j,l-Sn+1/2j,l)/τ/2+u/2hLxSn+1j,ls/h2δ2ySn+1/2j,l+W/2hLySn+1/2j,l水利论文(w&}.w*H @wuVh

(9)

E,i+N4o1z9`'@0

  可以看出,计算Sn+1j,l是由两步组成的,每一步仅是一个方向的隐式,故用两次追赶法即可。水利论文)|A*n-YEMF

3.3 精度分析

fV6l'KWO8l7R4J/z0

  现在,我们考虑Z-C格式的精度。先设法消去过渡值Sn+1/2j,l,为此,将(8)和(9)两式相加,可得水利论文\}tG |#G!qa8P0^

Sn+1j,l-Snj,l+τu/2hLx(Sn+1j,l+Snj,l)=τεs/2h2δ2ySn+1/2j,l+τW/4hLxSn+1/2j,l

;r Q#OW1g"k0

(10)水利论文-b#z9LDvaw

  将(8)和(9)两式相减,可得水利论文 ~Tt7l3Af SW

Sn+1/2j,l=(Sn+1j,l+Snj,l)/2+τεs/4h2δ2y(Snj,l-Sn+1j,l)+τW/8hLy(Snj,l-Sn+1j,l)水利论文~G4I&v g-}Mq

(11)水利论文#X)NR } {

把式(11)代入(10),变形整理,可得水利论文v_ WT+tez

(Sn+1j,l-Snj,l)(1-τuεs/8h3Lxδ2y2uW/16h2LxLy)=(Snj,l+Sn+1j,l)(τuεs/2h2δ2y+τW/4hLy-τu/4hLx)水利论文3V9OE0USJq+g&j0uOn

(12)水利论文A:C5] m7z VI*O$Ic

  设S(x,y,t)是(12)的精确解,并假定S(x,y,t)关于t三次连续可微,关于x,y四次连续可微,那么利用Taylor级数展开可得水利论文Qv|3dsE?*j

S(xj,yl,tn+1)-S(xj,yl,tn)/τ(1-τ2εs/8h3Lxδ2y2uW/16h2LxLy)-S(xj,yl,tn+1)-S(xj,yl,tn)/τ(τεs/2h22y+τW/4hLy-τu/4hLx)=O(τ2+h2)

s9K$as x mri(cgd0

(13)水利论文!Z)y,T G"?|;i

由此可见,Z-C格式具有二阶精度。水利论文7]-f}5x4w+G

3.4 稳定性分析

5[7e k4e5X[0

  现在,我们来讨论Z-C格式的稳定性。为此,把式(12)变形整理得水利论文C)d%pi*nN?]+~&@

(1-τεs/2h2δ2y-τW/4hLy)(1+τu/4hLx)Sn+1j,l=(1+τεs/2h2δy2+τW/4hLy)(1-τu/4hLx)Snj,l水利论文0m$kKL d8f%L5]

(14)水利论文$h7^nb'tn1V

由式(14)可得出过渡因子为水利论文R(w*[:T9o y$nu

G(τ,k)=[(1-2τεs/h2sin2k2h/2+iτW/2hsink2h)(1-iτu/2hsink1h)]/(1+2τεs/h2sin2k2h/2-iτW/2hsink2h)(1+iτu/2hsink1h)水利论文OV { lD^0y

(15)水利论文aM!v s+f

令a=2τεs/h2sin2k2h/2,b=τW/2hsink2h,c=τu/2hsink1h,则水利论文|C [3B)]

|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+b2]2=1-4a2+4a+4b2+4a3/[(1+a)2+b2]2

CS,Ou&\,m t9iCXV0

(16)

d3L5Ntk0G;J]6T!c0

显然,对于任意的τ,h,|G(τ,k)|2<1,所以Z-C格式是绝对稳定的。

E2Z D0`*o0

4 数值计算水利论文[v^Z@o

4.1 边界条件水利论文P&D^A{&mn4zo v8q2a

  我们考虑初边值问题。

VGL%@ V@^?0

(1)初始条件

iR|5U'd u0A{0
  用Rouse公式给出含沙量沿垂线分布

S(x,y,0)=Sa·5(a/h-a)z(h-y/y)z水利论文 O_M h]l%BS-vR2d4_

(17)

K9^z.kV:p0

式中z=ω/ku*为悬浮指标,Sa为近底含沙量,h为水深,一般取a=0.01~0.05h。水利论文H k ]*U7B;VK

(2)水面条件

990207e6.gif (1051 字节)(y=h)

3N"doH)q-]/tE$T:\;s0

(18)水利论文B0\W&H_8i+l4i

(3)底部边界条件

990207e7.gif (1087 字节)(y=a)水利论文2AIo9l7[~

(19)水利论文!wD)\%bEyR3kB

式中Sa*为近底挟沙力,即输沙平衡时的近底售沙量Sa.

-b/fI$r5\#O0X0

(4)近底含沙量计算近底含沙量在求解泥沙扩散方程时具有边界条件性质,它选取的正确与否,意味着所给边界条件是否正确。实际工程中一般缺乏实测资料,近底含沙量不易测定。这里,我们利用水流挟沙力和含沙量沿垂线分布公式来反求近底含沙量[4]水利论文(ce L.^ [9Z k

  已知断面平均挟沙力为

Uze)Ro1f%? t/S v0

S*=k(u3/ghω)m

,o'y}P6T+P"a0

(20)

|B#a+JB _?.S6C0
  假定

Sa*=αS*水利论文q G/x#Zl&XC {

(21)

v K\-dDl%d]:N0

  输沙平衡时,含沙量沿垂线分布用Rouse公式(17)表示,用(17)表达挟沙力的垂线分布 ,然后沿垂线积分得断面平均挟沙力为

_#r1Y"gz5G5\1wm D'E0

990207e8.gif (1567 字节)

$fM&`!w kf3hL0

(22)

h_4i#v d M\0
  将(22)与(21)比较,可得

990207e9.gif (1511 字节)

Fv1k9VN.P0

(23)水利论文vjC"` U4]

4.2 计算步骤水利论文I A'?1|O1yz

  为方便计算,将式(8)和(9)式变形整理,并对X,Y方向取不同的空间步长

VnW(gCJgw0

(τε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水利论文c^8gM{}`,H

(24)

n#V"D]K@tM0

-τu/4h1Snj-1,l+Snj,l+τu/4h1Snj+1,l=τεs/2h22ζ2ySn+1/2j,l+τW/4h2LySn+1/2j,l+Sn+1/2j,l

7q s!t#Alc_%b0

(25)水利论文){L(ae&f XM,];i

  在一个时间层(第n层)内,计算分两步进行:水利论文)Qx b0^"l i

  第一步,对式(24)用追赶法求第n+1/2层的过渡值。

j%VB+z#xM0

  令C1=-τu/4h1,C2=1,C3=τu/4h1,E1=τεs/2h22,E2=τW/4h2

HBX@+yCQ0

  D1=τεs/2h22-τW/4h2,D2=-τεs/h22-1,D3=τεs/2h22+τW/4h2

$nyK,uY]w0

  l=1时,D1Sn+1/2j,0+D2Sn+1/2j,1+D3Sn+1/2j,2=C3LxSnj,1-Snj,1水利论文(AY dd xP

  l=2时,D1Sn+1/2j,0+D2Sn+1/2j,1+D3Sn+1/2j,2=C3LxSnj,1-Snj,1

dBM!M} yo0

  ……………水利论文W|@3? L|!cc

  l=M时,D1Sn+1/2j,M-1+D2Sn+1/2j,M+D3Sn+1/2j,M+1=C3LxSnj,M-Snj,M水利论文+^fU IZ3KqT6m

   令Hl=-Snj,l+C3LxSnj,l(1<l<M)水利论文&t)p4J!S m5aB\#O

    H1=-Snj,1+C3LxSnj,1-D1Sn+1/2j,0

G}]|;@h{0

    HM=-Snj,M+C3LxSnj,M-D3Sn+1/2j,M+1水利论文Pn]XG(|!}*On

  其中Sn+1/2j,0和Sn+1/2j,M+1由边界条件给出,则用矩阵形式表示为

r1c/c%{ I3P,`0

990207e10.gif (2765 字节)

3v8@U7j#L,DN0

(26)水利论文c;c%c9P b!Z!u@;D

  第二步,再对(25)式用追赶法求第n+1层的值

\2?"C l^*ceL~0

  令Fj=E1δ2ySn+1/2j,l+E2LySn+1/2j,l+Sn+1/2j,l (1<j<N)

3q O&K-h5Sv:f |0

F1=E1δ2ySn+1/21,l+E2LySn+1/21,l+Sn+1/21,l-C1Sn+10,l水利论文(Kk%Q:qY"o^;{.e*{

FN=E1δ2ySn+1/2N,l+E2LySn+1/2N,l+Sn+1/2N,l-C3Sn+1N+1,l水利论文;b\.bbFN

  其中Sn+10,l和Sn+1N+1,l由边界条件给出,则同理可得矩阵方程水利论文2RL+n_-Ga(v*k

990207e11.gif (2631 字节)水利论文 bF0aX1vT

(27)

yn6F-_[J L5j0

  这样,按此步骤一层层地计算。水利论文q/Qx4J)L1D&x

4.3 数值模拟合理性分析水利论文 `}}_ Y&W2i;Y$V5I2_{

  受所掌握的实测资料的限制,目前尚无法对本文提出的算法与含沙量沿垂线分布的实测值进行对比。我们用库里·阿雷克沉沙池[5]的实测资料作了垂线平均值沿程变化的比较。该沉沙池的主要数据为:池深h=1.53m;平均流速u=0.12m/s;泥沙沉速ω=0.0176cm/s;悬浮指标Z=0.01。计算时取卡门常数κ=0.4,a=0.05h。表1给出了计算值和实测值,结果表明,计算值和实测值比较符合。

:gO5e oXY1?0

表1 断面平均含沙量验证(单位:kg/m3)

-~GHL }Wf0

Verifications of cross-sectional average sediment concentrations

7I4d'g _!@7gG E0

距离(m)

1l2\"coF*p0

0

x-lW5pj#Y6P)t@0

200水利论文#?@6o"q/W"RiFM

400

C,NGk7Q&P%@0

600水利论文*k4~@,M9CM1?OPU

800

!|Knd%z0

1000水利论文{Qk0i1bn,v9n/U

1200水利论文.d-P[%u;k4C^


计算值

Q} Q j%i:k-?$y0

3.012水利论文a&P ws I'A*s'y

2.573水利论文xoN0E[Ed

2.071水利论文.x1bOByi4nb@

1.639水利论文:|RQOw @ec

1.209

9g)sM*{(g3F9W!p0

0.849水利论文i$Sg[,t

0.539水利论文aX$SVP3gP

实测值水利论文T%s"l1V \:Y(Q#I

3.012水利论文*w"mOMh_f

2.650水利论文9L;C+o)v6wJ{ K:k~7[

1.810

H^]F4x/Hd)Q0

1.520水利论文Zt(a` dd2B.Mt

1.141

b(w9]rZ0j*x0

0.822水利论文1S3[G YhW D"K

0.561水利论文I`7g:Sw\"Xo


  为了进一步分析含沙量垂线分布计算结果的合理性,我们对另一组较粗的泥沙(ω=0.616cm/s,Z=0.4)进行了对比计算。图1\,图2分别为两组沙的计算结果。从图中可以看出,计算结果符合含沙量沿垂线分布的一般规律,粗沙分布不均匀,细沙分布较均匀;近底浓度相对较大,水面浓度相对较小,不存在Rouse公式中水面含沙量为0的缺陷;含沙量沿程衰减的特性较为明显。图3为较粗一组泥沙的相对含沙量沿垂线分布的沿程变化情况。图3表明,尽管进口断面按Rouse公式给出了含沙量沿垂线的分布,但由于该断面实际处于不平衡输沙状态,这种分布并不是稳定的。在距进口200m处,泥沙的分布调整到一种不平衡输沙状态,随着泥沙的沿程淤积,水流输沙向平衡方向发展,垂线平均含沙量趋向于水流挟沙力,而含沙量沿垂线分布向平衡时的分布状态(Rouse公式)发展。由于这种发展是趋向于稳定状态,因此愈接近下游,分布愈靠近Rouse公式。计算结果表明,本文提出的计算方法是合理可行的。

1HE0w%s'?E0

990207t1.gif (2561 bytes)水利论文;p4O;d6@&q`

990207t2.gif (2558 bytes)水利论文F TC)zn _

图1 含沙量垂线分布的沿程变化(Z=0.01)

4e.Su]1\]!B0

图2 含沙量垂线分布的沿程变化(Z=0.4)水利论文%WN9?N&v_

Changes of vertical distributions of sediment concentrations(Z=0.01)

D:o{8j2`1yR%t.n0

Changes of vertical distributions of sediment concentrations(Z=0.4)

[eTK'_ K;B)q0

5 结语水利论文%@\H5B_V6@%KFje

  本文建立了求解二维非恒定泥沙扩散方程的一种差分格式(Z-C格式)。这种格式具有如下特点:

&{I!^*c)po k0

  1.精度较高(具有二阶精度)。水利论文q A{m4NW;T}

  2.稳定性好(无条件稳定)。水利论文6UbhS"f7Z'L;J

  3.计算较方便(每一时段利用两次追赶法即可)。

&R s Y g#A`0

990207t3.gif (2330 bytes)水利论文m#R;gLANQN

图3 相对含沙量的垂线分布变化(Z=0.4)水利论文F9e?3J7CR

Changes of vertical distributions of relative sediment concentrations(Z=0.4)

&M2kd9@hMi#g0

参考文献水利论文Py p RZ&L%UV

1 陆金甫,关治。偏微分方程数值解法。清华大学出版社,1987年7月。水利论文.t7E5qYr?|$E

2 韦直林。二度恒定均匀流中泥沙淤积过程的研究。武汉水利电力学院研究生论文,1981年11月。水利论文` tjG!oGDQ

3 张瑞谨,谢鉴衡等。河流泥沙动力学。武汉水利电力学院,水利电力出版社,1989年6月。4 崔侠。剖面二维数学模型的研究。武汉水利电力学院研究生论文,1988年10月。

&Kd'_ fs9SZ0

5 韦直林。二度恒定均匀流中泥沙的淤积过程。武汉水利电力学院学报,1982年第4期。水利论文5C"n5k i`5CR+D

#|9Yv$OR06韦直林,傅小平。论泥沙连续方程的建立及相关的几个问题。武汉水利电力大学学报,199710月第30卷增刊。

*e3n ifmM zEn0

u/X f)F*ditd0 

YG:XO _0V3i(ABVb0
TAG: 泥沙 剖面 方程 吴卫民 张耀新
顶:13 踩:21
【已经有76人表态】
15票
极差
10票
很差
9票
较差
10票
稍差
11票
稍好
6票
较好
6票
很好
9票
极好
下一篇:空腔回流区水沙特性的计算分析(董耀华)
上一篇:金沙江流域的河流泥沙输移特性(潘久根)
查看全部回复【已有0位网友发表了看法】

广告投放

广告投放