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

1、水庫一維泥沙 淤積計(jì)算課程設(shè)計(jì)目錄一、目的與要求1二、基本原理11、基本方程12、方程離散13、公式補(bǔ)充2三、計(jì)算步驟2四、計(jì)算框圖3五、計(jì)算結(jié)果41、歷年輸沙量特征值 42、各年淤積總量53、各年水位庫容關(guān)系 64、水面線的變化 75、深泓變化86、壩前斷面變化 9六、結(jié)果分析121、剖面形態(tài)分析 122、庫容損失合理性分析 12七、 計(jì)算程序13、目的與要求具備一定并編制出完整通過課程設(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ù)方程:AQx0QQ2txAQQ2txAtSAcxy°Gb水流運(yùn)動(dòng)方程泥沙連續(xù)方程河床變形方程xgA0 thgAgA i。 ifx2 2Q n c g4-0AR3(S S*)(S S)推移質(zhì)平衡輸沙方程G=G水流挾沙力公式采用張瑞瑾公式，推移質(zhì)輸沙率公式采用Mayer-_Peter公式,MAYER-PETER公式中的能坡J按均勻流曼寧公式近似計(jì)算（每個(gè)斷面不同）。2、方程離散方程在恒定流情況下有q 0，離散為：Q

3、=constx方程變形為£_ QgA x AQ2n2a2r"3Q2x 2gA2 2Q n24A R3上式離散為1 Lj22gQj2A2Zj 1 Zj2xn(AjQj4R/3j 1方程S（4）去掉時(shí)間項(xiàng)得到Sj i S j i SjS*jexpS*jSji-1 exp xBQx對(duì)2號(hào)斷面以下，上式可以離散為:由方程（4-5 ）可得QS' y00B 0BG bBGb QS QSxx'Bjyoj(10)B j 1 y0j 1t一（S S） q該方程的解析解為：對(duì)于進(jìn)口斷面，推移質(zhì)不考慮，懸移質(zhì)采用單點(diǎn)離散方程（5）可離散為：y。(S SJ t103、公式補(bǔ)充S3m

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

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

6、床地形五、計(jì)算結(jié)果1、歷年輸沙量特征值1。（其中推移質(zhì)輸沙量是按照計(jì)算各年懸疑質(zhì)輸沙量、推移質(zhì)輸沙量、年均流量見表懸疑質(zhì)輸沙量的1.5%計(jì)算的）。表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年擁漲召年濟(jì)亂8毎儼泓迪年探泓1Q年水廚線200 11010QD020

18、0005000040000圖4深泓變化曲線6、壩前斷面變化該水庫為三角形淤積，最靠近壩前斷面淤積并不會(huì)特別嚴(yán)重，為使斷面變化情況表現(xiàn)的更加直觀，這里采用IP-17斷面作為典型斷面并分析其淤積情況。表6壩前斷面變化年起點(diǎn)距高程年起點(diǎn)距高程年起點(diǎn)距高程(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年起點(diǎn)距高程年起點(diǎn)距高程年起點(diǎn)距高程(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、縱斷面淤積情況分析：由深泓變化曲線圖可以看出，隨著時(shí)間的增加，河道深泓點(diǎn)高程成增高趨勢(shì)。該水庫運(yùn)行10年后，在距壩約10000米處，深泓點(diǎn)高程增加 40多米，計(jì)算結(jié)果表明，該水庫泥沙以 淤

25、積為主，致使水庫庫容逐年損失，且情況較為嚴(yán)重。B壩前斷面變化曲線分析：由壩前斷面變化曲線圖可以看出， 隨著時(shí)間的增加，壩前斷面各點(diǎn)的總體高程也成增高 趨勢(shì)。該水庫運(yùn)行10后，該斷面河床抬高 20多米，表明在該處河道逐年以淤積為主，且淤 積情況較為嚴(yán)重。2、庫容損失合理性分析由水位庫容曲線可以看出，水庫在運(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ì)算該

26、水庫庫容的泥沙損失時(shí),計(jì)算沖往下游的泥沙量不合理。在給定條件下，該水庫并沒有考慮水位的變化及壩前沖沙情況。該水庫的計(jì)算條件為正常蓄水位保持不變,在高水位運(yùn)行條件下， 必然會(huì)導(dǎo)致水庫的嚴(yán)重淤積， 而實(shí)際情況是，在洪水來臨時(shí)，壩前水位要降至防洪限制水位，以增大防洪庫容，該條件下,水流的速度及水流挾沙力較大，會(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,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='逐年累計(jì)淤積量.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!*計(jì)算初始水面線、深泓、庫容水位關(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)!*計(jì)算各個(gè)時(shí)段水庫的淤積情況*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)!*計(jì)算各斷面的水流挾沙力和推移質(zhì)輸沙率*CALL SGB( nn ,a,b,Q,Sx,Gb,rough,w, ndisp)!*計(jì)算各斷面的含沙量和沖淤厚度*CALL DEPO( nn ,mm,w,x,S,Sx,Gb,dt,B,Q,dx,dy,alow,alowO, npoi nt,faidy, ndisp)

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

33、/10000!累計(jì)流量!*輸出相關(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計(jì)算深泓*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!*泥沙沉降速度的計(jì)算，采用張瑞瑾公式*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水面線及各斷面水力要素計(jì)算函數(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計(jì)算斷面要素*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計(jì)算斷面要素*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二分法計(jì)算*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ì)輸沙率計(jì)算，懸移質(zhì)采用張瑞謹(jǐn)挾沙力公式唯移質(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!*斷面沖淤厚度計(jì)算*SUBROUTINE DEPO( nn ,mm,w,x,S,Sx,Gb,dt,B,Q,dx,dy,alow,alowO, np