河流模拟课程设计—水库一维泥沙淤积计算

合集下载
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

河流模拟课程设计—水库一维泥沙淤积计算
水库一维泥沙淤积计算课程设计
武汉大学水利水电学院
2013-3-15
目录
一、目的与要求 (1)
二、基本原理 (1)
1、基本方程 (1)
2、方程离散 (1)
3、公式补充 (2)
三、计算步骤 (3)
四、计算框图 (4)
五、计算结果 (5)
1、历年输沙量特征值 (5)
2、各年淤积总量 (5)
3、各年水位库容关系 (6)
4、水面线的变化 (7)
5、深泓变化 (8)
6、坝前断面变化 (9)
六、结果分析 (12)
1、剖面形态分析 (12)
2、库容损失合理性分析 (12)
七、计算程序 (13)
一、目的与要求
通过课程设计,初步掌握一维数学模型建立数学模型的基本过程和计算方法,具备一定的解决实际问题的能力。

以水流、泥沙方程为基础,构建恒定流条件下的河道一维水沙数学模型,并编制出完整的计算程序,并以某个水库为实例,进行水库泥沙淤积计算。

水流条件:恒定非均匀流。

泥沙条件:包括悬移质,推移质的均匀沙模型,推移质计算模式为饱和输沙,悬移质计算模式为不饱和输沙,水流泥沙方程采用非耦合解。

二、基本原理
1、基本方程
水流连续方程:
0=??+??x
Q
t A ①
水流运动方程
()f i i gA x h gA A
Q x t Q -=??+
+
02

或03
4222=+??+???? ????+??R
A n Q g x z gA A Q x t Q ③
泥沙连续方程
()())(*S S QS x
SA t --=??+??
αω④ 河床变形方程
)(*00
S S x
G t y b
--=??+??αωρ⑤ 推移质平衡输沙方程
G=G * ⑥
水流挟沙力公式采用张瑞瑾公式,推移质输沙率公式采用Mayer-_Peter 公式,MAYER-PETER 公式中的能坡J 按均匀流曼宁公式近似
计算(每个断面不同)。

2、方程离散
方程①在恒定流情况下有
0=??x
Q
,离散为:Q=const 方程③变形为
034222=+??+
R A n Q x z A Q
x gA Q 或023422222=+??+
R A n Q x z gA Q x 上式离散为
0)1((21343422
1212121222121=ψ-+ψ?+-+???? ??-++++++j
j j j j j j j j j j j R A Q R A Q xn z z A Q A Q g 方程(4)去掉时间项得到
)(*S S q
x S --=??αω 该方程的解析解为:
()()???? ?
--?-+
--+
=+++q x x q q x S S S S S S j j
j
j
j j αωαωαωexp 1exp 1
***1*1 由方程(4-5)可得
()()00'
0=??+??+??t
y B x QS x BG b ρ 对2 号断面以下,上式可以离散为:
()()()()0)1(1010'0
=
ψ+ψ-+?-+?-++t
y B y B x
QS QS x
BG BG j j j j b b ρ
对于进口断面,推移质不考虑,悬移质采用单点离散方程(5)可离散为: '
01*10)(ραωt
S S y ?-=
3、公式补充
m
gR u k S ?
=ω3*
K 取 0.124,m 取1.05,干密度'
0ρ取1.3 恢复饱和系数25.0=α
均匀沙粒径为d=0.041mm (悬移质),d=2 mm (推移质)
1、输入河床地形糙率等数据
求得断面面积与水位的关系(A ~Z ),进而求得断面平均流速
A Q
u =
,水力学半径χ
A R = 2、读入一个时段的水沙数据(特别注意,不要一次性将数据全部读入)读入第一时段(Q,S )值
3、计算水面线,同时得到各断面的水力要素
求得各个断面的河宽、断面面积、水深、平均流速等值计算前要注意在坝前输入水位,各断面均应对流量赋值4、计算悬移质水流挟沙力
m
gR u k S
=ω3*
K 取 0.124,m 取1.05。

5、计算推移质输沙率(采用mayer-peter 公式)
g
d hJ n n g s
s s b )(125.0)(047.0)(2/12
/32/3ρρ
ρργγγ---'=
6、计算各断面含沙量公式
()()???? ?
--?-+???? ???--+
=+++q x x q q x S S S S S S j j
j
j
j j αωαωαωexp 1exp 1
***1*1
7、计算各断面冲淤厚度进口断面'
1*10)(ραωt S S y ?-=?
8、修改水各断面水下河床高程 Y H H ?+=0
9、重新进入(2)进行下一循环
10、计算10年河床变形,计算时段为一天,单位为秒(s ) 11、淤积总量年输出一次,其余每两年输出一次计算结果
1、历年输沙量特征值
计算各年悬疑质输沙量、推移质输沙量、年均流量见表1。

(其中
推移质输沙量是按照悬疑质输沙量的1.5%计算的)。

表1龙开口坝址历年流量输沙量特征值
2、各年淤积总量
表2 历年淤积总量的变化
图1历年淤积总量的变化3、各年水位库容关系
表3各年水位库容关系
00.5
11.522.50
2
4
6
8
10
淤积量(
亿立方米)年份
图2各年水位库容关系曲线4、水面线的变化
表4水面线变化
240
2452502552602652702750 0.5 1
1.5
2
2.53水位(m )
库容(108 m3)
图3 水面线变化曲线5、深泓变化
表5逐年深泓变化
265 2662672682692702712722732740 100002000030000
40000
水位(m )距坝里程(m )
图4 深泓变化曲线
6、坝前断面变化
该水库为三角形淤积,最靠近坝前断面淤积并不会特别严重,为使断面变化情况表现的更加直观,这里采用IP-17断面作为典型断面并分析其淤积情况。

200
220
240
2602800
10000
20000
3000040000
高程
(m )
距坝里程(m )
表6 坝前断面变化
图5坝前断面变化曲线200 2202402602803003200 50
100
150
200
250300高程(m )
起点距(m )
0年2年4年6年8年10年
六、结果分析
1、剖面形态分析
A、纵断面淤积情况分析:
由深泓变化曲线图可以看出,随着时间的增加,河道深泓点高程成增高趋势。

该水库运行10年后,在距坝约10000米处,深泓点高程增加40多米,计算结果表明,该水库泥沙以淤积为主,致使水库库容逐年损失,且情况较为严重。

B、坝前断面变化曲线分析:
由坝前断面变化曲线图可以看出,随着时间的增加,坝前断面各点的总体高程也成增高趋势。

该水库运行10后,该断面河床抬高20多米,表明在该处河道逐年以淤积为主,且淤积情况较为严重。

2、库容损失合理性分析
由水位库容曲线可以看出,水库在运行10年后,经计算可知,在273m水位条件下,库容的损失率为77.89%,在241m的水位条件下,库容的损失率高达95.92%,计算结果表明,经过十年的运行,水库库容基本成全部淤满状态,所以在10年后将不能再继续使用。

结果表明,该水库的调节能力较差。

上述计算结果与现实不太相符,其原因是在计算该水库库容的泥沙损失时,计算冲往下游的泥沙量不合理。

在给定条件下,该水库并没有考虑水位的变化及坝前冲沙情况。

该水库的计算条件为正常蓄水位保持不变,在高水位运行条件下,必然会导致水库的严重淤积,而实际情况是,在洪水来临时,坝前水位要降至防洪限制水位,以增大防洪库容,该条件下,水流的速度及水流挟沙力较大,会造成库区泥沙处于冲刷状态,以增加水库库容。

七、计算程序
program main
parameter(nn=16,mm=60,ndisp=0,npxt=16,ND=3653,NY=
10,T=15)
paramete(rs=26.5,r=10.0,d=0.000041,failev=0.5,faidy=0.5,dt =86400)
dimension
x(mm,2,nn),rough(nn),npoint(nn),B(nn),A(nn),Gb(nn),dy(nn),Nda y(NY) dimension dx(nn),Xw(nn),dxa(nn),alow(nn),alow0(nn),Q(nn),Sx(nn),S(nn),Zle vel(nn)
data(Nday(m),m=1,10)/366,731,1096,1461,1827,2192,2557,2922 ,3288,3653/
!**************打开输出文件****************
open(10,file='河床地形.txt',status='old')
open(11,file='流量沙量.txt',status='old')
open(12,file='逐年累计淤积量.txt',status='unknown')
open(13,file='逐年深泓变化.txt',status='unknown')
open(14,file='逐年水面线变化.txt',status='unknown')
open(15,file='逐年坝前断面淤积变化.txt',status='unknown')
open(16,file='逐年水位库容关系变化.txt',status='unknown')
!**************读入地形数据****************
call Qbed(mm,npxt,x,dxa,alow0,npoint)
!**************计算沉降速度****************
W=FW(T,r,rs,d,0)
write(*,*)W
!**************计算断面间距****************
do i=1,npxt-1
dx(i)=dxa(i)-dxa(i+1)
enddo
!****计算初始水面线、深泓、库容水位关系****
dWQST=0.0
dWQS=0.0
dWGb=0.0
AVQ=0.0
do I=1,NPXT
rough(I)=0.035
enddo
Q(1)=295.
do I=2,NPXT
Q(I)=Q(I-1)
enddo
Zlevel(npxt)=267.0
CALL
level(x,rough,nn,zlevel,dx,q,npoint,b,a,xw,nn,mm,failev,NDISP) !初始水面线
!********输出到文件**********
CALL
FILE(12,13,14,15,16,0,x,npoint,Zlevel,alow,dxa,dx,mm,nn,dW QST,AVQ,dWQS,dWGb,ndisp) !*计算各个时段水库的淤积情况* DO 40,I=1,ND
WRITE(*,*)'读入天数:',I
!***读入流量、沙量数据***
READ(11,*)Q(1),S(1)
IF(Q(1)==0)THEN
GOTO 41
ENDIF
DO J=2,NPXT
Q(J)=Q(1)
ENDDO
!*****************************************
CALL
level(x,rough,nn,zlevel,dx,Q,npoint,b,a,xw,nn,mm,failev,0)
!**计算各断面的水流挟沙力和推移质输沙率**
CALL SGB(nn,a,b,Q,Sx,Gb,rough,w,ndisp)
!******计算各断面的含沙量和冲淤厚度******
CALL
DEPO(nn,mm,w,x,S,Sx,Gb,dt,B,Q,dx,dy,alow,alow0,npoint,faidy,nd isp)
!**************检验质量守恒**************
CALL CHECK(dy,B,S,Gb,dt,npxt,Zlevel,x,npoint,Q,dx,ndisp)
41 CONTINUE
!******计算悬移质、推移质年均流量********
WQS1=Q(1)*S(1)
Gb(1)=WQS1*0.015
WQSN=Q(nn)*S(nn)
dWQST=dWQST+WQS1*dt/10000000 !悬移质累计输沙量
dWQS=dWQS+(WQS1-WQSN)*dt/1300.0/100000000 !悬移质累计冲淤量
dWGb=dWGb+(Gb(1)-Gb(nn))*dt/1300.0/100000000 !推移质累计淤积量
AVQ=AVQ+Q(1)*dt/10000 !累计流量
!****输出相关文件***
DO m=1,NY
IF(I==Nday(m))then
CALL
FILE(12,13,14,15,16,m,x,npoint,Zlevel,alow,dxa,dx,mm,nn,dW QST,AVQ,dWQS,dWGb,ndisp) ENDIF
ENDDO
40 CONTINUE
close(10)
close(11)
close(12)
close(13)
close(14)
close(15)
close(16)
end
!**************读入地形数据****************** subroutine Qbed(mm,nn,x,dxa,alow0,npoint) dimension x(mm,2,nn),dxa(nn),alow0(nn),npoint(nn) do i=1,nn
read(10,*)
read(10,*) dxa(i),npoint(i)
read(10,*)
do j=1,npoint(i)
read(10,*) x(j,1,i),x(j,2,i)
enddo
!*******计算深泓*******
alow0(i)=x(1,2,i)
do m=2,npoint(i)
if(alow0(i)>x(m,2,i))then
alow0(i)=x(m,2,i)
endif
enddo
enddo
end
!*泥沙沉降速度的计算,采用张瑞瑾公式* FUNCTION FW(T,gama,gamas,d,ndisp)
IF(NDISP.EQ.-1)WRITE(*,*)'INTO FW'
call VISCOS(T,CMU)
A1=CMU*13.95/d
A2=1.09*9.8*d*(gamas-gama)/gama
FW=(A1**2.0+A2)**0.5-A1
RETURN
END
SUBROUTINE VISCOS(T,CMU)
X=1.775E-06
A=1+0.0337*T+0.000221*T*T
CMU=X/A
RETURN
END
!*************水面线及各断面水力要素计算函数*************** subroutine
level(x,rough,npxt,zlevel,dx,q,npoint,b,a,xw,nn,mm,failev,NDISP) dimension x(mm,2,nn),rough(npxt),dx(npxt),zlevel(npxt)
dimension q(npxt),npoint(npxt),b(npxt),a(npxt),xw(npxt)
IF(NDISP==-1)WRITE(*,*)'INTO LEVEL'
nc=1000
dz=0.5
dz1=0.1317
dz2=2.079
call
area(npoint(npxt),x(1,1,npxt),x(1,2,npxt),b(npxt),a(npxt),xw(npxt), zlevel(npxt),ndisp) do ip=npxt-1,1,-1
NR=0
zmin=zlevel(ip+1)+dz
30 call area(npoint(ip),x(1,1,ip),x(1,2,ip),b(ip),a(ip),xw(npxt),zmin,ndisp) if(A(ip)<=0)then
zmin=zmin+0.173
goto 30
endif
fmin=flevel(zlevel(ip+1),zmin,dx(ip),q(ip+1),q(ip),rough(ip),b
(ip+1),b(ip),a(ip+1),a(ip),failev,ndisp) if(fmin>0)then fr=(q(ip)/a(ip))**2.0*b(ip)/(9.8*a(ip))
if(fr<1)then。

相关文档
最新文档