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

一维泥沙淤积计算水库一维泥沙淤积计算课程设计目录一、目的与要求1二、基本原理11、基本方程12、方程离散13、公式补充1三、计算步骤1四、计算框图1五、计算结果11、历年输沙量特征值12、各年淤积总量13、各年水位库容关系14、水面线的变化15、深泓变化16、坝前断面变化1六、结果分析11、剖面形态分析12、库容损失合理性分析1七、计算程序1一、 目的与要求通过课程设计，初步掌握一维数学模型建立数学模型的基本过程和计算方法，具备一定的解决实际问题的能力。以水流、泥沙方程为基础，构建恒定流条件下的河道一维水沙数学模型，并编制出完整的计算程序，并以某个水库为实例，进行水库泥沙淤积计算。水流条件：恒定非均匀流。泥沙条件：包括悬移质，推移质的均匀沙模型，推移质计算模式为饱和输沙，悬移质计算模式为不饱和输沙，水流泥沙方程采用非耦合解。二、 基本原理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.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,-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(*,