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

一维泥沙淤积计算水库一维泥沙淤积计算课程设计武汉大学水利水电学院2013-3-15目录一、目的与要求1二、基本原理11、基本方程12、方程离散13、公式补充2三、计算步骤3四、计算框图4五、计算结果51、历年输沙量特征值52、各年淤积总量53、各年水位库容关系64、水面线的变化75、深泓变化86、坝前断面变化9六、结果分析121、剖面形态分析122、库容损失合理性分析12七、计算程序13一、 目的与要求通过课程设计，初步掌握一维数学模型建立数学模型的基本过程和计算方法，具备一定的解决实际问题的能力。以水流、泥沙方程为基础，构建恒定流条件下的河道一维水沙数学模型，并编制出完整的计算程序，并以某个水库为实例，进行水库泥沙淤积计算。水流条件：恒定非均匀流。泥沙条件：包括悬移质，推移质的均匀沙模型，推移质计算模式为饱和输沙，悬移质计算模式为不饱和输沙，水流泥沙方程采用非耦合解。二、 基本原理1、 基本方程水流连续方程： 水流运动方程 或 泥沙连续方程 河床变形方程 推移质平衡输沙方程G=G* 水流挟沙力公式采用张瑞瑾公式，推移质输沙率公式采用Mayer-_Peter公式，MAYER-PETER公式中的能坡J按均匀流曼宁公式近似计算（每个断面不同）。2、 方程离散方程 在恒定流情况下有，离散为:Q=const 方程 变形为 或 上式离散为方程（4）去掉时间项得到该方程的解析解为：由方程（4-5）可得 对2 号断面以下，上式可以离散为：对于进口断面，推移质不考虑，悬移质采用单点离散方程（5）可离散为：3、 公式补充K取 0.124，m取1.05，干密度取1.3恢复饱和系数 均匀沙粒径为d=0.041mm（悬移质），d=2 mm（推移质）三、 计算步骤1、输入河床地形糙率等数据求得断面面积与水位的关系（AZ），进而求得断面平均流速 ，水力学半径 2、读入一个时段的水沙数据（特别注意，不要一次性将数据全部读入）读入第一时段（Q,S）值3、计算水面线，同时得到各断面的水力要素求得各个断面的河宽、断面面积、水深、平均流速等值计算前要注意在坝前输入水位，各断面均应对流量赋值4、计算悬移质水流挟沙力K取 0.124，m取1.05。5、计算推移质输沙率（采用mayer-peter公式）6、计算各断面含沙量公式7、 计算各断面冲淤厚度进口断面 8、 修改水各断面水下河床高程9、重新进入（2）进行下一循环10、计算10年河床变形，计算时段为一天，单位为秒（s）11、淤积总量年输出一次，其余每两年输出一次计算结果四、 计算框图开始读入地形资料和糙率输出初始库容、深泓等i=i+1读入一个QYQ=0计算水流子程序计算S*和Gbj=j+1计算S(j)和dy(j)YDy(i)=0,重新计算S(j) 判断库尾冲刷与否YJnpnxt修改河床地形判断年份输出结果YIx(m,2,i)thenalow0(i)=x(m,2,i)endifenddo enddoend!*泥沙沉降速度的计算,采用张瑞瑾公式*FUNCTION FW(T,gama,gamas,d,ndisp) IF(NDISP.EQ.-1)WRITE(*,*)INTO FWcall VISCOS(T,CMU)A1=CMU*13.95/dA2=1.09*9.8*d*(gamas-gama)/gamaFW=(A1*2.0+A2)*0.5-A1 RETURN ENDSUBROUTINE VISCOS(T,CMU) X=1.775E-06 A=1+0.0337*T+0.*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,-1NR=0 zmin=zlevel(ip+1)+dz30 call area(npoint(ip),x(1,1,ip),x(1,2,ip),b(ip),a(ip),xw(npxt),zmin,ndisp) if(A(ip)0)thenfr=(q(ip)/a(ip)*2.0*b(ip)/(9.8*a(ip) if(fr0)thenzmax=zmax+dzgoto 20endifelsezmin=zmin+dz2NR=NR+1IF(NRNC)THENzmin=zmin-dz1write(*,*)the loop is death,pauseWRITE(*,*)nr,ip,zmin,fmin,nr,ip,zmin,fminread(*,*)endifgoto 30endif elsezmin=zmin-dz1goto 30 endif 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),a(ip+1),q(ip+1),dx(ip),rough(ip),failev) enddo return end FUNCTION flevel(zlelo,zlevel,dx,qlower,q,rough,blower,b,alower,a,failev,ndisp) IF(NDISP.EQ.-1)WRITE(*,*)INTO FLEVEL HLOWER=ALOWER/BLOWER H=A/B AA=ZLELO-ZLEVEL B1=failev*(QLOWER/BLOWER)*2.0/HLOWER*(10/3.0) B2=(1-failev)*(Q/B)*2.0/H*(10/3.0) BB=DX*ROUGH*2.0*(B1+B2) C1=(QLOWER/ALOWER)*2.0 C2=(Q/A)*2.0 CC=(C1-C2)/(2*9.8) FLEVEL=AA+BB+CC RETURN END!*计算断面要素*SUBROUTINE AREA(NPOINT,X,Z,B,A,XW,ZLEVEL,NDISP) DIMENSION X(NPOINT),Z(NPOINT) IF(NDISP=-1)WRITE(*,*)INTO AREA IF(NPOINT=2)THENWRITE(*,*)IN AREA THE INPUT DATA ARE FALSE N=,Npoint STOP ENDIF IF(Z(1)ZLEVEL.OR.Z(NPOINT)ZLEVEL)THEN WRITE(*,*)IN AREA THE WATER LEVEL IS TOO HIGH OR THE HEIGHT OF THE SECTION AT EDGE IS TOO LOW STOP ENDIF A=0.0 B=0.0 XW=0.0 DO I=1,NPOINT-1ZMIN=AMIN1(Z(I),Z(I+1) IF(ZMINZLEVEL)THENZMAX=AMAX1(Z(I),Z(I+1) IF(ZMAXZLEVEL)THENDB=X(I+1)-X(I)DH=ZLEVEL-0.5*(Z(I)+Z(I+1)DX=(Z(I)-Z(I+1)*2.0+DB*DB)*0.5 ELSEDB=(ZLEVEL-ZMIN)/(ZMAX-ZMIN)*(X(I+1)-X(I)DH=0.5*(ZLEVEL-ZMIN)DX=(2*DH)*2.0+DB*DB)*0.5 ENDIF IF(DBerr)thenff=f*fmin if(ff0)thenzmin=zlevel else zmax=zlevel endif ddz=zmax-zmin goto 10 endif return end!悬移质挟沙力和推移质输沙率计算，悬移质采用张瑞谨挟沙力公式!推移质采用mayer-peter公式SUBROUTINE SGB(nn,a,b,Q,Sx,Gb,rough,w,ndisp)DIMENSION A(nn),B(nn),Q(nn),Sx(nn),Gb(nn),rough(nn)IF(ndisp=-1)WRITE(*,*)INTO SGBDO I=1,nnH=A(I)/B(I)U=Q(I)/A(I)VA1=U*3.0VA2=9.8*H*WVA3=(VA1/VA2)*1.05Sx(I)=VA3*0.124 !悬移质挟沙力!推移质输沙率，断面输沙率Gb(I)=FGb(rough(I),10.0,H,Q(I),A(I),26.5,0.002,1.0,2.65,0)*B(i) ENDDOENDFUNCTION FGb(rn,r,H,Q,A,rs,d,rou,rous,ndisp)IF(ndisp=-1)WRITE(*,*)INTO FGbA1=(rn*2.0)*(Q*2.0)A2=(A*2.0)*(H*(4/3.0)RJ=A1/A2rnn=(d*(1/6.0)/26.0AA=(rnn/rn)*1.5)*rs*H*RJBB=0.047*(rs-r)*dCC=(rous-rou)/rous)*9.8*(rou*0.5)*0.125IF(AA-BB)x(L,2,J) alow(J)=x(L,2,J)ENDDOIF(alow(J)+dy(J)ABS(0.05*TOL1)THENWRITE(*,*)MASS NOT EQU