河流模拟课程设计--new.docx

1、键 入 公 司 地 址 武汉大学水利水电学院河流模拟课程设计班 级：港航二班姓 名： 学 号：指导老师：吴卫民河流模拟课程设计1目录1. 课程设计的目的 .22.课程设计的内容 .23 数学建模的原理 .23.1 基本方程 23.2 方程的离散 33.3 补充公式及参数选取 43.4 计算步骤 44 计算成果分析 .54.1 历年流量输沙量特征值 54.2 历年淤积总量 64.3 历年库容水位关系 64.4 深泓变化 74.5 坝前断面淤积变化 74.6 计算结果分析 85 源程序解释说明 .9河流模拟课程设计21. 课程设计的目的通过课程设计，初步掌握一维数学模型建立数学模型的的基本过程和计

2、算方法，具备一定的解决实际问题的能力。2.课程设计的内容以水流、泥沙方程为基础，构建恒定流条件下的河道一维水沙数学模型，并编制出完整的计算程序，并以某个水库为实例，进行水库泥沙淤积计算。水流条件：恒定非均匀流。泥沙条件：包括悬移质，推移质的均匀沙模型，推移质计算模式为饱和输沙，悬移质计算模式为不饱和输沙，水流泥沙方程采用非耦合解。3 数学建模的原理3.1 基本方程水流连续方程：(1)0xQtA水流运动方程(2)figAxhxt 02或 (3)3422RnQztQ泥沙连续方程（4）)(*SSxAt 河流模拟课程设计3河床变形方程（5）)(*0SxGtyb推移质平衡输沙方程G=G* （6）水流挟沙

3、力公式采用张瑞瑾公式，推移质输沙率公式采用 Mayer-Peter 公式：=()320.047()320.12512( )N 为曼宁糙率系数， 为河床平整情况下的沙粒曼宁糙率系数公式中的能坡 J 按均匀流曼宁公式近似计算（每个断面不同） 。3.2 方程的离散方程（1）在恒定流情况下有，离散为:Q=const0xQ方程（3）变形为 0342RAnxzgA或 2342Qx上式离散为 0)1(21 343421212 jjjjjjj RAxnzAQg方程（4）去掉时间项得到 )(*SqxS河流模拟课程设计4该方程的解析解为： qxxqqxSSSjjjjjj ep1ep1*1*由方程（4-5 ）可得

4、00tyBxQBGb对2号断面以下，上式可以离散为： 0)1(100 tyyxSxBGjjjbb为推移质输沙率，S为悬移质含沙量对于进口断面，推移质不考虑，悬移质采用单点离散方程（5）可离散为：01*0)(tSy3.3 补充公式及参数选取 mgRukS3*其中，k 取 0.124，m 取 1.05，恢复饱和系数 ，干密度 取 1.3，25.00均匀沙粒径为 d=0.041mm。3.4 计算步骤（1） 输入河床地形糙率等数据（2） 读入一个时段的水沙数据 （特别注意，不要一次性将数据全部读入）（3） 计算水面线，同时得到各断面的水力要素（4） 计算前要注意在坝前输入水位，各断面均应对流量赋值（5

5、） 计算水流挟沙力河流模拟课程设计5（6） 计算推移质输沙率（7） 计算各断面含沙量（8） 计算各断面冲淤厚度（9） 修改水各断面水下河床高程（10）重新进入（2）进行下一循环（11）计算10年河床变形，计算时段为一天，单位为秒（s）（12）淤积总量年输出一次，其余每两年输出一次计算结果4 计算成果分析4.1 历年流量输沙量特征值图表 1 历年流量输沙量特征值年份 悬移质输沙量(万 t) 推移质输沙量(万 t) 年径流量( 亿 m3)1984 1821.49 27.32235 48.45561985 3469.716 52.04574 61.48621986 11702.33 175.5349

6、5 86.63891987 3033.784 45.50676 50.39671988 1640.996 24.61494 29.58451989 2175.498 32.63247 31.7691990 6301.058 94.51587 65.84921991 4724.525 70.867875 56.60981992 1874.591 28.118865 32.81821993 2535.495 38.032425 37.9423河流模拟课程设计64.2 历年淤积总量图表 2 淤积总量随时间变化曲线4.3 历年库容水位关系图表 3 历年库容水位曲线0 0.5 1 1.5 2 2.5 3

7、2402452502552602652702750年2年4年6年8年10年库 容( 亿m3)水位(m)1982 1984 1986 1988 1990 1992 199400.511.522.5累计淤积总量( 万m3)年 份淤积量（m3）河流模拟课程设计74.4 深泓变化图表 4 历年深泓变化趋势图0 5 10 15 20 25 30 35 402002102202302402502602702800年2年4年6年8年10年距 坝 里 程(km)高程(m)4.5 坝前断面淤积变化图表 5 8#横断面历年淤积地形变化-100 -50 0 50 100 150 200 25023024025026

8、02702800年2年4年6年8年10年起 点 距(m)高程(m)河流模拟课程设计84.6 计算结果分析坝前呈三角洲淤积形态，选取坝前 CS-8#断面的历年冲淤变化作图，直观地显示了坝前逐年淤积的事实。由图表 2 也明显地看出淤积总量逐渐增加，随之库容逐渐减少，这可由图表 3 反映出来。另外，有了人为的限制条件，即允许河道冲刷已经淤积的泥沙，而不允许冲刷原始河床，故程序中事先设定了允许误差范围，由质量守恒关系计算冲淤的量的关系，入库泥沙与出库泥沙的差值理论上即为淤积在水库中的泥沙，但是考虑到计算过程中是以 267m 的正常蓄水位来计算的，尤其在库尾回水区可能会有高于 267m 出的淤积量没有计

9、入总淤积量，因此二者有一定的误差，但是只要控制在误差不大于 5%，即认为计算结果合理。河流模拟课程设计95 源程序解释说明PROGRAM MAIN PARAMETER (NPXT=16, MM=100, SFAI=0.85, FAI=0.50, IDAY=3653,IYEAR=10)PARAMETER (RS=26.5, RW=10.0, ROU0=1300,ROU=1000, ROUS=2650,DT=86400,&D=0.000041)DIMENSION X(MM,2,NPXT),DX(NPXT),DXA(NPXT),NPOINT(NPXT),DY(NPXT),&ALOW(NPXT),AL

10、OW0(NPXT)DIMENSION ROUGH(NPXT),A(NPXT),XW(NPXT),B(NPXT),H(NPXT),U(NPXT),&ZLEVEL(NPXT) DIMENSION Q(NPXT),S(NPXT),GB(NPXT),SX(NPXT)DIMENSION NDAY(IYEAR) DATA NDAY /366,731,1096,1461,1827,2192,2557,2922,3288,3653/ROUGH=0.035W=0.000909 !沉降速度,使用张瑞瑾泥沙沉降速度公式计算WRITE(*,*) 正在打开文件.OPEN(1,FILE=SECTION.TXT,STATU

11、S=OLD)OPEN(2,FILE=QS.TXT,STATUS=OLD)OPEN(3,FILE=淤积总量.TXT,STATUS=UNKNOWN)OPEN(4,FILE=纵剖面.TXT,STATUS=UNKNOWN)OPEN(5,FILE=坝前横断面.TXT,STATUS=UNKNOWN)OPEN(6,FILE=库容.TXT,STATUS=UNKNOWN)OPEN(7,FILE=水面线.TXT,STATUS=UNKNOWN)write(*,*) 已打开文件!输入部分河流模拟课程设计10WRITE(3,3) !输出“淤积总量“表头WRITE(4,4) !输出“纵剖面“表头WRITE(5,5) !输

12、出“坝前横断面“表头WRITE(6,6) !输出“库容“表头WRITE(7,7) !输出“水面线“表头3 FORMAT(2X,年,1X,悬移质输沙量 (万 t),1X,年径流量(万 t),1X,&总淤积沙量(万 m3)4 FORMAT(4X,年,5X,断面号 ,2X,距坝里程(km),3X, 深泓(m)5 FORMAT(4X,年,3X,断面号 ,6X,节点号,4X,起点距 m,3X,高程)6 FORMAT(4X,年,6X,水位 (m),7X,库容(亿 m3)7 FORMAT(4X,年,3X,距坝里程 (km),4X,水位(m)13 FORMAT(1X,I4,4X,F10.3,4X,F12.4,

13、4X,F10.3)14 FORMAT(1X,I4,6X,I4,4X,F8.2,4X,F8.2)15 FORMAT(1X,I4,6X,I4,6X,I4,4X,F8.2,4X,F8.2)16 FORMAT(1X,I4,4X,F10.2,4X,F10.3)17 FORMAT(1X,I4,4X,F8.2,4X,F8.2)DO J=1,NPXTREAD(1,*) DXA(J)read(1,*) NPOINT(J)rough(j)=0.035DO I=1,NPOINT(J)READ(1,*) X(I,1,J),X(I,2,J) !J 代表断面，I 代表节点IF(I.EQ.1).OR.(I.EQ.NPOIN

14、T(J) THENX(I,2,J)=X(I,2,J)+50ENDIFENDDOALOW0(J)=X(1,2,J)河流模拟课程设计11DO I=2,NPOINT(J)IF(ALOW0(J)=X(I,2,J) ALOW0(J)=X(I,2,J) !求初始的深泓（各断面） ,录入 ALOW0，后面求得的深泓不能小于该值ENDDO ENDDOALOW=ALOW0 !将深泓录入 ALOWWRITE(*,*) 已读入河床地形资料和糙率，并且得到初始深泓WRITE(*,*) 正在进行河床各断面的整理.DO J=1,(NPXT-1)DX(J)=DXA(J)-DXA(J+1) !计算断面间距ENDDODO J=

15、1,NPXTWRITE(4,14) 0,J,DXA(J)/1000.0,ALOW0(J) ENDDODO J=1,NPXTDO I=1,NPOINT(J)WRITE(5,15) 0,J,I,X(I,1,J),X(I,2,J) ENDDOENDDODO SKSW=273,240,-2 CALL ZV(NPOINT,NPXT,X,DX,SKSW,V,MM) WRITE (6,16) 0,SKSW,V ENDDOSKSW=267CALL ZV(NPOINT,NPXT,X,DX,SKSW,V,MM) WRITE (6,16) 0,SKSW,V 河流模拟课程设计12SXYZS=0.0 !累计悬移质输沙量

16、初始化SXYZ=0.0 !累计悬移质冲淤量初始化STYZ=0.0 !累计推移质淤积总量初始化AVQ=0.0 !累计年均流量初始化!-循环计算 -WRITE(*,*) 正在进行循环计算与输出 DO JDAY=1,IDAY !JDAY 用于计天数用IF(MOD(JDAY,200)=0) WRITE(*,*) 正在进行第,JDAY,次循环READ(2,*) Q(1),S(1)IF(Q(1)=0) GOTO 100Q=Q(1)ZLEVEL(NPXT)=267CALL LEVEL(X,ROUGH,NPXT,ZLEVEL,DX,Q,NPOINT,B,A,XW,& NPXT,MM,FAI,n) !计算各断面

17、水力要素和水面线DO J=1,NPXTU(J)=Q(J)/A(J)SX(J)=0.124*(U(J)*3.0)/(9.8*(A(J)/XW(J)*W)*1.05)H(J)=A(J)/B(J)GB(J)=ZMPM(ROUGH(J),RW,H(J),Q(J),A(J),RS,D,ROU,ROUS)*B(J)END DOGB(1)=Q(1)*S(1)*0.015!计算每个断面的 DYDY(1)=0.25*W*(S(1)-SX(1)*DT/ROU0 IF(ALOW(1)+DY(1)0.05)THENwrite(*,*) 质量不守恒,J,SUM1,SUM2PAUSEENDIFDO J=1,NPXTDO

18、I=1,NPOINT(J)IF(X(I,2,J)ZLEVEL(J) X(I,2,J)=ZLEVEL(J)!地形修改后，不能超过水面线ENDIFENDDOENDDODO J=1,NPXTALOW(J)=X(1,2,J)DO I=2,NPOINT(J) !修改地形后重新整理深泓数据（DY）IF(ALOW(J)=X(I,2,J) ALOW(J)=X(I,2,J)ENDDOENDDOXYZ1=Q(1)*S(1)!来沙量（入口，悬移质，kg/s）,此处 Q(1)指当天进口处来沙情况，非断面的意思。GB(1)=XYZ1*0.015!推移质来沙量按悬移质的 1.5%计算（kg/s ）SXYZS=SXYZS+

19、XYZ1*DT/10000000 !为累计悬移质输沙量 (万 t)SXYZ=SXYZ+(XYZ1-Q(1)*S(NPXT)*DT/(ROU0)/100000000 !为累计悬移质冲淤量 (万 m3),S(NPXT)表示出口处的泥沙含量STYZ=STYZ+(GB(1)-GB(NPXT)*DT/(ROU0)/100000000!为累计推移质淤积总量 (万 m3)AVQ=AVQ+Q(1)*DT/10000!为累计年均流量 (万 t)SYJ=SXYZ+STYZ!为累计年淤积总量 (万 m3)!输出部分河流模拟课程设计15100 DO JYEAR=1,IYEAR !JYEAR,计年数用IF(JDAY=N

20、DAY(JYEAR).AND.MOD(JYEAR,2)=0) THENDO J=1,NPXTWRITE(4,14) JYEAR,J,DXA(J)/1000.0,ALOW(J) !输出各个断面的深泓（每两年） ，到坝距离用 KM 计ENDDODO j=1,npxtDO i=1,npoint(j)WRITE(5,15) JYEAR,j,I,X(I,1,j),X(I,2,j)END DOEND DO !输出坝前横断面的坝距、高程（每两年）DO SKSW=273,240,-2 !SKSW 代表水库水位CALL ZV(NPOINT,NPXT,X,DX,SKSW,V,MM) WRITE (6,16) JY

21、EAR,SKSW,V !将每两年的库容曲线记录到 库容(亿 m3)（静库容）ENDDOSKSW=267CALL ZV(NPOINT,NPXT,X,DX,SKSW,V,MM) WRITE (6,16) JYEAR,SKSW,V !将每两年的 267m 水位下的库容曲线记录到后面DO J=1,NPXTWRITE(7,17) JYEAR,DXA(J)/1000.0,ZLEVEL(J)!输出水面线（每两年） （取最后一年）ENDDOENDIFIF(JDAY=NDAY(JYEAR) THEN !输出每年的淤积量河流模拟课程设计16WRITE(3,13) JYEAR,SXYZS,AVQ,SYJ !SYJ

22、为淤积总量SXYZS=0.0 !累计悬移质输沙量清零SXYZ=0.0 !累计悬移质冲淤量清零STYZ=0.0 !累计推移质淤积总量清零AVQ=0.0 !累计年均流量清零ENDIFENDDOENDDOWRITE(*,*) 循环完成，共进行了,JDAY-1,次循环WRITE(*,*) 质量守恒已验证!-循环结束 -WRITE(*,*) 文件已输出!-关闭文件- CLOSE(1)CLOSE(2)CLOSE(3)CLOSE(4)CLOSE(5)CLOSE(6)CLOSE(7)WRITE(*,*) 文件已关闭END!-子程序部分-FUNCTION ZMPM(ROUGH,R,H,Q,A,RS,D,ROU,

23、ROUS)!用 MAYER-PETER 公式求单宽推移质输沙率 gbROUT=ROU/1000ROUST=ROUS/1000河流模拟课程设计17!ROUT 为水的密度 (t/m3),ROUST 为沙粒密度(t/m3)RJ=(ROUGH*2.0)*(Q*2.0)/(A*2.0)*(H*(4.0/3) !曼宁公式求 JROUGHNN=(D*(1.0/6)/26.0 !河床平整情况下的沙粒曼宁糙率系数 NA0=(ROUGHNN/ROUGH)*1.50)*R*H*RJ-0.047*(RS-R)*DFENMU=(ROUST-ROUT)/ROUST)*9.8*(ROUT*0.5)*0.125IF(A0)0

24、) THENZMPM=0ELSEZMPM=1000*(A0)*1.5)/FENMU !单位转换为 kg/(s*m)ENDIFEND !-库容子程序 -SUBROUTINE ZV(NPOINT,NPXT,XX,DX,ZLEVEL,V,MM)! 计算给定条件下的相应某种水位下的库容( 亿 m3)(静库容)DIMENSION NPOINT(NPXT),XX(MM,2,NPXT),SA(NPXT),DX(NPXT)DO J=1,NPXTCALL AREA(NPOINT(J),XX(1,1,J),XX(1,2,J),B,A,XW,ZLEVEL,0)SA(J)=AENDDOV=0DO J=1,NPXT-1

25、V=V+(SA(J)+SA(J+1)*DX(J)/2/100000000ENDDORETURNEND!-过水断面面积子程序 -SUBROUTINE AREA(NPOINT,X,Z,B,A,XW,ZLEVEL,NDISP)DIMENSION X(NPOINT),Z(NPOINT)河流模拟课程设计18IF(NDISP.LE.-1)WRITE(*,*)INTO AREAIF(NPOINT.LE.2)THENWRITE(*,*)IN AREA THE INPUT DATA ARE FALSE N=,NpointSTOPENDIF IF(Z(1).LT.ZLEVEL.OR.Z(NPOINT).LT.ZL

26、EVEL)THENWRITE(*,*)IN AREA THE WATER LEVEL IS TOO HIGH OR THE HEIGHT & OF THE SECTION AT EDGE IS TOO LOWSTOPENDIF A=0.B=0.XW=0. DO 10 I=1,NPOINT-1ZMIN=AMIN1(Z(I),Z(I+1)IF(ZMIN.LT.ZLEVEL)THENZMAX=AMAX1(Z(I),Z(I+1)IF(ZMAX.LT.ZLEVEL)THENDB=X(I+1)-X(I)DH=ZLEVEL-0.5*(Z(I)+Z(I+1)DX=(Z(I)-Z(I+1)*2.+DB*DB)*

27、0.5ELSEDB=(ZLEVEL-ZMIN)/(ZMAX-ZMIN)*(X(I+1)-X(I)DH=0.5*(ZLEVEL-ZMIN)DX=(2*DH)*2.0+DB*DB)*0.5ENDIFIF(DB.LT.0.0)THENWRITE(*,*)THE DISTANCE FOR NODE ,I,AND,I+1,ARE FALSE河流模拟课程设计19WRITE(*,*)IN AREA,I,I+1,X(I),X(I+1) STOPENDIFDA=DB*DHB=B+DBA=A+DAXW=XW+DXENDIF10 CONTINUERETURNEND!-二分法求根子程序-SUBROUTINE BISE

28、C(ZMIN,ZMAX,FMIN,ZLEVEL,NPOINT,X,Z,B,A, XW,Q,& NDISP,ZLELO,BLOWER,ALOWER,XWLOWER,QLOWER,DX,ROUGH,FAI)DIMENSION X(NPOINT),Z(NPOINT)IF(NDISP.EQ.-1)WRITE(*,*)INTO BISECERR=0.00110 DDZ=ZMAX-ZMINZLEVEL=0.5*(ZMIN+ZMAX)CALL AREA(NPOINT,X,Z,B,A,XW,ZLEVEL,NDISP)F=FLEVEL(zlelo,zlevel,dx,qlower,q,rough,blower,

29、b,alower,a,fai,ndisp) IF(DDZ.GT.ERR)THENFF=F*FMINIF(FF.GT.0)THENZMIN=ZLEVELELSEZMAX=ZLEVELENDIF河流模拟课程设计20GOTO 10ENDIFRETURNEND!-水面线子程序 -SUBROUTINE level(x,rough,npxt,zlevel,dx,q,npoint,b,a,xw,nn,mm,fai,ndisp)! 恒定非均匀缓流水面线计算 DIMENSION X(MM,2,NN),ROUGH(NPXT),DX(NPXT),ZLEVEL(NPXT),& Q(NPXT),NPOINT(NPXT)

30、,B(NPXT),A(NPXT),XW(NPXT)IF(NDISP.EQ.-1)WRITE(*,*)INTO LEVELNC=1000 ! 最大循环次数DZ=3.DZ1=0.1317DZ2=0.179CALL AREA(NPOINT(NPXT),X(1,1,NPXT),X(1,2,NPXT),B(NPXT),& A(NPXT),XW(NPXT),ZLEVEL(NPXT),NDISP)DO 10 IP=NPXT-1,1,-1NR=0ZMIN=ZLEVEL(IP+1)+DZ30 CALL AREA(NPOINT(IP),X(1,1,IP),X(1,2,IP),B(IP), & A(IP),XW(N

31、PXT),ZMIN,NDISP)IF(A(IP).LE.0)THENZMIN=ZMIN+0.173GOTO 30ENDIFFMIN=FLEVEL(ZLEVEL(IP+1),ZMIN,DX(IP),Q(IP+1),Q(IP),& ROUGH(IP),B(IP+1),B(IP),A(IP+1),A(IP),FAI,NDISP)IF(FMIN.GT.0)THEN河流模拟课程设计21FR=(Q(IP)/A(IP)*2.0*B(IP)/(9.8*A(IP)IF(FR.LT.1)THENZMAX=ZMIN+DZ20 CALL AREA(NPOINT(IP),X(1,1,IP),X(1,2,IP),B(IP

32、),& A(IP),XW(IP),ZMAX,NDISP)FMAX=FLEVEL(ZLEVEL(IP+1),ZMAX,DX(IP),Q(IP+1),Q(IP), & ROUGH(IP),B(IP+1),B(IP),A(IP+1),A(IP),FAI,NDISP)IF(FMAX.GT.0)THENZMAX=ZMAX+DZGOTO 20ENDIFELSEZMIN=ZMIN+DZ2NR=NR+1IF(NR.GT.NC)THENZMIN=ZMIN-DZ1WRITE(*,*)THE LOOP IS DEATH,PAUSEWRITE(*,*),NR,IP,ZMIN,FMIN,NR,IP,ZMIN,FMINR

33、EAD(*,*)ENDIFGOTO 30ENDIFELSEZMIN=ZMIN-DZ1GOTO 30ENDIF CALL BISEC(ZMIN,ZMAX,FMIN,ZLEVEL(IP),NPOINT(IP),X(1,1,IP),& X(1,2,IP),B(IP),A(IP),XW(IP),Q(IP),NDISP,ZLEVEL(IP+1),B(IP+1),河流模拟课程设计22& A(IP+1),XW(IP+1),Q(IP+1),DX(IP),ROUGH(IP),FAI)10 CONTINUERETURNENDFUNCTION FLEVEL(ZLELO,ZLEVEL,DX,QLOWER,Q,ROUGH,BLOWER,B, & ALOWER,A,FAI,NDISP)IF(NDISP.EQ.-1)WRITE(*,*)INTO FLEVELHLOWER=ALOWER/BLOWERH=A/BAA=ZLELO-ZLEVELB1=FAI*(QLOWER/BLOWER)*2.0/HLOWER*(10/3)B2=(1-FAI)*(Q/B)*2.0/H*(10/3)BB=DX*ROUGH*2.0*(B1+B2)C1=(QLOWER/ALOWER)*2.0C2=(Q/A)*2.0CC=(C1-C2)/(2*9.8)FLEVEL=AA+BB+CCRETURNEND