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

2、。泥沙條件：包括懸移質(zhì)，推移質(zhì)的均勻沙模型，推移質(zhì)計(jì)算模式為飽和輸沙，懸移質(zhì)計(jì)算模式為不飽和輸沙，水流泥沙方程采用非耦合解。二、 基本原理1、 基本方程水流連續(xù)方程： 水流運(yùn)動(dòng)方程 或 泥沙連續(xù)方程 河床變形方程 推移質(zhì)平衡輸沙方程G=G* 水流挾沙力公式采用張瑞瑾公式，推移質(zhì)輸沙率公式采用Mayer-_Peter公式，MAYER-PETER公式中的能坡J按均勻流曼寧公式近似計(jì)算（每個(gè)斷面不同）。2、 方程離散方程 在恒定流情況下有，離散為:Q=const 方程 變形為 或 上式離散為方程（4）去掉時(shí)間項(xiàng)得到該方程的解析解為：由方程（4-5）可得 對(duì)2 號(hào)斷面以下，上式可以離散為：對(duì)于進(jìn)口斷面

3、，推移質(zhì)不考慮，懸移質(zhì)采用單點(diǎn)離散方程（5）可離散為：3、 公式補(bǔ)充K取 0.124，m取1.05，干密度取1.3恢復(fù)飽和系數(shù) 均勻沙粒徑為d=0.041mm（懸移質(zhì)），d=2 mm（推移質(zhì)）三、 計(jì)算步驟1、輸入河床地形糙率等數(shù)據(jù)求得斷面面積與水位的關(guān)系（AZ），進(jìn)而求得斷面平均流速 ，水力學(xué)半徑 2、讀入一個(gè)時(shí)段的水沙數(shù)據(jù)（特別注意，不要一次性將數(shù)據(jù)全部讀入）讀入第一時(shí)段（Q,S）值3、計(jì)算水面線，同時(shí)得到各斷面的水力要素求得各個(gè)斷面的河寬、斷面面積、水深、平均流速等值計(jì)算前要注意在壩前輸入水位，各斷面均應(yīng)對(duì)流量賦值4、計(jì)算懸移質(zhì)水流挾沙力K取 0.124，m取1.05。5、計(jì)算推移質(zhì)輸沙

4、率（采用mayer-peter公式）6、計(jì)算各斷面含沙量公式7、 計(jì)算各斷面沖淤厚度進(jìn)口斷面 8、 修改水各斷面水下河床高程9、重新進(jìn)入（2）進(jìn)行下一循環(huán)10、計(jì)算10年河床變形，計(jì)算時(shí)段為一天，單位為秒（s）11、淤積總量年輸出一次，其余每兩年輸出一次計(jì)算結(jié)果四、 計(jì)算框圖開始讀入地形資料和糙率輸出初始庫容、深泓等i=i+1讀入一個(gè)QYQ=0計(jì)算水流子程序計(jì)算S*和Gbj=j+1計(jì)算S(j)和dy(j)YDy(i)=0,重新計(jì)算S(j) 判斷庫尾沖刷與否YJ<npnxt修改河床地形判斷年份輸出結(jié)果YI<3653結(jié)束五、 計(jì)算結(jié)果1、歷年輸沙量特征值計(jì)算各年懸疑質(zhì)輸沙量、推移質(zhì)輸沙

5、量、年均流量見表1。（其中推移質(zhì)輸沙量是按照懸疑質(zhì)輸沙量的1.5%計(jì)算的）。表1龍開口壩址歷年流量輸沙量特征值代表年年平均懸移質(zhì)輸沙量 （萬t）年平均推移質(zhì)輸沙量 （萬t）年平均流量 （萬t）1984294144.1213301985412061.8017501986322148.3213601987482872.4217601988254338.1514401989642696.3919001990438565.7818001991597289.5818701992198229.7312901993645696.841900最大值645696.841900最小值198229.731290多年

6、平均值428764.3116402、各年淤積總量表2 歷年淤積總量的變化年懸移質(zhì)輸沙量(萬噸)年均流量(萬噸)總淤積量(億立方米)000011821.494845560.1425291.20510994180.387316993.5419658070.879420027.3424697751.085521668.3327656171.204623853.9830875691.368730155.0537460621.701834879.5943121571.987936754.1846403382.1171039289.6650197622.243 圖1 歷年淤積總量的變化3、各年水位庫容關(guān)系

7、表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.6280.2690.1792591.3841.1260.6680.5630.2370.1562571.2471.01

8、40.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 各年水位庫容關(guān)系曲線4、 水面

9、線的變化表4 水面線變化年距壩里程(m)水位(m)年距壩里程(m)水位(m)年距壩里程(m)水位(m)034920267.6434920267.88834920274.43033810267.16433810267.7683381027401431350267.53831350273.43029280267429280267.37829280272.93027000267427000267.21827000272.4024050267424050267.07824050271.5021846267421846267.02821846270.53019896267419

10、896267819896269.58018270267418270267818270268.92015790267415790267815790268.1013170267413170267813170267.41010170267410170267810170267.05078202674782026787820267049502674495026784950267021002674210026782100267002674026780267年距壩里程(m)水位(m)年距壩里程(m)水位(m)年距壩里程(m)水位(m)234920267011034920273.99

11、233810267.12633810271.751033810273.78231350267.06631350271.011031350273.49229280267.02629280270.371029280273.18227000267627000269.631027000272.79224050267624050268.641024050272.22221846267621846267.941021846271.71219896267619896267.391019896271.07218270267618270267.111018270270.482157902676157902671

12、015790269.72131702676131702671013170268.942101702676101702671010170267.882782026767820267107820267.19249502676495026710495026722100267621002671021002672026760267100267圖3 水面線變化曲線5、 深泓變化表5 逐年深泓變化年距壩里程(m)深泓(m)年距壩里程(m)深泓(m)年距壩里程(m)深泓(m)034920264.1434920264.77834920271.75033810261.7433810264.12833810270.

13、46031350253.3431350262.29831350268.95029280249.4429280261.92829280268.6027000246.6427000263.23827000269.21024050241.3424050261.99824050267.71021846239.2421846261.5821846267.87019896237.9419896261.498198962681418270252.83818270266.84015790230.1415790250.82815790264.62013170228.2413170240

14、.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.926313

15、50267.251031350269.37229280259.53629280266.71029280269.62227000259.48627000267.211027000270.09224050250.62624050265.391024050268.98221846247.11621846265.191021846269.24219896242941019896269.7218270242.84618270264.711018270268.61215790233.11615790258.711015790267.03213170232.41613170249.

16、521013170266.88210170226.49610170242.18101017026527820224.1167820234.07107820264.7924950220.6564950228.12104950245.6322100216.8662100222.71102100236.5420213.7860219.7100229.24圖4 深泓變化曲線6、 壩前斷面變化該水庫為三角形淤積，最靠近壩前斷面淤積并不會(huì)特別嚴(yán)重，為使斷面變化情況表現(xiàn)的更加直觀，這里采用IP-17斷面作為典型斷面并分析其淤積情況。表6 壩前斷面變化年起點(diǎn)距(m)高程(m)年起點(diǎn)距(m)高程(m)年起點(diǎn)距(

17、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.472

18、61.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

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

20、0.34293.1300.3年起點(diǎn)距(m)高程(m)年起點(diǎn)距(m)高程(m)年起點(diǎn)距(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.2210

21、44.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.

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

23、279.1304.526293.1300.38293.1300.310293.1300.3圖5 壩前斷面變化曲線六、 結(jié)果分析1、 剖面形態(tài)分析A、縱斷面淤積情況分析：由深泓變化曲線圖可以看出，隨著時(shí)間的增加，河道深泓點(diǎn)高程成增高趨勢。該水庫運(yùn)行10年后，在距壩約10000米處，深泓點(diǎn)高程增加40多米，計(jì)算結(jié)果表明，該水庫泥沙以淤積為主，致使水庫庫容逐年損失，且情況較為嚴(yán)重。B、壩前斷面變化曲線分析：由壩前斷面變化曲線圖可以看出，隨著時(shí)間的增加，壩前斷面各點(diǎn)的總體高程也成增高趨勢。該水庫運(yùn)行10后，該斷面河床抬高20多米，表明在該處河道逐年以淤積為主，且淤積情況較為嚴(yán)重。2、 庫容損失合理性分

24、析由水位庫容曲線可以看出，水庫在運(yùn)行10年后，經(jīng)計(jì)算可知，在273m水位條件下，庫容的損失率為77.89%，在241m的水位條件下，庫容的損失率高達(dá)95.92%，計(jì)算結(jié)果表明，經(jīng)過十年的運(yùn)行，水庫庫容基本成全部淤滿狀態(tài)，所以在10年后將不能再繼續(xù)使用。結(jié)果表明，該水庫的調(diào)節(jié)能力較差。上述計(jì)算結(jié)果與現(xiàn)實(shí)不太相符，其原因是在計(jì)算該水庫庫容的泥沙損失時(shí)，計(jì)算沖往下游的泥沙量不合理。在給定條件下，該水庫并沒有考慮水位的變化及壩前沖沙情況。該水庫的計(jì)算條件為正常蓄水位保持不變，在高水位運(yùn)行條件下，必然會(huì)導(dǎo)致水庫的嚴(yán)重淤積，而實(shí)際情況是，在洪水來臨時(shí)，壩前水位要降至防洪限制水位，以增大防洪庫容，該條件下，

25、水流的速度及水流挾沙力較大，會(huì)造成庫區(qū)泥沙處于沖刷狀態(tài)，以增加水庫庫容。七、 計(jì)算程序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)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),alow

26、0(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='逐年累計(jì)淤積量.txt',status='unknown')open(13,file='逐

27、年深泓變化.txt',status='unknown')open(14,file='逐年水面線變化.txt',status='unknown')open(15,file='逐年壩前斷面淤積變化.txt',status='unknown')open(16,file='逐年水位庫容關(guān)系變化.txt',status='unknown')!*讀入地形數(shù)據(jù)*call Qbed(mm,npxt,x,dxa,alow0,npoint)!*計(jì)算沉降速度*W=FW(T,r,rs,d,0)wri

28、te(*,*)W!*計(jì)算斷面間距*do i=1,npxt-1dx(i)=dxa(i)-dxa(i+1)enddo!*計(jì)算初始水面線、深泓、庫容水位關(guān)系*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

29、,x,npoint,Zlevel,alow,dxa,dx,mm,nn,dWQST,AVQ,dWQS,dWGb,ndisp)!*計(jì)算各個(gè)時(shí)段水庫的淤積情況*DO 40,I=1,NDWRITE(*,*)'讀入天數(shù):',I!*讀入流量、沙量數(shù)據(jù)*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)!*計(jì)算各斷面的水流挾沙力和推移質(zhì)輸沙率*CALL SGB(nn,a,b,

30、Q,Sx,Gb,rough,w,ndisp)!*計(jì)算各斷面的含沙量和沖淤厚度*CALL DEPO(nn,mm,w,x,S,Sx,Gb,dt,B,Q,dx,dy,alow,alow0,npoint,faidy,ndisp) !*檢驗(yàn)質(zhì)量守恒*CALL CHECK(dy,B,S,Gb,dt,npxt,Zlevel,x,npoint,Q,dx,ndisp)41CONTINUE!*計(jì)算懸移質(zhì)、推移質(zhì)年均流量*WQS1=Q(1)*S(1)Gb(1)=WQS1*0.015WQSN=Q(nn)*S(nn)dWQST=dWQST+WQS1*dt/10000000 !懸移質(zhì)累計(jì)輸沙量dWQS=dWQS+(WQS

31、1-WQSN)*dt/1300.0/100000000 !懸移質(zhì)累計(jì)沖淤量dWGb=dWGb+(Gb(1)-Gb(nn)*dt/1300.0/100000000 !推移質(zhì)累計(jì)淤積量AVQ=AVQ+Q(1)*dt/10000 !累計(jì)流量!*輸出相關(guān)文件*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)

32、close(14) close(15) close(16) end!*讀入地形數(shù)據(jù)*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!*計(jì)算深泓*alow0(i)=x(1,2,i)do m=2,npoint(i)if(alow0(i)>x(m,2,i)

33、thenalow0(i)=x(m,2,i)endifenddo enddoend!*泥沙沉降速度的計(jì)算,采用張瑞瑾公式*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.000221*T*T CMU=X/

34、A RETURN END!*水面線及各斷面水力要素計(jì)算函數(shù)*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

35、 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)

36、,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!*計(jì)算斷面要素*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(i

37、p+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-dz1goto 30 endif call bisec(zmin,zmax,fmin,zlev

38、el(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

39、 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!*計(jì)算斷面要素*SUBROUTINE AREA(NPOINT,X,Z,B,A,XW,ZLEVEL,NDISP) DIMENSION X(NPOINT),Z(NPOINT) IF(NDISP&l

40、t;=-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

41、.0 XW=0.0 DO I=1,NPOINT-1ZMIN=AMIN1(Z(I),Z(I+1) IF(ZMIN<ZLEVEL)THENZMAX=AMAX1(Z(I),Z(I+1) IF(ZMAX<ZLEVEL)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(DB<0.0)THENW

42、RITE(*,*)'THE DISTANCE FOR NODE ',I,'AND',I+1,'ARE FALSE'WRITE(*,*)'IN AREA',I,I+1,X(I),X(I+1) STOP ENDIFDA=DB*DHB=B+DBA=A+DAXW=XW+DXENDIF ENDDO RETURN END!*二分法計(jì)算*subroutine bisec(zmin,zmax,fmin,zlevel,npoint,x,z,b,a,xw,q,ndisp,zlelo,blower,alower,qlower,dx,rough,fail

43、ev) dimension x(npoint),z(npoint) IF(NDISP.EQ.-1)WRITE(*,*)'INTO BISEC' ERR=0.00110 ddz=zmax-zminzlevel=0.5*(zmin+zmax) call area(npoint,x,z,b,a,xw,zlevel,ndisp) f=flevel(zlelo,zlevel,dx,qlower,q,rough,blower,b,alower,a,failev,ndisp) if(ddz>err)thenff=f*fmin if(ff>0)thenzmin=zlevel els

44、e zmax=zlevel endif ddz=zmax-zmin goto 10 endif return end!懸移質(zhì)挾沙力和推移質(zhì)輸沙率計(jì)算，懸移質(zhì)采用張瑞謹(jǐn)挾沙力公式!推移質(zhì)采用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 SGB'DO I=1,nnH=A(I)/B(I)U=Q(I)/A(I)VA1=U*3.0VA2=9.8*H*WVA3=(VA

45、1/VA2)*1.05Sx(I)=VA3*0.124 !懸移質(zhì)挾沙力!推移質(zhì)輸沙率，斷面輸沙率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 FGb'A1=(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