河流模拟课程设计—水库一维泥沙-淤积计算(共29页).docx

精选优质文档-倾情为你奉上水库一维泥沙淤积计算课程设计武汉大学水利水电学院2013-3-15目录专心-专注-专业一、 目的与要求通过课程设计，初步掌握一维数学模型建立数学模型的基本过程和计算方法，具备一定的解决实际问题的能力。以水流、泥沙方程为基础，构建恒定流条件下的河道一维水沙数学模型，并编制出完整的计算程序，并以某个水库为实例，进行水库泥沙淤积计算。水流条件：恒定非均匀流。泥沙条件：包括悬移质，推移质的均匀沙模型，推移质计算模式为饱和输沙，悬移质计算模式为不饱和输沙，水流泥沙方程采用非耦合解。二、 基本原理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) 判断库尾冲刷与否YJ<npnxt修改河床地形判断年份输出结果YI<3653结束五、 计算结果1、历年输沙量特征值计算各年悬疑质输沙量、推移质输沙量、年均流量见表1。（其中推移质输沙量是按照悬疑质输沙量的1.5%计算的）。表1龙开口坝址历年流量输沙量特征值代表年年平均悬移质输沙量 （万t）年平均推移质输沙量 （万t）年平均流量 （万t）1984294144.1213301985412061.8017501986322148.3213601987482872.4217601988254338.1514401989642696.3919001990438565.7818001991597289.5818701992198229.7312901993645696.841900最大值645696.841900最小值198229.731290多年平均值428764.3116402、各年淤积总量表2 历年淤积总量的变化年悬移质输沙量(万吨)年均流量(万吨)总淤积量(亿立方米)000011821.490.1425291.2050.387316993.540.879420027.341.085521668.331.204623853.981.368730155.051.701834879.591.987936754.182.1171039289.662.243 图1 历年淤积总量的变化3、各年水位库容关系表3 各年水位库容关系水位（m）0年2年4年6年8年10年2732.5742.2071.5721.3130.7650.5692712.3852.0181.3841.1320.6160.4382692.2021.8351.2020.980.4930.3332672.0241.6611.0440.8590.3980.262651.8521.510.9280.7730.3460.232631.6881.3720.830.6980.3050.2042611.5321.2450.7450.6280.2690.1792591.3841.1260.6680.5630.2370.1562571.2471.0140.5940.5010.2070.1332551.1170.9070.5230.4420.1790.1122530.9950.8080.4570.3870.1540.0922510.880.7170.3990.3340.1310.0742490.7710.6320.3440.2870.110.0592470.670.5530.2930.2440.0890.0452450.5760.4780.2440.2040.0710.0342430.4920.4080.1990.1680.0530.0252410.4170.3470.1590.1370.0380.017图2 各年水位库容关系曲线4、 水面线的变化表4 水面线变化年距坝里程(m)水位(m)年距坝里程(m)水位(m)年距坝里程(m)水位(m)034920267.6434920267.88834920274.43033810267.16433810267.76833810274.17031350267.01431350267.53831350273.43029280267429280267.37829280272.93027000267427000267.21827000272.4024050267424050267.07824050271.5021846267421846267.02821846270.53019896267419896267819896269.58018270267418270267818270268.92015790267415790267815790268.1013170267413170267813170267.41010170267410170267810170267.05078202674782026787820267049502674495026784950267021002674210026782100267002674026780267年距坝里程(m)水位(m)年距坝里程(m)水位(m)年距坝里程(m)水位(m)234920267.14634920272.011034920273.99233810267.12633810271.751033810273.78231350267.06631350271.011031350273.49229280267.02629280270.371029280273.18227000267627000269.631027000272.79224050267624050268.641024050272.22221846267621846267.941021846271.71219896267619896267.391019896271.07218270267618270267.111018270270.482157902676157902671015790269.72131702676131702671013170268.942101702676101702671010170267.882782026767820267107820267.19249502676495026710495026722100267621002671021002672026760267100267图3 水面线变化曲线5、 深泓变化表5 逐年深泓变化年距坝里程(m)深泓(m)年距坝里程(m)深泓(m)年距坝里程(m)深泓(m)034920264.1434920264.77834920271.75033810261.7433810264.12833810270.46031350253.3431350262.29831350268.95029280249.4429280261.92829280268.6027000246.6427000263.23827000269.21024050241.3424050261.99824050267.71021846239.2421846261.5821846267.87019896237.9419896261.49819896268.14018270236.1418270252.83818270266.84015790230.1415790250.82815790264.62013170228.2413170240.94813170264.74010170224.4410170239.22810170262.5907820222.147820230.9887820253.7604950219.844950228.2984950237.2802100215.9242100219.9782100233.9900213.1740221.0480224.22年距坝里程(m)深泓(m)年距坝里程(m)深泓(m)年距坝里程(m)深泓(m)234920262.79634920269.961034920271.76233810262.41633810269.081033810270.93231350259.92631350267.251031350269.37229280259.53629280266.71029280269.62227000259.48627000267.211027000270.09224050250.62624050265.391024050268.98221846247.11621846265.191021846269.24219896242.17619896265.941019896269.7218270242.84618270264.711018270268.61215790233.11615790258.711015790267.03213170232.41613170249.521013170266.88210170226.49610170242.18101017026527820224.1167820234.07107820264.7924950220.6564950228.12104950245.6322100216.8662100222.71102100236.5420213.7860219.7100229.24图4 深泓变化曲线6、 坝前断面变化该水库为三角形淤积，最靠近坝前断面淤积并不会特别严重，为使断面变化情况表现的更加直观，这里采用IP-17断面作为典型断面并分析其淤积情况。表6 坝前断面变化年起点距(m)高程(m)年起点距(m)高程(m)年起点距(m)高程(m)00300.220300.240300.2010.5272.29210.5272.29410.5272.29015.8266.93215.8267415.8267.01018.2260.24218.2261.18418.2264.29018.9259.62218.9260.56418.9263.67019.4259.88219.4260.82419.4263.93024.8259.91224.8260.85424.8263.96026.4259.82226.4260.76426.4263.87044.2246.15244.2247.09444.2250.2061.2234.47261.2235.41461.2238.52078.6218.17278.6219.11478.6222.22091.5215.92291.5216.86491.5219.97097.3216.1297.3217.04497.3220.150106.8216.982106.8217.924106.8221.030121.6217.472121.6218.414121.6221.520132.4218.052132.4218.994132.4222.10145.6218.452145.6219.394145.6222.50173.8218.962173.8219.94173.8223.010191.7219.432191.7220.374191.7223.480203.6230.782203.6231.724203.6234.830216.5242.882216.5243.824216.5246.930228.7255.042228.7255.984228.7259.090241.1267.442241.1267.444241.1267.440253.7279.812253.7279.814253.7279.810265.3291.22265.3291.24265.3291.20279.1304.522279.1304.524279.1304.520293.1300.32293.1300.34293.1300.3年起点距(m)高程(m)年起点距(m)高程(m)年起点距(m)高程(m)60300.280300.2100300.2610.5272.29810.5272.291010.5272.29615.8267.01815.8267.011015.8267.01618.2267818.22671018.2267618.9266.41818.92671018.9267619.4266.67819.42671019.4267624.8266.7824.82671024.8267626.4266.61826.42671026.4267644.2252.94844.2264.221044.2267661.2241.26861.2252.541061.2255.09678.6224.96878.6236.241078.6238.79691.5222.71891.5233.991091.5236.54697.3222.89897.3234.171097.3236.726106.8223.778106.8235.0510106.8237.66121.6224.268121.6235.5410121.6238.096132.4224.848132.4236.1210132.4238.676145.6225.248145.6236.5210145.6239.076173.8225.758173.8237.0310173.8239.586191.7226.228191.7237.510191.7240.056203.6237.578203.6248.8510203.6251.46216.5249.678216.5260.9510216.5263.56228.7261.838228.726710228.72676241.1267.448241.1267.4410241.1267.446253.7279.818253.7279.8110253.7279.816265.3291.28265.3291.210265.3291.26279.1304.528279.1304.5210279.1304.526293.1300.38293.1300.310293.1300.3图5 坝前断面变化曲线六、 结果分析1、 剖面形态分析A、纵断面淤积情况分析：由深泓变化曲线图可以看出，随着时间的增加，河道深泓点高程成增高趋势。该水库运行10年后，在距坝约10000米处，深泓点高程增加40多米，计算结果表明，该水库泥沙以淤积为主，致使水库库容逐年损失，且情况较为严重。B、坝前断面变化曲线分析：由坝前断面变化曲线图可以看出，随着时间的增加，坝前断面各点的总体高程也成增高趋势。该水库运行10后，该断面河床抬高20多米，表明在该处河道逐年以淤积为主，且淤积情况较为严重。2、 库容损失合理性分析由水位库容曲线可以看出，水库在运行10年后，经计算可知，在273m水位条件下，库容的损失率为77.89%，在241m的水位条件下，库容的损失率高达95.92%，计算结果表明，经过十年的运行，水库库容基本成全部淤满状态，所以在10年后将不能再继续使用。结果表明，该水库的调节能力较差。上述计算结果与现实不太相符，其原因是在计算该水库库容的泥沙损失时，计算冲往下游的泥沙量不合理。在给定条件下，该水库并没有考虑水位的变化及坝前冲沙情况。该水库的计算条件为正常蓄水位保持不变，在高水位运行条件下，必然会导致水库的严重淤积，而实际情况是，在洪水来临时，坝前水位要降至防洪限制水位，以增大防洪库容，该条件下，水流的速度及水流挟沙力较大，会造成库区泥沙处于冲刷状态，以增加水库库容。七、 计算程序program mainparameter(nn=16,mm=60,ndisp=0,npxt=16,ND=3653,NY=10,T=15) paramete (rs=26.5,r=10.0,d=0.,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),Nday(NY)dimension dx(nn),Xw(nn),dxa(nn),alow(nn),alow0(nn),Q(nn),Sx(nn),S(nn),Zlevel(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-1dx(i)=dxa(i)-dxa(i+1)enddo!*计算初始水面线、深泓、库容水位关系*dWQST=0.0dWQS=0.0dWGb=0.0AVQ=0.0do I=1,NPXTrough(I)=0.035enddoQ(1)=295.do I=2,NPXTQ(I)=Q(I-1)enddoZlevel(npxt)=267.0CALL 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,dWQST,AVQ,dWQS,dWGb,ndisp)!*计算各个时段水库的淤积情况*DO 40,I=1,NDWRITE(*,*)'读入天数:',I!*读入流量、沙量数据*READ(11,*)Q(1),S(1)IF(Q(1)=0)THENGOTO 41ENDIFDO J=2,NPXTQ(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,ndisp) !*检验质量守恒*CALL CHECK(dy,B,S,Gb,dt,npxt,Zlevel,x,npoint,Q,dx,ndisp)41CONTINUE!*计算悬移质、推移质年均流量*WQS1=Q(1)*S(1)Gb(1)=WQS1*0.015WQSN=Q(nn)*S(nn)dWQST=dWQST+WQS1*dt/ !悬移质累计输沙量dWQS=dWQS+(WQS1-WQSN)*dt/1300.0/ !悬移质累计冲淤量dWGb=dWGb+(Gb(1)-Gb(nn)*dt/1300.0/ !推移质累计淤积量AVQ=AVQ+Q(1)*dt/10000 !累计流量!*输出相关文件*DO m=1,NYIF(I=Nday(m)thenCALL FILE(12,13,14,15,16,m,x,npoint,Zlevel,alow,dxa,dx,mm,nn,dWQST,AVQ,dWQS,dWGb,ndisp)ENDIFENDDO40 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,nnread(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)thenalow0(i)=x(m,2,i)endifenddo enddoend!*泥沙沉降速度的计算,采用张瑞瑾公式*FUNCTION FW(T,gama,gamas,d,ndisp) IF(NDISP.EQ.-1)WRITE(*,*)'INTO FW'call 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)thenzmin=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)thenfr=(q(ip)/a(ip)*2.0*b(ip)/(9.8*a(ip) if(fr<1)thenzmax=zmin+dz!*计算断面要素*20call area(npoint(ip),X(1,1,ip),x(1,2,ip),b(ip),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),failev,ndisp)if(fmax>0)thenzmax=zmax+dzgoto 20endifelsezmin=zmin+dz2NR=NR+1IF(NR>NC)THENzmin=zmin-dz1write(*,*)'the loop is death,pause'WRITE(*,*)'nr,ip,zmin,fmin',nr,ip,zmin,fminread(*,*)endifgoto 30endif elsezmin=zmin-dz1