1、水库一维泥沙 淤积计算课程设计目录一、目的与要求1二、基本原理11、基本方程12、方程离散13、公式补充2三、计算步骤2四、计算框图3五、计算结果41、历年输沙量特征值 42、各年淤积总量53、各年水位库容关系 64、水面线的变化 75、深泓变化86、坝前断面变化 9六、结果分析121、剖面形态分析 122、库容损失合理性分析 12七、 计算程序13、目的与要求具备一定并编制出完整通过课程设计，初步掌握一维数学模型建立数学模型的基本过程和计算方法, 的解决实际问题的能力。以水流、泥沙方程为基础，构建恒定流条件下的河道一维水沙数学模型, 的计算程序，并以某个水库为实例，进行水库泥沙淤积计算。水流
2、条件：恒定非均匀流。泥沙条件：包括悬移质，推移质的均匀沙模型，推移质计算模式为饱和输沙，悬移质计 算模式为不饱和输沙，水流泥沙方程采用非耦合解。二、基本原理1、基本方程水流连续方程:AQx0QQ2txAQQ2txAtSAcxyGb水流运动方程泥沙连续方程河床变形方程xgA0 thgAgA i。 ifx2 2Q n c g4-0AR3(S S*)(S S)推移质平衡输沙方程G=G水流挟沙力公式采用张瑞瑾公式，推移质输沙率公式采用Mayer-_Peter公式,MAYER-PETER公式中的能坡J按均匀流曼宁公式近似计算（每个断面不同）。2、方程离散方程在恒定流情况下有q 0，离散为：Q=const
3、x方程变形为_ QgA x AQ2n2a2r3Q2x 2gA2 2Q n24A R3上式离散为1 Lj22gQj2A2Zj 1 Zj2xn(AjQj4R/3j 1方程S（4）去掉时间项得到Sj i S j i SjS*jexpS*jSji-1 exp xBQx对2号断面以下，上式可以离散为:由方程（4-5 ）可得QS y00B 0BG bBGb QS QSxxBjyoj(10)B j 1 y0j 1t一（S S） q该方程的解析解为：对于进口断面，推移质不考虑，悬移质采用单点离散方程（5）可离散为：y。(S SJ t103、公式补充S3mk gRK取0.124 , m取 1.05，干密度0 取
4、 1.3恢复饱和系数0.25均匀沙粒径为d=0.041mm （悬移质），d=2 mm （推移质）二、计算步骤1、输入河床地形糙率等数据求得断面面积与水位的关系（A Z）,进而求得断面平均流速Q，水力学半径A2、读入一个时段的水沙数据（特别注意，不要一次性将数据全部读入）读入第一时段（Q,S）值3、计算水面线，同时得到各断面的水力要素求得各个断面的河宽、断面面积、水深、平均流速等值计算前要注意在坝前输入水位，各断面均应对流量赋值4、计算悬移质水流挟沙力gRK取 0.124，m取 1.05。5、计算推移质输沙率（采用 mayer-peter公式）3/2n 3/2（）hJ 0.047（ s ）d n
5、gb0.125 1/2（）gs6、计算各断面含沙量公式Sj1 S*j1 Sj Sj exp 2S*j Sj1expq x7、计算各断面冲淤厚度进口断面y。一（S1 SJ t08、修改水各断面水下河床高程H H。 丫9、重新进入（2）进行下一循环10、计算10年河床变形，计算时段为一天，单位为秒（ s）11、淤积总量年输出一次，其余每两年输出一次计算结果四、计算框图读入地形资料和糙率输出初始库容、深泓等YQ=0计算S和Gb计算水流子程序计算S(j)和dy(j)1r4i=i+11F读入一个判断库尾冲刷与否YJvnpnxt重新计算S(j)修改河床地形五、计算结果1、历年输沙量特征值1。（其中推移质输
6、沙量是按照计算各年悬疑质输沙量、推移质输沙量、年均流量见表悬疑质输沙量的1.5%计算的）。表1龙开口坝址历年流量输沙量特征值代表年年平均悬移质年平均推移质年平均流量输沙量（万t）输沙量（万t）（万t）1984294144.1213301985412061.8017501986322148.3213601987482872.4217601988254338.1514401989642696.3919001990438565.7818001991597289.5818701992198229.7312901993645696.841900取大值645696.841900最小值198229.7312
7、90多年平均值428764.3116402、各年淤积总量表2历年淤积总量的变化年悬移质输沙 量（万吨）年均流量（万吨）总淤积量（亿 立方米）000011821.494845560.1425291.20510994180.387316993.5419658070.879420027.3424697751.085521668.3327656171.204623853.9830875691.368730155.0537460621.701834879.5943121571.987936754.1846403382.1171039289.6650197622.243S 2 5 1 S 甲2 100 2
8、 4 6 8 10图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.6
9、680.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.13
10、70.0380.01727U2US26020245图2各年水位库容关系曲线4、水面线的变化表4水面线变化年距坝里程水位年距坝里程水位年距坝里程水位(m)(m)(m)(m)(m)(m)034920267.6434920267.88834920274.43033810267.16433810267.76833810274.17031350267.01431350267.53831350273.43029280267429280267.37829280272.93027000267427000267.21827000272.4024050267424050267.07824050271.502184
11、6267421846267.02821846270.53019896267419896267819896269.58018270267418270267818270268.92015790267415790267815790268.1013170267413170267813170267.41010170267410170267810170267.05078202674782026787820267049502674495026784950267021002674210026782100267002674026780267年距坝里程水位年距坝里程水位年距坝里程水位(m)(m)(m)(m)(m)
12、m)234920267.14634920272.011034920273.99233810267.12633810271.751033810273.78231350267.06631350271.011031350273.49229280267.02629280270.371029280273.18227000267627000269.631027000272.79224050267624050268.641024050272.22221846267621846267.941021846271.71219896267619896267.391019896271.072182702676182
13、70267.111018270270.482157902676157902671015790269.72131702676131702671013170268.942101702676101702671010170267.882782026767820267107820267.19249502676495026710495026722100267621002671021002672026760267100267图3水面线变化曲线5、深泓变化表5逐年深泓变化年距坝里程(m)深泓(m)年距坝里程(m)深泓(m)年距坝里程(m)深泓(m)034920264.1434920264.7783492027
14、1.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.14157902
15、50.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.416338
16、10269.081033810270.93231350259.92631350267.251031350269.37229280259.53629280266.71029280269.62227000259.48627000267.211027000270.09224050250.62624050265.391024050268.98221846247.11621846265.191021846269.24219896242.17619896265.941019896269.7218270242.84618270264.711018270268.61215790233.11615790258.
17、711015790267.03213170232.41613170249.521013170266.88210170226.49610170242.18101017026527820224.1167820234.07107820264.7924950220.6564950228.12104950245.6322100216.8662100222.71102100236.5420213.7860219.7100229.24280 r260220T0年潇泓T-丄年琛叛 亠4年擁涨召年济亂8毎儼泓迪年探泓1Q年水厨线200 11010QD0200005000040000图4深泓变化曲线6、坝前断面变
18、化该水库为三角形淤积，最靠近坝前断面淤积并不会特别严重，为使断面变化情况表现的更加直观，这里采用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
19、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.5201
20、32.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.
21、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.9267
22、619.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.2681
23、21.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.44102
24、41.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坝前断面
25、变化曲线分析：由坝前断面变化曲线图可以看出， 随着时间的增加，坝前断面各点的总体高程也成增高 趋势。该水库运行10后，该断面河床抬高 20多米，表明在该处河道逐年以淤积为主，且淤 积情况较为严重。2、库容损失合理性分析由水位库容曲线可以看出，水库在运行10年后，经计算可知，在273m水位条件下，库容的损失率为77.89%，在241m的水位条件下，库容的损失率高达95.92%，计算结果表明，经过十年的运行，水库库容基本成全部淤满状态，所以在10年后将不能再继续使用。结果表明，该水库的调节能力较差。上述计算结果与现实不太相符,其原因是在计算该水库库容的泥沙损失时,计算冲往下游的泥沙量不合理。在给定
26、条件下，该水库并没有考虑水位的变化及坝前冲沙情况。该水库的计算条件为正常蓄水位保持不变,在高水位运行条件下， 必然会导致水库的严重淤积， 而实际情况是，在洪水来临时，坝前水位要降至防洪限制水位，以增大防洪库容，该条件下,水流的速度及水流挟沙力较大，会造成库区泥沙处于冲刷状态，以增加水库库容。七、计算程序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.000041,failev=0.5,faidy=0.5,dt=86400)dime nsion
27、x(mm,2 ,nn ),rough( nn), npoint(nn ),B( nn ),A (nn ),Gb( nn ),dy (nn ),Nday(NY) dime nsion dx( nn ),Xw( nn ),dxa( nn ),alow( nn ),alowO( 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 /I*打开输出文件*open(10,file=河床地形.txt,status=old)open(11,fi
28、le=流量沙量.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(m m,n pxt,x,dxa,alow0 ,npoint)算沉降速*W=FW(T,r,rs,d,0)write
29、)W算断面间住巨 *do i=1, npxt-1dx(i)=dxa(i)-dxa(i+1)en ddo!*计算初始水面线、深泓、库容水位关系*dWQST=0.0dWQS=0.0dWGb=0.0AVQ=0.0do l=1,NPXT rough(l)=0.035 en ddoQ(1)=295.do I=2,NPXTQ(I)=Q(I-1)en ddoZlevel( npxt)=267.0CALL level(x,rough,nn,zlevel,dx,q,npoint,b,a,xw,nn,mm,failev,NDISP) !初始水面线输出到文件 *CALLFILE(12,13,14,15,16
30、0,x ,n poi nt,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)ENDDOI*CALL level(x,rough, nn ,zlevel,dx,Q, npoi nt,b,a,xw, nn ,mm,failev,0)!*计算各断面的水流挟沙力和推移质输沙率*CALL SGB( nn
31、 ,a,b,Q,Sx,Gb,rough,w, ndisp)!*计算各断面的含沙量和冲淤厚度*CALL DEPO( nn ,mm,w,x,S,Sx,Gb,dt,B,Q,dx,dy,alow,alowO, npoi nt,faidy, ndisp)I *验质量守丿恒 *CALL CHECK(dy,B,S,Gb,dt, npxt,Zlevel,x, npoi nt,Q,dx, ndisp)41CONTINUE!* 计算悬移质、推移质年均流量*WQS 仁Q(1)*SGb(1)=WQS1*0.015WQSN=Q( nn 广S( nn)dWQST=dWQST+WQS1*dt/10000000 !悬移质累计
32、输沙量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,NYIF(l=Nday(m)thenCALLFILE(12,13,14,15,16,m,x ,n poi nt,Zlevel,alow,dxa,dx,mm, nn ,dWQST,AVQ,dWQS,dWGb ,n disp)ENDIFENDDO40 CONTINUEclose(10)close(11)c
33、lose(12)close(13)close(14)close(15)close(16) endI *读入地形数据*subrouti ne Qbed(m m,nn, x,dxa,alowO ,npoint)dime nsion x(mm,2, nn ),dxa( nn ),alowO( nn), npoi nt( nn)do i=1, nnread(10,*)read(10,*) dxa(i), npoi nt(i)read(10,*)do j=1, npoin t(i)read(10,*) x(j,1,i),x(j,2,i)en ddo计算深泓*alow0(i)=x(1,2,i)do m=2
34、 npoin t(i)if(alow0(i)x(m,2,i)thenalow0(i)=x(m,2,i)en difen ddoen ddoend!*泥沙沉降速度的计算，采用张瑞瑾公式*FUNCTION FW(T,gama,gamas,d, ndisp)IF(NDIS 圧Q.-1)WRITE(*,*)INT0 FWcall VISCOS(T,CMU)A仁CMU*13.95/dA2=1.09*9.8*d*(gamas-gama)/gamaFW=(A1*2.0+A2)*0.5-A1RETURNENDSUBROUTINE VISCOS(T,CMU)X=1.775E-06A=1+0.0337*T+0.
35、000221*T*TCMU=X/ARETURNEND水面线及各断面水力要素计算函数*subroutine level(x,rough,npxt,zlevel,dx,q,npoint,b,a,xw,nn,mm,failev,NDISP) dime nsion x(mm,2, nn ),rough( npxt),dx( npxt),zlevel( npxt) dime nsion q(n pxt), npoint(n pxt),b( npxt),a( npxt),xw (n pxt) IF(NDISP=-1)WRITE(*,*)INTO LEVEL nc=1000dz=0.5dz1= 0.1317dz2=2.079call area( npoint(n pxt),x(1,1, npxt),x(1,2, npxt),b (n pxt),a (n pxt),xw (n pxt),zlevel( npxt), ndisp) do ip=n pxt-1,1,-1NR=0zmin=zlevel(ip+1)+dz30call area( npoin t(ip),x(1,1,ip),x(1,2,ip),b(ip),a(ip),xw (n pxt),zm in,n disp)if(A(i p) O)the nfr=(q(ip)/a(ip)*2.0*b(ip)/(
免费区 
