河流模擬課程設(shè)計—水庫一維泥沙淤積計算

1、水庫一維泥沙 淤積計算課程設(shè)計目錄一、目的與要求1二、基本原理11、基本方程12、方程離散13、公式補充2三、計算步驟2四、計算框圖3五、計算結(jié)果41、歷年輸沙量特征值 42、各年淤積總量53、各年水位庫容關(guān)系 64、水面線的變化 75、深泓變化86、壩前斷面變化 9六、結(jié)果分析121、剖面形態(tài)分析 122、庫容損失合理性分析 12七、 計算程序13、目的與要求具備一定并編制出完整通過課程設(shè)計，初步掌握一維數(shù)學模型建立數(shù)學模型的基本過程和計算方法, 的解決實際問題的能力。以水流、泥沙方程為基礎(chǔ)，構(gòu)建恒定流條件下的河道一維水沙數(shù)學模型, 的計算程序，并以某個水庫為實例，進行水庫泥沙淤積計算。水流

2、條件：恒定非均勻流。泥沙條件：包括懸移質(zhì)，推移質(zhì)的均勻沙模型，推移質(zhì)計算模式為飽和輸沙，懸移質(zhì)計 算模式為不飽和輸沙，水流泥沙方程采用非耦合解。二、基本原理1、基本方程水流連續(xù)方程:AQx0QQ2txAQQ2txAtSAcxy°Gb水流運動方程泥沙連續(xù)方程河床變形方程xgA0 thgAgA i。 ifx2 2Q n c g4-0AR3(S S*)(S S)推移質(zhì)平衡輸沙方程G=G水流挾沙力公式采用張瑞瑾公式，推移質(zhì)輸沙率公式采用Mayer-_Peter公式,MAYER-PETER公式中的能坡J按均勻流曼寧公式近似計算（每個斷面不同）。2、方程離散方程在恒定流情況下有q 0，離散為：Q

3、=constx方程變形為£_ QgA x AQ2n2a2r"3Q2x 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 QSxx'Bjyoj(10)B j 1 y0j 1t一（S S） q該方程的解析解為：對于進口斷面，推移質(zhì)不考慮，懸移質(zhì)采用單點離散方程（5）可離散為：y。(S SJ t103、公式補充S3m

4、k gRK取0.124 , m取 1.05，干密度0 取 1.3恢復飽和系數(shù)0.25均勻沙粒徑為d=0.041mm （懸移質(zhì)），d=2 mm （推移質(zhì)）二、計算步驟1、輸入河床地形糙率等數(shù)據(jù)求得斷面面積與水位的關(guān)系（A Z）,進而求得斷面平均流速Q(mào)，水力學半徑A2、讀入一個時段的水沙數(shù)據(jù)（特別注意，不要一次性將數(shù)據(jù)全部讀入）讀入第一時段（Q,S）值3、計算水面線，同時得到各斷面的水力要素求得各個斷面的河寬、斷面面積、水深、平均流速等值計算前要注意在壩前輸入水位，各斷面均應對流量賦值4、計算懸移質(zhì)水流挾沙力gRK取 0.124，m取 1.05。5、計算推移質(zhì)輸沙率（采用 mayer-peter公

5、式）3/2n 3/2（）hJ 0.047（ s ）d ngb0.125 1/2（）gs6、計算各斷面含沙量公式Sj1 S*j1 Sj Sj exp 2S*j Sj1expq x7、計算各斷面沖淤厚度進口斷面y。一（S1 SJ t08、修改水各斷面水下河床高程H H。 丫9、重新進入（2）進行下一循環(huán)10、計算10年河床變形，計算時段為一天，單位為秒（ s）11、淤積總量年輸出一次，其余每兩年輸出一次計算結(jié)果四、計算框圖讀入地形資料和糙率輸出初始庫容、深泓等YQ=0計算S和Gb計算水流子程序計算S(j)和dy(j)1r4i=i+11F讀入一個判斷庫尾沖刷與否YJvnpnxt重新計算S(j)修改河

6、床地形五、計算結(jié)果1、歷年輸沙量特征值1。（其中推移質(zhì)輸沙量是按照計算各年懸疑質(zhì)輸沙量、推移質(zhì)輸沙量、年均流量見表懸疑質(zhì)輸沙量的1.5%計算的）。表1龍開口壩址歷年流量輸沙量特征值代表年年平均懸移質(zhì)年平均推移質(zhì)年平均流量輸沙量（萬t）輸沙量（萬t）（萬t）1984294144.1213301985412061.8017501986322148.3213601987482872.4217601988254338.1514401989642696.3919001990438565.7818001991597289.5818701992198229.7312901993645696.841900取大

7、值645696.841900最小值198229.731290多年平均值428764.3116402、各年淤積總量表2歷年淤積總量的變化年懸移質(zhì)輸沙 量（萬噸）年均流量（萬噸）總淤積量（億 立方米）000011821.494845560.1425291.20510994180.387316993.5419658070.879420027.3424697751.085521668.3327656171.204623853.9830875691.368730155.0537460621.701834879.5943121571.987936754.1846403382.1171039289.6650

8、197622.243S 2 5 1 S 甲2 100 2 4 6 8 10圖1歷年淤積總量的變化3、各年水位庫容關(guān)系表3各年水位庫容關(guān)系水位（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.6

9、280.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.05

10、30.0252410.4170.3470.1590.1370.0380.01727U2US26020245圖2各年水位庫容關(guān)系曲線4、水面線的變化表4水面線變化年距壩里程水位年距壩里程水位年距壩里程水位(m)(m)(m)(m)(m)(m)034920267.6434920267.88834920274.43033810267.16433810267.7683381027401431350267.53831350273.43029280267429280267.37829280272.93027000267427000267.21827000272.4024050267

11、424050267.07824050271.5021846267421846267.02821846270.53019896267419896267819896269.58018270267418270267818270268.92015790267415790267815790268.1013170267413170267813170267.41010170267410170267810170267.05078202674782026787820267049502674495026784950267021002674210026782100267002674026780267年距壩里程水位年

12、距壩里程水位年距壩里程水位(m)(m)(m)(m)(m)(m)234920267011034920273.99233810267.12633810271.751033810273.78231350267.06631350271.011031350273.49229280267.02629280270.371029280273.18227000267627000269.631027000272.79224050267624050268.641024050272.22221846267621846267.941021846271.71219896267619896267.

13、391019896271.07218270267618270267.111018270270.482157902676157902671015790269.72131702676131702671013170268.942101702676101702671010170267.882782026767820267107820267.19249502676495026710495026722100267621002671021002672026760267100267圖3水面線變化曲線5、深泓變化表5逐年深泓變化年距壩里程(m)深泓(m)年距壩里程(m)深泓(m)年距壩里程(m)深泓(m)034

14、920264.1434920264.77834920271.75033810261.7433810264.12833810270.46031350253.3431350262.29831350268.95029280249.4429280261.92829280268.6027000246.6427000263.23827000269.21024050241.3424050261.99824050267.71021846239.2421846261.5821846267.87019896237.9419896261.498198962681418270252.8381

15、8270266.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.961

16、034920271.76233810262.41633810269.081033810270.93231350259.92631350267.251031350269.37229280259.53629280266.71029280269.62227000259.48627000267.211027000270.09224050250.62624050265.391024050268.98221846247.11621846265.191021846269.24219896242941019896269.7218270242.84618270264.711018270

17、268.61215790233.11615790258.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 11010QD020

18、0005000040000圖4深泓變化曲線6、壩前斷面變化該水庫為三角形淤積，最靠近壩前斷面淤積并不會特別嚴重，為使斷面變化情況表現(xiàn)的更加直觀，這里采用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

19、.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.

20、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.444

21、241.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.2267

22、618.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.

23、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.726

24、76241.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壩前斷面變化曲線六、結(jié)果分析1、剖面形態(tài)分析A、縱斷面淤積情況分析：由深泓變化曲線圖可以看出，隨著時間的增加，河道深泓點高程成增高趨勢。該水庫運行10年后，在距壩約10000米處，深泓點高程增加 40多米，計算結(jié)果表明，該水庫泥沙以 淤

25、積為主，致使水庫庫容逐年損失，且情況較為嚴重。B壩前斷面變化曲線分析：由壩前斷面變化曲線圖可以看出， 隨著時間的增加，壩前斷面各點的總體高程也成增高 趨勢。該水庫運行10后，該斷面河床抬高 20多米，表明在該處河道逐年以淤積為主，且淤 積情況較為嚴重。2、庫容損失合理性分析由水位庫容曲線可以看出，水庫在運行10年后，經(jīng)計算可知，在273m水位條件下，庫容的損失率為77.89%，在241m的水位條件下，庫容的損失率高達95.92%，計算結(jié)果表明，經(jīng)過十年的運行，水庫庫容基本成全部淤滿狀態(tài)，所以在10年后將不能再繼續(xù)使用。結(jié)果表明，該水庫的調(diào)節(jié)能力較差。上述計算結(jié)果與現(xiàn)實不太相符,其原因是在計算該

26、水庫庫容的泥沙損失時,計算沖往下游的泥沙量不合理。在給定條件下，該水庫并沒有考慮水位的變化及壩前沖沙情況。該水庫的計算條件為正常蓄水位保持不變,在高水位運行條件下， 必然會導致水庫的嚴重淤積， 而實際情況是，在洪水來臨時，壩前水位要降至防洪限制水位，以增大防洪庫容，該條件下,水流的速度及水流挾沙力較大，會造成庫區(qū)泥沙處于沖刷狀態(tài)，以增加水庫庫容。七、計算程序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,fa

27、idy=0.5,dt=86400)dime nsion 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=&#

28、39;河床地形.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='逐年壩前斷面淤積變化.t

29、xt',status='unknown')open(16,file='逐年水位庫容關(guān)系變化.txt',status='unknown')讀入地形數(shù)*call Qbed(m m,n pxt,x,dxa,alow0 ,npoint)算沉降速*W=FW(T,r,rs,d,0)write(*,*)W算斷面間住巨 *do i=1, npxt-1dx(i)=dxa(i)-dxa(i+1)en ddo!*計算初始水面線、深泓、庫容水位關(guān)系*dWQST=0.0dWQS=0.0dWGb=0.0AVQ=0.0do l=1,NPXT rough(l)=0.03

30、5 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,0,x ,n poi nt,Zlevel,alow,dxa,dx,mm, nn ,dWQST,AVQ,dWQS,dWGb, ndisp)!*計算各個時段水庫的淤積情況*DO 40,I=1,NDWRITE(*,*)'讀入天數(shù):',I!*讀入流量、沙

31、量數(shù)據(jù)*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)!*計算各斷面的水流挾沙力和推移質(zhì)輸沙率*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,alowO, npoi nt,faidy, ndisp)

32、I *驗質(zhì)量守丿恒 *CALL CHECK(dy,B,S,Gb,dt, npxt,Zlevel,x, npoi nt,Q,dx, ndisp)41CONTINUE!* 計算懸移質(zhì)、推移質(zhì)年均流量*WQS 仁Q(1)*SGb(1)=WQS1*0.015WQSN=Q( nn 廣S( nn)dWQST=dWQST+WQS1*dt/10000000 !懸移質(zhì)累計輸沙量dWQS=dWQS+(WQS1-WQSN)*dt/1300.0/100000000!懸移質(zhì)累計沖淤量dWGb=dWGb+(Gb(1)-Gb( nn )*dt/1300.0/100000000!推移質(zhì)累計淤積量AVQ=AVQ+Q(1)*dt

33、/10000!累計流量!*輸出相關(guān)文件*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)close(12)close(13)close(14)close(15)close(16) endI *讀入地形數(shù)據(jù)*subrouti ne Qbed(m m,nn, x,dxa,alowO ,npoint)dime nsion x(mm,2,

34、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, 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,

35、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.000221*T*TCMU=X/ARETURNEND水面線及各斷面水力要素計算函數(shù)*subroutine level(x,rough,npxt,zlevel,dx,q,npoint,b,a,xw,nn,mm,failev,ND

36、ISP) 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

37、), 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)<=0)the nzmin=zmin+0.173goto 30en diffmin=flevel(zlevel(ip+1),zm in, dx(ip),q(ip+1),q(ip),rough(ip),b(ip+1),b(ip),a(ip+1),a(ip),failev, ndisp)if(fmin> O

38、)the nfr=(q(ip)/a(ip)*2.0*b(ip)/(9.8*a(ip)if(fr<1)the nzmax=zm in+dz計算斷面要素*20callarea( npoin t(ip),X(1,1,ip),x(1,2,ip),b(ip),a(ip),xw(ip),zmax ,n disp)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)the nzmax=zmax+dzgoto 20en dif

39、elsezmin=zmin+dz2NR=NR+1IF(NR>NC)THENzmin=zmin-dz1write(*,*)'the loop is death,pause'WRITE(*,*)' nr,ip,zmi n,fmin', nr,ip,zmi n, fminread(*,*)en difgoto 30en difelsezmin=zmin-dz1goto 30en difcallbisec(zm in, zmax,fm in, zlevel(ip), npoin t(ip),X(1,1,ip),x(1,2,ip),b(ip),a(ip),xw(ip

40、),q(ip), ndisp,zlevel(ip+1),b(ip+1),a(ip+1),q(ip+1),dx(ip),rough(ip),failev)en ddoreturnendFUNCTION flevel(zlelo,zlevel,dx,qlower,q,rough,blower,b,alower,a,failev, ndisp)IF(NDISPEQ.-1)WRITE(*,*)'INTO FLEVEL'HLOWER=ALOWER/BLOWERH=A/BAA=ZLELO-ZLEVELB仁failev*(QLOWER/BLOWER)*2.0/HLOWER*(10/ 3.0)

41、B2=(1-failev)*(Q/B)*2.0/H*(10/3.0)BB=DX*ROUGH*2.0*(B1+B2)C仁(QLOWER/ALOWER)*2.0C2=(Q/A)*2.0CC=(C1-C2)/(2*9.8)FLEVEL=AA+BB+CCRETURNEND計算斷面要素*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 AR

42、EA THE INPUT DATA ARE FALSE N=',Npoi ntSTOPENDIFIF(Z(1)<ZLEVEL.OR.Z(NPOINT)<ZLEVEL)THENWRITE(*,*)'IN AREA THE WATER LEVEL IS TOO HIGH OR THE HEIGHT OF THE SECTIONAT EDGE IS TOO LOWSTOPENDIFA=0.0B=0.0XW=0.0DO l=1,NPOINT-1ZMIN=AMIN1(Z(I),Z(I+1)IF(ZMIN<ZLEVEL)THENZMAX=AMAX1(Z(I),Z(I+1)

43、IF(ZMAX<ZLEVEL)THENDB=X(I+1)-X (I)DH=ZLEVEL-0.5*(Z(l)+Z(l+1)DX=(Z(l)-Z(l+1)*2.0+DB*DB)*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<0.0)THENWRITE(*,*)'THE DISTANCE FOR NODE ',I,'AND',I+1,'ARE FALSE'WRITE(*,*)'

44、;IN AREA',I,I+1,X(I),X(I+1)STOPENDIFDA=DB*DHB=B+DBA=A+DAXW=XW+DXENDIFENDDORETURNEND二分法計算*subrouti nebisec(zmi n,zmax,fmi n,zlevel, npoi nt,x ,z, b,a,xw,q, ndisp,zlelo,blower,alower,qlower,dx,rough,failev) dime nsion x(npoin t),z( npoint)IF(NDISFEQ.-1)WRITE(*,*)'INTO BISEC'ERR=0.001 10 dd

45、z=zmax-zm inzlevel=0.5*(zm in+zmax)call area( npoin t,x,z,b,a,xw,zlevel, ndisp)f=flevel(zlelo,zlevel,dx,qlower,q,rough,blower,b,alower,a,failev, ndisp) if(ddz>err)the nff=f*fminif(ff>0)the nzmin=zlevelelsezmax=zlevelen dif ddz=zmax-zm in goto 10en difreturnend!懸移質(zhì)挾沙力和推移質(zhì)輸沙率計算，懸移質(zhì)采用張瑞謹挾沙力公式唯移質(zhì)采

46、用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(n disp=-1)WRITE(*,*)'INTO SGB'H=A(I)/B(I)U=Q(I)/A(I)VA仁 U*3.0VA2=9.8*H*WVA3=(VA1VA2)*1.05Sx(l)=VA3*0.124!懸移質(zhì)挾沙力!推移質(zhì)輸沙率，斷面輸沙率Gb(l)=FGb(rough(l),10.0,H,Q(l),A(l),26.5,0

47、.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 FGb'A1=(rn*2.0)*(Q*2.0)A2=(A*2.0)*(H*(4 /3.0)RJ=A1A2rnn=(d*(1 /6.0)/26.0AA=(rn n/rn)*1.5)*rs*H*RJBB=0.047*(rs-r)*dCC=(rous-rou)/rous)*9.8*(rou*0.5)*0.125IF(AA-BB)<O)THENFGB=0ELSEFGB=(AA-BB)*1.5)/CCENDIFRETURNEND!*斷面沖淤厚度計算*SUBROUTINE DEPO( nn ,mm,w,x,S,Sx,Gb,dt,B,Q,dx,dy,alow,alowO, np